Subduction termination through progressive slab deformation across Eastern Mediterranean subduction zones from updated P-wave tomography beneath Anatolia

The interpretation of fast seismic velocity anomalies as subducted lithosphere from images in the lower mantle using seismic tomography is frequently used to inform tectonic reconstructions of past subduction and infer properties of mantle dynamics such as slab sinking rates, lower-mantle viscosity, and ambient mantle flow (Van der Voo et al., 1999a; Hafkenscheid et al., 2006; Sigloch et al., 2008; van der Meer et al., 2010; Domeier et al., 2016; Wu et al., 2016). Such interpretations assume that subducted lithosphere remains stiff for up to hundreds of millions of years as it descends in the mantle before settling at the core-mantle boundary in the “slab graveyard” (Wysession, 1996; Rao and Kumar, 2014). However, slabs do not always remain cohesive, but instead can tear and segment as they descend into the mantle (Wortel and Spakman, 2000; Nolet, 2009; Biryol et al., 2011). This process is particularly notable in regions of continental collision through slab detachment (von Blanckenburg and Davies, 1995; Van der Voo et al., 1999b; Duretz et al., 2014).

The Eastern Mediterranean marks a major transition in the character of African-Arabian-Eurasian convergence from collision-dominated convergence in the east to subduction-dominated convergence in the west, making it an ideal region in which to investigate slab deformation during the latter stages of subduction. With this paper, we present some of the best-resolved seismic images to date of the Eastern Mediterranean subduction zones down to 900 km depth using a composite seismic array that includes data from the recent Continental Dynamics–Central Anatolian Tectonics (CD-CAT) seismic deployment in central Anatolia. Driven by these new images, we interpret the character of the modern subducted lithosphere beneath the entire Anatolian subcontinent in the context of Anatolian tectonics and discuss its implications for the fate of slabs in the mantle.

Anatolia formed from the amalgamation of major lithospheric blocks during a long history of Tethyan subduction that dates back at least to the Cretaceous (Şengör and Yilmaz, 1981; Menant et al., 2016; van Hinsbergen et al., 2016). Today, Anatolia is in a complex tectonic regime that can be described in three discrete neotectonic regions: the East Anatolian Contractional Province, the West Anatolian Extensional Province, and the Central Anatolian Province (Şengör et al., 1985). These neotectonic regions roughly correlate with the three convergent domains used later herein to describe our results: the Bitlis-Zagros domain, the Aegean domain, and the Cyprean domain (Fig. 1).

In eastern Anatolia, the Bitlis-Zagros suture zone marks the modern boundary between the Arabian and Eurasian plates (Fig. 1). The Bitlis-Zagros suture zone formed in the mid-Miocene when the Arabian and Eurasian plates collided upon closure of the southern branch of the Neo-Tethys ocean (Şengör and Yilmaz, 1981; Okay et al., 2010; Nouri et al., 2016). Around 11–10 Ma, following Arabian-Eurasian collision, steepening and delamination of the subducted Neo-Tethyan slab from beneath eastern Anatolia resulted in an overriding crust with a thin or missing mantle lithosphere (Şengör et al., 2003; Angus et al., 2006). Subsequent detachment opened a slab window, allowing for the upwelling of hot, fertile asthenosphere beneath what is now the high-elevation East Anatolian Plateau (Keskin, 2003, 2007; Şengör et al., 2003). Asthenospheric upwelling led to late Miocene to recent broad, domal uplift (Şengör et al., 2003) and enriched-mantle volcanism across the plateau (Pearce et al., 1990; Keskin, 2003) and continues to support the undercompensated crust of the plateau (Özacar et al., 2010). Today, convergence continues at a rate of 15 mm/yr (Reilinger et al., 2006), which is accommodated both by north-south–oriented shortening (Şengör et al., 2003) and the westward extrusion of Anatolia along the North Anatolian and East Anatolian faults (Barka and Reilinger, 1997).

The modern tectonic regime in western Anatolia is controlled by subduction along the Aegean Trench. Since at least 35–30 Ma, the Aegean Trench has migrated southward (Le Pichon and Angelier, 1979; Jolivet et al., 2013), contributing to extension in the overriding plate and thinning of the Aegean continental lithosphere. However, several mechanisms may contribute to extension, including trench migration, tectonic escape, and orogenic collapse (McKenzie, 1978; Şengör et al., 1985; Seyítoǧlu and Scott, 1996; Reilinger et al., 2006; Jolivet et al., 2013). Today, the western Anatolian/Aegean crust is ∼25–30 km thick (Zhu et al., 2006; Karabulut et al., 2013) after ∼50% extension in the Aegean (McKenzie, 1978). An active arc and continuous Wadati-Benioff zone down to ∼200 km depth (Papazachos et al., 2000) indicate continued active subduction along the Aegean Trench, while global positioning system (GPS) data indicate continued southward trench retreat (Reilinger et al., 2006).

Due to the Miocene collision of Arabia with Eurasia, the North Anatolian fault zone and conjugate faults in central-eastern Anatolia formed, accommodating westward extrusion and counterclockwise rotation of the Anatolian plate (Barka and Reilinger, 1997; Şengör et al., 2005). Deformation within central Anatolia is dominated by transtension, as expressed by increased exhumation rates and extensional basin formation related to the Central Anatolian fault zone (Koçyiǧit and Beyhan, 1998; Fayon and Whitney, 2007; Higgins et al., 2015), which have facilitated the voluminous Miocene to recent volcanism of the Central Anatolian Volcanic Province (Toprak and Göncüoǧlu, 1993; Toprak, 1998; Kuscu and Geneli, 2010) within the interior of the Central Anatolian Plateau. Conversely, the Central Taurus Mountains, which bound the southern margin of the Central Anatolian Plateau, have undergone little internal deformation since the mid-Miocene, despite ∼2 km uplift since 8 Ma (Cosentino et al., 2012; Schildgen et al., 2012). Active deformation in central Anatolia occurs in the context of ongoing subduction along the Cyprean Trench. The Cyprean trench connects to the Aegean Trench to the west via the Pliny-Strabo subduction-transform edge propagator (STEP) system (Govers and Wortel, 2005) and to the Bitlis-Zagros suture to the east through the diffuse plate boundary beneath the Adana Basin (Abgarmi et al., 2017). The subducting material along the Cyprean trench is itself complex. The western edge of Cyprus marks a boundary between subducting continental lithosphere to the east and oceanic lithosphere to the west (Fig. 1; Granot, 2016). The arrival of two attenuated fragments of Gondwanan continental lithosphere, the Eratosthenes and Anaximander Seamounts, complicate subduction further (Fig. 1). The Eratosthenes Seamount is a Gondwana-derived rifted continental fragment actively underthrusting beneath Cyprus (Robertson, 1998) and potentially contributing to the shallow slab dip angle of the Cyprean slab (Biryol et al., 2011; Bakırcı et al., 2012). The Anaximander Mountains are being underthrust offshore of Anatolia’s southwestern coast beneath the Isparta angle, which is a unique, largely undeforming, independently moving block separating the central Anatolian (Cyprean) and Aegean convergent domains (Barka and Reilinger, 1997; Zitter et al., 2003). While subduction of the Cyprean slab is a key control on the complexity of central Anatolian tectonics, a Wadati-Benioff zone that is largely restricted to the upper 100 km limits direct observations of Cyprean slab geometry.

Most attempts to understand the complex tectonics of Anatolia, as with the greater Mediterranean region, begin with African-Arabian subduction beneath Eurasia and its role in affecting mantle flow (Wortel and Spakman, 2000; Faccenna et al., 2014). For example, models proposed to explain the Eastern Anatolian and Central Anatolian Plateaus’ kinematic, magmatic, and/or uplift history invoke a variety of mechanisms, including propagating slab detachment (Wortel and Spakman, 2000; Cosentino et al., 2012; Schildgen et al., 2012, 2014), slab rollback–related delamination (Göğüş and Psyklywec, 2008; van Hinsbergen et al., 2010; Bartol and Govers, 2014), multiple contemporaneous subduction zones (van Hinsbergen et al., 2016; Gürer et al., 2016), and slab-induced mantle convection (Faccenna and Becker, 2010; Reid et al., 2017). Each proposed mechanism makes predictions about the current structure of subducted Tethyan lithosphere and the surrounding mantle structure that can only be probed with seismic imaging techniques.

Several studies have investigated the seismic structure of the mantle beneath Anatolia at a variety of scales. Studies of the uppermost mantle show velocities slower than the global average across the entire region (Gans et al., 2009; Mutlu and Karabulut, 2011; Salaün et al., 2012; Delph et al., 2015a; Govers and Fichtner, 2016). Evidence from Pn tomography (Gans et al., 2009; Mutlu and Karabulut, 2011), surface wave tomography (Bakırcı et al., 2012; Salaün et al., 2012; Delph et al., 2017), teleseismic tomography (Biryol et al., 2011; Salah, 2017), ambient noise tomography (Warren et al., 2013), full-waveform tomography (Fichtner et al., 2013; Govers and Fichtner, 2016), and joint receiver function–surface wave dispersion inversion (Gök et al., 2007; Delph et al., 2015b) all suggests that velocities are particularly slow beneath eastern Anatolia throughout the upper mantle, supporting the hypothesis of rising hot asthenosphere through a slab window (Keskin, 2003; Şengör et al., 2003). Slow velocities as far west as the Isparta angle have been included with those beneath the East Anatolian Plateau to represent inflowing asthenosphere resulting from recent large-scale slab delamination and rollback beneath central and eastern Anatolia (Bartol and Govers, 2014; Govers and Fichtner, 2016).

The subducting Aegean and Cyprean slabs have also been imaged at a variety of scales beneath Turkey. The shallow slabs are thought to be separated by a vertical STEP-related tear approximately below the Isparta angle (de Boorder et al., 1998; Biryol et al., 2011; Salaün et al., 2012). While the Aegean slab appears to dip at a typical subduction angle at shallow depths (Salaün et al., 2012), the Cyprean slab appears to subduct subhorizontally as far north as 38°N (Bakırcı et al., 2012; Abgarmi et al., 2017; Delph et al., 2017). Regional- and global-scale body wave tomography images show both slabs continuing through the mantle transition zone, where they appear to connect, thicken, and temporarily flatten (Bijwaard et al., 1998; Piromallo and Morelli, 2003; Biryol et al., 2011). However, while the Aegean slab appears continuous with a near-constant dip through the upper mantle in most models, the Cyprean slab anomaly appears more diffuse and irregular, suggesting the slab may be actively deforming and segmenting as it subducts (Biryol et al., 2011). Upon flattening in the mantle transition zone (Biryol et al., 2011), both slabs appear to continue into the lower mantle (Bijwaard et al., 1998; Piromallo and Morelli, 2003) as a result of the accumulation of a critical mass of oceanic lithosphere (Faccenna et al., 2003).

The abundance of seismic research in Anatolia provides an excellent framework within which to view the general character of Eastern Mediterranean slabs from the surface into the lower mantle, but methodological and data limitations in these studies inhibit their ability to resolve the finer-scale structure of the various slab segments beneath Anatolia, which are critical to evaluating the tectonic and geodynamic evolution of the Eastern Mediterranean system. More recent regional body wave tomography has improved upon this with the incorporation of frequency-dependent sensitivity kernels and a denser seismic network, providing images of tears within a once-continuous Tethyan slab and revealing the segmented nature of the subducting Cyprean slab (Biryol et al., 2011). However, this model lacks the resolution in central Turkey to delineate the structure, deformation, and segmentation of the Cyprean slab due to relatively sparse station density in the region at the time of the study. The model presented in this paper builds on the work of Biryol et al. (2011) through the addition of a dense data set in central Turkey and incorporation of data from the largely expanded national seismic network in Turkey (Fig. 2), allowing us to obtain a detailed image of the structure of the Cyprean slab and surrounding slabs down to 900 km depth.

The tomography model presented here serves as a significant update on the tomography model of Biryol et al. (2011). Next, we briefly review the data and methodology used in these studies while highlighting the relevant modifications made for the current study. For a more detailed description of the methodology and its limitations, see Biryol et al. (2011) and Schmandt and Humphreys (2010).

For this study, we utilized the data set of Biryol et al. (2011), which consists of 33,990 direct P and PKIKP phase arrival times from 409 unique earthquakes recorded at a variety of temporary (Sandvol, 1999; Mitchell and Zhu, 2002; Beck and Zandt, 2005) and permanent (Albuquerque Seismological Laboratory [ASL]/USGS, 1988; GEOFON Data Centre, 1993; Bogazici University Kandilli Observatory and Earthquake Research Institute, 2001) seismic networks. Arrival times in that study were picked on seismograms filtered in three frequency bands with corner frequencies of 0.2–0.8 Hz, 0.1–0.4 Hz, and 0.04–0.16 Hz for broadband stations and one frequency band with corner frequencies of 0.5–1.5 Hz for short-period stations.

We added data collected between 2013 and 2015 across two broadband seismic networks: the recent CD-CAT Project (72 stations; Sandvol, 2013) seismic array, and the Kandilli Observatory Digital Broadband Seismic Network (KOERI; 121 stations; Bogazici University Kandilli Observatory and Earthquake Research Institute, 2001; see Fig. 2). We picked arrivals from earthquakes with magnitude 5.0 and greater at distances between 30° and 90° from the stations for direct P phases (509 earthquakes) and between 155° and 180° from the stations for PKIKP phases (18 earthquakes; Fig. 3A). An excess of earthquakes from the west Pacific subduction zones led to an uneven back-azimuthal distribution (Fig. 3B), so to more evenly distribute the back-azimuthal distribution of rays, we excluded from our data set ∼47,000 rays from earthquakes in the west Pacific by preferentially removing earthquakes with fewer recovered residuals. Arrival times were picked on seismograms filtered in four frequency bands with corner frequencies of 0.5–1.5 Hz, 0.2–0.8 Hz, 0.1–0.4 Hz, and 0.04–0.16 Hz for all stations. After combining data sets, we used 66,331 arrival times (12% from bandpass 0.5–1.5 Hz, 34% from bandpass 0.2–0.8 Hz, 31% from bandpass 0.1–0.4 Hz, and 23% from bandpass 0.04–0.16 Hz), consisting of 65,352 direct phases and 979 core phases from 936 earthquakes over the time period between 1999 and 2015.

To obtain traveltime residuals, we used a multichannel cross-correlation algorithm (VanDecar and Crosson, 1990; Pavlis and Vernon, 2010). Relative traveltimes were determined with respect to the IASP91 one-dimensional (1-D) P-wave velocity model (Kennett and Engdahl, 1991) and de-meaned for each event. We additionally made a correction to the traveltime residual for crustal thickness beneath each station by assuming a vertical ray and constant velocity of 6.2 km/s in the lower crust. Crustal thicknesses were determined using the Moho map of Vanacore et al. (2013) and adjusted with more recent constraints from Abgarmi et al. (2017); see Supplemental Figure S1¹. Outside of these models, crustal thickness was manually extrapolated. Notably, the crustal correction had negligible effect on the resulting model. Average traveltime residuals for each station are summarized in Figure 4.

The method used in this study was the finite-frequency teleseismic tomography method (Dahlen et al., 2000; Schmandt and Humphreys, 2010). The finite-frequency method accounts for the frequency-dependent volumes surrounding the geometric ray path that are sampled by any given ray by approximating sensitivity with the Born theoretical “banana-doughnut” kernels for the first Fresnel zone (Dahlen et al., 2000). The inversion is a smoothed, damped least-squares inversion to account for uncertainty in ray geometry and to regularize the inversion. Trade-off analysis indicates that a range of damping and smoothing parameters are justified (Fig. S2 [^{footnote 1}]). We chose to adopt damping and smoothing parameters of 6 and 7, respectively, in accordance with the study of Biryol et al. (2011), for consistency and to facilitate comparison. Station and event terms were additionally included in the inversion to account for small-scale heterogeneity unaccounted for by our model geometry, such as local site effects.

The most significant changes to the Biryol et al. (2011) approach were in model space parameterization. The modeled region is 3300 km (E-W) by 2300 km (N-S) centered on 39°N, 35°E, with horizontal node spacing at the surface ranging from 35 km in the center to 56 km on the edges and dilating with depth to a maximum spacing ranging from 55 to 90 km at the base of the model space. Vertically, the model space extends from 60 km depth to 1010 km depth, with spacing ranging from 35 km near the surface to 70 km at the bottom. This is in contrast to Biryol et al. (2011), which extended from 35 km to 655 km. Model parameterization was chosen to optimize crossing-ray density as determined by hit quality, which quantitatively assesses the quantity and back-azimuthal distribution of rays sampling a given node (Schmandt and Humphreys, 2010; Biryol et al., 2011). The inclusion of 32,341 additional arrival times, primarily from the CD-CAT seismic array in central Anatolia, served to most significantly improve sampling of the upper mantle beneath central Anatolia and the transition zone and uppermost lower mantle across the study area, while having a minor effect on sampling of the upper mantle beneath western and eastern Anatolia (Fig. S3 [^{footnote 1}]). This new data allowed our modeled region to have sufficient hit quality (hit quality > 0.5) within the lower mantle to justify extending the model space to 1010 km depth (Fig. S3G [^{footnote 1}]). However, various tests suggest that the deepest two layers of the model space exhibit artificially high-amplitude velocity perturbations (Fig. S4 [^{footnote 1}]); thus, we only interpreted results to 900 km depth.

Data variance reduction is a commonly used metric to assess the fit of a model to the data for the tomographic inverse problem. For this study, we obtained a variance reduction of 70%. However, while variance reduction indicates fit to the data, there is no such quantitative indicator for uncertainty. Instead, we performed a series of synthetic anomaly recovery tests to qualitatively assess uncertainty and the model’s ability to recover structure.

First, we performed a standard checkerboard analysis, using alternating square synthetic input anomalies of ±5% separated by neutral zones of the same size (Figs. S5 and S6 [^{footnote 1}]). On average, the checkerboard structure is well recovered beneath Turkey, where up to 50% of checker amplitude is recovered. Nonzero amplitude recovery in neutral input layers (Fig. S5C [^{footnote 1}]) implies some degree of vertical smearing, while amplitude recovery between checkers indicates some lateral smearing induced by the smoothness constraint. Outside of Turkey (including beneath the Aegean Sea, the Black Sea, the Mediterranean Sea, and the Arabian plate), the checkerboard recovery drops significantly. In general, we found that anomaly geometry was recovered within the 0.5 hit quality contour, so we used this quantity to outline the well-resolved portion of our model in the results. However, this outline is not a deciding factor for model resolution; it was instead used as a proxy for the well-sampled portions of the model space.

Finally, in addition to the standard checkerboard analysis, isolated anomaly recovery tests were performed to assess the potential of our model to resolve synthetic structures and their geometries within the well-sampled portion of the model space (Figs. S7, S8, S9, and S10 [^{footnote 1}]). These recovery tests are described in the discussion section.

The bulk character of the tomography model (Fig. 5) indicates dominantly slow velocity perturbations throughout the upper mantle, consistent with previous imaging efforts (Gans et al., 2009; Salaün et al., 2012; Delph et al., 2015a; Govers and Fichtner, 2016). These slow velocity perturbations grade into dominantly fast velocity perturbations below 500 km depth. Images of the tomography model are shown in Figures 6 and 7, as well as with a discrete color scale in Figures S11 and S12 (^{footnote 1}). Several prominent velocity anomalies stand out in the upper mantle: a fast velocity anomaly parallel to the Aegean trench (labeled Aegean slab, or AS, in figures), a fast velocity anomaly north of the Cyprean trench (labeled WC and EC in figures), a slow velocity anomaly (labeled East Anatolian slow, or EAS, in figures) beneath the high elevations of eastern Anatolia, and a fast velocity anomaly (labeled as İstanbul zone, or IZ, in figures) underlying the İstanbul zone and the Pontide block in northern Anatolia. These prominent anomalies give way to interconnected low-amplitude fast velocity anomalies in the lower mantle across the whole region. To more clearly illustrate the three-dimensional (3-D) geometries of each of these anomalies, we present further descriptions of the model in three separate convergent domains, here listed from west to east: the Aegean domain, the Cyprean domain, and the Bitlis-Zagros domain.

In the Aegean domain, the most prominent anomaly is the fast velocity anomaly coincident with the Aegean Wadati-Benioff zone, anomaly AS (Figs. 6 and 7). However, while Wadati-Benioff zone earthquakes terminate at ∼200 km depth (Papazachos et al., 2000), anomaly AS extends to the base of the model domain at 900 km depth, with a regular dip of ∼60° toward the northeast (Fig. 7A). Anomaly AS has a reduced amplitude (average velocity perturbation ∼1%) between 300 and 400 km depth before widening in the transition zone and uppermost lower mantle.

The Cyprean domain is dominated by a fast velocity anomaly to the north of the Cyprean trench that is distinct from the fast velocity anomaly in the Aegean domain, separated by vertically extensive slow velocities (labeled STEP tear in Fig. 6B). We refer to this as the Cyprean slab anomaly (labeled as EC and WC in Figs. 6, 7B, and 7C). The character of the Cyprean slab anomaly varies from west to east. The western part of the Cyprean anomaly (anomaly WC) coincides with Wadati-Benioff zone earthquakes that extend to ∼150 km depth. Anomaly WC dips regularly to the northeast with high amplitude (average velocity perturbation ∼3%) in the upper mantle (Fig. 7B). Upon reaching the transition zone, the amplitude of the anomaly is significantly reduced to ∼2%, and the anomaly widens to the north by ∼200 km. Low-amplitude fast velocity anomalies in the uppermost lower mantle beneath the Cyprean slab appear to be connected to the anomaly in the transition zone. Conversely, the eastern part of the Cyprean fast velocity anomaly (anomaly EC) directly north of Cyprus extends northward to ∼37.5°N and in depth from the model surface to 200 km at a subhorizontal angle (Fig. 7C). Approximately 50 km to the north, there is a high-amplitude (average velocity perturbation ∼3%), subvertically dipping, fast velocity anomaly (anomaly EC2) that extends from 200 km to 700 km depth. In cross section, this block appears separated from anomaly EC by low-amplitude slow velocities. While the amplitude of the intervening slow velocities is low (<1%), synthetic recovery tests indicate that significant separation between EC and EC2 is necessary to allow for intervening slow velocities, suggesting the separation is robust (Fig. S8 [^{footnote 1}]). Within the mantle transition zone, this large, vertically dipping anomaly appears connected to a more disparate, low-amplitude fast velocity anomaly (anomaly EC3; average velocity perturbation ∼1.5%) extending to the north, as well as into the lower mantle. In the lower mantle, smaller fast velocity anomalies are scattered below the larger blocks. Notably, while the western and eastern portions of the Cyprean slab anomaly display distinct character, they appear to be connected.

The Bitlis-Zagros domain (Fig. 7D) is almost completely dominated by slow velocity perturbations (anomaly EAS), particularly at depths above ∼200 km (average velocity perturbation ∼–3%). One prominent exception is a small fast velocity anomaly (anomaly EAF), centered at 42°E, 40°N, that is ∼100 km in diameter and extends from ∼100 to 300 km depth (average velocity perturbation ∼2%). Below 200 km depth, the amplitudes of slow velocity perturbations are significantly lower (average velocity perturbation ∼–1%; Fig. 7D). Below ∼300 km, additional small, low-amplitude fast velocity anomalies (average velocity perturbation ∼1%) appear scattered across the domain. Interconnected low-amplitude fast velocity anomalies (average velocity perturbation ∼1.5%) are prevalent near the base of the model (Fig. 6G).

Across all three convergent domains at shallow depths, there is a prominent fast velocity anomaly (anomaly IZ) that extends below the İstanbul zone, the Pontide block, and the southern Black Sea (Fig. 6). Anomaly IZ (average velocity perturbation ∼3%) exhibits variable amplitude and shape, but it extends vertically to ∼200 km depth. Widespread slow velocity perturbations separate anomaly IZ from the shallow slab anomalies.

Before discussing the implications of the new features imaged with this model, we confirm the consistency of our model with other published seismic models in the literature by identifying expected velocity anomalies in the model domain that have previously been deemed robust. The most obvious expected anomalies are fast velocity anomalies associated with Wadati-Benioff zone earthquakes. Wadati-Benioff zones extend within our study region downdip from the Aegean Trench to ∼200 km depth and from the Cyprean trench to ∼100 km depth. As expected, earthquakes in both Wadati-Benioff zones fall within high-amplitude fast velocity anomalies in the tomography model (Fig. 6A) that we interpret as the Aegean and Cyprean slabs, respectively. From prior seismic tomography (Bijwaard et al., 1998; Piromallo and Morelli, 2003), the Aegean slab anomaly has been known to extend into the lower mantle to at least 1000 km depth (Faccenna et al., 2003), and the Cyprean slab has been seen to extend at least into the mantle transition zone, indicating that both slabs are aseismic at depth. Both observations are confirmed with our new model.

A second expected feature in the tomography model is the presence of widespread slow velocity anomalies throughout the uppermost mantle beneath Anatolia, with particularly slow velocities beneath eastern Anatolia (Gök et al., 2007; Gans et al., 2009; Biryol et al., 2011; Mutlu and Karabulut, 2011; Bakırcı et al., 2012; Salaün et al., 2012; Fichtner et al., 2013; Warren et al., 2013; Delph et al., 2015a; Govers and Fichtner, 2016; Salah, 2017). At the regional scale, shallow slow velocities are typically explained by mantle upwelling in response to slab detachment in eastern Anatolia (Keskin, 2003; Şengör et al., 2003), drip-style lithospheric delamination in central Anatolia (Göğüş et al., 2017), or rollback-style lithospheric delamination across central and eastern Anatolia (Bartol and Govers, 2014) and western Anatolia (van Hinsbergen et al., 2010). This manifests at the surface in a number of young volcanic fields that show an asthenospheric contribution to their geochemical signatures, such as in the widespread volcanism across the East Anatolian Plateau (Pearce et al., 1990; Şengör et al., 2008), mafic volcanism in the Central Anatolian Volcanic Province (Reid et al., 2017), the Miocene–Pliocene Kırka-Afyon-Isparta volcanic field (Francalanci et al., 2000), and the Pliocene–Pleistocene Kula volcanic field (Richardson-Bunbury, 1996; Savaşçın and Oyman, 1998; Tokçaer et al., 2005; Dilek and Altunkaynak, 2009). Our model confirms the presence of slow velocities at shallow depths across the region (Fig. 6A), with particularly slow velocities beneath eastern Anatolia (Fig. 4).

A third expected feature is the shallow fast velocity anomaly observed at the northern edge of the model (anomaly IZ; Figs. 6A and 7A). A similar anomaly was present in previous regional P-wave (Biryol et al., 2011; Yanovskaya et al., 2016) and S-wave (Salaün et al., 2012; Bakırcı et al., 2012; Fichtner et al., 2013) tomographic velocity models. This is typically interpreted as the strong mantle lithosphere of the Black Sea (Jiménez-Munt et al., 2003; Yegorova et al., 2013), the İstanbul zone near the western Pontides (Biryol et al., 2011), and the Sakarya zone (Fig. 1) near the central Pontides (Fichtner et al., 2013). Alternatively, correlation of the anomaly’s southern boundary with the Izmir-Ankara-Erzincan suture zone has led to the interpretation of the anomaly as a remnant of Neo-Tethys subduction (Pourteau et al., 2010; Salaün et al., 2012) and/or partly related to underthrusting of the crystalline Kırşehir block (Espurt et al., 2014). However, synthetic recovery tests (Fig. S7 [^{footnote 1}]) indicate that our model is unable to distinguish additional details of the structure of the anomaly, so we refrain from interpreting its origin.

As a whole, the present model accordingly corroborates each of the aforementioned studies, showing fast velocity anomalies corresponding to the locations of expected subducting slabs, prevalent slow velocities throughout most of the upper mantle beneath eastern Anatolia and recently active volcanism in central and western Anatolia (Fig. 7D), and fast seismic velocities at shallow depths north of the Izmir-Ankara-Erzincan suture zone.

Significant fast velocity anomalies are observed in the upper mantle with varying intensity and structure in each of the three convergent domains: the Aegean domain, the Cyprean domain, and the Bitlis-Zagros domain. Here, we present our interpretations of these fast velocity anomalies and their relationship to subducted African-Arabian lithosphere. Our interpretation of the subducted slab structure is summarized in Figure 8.

Fast velocity anomaly AS is coincident with Wadati-Benioff zone earthquakes at shallow depths and extends to the base of our model at 900 km depth (Figs. 6 and 7). The slab appears as a high-amplitude fast velocity anomaly that is continuous from the shallowest layers in our model to the lower mantle. Besides the concave curvature of the Aegean slab at shallow depths, which is related to the shape of the trench, the slab appears to have an otherwise simple character, with a consistent dip through the mantle. In contrast, the subducting Cyprean slab can be described in two parts: the western Cyprean slab (Fig. 7B) and the eastern Cyprean slab (Fig. 7C). Along the western Cyprean trench, where oceanic lithosphere is entering the trench (Granot, 2016), the slab exhibits a simple geometry similar to the Aegean slab. Here, the slab appears to subduct continuously through the upper mantle with a dip of ∼45° to the NE, coinciding with a clear Wadati-Benioff zone. However, along the eastern portion of the Cyprean trench, the complex geometry of the eastern slab may be a result of the underthrusting of attenuated continental lithosphere at the trench (Delph et al., 2017), such as the Eratosthenes Seamount (Robertson, 1998; Schattner, 2010). Not only did the impingement of the Eratosthenes Seamount on the trench contribute to recent uplift on Cyprus (Robertson, 1998), but the subduction of attenuated continental material along the Cyprean trench is well correlated with a region of relative aseismicity in the shallow slab and northward-directed subhorizontal subduction (anomaly EC) between longitudes 32°E and 34°E (Fig. 7C). While relatively poor resolution near the surface hinders our ability to image the shallow slab segment in detail, this result is consistent with previous seismic imaging studies sensitive to uppermost mantle velocity structure (Bakırcı et al., 2012; Abgarmi et al., 2017; Delph et al., 2017).

Our interpretation of the eastern Cyprean slab supports a buoyancy-driven slab deformation mechanism. The shallowly subducting slab segment is separated from a large vertically dipping block of the Cyprean slab (anomaly EC2) that appears to be connected to the western Cyprean slab in the mantle transition zone. If the shallow segment consists of subducted continental lithosphere while the separated block consists of the leading oceanic lithosphere, we would expect a discrepancy in buoyancy, leading to a difference in the dip of the slab segments. Separation of the two slab segments (Fig. 6C; Fig. S8 [^{footnote 1}]) further implies that this buoyancy contrast recently induced horizontal tearing of the slab. We infer the tearing to be propagating from east to west, where it ends near the edge of the subducting continental material, in a similar fashion to what is observed elsewhere in the Mediterranean system (Wortel and Spakman, 2000). Such a tear, coupled with slab rollback, likely alters mantle flow, allowing for upwelling through the tear and possibly contributing to the observed slow upper-mantle velocities and active volcanism in the Central Anatolian Volcanic Province above the vertical block. Additionally, detachment of the vertical block and potential rebound of the buoyant continental material have previously been interpreted to contribute to recent uplift of the Central Taurus Mountains (Schildgen et al., 2014; Abgarmi et al., 2017; Delph et al., 2017).

To the east of the Cyprean trench, there is the Bitlis-Zagros suture zone, which demarcates the former subduction margin of the once-expansive Neo-Tethys Ocean (Şengör and Kidd, 1979; Bozkurt, 2001; Okay et al., 2010). This segment bridges the gap between the active subduction zone at the Cyprean trench to the west and the Zagros collision front to the east. We observe the eastern edge of the Cyprean slab at longitude 35°E in our tomography model. The Zagros slab to the east has been well imaged with teleseismic tomography (Alinaghi et al., 2007), but poor resolution on the edge of that model prevented a clear demarcation of the western edge of the Zagros slab. Similarly, our model has poor resolution of the Zagros slab, but we do observe a fast velocity anomaly just outside our well-resolved model space (labeled as “Zagros slab?” in Fig. 6). This confirms that the Zagros slab fast anomaly is not continuous with the Cyprean slab anomaly across the Bitlis-Zagros suture zone. Instead, the upper mantle is dominated by slow velocity anomalies—a marked contrast to the Aegean and Cyprean domains, which each contain prominent fast velocity anomalies in the upper mantle. Since collision along the Bitlis-Zagros suture zone, the subducted Neo-Tethyan lithosphere in this region, which we refer to as the Bitlis slab, steepened and detached (Keskin, 2003; Şengör et al., 2003; Keskin, 2007) and likely sank into the mantle. Sinking of the Bitlis slab and the resultant upwelling of asthenosphere are often interpreted to be the source of slow upper-mantle velocities beneath eastern Anatolia (Lei and Zhao, 2007; Zor, 2008). Thus, our results are consistent with slab detachment beneath eastern Anatolia, but they raise the question of the current whereabouts of the detached slab.

One approach to predicting the current lateral location of the detached Bitlis slab is to assume steady plate kinematics since the time of detachment ca. 10 Ma (Şengör et al., 2003). If we assume that (1) Arabian plate velocity has been nearly constant since the early Miocene (ArRajehi et al., 2010) with absolute velocities of ∼35 mm/yr WNW, (2) the slab detached at the reconstructed location of the Bitlis-Zagros suture zone, and (3) the slab sank vertically since detachment, then we predict that the detached Bitlis slab should be roughly 350 km to the ESE of the Bitlis-Zagros suture zone, in the same location as the imaged Zagros slab (Alinaghi et al., 2007). This would imply that the poorly resolved fast anomaly we see in the southeast portion of our model could be the Bitlis slab rather than the Zagros slab. However, this interpretation does not allow for the observed fast velocity anomalies to comprise the total volume of subducted lithosphere, which we would expect to span the entire margin. If the poorly resolved fast velocity anomaly is in fact the Bitlis slab, then there is no space in the mantle to accommodate the subduction of the Zagros slab to the east, suggesting the anomaly is unlikely to be the Bitlis slab. This further implies that a cohesive Bitlis slab is not present at these depths (∼350 km and below), or it has migrated laterally outside of our model domain, presumably due to ambient mantle flow.

Another possible anomaly that could account for the Bitlis slab is anomaly EAF (Figs. 6A and 6B). Previous authors have interpreted this anomaly to represent the detached Bitlis slab (Lei and Zhao, 2007; Bakırcı et al., 2012). However, the observed anomaly accounts for only a small portion of the expected subducted slab volume. If it represents the detached Bitlis slab, its presence at 100 km also implies the body was dense enough to detach from the Arabian plate at the surface, but it must have quickly become neutrally buoyant, allowing it to remain at shallow depths from the mid-Miocene to the present.

Instead of these two possibilities, neither of which can account for the assumed once laterally continuous subducted slab, we consider it more likely that the Bitlis slab did not detach as a large, intact body, but rather the Bitlis slab fragmented as it detached. Around 11–10 Ma, collision and tearing of the leading oceanic lithosphere allowed for the opening of a slab window and subsequent asthenospheric upwelling (Keskin, 2003; Şengör et al., 2003). Following this initial opening, further tearing and deformation caused the subducted slab to fragment, leading to the dispersed, small, and low-amplitude fast velocity anomalies seen in our images (Fig. 7D). These anomalies, and thus the fragmented Bitlis slab, are located beneath eastern Anatolia between the Cyprean and Zagros slabs, indicating that detachment, deformation, and fragmentation are not instantaneous, but rather represent a multistage process (Wortel and Spakman, 2000; Govers and Wortel, 2005). This in turn prevents the detaching Bitlis slab from decoupling from the surrounding Cyprean and Zagros slabs and sinking independently of plate motion. The Bitlis slab fragments have not previously been imaged as seen in this study. However, analyses of conversion amplitudes from the 410 and 660 km discontinuities and transition zone thickness beneath eastern Anatolia using the receiver function technique suggest that cold temperatures influence phase conversions in this region (Ozacar et al., 2008). Ozacar et al. (2008) interpreted this as the presence of fragments of delaminated Eurasian lithosphere, rather than African-derived lithosphere. As noted in that study, however, the presence of cold slab fragments is equally consistent with the observed variations in transition zone thickness.

Given the wide range of upper-mantle slab sinking rates in the literature (Billen, 2010; van der Meer et al., 2010; Butterworth et al., 2014; Ballmer et al., 2015), it is difficult to determine if the depth of the Bitlis slab fragments matches the inferred mid-Miocene timing of slab break-off. One common way to investigate the appropriate depth is to assume that slab sinking speed is equal to convergence rate (Schellart et al., 2009). With an Arabia-Eurasia convergence rate of 15 mm/yr, this amounts to a slab depth of ∼150 km after a 10 Ma detachment. However, detachment and fragmentation of the Bitlis slab make this a minimum estimate because free sinking would likely be faster than the relatively slow Arabia-Eurasia convergence rate. The high end of estimated sinking rates (up to 150 mm/yr) would place the slab well into the lower mantle and possibly below our model domain, even when considering a sinking deceleration in the lower mantle (Lithgow-Bertelloni and Richards, 1998). While this would provide an alternative explanation for the lack of a continuous fast velocity anomaly, these rates are considered outliers, as sinking rates >60 mm/yr only apply to continuous slabs that have stalled in the mantle transition zone (Billen, 2010). Thus, having the shallowest Bitlis fragments visible within our model domain in the transition zone and uppermost lower mantle is consistent with our interpretation.

Interpreting the disparate fast velocity anomalies beneath eastern Anatolia as slab fragments requires that the small, low-amplitude fast velocity anomalies we image are robust. While it is difficult to prove that such low-amplitude anomalies are not artifacts of the tomographic inversion, we show that anomalies of that magnitude beneath eastern Anatolia are unlikely to appear without true fast velocities. This is made clear with synthetic recovery tests of slab fragments with varying input amplitude (Fig. S9 [^{footnote 1}]). In such tests, only true velocity anomalies with high amplitude (at least 3%) yield recovered amplitudes comparable to the observed velocities. Last, synthetic tests with induced data noise (Rawlinson et al., 2014) indicate that unreasonably high amounts of noise (>100% of model misfit) are required to produce anomalies as strong as the interpreted fragments beneath eastern Anatolia (Fig. S10 [^{footnote 1}]). These tests are described in more detail in the Supplemental Material (see ^{footnote 1}).

One key feature of the Biryol et al. (2011) model is that each of the subducted slabs beneath Anatolia appears to flatten in the mantle transition zone, at least temporarily. However, this may be an artifact resulting from a model domain that does not extend into the lower mantle. As we find in our model, significant amplitude is artificially placed in the deepest layers of the model during inversion (Fig. S4 [^{footnote 1}]). It is likely that a similar effect occurred during inversion for the model of Biryol et al. (2011), such that artificially fast velocity anomalies were seen to be continuous and flat-lying in the mantle transition zone.

Our new model differs from the model of Biryol et al. (2011) in that the model domain extends into the lower mantle, allowing for the assessment of the subducting African lithosphere’s interactions with the 660 km discontinuity. Seismic tomography (Fukao et al., 2001; Scire et al., 2017) and geodynamic modeling (Čížková et al., 2007; Billen, 2010; Stegman et al., 2010; Garel et al., 2014) studies indicate that the viscosity jump observed at the base of the mantle transition zone (Quinteros et al., 2010; Rudolph et al., 2015) may serve to inhibit slab descent through stagnation and/or deformation. The way in which the subducting plate interacts with the mantle transition zone has broad implications for processes such as mantle convection (van der Hilst et al., 1997), overriding plate volcanism (Richard and Iwamori, 2010), and possibly for orogenic plateau development (Faccenna et al., 2013).

In the Mediterranean region, southward trench migration from slab rollback is generally thought to cause stagnation in the mantle transition zone as the 660 km discontinuity serves to block slab descent into the lower mantle (Wortel and Spakman, 2000; Piromallo and Morelli, 2003). This constrains mantle convection to the upper mantle, which may have a significant impact on surface kinematics and deformation (Faccenna et al., 2013, 2014). Previous seismic imaging efforts differ on how slabs beneath Anatolia relate to this model. While the regional tomography model of Biryol et al. (2011) shows the subducted lithosphere lying flat in the mantle transition zone, global tomography models show that the subducted lithosphere beneath Anatolia may temporarily stagnate, but eventually sinks into the lower mantle (Bijwaard et al., 1998; Piromallo and Morelli, 2003). In our new model, fast seismic velocity perturbations do not dominate the mantle transition zone (Figs. 5, 6, and 7), appearing more disjointed, but they nonetheless continue into the lower mantle throughout the study area (Figs. 6E and 6F).

In the Aegean domain, the Aegean slab anomaly is continuous from the surface into the lower mantle, appearing to traverse the transition zone relatively uninhibited (Fig. 7A). It is possible that the thickness of the Aegean slab fast velocity anomaly in the transition zone (∼300 km) indicates deformation associated with slab penetration into the lower mantle (Stegman et al., 2010; Garel et al., 2014). However, due to a number of factors, including the relatively constant thickness of the slab anomaly and continuous penetration into the lower mantle, we consider the slab thickness to more likely result from a combination of tomographic smearing in the Aegean region and the 3-D concave structure of the Aegean slab. In the Cyprean domain, the vertical limb of the Cyprean slab anomaly cohesively extends to the base of the transition zone without penetrating to the lower mantle (anomaly EC2; Fig. 7C). However, north of the vertical block, lower-amplitude fast velocity perturbations extend from the transition zone into the lower mantle (anomaly EC3). While the vertical block itself appears continuous in its descent, one interpretation is that, contrary to the Aegean slab, the vertical Cyprean slab block buckles upon interacting with the 660 km discontinuity. Weaker fast velocities in front of the separated block may be fragments of older Cyprean slab that have broken off, piled up, and eventually sunk into the lower mantle. Thus, in central Anatolia, the 660 km discontinuity induces slab deformation but allows for slab penetration to the lower mantle in lieu of stagnation. In the Bitlis-Zagros domain (Fig. 7D), disparate fast velocity bodies are scattered within the transition zone and uppermost lower mantle. We similarly interpret these fast velocity anomalies as fragments of the detached Bitlis slab, though it is unclear whether fragmentation was associated with the slab’s interaction with the 660 km discontinuity or took place at shallower depths. As with the Cyprean slab, fast velocity perturbations in the uppermost lower mantle indicate that the 660 km discontinuity does not prevent the descent of the deformed Bitlis slab fragments into the lower mantle. Across each domain, contrary to the greater Mediterranean system, the viscosity boundary associated with the transition from the upper to lower mantle does not appear to cause long-term stagnation of Eastern Mediterranean slabs in the mantle transition zone, but instead likely induces at least some slab deformation.

The traditional view of subduction is one in which denser lithosphere descends into the upper mantle, possibly stalling in the transition zone due to a viscosity contrast at the 660 km discontinuity (Fukao et al., 2001; Goes et al., 2017), before coming to rest near the core-mantle boundary (Wysession, 1996; Rao and Kumar, 2014). However, evidence is mounting that many slabs do not follow such a simple path (Nolet, 2009). In particular, Mediterranean slabs are seen to significantly deform prior to reaching the 660 km discontinuity (von Blanckenburg and Davies, 1995; Wortel and Spakman, 2000; Biryol et al., 2011) and often do not readily continue their descent into the lower mantle (Piromallo and Morelli, 2003). Eventually, continental impingement on a trench may lead to detachment of the slab and a transition from subduction to collision, as inferred from seismic tomography in the India-Asia collision (Van der Voo et al., 1999b; Koulakov and Sobolev, 2006).

Anatolia marks a major transition from primarily oceanic subduction in the west to continental collision in the east, separated by a transitional segment in central Anatolia where attenuated continental lithosphere is subducting. In Anatolia, we suggest that we can relate the variable morphology of subducting lithosphere across three distinct convergent domains with the temporal evolution of subduction termination. Thus, Figure 8 can be viewed not only as an interpretation of the modern mantle structure beneath Anatolia, but additionally as a depiction of the temporal evolution of the subduction termination process when viewed from west to east. With this model, the Aegean domain represents the continuous subduction of oceanic lithosphere beneath a continental upper plate, which occurs throughout the majority of the lifetime of a subduction zone. In this domain, we observe a Wadati-Benioff zone down to ∼200 km depth and a continuous, relatively undeformed slab extending from the surface down into the lower mantle. East of the Aegean domain, the Cyprean domain represents subduction actively transitioning to termination. Here, slab character is variable, in that the western Cyprean slab appears continuous with a consistent dip and a Wadati-Benioff zone extending to ∼150 km depth, while the eastern Cyprean slab exhibits signs of severe deformation and fragmentation, with a horizontal tear separating slab segments of distinct dip. The variation in subduction character within the Cyprean domain is likely determined by the character of material entering the trench on the downgoing plate, with oceanic lithosphere subducting at the western Cyprean trench and thinned continental lithosphere further to the east (Figs. 1 and 8). The severe, shallow slab deformation associated with the impingement of continental lithosphere on the trench is evidence in the mantle that the Cyprean domain is both a spatial transition from the Aegean domain and a transition in subduction stage. Finally, the Bitlis-Zagros domain east of the Cyprean domain represents the final stage in the transition from subduction to the collision of two continental plates and the complete removal of the subducting lithosphere. Subduction is thought to have ended ca. 10 Ma in this region following collision and slab detachment (Keskin, 2003; Şengör et al., 2003), such that the most recently subducted slab has likely descended to at least the mantle transition zone. What remains is a collection of torn, deformed fragments of Bitlis slab with a weakened seismic signature. This evolutionary model predicts that the mantle character of the eastern Cyprean domain is approximately what subduction in the Bitlis-Zagros domain may have looked like ca. 10 Ma, and what the Aegean domain will look like when the African continental margin begins to enter the Aegean Trench. This model additionally has significant implications for the fate of subducted material. While some slabs may be able to cohesively traverse the mantle to the core-mantle boundary (Wysession, 1996), our images suggest that in other cases, slab destruction and assimilation are catalyzed at much shallower depths by tearing and fragmentation upon continental collision.

We used finite-frequency teleseismic P-wave tomography to obtain an improved seismic velocity model of the Eastern Mediterranean mantle. Our new P-wave tomography model reveals the complex, segmented structure of subducting African-Arabian lithosphere beneath Anatolia and its interactions with the surrounding mantle and the mantle transition zone in three separate convergent domains. Three main conclusions can be drawn across convergent domains:

• (1) Our model confirms the work of several well-established features of Anatolian mantle velocity structure shown in prior seismic studies, including: (a) a continuous fast velocity anomaly associated with the Aegean slab extending from the surface into the lower mantle, (b) widespread slow velocities in the shallow mantle beneath Anatolia, with particularly slow velocities throughout the upper mantle beneath eastern Anatolia, and (c) fast seismic velocities in the shallow mantle north of the Izmir-Ankara-Erzincan suture zone.

• (2) Slabs in the Eastern Mediterranean interact with the 660 km discontinuity differently than is generally observed in the greater Mediterranean system. Rather than flattening in the mantle transition zone, slabs beneath Anatolia all appear to descend into the lower mantle, with the Aegean slab crossing the discontinuity unimpeded, the Cyprean slab continuing into the lower mantle with signs of significant deformation in the mantle transition zone, and a significantly fragmented Bitlis slab appearing both above and below the 660 km discontinuity.

• (3) The distinct character of subduction within each convergent domain is correlated with the material entering the trench on the downgoing plate, varying from oceanic lithosphere subducting at the Aegean Trench to continental lithosphere colliding with eastern Anatolia at the Bitlis-Zagros suture zone. This allows us to interpret spatial variations in slab structure as different temporal stages of subduction termination. The Aegean slab is undergoing little to no deformation as it dips continuously into the lower mantle, implying it is in its main subduction stage. The Cyprean slab, which is severely deformed through tearing and fragmentation, is the middle stage between those of the Aegean and Bitlis slabs, when continental lithosphere is beginning to underthrust at the trench, likely blocking subduction, inducing tearing, and initiating the breakup of the downgoing Cyprean slab. Thus, the Cyprean slab exhibits subduction termination in progress. The Bitlis slab is completely fragmented within and below the mantle transition zone and shows no continuity with the surface plate, indicating a complete transition from subduction to continental collision in the Bitlis-Zagros domain. This model has significant implications for our view of a slab’s life cycle, suggesting that some slabs may significantly fragment at the end of subduction rather than cohesively sink through the mantle after detachment.

We would like to thank the extraordinary efforts made by the entire Continental Dynamics–Central Anatolian Tectonics (CD-CAT) team, both in collecting our critical seismic data set and in invaluable workshops that have made this article possible. We would also like to thank Brandon Schmandt and Alissa Scire for thoughtful consultations on the finite-frequency tomography method. Last, we would like to thank two anonymous reviewers for thorough, thoughtful revisions that helped to greatly improve this manuscript. This research was supported by National Science Foundation grant EAR-1109336 (CD-CAT). Delph was supported by the Wiess Postdoctoral Research Fellowship at Rice University. The tomography model from this study will be made available online through the IRIS Earth Model Collaboration (http://ds.iris.edu/ds/products/emc/) with the label ANA2_P_2018.