Geochemical and Isotopic Evolution of Late Oligocene Magmatism in Quchan, NE Iran

Magmatic activity that accompanied the collision between Arabia and Eurasia at ∼27 Ma, provides unique opportunities for understanding the triggers and magma reservoirs for collisional magmatism and its different styles in magmatic fronts and back‐arcs. We present new ages and geochemical‐isotopic results for magmatic rocks that formed during the collision between Arabia and Eurasia in NE Iran, which was a back‐arc region to the main magmatic arcs of Iran. Our new zircon U‐Pb ages indicate that collisional magmatism began at ∼24 Ma in the NE Iran back‐arc, although magmatism in this area started in the Late Cretaceous time and continued until the Pleistocene. The collisional igneous rocks are characteristically bimodal, and basaltic‐andesitic and dacitic‐rhyolitic components show significant isotopic differences; εNd(t) = +4.4 to +7.4 and εHf(t) = +5.4 to +9.5 for mafic rocks and εNd(t) = +0.2 to +8.4 and εHf(t) = +3.4 to +12.3 for silicic rocks. The isotopic values and modeling suggest that fractional crystallization and assimilation‐fractional crystallization played important roles in the genesis of felsic rocks in the NE Iran collisional zone. Trace element and isotopic modeling further emphasize that the main triggers of the magmatism in NE Iran comprise a depleted to the enriched mantle and the Cadomian continental crust of Iran. Our results also emphasize the temporal magmatic variations in the NE Iran back‐arc from Late Cretaceous to Pleistocene.

slab and continental crust (Miller et al., 1994;Straub et al., 2010). The relative importance of these contributions is assumed to vary during the lifespan of the continental arcs and can also change after subduction ceases and magmatism is transferred to collisional and postcollisional regimes. In collisional regimes, active subduction has ceased but pre-existing subduction can leave a metasomatized, young sub-continental lithospheric mantle, such as the E-MORB-like mantle that may have been the source of postcollisional magmas in Central Anatolia (Reid et al., 2017). The metasomatized mantle in both subduction-and collision-related settings can be melted to produce a series of magmatic rocks with high concentrations of incompatible elements along with depletion in Nb-Ta, such as shoshonites and/or ultrapotassic rocks (Li et al., 2002;Liu et al., 2014;Tian et al., 2017) and even magmatic rocks without depletion in Nb-Ta such as the igneous rocks from the Argentinian back-arc region (in contrast to the Chilean volcanic front) of the southern Andean magmatic system (e.g., Jacques et al., 2013Jacques et al., , 2014. Furthermore, the role of the continental crust will increase in arc magmatism in collisional settings with a thickened lithosphere. The role of refertilized lithospheric mantle, and especially continental crust, is expected to increase from the magmatic front to the back-arcs and/or rear-arcs, as observed in central to South America (e.g., Aragon et al., 2013;Sadofsky et al., 2007) and/or Tibet (e.g., Alexander et al., 2019).
There are several ambiguities regarding magmatism in collisional settings such as: (a) What is the relative importance of each compositional reservoir, including mantle and crust, on collisional magmatism (e.g., Chung et al., 2005;Huang et al., 2019;Seghedi et al., 2019)? (b) Do geochemical and isotopic signatures of continental crust overprint the mantle-derived magmas (e.g., Davidson, 1996;Jahn et al., 1999)? Answering these questions is crucial for understanding the role of continental collision in defining the geochemical and isotopic characteristics of erupting magmas. Collision-related igneous rocks can have variable geochemical signatures, ranging from adakitic (high Sr/Y, La/Yb rocks) to shoshonitic and ultrapotassic and rarely, to Oceanic Island Basalt (OIB)-type alkaline rocks (e.g., Chung et al., 2005). The continental collision can also produce bimodal magmatic suites, in which silicic rocks are usually accompanied by basaltic rocks (e.g., Gómez-Tuena et al., 2018;Kelemen & Behn, 2016;Maunder et al., 2016;Mayen et al., 2017).
Recent advances in numerical modeling support the generation of silicic magmas by assimilation and fractional crystallization (AFC) accompanying crustal melting during the intrusion of basaltic melts from the upper mantle, particularly in regions with thickened continental crust (Annen et al., 2006a;Meade et al., 2014;Patchett, 1980). The AFC processes are considered to be the main mechanisms for the formation of silicic magmas with adakitic (high Sr/Y and La/Yb) signatures (Stern & Kilian, 1996;J.-F. Xu et al., 2002; W.-C. Xu et al., 2015). Assimilation and fractional crystallization processes occur at deep crustal levels (∼20-30 km deep)-mostly associated with large intrusions of intermediate to felsic composition in the lower crust; in "deep crustal hot zones" (Annen et al., 2006a;Wan et al., 2018) and/or in mafic zones (Walker et al., 2015), where differentiating mantle-derived magmas interact with pre-existing crust, in MASH (combined Mixing, Assimilation, Storage, and Homogenization) zones.
The close association of the similarly aged basaltic-intermediate and felsic rocks-such as the magmatic calc-alkaline to adakitic rocks recorded from south Quchan in NE Iran-can be used to investigate the role of FC and AFC processes in the formation of felsic and/or adakitic rocks in a collisional setting. The basaltic lavas can be further used to unravel the mantle source of the magmatic precursors. NE Iran is a back-arc magmatic belt that contains Late Cretaceous to Cenozoic magmatic rocks. Silicic magmatic rocks from south Quchan (Figure 1a) occupy an area of ∼2,000-3,000 km 2 in NE Iran. These rocks are suggested to show a long-lived magmatic activity from Oligocene to Pleistocene (K-Ar ages; Bauman et al., 1983;Shabanian et al., 2012) and are associated with Late Oligocene (24-23 Ma) basaltic-andesitic magmatism in the region. Mafic-intermediate volcanic rocks have been geochemically classified as high-Mg andesites to basalts as well as high-Nb hawaiites and mugearites (Ahmadi et al., 2017), whereas silicic rocks display adakitic signatures (Shabanian et al., 2012). Both basaltic-andesitic and felsic magmatism have occurred in a subduction-related context from Oligocene to Pleistocene time but after the collision between Iran and Arabia that commenced at ∼27 Ma (McQuarrie & van Hinsbergen, 2013). Therefore, these rocks present a "collisional" perspective, and we prefer to use the term "collision-related" (vs. "postcollisional") for describing these rocks, as the collision between Arabia and Iran is still ongoing and these rocks do not occur in a "relaxation" period after the main collisional event ceased. In this study, we present the first systematic geochemical, isotopic and geochronological (zircon U-Pb) study of Late Oligocene (∼24 Ma) collision-related bimodal magmatism (silicic and mafic rocks) from NE Iran back-arc (south Quchan). The primary goals of this paper are (a) to characterize-based on their geochemical and isotopic signatures-the sources and origin of collision-related mafic magmatic rocks from NE Iran; (b) to understand the role and relationship between fractional crystallization (FC), AFC and recharge, evacuation, and fractional crystallization (REFC) in the genesis of felsic and/or adakitic rocks via their mafic precursors; and (c) to establish in a broader sense the relationships between magmatism, crustal growth and the role of MASH zones in generating silicic and adakitic rocks within continental arcs. This study has the potential to bring about a step-change in assessing the causes of the geochemical-isotopic variations in collisional magmatism of NE Iran, which can be tapped from the lithospheric mantle and the Late Neoproterozoic-Early Paleozoic continental crust of Iran.

Geological Setting
The backbone of the Arabia-Eurasia collision zone, the Zagros orogen, was built through Neotethyan subduction, magmatism and collision between several continental blocks including Iran and Anatolia (Berberian & King, 1981). These continental blocks have an Ediacaran-Early Paleozoic (Cadomian) crust and have been rifted from Gondwana in Late Paleozoic time and accreted to Eurasia during Triassic-Jurassic time (Stampfli & Borel, 2002). The Cadomian magmatic rocks exposed in several segments of Iran ( Figure 1a) have ages of 620-500 Ma . The rifting of Gondwana in the Late Paleozoic (Permian) caused the opening of the Neotethyan ocean in western parts of the Iranian blocks and the consumption of Paleotethys in the northeastern parts (Berberian & King, 1981). The collision between the Iranian continental blocks and Eurasia during Triassic time led to the closure of the Paleotethyan ocean. The subduction of the Paleotethyan oceanic lithosphere beneath Eurasia (Turan platform) and then collision between the Iranian continental blocks and Eurasia are marked by the presence of a Late Paleozoic suture zone (Late Paleozoic ophiolites through Mashhad to Fariman), Carboniferous-Permian arc magmatism and Late Triassic collisional I-and S-type granites in NE Iran (Alavi, 1991;Mirnejad et al., 2013;Zanchetta et al., 2013). There was a passive margin along western Iran from Late Triassic to Early Cretaceous time, with accommodation of Jurassic continental and Early Cretaceous marine sediments leading to collapse of this margin (Azizi & Stern, 2020;Azizi et al., 2018), and then subduction initiation in Late Cretaceous. The subduction initiation along western Iran (Main Zagros Thrust, Figure 1) led to the opening of several back-arc basins such as the Late Cretaceous Sabzevar-Torbat-e-Heydarieh basin.
Igneous activity in Iran has been continuous from Late Cretaceous (∼110 Ma) up to the Pleistocene as a result of the subduction of the Neotethyan oceanic crust beneath Iran from Late Cretaceous to Oligocene and then through the collision between Arabia and Iran from ∼27 Ma onward (McQuarrie & van Hinsbergen, 2013). The subduction was accompanied by the development of a continental (Andean-type) magmatic arc, parallel to but ∼150-200 km from the Zagros suture zone (Main Zagros Thrust), known as the Urumieh-Dokhtar Magmatic Belt (UDMB) (Figure 1a). The UDMB thus outlines the magmatic front (MF) of the Neotethyan subduction. Magmatic rocks are also abundant in northwestern to northeastern parts of the UDMB and are considered as the products of a back-arc (BA) magmatism. Previous geochronological data suggest the UDMB experienced magmatic flare-ups in the Eocene (∼54 Ma until 37 Ma, subduction-related magmatism) and the Miocene (20-5 Ma, collision-related magmatism) (Chiu et al., 2013;Verdel et al., 2007). However, recent data show that the high magmatic fluxes did not occur at the same time in the magmatic front and back-arcs. The UDMB (the magmatic front) magmatism seems to be subdivided into two distinct episodes: ∼80-70 Ma and ∼50 Ma present, with protracted magmatic activity lasting from 40 Ma to the present (Sepidbar et al., 2019), while the back-arc magmatic flare-ups occurred at 110-80 and 75-35 Ma in the NE Iran back-arc and 55-35 Ma in the NW Iran back-arc .
Below we subdivide the Cenozoic Iran arc into the magmatic front and the NE Iran back-arc and discuss these sub-provinces separately.

Magmatic Front
The UDMB is a 50-80 km wide volcano-plutonic belt that defines the magmatic front and trends NW-SE for >1,000 km across Iran between 28° and 39°N ( Figure 1). The UDMB evolved for ∼80 m.yr. from when subduction began underneath the Iranian Plateau in Late Cretaceous time, continuing as a mature arc in Paleogene with a magmatic flare-up during the Eocene (Verdel et al., 2011), before swapping magmatic style into collisional high magmatic fluxes following the collision with Arabia during Late Oligocene (∼27 Ma). The UDMB arc started with the eruption of low-K tholeiitic and calc-alkaline magmas. A magmatic lull occurred during the Early Paleocene time, accompanying the uplift of the Zagros forearc ophiolites (Moghadam & Stern, 2011). The uplift is evidenced by the presence of lower Paleocene conglomerates, rich in ophiolitic fragments, which cover the ophiolites and are also present in the UDMB. The magmatism restarted in the Middle-Upper Paleocene with the eruption of pyroclastic rocks and lavas. Middle to Late Paleocene igneous rocks have calc-alkaline characteristics and most of them show mantle-like juvenile isotopic signatures . Eocene igneous rocks are abundant and dominantly high-K calc-alkaline. Late Miocene to Pleistocene UDMB igneous rocks have adakitic and ultrapotassic geochemical signatures (Pang et al., 2015(Pang et al., , 2016. Miocene adakitic lavas are mainly present in the UDMB and are interpreted as products of melting of the eclogitic parts of a thickened lower continental crust in a postcollisional setting (e.g., Deng et al., 2018;Wan et al., 2018).

Back-Arc Magmatism in NE Iran
Back-arc igneous activity occurred along an NW-NE arcuate belt >1,200 km long (Figure 1a) during the latest Cretaceous to Pleistocene, with the eruption of marine to subaerial magmatic rocks, including thick sequences of marine pyroclastic rocks (Ballato et al., 2011;Verdel et al., 2011). Magmatism in the NE Iran back-arc started in the latest Cretaceous and continued through the Paleogene-Neogene to the Pleistocene. The Late Cretaceous to Neogene magmatic rocks are distributed north and south of the Sabzevar-Torbat-e-Heydarieh ophiolitic belt, but also intrude the ophiolites to the north and south (Figure 1b). The Sabzevar-Torbat-e-Heydarieh ophiolites are considered as fragments of a Late Cretaceous back-arc oceanic basin associated with subduction initiation along what is now the Zagros orogen. Late Cretaceous igneous rocks with zircon U-Pb ages of 104-92 Ma are abundant in the NE Iran back-arc (Figure 1b) (Alaminia et al., 2013;Kazemi et al., 2019). The Late Cretaceous sequence is largely marine, comprising low to high-K calc-alkaline rhyolitic to dacitic-andesitic lavas, felsic tuffs and radiolarites. Subaerial lavas including alternating basalts-andesites and dacites are also common. These are overlain by both Middle Paleocene-Eocene terrigenous sediments (Oryan basin; Figure 1b) and Mid-Late Paleocene-Eocene acidic to intermediate pyroclastic rocks with intercalated basaltic to andesitic lavas. In the NE back-arc, Late Cretaceous intrusions (104-76 Ma) were emplaced into thick sequences of Cretaceous terrigenous sediments as well as pyroclastic and volcanic rocks (Mazhari et al., 2019).
A thick pile of magmatic rocks in the NE Iran back-arc occurred after the emplacement and exhumation of the Sabzevar-Torbat-e-Heydarieh ophiolites in the Latest Cretaceous-Middle Paleocene time. The postophiolite subduction-related magmatism started in the Mid-Late Paleocene with the eruption of pyroclastic rocks and interlayered lavas and continued to the Early Eocene with the outpouring of basaltic to andesitic melts. Intermediate to felsic plutonic rocks with Latest Cretaceous-Paleocene (66-58 Ma) ages also occur in NE Iran back-arc ( Figure 1b). Latest Cretaceous-Paleocene plutonic rocks have low-K tholeiitic to calc-alkaline geochemical signatures, and their radiogenic bulk rock Nd and zircon Hf isotope characteristics suggest derivation from a depleted mantle . Eocene volcanic and plutonic rocks (47-40 Ma) are mostly high-K calc-alkaline, but adakitic rocks are also common. These rocks show traces of AFC processes within the middle-upper crust .
Collisional magmatism in the NE Iran back-arc began at 27 Ma and continued to the Late Miocene (∼10 Ma) and even to the Pleistocene (Kheirkhah et al., 2015). Most of the collision-related mafic to felsic volcanic rocks in the back-arc region of NE Iran occur south of Quchan and also in W-NW Sabzevar (Rostami-Hossouri 10.1029/2021GC009973 6 of 40 et al., 2020) ( Figure 1a). Pleistocene subaerial andesites also occur north Torbat-e-Heydarieh. The felsic rocks in south Quchan appear both as volcanic domes 1-4 km wide and as lava flows with K-Ar ages of 19.5 ± 0.5 to 2.3 ± 0.1 (Ghasemi et al., 2010). The felsic rocks are accompanied by andesitic to basaltic lava flows. The volcanic domes are intruded into Eocene, but occasionally Miocene, volcano-sedimentary sequences including sandstones, marls, tuffites and shales ( Figure 2). Lava flows rest on Miocene red shales and marls. There are xenoliths of sedimentary hosts within the felsic rocks, which suggest assimilation of the crustal rocks and contamination of the magma. Volcanic rocks are also covered by Pliocene-Pleistocene un-sorted conglomerates. Geochemically the felsic rocks have adakitic signatures and were generated during local tectonic activity in NE Iran-Kopet Dagh (Bauman et al., 1983;Shabanian et al., 2012). Mafic rocks were suggested to include high-Mg andesites to basaltic as well as high-Nb hawaiites and mugearites and show 40 Ar-39 Ar ages of 23.1 ± 0.3 and 22.9 ± 0.5 Ma (Ahmadi et al., 2017). Dacitic lava flows yielded ages of 24.1 ± 0.4 and 23.2 ± 0.5 Ma (Ahmadi et al., 2017).

Analytical Procedures
This paper deals with the Late Oligocene volcanic (both mafic and felsic) rocks from south Quchan. We also compare our results with compiled data on the Late Cretaceous-Pleistocene magmatic rocks of NE Iran . We have used six main analytical procedures: (a) JEOL wavelength-dispersive electron probe X-ray micro-analyzer (JXA 8800R) to determine the composition of minerals at Institute of Geology and Geophysics (IGG-CAS), China; (b) X-Ray Fluorescence (XRF) and Inductively Coupled Plasma-Mass Spectrometry (ICP-MS) for whole-rock major-and trace-element analyses at Institute of Geosciences, Kiel University, Germany; (c) Cathodoluminescence imaging of zircons at IGG-CAS; (d) Secondary-ion microprobe (SIMS) analysis for U-Pb zircon ages at IGG-CAS; (e) Multi-Collector Inductively Coupled Plasma-Mass Spectrometry (MC-ICP-MS) to analyze zircon Lu-Hf isotopes at IGG-CAS; (f) The Sr-Nd-Hf-Pb whole-rock isotope composition of the bulk rocks were analyzed using a MC-ICP-MS at University of Bonn/Köln. We have studied >70 samples petrographically; 26 were used for whole-rock chemical analysis, 4 for SIMS U-Pb zircon ages and zircon Lu-Hf analysis and 21 for Sr-Nd-Pb-Hf isotopes. We have selected fresh samples for whole-rock geochemistry. We sampled different geographic locations to cover all rock units from all parts of the studied area. Analytical details are described in Supporting Information S1.

Sample Descriptions
Our Late Oligocene samples from NE Iran include basalts, olivine basalts, andesites and dacites-rhyolites ( Figure 2). Olivine basalts contain clinopyroxene phenocrysts and large crystals of olivine (∼2-3 mm), with interstitial plagioclase laths (Figures S1a and S1b), whereas basalts comprise clinopyroxene and plagioclase with fine-grained olivine in the groundmass. The groundmass of the investigated samples is dominated by clinopyroxene + plagioclase lath + Fe-Ti oxides + glass ± olivine. Olivine in the basaltic rocks is fresh in most cases and contains inclusions of alumina-rich spinel. Andesitic rocks contain more plagioclase and clinopyroxene microphenocrysts than basalts, without olivine. Iron-titanium oxides occur in small amounts (<2 vol.%) in andesites. These mafic rocks show holocrystalline to intergranular texture. Some andesites contain plagioclase + sanidine microlites and clinopyroxene microphenocrysts and show trachytic textures. Amphibole is absent in olivine basalts but occurs as rare needles in the groundmass of andesites and rarely in basalts.
Most felsic rocks have dacite and rhyolitic composition. Dacites contain plagioclase and sanidine phenocrysts with rounded quartz. Modal abundances of quartz and sanidine phenocrysts are higher in the rhyolites than in the dacites. Plagioclases are zoned and show sieve texture in the cores (Figures S1c and S1d). Zoned amphiboles with resorbed rims are also present in dacites (Figures S1e and S1f). The groundmass of dacites contains glass, plagioclase microlites, sanidine and quartz. Dacites and rhyolites show hyaloporphyritic and microlitic to cryptocrystalline textures.

Mineral Geochemistry
We have analyzed olivine, pyroxene, plagioclase, amphibole and oxides from Late Oligocene basalts and dacites. Olivine phenocrysts in basalts have forsterite contents of 69.7%-85%. Olivine shows zonation; the cores 10.1029/2021GC009973 7 of 40 of olivine crystals have higher forsterite content (Fo ∼85%) than their rims (Fo ∼69.7%). Their NiO content varies from 0.04 to 0.15 wt%. Clinopyroxene in basalts is diopside (Figure 3a). The cores of the plagioclase phenocrysts in basalts have compositions from bytownite to labradorite (An = 50.7-81.5) (Figure 3b), except for some grains which show andesine (An = 47.1-48.4) composition toward the rims. Some plagioclase from sample MB12-1 also shows more sodic compositions toward the rims (An = 17.6-26.1) with high K 2 O Amphiboles also show oscillatory zonation and Al 2 O 3 content varies in different bands ( Figure S2). There seem to be two different amphiboles in dacites; one group has high-Mg cores, while the second type shows low-Mg cores; both show oscillatory zoning toward the rim ( Figure S2). Oxide micro-phenocrysts in basalts are compositionally similar to titanohematite. Spinel in basalts occurs as inclusions in olivine phenocrysts and has high Al contents (Cr#; 0.21-0.23). In a TiO 2 versus Al 2 O 3 diagram, the spinels have similarities to spinel from mid-oceanic ridge basalts (MORB) (Figure 3e).

Zircon U-Pb Ages
We have dated four samples from dacitic domes and lavas from the NE Iran back-arc. Sample CG12-2 is from a dacitic dome. Zircons from this sample are medium-to coarse-grained (∼80-100 μm) and show oscillatory zoning. Seventeen analyses from this sample show an intercept age of 23.7 ± 0.4 Ma (Figure 4). Two inherited cores give ages of 98 and 224 Ma. These inherited cores show the assimilation of older rocks by the Quchan Oligocene magmas. The common-lead content of the zircons is f 206 < 2.4%.   (Morimoto, 1988) for clinopyroxenes in Quchan basalts. (b) Triangular plot for classification of the Quchan feldspars (Deer et al., 1992). (c) Mg/Mg + Fe +2 versus Si cation binary plot (Leake, Woolley, Arps, et al., 1997) for the classification of Quchan amphiboles. (d) Triangular TiO 2 -FeO-Fe 2 O 3 plot for the classification of Quchan Fe-Ti oxides (Deer et al., 1992). (e) TiO 2 versus Al 2 O 3 diagram (Kamenetsky et al., 2001) for discrimination of spinels from Quchan basalts. (f) Compositional profiles for anorthite content in plagioclase (upper panel) and Al 2 O 3 in amphibole (lower panel) from Quchan dacites.
Sample ZO12-1 is taken from a dacitic lava flow. Zircons from this sample are coarse-grained (∼100-150 μm). Most zircons show oscillatory zonation, although some of them are weakly zoned. Fifteen analyses were obtained from this sample. Common-lead content is low in most grains; with f 206 < 1.5%. Some grains have a higher common lead with f 206 = 2.6-5.4. Zircons from this sample show an intercept age of 23.7 ± 0.3 Ma ( Figure 4).
Sample KN12-4 is from a dacitic dome. Zircons from this sample are medium-to coarse-grained (∼80-130 μm) with oscillatory zoning, although some zircons have only weak oscillatory zonation. Seventeen analyses from magmatic zircons show an intercept age of 23.8 ± 0.2 Ma ( Figure 4). The common-lead content of magmatic zircons from this sample is f 206 < 2.8.
Sample SM12-2 is taken from a dacitic dome. Fifteen analyzed magmatic zircons with oscillatory zonation show an intercept age of 24.5 ± 0.3 Ma, which is interpreted to be the crystallization age of this sample. The zircons have a common lead with f 206 < 5.4.
The new ages obtained in this study range from ∼25 to 24 Ma (Late Oligocene), which is in accordance with the reported 40 Ar-39 Ar ages (23-24 Ma; Ahmadi et al., 2017) for mafic and felsic rocks. However, some felsic volcanic rocks/domes seem to have younger K-Ar ages of 19.5 ± 0.5 to 2.3 ± 0.1 (Ghasemi et al., 2010).
We have normalized the rare earth-and trace-elements data from Quchan rocks to chondrite and N-MORB.

Radiogenic Isotopes
We analyzed Late Oligocene volcanic rocks for Sr-Nd-Pb-Hf isotopes ( good negative correlation, although two non-adakitic samples (dacite sample CG12-1 and rhyolite sample AG12-1) have higher 87 Sr/ 86 Sr (t) at a given 143 Nd/ 144 Nd and could reflect late-stage processes, e.g., alteration and/or assimilation of altered host rocks such as tuffites or marls and/or biotite-rich Cadomian paragneisses which are abundant in NE Iran. Assimilation of the biotite-rich Cadomian rocks could explain the radiogenic Sr-isotope compositions of some Quchan felsic rocks, but biotite does not have high U/Pb and thus cannot explain the lack of disturbance in the Pb isotopic composition of these rocks. However, sample CG12-1 has higher CaO and LOI contents compared to other dacites, whereas sample AG12-1 has higher LOI relative to other rhyolites, which would be consistent with their high 87 Sr/ 86 Sr (t) being related to alteration. Sample AD12-1 has higher 87 Sr/ 86 Sr (t) and lower εNd(t) than other dacites and rhyolites and  plots toward the Ediacaran-Early Cambrian (Cadomian) crust of Iran. All these non-adakitic samples (CG12-1, AG12-1, and AD12-1) have low Rb (2.1-10.8 ppm) and high Rb/Sr (0.006-0.03) compared to other felsic samples. Rhyolite sample AG12-1 has higher εNd(t) and εHf(t) than mafic rocks. Most Quchan rocks plot in a trend defined by Paleocene-Eocene to Pleistocene magmatic rocks of NE Iran (Figure 8a). This trend runs between a depleted mantle (DM), similar to the mantle source of arc tholeiites and the Cadomian continental crust of Iran. The εHf(t) values for basalts and andesites vary between +5.4 and +9.5. These samples fall both along and below the mantle array in a plot of εHf(t) versus εNd(t) (Figure 8b). The εHf(t) values for felsic rocks range from +3.4 to +12.3. They plot parallel to, but mostly below, the mantle array, in a trend between the depleted mantle and the Cadomian continental crust of Iran. In the thorogenic-Pb isotope plot (Figure 8c), the Quchan igneous rocks fall in a narrow linear array, well above the Northern Hemisphere Reference Line (NHRL) of (Hart, 1984), with ∆8/4 of 75-56 (Figure S3). The samples define a trend between the depleted mantle and Cadomian crust (=CC) reservoirs (Figures 8c and S3). On the uranogenic Pb diagram (Figure 8d), the investigated samples plot above the NHRL (∆7/4 to 6-12, Figure S3) and overlap with the Paleocene-Eocene igneous rocks of NE Iran. High uranogenic-Pb isotope ratios of Quchan rocks may show the involvement of an enriched mantle. This enriched mantle could be the result of contamination and/or re-fertilization of a depleted mantle source with marine sediments through the subduction of the Neotethyan oceanic lithosphere beneath Iran, creating a metasomatized lithospheric mantle. This type of enriched lithospheric mantle has also been suggested for the formation of high-K rocks from the Mediterranean region (e.g., Conticelli et al., 2009;Cvetkovic et al., 2013). A metasomatized enriched mantle with radiogenic Pb (∆8/4 > 45) is also suggested as the source of postcollisional (Pleistocene-Holocene) volcanism in Central Anatolia (Reid et al., 2017). Sample AD12-1 has higher thorogenic and uranogneic Pb isotope ratios whereas sample AG12-4 has lower values of thorogenic and uranogneic Pb. Both samples are non-adakitic dacites with low Sr/Y and  La (n) /Yb (n) . The radiogenic Pb isotopes of some felsic rocks can also show assimilation of Cadomian continental crust during magma ascent.
Zircons from dacites and rhyolites display εHf(t) from +13.7 to −6.5. The highly variable εHf(t) values for these rocks again could attest to contamination with the Cadomian continental crust.

Discussion
The composition of magmas erupted in subduction-related continental arcs can reflect an interplay between the mantle source composition and processes occurring during melt generation and ascent to the arc surface. These processes may take place in the mantle (e.g., mixing of sediment/slab melts and fluids into a depleted mantle) or crust (e.g., MASH or assimilation-fractional crystallization which can incorporate lower-or upper continental crust during magma ascent and storage). Heterogeneities within the mantle and/or lower-crustal reservoirs can result in geochemical variations in magma composition (e.g., Rapp et al., 2008;Stracke, 2012;Willbold & Stracke, 2006. In addition to source heterogeneities, continental middle to the upper crust, especially in the case of collisional systems, can contribute to final magma compositions. Below, we first explore the characteristics of the mantle source in the Quchan area using basaltic rocks, and then we use the basaltic rocks to understand the origin and evolution of Quchan felsic and/or adakitic rocks. To unravel the genesis of the felsic/adakitic rocks, we first focus on the geochemical signatures of these rocks as well as the composition of their amphiboles and then use the FC, AFC, and REFC models to show the evolution of the major-and trace-elements from basaltic rocks into intermediate and then into felsic/adakitic rocks. Finally, we will use the Sr, Nd, and Pb isotopes to see if the isotopic signatures of the Quchan rocks are solely controlled by a depleted mantle and its interaction with crustal components via assimilation, or whether an enriched mantle source is required.

Source Characteristics
Postcollisional Oligocene magmas in NE Iran could come from the melting of a metasomatized mantle, with traces of the previous subduction, from Late Cretaceous to Oligocene, before the continental collision occurred. Therefore, several complex processes, including the interaction of subducting sediments/ sediment melts with a depleted mantle, could enrich the mantle wedge or the sub-continental lithospheric mantle.
The Quchan volcanic rocks isotopically plot in a trend between the depleted mantle and the Cadomian continental crust (Figure 8) and their isotopic and trace-element signatures rule out a depleted MORB-type mantle as a unique source for the genesis of these rocks. The presence of spinels with high TiO 2 and Al 2 O 3 in near primitive basaltic rocks (Figure 3e) could also eliminate a depleted mantle as the source of Quchan mafic rocks. Although crustal processes such as AFC could lead to the generation of felsic rocks, these processes are inadequate to cause the range of trace elements and isotopic compositions observed in basalts (see Section 5.3.4). Instead, we suggest some aspects of the wide range in trace elements and Sr, Nd, Hf and Pb isotopic composition of the Quchan volcanic rocks can be inherited from heterogeneity within the mantle. Both source heterogeneities and contamination with the Cadomian continental crust of Iran could generate such trace element and isotopic variations.
Most volcanic rocks from Quchan are fractionated silicic rocks and the primary mantle-derived basaltic rocks needed to unravel the mantle composition are rare. However, some basalt samples with high Mg# (57%-59%) and with εNd(t) = 5.4-4.4; εHf(t) = 9.5-5.4 have high Nb/U (∼9-26) and low Zr/Nb (∼9-15), which suggests that these samples could allow the extrapolation of the composition of the primary man- Figure 6. Na 2 O + K 2 O and K 2 O versus SiO 2 discrimination diagrams (a and b) for the classification of Quchan volcanic rocks (modified after Lebas et al. [1986]). La (n) /Yb (n) versus Yb (n) (c) and Sr/Y versus Y diagrams for discriminating Quchan volcanic rocks (compositional domains for adakite and normal arc rocks are according to Castillo [2012] and Defant and Drummond [1990]). Geochemical data for Late Cretaceous to Pleistocene magmatic rocks are from .
tle melts. As such, we have applied a simple mixing model and subsequent partial melting considering a Depleted Mantle (DM) that was metasomatized by bulk mixing with altered oceanic crust and terrigenous sediments (Figure 9). Our model is not a unique solution to explain the genesis of all basaltic rocks, but this model is geologically reasonable for a region that has been experienced long-time subduction from Late Cretaceous to Late Oligocene. For modeling, we have used three basaltic samples with MgO ∼5-7 wt%, that is, samples MB12-1, MB12-2, and DA12-1. In order to examine the source mantle composition and processes, we used the PRIMACALC2 spreadsheet, which allows an estimation of a primary melt in equilibrium with the upper mantle (Kimura & Ariskin, 2014). Normalization data are from Sun and McDonough (1989). Geochemical data for Late Cretaceous to Pleistocene magmatic rocks are from . The Palaeocene-Eocene geochemical data are heterogenous, and describing these heterogeneities is beyond the scope of this paper, but these data support the case for complex magmatic plumbing systems that evolve in continental arcs.
Data for the altered oceanic crust and terrigenous sediments, DM, mantle modal mineralogy and bulk Kd, come from (Iveson et al., 2018;Jacques et al., 2014;Kimura et al., 2009;Nandedkar et al., 2016;Saginor et al., 2013;Workman & Hart, 2005;Zhang et al., 2019). The bulk partition coefficients assume mantle mineralogy of 45% olivine, 19% clinopyroxene, 23% orthopyroxene, 10% amphibole and 3% garnet. Subsequently, we assumed aggregated fractional melting of the mixed source (DM + altered oceanic crust and trench sediments) with variable degrees of melting. All modeling parameters are presented in Table S6. Since both the degree of melting and ratio of DM to altered oceanic crust and trench sediments are unknown, we applied an iterative method whereby both unknowns have been changed until the composition of the modeled melt resembled the selected Quchan basalts. Our model (modeled melt (1) in Figure 9) indicates that a 90:6:4 mixture of DM and terrigenous trench sediment and altered oceanic crust, with 5% aggregated fractional melting, will closely match the trace-element abundances of our target basalts (Figure 9). Higher proportions of altered oceanic crust and trench sediments or different melting percentages give patterns that do not match our samples. However, the modeled melt has slightly lower Sr (±Nb) than the target samples, but also anomalies in Pb, U, and Rb. The modeled melt has a positive anomaly in Pb compared to our basalts. The Pb content in the subduction-related magmas is very dependent on the amount and the composition of subducting sediments. The target basalts instead have a peak in Sr and sample MB12-2 does not show depletion in Nb-Ta, which is different from the modeled Nb-Ta depleted melt. We thus emphasize that two basaltic samples with Nb-Ta depletion could come from melting of a metasomatized mantle (i.e., a mixture of DM and terrigenous trench sediment + altered oceanic crust), similar to the modeled melt. Sample MB12-2 without Nb-Ta depletion is alike to the enriched MORBs and could come from partial melting of an enriched mantle. Similar mantle sources have been considered for the formation of postcollisional magmas  (Zindler & Hart, 1986 in Central Anatolia (Reid et al., 2017) and Quaternary basaltic lavas from the Argentinean back-arc region of the southern Andean subduction system (e.g., Jacques et al., 2013Jacques et al., , 2014. This basalt sample MB12-2 has radiogenic Pb isotopic composition-with ∆7/4 = 6.98 and ∆8/4 = 60.08 ( Figure S3)-consistent with formation in a subduction-related setting.
It is possible that the altered oceanic crust and trench sediments and/or the "ambient mantle" composition we selected for our model are different from those that prevailed in the magmatic source of Quchan rocks. Unlike recent subduction systems, it is hard to evaluate the compositions of the source mantle and subducted materials for an ancient subduction system like the Cenozoic arcs of Iran. For further testing of our model, we have used the average composition of Eastern Mediterranean Sea sediments (Klaver et al., 2015), to mix with a depleted mantle. Our new model (modeled melt (2)) support that an 89:11 mixture of DM and Eastern Mediterranean Sea sediments-with 6% aggregated fractional melting-again closely match with the trace-element composition of our target basalts, excluding positive anomalies in Pb, Rb, Ba, Th, and U and less negative anomalies in Nb-Ta than our target samples. The positive anomalies in the above-mentioned elements could be related to the nature of the presumed subducting sediments (the Eastern Mediterranean Sea sediments) which have marl and mudstone composition.

Amphibole as a Precursor
Experimental studies emphasize the role of amphibole fractionation during the crystallization and geochemical evolution of subduction-related, H 2 O-rich magmas (e.g., Carmichael, 2002;Grove et al., 2002Grove et al., , 200 3, 2006Grove et al., , 2012. These studies attest to the importance of amphibole crystallization and its effect on the HREE to MREE patterns and probably it is a key for producing adakitic signatures in some magmatic rocks from continental arcs (e.g., Davidson et al., 2007;. Extreme amphibole fractionation could be considered as a mechanism for the formation of spoon-like HREE to MREE patterns in Quchan felsic samples and also can reconcile some elemental ratios such as La/Nb, Nb/U, Dy (n) /Yb (n) , La (n) /Sm (n) , and Dy/Dy* which are keys for recognizing the adakitic signature (see next section). We believe the composition of amphibole phenocrysts from Quchan felsic rocks can give a representative indicator of the magmatic processes Figure 9. Sample/primary-mantle normalized diagram for the primary Quchan melts (basaltic samples including MB12-1, MB12-2 and AD12-1) in equilibrium with mantle-calculated using PRIMACALC2 (Kimura & Ariskin, 2014). Modeled melts (melts [1 and 2]) assumes 5% and 6% aggregated fractional melting of the mixed (90% DM + 6% subducted sediments + 4% altered oceanic crust, AOC and/or 89% DM + 11% Eastern Mediterranean Sea sediments) mantle sources. For compression, we also show the modeled melt after 1% aggregated fractional melting of the mixed mantle (1). Data and parameters for these calculations are provided in Table S6. Data on Argentinian back-arc average basalt is from Jacques et al. (2013Jacques et al. ( , 2014. during the formation of adakitic-like geochemical signatures. For this purpose, the mineral formulae for amphiboles were re-calculated following the spreadsheet supplied by , and re-interpreted by . In this spreadsheet, the Fe +3 /Fe +2 ratio is re-calculated by charge balance after correcting the tetrahedral (Si, Al, and Ti) plus octahedral (Al, Ti, Cr, Fe, Mn, and Mg) cations to 13 . This calculation shows that amphiboles from Quchan felsic rocks mostly have Tschermakitic pargasite composition, with ∼0.02-0.4 apfu Na cation (atoms per formula unit = apfu) in site A, but with 1.5-1.8 Ca and 0.2-0.5 Na cations in site B and 0.08-0.2 Ti in Site C, and were in equilibrium with calc-alkaline melts. Site A also contains a noticeable amount of K (0.04-0.2).
In a plot of Al (T) (total Al = [6] Al + [4] Al) versus [4] Al, considering the Al# content (= [6] Al/ [6] Al + [4] Al) ( Figure S4), all amphiboles from Quchan felsic rocks follow a common compositional trend. The Al (T) of these amphiboles ranges between 1.6 to 2.2, with Al# between 0.12 and 0.28. Our data show that these amphiboles have higher Al# content than amphiboles from shallow volcanic rocks (with Al# = 0) and even some analyzed spots from amphiboles (toward the cores of grains) have an Al# content comparable to high-P crustal and/or mantle-derived (experimental) amphiboles (with Al# = 0.21) . Thermobarometric calculations record pressures of 4-6 Kbar and temperatures ranging from 834 to 900°C. Oxygen fugacity can be calculated from the amphibole composition  using the amphibole magnesium index (Mg*): ∆NNO = 1.644 Mg* -4.01 (where Mg* = Mg + Si/47 - [6] Al/9 -1.3 [6] Ti + Fe +3 /3.7 + Fe +2 /5.2 -B Ca/20 -A Na/2.8 + A [t]/9.5). The calculated values range from log-fO 2 = −13.7 to −10.5, and ∆NNO = −0.6 to +1.8. For comparison, the ∆NNO for calc-alkaline magmas is suggested to be in the range of −1 to +3 (Carmichael, 1991;. The relationships between temperature (834-900°C) and H 2 O melt (7-9.5; inferred from amphibole composition) indicates maximum stability for the formation and fractionation of Quchan amphiboles in the middle-lower crust (e.g., Müntener et al., 2001). Finally, these P, P H2O , and T estimates on amphiboles from different Quchan rocks could suggest that differentiation of the Quchan melts occurred at ca 14-20 km depth (assuming an average crust density of 2,700 kg/m 3 [e.g., Rossetti et al., 2017] in the middle crust), although amphiboles with high Al# (Al# >0.2) may suggest greater depths. Therefore, our data are consistent with the observation that parts of amphiboles (or amphibole cores) crystallized earlier, in the deeper crust-and probably in deep crustal hot zones. Crystallization of H 2 O-rich magmas similar to the parental melts of Quchan lavas (with H 2 O melt = 7 to 9.5) in the middle to lower crust (at >∼7 km) will lead to amphibole appearing as the first mineral on the liquidus of basaltic to andesitic melts (Carmichael, 2002;Moore & Carmichael, 1998;Rooney et al., 2011). However, amphibole fractionation in the lower-middle crust is supposed to follow the earlier crystallization of olivine and pyroxene in the upper mantle-lower crust (e.g., Rooney et al., 2011), which is also consistent with our fractional-crystallization modeling (see below).

Geochemical Evaluation of Quchan Adakites
Several models have been suggested for the origin of adakites including melting of the eclogitic parts of subducting oceanic crust and/or eclogitic segments of the lower continental crust (e.g., Castillo, 2006;Chiaradia et al., 2009;Kay & Kay, 2002;Kolb et al., 2013;Rapp & Watson, 1995). Garnet and/or amphibole fractionation along with AFC processes at the base of the thickened continental crust has been also considered as a candidate for the genesis of adakites (e.g., Chiaradia, 2009;Hidalgo & Rooney, 2014;. AFC and/or MASH processes have been suggested as the main trigger for the formation of both silicic and adakitic magmas in the lower continental crust, in deep crustal hot zones, where mantle-derived magmas are being modified (Annen et al., 2006a).
Our new geochemical data show that the Late Oligocene Quchan felsic volcanic rocks show both normal arc (±mafic volcanic rocks) and adakitic compositions. Fractionated Quchan felsic rocks show shallower slopes in plots of Th/Yb, Ba/Yb, and La/Yb versus Nb/Yb (Figures 10a, 10b, and 10e), which can reflect enrichment in Nb and depletion in Yb and Y because of crystal fractionation. The Quchan adakitic rocks have higher Nb/Yb and Th/Yb than basalts and andesites, but some non-adakitic samples (AD12-1, AG12-1, AG12-4, and CG12-1) have similar ratios to basalts and andesites. The Quchan Oligocene adakites have Nb/Yb, La/Yb, Th/Yb, and Ba/Yb similar to highly fractionated Eocene adakitic-like dacitic-rhyolitic domes (barren adakites) from NW Sabzevar and are geochemically distinct from Eocene, fertile (Cu-Au-  . For clarity, the domains of Eocene highly fractionated barren silicic rocks with adakitic signatures and fertile (Cu-Au-bearing) adakitic andesites, monzonites and monzo-diorites are shown.
To test the role of amphibole versus garnet fractionation (and/or garnet as residue) in the formation of the Quchan Oligocene adakites we used a series of plots, for example, La/Yb versus SiO 2 (Figure 10c), La/ Nb versus Nb/U (Figure 10d), Dy (n) /Yb (n) versus La (n) /Sm (n) (Figure 10g) and Dy/Dy* versus Dy/Yb (Figure 10h). Amphiboles should have high Nb/U ratios and therefore a decreasing Nb/U along with increasing La/Nb can be considered as an indication of amphibole fractionation. Also, trends observed in Dy/Dy* versus Dy/Yb and Dy (n) /Yb (n) versus La (n) /Sm (n) can be considered as reflecting amphibole fractionation (Davidson et al., 2013). Garnet fractionation would increase La/Yb and Dy (n) /Yb (n) in plots of Dy (n) /Yb (n) versus La (n) /Sm (n) and La/Yb versus SiO 2 , which is not the case for Quchan adakites but instead can be seen in fertile adakites from S-SW Neyshabour. There is also a fractionation trend extending from mafic rocks toward dacites-rhyolites in plots of Zr/Hf versus Sm/Hf (Figure 10f), which can show fractionation of clinopyroxene and amphibole. Since amphibole and pyroxene have higher Kd for Hf than for Zr, contributions of trace minerals such as zircon, allanite and xenotime with Kd(Hf)/Kd(Zr) < 1 (Bea et al., 2006) are necessary to explain the fractionation of the Zr/Hf ratio.
There are few geochemical-isotopic data for the Quchan adakites, but available data suggest they were produced by partial melting of the lower continental crust of NE Iran, which has been thickened enough (40-50 km) to generate such melts and/or facilitate high-pressure differentiation of basaltic melts in the deep crust (e.g., Shabanian et al., 2012). However, high-pressure fractionation trends involving garnet cannot be reconciled with our geochemical data, especially the N-MORB-normalized spoon-shaped MREE-HREE patterns for dacites and rhyolites (Figure 7). Instead, we propose that the elemental ratios discussed above suggest that the adakitic signatures of the felsic adakitic rocks (except non-adakitic samples AD12-1, AG12-1, AG12-4, and CG12-1) could have originated from crustal assimilation and extreme fractionation of amphibole ± clinopyroxene ± plagioclase and accessory minerals from a mafic melt similar to the Quchan basaltic melts and/or modeled melt obtained from Section 5.1. In addition, the elemental patterns such as La/Yb versus SiO 2 and/or Dy (n) /Yb (n) versus La (n) /Sm (n) could reflect variable fractionating assemblages, changes in pressure of fractionation and/or the input of new pulses of magma into a pre-existing chamber.

FC-AFC-REFC Processes and Silicic Magmatism
Our elemental ratios and isotopic data as well as results obtained for amphibole chemistry indicate that crystal fractionation (amphibole dominant) and assimilation of continental crust (AFC) by the Quchan basaltic melts could play an important role in the genesis of silicic magmas (both adakites and non-adakitic rocks) in Quchan. Field, petrography, and isotopic data show evidence for assimilation of crustal host rocks including Cadomian biotite-rich paragneisses. The evidence includes xenoliths in felsic rocks, rounded quartz grains in dacites and rhyolites and variable bulk-rock Sr-isotope compositions as well as zircon εHf(t) values. Major-and trace-element data could also attest to assimilation, which we discuss below. We suggest that the Oligocene Quchan collisional volcanic rocks might represent melts that were extracted from a metasomatized mantle reservoir (see Section 5.1) and evolved via fractional crystallization and assimilation-fractional crystallization in a magmatic plumbing system developed in the Cadomian crust of NE Iran. In order to assess this hypothesis, we tested models invoking Rayleigh FC (e.g., Moghadam et al., 2016;Wanless et al., 2010;, coupled with AFC (e.g., DePaolo, 1981; and REFC in mafic magmatic reservoirs . The complete workflow of the FC-AFC-REFC models including, (a) the possible magmatic parental melts and assimilated crustal rocks, (b) calculation of trace-element partition coefficients, and (c) presentation of the FC-AFC-REFC equations is explained in Supporting Information S1. All the parameters and results obtained from FC-AFC-REFC modeling are presented in Tables 4 and 5 and Tables S4 and S5.

Fractional Crystallization (FC)
Major elements mass-balance modeling has been used to simulate the role of FC and AFC processes in the formation of Quchan intermediate and felsic (both adakites and non-adakitic-) rocks from basaltic pre-cursors. The quality of the calculated models is evaluated through the minimization of the sum of the squared residual (Σr 2 ) of the Parental Melt versus Daughter Liquid + Fractionating Mineral Assemblage relationship. Major element mass-balance models are considered acceptable when Σr 2 <1.0. The calculated models show that basaltic samples such as CG12-4 (SiO 2 49.1 wt%) and DA12-1 (SiO 2 53.4 wt%) reflect ca 30%-45% FC of an average basaltic parental melt (MB-av; average composition of samples MB12-1 and MB12-2 with MgO values of 7.1 and 6.8 wt%, respectively), with crystallizing assemblages of hornblende + plagioclase together with Cr-spinel and olivine (Hbl 14-59 + Pl 18-67 ± Spl 0-28 ± Ol 0-17 ; hereafter percentages are given for the total crystallized assemblage; see Tables 4 and 5). However, the resulted crystallizing assemblage fails (Σr 2 ∼3.6-16.2), mainly due to the large r 2 for Na, to directly reproduce the differentiation trend from the basaltic source to the Quchan andesites (see representative models 4, 9, and 10 and replicated models 31-33 using plagioclases with different Na 2 O contents, in Table S4), and to some dacites and rhyolites (see representative models 34 to 41, in Table S4).

Assimilation-Fractional Crystallization (AFC)
The failure of fractionation of basalts to yield andesites, dacites, and rhyolites strongly suggest the possible contribution of crustal materials through AFC; mass-addition, mainly during the genesis of fractionated melts (e.g., McBirney et al., 1987;. To test this hypothesis, we applied AFC models based on the interaction between Quchan collisional mafic melts (i.e., basalts and andesites) and the Cadomian felsic crust that outcrops in NE Iran and hosted the Quchan plumbing system. Maximum mass-addition (assimilation) trends calculated through the lever-rule method are reported in Table 5 and Figures 11 and 12. In agreement with the method proposed by , if a pure crystal-fractionation (FC) trend corresponds to the null-assimilation value, then the compositional space between the FC-and the maximum mass-addition (assimilation) trends correspond to the possible compositions of residual liquids produced by AFC-processes.  .

Table 4 Bulk Compositions Used in the FC-AFC Models
Results obtained from AFC modeling suggest that the Quchan andesites could represent basaltic melts that have undergone ca 50% coupled assimilation and fractional crystallization processes. On the other hand, the direct effect of assimilation of the Cadomian felsic rocks appears to be limited in the genesis of dacites and rhyolites from the andesitic melts.
Based on these results, we further investigated the evolution of Quchan melts through FC-AFC modeling applied to representative trace-(Sr and Y) and rare-earth (Dy and Yb) elements (results are in Table S5 for FC and AFC models). Solutions for these FC-AFC models are consistent with those for major-element models, supporting the scenario of Quchan magmatism controlled by FC coupled with assimilation processes. Rayleigh FC models calculated for the Bas1 assemblage generally fail to produce andesites and felsic melts (Figure 11), whereas a combination of fractional crystallization (FC ∼ 40-60%) and assimilation up to ca 30% (proportion of assimilant to fractionates r = 0.3) of Cadomian felsic crust seems capable of generating andesites from Quchan basalts. Dacites and rhyolites (including adakites) are then the residual liquids after 50%-70% FC of an andesitic parental melt ( Figure 12). Higher rates of assimilation (>30%) of Cadomian felsic crustal rocks are still recognizable. However, at this stage, we cannot rule out the inheritance of Note. C 0 stands for the initial concentration of an element in the parental liquid, C A is the concentration of element in assimilant, C L is the concentration of element in the generated liquid. a Proportion of assimilated rock.

Table 5
Major-Element AFC Model for MB-av and SB-av to the Cadomian Crust this signature from the previous AFC processes which were responsible for the genesis of andesites, which can suggest there is no prerequisite for high-degree assimilation of Cadomian rocks.

Recharge, Evacuation, and Fractional Crystallization (REFC)
The wide compositional range of the Oligocene lavas in NE Iran, and the fact that some samples do not follow the AFC trends (Figure 12), together with the complex FC-AFC processes involving different fractionating assemblages and various degrees of crustal assimilation, suggests a magmatic build-up in deep "hot zones" and the persistence of an extensive magmatic feeding system maintained by a continuous recharge and evacuation of mafic magmas at the time of collisional magmatism. This idea can be also supported by the complex zoning of plagioclase and amphibole phenocrysts. A closed magmatic system without a magmatic recharge, in which the evolution of melts is driven only by FC-processes, will produce progressively smaller volumes of more residual liquids (e.g., Albarède, 1996; with up to 50% of melt cooling and crystallization, <2500 years as a result of the interaction with cold host rocks . In contrast, a continuous magmatic recharge in deep crust reservoirs will decrease both cooling and crystallization rates, and thus will sustain the mass of the residual liquids and can trigger a repeated Figure 11. Fractional crystallization (FC) and assimilation-fractional crystallization (AFC) major element modeling of Quchan volcanic rocks. Major elements have been recalculated to 100% anhydrous, in the system SiO 2 -TiO 2 -Al 2 O 3 -FeO*-MnO-MgO-CaO-Na 2 O-K 2 O. Orange and yellow stars represent the presumed parental melts; MB-av basalt (Bas1 and Bas2 parental melts) and SB-av andesite, respectively. Blue and black lines represent Bas1 and Bas2 fractionation trends respectively, as calculated for MB-av basaltic melt. The Red line represents And fractionation trend as calculated for SB-av andesitic melt. Green and pink lines represent maximum assimilation trends as calculated for the suspected parental SB-av andesitic and MB-av basaltic melts, respectively. The domains of possible AFC processes are presented in green and purple shaded areas (for SB-av andesitic and MB-av basaltic melts, respectively) situated between maximum assimilation/mixing-and FC-trends. The compositions of the Cadomian assimilant and the NE Iran Cadomian crust are taken from . The detailed description of the FC-AFC modeling using major elements is presented in Supporting Information S1. ascent of melts into the upper parts of the feeding system (Annen et al., 2006a;. Furthermore, it has been demonstrated through numerical modeling that deep-crustal chambers undergoing prolonged or repeated recharge events are less prone to magma-evacuation concerning the shallower stagnation layers . Such recharge can produce growing mafic to intermediate magma reservoirs by deforming and assimilating the surrounding lower-to middle crustal rocks (e.g., . To evaluate the effect of the growing deep reservoir(s) on Quchan magmatism, we examined the REFC model of . However, since the literature contains few well-documented examples of REFCtype magmatic systems (e.g., Portnyagin et al., 2015 and references therein), we approached the mode- The percentages indicate the amounts of fractionating assemblages. Equations and parameters used are presented in Supporting Information S1. We also tested the recharge, evacuation and fractional crystallization (REFC) model of , using the same partition coefficients as for the FC-AFC trends. The REFC modeling was applied to both basaltic (MB-av) and andesitic (SB-av) parental melts and has been calculated to 10 times (10×) of the initial mass of the magma chamber (M 0 ch = 1). See text for explanation and Supporting Information S1 for details of the calculations.
ling through a conservative hypothesis assuming (a) the recharge rate (dM r = 0.4) is higher than the sum of crystallization and evacuation rates (dM x + dM e = −0.3), and this corresponds to a magma chamber growth of 10% (DM ch = +0.1) at each overturn; and (b) possible large fractionation between incompatible trace and REE elements can be achieved as a function of the crystallizing phases involved (e.g., Lee et al., 2006;Portnyagin et al., 2015). Using the same partition coefficients as in our FC-AFC models, the REFC modeling was applied to both basaltic (MB-av) and andesitic (SB-av) reservoirs and calculated up to 10 times (10×) the initial mass of the magma chamber (M 0 ch = 1). Results are reported in Table S5 and presented in Figure 12. Compared to pure FC, the REFC modeling for the basaltic reservoir suggests an interesting scenario characterized by (a) a progressive depletion of Sr in residual liquids compared to the Sr enrichment in the FC-trend, (b) a stronger depletion of Y and (c) less fractionation of HREEs. We emphasize that the progressive growth of the basaltic reservoir seems to have no direct influence on the HREE budget in the evolved melts. On the other hand, nearly all Quchan volcanic rocks are included in the compositional space between FC and REFC trends, when departing from the same MB-av parental melt in the Sr-Y plot (Figure 12a).
The complexity behind the evolution of magmas in the Quchan volcanic system can be described by three samples from the same AG locality: AG-12-3 andesite, AG12-4 dacite and AG12-1 rhyolite. These three samples (a) show Sr-Y signatures close to the Cadomian assimilants ( Figure 12a); (b) can be interpreted as daughter liquids after ca 50% fractionation of the andesitic SB-av melt; (c) can be interpreted as fractionated melts extracted by a progressive growth (1.3×) of the deeper basaltic reservoir; (d) can be produced by 60%-80% AFC of a pristine basaltic melt; or (e) can represent the evolved melts extracted from a growing (1.1-1.2×) andesitic reservoir (Figure 12a). We suggest that all these scenarios may have occurred during the genesis of AG rocks, as might be expected in a real magmatic plumbing system. However, since non-adakitic sample AG12-1 has higher εNd(t) and εHf(t), we suggest this sample can be considered as a fractionated melt, extracted during the progressive growth of deeper basaltic reservoirs with more juvenile isotopic signatures than the other rocks.
The FC-AFC-REFC models presented in this study highlight the general perspective of magmatic plumbing systems toward an innovative vision of magmatic feeding systems made up of multiple types of transport of basaltic magmas and storage chambers distributed within the crust. In these settings, multi-stage open-system processes such as the FC-AFC processes can lead to the evolution of magmas to yield felsic rocks and/or via AFC and extreme amphibole fractionation to produce adakites.

Isotope Modeling of Subducted Trench Sediments Added to DM
We have used Sr-Nd-Pb isotope modeling ( Figure 13) to further investigate the role of mantle heterogeneity as well as crustal contamination (or AFC) in the genesis of the Quchan felsic rocks. However, the other ambiguity is that the Quchan basalts show depletion in Nb-Ta except for sample MB12-2 which lacks such depletion; all these basalts have Nb/La∼0.6 to 1. Sample MB12-2 with less of a subduction signature has REE-and trace-element patterns similar to enriched MORBs and shows derivation from an enriched mantle. However, this sample with εNd(t) = 4.4 and εHf(t) = 9.5, still has radiogenic Pb-isotope ratios with ∆8/4-60.1, which show the involvement of components from subducted sediments in their mantle source. Using trace elements, we suggested that a 90:6:4 mixture of DM, terrigenous trench sediment and altered Figure 13. Mixing mass balance modeling of the initial 143 Nd/ 144 Nd (a) and 87 Sr/ 86 Sr (b) isotopic ratios versus 206 Pb/ 204 Pb, using the depleted mantle (DM) and Cadomian continental crust endmembers to unravel the isotopic signatures of the Quchan magmatic rocks. Mixing mass balance modeling is also presented between mixed mantle (1) (MM (1), i.e., 90% depleted mantle + 6% subducted sediments + 4% altered oceanic crust, AOC) and the Cadomian continental crust of Iran. Diamond marks on curves show 1%-10% mixing between endmembers. Data for Pacific MORBs and Pacific arc tholeiites are from EarthChem (https://www. earthchem.org). Geochemical data for Late Cretaceous to Pleistocene magmatic rocks are from . All data including those for MORBs and arc tholeiites have been corrected for 24 Ma radiogenic growth. oceanic crust, following 5% aggregated fractional melting (and/or 89:11 mixture of DM and Eastern Mediterranean Sea sediments, following 6% aggregated fractional melting) could produce modeled melts that are geochemically close to the composition of the Quchan basalts (except for sample MB14-2 which needs an enriched mantle with less input from the subducting sediments). Again, we clarify that our trace elements model is not a unique solution to explain the genesis of all basaltic rocks, but this model is geologically reasonable. To better understand the role of mantle heterogeneity and crustal contamination, we have further used mass-balance modeling to calculate the Sr, Nd and Pb isotopic composition of a mixed mantle (1) and (2) (MM(1) and MM(2) in Figure 9) sources using both mixtures of 90:6:4 of DM, terrigenous trench sediment, altered oceanic crust, and 89:11 of DM and Eastern Mediterranean Sea sediment-probably similar to the depleted mantle-sediment melange model suggested by Marschall and Schumacher (2012 (Figure 13). Then the mixing trends between a depleted mantle (DM) and Cadomian continental crust and also between the mixed mantle (1) (MM(1)) and the Cadomian continental crust have been modeled. All parameters for modeling are presented in Table S6.
Our modeling displays two features. (a) The Quchan samples show two mixing trends (Figure 13), one between the depleted mantle and Cadomian continental crust and the second one between the mixed mantle (1) and the Cadomian continental crust. These trends are better recognized in the Nd versus Pb isotope plot. (b) The mixed mantle (1) can produce the Nd-Pb isotopic signatures of the basalt samples MB12-1 and MB12-2, although the 87 Sr/ 86 Sr (t) value in sample MB12-2 is a bit higher than in sample MB12-1. This mixed mantle (1) can melt to produce the Quchan basalts, although it should produce basalts with Nb-Ta depletion. Since basalt sample MB12-2 shows no clear Nb-Ta depletion (and together with sample MB12-1 displays less enrichment in Rb, Ba, U and Th compared to the modeled melt), it could be derived from a distinct, but enriched source-compared to the other Quchan basalts-possibly with the addition of lesser amounts of subducting sediments. Therefore, the isotopic data require at least three components, which are also recognized in the trace elements of basalts and felsic rocks using our modeling for the mantle source (Section 5.1) and FC-AFC processes (Section 5.3). Our modeling suggests that the isotopic signatures of the Quchan lavas most likely reside in both the lithospheric mantle and in the Cadomian continental crust of Iran.

Formation of Quchan Rocks
Our results are graphically summarized in Figure 14 and show that the geochemical-isotopic variations in Quchan lavas can be provided first by both sub-continental lithospheric mantle and then by the Cadomian continental crust of Iran. However, since the Quchan rocks have mostly high εNd(t)>+4 and εHf(t)>+6, it can be assumed that these isotopic values are not consistent with the presence of a lithospheric mantle (unless that lithospheric mantle was itself a previous subducted slab), but instead could attest to a variably depleted mantle beneath an oceanic back-arc or an intraoceanic arc system. We must emphasize that there is neither any trace of an oceanic back-arc basin during the Late Oligocene in NE Iran nor evidence for an intraoceanic arc within the narrow Sabzevar-Torbat-e-Heydarieh fossilized back-arc basin. However, the radiogenic Nd and Hf isotopes for the Quchan rocks can be inherited from the pre-existing depleted mantle which prevailed beneath the Late Creta- Figure 14. Schematic representation of the magmatic processes to yield the Quchan mafic and felsic rocks. Melting of a depleted to enriched (refertilized) lithospheric mantle and further fractional crystallization of olivine (Ol) and pyroxene (Py) in lower crust magma chambers produce the more evolved melts. These evolved melts further undergone fractional crystallization and assimilation-fractional crystallization in lower to middle magma chambers (hot zones) along with extreme amphibole crystallization to produce felsic and/or adakitic rocks. Black dikes display the melt conduits for primitive basaltic melts to enter the Melting, Assimilation, Storage, and Homogenization (MASH) zones in the deep lower crust, whereas colored dikes show the migration of evolved magmas into shallow magma chambers and/or into the surface.
ceous-Paleocene Sabzevar-Torbat-e-Heydarieh back-arc basin. According to our major-and trace-element and Sr-Nd-Pb modeling, we suggest a four-stage scenario for the formation of the Quchan lavas. (a) Melting of a metasomatized mantle, a DM re-fertilized by altered oceanic crust, trench sediments and/or sediment melts, to produce the original basaltic melts. The refertilization of the mantle can have occurred during the previous subduction, that is, from Late Cretaceous to Oligocene (ca 27 Ma). However, our isotopic modeling shows that both depleted and refertilized mantle sources seem to be involved during the magma genesis.
(b) Early fractionation of olivine and pyroxene from basaltic melts in lower-to middle-crust hot zones to produce more-evolved (andesitic) melts with higher H 2 O contents. (c) Extensive fractionation of amphibole and assimilation of continental crust (AFC) in the middle crust and shallow magma chambers, to produce the adakitic magmas with spoon-like MREE-HREE patterns. This process would further perturb the isotopic ratios of the Quchan magmas. (d) Recharge of magma chambers to perturb the geochemical-isotopic systematics of fractionated magmas, which can further suffer AFC processes to generate the felsic rocks. Empirical models suggest that hydrous basaltic magmas can stall in the lower-middle continental crust, in hot zones, and fractionate to produce crystal mush (e.g., Annen et al., 2006aAnnen et al., , 2006bBachmann & Bergantz, 2008;Hidalgo & Rooney, 2014;Hidalgo et al., 2007). The crystal mush or the hot zones can be replenished by new pulses of basaltic magmas. These processes that is, replenishment and thus magmatic fractionation, can explain the oscillatory zoning in amphibole and plagioclase and the variable isotopic signatures of the fractionated rocks. Extreme amphibole fractionation could be expected to generate amphibole cumulates in the middle-lower crust. Outcrops of these cumulates are lacking in the Quchan area, which is covered by Oligocene magmatic rocks and younger volcanic-sedimentary rocks, but there are 1-2 m thick layers of amphibole (±<10 clinopyroxenes) cumulates in other outcrops in NE Iran that accompany deep-seated, deformed Eocene granodiorites and diorites. These cumulates are similar to mono-mineralic amphibole cumulates reported from continental arcs worldwide, for example, Antarctica (Tiepolo & Tribuzio, 2008), south China (Sun & Zhou, 2008) and from southern and central American volcanic arcs. Crystallization of amphibole cumulates in the middle crust, along with assimilation of continental crust, could rapidly increase the SiO 2 content of the residual magmas (e.g., Carmichael, 2002), which could produce Quchan dacites and rhyolites upon eruption. Amphibole fractionation and AFC processes further can generate residual magmas with adakitic geochemical compositions (Castillo et al., 1999).

Magmatic Triggers in NE Iran
Zircon U-Pb, K-Ar, and Ar-Ar ages imply that magmatism in the NE Iran back-arc started in the Mid-Late Cretaceous (∼110 Ma) and continued until the Pleistocene (e.g., Alaminia et al., 2013;Ghasemi et al., 2010;Jamshidi et al., 2015;Kazemi et al., 2019;Rostami-Hossouri et al., 2020;Shabanian et al., 2012, among others). Late Cretaceous magmatism in the NE Iran back-arc is consistent with subduction initiation along the Zagros suture zone or the Main Zagros Thrust (MZT) at ca 110-100 Ma and along the Makran zone in south Iran (e.g., Barbero et al., 2020;Burg, 2018;Esmaeili et al., 2020;Moghadam & Stern, 2011;Moghadam et al., 2010;Moghadam, Khedr, et al., 2014;Monsef et al., 2018) (Figure 1a). Subduction initiation drove extension in the overlying plate of the Iranian plateau, caused exhumation of high-pressure rocks (blueschists) along with the Zagros ophiolites (Angiboust et al., 2016; and triggered the formation of back-arc oceanic basins within the Iranian plateau, such as the Sabzevar-Torbat-e-Heydarieh oceanic back-arc in NE Iran (Moghadam, Corfu, et al., 2014). This extension was also responsible for a high rate of magmatism in the Iranian plateau, mainly in the NE segment. This Mid-Late Cretaceous magmatism in the NE Iran back-arc is represented by arc tholeiites and calc-alkaline rocks ( Figure 6).
The NE Iran magmatism entered a waning stage during Paleocene to Early Eocene time coeval with the exhumation of blueschists and thus the closure of the Sabzevar-Torbat-e-Heydarieh back-arc oceanic basin (Bröcker et al., 2021). The Paleocene to Early Eocene magmatism is characterized by highly radiogenic arc-tholeiitic to calc-alkaline igneous rocks. Magmatism began to increase during the Middle Eocene (48-40 Ma), coeval with flare-up magmatism in other segments of both the magmatic front and back-arcs in Iran van der Boon et al., 2021). High magmatic fluxes during the Eocene were accompanied by extension throughout the Iranian plateau driven by the roll-back of the Neotethyan slab (Chiu et al., 2013;Verdel et al., 2007). The extension was accompanied by core-complex formation, basin deepening and subsidence, leaving a thick sequence of Nummulite-bearing limestones and green marine pyroclastic rocks, which are accompanied by and/or interlayered with Eocene volcanic rocks ( Eocene magmatism in NE Iran shows both calc-alkaline and adakitic geochemical signatures ( Figure 6). Magmatism seems to have waned during the Late Eocene and restarted during the Mid-Late Oligocene, after the collision between Arabia and Eurasia at ∼27 Ma. Magmatism continued during Miocene and Pleistocene time in NE Iran with the eruption of subaerial andesites from NE Iran (NW Torbat-e-Heydarieh). Oligocene and younger magmatic rocks from NE Iran show enrichment in large ion lithophile elements and depletion in high field strength elements, and these magmas are inferred to be mainly derived from a subduction-modified mantle source, refertilized by subduction components from subducting Neotethyan oceanic lithosphere beneath Iran, before Iran and Arabia collided at ca 27 Ma.
Since there is ∼600 km between the magmatic front (Urumieh-Dokhtar magmatic belt, Figure 1) and backarc region in NE Iran, some believe that the Cenozoic magmatism in NE Iran could be related to the subduction and thus closure of the Neotethyan branches or the Sabzevar-Torbat-e-Heydarieh oceanic back-arc basin in NE Iran (e.g., Moghadam, Khedr, et al., 2015). We would argue that (a) the distance between the magmatic front and back-arcs is controlled mainly by the dip of the subducting slab; (b) in most retreating arcs, the back-arc magmatic belt is mostly restricted to ∼100-300 km from the magmatic front (as in the southern Andean subduction system); (c) in some subduction systems-for example, extensional ones-the back-arc magmatism is located far away from the trench. For example, the Cenozoic back-arc magmatism in the Oligocene San Juan Volcanic Field (USA) could be traced ∼1,000 km away from the trench (Lake & Farmer, 2015).
We cannot rule out the subduction of the Neotethyan back-arc lithosphere in NE Iran, but the polarity of this subduction has been challenged. Many argue that the back-arc lithosphere has been subducted toward the north-beneath the Turan Plate-and thus has produced the Quchan Oligocene rocks (e.g., Ghasemi et al., 2010), but without a trace of Late Cretaceous-Eocene magmatism. However, the "collisional" or "collision-related" terms we use in this study refer to the magmatism occurring after the collision between Iran and Arabia at ∼27 Ma, which is abundant in all parts of Iran. The tectonic forces needed for continental collision-that is, those related to the opening of the Red Sea and convergence between Arabia and Iran-could also be a trigger for the closure of the Sabzevar-Torbat-e-Heydarieh oceanic basin in NE Iran, although the exhumation of the Paleocene-Early Eocene blueschists in the Sabzevar-Torbat-e-Heydarieh oceanic basin (Bröcker et al., 2021) and the deposition of Paleocene-Eocene basal conglomerates (+terrigenous sediments upwards the section) to seal these ophiolites show that this Neotethyan back-arc basin was closed earlier than Late Oligocene.

Conclusions
Late Cretaceous to Pleistocene magmatic rocks are common in NE Iran and are related to the subduction initiation beneath the Iranian Plateau during Late Cretaceous and then to a magmatic flare-up in Eocene time due to the Neotethyan rollback, and finally to the collision between Iran and Arabia at ca 27 Ma. Late Oligocene volcanic rocks including basalts, andesites, dacites and rhyolites are abundant in NE Iran (south Quchan), with the felsic rocks showing zircon U-Pb ages of 25-24 Ma. These rocks have subduction-related geochemical signatures and are characterized by variable Sr-Nd-Pb-Hf isotope compositions. The geochemical compositions of Quchan basalts are generally in agreement with partial melting of a refertilized mantle, generated by the interaction of a depleted mantle with subducting altered oceanic crust and overlying sediments. However, the Sr-Nd-Pb isotopic ratios and modeling could highlight that both a depleted and an enriched mantle were involved in the genesis of the Quchan rocks. Fractional crystallization and assimilation-fractional crystallization models imply that basaltic melts could pool in the lower crust and fractionation of early olivine and clinopyroxene can yield fractionated, H 2 O-rich melts. These evolved melts can foster extreme amphibole fractionation in the lower-middle crustal hot zone, along with assimilation of the Cadomian continental crust, to produce the Quchan felsic rocks and/or adakites.

Data Availability Statement
All data underlying the finding of this paper can be accessed from both Supporting Information S1 and http://dx.doi.org/10.17632/8zm7zkrrnb.1.