Evolution of mafic lavas in Central Anatolia: Mantle source domains

Continental assembly and breakup are key components of crustal evolution, and their effects are seen in the rocks of orogenic belts across the globe. The Alpine-Himalayan orogen is one such region, in which complex convergence of continental masses has been accompanied by the growth of major mountain belts and followed locally by periods of extensional collapse. Central and Western Anatolia are regions of post-collisional volcanism within this major orogen where both calc-alkaline and alkaline mafic lavas have erupted episodically from the early Miocene to the present. The transition from convergence to extension, and the complex evolution of the continental lithosphere and sublithosphere that accompanied this change, are manifest in the geochemical and isotopic signatures of Anatolian mafic lavas that are the focus of this work.

The thermal effects of rifting or post-orogenic crustal relaxation can initiate melting of hydrated lithosphere that was metasomatized during previous subduction events, given that metasomatic fluids in suprasubduction zones can produce abundant veins rich in fluid-mobile trace elements within the lithospheric mantle. In subduction zones, both slab-derived fluids and sediment-mélange packages may mix with mantle wedge material and contribute to the overall chemistry of continental volcanics (e.g., Behn et al., 2011; Nielsen and Marschall, 2017). At a larger scale, lithospheric thickening can occur in any magmatic environment through the underplating of cumulates, and in subduction settings through the accretion of slabs to the underside of an existing crustal section. Given the propensity of continental material to rift along pre-existing zones of convergent and transverse movement, intraplate extensional mafic lavas likely interact with lithospheric materials that record earlier subduction processes and thus the geochemical consequences of continental evolution.

Recent geophysical studies have resolved a complex network of subducted materials beneath Anatolia. Tomographic images reveal fast seismic anomalies beneath Anatolia that correspond to the subducted portion of the African lithosphere along the Cyprean and Aegean trenches (Biryol et al., 2011; Bartol and Govers, 2014; Govers and Fichtner, 2016). The two slabs are clearly separated south of Anatolia; both dip to the north, and they merge at depth beneath mid-latitude Anatolia. It is significant that tomograms from depths below 500 km show that the separate fast anomalies (Cyprean and Aegean slabs) merge to form a continuous belt that extends beneath all of Central Anatolia and flattens out along the mantle transition zone (∼660 km depth). The implication is that this feature is a single slab that sank beneath the present-day İzmir-Ankara-Erzincan suture, and Bartol and Govers (2014) suggested it represents delaminated lithospheric mantle (i.e., the northern Neotethys slab) that sank into the deeper mantle in the early Miocene (ca. 20–15 Ma). Delamination of this slab was followed by slab rollback to the present location south of Cyprus, while in Western Anatolia, it has resulted in extension both on land and in the eastern Aegean Sea. Slabs beneath both Central and Western Anatolia show geophysical evidence of tearing that allows warm asthenosphere to penetrate to shallow levels both through and around the subducted materials.

In this study, we report new major and trace element abundances and Sr-Nd-Pb-Hf isotopic data for late Miocene to recent mafic volcanic rocks from Central Anatolia and interpret them within the geodynamic context of complex slab configurations beneath the region; our new data are considered in the context of published data for Miocene and younger lavas from both Central and Western Anatolia (Alıcı et al., 1998, 2002; Kürkçüoğlu et al., 1998, 2001; Parlak et al., 2001; Şen et al., 2004; Chakrabarti et al., 2012; Aldanmaz et al., 2015; Reid et al., 2017; Di Giuseppe et al., 2018; Kocaarslan and Ersoy, 2018). Our goal is to provide a better understanding of the mantle source domains that contribute to volcanism in these areas, and specifically to explore relationships between major lithospheric and sublithospheric features as they relate to the tectonic evolution of Anatolia. Recent isotopic studies of Anatolian mafic lavas have attributed their genesis to contributions from depleted mid-ocean-ridge (MORB)–like mantle and a subduction modified oceanic island basalt (OIB)–like source domain containing fluids or melts associated with subducted sediments (e.g., Chakrabarti et al., 2012; Reid et al., 2017). Here, we extend the geographic and temporal scope of mafic lavas under consideration, allowing us to document source evolution and distribution in a complex tectonic environment affected by regional-scale lithospheric delamination subsequent to the closure of Neotethys and the onset of southwestward movement of the Anatolian microplate (Bartol and Govers, 2014; Göğüs et al., 2017). Ultimately, the insights provided through the geochemical characteristics of Anatolian mafic volcanism allow us to unravel the geochemical and tectonic development of dynamic orogenic and post-orogenic regions.

Anatolia is part of the Alpine-Himalayan orogenic belt, formed by tectonic juxtaposition of pieces of continental lithosphere during the Cenozoic closure and terminal suturing of the northern and southern branches of the Neotethys Ocean. Between the Cretaceous and late Eocene (ca. 100–34 Ma), the majority of the convergence between Africa and Arabia-Eurasia was accommodated by northward subduction of the northern Neotethys Ocean along the İzmir-Ankara-Erzincan suture zone (IAESZ; Fig. 1; Şengör and Yılmaz, 1981; Okay and Tüysüz, 1999; Bartol and Govers, 2014). By the late Paleocene, the northern portion of the ocean was closed along this suture, allowing for the terminal collision of the Eurasian plate with the Anatolian-Iranian Platform, which had approached from the south. The North Anatolian fault zone follows this northern Neotethyan suture, while more complex subduction of the southern Neotethys is manifest in the Aegean and Cyprean trenches immediately south of Anatolia.

Subduction of the African plate and closure of the southern Neotethys began south of Cyprus by the early Miocene, continuing through the mid-Miocene in the east and the early Quaternary in the west (Dewey et al., 1989). Final collision along the Bitlis-Zagros suture zone (Fig. 1) halted northward movement of the Arabian plate with respect to Africa and thickened the crust across the Anatolian-Iranian Platform by as much as 2 km (e.g., Yılmaz et al., 1993). Northward subduction of African lithosphere occurs at present along the Aegean and Cyprean trenches, which are separated by the Isparta Angle in southern Anatolia. Subduction along the Cyprus section of the trench is influenced by the collision of the Eratosthenes Seamount, which obstructs descent of the African plate, leading to steep subduction and tearing of the downgoing slab (Wortel and Spakman, 1992).

The interplay between northern Neotethyan slab delamination and southern Neotethyan subduction sets the stage for modern volcanism across Anatolia. Subsequent to the regional-scale delamination event, an influx of asthenosphere into the uppermost mantle resulted in broad heating and uplift of the Anatolian plateau, manifest today in an anomalously thin lithospheric mantle (30–40 km; Govers and Fichtner, 2016) and high heat flow (e.g., Keskin, 2003). Although delamination may have occurred at different times across Anatolia, there is no compelling evidence to support an age progression (Bartol and Govers, 2014). In Western Anatolia, post–mid-Miocene extension and volcanism have been broadly associated with slab rollback along the Hellenic Trench since ca. 12 Ma (Meulenkamp et al., 1988) and lithospheric delamination (Biryol et al., 2011; Bartol and Govers, 2014). Slab rollback across the Aegean extensional province reflects the differential convergence rates between Africa and Europe (∼10 mm/yr) and along the Hellenic Trench (∼40 mm/yr; McClusky et al., 2000). Thus, early Miocene mafic volcanism in Western Anatolia followed almost immediately upon the heels of delamination, while the neotectonic environment of Western Anatolia is dominated by extension related to rollback of the Hellenic-Aegean slab.

Mio-Pliocene through recent faulting and volcanism in Central Anatolia are consequences of internal extension governed primarily by “escape tectonics” as the microplate moves westward at a rate of ∼15 mm/yr along major strike-slip fault systems (Reilinger et al., 1997) with ascent of mafic melts along crustal faults (e.g., Notsu et al., 1995). The most recent phase of volcanism in Central Anatolia ensued several million years following regional slab delamination, presumably related to slow heating of the lower lithosphere, and post-dates two pulses of ignimbrite activity dated at 9–8 Ma and 7–5 Ma in the Cappadocian volcanic province (Innocenti et al., 1975; Temel et al., 1998; Aydar et al., 2012). This phase of activity began with growth of multi-stage Miocene through Quaternary calc-alkaline andesitic stratovolcanoes associated with transtensional pull-apart basins including Erciyes and Hasandağ (Kürkçüoğlu et al., 1998; Şen et al., 2004; Fig. 1) and continued through the Quaternary at scattered mafic and felsic vents with an increasing representation of alkali olivine basalts (e.g., adjacent to Hasandağ and Karapınar, 0.6–0.2 Ma; Reid et al., 2017).

We present new geochemical and isotopic analyses of the mafic volcanic rocks from several volcanic areas including Sivas, the northernmost magmatic province in Central Anatolia (Fig. 1). Parlak et al. (2001) referred to Sivas volcanism overall as Plio-Pleistocene and the basalts as Plio-Quaternary, but K-Ar ages of 3.3 Ma (Türkecan et al., 2000) and 4.0 ± 1.2 to 5.88 ± 0.77 Ma (Platzman et al., 1998) suggest volcanism developed in latest Miocene to early Pliocene time. Sivas lacks a stratovolcanic edifice but rather comprises numerous small-volume mildly alkaline mafic flows (Parlak et al., 2001; Kocaarslan and Ersoy, 2018). Primitive basanites appear to be restricted to the northwestern portion of the study area, while basaltic through trachyandesitic lavas are common in the southeast; the two eruptive groups are separated by the Kızılırmak River, which overlies a major strike-slip fault associated with the IAESZ (Kürkçüoğlu et al., 2015). Compositional differences among Sivas lavas are variously attributed to crustal contamination (Parlak et al., 2001; Kocaarslan and Ersoy, 2018) or distinct thermal and magmatic environments across the suture (Kürkçüoğlu et al., 2015).

Our data set also includes two new basalts from the youngest Plio-Quaternary phase of activity at Erciyes stratovolcano (Notsu et al., 1995) and undated mafic lavas from Erkilet and Develidağ (Fig. 1). The Erkilet and Develidağ volcanic complexes are commonly referred to as Mio-Pliocene in age (e.g., LePennec et al., 1994), although no comprehensive age determinations have been published. Dönmez et al. (2003) reported a K-Ar age of 5.3 Ma on a Develidağ andesite, and Türkecan et al. (1998) reported a K-Ar age of 3.1 Ma on an Erkilet lava but provided no compositional or stratigraphic metadata.

In this paper, we use the Sr-Nd-Pb-Hf radiogenic isotope compositions of early Miocene through Quaternary mafic lavas from Central and Western Anatolia to determine the origin and geodynamic evolution of their mantle sources and to relate their occurrence to the regional tectonic history of continental assembly and post-orogenic extension. These results enable us to compare mantle source evolution across Anatolia, and in particular, to explore the contribution from subducted sediments and to identify the geochemical and isotopic end members involved.

Samples of visibly fresh lavas were cut into slabs with any weathered portions removed and reduced to ∼5–7 mm pieces using an alumina ceramic mini-crusher. Powders were prepared by grinding in a tungsten carbide disc mill for 30–90 s. Whole-rock analyses for major and minor elements (including Ba and Sr) were obtained by direct current plasma spectrometry on an ARL-Fisons Spectraspan 7; remaining trace elements were analyzed by inductively coupled plasma mass spectrometry using a VG PlasmaQuad-3 at Duke University (Durham, North Carolina, USA). These data are presented and discussed in Kürkçüoğlu et al. (2015).

Clean chips of selected whole-rock samples, prepared with stainless steel tools, were separated for Sr-Nd-Pb-Hf isotopic analysis at San Diego State University (SDSU; San Diego, California, USA). Rock chips were cleaned first in ultrapure water (18 MW H₂O) followed by distilled 1N HBr, H₂O rinse, then distilled 2.5N HCl, and a final rinse in H₂O. They were then dissolved in a distilled HF-HNO₃ solution in sealed radius-bottom beakers perfluoroalkoxy (PFA) Teflon beakers on hot plates for a period of three days. Techniques for Pb, Sr, Hf, and Nd separation using anion exchange chromatography are described in Hanan et al. (2004, 2008). Sr isotopes were measured on a VG Sector 54 seven-collector thermal ionization mass spectrometer, and Hf, Nd, and Pb isotopes were analyzed on a Nu 1700 multi-collector inductively coupled plasma mass spectrometer following the methods of Hanan et al. (2004, 2008).

Measured ⁸⁷Sr/⁸⁶Sr isotope ratios were normalized to ⁸⁸Sr/⁸⁶Sr = 0.1194 to correct for mass fractionation, and the Sr isotopic data are reported relative to the NIST SRM 987 standard (⁸⁷Sr/⁸⁶Sr = 0.710250). The Nd and Hf isotope ratios were corrected for instrumental mass fractionation and machine bias by application of a discrimination factor determined by bracketing sample runs with analyses of the SDSU AMES Nd and the JMC 475 standards, respectively (every two or three samples). Hf isotope sample data were normalized to ¹⁷⁹Hf/¹⁷⁷Hf = 0.7325 and are reported relative to JMC 475 (¹⁷⁶Hf/¹⁷⁷Hf = 0.282162); Nd isotope sample data were normalized to ¹⁴⁶Nd/¹⁴⁴Nd = 0.7219 and are reported relative to the lab value of ¹⁴³Nd/¹⁴⁴Nd = 0.512130 for the SDSU AMES Nd standard. The measured value of ¹⁴³Nd/¹⁴⁴Nd of the La Jolla Nd standard at SDSU is 0.511844. Pb isotope ratios were monitored for isotope fractionation using the NIST SRM 997 Tl standard (White et al., 2000) and were corrected for instrumental mass fractionation and machine bias by applying a discrimination factor determined by bracketing sample analyses with analyses of NIST standard SRM 981 (every one or two samples) using values determined by Todt et al. (1996). The external reproducibility for ⁸⁷Sr/⁸⁶Sr and ¹⁴³Nd/¹⁴⁴Nd is 2 and 3 ppm, respectively; for ²⁰⁶Pb/²⁰⁴Pb, ²⁰⁷Pb/²⁰⁴Pb, and ²⁰⁸Pb/²⁰⁴Pb it is 21, 36, and 44 ppm, respectively; and for ¹⁷⁶Hf/¹⁷⁷Hf it is 2 ppm, based on replicate analyses of standards. Analytical results for procedural blanks were: <250 pg Sr, <200 pg Nd, <90 pg Pb, and <20 pg Hf. No blank corrections were applied to the data because they would have been insignificant.

Young Sivas lavas comprising basalts through basaltic trachyandesites and subordinate basanites (Fig. 2) have been analyzed previously for their major and trace element compositions (Parlak et al., 2001; Kürkçüoğlu et al., 2015); five new analyses are included in Table 1. Incompatible trace element abundances of Sivas basanites show elevated Ba, Th, and U relative to the basalts, corresponding negative Nb, Ta, and K anomalies, steeper rare earth element (REE) profiles (La/Sm_n ∼4, where n refers to the chondrite-normalized abundance [normalizing values of Sun and McDonough, 1989]), and enrichment in REEs over high-field-strength elements of similar compatibility (Fig. 2).

Table 1 presents 11 new analyses of late Miocene–Pliocene basalts through basaltic trachyandesites from Erkilet volcano in north-central Anatolia (Fig. 1). The Erkilet lavas include seven geochemically homogeneous basalts with ∼7 wt% MgO; their CaO/Al₂O₃ (∼0.6) values overlap those of Sivas basalts with similar MgO, suggesting a similar crystallizing assemblage. Erkilet basalts have lower incompatible trace element abundances than Sivas basalts and show much less geochemical variability overall. Their primitive mantle–normalized incompatible trace element patterns have positive Ba anomalies and marked negative Nb-Ta anomalies.

From Central Anatolia, we analyzed the Sr-Nd-Pb-Hf isotopic compositions of 11 Sivas lavas including basanite, basalt, trachybasalt, and basaltic trachyandesite (Kürkçüoğlu et al., 2015). The new isotopic data presented in Table 2 also include 11 Miocene and younger lavas from Central Anatolia, specifically: five Erkilet lavas, four late Miocene–Pliocene Develidağ basalts, and two Quaternary Erciyes lavas; and new analyses of two basaltic andesites from Tendürek in Eastern Anatolia. The new Anatolian data have a wide range in isotopic values, with ²⁰⁶Pb/²⁰⁴Pb of 18.667–19.320, ²⁰⁷Pb/²⁰⁴Pb of 15.638–15.688, ²⁰⁸Pb/²⁰⁴Pb of 38.725–39.354, ⁸⁷Sr/⁸⁶Sr of 0.70366–0.70590, ¹⁴³Nd/¹⁴⁴Nd of 0.51257–0.51282, and ¹⁷⁶Hf/¹⁷⁷Hf of 0.28000–0.28302 (Table 2).

The data show a striking temporal variation in the binary ε_Nd versus ε_Hf plot (Fig. 3A) from Miocene through Quaternary time. Quaternary lavas from Central (Develidağ, Erciyes, Sivas, Hasandağ cinder cone province, Karapınar), Eastern (Tendürek), and Western Anatolia (Kula) and Central Anatolian late Miocene to Quaternary lavas from Erkilet define a broad field with intermediate ε_Hf (+4.11 to +9.50) and ε_Nd (−1.34 to +5.74) that extends above the terrestrial array to higher ε_Hf for a given ε_Nd relative to the Miocene data trend which falls about the terrestrial array (Fig. 3A; J. Blichert-Toft, personal commun., 2020). Two Sivas basanites plot along the terrestrial array, substantially closer to the field of northeastern Mediterranean sediments than the other Central Anatolian samples.

The ⁸⁷Sr/⁸⁶Sr versus ε_Nd plot (Fig. 3B) shows a temporal variation relationship analogous to that shown for ε_Nd versus ε_Hf (Fig. 3A). Our new samples collectively have Sr-Nd isotopic values that fall within the previously published range of Miocene and younger Central Anatolian mafic lavas (Wilson et al., 1997; Kürkçüoğlu et al., 2001; Parlak et al., 2001; Alici et al., 2002; Reid et al., 2017; DiGiuseppe et al., 2018; Gall et al., 2021). The Sr-Nd isotopic ranges for each volcanic area are restricted and define strong correlations (Fig. 3B). Sivas ⁸⁷Sr/⁸⁶Sr values range from 0.70399 to 0.70587 and ¹⁴³Nd/¹⁴⁴Nd from 0.51257 to 0.51280 (ε_Nd +3.12 to −1.34; Fig. 3); high-MgO basanites tend toward more-enriched Sr values. Develidağ and Erkilet basalts extend to slightly lower ⁸⁷Sr/⁸⁶Sr (0.7037) and higher ε_Nd (+3.61) values, while Erciyes basalts have ε_Nd between +0.5 and +1.6 and ⁸⁷Sr/⁸⁶Sr from ∼0.704 to 0.705.

The Pb isotope signatures are consistent with those of the Sr, Nd, and Hf isotopes, showing temporal and spatial variations on Pb-Pb binary plots (Figs. 3C, 3D). In ²⁰⁶Pb/²⁰⁴Pb versus ²⁰⁷Pb/²⁰⁴Pb or ²⁰⁸Pb/²⁰⁴Pb space, the new data on late Miocene through Quaternary mafic lavas from Central Anatolia define a trend with a positive slope that is roughly parallel to the Northern Hemisphere reference line (NHRL) but at higher ²⁰⁷Pb/²⁰⁴Pb (15.597–15.657) and ²⁰⁸Pb/²⁰⁴Pb (38.503–39.334) for a given ²⁰⁶Pb/²⁰⁴Pb (18.667–19.320).

Major element variations suggest that evolution of the Quaternary Sivas basaltic suite is controlled by as much as ∼65% fractional crystallization of olivine (∼16%), clinopyroxene (∼35%), plagioclase feldspar (∼40%), and Fe-Ti oxides (∼9%) from a parental liquid with ∼9 wt% MgO and ∼47% wt% SiO₂ (Kürkçüoğlu et al., 2015). These lavas have primitive mantle–normalized incompatible trace element patterns with prominent negative anomalies for Nb, Ta, and K (Fig. 2) and steep REE profiles with La/Sm_n between 2.5 and 3.7. Pliocene Sivas basanites are not related to one another or to the basaltic lavas through crystal fractionation, i.e., P₂O₅ and K₂O contents of primitive basanites are elevated relative to those of the basalts and correlate positively with MgO (Kürkçüoğlu et al., 2015). Kürkçüoğlu et al. (2015) interpreted these samples as individual small-volume melts with significant contributions from a lithospheric source region metasomatically modified by fluids. Petrogenesis of the Erkilet lavas is not investigated here in detail. Their geochemical overlap with Sivas lavas suggests a comparable fractionating assemblage, while the shallow slope of their REE profiles (La/Sm_n ∼2.0; Fig. 2) suggest they derive from melting within the spinel stability field.

Mildly alkalic mafic lavas from Central Anatolia were previously interpreted as asthenosphere-derived primary basaltic liquids that are heavily overprinted by crustal contamination (e.g., Parlak et al., 2001; Kocaarslan and Ersoy, 2018) or by interaction with mantle that has been modified by input of subduction-derived fluids (Notsu et al., 1995; Reid et al., 2017). Our data are not consistent with crustal contamination among parental mafic lavas, given that the Sr isotopic values of Central Anatolian lavas correlate poorly with MgO: Sivas lavas with >9 wt% MgO have the highest ⁸⁷Sr/⁸⁶Sr values (Fig. 4). This observation speaks against significant crustal contamination but rather supports an amended version of the latter model in which melts (and fluids?) of recycled continent-derived material located within the mantle dominate the isotopic and trace element signatures of the mafic lavas. Here we explore the mantle source domains and the temporal evolution of their contributions to Central Anatolian magmatism.

The entire range of Central Anatolian Sr-Nd-Hf isotopic data falls within the broad array defined by Miocene and younger lavas from Western Anatolia, which encompasses widespread early and mid-Miocene lavas as well as late Miocene through Pliocene samples from Kula and Thrace (Aldanmaz et al., 2000, 2006, 2015; Alıcı et al., 2002; Chakrabarti et al., 2012). The ages of mafic volcanism in Western Anatolia are better constrained than those from Central Anatolia; samples reported here are categorized following whole-rock K-Ar determinations by Aldanmaz et al. (2000, and references therein). The early Miocene lavas have the highest ⁸⁷Sr/⁸⁶Sr (0.70835–0.70860) and the lowest ε_Nd (−6.24 to –5.89). Middle Miocene lavas range to less-radiogenic ⁸⁷Sr/⁸⁶Sr (0.70789–0.70868) and more-radiogenic ε_Nd (−5.19 to −4.68), while the late Miocene lavas cluster tightly at much-less-radiogenic ⁸⁷Sr/⁸⁶Sr (0.70311–0.70354) and more-radiogenic ε_Nd (+5.68 to +6.63). The late Miocene to Quaternary lavas from Kula (Aldanmaz et al., 2006; Chakrabarti et al., 2012) define a steeper trend with more variable ⁸⁷Sr/⁸⁶Sr values between 0.70303 and 0.70619 (Fig. 3B). In Western Anatolia, the early Miocene lavas of Foça and Yenişakran have lowest ε_Hf (−5.06 to −2.02) and ε_Nd (−6.24 to −5.89), increasing in the middle Miocene lavas of Foça, Ayvacık, and Bergama to ε_Hf of −1.52 to +0.46 and ε_Nd of −5.19 to −4.68), and jumping to very high ε_Hf (+8.84 to +9.85) and ε_Nd (+5.72 to +6.63) in the late Miocene lavas of Čanakkale, Thrace, Ayvacık, and Taştepe, similar to that of depleted mantle sources.

The Anatolian data set for Pb isotopic variations shows complex variations, given that the data do not define a single linear field or array in binary Pb-Pb space (Fig. 3). Specifically, late Miocene and younger lavas from Central Anatolia, Quaternary Western Anatolia mafic lavas, and Tendürek lavas in Eastern Anatolia define a trend in Pb-Pb isotopic space that is oblique to that of early Miocene through Pliocene Western Anatolian basalts. These Miocene lavas from Western Anatolia plot along negatively sloped trends that extend from the late Miocene data which plot near the NHRL (²⁰⁶Pb/²⁰⁴Pb 18.910–19.102; ²⁰⁷Pb/²⁰⁴Pb 15.580–15.616; ²⁰⁸Pb/²⁰⁴Pb 38.649–38.753), through the middle Miocene data (²⁰⁶Pb/²⁰⁴Pb 18.616–18.712; ²⁰⁷Pb/²⁰⁴Pb 15.686–15.678; ²⁰⁸Pb/²⁰⁴Pb 39.004–39.016), to the early Miocene data which plot at ²⁰⁷Pb/²⁰⁴Pb (15.689–15.693) and ²⁰⁸Pb/²⁰⁴Pb (39.005–39.048) and ²⁰⁶Pb/²⁰⁴Pb (18.589–18.675). We restrict our modeling and detailed interpretations to these high-precision data and published isotopic results for Western Anatolia lavas that have been comprehensively analyzed at San Diego State University using the procedures outlined herein (Aldanmaz et al., 2015).

The trends are seen most clearly in ternary Pb isotope space (Fig. 5) but are apparent in binary plots as well (Fig. 3). We recognize that the small number of samples from Eastern Anatolia makes a broad interpretation impossible but note that the high-precision data presented here are consistent with our geodynamic model for late Miocene through Quaternary lavas from across Anatolia. The positive Δ7/4Pb and Δ8/4Pb (i.e., the displacement of the samples from the NHRL; Hart, 1984), high ⁸⁷Sr/⁸⁶Sr, and low ε_Nd observed for these lavas require a source component with older continent-derived Pb, Sr, and Nd isotope contents, i.e., higher values of Rb/Sr and lower Sm/Nd than the convecting mantle. In Sr-Nd isotope space, variations observed in Central Anatolia are analogous to those documented in the Gulf of Aden (the Sheba Ridge and its extension into the Gulf of Tadjoura; Fig. 3) where decreasing ε_Nd (from the Sheba Ridge to the Tadjoura Trough, Asal Rift, and Tadjoura Rift) reflects variable pollution of the upper mantle MORB source by recycled continent-derived material originating from the subcontinental lithospheric mantle (SCLM) and recycled oceanic lithosphere (Schilling et al., 1992; Rooney et al., 2011).

We interpret the Western Anatolian Miocene ²⁰⁶Pb/²⁰⁴Pb-²⁰⁷Pb/²⁰⁴Pb isotopic data array (Fig. 3) as a mixing trend and note that the intersection of the trend with the Gulf of Aden–Sheba Ridge MORB field and the NHRL characterizes the ambient upper-mantle isotope signature beneath Anatolia. This composition has distinctly more-radiogenic Pb isotopic values (and less-radiogenic Sr isotopic values) than the Quaternary Kula basalts that were interpreted by Chakrabarti et al. (2012) and Reid et al. (2017) to represent the depleted ambient end member. The Miocene trend can therefore be described as a mixture of two mantle sources, i.e., a grossly homogeneous ambient MORB source with ²⁰⁶Pb/²⁰⁴Pb of ∼19, ²⁰⁷Pb/²⁰⁴Pb of ∼15.6, and ²⁰⁸Pb/²⁰⁴Pb of ∼38.6 and a heterogeneous lithospheric domain with lower ²⁰⁶Pb/²⁰⁴Pb but higher ²⁰⁷Pb/²⁰⁴Pb and ²⁰⁸Pb/²⁰⁴Pb that was modified by subduction-related contributions including direct melts of sediments and/or mélange.

Multi-component mixing is required to explain the isotopic signatures in late Miocene through Quaternary time as represented by lavas of these ages from Central and Eastern Anatolia and Quaternary basalts from Kula in Western Anatolia. Although the model parameters (Supplemental Table S1¹) and results (Fig. 6; Fig. S1) are not unique, the end product is consistent with the available data and provides a geologically reasonable framework within which to interpret the genesis of the Anatolian suites. It takes at least three end-member components to bracket the compositional range of the Anatolian isotope data. We focus our discussion in Pb isotope space because variations of Pb isotopic values are very sensitive to mixing involving continent-derived recycled material. The identification of the sublithospheric ambient mantle as one end member requires a geometry where the remaining two end members have high Δ7/4Pb and Δ8/4Pb and plot in the Pb isotope region of continent-derived materials. The join between the two end members in Pb isotope space represents the range of heterogeneous source components (which vary in isotope composition and age) that have mixed with material derived from the ambient mantle to form the Anatolian lavas. One end member is characterized by radiogenic Pb isotopes with ²⁰⁶Pb/²⁰⁴Pb >19 similar to that of modern oceanic sediment, the other with ²⁰⁶Pb/²⁰⁴Pb ∼18 similar to that of ancient sediment and/or continental lithospheric mantle.

We interpret the Western Anatolian late Miocene lavas with ²⁰⁶Pb/²⁰⁴Pb of ∼19 to represent a sublithospheric component located in the upper mantle. For ease in modeling and visualization, we have chosen a single point located in the region of overlap between the Anatolian late Miocene data, Gulf of Aden–Sheba Ridge data, and NHRL to represent the ambient mantle isotope composition. Each of the Quaternary basalts represents a binary mix of a specific lithospheric composition along the join between the two continent-derived end members and this sublithospheric ambient mantle composition (see Table S1 [^{footnote 1}]). The Anatolian lithosphere is envisioned as being isotopically heterogeneous, having undergone subduction-related enrichment resulting in pollution of the ancient lithosphere with younger recycled material, e.g., recycled sediment. Low-degree melts from subduction-modified lithosphere are known to be significantly enriched in Pb, Sr, Nd, and Hf (Hanan et al., 2008; Jean et al., 2014) relative to sublithospheric melts by factors of tens to hundreds. Multi-component mixing involving sediment with relatively high Pb (Sr, Nd, Hf) abundances can produce linear-like arrays in Pb isotope space, similar to the overall late Miocene–Quaternary array, that resemble binary mixing trends (Douglass and Schilling, 2000; Hanan et al., 2008). The isotope compositions of basalt derived by mixing between the MORB-like sublithosphere and Anatolian continental lithosphere components would be dominated by the continental lithosphere contribution because of its greater mass proportion of Pb, Sr, Nd, and Hf. This mixing scenario explains the origin of the Quaternary trend, with positive Δ7/4Pb and Δ8/4Pb subparallel to the MORB array, and the trends in ⁸⁷Sr/⁸⁶Sr and ε_Hf versus ε_Nd space extending between the sublithospheric ambient mantle (represented by the Western Anatolia early Miocene lavas) and oceanic sediment, including northeastern Mediterranean sediments and the Ulukişla ultrapotassic lavas derived from continental lithosphere. A similar argument can be made for the Miocene trend. For these samples, the heterogeneity in the continent-derived components is much less than in the Quaternary, so the trend of Miocene samples clusters about a single tie line between ambient mantle values and label 4 on the Pb isotopic model (Fig. 6).

The mixing model (Fig. 6) shows that the Miocene lavas have source contributions of continent-derived material ranging from <1% to <10%. The late Miocene to Quaternary Anatolian lavas, most of which are from Central Anatolia, can be modeled with contributions from lithospheric source domains of ∼1%–3%. While the relationships are most pronounced in Pb-Pb plots, ε_Hf versus ε_Nd and ⁸⁷Sr/⁸⁶Sr versus ε_Nd isotope plots are consistent with the Pb model (Fig. S1 [^{footnote 1}]).

We note that the Miocene and Quaternary mixing trends defined in Pb isotope space by the Anatolian samples are similar to those observed in basalts from the Gulf of Aden, including the Sheba Ridge and its extension into the Gulf of Tadjoura (Schilling et al., 1992; Rooney et al., 2011), and suggest that these suites have similar source components. The Sheba Ridge data parallel the NHRL, and their isotope signature is dominated by mixing between the upper-mantle MORB source and a young C-like recycled mantle component (Schilling et al., 1992; Hanan and Graham, 1996; Rooney et al., 2011). The data field for the Sheba Ridge westward extension into the Gulf of Tadjoura (Asal Rift, Tadjoura Ridge, Gulf of Tadjoura) extends away from the Sheba Ridge in Pb isotope space toward high Δ7/4Pb and Δ8/4Pb, indicative of significant influence from ancient SCLM. The older Anatolian Miocene basalts are similar in isotope composition to the Asal Rift, Tadjoura Ridge, and Gulf of Tadjoura basalts due to the influence of older Pb presumably derived from the SCLM, whereas the Western Anatolian late Miocene and the Quaternary basalts resemble those of the Sheba Ridge in being more influenced by younger recycled Pb isotope signatures.

Evidence from the Pb isotopes requires that mantle source regions of the Anatolian lavas have experienced pollution involving ancient crust and/or sediment and subduction products (e.g., fluids and/or melts). Specifically, early to mid-Miocene lavas from Western Anatolia all have Pb isotope compositions with high positive Δ7/4Pb and Δ8/4Pb that touch upon or overlap the field for eastern Mediterranean sediments and Ulukişla ultrapotassic rocks (Figs. 3, 5). We note that sediments are enriched in Pb, Sr, Nd, and Hf over asthenospheric source components. The Late Cretaceous to middle Eocene Ulukişla suite is interpreted to be derived from a continental lithosphere mantle source that was metasomatized by fluids and/or melts generated during subduction (Alpaslan et al., 2004, 2006). We suggest here that the same geochemical characteristics could reflect incorporation of sediment and/or mélange into the sublithospheric source region (Behn et al., 2011; Nielsen and Marschall, 2017).

Support for this interpretation comes from the incompatible trace element abundances of early to mid-Miocene Western Anatolian lavas. When normalized to global oceanic sediment (GLOSS II; Plank, 2014), these samples exhibit essentially flat patterns anchored at unity, indicating the extent to which sediment contributions dominate the trace element signatures; the Ulukişla lavas show this same feature (Fig. 7). The trace element signatures of all other samples from Central Anatolia and the Quaternary Kula suite from Western Anatolia also have substantial overlap with GLOSS II but show differential enrichments and depletions that are likely inherited from sublithospheric sources. Specifically, both suites are depleted in large ion lithophile elements (including Pb) relative to GLOSS II, and the Kula lavas are enriched in Nb and Ta (Fig. 7). These features indicate the significance of sediment inputs and confirm the change in the source region(s) that we identified on the basis of Pb isotopic data (Fig. 5).

Contributions from sediments and subducted mélange material can be explored further through the abundances and isotopic signatures of Hf and Nd. These elements are effectively immobile in hydrous fluids (e.g., Kessel et al., 2005; Handley et al., 2011), and their similar compatibility in mantle phases means they are unlikely to fractionate during mantle melting or low degrees of crystallization in mafic melts. As a result, Hf-Nd fractionation is most likely to occur during melting of sediments or mélange (we use the term mélange following Nielsen and Marschall [2017] to indicate a physical mixture of sediments, altered oceanic crust, and hydrated mantle). Hf isotope values in oceanic sediment also vary independently from the Nd isotope values. This variation results from the “zircon effect” (Patchett et al., 1984), where sediments deposited close to their source have very low ε_Hf reflecting their detrital old zircon content; ocean sediments deposited further from the continent lack zircon and have relatively high ε_Hf controlled by the fine clay–rich continental material. Both Central and Western Anatolian lavas trend toward the high-ε_Hf field of oceanic sediment (Fig. 3). Note that the Central Anatolian late Miocene through Quaternary lavas trend to higher ε_Hf relative to contemporaneous Western Anatolian lavas, consistent with the idea that in the Central Anatolia region, where the Pb-Pb trends are subparallel to the MORB array, the mantle source for younger magmatism has experienced more subduction modification through incorporation of recycled oceanic sediments or sediment-derived fluids.

Variation in Nd isotopes as a function of Hf/Nd (Fig. 8) further delineates these mixing processes. Central Anatolian lavas define a broad mixing trend between mantle values and those of global and Mediterranean pelagic sediments (Nielsen and Marschall, 2017). This relationship indicates a role for direct contributions from sedimentary materials into the mantle source region. In detail, lavas from the Miocene and younger volcanic systems present a series of near-horizontal arrays of varying Hf/Nd at near-constant Nd isotopic compositions, interpreted as representing variable-degree melts of sediment and/or mélange + mantle (i.e., mantle that has been contaminated by mixing prior to melting; Nielsen and Marschall, 2017). This relationship cannot reflect incorporation of slab-derived fluids in which both Hf and Nd are immobile. Compositions of Ulukişla ultrapotassic lavas approach those of pelagic sediments in this representation, consistent with their Sr-Nd-Pb isotopic systematics which suggest they represent infiltration of a sediment-derived melt into the lithosphere prior to magma formation. Western Anatolian lavas show consistent patterns that correspond to the age of magmatism. Early- to mid-Miocene lavas define distinct arrays that approach the Ulukişla and/or pelagic sediment values and extend to higher Hf/Nd reflecting variable degrees of mélange melting. Late Miocene lavas define a tight array, with compositions that lack significant contribution from sediment and/or mélange. Quaternary Kula basalts from Western Anatolia plot within and adjacent to the Central Anatolia array at intermediate values. These geographic and temporal relationships are consistent with the mixing inferred from Pb isotope relationships, in which a heterogeneous sublithospheric mantle interacts with much more heterogeneous contributions from a mixed lithospheric and/or mélange source domain with variable contributions from bulk sediments.

We emphasize that the proportion of sediments and/or mélange contributing to Anatolian magmatism is likely to be quite small (e.g., Klaver et al., 2015). Ancient continental lithosphere, rejuvenated by subduction processes, can contain high concentrations of Pb, Sr, Nd, and Hf relative to the ambient MORB sublithosphere mantle source (e.g., the Pb concentration in subduction-related fluids and metasomatized lithosphere can be tens to hundreds of times higher than in depleted mantle).

The picture that emerges from the geochemical signatures is one in which an abrupt change in source material occurs across all of Anatolia in the latest Miocene. Western Anatolia Miocene lavas were derived from a mixed source of subduction modified lithosphere and/or sediment and/or mélange and a MORB-like source region with ²⁰⁶Pb/²⁰⁴Pb >19. The ²⁰⁶Pb/²⁰⁴Pb of ∼19 in the MORB-like late Miocene lavas suggests recent subduction enrichment in U, given that these samples plot near the NHRL. We interpret this to be due to the enhanced radiogenic growth of ²⁰⁶Pb from ²³⁸U decay relative to the decay ²³⁵U, the parent of ²⁰⁷Pb, which has a much shorter half-life. In contrast, latest Miocene through Quaternary lavas erupted in Western, Central, and Eastern Anatolia sample a distinct sublithospheric source that is relatively more isotopically heterogeneous due to subduction-process pollution by a young sediment component with radiogenic Pb isotopes like that of modern oceanic sediment. This is evidenced by the more-radiogenic ²⁰⁶Pb/²⁰⁴Pb and ε_Nd-ε_Hf values of the Central Anatolian trend relative to that of young Western Anatolian lavas. This shift marks a profound change in the mantle source regions for mafic volcanism.

The onset of slab rollback in the mid-Miocene appears to have been associated with an abrupt change in the isotopic characteristics of mafic lavas (Seyitoğlu et al., 1997; Aldanmaz et al., 2000; Agostini et al., 2007; Chakrabarti et al., 2012; Ersoy and Palmer, 2013). Late Miocene Western Anatolian mafic lavas with high-precision Pb isotopic data (Figs. 3C, 3D) record the greatest contribution from the sublithosphere, presumably reflecting the flow of asthenospheric material following slab rollback, accompanied by rotation or tearing, in this region (Innocenti et al., 2005; Dilek and Altunkaynak, 2009; van Hinsbergen and Schmid, 2012; Ersoy and Palmer, 2013). Similarly, Agostini et al. (2008) found Li and B isotope variations in Biga Peninsula basalts that define an abrupt change from subduction-related to rift-related signatures prior to eruption of the Kula basalts, i.e., in the latest Miocene, a scenario that corresponds well to our interpretation of Sr-Nd-Pb-Hf isotopic data. The occurrence of 14–17 Ma post-collisional lamproites and ultrapotassic mafic lavas near the Kula eruptive area (Innocenti et al., 2005; Ersoy and Palmer, 2013) marks a turning point in the tectono-magmatic evolution of Western Anatolia, given that it was followed by mafic magmatism derived from source domains with increasing representation of less-metasomatized (i.e., sublithospheric) contributions. This scenario suggests that subduction of the Aegean or Tethyan seafloor took place throughout the region until the mid-Miocene (ca. 16–17 Ma), after which time rapid lithospheric thinning released subduction-related fluids that led to the eruption of these highly alkaline lavas (Agostini et al., 2008) and allowed for ascent of hotter sublithospheric material.

Isotopic and geochemical variations among Central Anatolian mafic lavas reflect spatially distinct melting scenarios of a shared heterogeneous mantle source in late Miocene through Quaternary time. Recent tomographic images suggest the presence of a delaminated Neotethyan slab beneath this region at a depth approaching that of the mantle transition zone, and mantle shear-wave splitting results suggest asthenosphere ascends beneath the Sivas area (Gans et al., 2009; Biryol et al., 2011; Bartol and Govers, 2014). The most recent phase of volcanism, e.g., in Sivas and Erkilet as described herein, occurred ∼5 m.y. subsequent to delamination of the older slab and thus likely reflects melting of sediments and/or mélange along the top of the slab. We note that the greater depth of melting at Sivas (i.e., within the garnet stability field; Kürkçüoğlu et al., 2015) relative to more southern volcanic centers is consistent with this picture, given that the imaged depth to the top of the slab increases from south to north across Central Anatolia. We attribute the incorporation of bulk sediment and/or mélange (indicated by incompatible trace element abundances and Nd isotopic systematics; Fig. 8) as reflecting melting along the top of the descending, delaminated Tethyan slab in response to the ascent of warm asthenosphere (ambient mantle) through a tear in the subducted slab as imaged beneath Central Anatolia (Gans et al., 2009; Biryol et al., 2011). Farther south at Hasandağ, the Aegean slab displays multiple tears and fractures through and around which upwelling circum-Mediterranean asthenospheric mantle rises; heating along the base of the subcontinental lithosphere supports delamination melting at shallow levels (Reid et al., 2017; Gall et al., 2021) and may also affect sediments remaining along the delaminated slab itself at greater depth. Tectonic extension in the overlying plate likely also causes the asthenospheric material to rise to shallower depths and melt upon decompression.

New Sr-Nd-Pb-Hf isotopic results for mildly alkaline Anatolian mafic lavas reveal consistent spatial and temporal variations in mantle source mixing between lithospheric and sublithospheric source domains. The Pb isotopic signatures of Anatolian mafic lavas indicate a fundamental role for lithosphere that has interacted with mélange and/or sediments in the genesis of Pleistocene–recent volcanism. Despite the presumably volumetrically small contribution of sediment-derived melts and fluids, this material exerts a disproportionate control on the Pb (and Nd, Sr, and Hf) isotopic composition of the mafic silicate lavas. In Central Anatolia, latest Miocene through Quaternary lavas from Sivas, Develidağ, Erkilet, and Erciyes plot together with Quaternary lavas from Western and Eastern Anatolia in a region of Pb isotope space that is distinct from that of Miocene mafic lavas from Western Anatolia. The different isotopic groupings correspond to distinct mantle source regions, and this abrupt shift in contributing domains reflects a fundamental change in the geodynamics of the area.

We note that the isotope relationships and geochemical data support geophysical and tomographic models of the evolving upper-mantle structure beneath Anatolia. Delamination of the Neotethyan slab is inferred to have taken place in the early to mid-Miocene, and interaction between a sublithospheric component that represents the regional subduction modified ambient upper-mantle MORB source and slab-related materials such as sediments and mélange drives the isotopic composition of mafic lavas erupted across the region. One important implication of this result is that the subduction modification of mantle source regions contributing to Anatolian volcanism is not related solely to Miocene closure of Africa and Eurasia, but rather includes as well an older feature represented in sediments that melt along the top of the delaminated slab. Sampling of these domains occurs as upwelling sublithospheric subduction-modified MORB source material moves around and through the broken Aegean slab that dips steeply beneath southernmost Anatolia.

This work was supported by grants from the U.S. National Science Foundation (EAR-0738963, DRL-0962792 to T.F.) and TUBİTAK (107Y215 to B.K.), and by the Pennsylvania State University Department of Geosciences Hiroshi and Koya Ohmoto Fellowship in Geosciences, Paul D. Krynine Memorial Fund, and Chevron Research Travel Award to M.P.S. B.B.H. gratefully acknowledges funding for the SDSU Mass Spectrometry Laboratory that has been provided by the MRI and OCE divisions of the U.S. National Science Foundation and the W.M. Keck Foundation. Comments from L. Farmer, A. Hofmann, J. Konter, and several anonymous reviewers helped to greatly clarify our thinking around the role of sediment in subduction zones.