The effects of subduction termination on the continental lithosphere: Linking volcanism, deformation, surface uplift, and slab tearing in central Anatolia

The tectonic implications of ocean-basin closure are ill defined, despite being a crucial part of the Wilson Cycle in plate tectonics (Wilson, 1968). This process is currently taking place in the Mediterranean region, where convergence between the African and Eurasian plates is leading to the termination of subduction and final closure of the long-lived Tethyan Ocean as part of the ∼15,000-km-long Alpine-Himalayan orogenic belt (Lister et al., 2001). The Anatolian region is a complex zone of deformation that is centrally located in this extensive orogenic belt. In the east, complete oceanic closure was accomplished at the Arabia-Eurasia collision zone in the early Miocene (Okay et al., 2010), while continuous subduction of old Neotethyan lithosphere is still taking place along the Aegean trench in the west (Granot, 2016; Fig. 1). Central Anatolia represents a transition between these end-member processes, making it an ideal place to study the tectonic interactions between the overriding and subducting plates during the terminal stages of subduction.

Deformation within the Anatolian plate is intimately tied to the behavior of subducting African lithosphere along its southern margin. In the west, slab rollback along the Aegean trench is imposing some of the highest observed extensional strain rates in the world (Le Pichon and Kreemer, 2010) above a seemingly coherent subducting plate (Spakman et al., 1993; Piromallo and Morelli, 2003; Biryol et al., 2011), while at the Arabia-Eurasia collision zone in the east, evidence for an underthrusting slab is largely absent below the third largest orogenic plateau on Earth (the East Anatolian Plateau, Keskin, 2003; Şengör et al., 2003). Between these regions lies central Anatolia, where the character of the subducting slab is enigmatic due to a lack of a laterally continuous Wadati-Benioff zone and the fragmented appearance of the slab from tomographic images (Biryol et al., 2011). This subduction is further complicated by the arrival of attenuated Gondwanan continental lithosphere at the Cyprean trench (Schattner, 2010; Granot, 2016). These observations indicate that a recent transition from oceanic subduction to passive margin underthrusting has occurred in central Anatolia.

To better understand the relationship between the tectonics of central Anatolia and the character of subduction along its southern margin, a dense array of 64 broadband seismic instruments was deployed at 72 locations as part of the Continental Dynamics–Central Anatolian Tectonics (CD-CAT) project. These stations were operational for a two-year period (May 2013–May 2015) and span the boundaries between the Anatolian, Eurasian, and westernmost Arabian plates, providing unprecedented spatial sampling throughout the region. In this study, we jointly invert P-wave receiver functions (Abgarmi et al., 2017) with Rayleigh-wave dispersion data derived from ambient noise cross correlations and earthquakes to obtain a high-resolution model of the shear-wave velocity structure beneath central Anatolia down to ∼150 km. These new seismic images will provide information about the fine-scale interactions between the underthrusting African lithosphere and the overriding lithosphere of the Anatolian plate.

The Anatolian region comprises various continental fragments that have undergone a complex history of amalgamation related to the closure of the Tethyan Ocean system (Şengör and Yilmaz, 1981; Bozkurt, 2001; Whitney and Hamilton, 2004). Collision of the Arabian plate in the east and slab rollback along the Aegean trench in the west have led to the formation of the conjugate strike-slip North and East Anatolian fault zones, which exemplify “tectonic escape” (Fig. 1) (Şengör et al., 1985; Burke and Şengör, 1986). These plate-bounding fault zones allow the Anatolian plate to translate westward and rotate counterclockwise (with respect to fixed Eurasia) as a coherent unit in a regional sense (Reilinger et al., 2006). The apparent coherency of this tectonic block, however, does not imply simplicity in the tectonic structures that have been imposed by the long-lived evolution of the Tethyan Ocean and the recent dynamics of the convergent margin.

The Late Cretaceous–Early Tertiary evolution of central Anatolia can be summarized by the closure of the northern branch of the Neotethys Ocean along the south-vergent Late Cretaceous–Paleocene Izmir-Ankara-Erzincan suture zone (IAESZ) and Inner-Tauride suture (ITS) (Şengör and Yilmaz, 1981; Boztuǧ et al., 2009; Parlak et al., 2013). The accretion of continental material (namely, the Kirşehir and Anatolide-Tauride blocks; Fig. 2 inset) onto the southern margin of Eurasia as crustal-scale nappe complexes was accommodated by detachment and decoupling of the crust from the underlying mantle lithosphere (Bartol and Govers, 2014; van Hinsbergen et al., 2016). This has led to a continuous southward migration of the trench from the Eocene to present (Menant et al., 2016; van Hinsbergen et al., 2016). However, the mantle lithosphere of the subducting plate may have continued to underthrust at a relatively shallow angle during the Eocene–Oligocene (Bartol and Govers, 2014). Shortly thereafter, the subducting slab began to founder from the overriding plate, leading to extension in the upper plate as evidenced by the early-mid Miocene unroofing of metamorphic core complexes in the Kirşehir Block (Fayon and Whitney, 2007). Foundering was likely caused by a decrease in the convergence rate between Africa and Eurasia, as indicated by the development of backarc basins and trench migration throughout the entire Mediterranean region (Jolivet and Faccenna, 2000). The presumably thin lithosphere of the Central Anatolian region was then exposed to hot upwelling asthenospheric material that adiabatically melted and initiated large ignimbritic eruptions beginning in the mid-Miocene (Le Pennec et al., 1994; Aydar et al., 2012). The subducting slab appears to have foundered in the southwestward direction in central Anatolia, as indicated by the southwest-younging of the initiation of volcanism in the Central Anatolian Volcanic Province (Schleiffarth et al., 2015) and NE-SW–oriented shear-wave splitting orientations in Anatolia (Paul et al., 2014 and references therein). In the late Miocene, the initiation of the uplift of nearly flat-lying carbonates in the Central Taurus Mountains is interpreted to result from slab fragmentation and subsequent rebound in the downgoing plate (Duretz et al., 2011; Schildgen et al., 2014). In Pliocene to younger times, subduction has been further complicated by the arrival of Gondwanan-derived continental fragments along the southern margin of Anatolia (e.g., Schattner, 2010; Delph et al., 2015a). This underthrusting of attenuated continental material has led to an acceleration in the uplift rate of the Central Taurus Mountains (Schildgen et al., 2012) and an apparent change in dip of the Cyprean slab from W to E (Biryol et al., 2011); this change in dip corresponds well with the inferred boundary between subducting oceanic and continental material in the Eastern Mediterranean (Granot, 2016).

Two main geologic terranes comprise our study area: the Kirşehir Block and the Anatolide-Tauride Block (Fig. 2). The basement of the Kirşehir Block consists primarily of crystalline metamorphic and plutonic rocks of island-arc affinity as well as obducted ophiolitic fragments (Whitney et al., 2001; Lefebvre et al., 2015). This domain is bounded in the north and south by the Izmir-Ankara-Erzincan suture zone (IAESZ) and Inner-Tauride suture (ITS) (Okay and Tüysüz, 1999; Parlak et al., 2013). The Inner-Tauride suture is thought to mark the closure of a small Neotethyan ocean basin that separated the Anatolide-Tauride Block from the Kirşehir Block (van Hinsbergen et al., 2016) and correlates with the surficial location of a drastic change in geophysical properties, including Bouguer gravity anomalies (Ates et al., 1999), Pn velocity (Gans et al., 2009; Mutlu and Karabulut, 2011), phase velocities from ambient noise tomography (Warren et al., 2013; Delph et al., 2015a), and crustal thickness (Abgarmi et al., 2017). The Anatolide-Tauride Block originated on the Gondwanan passive margin. In central Anatolia, it generally consists of unmetamorphosed marine sedimentary rocks where exposed (Cosentino et al., 2012) and is separated from the Kirşehir Block by the sedimentary basins of the Tuz Gölü, Konya, and Ulukişla basins (hereafter referred to as the central Anatolia Basin based to their similarities in velocity structure; however, this grouping does not imply that the evolution of these basins is related). The Anatolide-Tauride Block becomes more heavily deformed in the Eastern Taurus Mountains across the Central Anatolian fault zone due to stresses imposed by the northward motion of the Arabian plate.

The ∼730 km Central Anatolian fault zone is a seismically active, transtensional structure dominated by low-magnitude events and low rates of lateral displacement (no more than 80 km total and only ∼2–3 km since 13 Ma; Koçyiğit and Beyhan, 1998) that demarcates the transition between transpressional deformation to the east and transtension to the west and largely exploits the lithospheric weakness associated with the eastern Inner-Tauride suture. Deformation in the central portion of the Central Anatolian fault zone (Ecemiş corridor) has become largely extensional since the late Pliocene, leading to the formation of a pull-apart basin that is home to the northern portion of the Central Anatolian Volcanic Province (Toprak, 1998; Higgins et al., 2015).

The Central Anatolian Volcanic Province is a NE-SW–trending, long (∼300 km) and broad (40–100 km) volcanic region comprising voluminous felsic-intermediate ignimbrite deposits (>1000 km³; Le Pennec et al., 1994), large polygenetic stratovolcanoes, and smaller monogenetic cones and domes that correlate with inherited deformational structures near the Central Anatolian fault zone and Inner-Tauride suture (Toprak and Göncöoğlu, 1993; Toprak, 1998; Aydin et al., 2014). Ignimbrite volcanism in the region began in the mid-late Miocene (Aydar et al., 2012) and is concentrated near the projected junction of the Tuz Gölü fault and Central Anatolian fault zone at the buried transition from the Kirşehir Block to the Anatolide-Tauride Block. This large-volume ignimbrite volcanism gives way to smaller, more mafic eruptions to the south in the Anatolide-Tauride Block. The initiation of volcanism in the Central Anatolian Volcanic Province also appears to show a southwestward-younging trend (Schleiffarth et al., 2015).

Bounding the southern edge of the Central Anatolian Volcanic Province and Central Anatolian Basin are the Central Taurus Mountains. The Central Taurus Mountains have experienced multiple phases of uplift since the late Miocene: the first ca. 8 Ma, interpreted as a result of slab break-off, and the second beginning ca. 1.6 Ma, contemporaneous to the arrival of the Eratosthenes seamount at the Cyprean trench (Cosentino et al., 2012; Schildgen et al., 2012). These processes are interpreted to have caused a cumulative uplift of ∼2 km in the absence of significant compressional deformation as evidenced by nearly flat-lying Miocene shallow-marine carbonates that lie atop much of the Central Taurus Mountains. Upwelling asthenosphere as a result of slab break-off has been proposed to be responsible for this dramatic uplift (Cosentino et al., 2012); however, recent geophysical imaging of upper-mantle structure in central Anatolia indicates that the Cyprean slab may be present as far north as the northern margin of the Central Taurus Mountains (Bakırcı et al., 2012; Abgarmi et al., 2017).

The Bitlis-Zagros suture formed as a result of the final closure of the Neotethyan Ocean between Arabia and Eurasia in the early Miocene (Okay et al., 2010). Today, it represents the boundary between the Anatolian and Arabian plates in the west and Eurasian and Arabian plates in the east. The western portion of the Bitlis-Zagros suture roughly delineates the location of the long (∼500 km), broad (∼20 km), and very seismically active East Anatolian fault zone (EAFZ) (Bulut et al., 2012). In the Adana Basin, the EAFZ branches into multiple strands and transitions from dominantly left-lateral strike-slip motion to transtension, marking a diffuse plate boundary resulting from a compatibility gap at a continental triple junction (Dewey et al., 1986). This has led to drastic crustal thinning beneath the Adana Basin, a steep gradient in crustal thicknesses (Abgarmi et al., 2017), and differential uplift rates across the Kozan fault, which separates the Adana Basin from the Central Taurus Mountains (Radeff et al., 2015).

In this study, we add an additional two years of seismic data (May 2013–May 2015) from ∼170 stations from the recent CD-CAT seismic deployment and the Kandilli Observatory and Earthquake Research Institute (KOERI) to the data set of Delph et al. (2015a) (Fig. 1). The seismic portion of the CD-CAT project involved the deployment of 64 broadband seismic stations (STS-2 and CMG-3T sensors) in 72 locations throughout central Anatolia.

We perform ambient noise tomography on this new data set following the method of Bensen et al. (2007). The addition of the CD-CAT data set nearly doubles the number of raypaths that sample the Anatolian region from Delph et al. (2015a) while adding significant data coverage to the central portion of Anatolia. We first remove the trend, mean, and instrument responses from all vertical-component data and filter between the periods of 5–150 seconds so that consistent measurements of ground motion can be cross correlated for all contemporaneously recording seismic stations. These interstation cross correlations are then stacked over the duration of operation to increase their signal-to-noise ratio (SNR), and the positive and negative time lags of the cross correlations are averaged to obtain waveforms that represent surface (Rayleigh) waves traveling between stations. Waveforms with a SNR ≥ 10 are analyzed using frequency-time analysis (FTAN) to estimate interstation Rayleigh-wave phase velocities (Bensen et al., 2007). Recent research has shown that accurate Rayleigh-wave measurements can be obtained at a station spacing ∼1 wavelength, contrary to the assumption that an interstation spacing of at least 3× wavelength was necessary to accurately recover surface-wave velocities (Luo et al., 2015). This significantly increases the number of useable interstation phase velocities recorded by the CD-CAT array. These interstation phase velocities are then inverted for phase-velocity maps as a function of period (from 8 to 50 seconds) using the tomographic method of Barmin et al. (2001) with damping and smoothing parameters of 300 and 50 km. These parameters are chosen based on a trade-off plot of the root-mean-square (RMS) travel-time residual and the trace of the smoothing term in the inversion matrix normalized by the size of the smoothing matrix (Fig. 3A). The smoothing parameter is slightly smaller than that used in Delph et al. (2015a) as justified by the trade-off plot and the very dense ray coverage after the introduction of the CD-CAT data (Fig. 3A). This inversion also incorporates a ray-density–dependent weighting parameter that controls velocity variations in areas of poor ray density. Variation of this parameter has a relatively small effect on the resulting phase-velocity maps when applied to this ambient noise data set due to the very dense sampling of our study area, and we use a value of 100 for this study.

Two plane-wave tomography (TPWT; Forsyth and Li, 2005) has proven to be a powerful way to obtain surface-wave dispersion information with a distributed array. TPWT models an incoming surface wave as two interfering plane waves described by amplitude, phase, and direction to account for lateral velocity heterogeneities and scattering outside of the array, which can cause the incoming ray path to deviate from the great-circle path between the earthquake and array. To fully describe the different lateral sensitivities of the Rayleigh waves used in this study, finite frequency sensitivity kernels using the Born single-scattering approximation are used; these kernels allow for the accurate recovery of anomalies as small as ∼1 wavelength at the period of interest (Zhou et al., 2004; Yang and Forsyth, 2006). Since this research focuses on central Anatolia, we perform TPWT only on the CD-CAT stations.

For this data set, we inspect M_w ≥ 6 earthquakes at 25°–125° epicentral distances from our station locations on vertical 1 sample/second waveform data. First, we removed the trend, mean, and instrument responses of all stations to ensure data are in a common unit (displacement). After discarding all traces where an earthquake was not easily observable, we perform a narrow bandpass filter with a 10 mHz filter width centered at 13 frequencies with periods corresponding to 25–143 seconds. High-quality Rayleigh-wave arrivals were then windowed and used for TPWT. Periods <40 seconds returned an insufficient number of high-quality Rayleigh waveforms and were omitted from further analysis, as most of the higher frequency content of the Rayleigh waves had attenuated due to the large epicentral distances of the earthquakes in this study (>70°) (Fig. 4). Some longer periods were also excluded (125 and 143 seconds) because the wavelength of the Rayleigh waves at these periods surpassed the spatial dimensions of the seismic array, compromising the accuracy of the phase-velocity measurements (Fig. 5) (Yang and Forsyth, 2006). However, the 111 second phase-velocity measurement was retained because it could still recover the pattern of anomalies, albeit with underestimated amplitudes (Supplemental Fig. S1¹). The TPWT model region extends beyond our array so that nodes outside of the array can absorb velocity heterogeneities from scattering outside of the footprint of the seismic array.

Because the two plane-wave inversion is a smoothed, damped least-squares inversion, it is important to have an initial phase-velocity model that is close to the real phase-velocity structure of the study region (Rau and Forsyth, 2011). To account for vertical and lateral variations in crustal-velocity structure, we follow the method described by Delph et al. (2015a) to obtain a simple 3D shear-velocity model from ambient noise tomography (ANT)–derived dispersion data using an updated compilation of crustal thickness estimates for the Anatolian region (Supplemental Fig. S2 [see ^{footnote 1}]) (Özacar et al., 2010; Vanacore et al., 2013; Abgarmi et al., 2017). Below the crust, we append a constant-velocity mantle of Vs = 4.3 km/s based on the average shear velocity in the upper mantle beneath Anatolia (Fichtner et al., 2013; Delph et al., 2015a). Dispersion data from this 3D shear-velocity model was then calculated to produce a starting model of Rayleigh-wave phase velocities at different periods (40, 45, 50, 59, 67, 77, 91, 100, and 111 seconds) for the two plane-wave inversions. Other mantle velocities were tested as well, but these did not change the inverted phase-velocity maps by more than ∼50 m/s.

The two plane-wave inversion is dependent on a spatial smoothing and damping parameter to stabilize anomalies within the model area. For spatial smoothing, we use a similar Gaussian pulse width as used in the ambient noise tomographic inversion (50 km), which acts to smooth the sensitivity kernels estimated from the Born approximation. This smoothing parameter has a much greater effect on short-period kernels than those of longer periods (see Li, 2011, fig. 3, for example). We choose a damping parameter of 0.2 because the rank of the resolution matrix stabilizes at this value (Fig. 3B). Our grid spacing is set at 0.2° × 0.2° to match that of the ambient noise tomography grid.

This inversion is a two-step process for which we perform three iterations at each step. In the first step, the inversion solves for the amplitude, phase, and direction of two interfering plane waves that best reproduce the windowed waveforms from each event using a simulated annealing approach. In the second step, phase velocities are estimated using the observed and predicted arrival times of the waveforms based on the best-fitting wave parameters from the first step, which can be further adjusted in this step to better fit the waveforms (Forsyth and Li, 2005; Yang and Forsyth, 2006). With the data and inversion parameters used in this study, we were able to resolve anomalies of 1° × 1° at all periods investigated based on checkerboard tests (Supplemental Fig. S1 [see ^{footnote 1}]). Due to the increasing wavelength as a function of period, shorter period measurements are able to resolve smaller-scale structures enabling accurate recoveries of lateral gradients in phase-velocity structure, while longer periods (>91 seconds) accurately recover the pattern of the anomalies smaller than 1° × 1° but with diminished amplitudes.

By combining receiver functions derived from common conversion point (CCP) stacking (Dueker and Sheehan, 1997) with ANT-derived dispersion data, a high-resolution shear-wave velocity model can be obtained for the crust and uppermost mantle (Julià et al., 2000; Delph et al., 2015b; Delph et al., 2017). In this study, we build on this approach by combining longer period dispersion data obtained through TPWT to obtain shear-wave velocity results deeper into the upper mantle. Rayleigh waves at periods obtained from earthquake sources (generally >20 seconds) have very little sensitivity to upper-crustal structure, and thus it is important to pair them with shorter period data from ANT to better constrain crustal-velocity structure. By utilizing the joint inversion of receiver functions and ANT-derived dispersion data, we can obtain a crustal model that mitigates the need to assume an a priori crustal model that may be weakly constrained.

We follow the methodology described by Delph et al. (2015b) and Delph et al. (2017) to create CCP-derived receiver functions. First, we re-compute the high-quality receiver functions obtained by Abgarmi et al. (2017) using a Gaussian low-pass filter of ∼1.4 Hz (Gaussian pulse with a 2.8 alpha parameter). We then compute a CCP volume with bin spacing of 0.2° and a minimum bin width of 0.3° that can dilate until a minimum of 10 rays is incorporated into the bin or a bin width of 1° is reached. This approach creates a laterally smoothed and gridded 3D model of conversion amplitudes in the Earth. The receiver functions in this analysis are migrated with a two-layer velocity model consisting of crustal and mantle shear velocities of 3.4 km/s and 4.2 km/s with a Vp/Vs ratio of 1.78, which is representative of the average structure in central Anatolia (Delph et al., 2015a). This velocity model has a minimal effect on the resulting CCP-derived receiver functions, as the same velocity model is used to migrate the newly stacked receiver functions back to time (Delph et al., 2015b). A Gaussian filter is then used on the resulting CCP profiles to reduce the effects of spurious and unrealistic amplitude artifacts, which result from introducing new data into a CCP bin (see Delph et al., 2017, for details). The end result is equivalent to a P-wave receiver function filtered with a 1.2 Hz Gaussian (2.5 alpha parameter, ∼1 km vertical resolution) at each grid point. This receiver function can then be paired with a surface-wave dispersion measurement for the joint inversion.

To prevent inverting receiver-function multiples generated from the crust-mantle boundary (PpPs and PsPs + PpSs, respectively), we invert only the first 8 seconds of the receiver function (sensitive to depths <65 km). At depths greater than 65 km, only Rayleigh-wave dispersion data are used in the inversion for shear velocities. Because of the broad sensitivity kernels of Rayleigh waves (e.g., Fig. 6B), we increase our model layer thickness so that the diagonal of the resolution matrix is >0.1 with respect to the dispersion data. With the dispersion data used in this study, we can obtain shear-wave velocities down to ∼150 km (Figs. 6B and 7A). In order to create a smoothed, continuous volume of shear-wave velocities in the upper mantle where layers are thick, we perform a cubic interpolation between layer velocities for every depth profile. The resulting profile honors all layer velocities and does not significantly degrade fit to the dispersion data. Some of the stations of Abgarmi et al. (2017) were omitted from this analysis due to complicated waveforms at times corresponding to shallow crustal depths that cannot be adequately modeled by the joint inversion algorithm. Gaps in the resulting 3D model were filled using a linear interpolation between nearby data points for each layer. There were no gaps more than 0.4° away from any data point constrained by the joint inversion, and these gaps were largely focused in the western Central Anatolian Basin (see Supplemental Fig. S3 [^{footnote 1}] for fits to data sets and gap locations).

The new ambient noise tomography model that incorporates two years of CD-CAT and KOERI seismic data is similar to the ambient noise tomography data set of Delph et al. (2015a), with slightly faster phase velocities beneath the Central and Eastern Taurus Mountains at longer periods (>30 seconds). At periods where both ANT and TPWT obtain phase-velocity results (40, 45, and 50 seconds), the two methods return values consistent with each other (within ∼0.05 km/s on average; Table 1).

Phase velocities from TPWT and ANT at all periods are slow throughout most of central Anatolia relative to the global model AK135 (Fig. 6; Fig. 8A, anomaly 1). Dispersion curves from ANT and TPWT in the Arabian plate (AP), the Eastern Taurus Mountains, the Central Taurus Mountains, and the Central Anatolian Basin are shown in Figure 6 (locations are shown as colored stars in Fig. 8D). Beneath the northwestern Arabian plate, phase velocities are indicative of a thin crust and fast uppermost mantle as seen by the steep increase in phase velocities at short periods (green) (Fig. 6A; Supplemental Fig. S4 [see ^{footnote 1}]). At longer periods, phase velocities beneath the northern Arabian plate are similar to those in the Eastern Taurus Mountains and the Central Anatolian Basin, indicating that deeper mantle shear velocities in both regions may be similar. The slow phase velocities in the shortest periods (<20 seconds) for the Central Anatolian Basin are most likely due to the presence of thick sedimentary deposits. These phase velocities remain slow at periods sensitive to upper-mantle shear-wave velocities (>20 seconds) despite the region having a thin crust. These results for the Central Anatolian Basin, along with numerous other studies (e.g., Bakırcı et al., 2012; Salaün et al., 2012; Fichtner et al., 2013; Delph et al., 2015a; Govers and Fichtner, 2016), provide consistent evidence that the shear-wave velocities of the uppermost mantle in central Anatolia are slower than typical upper-mantle velocities. The slightly faster phase velocities beneath the Eastern Taurus Mountains at intermediate periods (30–59 seconds, orange curve in Fig. 6A) indicate that the lower crust and/or uppermost mantle of the Eastern Taurus Mountains must contain seismically faster material than the nearby Central Anatolian Basin, as the crust is ∼15 km thicker beneath the Eastern Taurus Mountains (Abgarmi et al., 2017). Lastly, dispersion data from the Central Taurus Mountains are consistent with a relatively thicker crust based on its low gradient at short periods. These dispersion data become much faster at intermediate and long periods (>30 seconds) compared to most of Anatolia, indicating the presence of fast lower-crustal and/or mantle shear velocities extending to significant depths. In fact, a fast phase-velocity anomaly underlies much of the Central Taurus Mountains (Fig. 8, anomaly 2). Abgarmi et al. (2017) proposed that the African lithosphere (or the “Cyprean slab”) was present beneath central Anatolia as far north as the northern margin of the Central Taurus Mountains based on intermediate depth (∼50–70 km) conversions in receiver functions. This is consistent with fast (>4.4 km/s) shear velocities obtained from the inversion of Rayleigh-wave dispersion data (Bakırcı et al., 2012) and full-waveform tomography (Govers and Fichtner, 2016). On a regional scale, longer-period phase velocities (≥100 seconds) show very little lateral variation, with a low-magnitude relative fast velocity anomaly underlying much of central Anatolia (Fig. 8D, anomaly 3).

Shear-wave velocities in the crust and uppermost mantle throughout the study area are generally slow (Fig. 7B), consistent with previous studies of shear velocity in the region (Bakırcı et al., 2012; Salaün et al., 2012; Fichtner et al., 2013; Delph et al., 2015a; Govers and Fichtner, 2016). The Arabian plate can be distinguished from the Anatolian plate by its thin crust and relatively fast uppermost-mantle shear-wave velocities (>4.4 km/s) compared to the thicker crust and seismically slower uppermost mantle of the Anatolian lithosphere (Fig. 9B). Crustal thickness estimates in the Arabian plate (30–35 km) agree well with a sharp increase in shear-wave velocities from ∼3.5 km/s to >4.4 km/s. These fast shear velocities extend down to ∼50–60 km where they are underlain by a large negative gradient, eventually reaching shear velocities that are much slower than expected for mantle lithosphere (<4.2 km/s at ∼60–110 km below the surface) (Figs. 9C, 9D, 10B, and 10D).

Farther west, there is a sharp transition from the fast velocities observed in the crust and mantle of the Arabian plate to very slow crustal shear velocities in the Adana Basin (<3.2 km/s) extending down to ∼17 km (Figs. 9A and 10A). Faster shear velocities at depths ∼30–60 km consistent with lithospheric mantle wave speeds extend below the eastern portion of the Adana Basin as well as locally beneath the Eastern Taurus Mountains. These fast velocities terminate to the west near the center of the Adana Basin, transitioning to slower velocities consistent with the lower crust of the Anatolian plate (Fig. 10B).

Beneath the Central Taurus Mountains and westernmost Adana Basin, we observe vertically extensive fast shear velocities (up to 4.8 km/s) beginning at a depth of ∼60 km (Figs. 9C, 9D, 10B, and 10C). These fast velocities appear to underlie much of the Central Taurus Mountains and extend as far north as the Central Anatolian Basin, where they transition rather abruptly to slow upper-mantle shear velocities. To the east, these fast velocities are separated from the fast velocities associated with the upper mantle of the Arabian plate by a narrow band of slow shear velocities that indicates that these are separate features (Fig. 10B).

Shear velocities in the upper crust of the Central Anatolian Basin are generally low (<3.4 km/s) due to the thick sedimentary deposits in the region (Fig. 9A). In the mantle, we observe the slowest shear velocities in the study area near the western edge of the Central Anatolian fault zone. These slow velocities show a NE-SW trend similar to that of the Central Anatolian Volcanic Province and reach velocities as low as 3.8 km/s beneath Mount Erciyes, the largest Quaternary stratovolcano in the Central Anatolian Volcanic Province (Figs. 9C and 10A). The spatial footprint of these slow mantle shear-wave velocities extends across a much broader region than just the Central Anatolian Volcanic Province, encompassing nearly the entire Central Anatolian Basin, the southern Kirşehir Block, and even the Eastern Taurus Mountains near the Central Anatolian fault zone.

The lithosphere-asthenosphere boundary (LAB) is a rheological boundary between the rigid lithosphere and ductile asthenosphere. In seismic studies, the LAB is generally associated with a gradational decrease in shear-wave velocities corresponding to a decrease in viscosity and possibly a small amount of melt (Eaton et al., 2009; Fischer et al., 2010). In this study, we investigate the character of the upper mantle using two separate approaches to provide evidence (or lack thereof) for the presence of subcontinental mantle lithosphere beneath central Anatolia.

In the first approach, we define the LAB as the depth of the maximum negative velocity gradient after a shear velocity of 4.4 km/s is reached, a velocity that we infer to be representative of lithospheric mantle. The lithospheric thickness of the Anatolian plate was found to be at a relatively constant depth of ∼80 km based on a recent regional-scale S-wave receiver-function study (Kind et al., 2015); however, at the finer scale resolved in this study, we observe a much more variable lithospheric thickness ranging from <50 km in the Arabian plate to ∼80 km in the Eastern Taurus Mountains and northern Kirşehir Block (Fig. 11). The thin lithosphere of the Arabian plate in this region is consistent with previous S-wave receiver-function results farther east (Angus et al., 2006). In Figure 11, the white regions indicate that the depth to the LAB is ambiguous due to shear velocities that do not reach our defined cutoff of 4.4 km/s as a proxy for lithospheric mantle (e.g., the Central Anatolian Volcanic Province) or as the result of a nearly monotonically increasing velocity-depth profile such that a negative velocity gradient isn’t obtained at a realistic depth for the LAB (e.g., the Central Taurus Mountains). It is important to note that if no mantle lithosphere exists beneath a region, the crust-mantle transition (Moho) would also represent the LAB. In this case, the negative gradient algorithm described above will not be useful in constraining the depth to the LAB, as the transition from crustal to asthenospheric mantle would be a velocity increase (no negative gradient) as opposed to the decrease that we usually associate with the LAB. This may be the case in some of central Anatolia.

Because shear-wave velocities in the uppermost mantle beneath central Anatolia can be very slow (Fig. 10C), the “negative-gradient” approach described above does not return results for the entire region. Thus, to provide evidence for the presence of lithospheric mantle throughout the entire region, we compute the average shear-wave velocity in the top 30 km of the mantle relative to crustal thickness estimates derived from receiver functions (Özacar et al., 2010; Vanacore et al., 2013; Abgarmi et al., 2017) (Supplemental Fig. S2 [see ^{footnote 1}]). We use a lower shear-velocity limit (4.3 km/s) as indicative of the presence of mantle lithosphere in this analysis to allow for uncertainties in crustal thickness estimates (which could incorporate crustal velocities into the mantle average if our estimated crustal thickness is shallower than the true crustal thickness) and the possibility that the lithospheric mantle is <30 km thick (Fig. 12A).

Based on the combination of these two analyses, there is very little evidence for the presence of mantle lithosphere beneath the Central Anatolian Volcanic Province. Thus, this thin to negligible mantle lithosphere allows for hot asthenospheric material to come into close proximity with the base of the Central Anatolian crust. Beneath the Central Taurus Mountains, the lithospheric thickness of the overriding plate is difficult to differentiate from the fast, vertically extensive shear velocities associated with the underthrusting African lithosphere (cross-hatched area, Fig. 11). Regardless, the presence of this fast-velocity material would protect the lithosphere of the Central Taurus Mountains from being exposed to asthenospheric material, making the interpretation that the Central Taurus Mountains are supported by shallow asthenosphere unlikely (e.g., Cosentino et al., 2012). Beneath the Arabian plate, the Eastern Taurus Mountains, and the northwestern Kirşehir Block, there is evidence from both approaches for the presence of a thin and seismically fast lithospheric mantle keel.

The strikes of neotectonic fault zones and structures appear to spatially correlate with boundaries between relatively thick and thin and/or negligible mantle lithosphere as inferred from the lateral gradient of the average upper mantle shear velocities (Fig. 12B). In particular, a semi-circular pattern of faulting is located in the Eastern Taurus Mountains at the boundary between slow and fast average upper-mantle shear velocities (Fig. 12, anomaly 1). Also, the N-S–trending thicker lithospheric mantle of the Eastern Taurus Mountains appears to correlate well with areas that show the highest density of neotectonic faulting (Fig. 12, anomaly 2). A similar pattern can be seen southwest of the Central Anatolian Basin, where NW-SE–striking fault orientations bound a thin zone of faster mantle shear velocities (Fig. 12, anomaly 3). The correlation between neotectonic structures and the transition from relatively thick to thin mantle lithosphere (as inferred from shear velocities) could imply that variations in lithospheric thickness play a crucial role in the location and concentration of deformation. Thus, the thickness of mantle lithosphere in Anatolia could contribute to the strength of the lithosphere as a whole, similar to what has been proposed in the western United States (Levander and Miller, 2012).

The magmatic products of the Central Anatolian Volcanic Province are thought to utilize regional tectonic structures that align with the Central Anatolian fault zone and Tuz Gölü fault (Toprak and Göncöoğlu, 1993; Toprak, 1998) (Fig. 2). As with deformation, volcanism throughout central Anatolia also seems to correlate with the transition between slow and fast upper-mantle velocities. This may be because strain in the crust and upper mantle is concentrated at these boundaries, leading to zones of weakness through which melt from the mantle is able to migrate relatively efficiently. In particular, Mount Erciyes, the largest stratovolcano in the Central Anatolian Volcanic Province, is located near the Inner-Tauride suture at the edge of a step at a releasing bend in the Central Anatolian fault zone (Higgins et al., 2015). The depth projection of the Inner-Tauride suture corresponds well with ∼40 km band of slow velocities that separate the fast upper mantle beneath the northern Kirşehir Block and Eastern Taurus Mountains (Fig. 10A, ITS/CAFZ). Thus, the relict lithospheric weakness associated with the Inner-Tauride suture (and now utilized by the Central Anatolian fault zone) represents an efficient pathway for mantle-derived melt to propagate to the surface, and may control the distribution of volcanism in the Central Anatolian Volcanic Province.

The very slow (<4.2 km/s) shear velocities observed beneath the Central Anatolian Volcanic Province support previous studies that indicate the presence of a hot upper mantle beneath the region (Gök et al., 2003; Biryol et al., 2011; Mutlu and Karabulut, 2011; Bakırcı et al., 2012; Salaün et al., 2012). The presence of melt beneath the Central Anatolian Volcanic Province is harder to constrain, as the effects of high temperatures versus the presence of melt on shear-wave velocity are difficult to separate. Through the numerical modeling of an upper-mantle system, Hammond and Humphreys (2000) estimated a decrease in shear velocity of ∼8% per 1% melt. Assuming the mantle lithosphere and asthenosphere are compositionally similar and that the shear velocities in lithospheric mantle would be 4.5 km/s at typical pressures and temperatures, our shear velocities of ∼4.1 km/s suggest a melt percentage of >1% beneath the Central Anatolian Volcanic Province (9% velocity decrease). Alternatively, neglecting the temperature effects and considering only the poroelastic effects of interstitial melts (Yoshino et al., 2005), ∼4%–6% melt could be present assuming a peridotitic upper-mantle composition based on geochemical results from basaltic magmas in central Anatolia (Delph et al., 2017; Reid et al., 2017). While the effects imposed by the combination of temperature and melt are not well quantified, recent experimental data on polycrystalline aggregates with analogous characteristics to the mantle olivine + basalt melt system found that it is difficult to obtain shear velocities below ∼4.2 km/s without temperatures crossing the solidus and presumably initiating melting (McCarthy et al., 2011). The anomalously slow shear velocities observed in the upper mantle beneath the Central Anatolian Volcanic Province, along with the aforementioned studies on the relationships between seismic velocities and melts, provide strong evidence that the upper mantle beneath central Anatolia likely contains some fraction of melt.

Evidence for the presence of the Cyprean slab below the southern margin of the Anatolian plate is accumulating based on upper-mantle seismic velocities and impedance contrasts beneath the Central Taurus Mountains (Bakırcı et al., 2012; Salaün et al., 2012; Govers and Fichtner, 2016; Abgarmi et al., 2017). The fast (>4.5 km/s) shear-wave velocities below the Central Taurus Mountains observed in this study corroborate this interpretation. Abgarmi et al. (2017) recently proposed that the rapid ∼2 km uplift of the Central Taurus Mountains was not due to the influx of asthenosphere beneath the region (Cosentino et al., 2012) but rather slab break-off and subsequent rebound after the detachment of the subducting oceanic lithosphere (Duretz et al., 2011; Schildgen et al., 2014). The removal of this load would allow for the uplift of the Central Taurus Mountains without requiring internal compressional deformation, preserving the near horizontal layering of late Miocene sedimentary rocks (Duretz et al., 2011; Cosentino et al., 2012). Although the appearance of the fast anomaly resolved in this study does not delineate the structure of the subducting slab in detail, it appears that the slab is present until the northern edge of the Central Taurus Mountains near the transition to the Central Anatolian Basin, where it is replaced by very slow shear velocities beneath the Central Anatolian Basin and the Central Anatolian Volcanic Province. To the east, the slab anomaly is replaced with anomalously slow shear-wave velocities in the upper mantle beneath the Adana Basin. These slow velocities appear to merge with the slow and shallow (∼60 km) shear-wave velocities beneath the Arabian plate (Fig. 10B). The discontinuous appearance of the fast anomaly along the southern margin of central Anatolia is consistent with the interpretation that the slab is actively deforming and fragmenting, leading to the development of slab tears and gaps (Biryol et al., 2011). The presence of attenuated continental material underthrusting beneath the southern margin of central Anatolia was recently verified by magnetic anomalies (Granot, 2016); thus, this slab deformation may be due to a recent transition from the underthrusting of old Tethyan oceanic lithosphere to the attenuated continental lithosphere of the African passive margin. The underthrusting of continental fragments appears to be common throughout the eastern Mediterranean (Schattner, 2010; Biryol et al., 2011; Delph et al., 2015a; Delph et al., 2015b), with the Anaximander Mountains and Eratosthenes Seamount representing un-subducted analogues.

While our seismic results only lend insight into the current character of the crust and upper mantle beneath Anatolia, they must be interpreted in light of the recent tectonic history of the region. In Figure 13, we interpret a cross section of our model based on a proposed model for the tectonic history of Central and Eastern Anatolia involving the rollback of a shallow dipping slab initiating in the Miocene (Bartol and Govers, 2014).

During the accretion and transfer of the continental Anatolide-Tauride Block (ATB) and Kirşehir Block (KB) from the downgoing plate to the overriding Eurasian margin, decoupling between their crust and mantle lithosphere may have allowed for the shallow subduction of mantle lithosphere as far north as the Izmir-Ankara-Erzincan suture zone (Bartol and Govers, 2014). Eocene–Oligocene contractional deformation is documented in the Kirşehir and Anatolide-Tauride blocks, indicating a period of compression related to this accretion (Çemen et al., 1999; Gürer et al., 2016). In the Miocene (Fig. 13A), contraction gave way to extension as seen by extensional exhumation in Kirşehir Block, which initiated by at least 17 Ma (Fayon and Whitney, 2007) and the rejuventation of deposition in the Tuz Gölü basin ca. 11 Ma (Fernandez-Blanco et al., 2013). Voluminous ignimbritic volcanism in the Central Anatolian Volcanic Province also initiated during this time. Felsic-intermediate ignimbrites dominate the northern extent of the Central Anatolian Volcanic Province in the Kirşehir Block while smaller mafic-intermediate composition volcanics are more prevalent in the Anatolide-Tauride Block to the south. A possible explanation for this spatial variation in composition is that more felsic material from the underthrusting of some of the Anatolide-Tauride Block’s crust beneath the Kirşehir Block along the Inner-Tauride suture was incorporated into mantle melt after slab delamination. Alternatively, the compositional variation may be due to the amount of time that the crust has been exposed to shallow, hot asthenosphere, as volcanism appears to have initiated earlier in the northern Central Anatolian Volcanic Province than in the south. Regardless, exhumation, extension, and volcanism indicate that the delamination of the mantle lithosphere from the accreted terranes initiated by at least mid-Miocene times, presumably foundering to the southwest as indicated by SKS splitting directions in central Anatolia (Paul et al., 2014) and the southwestward younging of the initiation of volcanism (Schleiffarth et al., 2015). This delamination event may also be largely responsible for the thin mantle lithosphere observed beneath most of the Kirşehir and Anatolian-Tauride blocks. As these blocks are thought to have decoupled from their mantle lithosphere during their accretion onto the Eurasian margin (Bartol and Govers, 2014; van Hinsbergen et al., 2016), the majority of the mantle lithosphere associated with these blocks was likely part of the delaminating slab. The removal of the slab led to the influx of asthenospheric mantle to very shallow depths (<50 km) beneath central Anatolia, as seen by the very slow shear-wave velocities imaged in this study.

Along the southern margin of central Anatolia, the initiation of uplift of the Central Taurus Mountains is attributed to slab break-off at ca. 8 Ma (Cosentino et al., 2012; Schildgen et al., 2012; Radeff et al., 2015) (Fig. 13B). This break-off may have been caused by the transition from the subduction of oceanic to continental material at the Cyprean trench. The subduction of this buoyant continental material resulted in the shallow subduction of the African lithosphere as far as the northern margin of the Central Taurus Mountains (Figs. 13B and 13C). Finally, uplift along the Central Taurus Mountains accelerated ca. 1.6 Ma as a result of the arrival of a relatively thicker fragment of continental material (the Eratosthenes Seamount) into the Cyprean trench (Schattner, 2010; Schildgen et al., 2012). An alternative possibility for explaining the late Miocene to recent uplift is the initiation of slab necking ca. 8 Ma caused by the introduction of continental material into the Cyprean trench ca. 8 Ma, followed by complete slab break-off ca. 1.6 Ma which led to the later phase of accelerated uplift in the Central Taurus Mountains.

The possibility of the presence of oceanic lithosphere between the Anatolide-Tauride Block and the African continental margin along the southern margin of central Anatolia, as seen further to the west along the Aegean and westernmost Cyprean trenches and inferred by the presence of ophiolites along the Bitlis-Zagros suture, could make it difficult to invoke a single-slab delamination model. If the lithosphere of the spreading ridge that likely separated the Anatolide-Tauride Block and African plate in the southern branch of the Neotethys was extinct for long enough to become rigid and act as a single plate, then the continuous subduction and subsequent delamination of a single slab throughout Cenozoic may be possible. However, if this was not the case, oceanic subduction could have initiated along the relict lithospheric weakness associated with the spreading ridge, and Fig. 13A could involve the presence of two slabs: the delaminating slab composed of Anatolide-Tauride and Kirşehir mantle lithosphere (shown in Fig. 13A), and a new slab composed of African lithosphere (not shown). An investigation of the deeper seismic structure interpreted alongside constraints from the surface geology is needed to better infer the Miocene evolution of subduction.

The CD-CAT seismic deployment, along with the addition of >100 stations from the Kandilli Observatory and Earthquake Research Institute at Boǧaziçi University, allows for unprecedented spatial sampling of the crust and upper mantle in central Anatolia. We utilize a recently developed approach to the joint inversion of P-wave receiver functions, ambient noise-derived dispersion data, and earthquake-generated dispersion data to obtain a 3D shear-wave velocity model down to 150 km to understand the interplay between the recent tectonic history of the Central Anatolian plate and the character of subduction along its southern margin.

The main findings of this research are:

• 1. The thickness of the lithosphere beneath Anatolia is variable but generally thin. Beneath the Arabian plate, Eastern Taurus Mountains, and northwestern Kirşehir Block, lithospheric thickness ranges from ∼50–80 km, while very little evidence for mantle lithosphere based on shear-wave velocities is found beneath the Central Anatolian Volcanic Province. As the Kirşehir and Anatolide-Tauride blocks were initially attached to the downgoing plate, most of their mantle lithosphere was most likely removed as part of a Miocene slab delamination event. Gradients in lithospheric thickness also correlate with neotectonic fault density and Miocene to recent volcanism, implying that lithospheric thickness may be playing a crucial role in controlling the distribution of neotectonic structures and volcanism at the surface.

• 2. The fastest shear-wave velocities in the region are located beneath the Central Taurus Mountains. The presence of these fast velocities provides further evidence that the subducting Cyprean slab underlies the Central Taurus Mountains, which have been uplifted ∼2 km in the past 8 m.y. The presence of this seismically fast material makes the interpretation that uplift was driven by the influx of low-density asthenosphere unlikely. Rather, we propose that the majority of the uplift of the Central Taurus Mountains is due to slab rebound after the detachment of the Cyprean slab beneath Anatolia, and we interpret that slab necking/break-off was caused by a transition in the composition of the material arriving at the Cyprean trench from oceanic to attenuated continental in the late Miocene.

• 3. The fast seismic velocities interpreted as the underthrusting Cyprean slab appear to terminate at the northern margin of the Central Taurus Mountains. These fast velocities give way to a NE-SW trend of very slow upper-mantle shear-wave velocities (<4.2 km/s) that correlate well with the Central Anatolian Volcanic Province. These slow velocities are interpreted to be shallow asthenosphere in which the presence of melt is required. This shallow asthenosphere combined with the lithospheric-scale weakness of the Inner-Tauride suture is responsible for the distribution of volcanism in the Central Anatolian Volcanic Province. Also, the change in composition of the volcanism in the Central Anatolian Volcanic Province from large volume silicic ignimbrites in the north to smaller mafic to intermediate volcanics in the south, could be due to the melting of highly silicic re-laminated material underthrust beneath the Kirşehir Block. Alternatively, this compositional variation could reflect the amount of time that the overriding plate has been exposed to the hot, shallow asthenosphere.

This research was supported by National Science Foundation grant 1109336. J.R.D. also received support for this project through the Weiss Post-Doctoral Fellowship at Rice University. The authors want to thank Rob Govers and an anonymous reviewer for their helpful comments, which improved this manuscript, as well as the CD-CAT collaborators for useful and insightful discussions. The success of the seismic deployment could not have been accomplished without the assistance of Angela “Mouse” Reusch (The Incorporated Research Institutions for Seismology [IRIS] Portable Array Seismic Studies of the Continental Lithosphere [PASSCAL]) and our Turkish collaborators C. Berk Biryol (University of North Carolina [UNC]), Uğur Teoman (Kandilli Observatory and Earthquake Research Institute [KOERI]), Metin Kahraman (KOERI), Zafer Kaplan (Middle East Technical University), and Savaş Ceylan (Eidgenössische Technische Hochschule Zürich). The authors want to thank Savaş Ceylan for providing a code for surface-wave windowing, as well as Lara Wagner (Carnegie Department of Terrestrial Magnetism) and Sanja Antonijevic (UNC) for providing insight into the two plane-wave tomography codes. The facilities of IRIS Data Services, and specifically the IRIS Data Management Center, were used for access to waveforms and related metadata used in this study. IRIS Data Services are funded through the Seismological Facilities for the Advancement of Geoscience and EarthScope (SAGE) Proposal of the National Science Foundation under Cooperative Agreement EAR-1261681. Waveform data were also obtained from the GeoForshungsZentrum (GFZ) German Research Center for Geosciences and KOERI. The figures in this paper were created using Generic Mapping Tools (GMT; Wessel et al., 2013), and data were analyzed using Seismic Analysis Code (SAC; Goldstein and Snoke, 2005).