Melting of Subducted Continental Crust in the Cimmerian Collisional Orogenic Belt, NE Iran Open Access

Continental crust is generally thought to be resistant to subduction during or following continental collision because its more felsic composition results in lower density relative to mafic oceanic crust. However, the identification of (ultra)high-pressure minerals such as phengite, coesite, and diamond in eclogite from high-pressure metamorphic massifs provides evidence for the recycling of continental crust material into the deep mantle [1, 2]. The similar density of the upper continental crust and ambient mantle has been revealed by experimental studies performed at mantle pressure conditions [3, 4], suggesting that nonmafic crustal material can also descend into the Earth’s upper mantle.

Subduction of continental crust can occur after the oceanic crust between two continental terranes has been consumed, as long as the drag from the subducting basaltic plate continues [1, 5]. Melts derived from the underthrust continental crust could modify the orogeny lithospheric mantle, as indicated by postcollisional mafic rocks [6, 7]. Nevertheless, continental crust subduction is generally thought to have taken place where (ultra)high-pressure metamorphism has been documented, such as the Dabie-Sulu and Qilian orogenic belts in China, the Kokchetav complex in Kazakhstan, and the Bohemian Massif, Germany [1]. However, it is unclear whether the subduction and partial melting of continental crust operated in ancient collisional belts where no collision-related ultrahigh-pressure metamorphism has been documented [8]. This limits the full evaluation of the recycling of continental crust and its possible contributions to mantle heterogeneity.

It is challenging to recognize continental crust subduction in orogenic belts where no diagnostic ultrahigh-pressure mineral phase occurs in metamorphic rocks of continental origin. Fortunately, because of its common occurrence and chemically resistant and refractory properties, zircon from magmatic rocks provides potential tools to explore the recycling of continental materials [6, 9]. For example, recycling of continental crust via delamination of thickened continental lithosphere was suggested by zircons from garnet/spinel pyroxenite veins within mantle xenoliths from the Hannuoba basalts, North China [10]. Inherited zircon cores with magmatic overgrowths found in granitic rocks have generally been interpreted to be the result of crustal contamination and can also be used to trace their source provenance [9, 11]. For example, Eocene to Miocene andesitic to rhyolitic arc rocks in the Southern Mountains Arc, Java contain a range of Cambrian and older inherited zircons, which have been interpreted to be the result of partial melting of a Gondwana continental fragment beneath the arc [11]. More generally, a combination of numerical modeling of crustal rock melting and zircon solubility has been used to argue that inherited zircons in granitic magmas were mostly derived from the source region [12].

Late Triassic intermediate to felsic magmatic rocks exposed in NE Iran have been interpreted to be the magmatic response to the closure of the Paleo-Tethyan Ocean and the subsequent collisional event. The genesis of these igneous rocks has been studied previously, but their origin is still debated [13-16]. More importantly, the genetic relationship between them and the central Iran continent has not been well explained, and they may contain valuable information for deciphering the subduction and partial melting of the continent crust.

In the present study, we report zircon U-Pb and Lu-Hf isotopes as well as bulk-rock geochemical data for twenty-two samples from three Late Triassic granitic plutons in Mashhad, NE Iran. Magmatic rocks in the northwest and southeast portions of Mashhad contain zircons with relict cores rimmed by magmatic overgrowth zones. These data combined with whole-rock geochemical data were used to examine their source and shed new insights on the subduction and melting of Central Iran’s continental crust during the Cimmerian orogen.

The Cimmerian orogen in northeast Iran formed from the collision between the Turan terrane in the north and the Central Iran terrane in the south (Figure 1(a)). The Paleo-Tethyan suture zone (Figure 1(b)) between them is characterized by mafic-ultramafic-chert complexes in the Alborz and Kopet Dag mountains [17-19].

South of the Paleo-Tethyan suture, the Central Iran terrane is characterized by extensive and widespread Neoproterozoic to Early Cambrian magmatism with ages ranging from ca. 640 to 520 Ma [20, 21]. The calc-alkaline signature of these Neoproterozoic magmas has been linked to south-dipping subduction beneath north Gondwana [21]. In contrast, Late Cambrian to Early Devonian basalt to andesites, alkali gabbro, and A-type magmas [17] were related to the opening of the Paleo-Tethyan Ocean.

During the Late Paleozoic to early Mesozoic, consumption of the intervening Paleo-Tethyan oceanic slab through northward subduction beneath the Turan terrane (southern margin of the Eurasia continent) resulted in the development of an active Andean-type continental margin and an extensive magmatic arc [17, 19, 22]. This is supported by emplacement of Middle-Upper Devonian boninitic rocks, calc-alkalic basaltic to intermediate volcanic rocks, and ca. 350 Ma eclogite related to oceanic slab subduction [19, 23-25]. At this time, the Central Iran terrane was a passive margin dominated by platform-type deposits [20]. The closure of the Paleo-Tethyan ocean and the final collision of the Central Iran and Turan terranes occurred in the Late Triassic-Early Jurassic [19, 26]. In northeastern Iran, the Cimmerian orogeny is represented by syn- to postcollisional clastic successions of the Shemshak Group conglomerate and sandstone [20, 27].

Metamorphism related to the Cimmerian orogeny in NE Iran varies from lower greenschist to amphibolite facies [28, 29]. The Permian Mashhad metamorphic complex exposed in the north of the Binalud mountains consists of schist, phyllite, marble, and amphibolite, which has been interpreted as an accretionary wedge [26]. The Triassic to Early Jurassic Mashhad phyllite, occupying most of the Binalud mountains, is a greenschist facies sequence including phyllite, slate, meta-graywacke, quartzite, and sandstone [28].

In the study area, Late Triassic magmatic rocks from northwest to southeast are composed of diorite to granodiorite in the Vakilabad and Kuhsangi area (NW Mashhad pluton in following sections), granodiorite and two-mica granite in southeast Mashhad (Figure 1(c), refer to SE Mashhad pluton), and the Torbat-e-Jam pluton north of the city of Torbat-e-Jam (Figure 1(d)). Intermediate to felsic rocks in NW Mashhad intruded into the Mashhad ophiolite and flysch deposits, both of which are characterized by the presence of ilmenite and low magnetic susceptibility [15]. The Torbat-e-Jam pluton intruded into metamorphosed sedimentary rocks of the Upper Triassic Miankuhi Formation [14]. Generations of these magmas were related to the collision of the Iran and the Turan terranes or to the subduction of the Paleo-Tethyan oceanic lithosphere during the Late Triassic [13, 14, 16, 19].

The Torbat-e-Jam pluton intruded into the Triassic Miankuhi Formation generating contact with metamorphic aureoles a few meters wide [14]. The pluton consists of a range of rock types including gabbro, diorite, granodiorite, and granite, but more evolved syenogranite and aplite are also found. Elongate metasedimentary xenoliths are present near the contact with the wall rocks [14]. The largest portion of the Torbat-e-Jam pluton is occupied by the granodiorite. The medium- to coarse-grained granodiorite (online supplementary Figure S1) is composed mainly of plagioclase, K-feldspar, and quartz, with biotite and amphibole. Brick-red syenogranite shows gradual to sharp contact with the granodiorite. The syenogranite is composed of K-feldspar (>30%), quartz (25%–35%), plagioclase (20%–25%), biotite (5%–10%), and rare amphibole (<3%). Apatite, zircon, and titanite are common accessory phases. Magmatic microgranular enclaves (MMEs) in the pluton are fine-grained and composed of plagioclase, biotite, amphibole, quartz, and K-feldspar, with rare pyroxene (online supplementary Figure S1).

The SE Mashhad pluton is elongated with an NW-SE strike and consists mainly of granodiorite and two-mica granite (online supplementary Figure S1). MMEs (online supplementary Figure S1) are present in the granodiorites but are rare in the two-mica granite [13, 16]. Granodiorite comprises the southeast part of the pluton and is medium to coarse grained (online supplementary Figure S1), with a massive structure with no deformation. The rock consists mainly of plagioclase (25%–35%), K-feldspar (20%–30%), quartz (15%–25%), amphibole (<5%), and biotite (5%–10%) with accessory phases including zircon, apatite, and titanite. The two-mica granite (online supplementary Figure S1) crops out in the northwest part of this pluton and is composed of quartz (25%–30%), K-feldspar (20%–30%), plagioclase (15%–25%), biotite (<5%), muscovite (5%–10%), and minor garnet (<3%). Accessory minerals include apatite and zircon. Field observation shows that the gray to reddish two-mica granite intruded into the dark gray granodiorite.

In the Vakilabad and Kuhsangi area, northwestern Mashhad, Late Triassic magmatic rocks occur as small plutons. The granodiorite is grayish in color with a medium- to fine-grained texture (online supplementary Figure S1). It contains plagioclase (30%–40%) and K-feldspar (10%–20%), with varying amounts of quartz (15%–30%), but the only mafic mineral is biotite (5%–10%). Trace ilmenite occurs locally, with accessory minerals including apatite and zircon.

Nine samples were collected from the Torbat-e-Jam pluton, including two microgranular enclaves, six granodiorites, and one syenogranite. Two diorites, four granodiorites, and two two-mica granites were sampled from the SE Mashhad pluton, with one of the granodiorite samples (MH23) showing gneissic fabrics defined by the oriented alignment of biotite and feldspar. Five granodiorites were collected from granitic plutons exposed in NW of Mashhad, two are from outcrops near Vakilabad village, and the remaining three are from the area near the Kuhsangi Park (Figure 1c).

Zircons were separated with conventional liquid and magnetic techniques and hand-picked under a binocular microscope. Zircon grains were mounted in epoxy and photographed under a JSM-7800F thermal field emission scanning electron microscope equipped with a Gatan MonoCL4 cathodoluminescence spectroscope. All chemical analyses were carried out at the State Key Laboratory of Ore Deposit Geochemistry, Institute of Geochemistry, Chinese Academy of Sciences, Guiyang, China.

Zircon U-Pb isotopic compositions were analyzed with a Geolas Pro 193 nm laser ablation system coupled with Agilent 7500× inductively coupled mass spectrometry (ICP-MS). During analysis, the laser spot diameter was 32 or 24 µm with a pulse rate of 5 Hz. Standard zircon 91500 was used to monitor elemental fractionation and instrumental bias correction. These data are presented in online supplementary Table S1. The secondary standard zircons Qinghu and Plésovice yielded a mean age of 156.9 ± 1.1 Ma (2σ, n = 32) and 338.1 ± 1.4 Ma, respectively, which are consistent with previously reported values [30, 31].

Whole-rock major oxide compositions were determined using an ARL Perform X4200 X-ray fluorescence spectrometer on fused discs. Trace element concentrations were measured by PlasmaQuant MS Elite ICP-MS, and powders were dissolved in an HF-HNO₃ mixture. Whole-rock geochemical results are listed in online supplementary Table S2.

Whole-rock Sr and Nd isotopic compositions were determined on aliquots of the same powders analyzed for major and trace elements. Cation ion-exchange resin AG50W-X12 was used for Sr and rare earth element (REE) separation and purification, and Nd was separated from other REE fractions through P507 extraction chromatography resin. A Triton thermal ionization mass spectrometer was used for Sr and Nd isotopic measurements. For mass fractionation correction, ⁸⁷Sr/⁸⁶Sr and ¹⁴³Nd/¹⁴⁴Nd ratios were normalized to 0.1194 and 0.7219, respectively. The measured values of ⁸⁷Sr/⁸⁶Sr and ¹⁴³Nd/¹⁴⁴Nd for USGS standard reference material BCR-2 are 0.705038 ± 0.000004 (1δ) and 0.512653 ± 0.000003 (1δ), respectively.

Zircon Hf isotopic compositions were measured using a Neptune Plus III multicollector combined with a RESOlution S-155 193 nm excimer laser-ablation system attached to an Agilent 7700× ICP-MS. During the analytical session, standard zircon Penglai (¹⁷⁶Hf/¹⁷⁷Hf = 0.282906 ± 0.000010) was used for fractionation correction. Results of the Lu-Hf isotopes and the correlation between Hf and U-Pb analyses are shown in online supplementary Tables S3 and S4. The secondary standard zircons GJ-1 and Plešovice yielded mean ¹⁷⁶Hf/¹⁷⁷Hf values of 0.282013 ± 0.000009 (2SD; n = 7) and 0.282474 ± 0.000014 (2SD; n = 11), respectively, consistent with previously reported values [31, 32].

Zircon crystals are generally euhedral, transparent, and colorless. They are characterized by oscillatory zoning (Figure 2), and their high Th/U ratios (mostly >0.1) indicate a magmatic origin. Granitic rocks in the NW and SE Mashhad areas contain variable volumes of inherited zircon cores as seen by their sharp contacts in cathodoluminescence (CL) images (Figure 2). The zircons with inherited cores are commonly euhedral and prismatic in morphology, similar to those without inherited cores. The inherited cores are generally identified by the contrasting luminous character compared with the magmatic growth zones and their curved and sharp contacts. The inherited cores vary in sizes from ~10 to 100 µm, with some displaying oscillatory zoning, suggesting a magmatic origin. A few of them lack visible zoning or weakly planar/oscillatory zoning, possibly representative of relicts after partial melting of source rock [9].

Various rocks including microgranular enclaves, granodiorite, and syenogranite samples from the Torbat-e-Jam pluton were dated. Zircons from the microgranular enclave (MH13E) show variable U concentrations ranging from 158 to 3677 ppm. All analyses provided a lower intercept age of 218.5 ± 3.5 Ma (n = 30; Figure 3(a)), and ten spots with U <1200 ppm yielded a Concordia age of 217.8 ± 1.2 Ma (Mean Square of Weighted Deviates (MSWD) = 0.2) and a weighted mean age of 217.8 ± 2.4 Ma (MSWD = 0.2). Zircons from granodiorite sample MH12 have a weighted mean age of 219.9 ± 1.6 Ma (n = 30; MSWD = 0.7; Figure 3(b)) and an intercept age of 218.8 ± 2.1 Ma. The two granodiorite samples MH14 and MH15 yielded mean ages of 217.3 ± 1.2 Ma (n = 45; MSWD = 0.7; Figure 3(c)) and 218.9 ± 1.3 Ma (n = 45; MSWD = 0.6; Figure 3(d)), respectively. A large population of zircon grains from the brick-red syenogranite (MH16) is fractured with low CL response. Their apparent ages vary from 164.4 to 323.5 Ma with four analyses showing concordance less than 90%. After excluding those spots with concordance <90% and four spots with ages younger than the main population, the remaining analyses yielded a mean age of 217.1 ± 1.3 Ma (MSWD = 0.6, n = 37; Figure 3(e)). The analyses with low U contents (<1000 ppm) yielded a Concordia age of 217.8 ± 0.8 Ma (MSWD = 0.4, n = 27; Figure 3(f)).

A total of forty-seven analyses were completed on synmagmatic zircons from the granodiorite (MH10), with forty-two spots yielding a weighted mean ²⁰⁶Pb/²³⁸U age of 216.0 ± 1.3 Ma (MSWD = 1.2; Figures 3(g) and 3(h)). For the granodiorite from the Kuhsangi park area, two samples (MH35 and MH36) were dated at 216.4 ± 1.7 Ma (MSWD = 0.7, n = 38; Figures 3(i) and 3(j)) and 217.3 ± 1.3 Ma (MSWD = 1.0; n = 41; Figures 3(k) and 3(l)).

Magmatic zircons from the granodiorite sample MH28 in the SE Mashhad area yielded a weighted mean ²⁰⁶Pb/²³⁸U age of 199.9 ± 1.2 Ma (MSWD = 0.9; n = 29; Figures 3(m) and 3(n)). Another sample (MH34) gave a mean age of 200.0 ± 1.3 Ma (MSWD = 0.6; n = 25), consistent with the concordia age of 200.1 ± 0.7 Ma (Figures 3(o) and 3(p)).

The granitic rocks cropping out in the NW and SE Mashhad areas are characterized by abundant inherited zircon cores based on CL images (Figure 2). ²⁰⁷Pb/²⁰⁶Pb ages for zircon with ages >1500 Ma and ²⁰⁶Pb/²³⁸U ages for younger zircon were used in plots and discussion [33]. Fifteen analyses of inherited cores from sample MH10 from NW Mashhad yielded apparent ages ranging from 420 ± 11 to 1939 ± 83 Ma. Three grains show concordance lower than 80% which could be the result of variable Pb loss. Fourteen inherited cores from granodiorite sample MH35 in the NW Mashhad area have ages varying from 615 ± 17 Ma to 2794 ± 59 Ma, whereas ten analyses for inherited cores of sample MH36 yielded apparent ages ranging from 416 ± 6 to 1913 ± 86 Ma.

Granodiorites (MH28 and MH34) from the SE Mashhad pluton also contain inherited zircon cores with Devonian to Archean ages (Figures 4(a)–4(d)). Twenty-two analyses for inherited zircon cores from sample MH28 provided apparent ages ranging from 414 ± 38 to 2426 ± 45 Ma. Twenty-six analyses were measured on inherited zircon cores from sample MH34, and twenty-four spots with concordance >80% yielded apparent ages from 401.7 ± 7 to 2599 ± 47 Ma.

The granodiorite from the Torbat-e-Jam pluton has high SiO₂ contents ranging from 64 wt.% to 71 wt.%, with K₂O varying from 2.94 to 4.63 wt.% (Figure 5(a)). They are metaluminous to peraluminous with A/CNK values (molar ratio Al₂O₃/ CaO+Na₂O+K₂O) ranging from 0.98 to 1.1 (Figure 5(b)). Their Al₂O₃, CaO, MgO, Fe₂O₃t, and P₂O₅ contents decrease with increasing SiO₂, whereas Mg# [molar ratio Mg/(Fe+Mg)] decreases steadily with increasing SiO₂ till 72 wt.% (Figures 6(a)–6(f)). The syenogranites are more evolved with higher SiO₂ (75, 77 wt.%) but lower MgO (0.94, 0.06 wt.%) and Fe₂O₃t contents than those of the granodiorite. The two microgranular enclave samples have intermediate SiO₂ of 57 and 59 wt.% and MgO contents of 3.42 and 2.45 wt.%, respectively.

Granodiorites in SE Mashhad have higher SiO₂ (66, 70 wt.%) than those in the NW with high K₂O (2.4, 4.2 wt.%) plotting within the high-K₂O field (Figure 5(a)). They have high Al₂O₃ (>15.0 wt.%) but low MgO (<3.0 wt.%), Cr (2.41, 100 ppm), and Ni (2.18, 18.3 ppm). Granodiorites in SE Mashhad have higher Sr (474, 1234 ppm) with the exception of MH23 (Sr: 15.9 ppm and La: 68.4 ppm) and La (52.8, 94.6 ppm) than the felsic rocks in the NW Mashhad area (Sr: 335 to 564 ppm and La: 20.8 to 55.5 ppm). Consequently, the SE Mashhad granodiorites show higher Sr/Y and more fractionated (La/Yb)_N ratios (Figure 7). Two diorites hosted in granodiorite from the SE Mashhad pluton have SiO₂ contents of 53 wt.% and 55 wt.%. The two-mica granites from the northern part of the pluton have high SiO₂ (75 and 76 wt.%) contents and are strongly peraluminous with A/CNK values >1.1 (Figure 5).

The various rocks emplaced in Triassic units in the Mashhad area generally show moderately fractionated REE patterns (Figures 8(a)–8(d)) as well as enrichment in large ion lithospheric elements (LILEs) and variable depletion in high-field strength elements (HFSEs, Figures 8(e)–8(f)). Granodiorites from the Torbat-e-Jam pluton have fractionated REE patterns with (La/Yb)_N ranging from 13.0 to 21.9 and moderate negative Eu anomalies (Eu/Eu* = 0.37 to 0.74, Eu/Eu* = Eu/[Sm_N × Gd_N]^1/2, the subscript N refers to chondrite normalization). Granodiorites from the SE Mashhad pluton show enrichment in light rare earth element (LREE) relative to heavy rare earth element (HREE) with (La/Yb)_N varying from 31.8 to 42.9, with weak Eu anomalies (Eu/Eu* = 0.74 to 0.78). The granodiorites from the NW Mashhad pluton have slightly lower LREEs (La: 27.3 to 47.6 ppm) relative to the SE Mashhad pluton (La: 64.2 to 70.8 ppm). They have (La/Yb)_N ranging from 15.0 to 24.2 and Eu/Eu* from 0.78 to 1.06. The trace element patterns of granodiorites from the SE and NW Mashhad plutons are broadly similar to enrichment in LILEs and depletion in HFSEs.

Two granodiorite samples from the Torbat-e-Jam pluton have initial ⁸⁷Sr/⁸⁶Sr ratios of 0.70418 and 0.70596 and εNd(t) values of −0.46 and +0.33, and one syenogranite sample shows a similar Sr isotope composition of 0.70525 and εNd(t) of −0.18. Two microgranular enclaves hosted in the granodiorite have isotopic compositions comparable to Torbat-e-Jam granodiorites with ⁸⁷Sr/⁸⁶Sr of 0.70498 and 0.70669 and εNd(t) of −0.26 and +0.33.

Two-mica granites from the SE Mashhad pluton have radiogenic Sr (initial ⁸⁷Sr/⁸⁶Sr: 0.71045 and 0.71396) and unradiogenic Nd isotopic compositions with εNd(t) values of −6.96 and −6.26 (Figure 9). Granodiorites of the SE Mashhad pluton have initial ⁸⁷Sr/⁸⁶Sr values varying from 0.68050 to 0.70600 and εNd(t) values varying from −2.44 to −2.45, but the lowest initial ⁸⁷Sr/⁸⁶Sr of sample MH23 (0.68050) results from its abnormally high Rb/Sr (13.6) that could cause excessive time correction during calculation of initial ⁸⁷Sr/⁸⁶Sr values and might be geologically meaningless. Two diorites in this pluton are slightly more depleted in their Sr-Nd isotopic compositions (initial ⁸⁷Sr/⁸⁶Sr: 0.70560 and 0.70592, εNd(t): −1.18 and −1.73) than the granodiorites. Granodiorites from the NW Mashhad area have comparable Sr-Nd isotopic compositions (initial ⁸⁷Sr/⁸⁶Sr = 0.70800 to 0.70840; εNd(t) = −4.93 to −6.01) to the two-mica granites but are more enriched than granodiorites from the SE Mashhad pluton.

Zircon grains from five samples including one MME and four host granitic rocks from the Torbat-e-Jam pluton exhibit high similarity in Hf isotope (Figure 10 and online supplementary Figure S2). Zircon grains from the MME sample yielded εHf(t) values ranging from −2.9 to 5.5, with a peak of +3.0. Zircon crystals from the four host granitic rocks have εHf(t) values varying from −4.9 to 5.5, mostly ranging from +1 to +5 with a peak of +2.3.

Syn-magmatic zircons in granodiorite sample MH10 from NW Mashhad yielded ¹⁷⁶Hf/¹⁷⁷Hf values from 0.282330 to 0.282576, with εHf(t) values varying from −10.9 to −2.3 with a peak at −4.5 (online supplementary Figure S2). Four outliers yielded negative εHf(t) values less than −12, which could be mixtures due to the size of the laser beam given the zoning revealed in CL images. Two granodiorites from near Kuhsangi Park in the NW Mashhad show εHf(t) values ranging from −5.5 to +0.3 with a peak of −3.3 and an outlier yielding an εHf(t) value of −11.2. Granodiorites from the SE Mashhad pluton have magmatic zircon ¹⁷⁶Hf/¹⁷⁷Hf values from 0.282382 to 0.282633 and εHf(t) values from −9.5 to −0.7, and most of them were concentrated in a range from −3.3 to −1.8.

Granodiorites and granites dominate the Torbat-e-Jam pluton, and they are characterized by high SiO₂ and low MgO contents. These geochemical features are against the direct partial melting of mantle peridotite and rather suggest the partial melting of basaltic rocks in the crust. The one sample with a higher MgO content of 4.38 wt.% indicates some mantle melt involvement through magma mixing in its origin, as indicated by the occurrence of MMEs. The presence of amphibole and the metaluminous character of most of the intrusions are consistent with I-type magmas [34].

The I-type magmas could be generated by magma mixing or fractional crystallization of mantle-derived magmas with or without crustal assimilation (AFC) assimilation-fractional crystallization, but both these seem unlikely for the granodiorite and granite from the Torbat-e-Jam pluton as explained below. The mafic endmember involved in magma mixing could be represented by gabbro from this pluton on the basis of their close relationship in spatial and temporal. Magmas formed by magma mixing could fall on a linear array defined by mafic and felsic endmembers; however, on the P₂O₅ versus SiO₂ plots (Figure 6(f)), granodiorite and granite show a coherent array, but gabbro displays a more scattered pattern. Moreover, the small variation in zircon εHf(t) of the different granitic rock samples is not consistent with magma mixing, as a large variation (up to 10 units) is expected [35]. Crystallization experiments on basaltic magmas have shown that granitic liquids can be obtained by fractional crystallization along with the formation of large amounts of ultramafic to gabbroic cumulate plutonic rocks [36]. However, the granodiorites and granites are significantly larger than the gabbro at outcrop scale (Figure 1(d)), and this observation conflicts with a fractional crystallization model. Furthermore, the broadly constant whole-rock εNd(t) values over the large range in SiO₂ and Th/La values are also inconsistent with the AFC hypothesis (Figures 11(a) and 11(b)). In addition, the low efficiency of mixing between mafic and felsic magmas appears to preclude the generation of large volume of homogeneous high-K, I-type granitoids [37].

Partial melting of intracrustal mafic to intermediate magmatic rocks has been suggested to explain I-type magmas [38]. Liquids with silicic to intermediate composition and low Mg# have been generated by melting experiment studies on basaltic rocks with granulite residues at 8–12 kbar [39]. Experiments testing partial melting of crustal rocks suggest that high-K, I-type magmas can be derived from mafic to intermediate rocks in the crust with high-K compositions [37]. Granodiorites and granites from the Torbat-e-Jam pluton mostly have K₂O/Na₂O >1 and Mg# <45, which is consistent with the partial melting of basaltic rocks of medium to high-K basaltic compositions. Formation of granodiorite and granite of the Torbat-e-Jam pluton by high-pressure melting (>12 kbar) is unlikely, supported by their high and flattened HREE patterns (Yb >2 ppm, Dy_N/Yb_N varying from 0.9 to 1.1), which is indicative of no garnet residue and thus relatively a low pressure for the magma source region [39]. Given that their slightly depleted zircon εHf(t) values (+0.5 to +5.5) are similar to those of the gabbro from this pluton [14], and their Sr-Nd isotopic characteristics are moderately enriched, partial melting of juvenile crust could explain the formation of the granodiorite and granites of the Torbat-e-Jam pluton.

The generation of the enclaves in the pluton was most likely the result of magma mixing. They have zircon U-Pb ages and εHf(t) values similar to those of host granitic rocks (Figure 10), which is inconsistent with a melting residue origin for these enclaves because relicts should have an age older than magma derived from it. The enclaves are characterized by fine-grained minerals (online supplementary Figure S1), which is not a typical feature of cumulate plutonic rocks. Some enclaves show curved and cuspate margins, as well as quartz ocelli, plagioclase megacrysts with biotite- and amphibole-rich corona, and acicular apatite with large aspect ratios (online supplementary Figure S1). These observations are evident of dynamic interaction of two or more magmas [38]. This is in line with the linear arrays and their intermediate major element composition between the gabbro-diorite and granodiorite-granite assemblages (Figure 6).

Syenogranites and aplite sills contain higher SiO₂ (Figure 6) than granodiorites but show similar Sr-Nd (Figure 9) and zircon Hf isotopic compositions (Figure 10), suggesting fractionation from parental magmas similar to those of the less evolved granodiorites. This is consistent with their gradational contacts with the granodiorite [14], as well as the larger negative Eu anomalies (Figure 8(b)) and lower Sr, Ba, and Ti concentrations than the granodiorites (Figure 8(f)).

Two-mica granites from the SE Mashhad pluton contain elevated SiO₂, Al₂O₃, and CaO, as well as relatively high Ba and Sr, but moderate to high Rb/Sr ratios (1.2 and 13.6). These chemical characteristics are similar to Malashan-Gyirong leucogranites in southern Tibet, which have been suggested to have formed from melts produced by the fluid-fluxed breakdown of metapelite [40]. This is consistent with the radiogenic initial Sr isotopes (0.7074 and 0.7140) and unradiogenic Nd isotope compositions (−5.9 and −7.0; this study and Reference 16) of the studied two-mica granites.

The two-mica granites are unlikely to represent melts formed after protracted fractionation of less evolved magmas similar to the SE Mashhad adakitic granodiorites. Highly fractionated granites are generally characterized by high Rb/Sr, low REE, and prominent negative Eu anomaly, as well as low Nb/Ta ratios [41]. The two-mica granites possess higher REE contents, Eu/Eu*, and Nb/Ta but lower Rb/Sr ratios than ca. 8 Ma highly differentiated garnet granites from Himalaya [42]. The potential parental magmas for the two-mica granites are most likely the coeval granodiorites in the southeastern part of this pluton, as indicated by their similar age and contact relation. However, the different evolution trends of the two-mica granites and granodiorites as illustrated on plots of Ba and Al₂O₃ versus Eu/Eu* (online supplementary Figure S3) imply that the former did not evolve from magmas represented by the granodiorites [16]. In addition, the outcrop area of the two-mica granites is equivalent to or greater than that of the granodiorites (Figure 1(c)); this observation did not support the fractionation hypothesis as fractional crystallization could generate a significantly larger volume of cumulate relative to extracted liquids.

The granodiorites of the SE Mashhad pluton exhibit moderately high Sr/Y and (La/Yb)_n ratios along with low HREE, similar to adakitic rocks [43]. Magmas with adakitic geochemical signatures can form from a variety of sources and processes [8, 43-46], such as (a) partial melting of an oceanic slab, (b) fractionation of basaltic parental magmas, and (c) melting of thickened continental crust. The possibility of these models will be evaluated in the following sections.

Oceanic slab melting in a subduction setting is unlikely. Tectonic and sedimentary analyses indicate that Gondwana-derived Central Iran terrane was amalgamated with the Turan terrane before 220 Ma [26, 27, 47], predating the ca. 200 Ma SE Mashhad adakitic magmatism. Second, most of the granodiorites of the SE Mashhad pluton have K₂O/Na₂O values >1, unlike typical adakites derived from a basaltic slab which are commonly sodium-rich [43]. Furthermore, the granodiorites exhibit Sr-Nd isotope compositions lower than Tethyan midocean ridge basalt of NE Iran ([17]; Figure 9).

High-pressure fractional crystallization of basaltic magma [45] is also not likely as their La/Yb and Sr/Y show no systematic positive correlation with whole-rock SiO₂ contents (Figures 12(a) and 12(b)). High-pressure fractionation in the presence of garnet would deplete evolved melts in heavy REEs as garnet has a high partition coefficient for HREE relative to middle REE, but the expected trend of increasing Dy/Yb with increasing of SiO₂ is not observed (Figure 12(d)). In addition, the genesis of adakitic magmas through fractional crystallization from basaltic parental magmas is commonly associated with a large volume of contemporaneous mafic rocks [46], which is not the fact in the SE Mashhad.

Adakitic rocks in SE Mashhad could not be generated by the melting of the thickened crust. Crustal thickening during orogen development would cause metamorphic eclogitization and densification of the lowermost crust, leading to delamination of the lower crust [48]. Melts derived from delaminated crust that was contaminated by mantle peridotitic rocks should have high Cr, Ni, and MgO [49], but Late Triassic adakitic magmas in SE Mashhad show low Mg#, Cr, and Ni, which suggests a negligible contribution from the mantle. Delamination as a consequence of significantly thickened lithosphere could induce intensive magmatism [49], but this was not the case in NE Iran during the Late Triassic period. Melting at high-pressure conditions in significantly thickened crust where garnet is stable and plagioclase breaks down would generate magmas with extremely high Sr and low Y and HREE. For instance, the Cenozoic adakitic rocks from the Tibetan Plateau have unusually high Sr/Y (even >200 [50]) and La/Yb (up to 149 [44]) that are significantly higher than those of Late Triassic granodiorites in SE Mashhad area. The moderate Sr/Y and La/Yb values of SE Mashhad granodiorites do not support partial melting of the significantly thickened crust.

The microgranular enclaves in this pluton have been suggested to be the result of magma mixing [13], as indicated by the reverse zoning of plagioclase [13] and the presence of needle-like apatite (online supplementary Figure S1). However, this could not explain the formation of the granodiorite of this pluton [13]. The only candidates of crust- and mantle-derived magmas involved in magma mixing are possibly two-mica granites and mafic microgranular enclaves, respectively. On the Sr versus SiO₂ plots (Figure 12(d)), the granodiorites plot away from the straight linear arrays, which is typical of magma mixing, between the coeval two-mica granites and microgranular enclaves. More importantly, the mean Sr/Y ratio (~50) of adakitic rocks is higher than both of the two-mica granites (~20) and enclaves (~40). In general, adakitic rocks produced by magma mixing are characterized by high MgO values [51], but this is not the case here (Figures 6(c) and 6(d)).

Partial melting of the subducted continental crust [8] could be responsible for the formation of the Late Triassic adakitic rocks in the SE Mashhad area. This is supported by (1) their relatively high K₂O/Na₂O values, which indicate that their source region was K-enriched as suggested by experimental works [37, 52]. (2) Their slightly higher Mg# (up to 51) than experimental melts of basaltic rocks at 1–3 GPa [39] implies mild interaction of the continent-derived melts with mantle peridotite [8]. (3) Enrichment in LREEs and depletion in HFSEs are indicators of source rocks similar to continental crust. (4) These rocks have Nb/Ta (7.8, 15.7, and mean of 12.8) and Th/La (0.28, 0.43, and mean of 0.33) values largely overlapping with continental crust (bulk continental crust: 11.4, lower continental crust: 8.3, and upper continental crust: 13.3 [53]). (5) Their Sr-Nd isotope compositions are more enriched relative to granitic rocks of the Torbat-e-Jam pluton that formed by the melting of a juvenile crust. Moreover, the inherited zircon cores from these rocks have Neoproterozoic to Early Paleozoic ages consistent with major magmatism of the Central Iran terrane [20, 21, 54] but different from late Paleozoic–Early Mesozoic magmas of the Turan active continental margin (detailed in section 5.2 [17, 55]). Taken together, partial melting of subducted continental crust is applicable for these magmas production.

Circa 217 Ma granodiorites in the NW Mashhad area have comparable SiO₂ and Mg# to the SE Mashhad adakitic rocks (Figure 6(d)). The NW Mashhad granodiorites have slightly lower LREEs than granodiorites in SE Mashhad but similar HREE patterns. A continent crust origin has been suggested for the SE Mashhad adakitic granodiorites, and the NW samples exhibit whole-rock εNd(t) values even lower than those of the granodiorites in SE Mashhad. It is reasonable to infer that the model of ancient continental crust remelting is also applicable to the NW Mashhad granodiorites. On the Sr/Y versus Y and La/Yb versus Yb binary plots (Figure 7), granodioritic rocks in NW Mashhad straddle the field between the common arc and adakitic rocks. Their Dy/Yb ratios show negligible variation over the range of SiO₂ from 66.0 to 70.2 wt.% (Figure 12(d)), indicating residue of garnet in the source [39, 56]. Because the stability of garnet is pressure sensitive, and the Dy/Yb ratios of NW Mashhad plutons are comparable to the SE Mashhad adakitic rocks, which indicate similar melting pressure conditions for the two intrusions. Granodiorites exposed in both NW and SE Mashhad areas contain variable amounts of zircon with ancient inherited cores. Furthermore, their nearly identical zircon εHf(t) values (online supplementary Figure S2) denote similarity in their source character. On the basis of these observations, granodioritic magmas in the NW Mashhad area could (at least partly) be produced from similar sources of adakitic rocks in the SE Mashhad.

Late Triassic granodiorites in NW and SE Mashhad areas suggest the partial melting of the heterogeneous crust. These rocks are geochemically transitional between I- and S-type magmas. Both lack diagnostic minerals, such as amphibole, muscovite, garnet, and cordierite, which makes their source nature of protoliths enigmatic. An I-type magma affinity is indicated by (a) their major element compositions, which are similar to experimental melts produced by partial melting of amphibolites (Figure 13) and (b) the trend of decreasing P₂O₅ with increasing SiO₂ (Figure 6(f), [34]). (c) In addition, they show meta- to peraluminous geochemical signatures (A/CNK values of SE Mashhad pluton: 0.86–1.16; NW pluton: 0.94–1.15). However, their enriched Sr-Nd isotope composition and abundance of inherited zircon cores are consistent with S-type magmas produced by the melting of metasedimentary rocks [34]. The transitional compositions are similar to granitic rocks of the Harcourt batholith in southeastern Australia, where heterogeneous source rocks varying from meta-igneous to mildly aluminous metasedimentary rocks [57] were suggested. Alternatively, the origin of I-type magma has been explained by hybridization between melts from metasedimentary rocks and interlayered mafic horizons [58] Therefore, partial melting of heterogeneous crustal rocks could be applicable to the generation of Late Triassic granodiorites in the NW and SE Mashhad areas.

For zircons from three granodiorites from the NW Mashhad, ten spots on syn-magmatic growth zones yielded ages ranging from 235.0 ± 5.5 to 296.2 ± 10.8 Ma. Similarly, eleven spots on magmatic growth zones on zircons from two granodiorites of the SE Mashhad pluton yielded apparent ages ranging from 214.8 ± 3.1 to 293.3 ± 10 Ma and one analysis of 413.9 ± 38 Ma. For these outliers, however, five repeated analyses on syn-magmatic growth zone for five zircon yielded apparent ages that agree well with their corresponding sample’s weighted mean age within uncertainty. In addition, for the spot with age of 413.9 ± 38 Ma from sample MH28, the laser sampling location straddles the new magmatic growth zone and inherited core with age of 752.3 ± 21 Ma (Figure 2). Based on the CL images and repeated analyses results, analyses on magmatic zircon rims with apparent age older than the corresponding mean age are most likely result from sampling overlap between different zones and thus meaningless. So, they have been excluded from the following discussion.

Inherited zircon cores with syn-magmatic overgrowths are common in the granodiorites from the NW and SE Mashhad plutons. The youngest ages for the inherited zircon cores from each sample are Devonian or older. The inherited cores could be derived either from (a) crustal contamination during emplacement or (b) the source from which the host magma is generated.

Lines of evidence indicate that crustal contamination during magma ascending through the crust is unlikely. First, granodiorites from the NW and SE Mashhad intruded directly into slate with a whole-rock εNd(t) value of −15 [15]. Significant crustal contamination would result in higher Th/La and lower εNd(t) with increasing SiO₂ for melts, which is not present in these samples (Figures 11(a) and 11(b)). In general, a melt should display lower Nb/La and Nb/Ta with increasing initial ⁸⁷Sr/⁸⁶Sr values as crustal contamination occurs, but the expected trend is absent (Figures 11(c) and 11(d)). Similarly, the relatively homogeneous zircon εHf(t) values for individual samples are not consistent with significant crustal contamination, and the moderate variation in zircon εHf(t) was similar to those of zircons from the leucosome of the Hidaka migmatitic mafic granulites, in northern Japan, which reflect heterogeneity in source regions [58]. Furthermore, the Late Triassic granitoids in Mashhad intruded the Paleo-Tethyan suture, possibly southern margin of the Turan terrane, which could not provide the dominant components with the Early Paleozoic to Late Neoproterozoic ages. Because Early Mesozoic sedimentary rocks sourced from the Turan terrane were characterized by abundant detrital zircon age populations of ca. 1800–1900 Ma (Figure 4(f), [47]), this age population is rare in our samples. To the north of the Paleo-Tethyan suture zone, magmatism ranging from ca. 400 to 240 Ma is present in the Turan terrane [19]. The Triassic Nakhlak Group sandstones were sourced from magmatic arc related to the subduction of the Paleo-Tethyan slab, which has a detrital zircon prominent age group of 240–280 Ma [55]. However, no inherited zircon component with such ages has been revealed from these granitic rocks. In summary, crustal assimilation during the generation of Late Triassic magmas in Mashhad is negligible, although cannot be excluded completely.

Given the weak evidence for crustal contamination, we propose that the inherited zircon cores with syn-magmatic overgrowth zone were derived mostly from the magma source region. Inherited zircons from intermediate to felsic rocks have been suggested to be derived predominantly from their source materials [12]. For example, the Eocene to Miocene Southern Mountains Arc, eastern Java, consists of andesite to dacite, and Archean to Cambrian zircons in those rocks have been explained by partial melting of underlying continental crust beneath the arc [11]. A Cambrian source has also been suggested for ca. 430 Ma granites in southeastern Australia supported by the presence of a large number of Cambrian inherited zircons [59].

Major age peaks of the inherited zircon cores in granitic rocks in the Mashhad area largely reflect magmatism in the Central Iran terrane. The dominant age peak at ca. 540 Ma is consistent with the age of Late Neoproterozoic–Early Paleozoic Cadomian magmatic arcs in Iran and Anatolia, and they were formed as a result of southward subduction of Proto-Tethys underneath the northern margin of Gondwana [21]. On the zircon εHf(t)-age plot (Figure 10), Neoproterozoic-inherited zircon cores have εHf(t) values similar to ca. 540 Ma magmatic rocks in the central Iran terrane and detrital zircons, and this implies that the inherited zircon cores were derived principally from Neoproterozoic–Early Paleozoic magmatism. In addition, Triassic zircon overgrowth and their inherited cores in NW and SE Mashhad intrusions overlap the zircon Hf isotopic evolution range of detrital zircons from Neoproterozoic–Cambrian formations [54] and Neoproterozoic magmatic rocks [60, 61]; this suggests that the parent magmas were produced from a common continental source region. While igneous rocks with ages older than ca. 600 Ma are scarce in the Central Iran terrane, only tuffaceous rocks from the Tashk Formation have zircon ages ranging from ca. 624 to 805 Ma [62]. The sandstones of the Lower Cambrian Zaigun Formation exposed in the northwestern part of the Central Iran terrane exhibited a zircon age cluster at 700–900 Ma [54] and could be sourced from exhumation and erosion of the early Neoproterozoic magmatic rocks in the Central Iran terrane [20]. In addition, detrital zircons of roughly similar ages are abundant in Early Paleozoic clastic rocks [20, 61]. The available data suggest that the Central Iran terrane is a suitable source region to contribute to the production of Late Triassic magmas in Mashhad area.

Lines of evidence below indicate that Paleo-Tethyan oceanic lithosphere did not detach from the Central Iran terrane. In NE Iran, magmatism following the Late Triassic pulse is negligible [63]. This argues against the detachment of the Central Iran continental lithosphere from the early oceanic slab [64]. The low-grade metamorphism and deformation of the Mesozoic rocks in NE Iran [27, 28] suggest continued convergence between the Turan and the Central Iran terrane, which likely resulted from drag of the early subducted Paleo-Tethyan lithosphere. Thus, the subduction of the Central Iran continental crust was supported by pull force from the descending oceanic slab (Figure 14(a)).

Enrichment of the Paleo-Tethyan mantle in the Late Triassic supports the subduction of the Central Iran continental crust into the mantle. The least evolved mafic microgranular enclave (MgO of 5.7 wt%) in the SE Mashhad granodiorite has a whole-rock εNd(t) value of −1.4 [13] lower than midocean ridge basalts of the Iranian Paleo-Tethys [17]; this indicates that the Late Triassic lithospheric mantle has been enriched by subduction materials. The Devonian Darrehanjir plagiogranite has arc-like trace element patterns and a high zircon εHf(t) value (mean of +15) and is considered to be the melts of depleted mantle [17]. In contrast, zircon εHf(t) values of the Late Triassic gabbros from the Torbat-e-Jam pluton are lower than the Devonian plagiogranite, suggesting that it formed by melting subarc mantle modified by subduction components [14]. These observations indicate that the postcollisional Paleo-Tethyan mantle was at least partly enriched by subducted crustal materials, for example, the Central Iranian continental crust. Heat required for partial melting of the subducted Central Iran continental crust might be provided by the enhanced asthenosphere corner flow induced by roll-back of the Paleo-Tethyan lithosphere ([23, 65]; Figure 14), as indicated by Triassic alkaline basaltic rocks of northern central Alborz [66].

The moderately high La/Yb and Sr/Y values of granodiorites in SE and NW Mashhad areas indicate partial melting under relatively high-pressure conditions (~0.8–1.5 GPa) approximate to the stability of garnet [67]. There is no geological or geophysical evidence for crustal thickening before and during the continental collision [13, 16], and thus, this pressure condition might reflect the depth of lower crust to the crust-mantle boundary, possibly the upper section of the subducted Central Iranian continental lithosphere. Furthermore, (a) granodiorites in NW and SE Mashhad show whole-rock εNd(t) values lower than granodiorites of the Torbat-e-Jam pluton, which were considered as melts derived from the lower crust of the Turan terrane (this study [14]) and (b) these granodiorites intruded the Turan terrane, north of the Paleo-Tethyan suture. Collectively, this suggests that they were probably derived from partial melting of the subducted Central Iran continental crust.

The occurrence of granitic rocks intersecting terrane boundaries is not unusual, generally known as stitching plutons [68-70]. Granodiorites in NW and SE Mashhad areas are stitching plutons in NE Iran [47]. In general, magma ascent is facilitated by positive buoyancy. And thus, magmas derived from the subducting continental crust ascending upward, crossing the overriding crust, and even the suture of amalgamated terranes could be a more straightforward and suitable explanation. Therefore, like other stitching plutons in the world, the formation of Late Triassic magmas in the Mashhad area was associated with partial melting of the underthrusting continental crust, which could be more common than previously thought.

Melts of continental crust have been invoked for the genesis of postcollision mafic magmas that are enriched in elements that are incompatible in the mantle (such as Rb, Th, and U), LREE, and Sr-Nd-O isotopes [6, 7]. However, continental crust subduction was generally only recognized in collisional orogeny associated with (ultra)high-pressure metamorphism [1], whereas the presence of high-pressure metamorphic event is rare and limits the role played by continental crust subduction to mantle enrichment. Our study provides evidence that continental crust subduction occurs in collisional orogen without high-pressure metamorphism to data. Taking into account the larger length of the common collisional belt relative to the recorded (ultra)high-pressure metamorphism [1, 2], this would suggest that subduction and melts of continental crust contribute to mantle metasomatism at a larger spatial scale.

Late Triassic postcollision (ca. 217–200 Ma) magmas in NE Iran provide new insights into magma generation and continental crust recycling.

1. The inherited zircon cores in NW and SE Mashhad plutons likely survived melting of their source region.

2. Adakitic rocks in the NW and SE Mashhad areas were likely produced by partial melting of the subducted Central Iran continental crust.

3. Partial melting of subducted continental crust can take place at collisional orogenic belt without occurrence of ultrahigh-pressure metamorphism.

The authors declare that they have no conflicts of interest.

We are benefited greatly by incisive comments and helpful suggestions from three reviewers and an editor on early version of this manuscript. We thank Mr. Hosein Hadizadeh from the Geological Survey of Iran for his assistance in the field expedition. We appreciate Yang Shuqin, Tang Yanwen, Hu Jing, Xiao Fang, Dai Zhihui, and Dong Shaohua for their help during laboratory analyses. Financial support was conjointly provided by the National Natural Science Foundation of China (42222206, 42073047, 91955209, and 42121003) and the projects of “Innovative Team of One Belt and One Road,” K.C. Wong Education Foundation (GJTD-2020-13), and Hundred Talent Plan of the Chinese Academy of Sciences.