Metamorphic and intrusive history of the Hindu Raj region, northern Pakistan

Abstract

The Hindu Raj region of northern Pakistan is situated between the Karakoram to the east and the Hindu Kush to the west. Both the Karakoram and the Hindu Kush are better studied and have well-documented, distinct geological histories. Investigation of the Hindu Raj region has been mainly limited to reconnaissance exploration and as such little is known about its tectonometamorphic history and whether that history is similar to its neighbouring areas. Analysis of new specimens collected along the Yasin Valley within the Hindu Raj region outline mid-to-Late Cretaceous pluton emplacement (ca. 105 and 95 Ma). Some of those plutonic rocks were metamorphosed to ∼750 ± 30 °C and 0.65 ± 0.05 GPa during the ca. 80–75 Ma docking of the Kohistan arc. A record of this collisional event is well-preserved to the west in the Hindu Kush and variably so to the east in the Hunza Karakoram. A subsequent, ca. 61 Ma, thermal event is partially preserved in Rb–Sr geochronology from the Hindu Raj, which overlaps with sillimanite-grade metamorphism in the Hunza portion of the Karakoram region to the east. Finally, apatite U–Pb and in situ Rb–Sr both record a late Eocene thermal/fluid event likely related to the India-Asia collision. These new data outline a complex geological history within the Hindu Raj, one that shares similarities with both adjacent regions. The information about the tectonometamorphic development of the Hindu Raj is important to gaining a detailed view of the geological characteristics of the southern Asian margin prior to the India-Asia collision.

1. Introduction

Accretionary orogens are among the most complex tectonic environments on Earth. They are structurally segmented (e.g. van Staal et al. 2012), lithotectonically heterogeneous (e.g. Monger & Price, 1979), often host to multiple intrusive suites (e.g. Li et al. 2016) and record various deformation ± reactivation events (e.g. Piette-Lauzière et al. 2020). The geological histories of such environments are further obfuscated when overprinted by a terminal continent–continent collision (e.g. Hatcher, 2010).

The southern Asian margin developed as a long-lived accretionary orogen throughout much of the Mesozoic and earliest Cenozoic eras following the breakup of Gondwana and prior to the final closure of the Neo-Tethys ocean (Green et al. 2008) and formation of the Himalayan orogen (Yin & Harrison, 2000). Rocks of Asian affinity that now dominantly comprise the Tibetan plateau and adjacent portions of the Pamir and Northern Pakistan record a history of terrane accretion along and between a series of sutures ranging from early Paleozoic (Qilian) to Jurassic (Bangong) age (Tapponnier et al. 2001; Zanchi et al. 2001). The southern margin of the Asian continent was also the locus of subduction from Late Jurassic to Eocene time where the large-scale, Andean-style continental Ladakh-Gangdese batholith and calc-alkaline volcanic arc developed over 1500 km along the margin (e.g. Schärer et al. 1984). This protracted history of accretion, deformation and magmatism preconditioned the margin prior to its collision with India, laying the foundation for the lithospheric response to subsequent orogenesis. It follows then, that a full understanding of the history of southern Asian rocks is a prerequisite for interpretations about the orogen commonly viewed as the type example of continent–continent collision.

Geological research has a long history across the Tibetan plateau, the Karakoram (Searle et al. 1990; Fraser et al. 2001; Rolland et al. 2001; Rolland et al. 2006a, 2006b; Searle et al. 2010; Zanchi & Gaetani, 2011; Palin et al. 2012) and the Hindu Kush (Gaetani et al. 1996; Hildebrand et al. 1998; Zanchi et al. 2000; Hildebrand et al. 2001; Faisal et al. 2014; Faisal et al. 2016; Faisal et al. 2018; Soret et al. 2019) and has provided much insight into the pre-Himalayan (>55 Ma) history of the south Asian margin. There remain, however, some portions of the former margin that have not yet been investigated in detail leaving significant gaps in the continuity of knowledge about the southern Asian antecedent. One such place is the Hindu Raj region of north central Pakistan, which encompasses the Yasin and Ishkuman valleys, situated between the Hindu Kush region to the west and the Karakoram region to the east (Fig. 1a). The main geological reports from the area comprise reconnaissance mapping commissioned through the Geological Survey of Pakistan, completed in 1977 and published in 1988 (Ebblin, 1988), preliminary geological observations to the immediate west of the main Yasin valley by Le Fort and Gaetani (1998), local geological observations across northern Pakistan (Jan et al. 1981) and initial K-Ar (Casnedi et al. 1978) and whole Rb–Sr and geochemical data from the broader region (Debon et al. 1987).

Figure 1. (a) Simplified regional tectonic map after Robinson (2015) showing the distribution of tectonometamorphic blocks/terranes (i.e. Kohistan, Karakoram, Southern Pamir/Hindu Kush, Central Pamir, North Pamir and Kunlun) and intervening structures. Inset shows approximate location of the map within the larger geographic area. Dashed black line with white background shows the approximate location of the map in B. HK – Hindu Kush region; NZ – Hunza region; BT – Baltoro region. (b) Simplified geological map of northern Pakistan modified after Searle and Khan (1996) showing sampling locations and structural measurements.

In this contribution, we present the results of a new geological investigation of specimens collected along the Yasin valley in the Hindu Raj region. This work includes the first monazite U(-Th)-Pb petrochronology, phase equilibria models, apatite U–Pb geochronology and biotite in situ Rb–Sr geochronology from metamorphic rocks the area and the first U–Pb zircon geochronology and trace element geochemical information from plutonic rocks. These new data help close a gap in our understanding of the development of southern Asia in northern Pakistan, between the better studied Hindu Kush and Karakoram regions, and provide additional detailed information about the evolution of the region prior to – and during – India-Kohistan-Asia collision.

2. Geological background

Here we provide a brief summary of previous work in the Hindu Kush, Karakoram and Pamir regions that surround the present study area.

2.a. Hindu Kush region

The Hindu Kush region is mainly underlain by rocks of the Eastern Hindu Kush block (Zanchi et al. 2000; Zanchi & Gaetani, 2011), also variably referred to as the Wakhan zone (Gaetani et al. 1996) or part of the Southern Pamir/Hindu Kush terrane (Fig. 1a; Robinson, 2015), which was accreted to the south Asian margin in the Late Triassic Epoch (Faisal et al. 2014). The rocks within the Eastern Hindu Kush block are of Gondwanan affinity (Angiolini et al. 2013) and typically comprise argillaceous and calcareous protoliths (Zanchi & Gaetani, 2011) that have been variably metamorphosed to sillimanite-bearing migmatite (Hildebrand et al. 2001). The tectonometamorphic record preserved in the rocks of the Eastern Hindu Kush block includes ca. 185 Ma accretion of the Karakoram block, which reached pressure-temperature (P-T) conditions of 0.3 GPa and 400–500 °C and the Late Cretaceous docking of the Kohistan island arc (0.6 GPa and 500–600 °C) (Faisal et al. 2014; Soret et al. 2019). Cenozoic collision-related sillimanite grade rocks reached P-T conditions of 0.7–0.8 GPa and 700–750 °C in the late Eocene-early Oligocene while those at staurolite grade experienced 0.6 GPa and 580–600 °C at the same time (Soret et al. 2019). The region is also host to a series of rift-related structures and intrusions (ca. 500–480 Ma – Kafiristan) and pre-collision subduction-related granodiorites (ca. 123–127 Ma – Tirich Mir). Localised post-collision anatectic leucogranites are seen in the ca. 24 Ma Gharam Chasma biotite + muscovite ± garnet ± tourmaline leucogranite (Hildebrand et al. 1998; Hildebrand et al. 2001; Heuberger et al. 2007; Faisal et al. 2016). Gharam Chasma crystallization is synchronous with the early stages of the large Baltoro granite batholith in the Central and Eastern Karakoram (Rex et al. 1988; Searle, 1991; Searle et al. 1992; Searle et al. 2010).

The Eastern Hindu Kush block is separated from the adjacent Karakoram block to the south by the Tirich Mir boundary zone (Hildebrand et al. 2000), which is marked by peridotite and metagabbro (Zanchi et al. 2000). The rocks of the Karakoram block to the south-east of the Tirich Mir boundary zone in the Hindu Kush region can be divided into the Karakoram Northern sedimentary belt and Southern sedimentary belt, which comprise low metamorphic grade slate + limestone and argillaceous + siliciclastic + carbonate-rich phyllite/chert/breccia, respectively (Hildebrand et al. 2000). Like those in the Eastern Hindu Kush block, the rocks of the Karakoram block also host 110–108 Ma subduction-related plutons (110–108 Ma Buni Zom; Faisal et al. 2016; ca. 104 Ma Phargam; Heuberger et al. 2007).

2.b. Karakoram region

The Karakoram region to the east of the Hindu Raj contains rocks of the Karakoram block. The Hunza and Baltoro (Fig. 1a) portions of the Karakoram region record complex, protracted and spatially distinct histories of high-grade metamorphism and magmatism. The broader region is host to the >700 km long Karakoram batholith, which comprises pre-collisional subduction-related granodiorites and granites, some metamorphosed to amphibolite-facies, that range in age from ca. 200 to 80 Ma (Searle et al. 1990 and references therein; Crawford & Searle, 1992). In the Baltoro region, U–Pb dating of granitoids indicates that the Baltoro granite component of the Karakoram batholith intruded between ca. 27 and 13 Ma (Searle et al. 2010). Regional low-pressure andalusite-biotite ± sillimanite metamorphism associated with pluton emplacement is interpreted to have occurred in the southern Baltoro Karakoram and Hunza valley region at ca. 106 Ma, followed by regional sillimanite-grade metamorphism at ca. 82–44 Ma (Fraser et al. 2001; Foster et al. 2002; Foster et al. 2004) and a kyanite-grade metamorphic overprint at ca. 28 Ma (0.78 GPa and 645 °C) (Palin et al. 2012). This contrasts with post ca. 21.8 Ma kyanite-grade (0.74–0.80 GPa and 640–660 °C) metamorphism recorded in the Baltoro region that was driven by conductive heating following post-collisional crustal thickening (Searle & Tirrule, 1991; Palin et al. 2012).

The Dassu-Askole region of the Karakoram, to the south and east of the Baltoro region (Fig. 1a), contains amphibolite facies metamorphism (0.7 GPa and 700 °C) associated with Paleogene (55–37 Ma) crustal thickening (Rolland et al. 2001). A younger, higher temperature but lower pressure (∼775–850 °C and 0.5–0.6 GPa; Rolland et al. 2001; Rolland et al. 2006b) event is also recorded in domal culminations (Dassu and Askole) in the area (Searle, 1991). Monazite from sillimanite-bearing orthogneiss in the core of the Dassu dome are as young as ca. 5.4 Ma (Fraser et al. 2001; Searle et al. 2010), while zircon and monazite from a tourmaline pegmatite are even younger at ca. 3.5 Ma (Searle et al. 2010).

2.c. Pamir

The rocks of the Southern Pamir appear to be related to those in the Karakoram (Angiolini et al. 2013) and are separated from the North and Central Pamir by the Rushan – Pshart suture zone (Fig. 1a; Schwab et al. 2004; Robinson et al. 2012; Hacker et al. 2017; Searle & Hacker, 2019). A series of thrust sheets and metamorphic domes, such as the large Shakhdara dome, record crustal thickening and regional metamorphism from 37 to 20 Ma, concomitant with thickening and metamorphism along the Karakoram (see Hacker et al. 2017; Searle & Hacker, 2019 and references therein). Late strike-slip faults, such as the Kilik fault (Zanchi & Gaetani, 2011), the possible eastern extension of the Tirich Mir boundary zone (Hildebrand et al. 2000; Zanchi et al. 2000; Robinson, 2015) cut through the Hindu Kush and southern Pamir during Cenozoic indentation of the Pamir.

The rocks of the Southern Pamir, like those of the Hindu Kush and the Karakoram, comprise a basement of Gondwanan affinity (Angiolini et al. 2013) that host subduction-related Cretaceous plutonic rocks (Schwab et al. 2004; Stübner et al. 2013). The metamorphic record in the region is distinct and dominated by rocks deeply exhumed within large domal culminations where conditions of 0.7–1.5 GPa and 700–800 °C were reached (Schmidt et al. 2011; Hacker et al. 2017). Titanite dates from amphibolite within the Shakdhara dome, the largest dome in the southern Pamir, range between 34 ± 4.6 and ca. 10 Ma, which overlaps with garnet Lu-Hf dates from the same lithology (ca. 37 Ma) and monazite from proximal paragneiss and leucocratic gneiss (30–18 Ma)(Stearns et al. 2015; Hacker et al. 2017).

2.d. Hindu Raj

The Hindu Raj range is situated between the Hindu Kush region to the west and Karakoram region to the east. Two major valleys transect the range, the Yasin River in the west which rises from the Darkot pass and the Ishkuman – Karambar River which extends north as far as the Wakhan corridor (Afghanistan) (Fig. 1b). A northern sedimentary belt is dominated by a thick, mainly carbonate sequence ranging from Permian to Cretaceous in age (Casnedi, 1979; Ebblin, 1988; Zanchi et al. 2001). A central granitoid belt extends from the Darkot pass east to the Hunza Karakoram and comprises both pre-collision granodiorites and monzogranites and post-collision leucogranite dykes (Debon et al. 1987). A distinct alkali granite, including amphibole and pyroxene-bearing syenite and quartz syenite, approximately 20 km, long has been mapped along the Karakoram batholith in the Karambar valley (Koz Sar complex; Debon & Khan, 1996) (Koz Sar complex; Debon & Khan, 1996). These rocks are similar to the Hemasil syenite in the southern Baltoro Karakoram (Mahéo et al. 2009) and are geochemically related to a suite of lamprophyre dykes intruding the Baltoro Karakoram region (Rex et al. 1988; Searle, Crawford & Rex, 1992; Searle et al. 2010). The Koz Sar monzonites have a Rb–Sr isochron age of 88 ± 4 Ma (Debon & Khan, 1996). If this age marks an intrusive age, then alkaline magmatism accompanied the calc-alkaline granites of the Hunza plutionc unit. The Koz Sar syenites are also petrologically similar to a suite of alkaline intrusives along the southern Pamir region (Searle & Hacker, 2019). These syenites and lamprophyres in the Baltoro and Pamir regions have been interpreted as mantle-derived melts intruding at approximately the same time as the Baltoro-type crustal melt granites (Searle et al. 1992; Searle et al. 2010).

An additional granitic batholith, the subduction-related Ghamu Bar pluton (Debon et al. 1987) extends into the Hindu Raj region from the west (Fig. 1b) and occurs within variably metamorphosed sedimentary rocks. Metamorphism within these zone reaches sillimanite grade (Ebblin, 1988). The southern margin of the Hindu Raj is marked by the Shyok suture zone which is the northern margin of the Kohistan island arc (Fig. 1b).

The absolute timing of metamorphism in the area is unknown. High-metamorphic grade rocks are observed to be spatially associated with granodiorite, also of unknown age, but inferred to be related to similar rocks farther to the south-west that yield zircon dates between ca. 110 and 102 Ma (Heuberger et al. 2007; Faisal et al. 2016). At least two regional folding events are reported with F1, which records N-S maximal finite strain, dominant over F2 with a near vertical maximal finite strain (Ebblin, 1988).

3. Methods

3.a. Field mapping

Six days of vehicle-supported field mapping along the Yasin and Ishkuman valleys was carried out in the Spring of 2018. Spatial relationships between geological lithologies were recorded and documented along with measurements of geological structures. A suite of representative specimens was collected for further laboratory analyses.

3.b. EPMA

Mineral compositions were obtained at the Fipke Laboratory for Trace Element Research (FiLTER; University of British Columbia Okanagan) using a Cameca SXFiveFE electron microprobe analyser (EMPA). Operating conditions include an excitation voltage of 15 kV, a beam current of 20 nA, a peak count time of 30 s with a background count time of 15 s and a focused spot size of 5 μm. Unfiltered data including the standards, X-ray lines and crystals used are provided in Tables S1 and S2. For all applicable phases, Fe³⁺was calculated on the basis of charge-balance using the software AX (Holland & Powell, 2009).

3.c. Automated mineralogy

The modal abundances of phases in thin section were determined by an automated mineralogy routine that integrated backscattered electron (BSE) and energy-dispersive X-ray (EDS) signals (Gottlieb et al. 2000). Thin sections were analysed at Vidence Inc., Burnaby, Canada, using a Hitachi SU3900 scanning electron microscope fitted with a single large area (60 mm²) Bruker SDD energy-dispersive spectrometer and running the AMICS automated mineralogy package. BSE brightness was calibrated against quartz, gold and copper standards.

3.d. P-T determination

Phase diagrams were calculated using the petrological modelling software Theriak–Domino (de Capitani & Petrakakis, 2010). Sample specific bulk compositions (Appendix S2) were constructed by integrating modal phase abundances for the thin sections with mineral chemistry (Palin et al. 2016). Modelling was undertaken in the MnO–Na₂O–CaO–K₂O–FeO–MgO–Al₂O₃–SiO₂–H₂O–TiO₂–O₂(MnNCKFMASHTO) chemical system using the internally consistent thermodynamic dataset ds-55 of Holland and Powell (1998) (updated to August 2004) and activity–composition relations for melt, biotite and ilmenite (White et al. 2007), feldspar (Holland & Powell, 2003), white mica (Coggon & Holland, 2002) and garnet (White et al. 2007). Dataset ds-55 was chosen over the more recent ds-62 because it better predicted the observed mineral assemblages, and in particular, the composition and stability of garnet in our samples (c.f., Dyck et al. 2020; Dyck et al. 2021). The amount of water was set to minimally saturate the assemblage in the immediate sub-solidus at 0.7 GPa. Bulk rock X_Fe³⁺ = Fe³⁺/(Fe_total) was set for each sample and is a function of the Fe³⁺ determined for all ferric minerals by stoichiometric charge balance. Garnet isomodes and end-member isopleths were generated from the output of a 50 × 50 gridded Gibb’s free energy minimization calculation (pixelmap) and contoured in Matlab. Uncertainty on the absolute positions of phase assemblage boundaries in pressure–temperature space for our choice of dataset and activity–composition models generally do not exceed ±20–30 °C and ± 0.05 GPa at 1σ (Powell & Holland, 2008; Palin et al. 2016).

We used the average PT (avPT) function in thermocalc (versions tc333 with ds-55; Powell & Holland, 1994) and the Ti-in-biotite thermometer of Henry et al. (2005) to independently verify the metamorphic conditions determined using the phase diagrams results. The 1σ uncertainty for avPT calculations is expressed as an error ellipse, the dimension of which scales with the collective fit of the end-member reactions (Powell & Holland, 1994). We report all Ti-in-biotite temperatures with ± 24 ºC uncertainty (Henry et al. 2005).

3.e. U(-Th)-Pb geochronology

3.e.1. Metamorphic specimens

3.e.1.a. Monazite

Monazite in selected metamorphic rocks collected during fieldwork was analysed in situ for U(-Th)-Pb geochronology and trace element concentrations in the FiLTER facility via laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) using a Photon Machines Analyte 193 Excimer laser paired with an Agilent 8900 triple quadrupole instrument run in single quadrupole mode. ICP-MS signal was optimized for maximum signal using the glass reference material NIST610 (Jochum et al. 2011) and minimal oxide and doubly charged ion production (<1%), while maintaining ²³⁸U/²³²Th within 3% of target values. A laser spot size of 12 microns was used with a repetition rate of 6 Hz and an estimated fluence of 3.5 J/cm². All spots were pre-ablated with two laser pulses followed by 35 s of background signal collection and 25 s of ablation. Instrument drift through analytical sessions and downhole fractionation was monitored through analysis of matrix-matched reference monazite bracketing every ∼10 unknown laser spots. Analysis and correction for both were performed using Iolite v. 4.5 (Paton et al. 2010, 2011). Excess dispersion for each radiogenic ratio of interest was calculated based on assumed homogeneity of secondary reference monazites that were also analysed bracketing the unknowns. Calculated excess dispersion, typically <1%, was propagated in quadrature through the ratios of all analyses for the associated run. All data are uncorrected for common Pb. The full datasets collected are presented in Tables S2 and S4.

Monazite “44069” – ²⁰⁶Pb/²³⁸U date of 424.9 ± 0.4 Ma (Aleinikoff et al. 2006) – was used as the primary reference during the two monazite analytical sessions with ‘Bananeira’ and ‘554’ monazite analysed as secondary reference materials. Bananeira, also referred to as ‘Stern’, returned ²⁰⁷Pb corrected (Stacey & Kramers, 1975) ²⁰⁶Pb/²³⁸U dates of 507 ± 4 Ma (mean squared weighted deviation – MSWD = 0.39, n = 8/10) and 514 ± 6 Ma (MSWD = 2.1, n = 16/16), which overlap reported dates of 512 ± 2 Ma (Palin et al. 2013) and 508 ± 1 Ma (Kylander-Clark et al. 2013). The ²⁰⁸Pb/²³²Th dates returned from the two analytical runs of Bananeira are 493 ± 4 Ma (MSWD = 1.8, n = 8/10) and 508 ± 4 Ma (MSWD = 1.5, n = 14/16), respectively. Both dates are within the estimated reproducibility of the laser ablation method (2%; Horstwood et al. 2016) of a reported, but overdispersed ²³²Th/²⁰⁸Pb date (LA-ICP-MS) of 497 ± 2 Ma (MSWD = 6.1; Kylander-Clark et al. 2013). Monazite 554 returned ²⁰⁷Pb corrected (Stacey & Kramers, 1975) ²⁰⁶Pb/²³⁸U dates of 45 ± 3 Ma (MSWD = 1.4, n = 9/10) and 49.8 ± 1.8 Ma (MSWD = 1.6, n = 15/16) and ²⁰⁸Pb/²³²Th dates of 46.9 ± 0.4 (MSWD = 0.31, n = 10/10) and 46.2 ± 0.3 Ma (MSWD = 2.5, n = 14/16). All dates fall within the range of published dates including a ²⁰⁸Pb/²³²Th date of 45.0 ± 1.3 Ma (M. Tatsumoto, as cited in Harrison et al. 1999), and ²⁰⁶Pb/²³⁸U dates of ∼44.3–50.8 Ma (multi-step thermal ionization mass spectrometry; Peterman et al. 2012).

Monazite trace element data were collected from the same ablated volume as the U-Th-Pb isotopes. Bananeira was used as the matrix-matched primary reference material for trace element concentrations using the values reported in Kylander-Clark et al. (2013) with stoichiometric P employed as the internal calibrant for the unknown analyses. Previous studies have estimated a reproducible accuracy of ∼5% for repeat analyses of secondary reference monazites (Cottle et al. 2019). In the present work, we do not use the absolute values of trace elements for interpretations, but instead examine the relative changes between analyses

3.e.1.b. Apatite

Analysis of apatite for U–Pb geochronology was also carried out in FiLTER using the same instrumentation as the monazite analysis. It follows the basic procedure outlined above for instrument tuning, unknown bracketing and drift/downhole fractionation corrections. A laser spot size of 40 microns at a repetition of 8 Hz and an estimated fluence of 4 J/cm² was used for all analyses. The apatite MCR1 (153.3 ± 0.2 Ma isochron date; Apen et al. 2022) was used as the primary reference apatite, with MAD (Thomson et al. 2012; 487 ± 1 to 467 ± 9 Ma; Apen et al. 2022) and Mount McClure (523.5 ± 1.5 Ma; Schoene & Bowring, 2006) apatite analysed as matrix-matched secondary reference materials. Repeat analysis of MAD returned a ²⁰⁷Pb corrected (Stacey & Kramers, 1975) ²⁰⁶Pb/²³⁸U date of 481 ± 2 Ma (MSWD = 1.4, n = 12/12), while Mount McClure yielded an overdispersed Tera-Wasserburg intercept date of 525 ± 8 Ma (MSWD = 11, n = 11/12). The full apatite U–Pb dataset is presented in Table S5.

3.e.2. Plutonic specimens

A total of four plutonic specimens, Y7, 9, 10 and 12, from the study area were crushed and separated following standard procedures for accessory mineral geochronology and geochemistry. Of those specimens, Y7 returned titanite and zircon (only zircon was targeted), Y9 yielded monazite and zircon and Y10 produced only zircon. No datable accessory minerals were extracted from specimen Y12.

Zircon U–Pb dates and trace element data were acquired at the University of California Santa Barbara using a laser ablation split stream system consisting of a Photon Machines Excimer 193 nm laser ablation unit coupled to a Nu Instruments, ‘Nu Plasma 3D’ multi-collector ICP-MS (U, Th and Pb isotopes) and an Agilent 7700S quadrupole ICP-MS (trace elements), using methods of Kylander-Clark et al. (2013) with modifications as outlined in McKinney et al. (2015). Samples were analysed for 20 s using a fluence of 1.5 J cm^–2, a frequency of 4 Hz and spot size of 20 μm diameter resulting in crater depths of ∼9 μm. Using a standard-sample bracketing technique, reference zircons with known isotopic compositions were analysed before and after each set of eight unknown analyses. Data reduction, including corrections for baseline, instrumental drift, mass bias, downhole fractionation, age and trace element concentration calculations were carried out using Iolite v. 4.1 (Paton et al. 2010). Zircon U–Th/Pb and trace element data were normalized to 91500 zircon (1062.4 Ma 206Pb/238U ID-TIMS age; Wiedenbeck et al. 2004 – trace element values from GeoREM as of December 2022). A ²⁰⁷Pb-based correction method was utilized using a common-Pb composition derived from the single-stage model of Stacey and Kramers (1975) at the inferred crystallization age. The uncertainty on the ²⁰⁷Pb corrected age incorporates uncertainties on the measured ²⁰⁶Pb/²³⁸U and ²⁰⁷Pb/²⁰⁶Pb ratios as well as a 1% uncertainty on the assumed common lead composition. Plešovice (Sláma et al. 2008), with an expected ²⁰⁶Pb/²³⁸U date of 337.13 ± 37 Ma, and GJ1 (Jackson et al. 2004), with an expected ²⁰⁶Pb/²³⁸U date of 601.86 ± 0.37 Ma (Horstwood et al. 2016), were analysed as secondary reference zircons, and returned weighted mean ²⁰⁶Pb/²³⁸U dates of 337 ± 1 Ma (MSWD = 0.62, n = 17/18) and 601 ± 2 Ma (MSWD = 0.75, n = 18/18), respectively. Trace element concentrations for GJ-1 are accurate to within 5% of the reference values. All uncertainties are quoted at the 95% confidence or 2 sigma level and include contributions from the external reproducibility of the primary reference material for the ²⁰⁷Pb/²⁰⁶Pb and ²⁰⁶Pb/²³⁸U ratios. The full zircon U–Pb dataset is provided in Table S6.

3.f. ⁸⁷Rb/⁸⁷Sr geochronology

In situ Rb–Sr geochronology was carried out in the FiLTER facility at UBCO following the procedures outlined in Hogmalm et al. (2017), Redaa et al. (2021), and Rösel and Zack (2022) and those detailed in Larson et al. (2023). An ESL NWR193 laser with a TwoVol3 ablation cell paired with an Agilent 8900 triple quadrupole inductively coupled plasma mass spectrometer (QQQ-ICP-MS) was used for the analyses. The collision cell in the QQQ-ICP-MS was charged with N₂O gas (∼25 ml/min). The ablation cell He carrier gas (0.35 l/min) was mixed with the Ar sample gas (0.9 l/min) through an in-house signal smoothing device before entering the plasma. A laser spot-size diameter of 50 microns was used for all analyses with a repetition rate of 10 Hz and a ∼4 J/cm² fluence; spots were ablated for 25 s followed by a 20 s washout. All spot locations were pre-ablated with two laser pulses to ensure a clean surface.

A specimen bracketing technique was employed to monitor and correct for instrument drift over the analytical session. The pressed nano-powder, Mica Mg (Zack & Hogmalm, 2016; Hogmalm et al. 2017), was analysed as the primary standard every ∼10 unknowns, while analyses of NIST610 glass were used to correct final ⁸⁷Sr/⁸⁶Sr. Two phlogopite specimens (1B and 1O) and the nano-pressed powder Mica Fe (Rösel & Zack, 2022) were analysed as unknowns. Phlogopite 1B (initial ⁸⁷Sr/⁸⁶Sr intercept of 0.7035; Camacho et al. 2020) yields an isochron date of 978 ± 15 Ma (MSWD = 0.4, n = 7/11). Phlogopite 1O was collected from the same outcrop as 1B and returned an isochron date of 986 ± 18 (MSWD = 0.1, n = 11/11). Both dates overlap the two-point Rb–Sr isochron of 986 ± 5 Ma reported in Camacho et al. (2020) for 1O and 990 ± 5 Ma reported in Camacho et al. (2012) for 1B. Analyses of the pressed nano-powder Mica Fe yield an isochron date of 303 ± 3 Ma (initial ⁸⁷Sr/⁸⁶Sr intercept of 0.7260; Rösel & Zack, 2022 – MSWD = 0.26, n = 10/10), which overlaps with a reported weighted mean of Rb–Sr spot dates of 305 ± 2 Mica Fe (Rösel & Zack, 2022). The full Rb–Sr dataset is included in Table S7.

4. Results

4.a. Field observations

The northern most part of the study area is dominated by intercalated phyllite and psammite (Fig. 2a). Metamorphic grade increases to the south towards the Ghamu Bar pluton (Fig. 1b) where the rock outcrops comprise orthogneiss, quartzofeldspathic gneiss (Fig. 2b) and amphibolite. Metamorphism near pluton contact reaches sillimanite-grade (see description of specimens Y5 and Y3DD below).

Figure 2. Representative field photographs from the Yassin valley. (a) Intercalated phyllite and meta-sandstone. Geology hammer circled for scale. (b) Quartzo-feldpathic gneiss from just north of the Ghamu Bar pluton. Geology hammer circled for scale. (c) Outcrop of the Ghamu Bar pluton. Geologist outlined for scale. Representative photographs of the Ghamu Bar pluton showing unfoliated material (top) from the pluton interior and foliated material (bottom) from near the pluton margin. Coin circled in the top photo for scale. Thumb visible in the lower image for scale. (e) An outcrop of slate from between the Ghamu Bar pluton and Shyok suture; geologist outined for scale. (f) Low metamorphic grade sedimentary rocks between the Ghamu Bar pluton and Shyok suture folded into an antiform and cut by a south-verging reverse fault.

Near its contact with the country rocks, the Ghamu Bar pluton is massive in character (Fig. 2c) with a moderate foliation (Fig. 2d bottom). Farther towards the pluton’s center the rocks are not foliated (Fig. 2d top). Large (centimetre-sized) K-feldspar phenocrysts are found in both the foliated margin and undeformed centre of the pluton (Fig. 2d).

Rocks south of the Ghamu Bar pluton comprise a sedimentary succession affected by sub-greenschist facies metamorphism. The lithologies are dominated by slate (Fig. 2e), quartzite and mixed siliciclastic packages. Moreover, the structural characteristics of the rocks south of the pluton are dominated by meso-scale folding and faulting; both of which typically verge to the south (Fig. 2f).

The geology south of the Shyok suture (Fig. 1b) was not examined in detail. A few specimens of plutonic rocks (see descriptions below) were collected to compare with those collected from the Ghamu Bar pluton.

4.b. Petrology

4.b.1. Metamorphic Specimens

Two specimens, Y5 and Y3DD, were chosen to determine the metamorphic history of the region based on their higher variance mineral assemblages. Specimen Y5 is an orthogneiss with plagioclase (∼69 %), biotite (∼20 %), quartz (∼5 %), garnet (∼ 4%) and minor sillimanite, ilmenite, Fe-oxide, cordierite/white mica and chlorite (Fig. 3a). Based on a 11-oxygen calculation, biotite contains Ti = 0.19 apfu (atoms per formula unit) and has a XMg = Mg/(Mg+Fe²⁺) = 0.42 (Table 1). Garnet are subhedral and contain inclusions of quartz that do not define an internal fabric (Fig. 3a, c). EMPA line profiles across garnet reveal a flat chemical profile, with a composition of Alm_0.76Prp_0.18Sps_0.03Grs_0.03 (Table 1), for which Alm = almandine = Fe²⁺/(Fe²⁺+Mg+Ca+Mn), Prp = pryrope = Mg/(Fe²⁺+Mg+Ca+Mn), Sps = spessartine = Mn/(Fe²⁺+Mg+Ca+Mn) and Grs = grossular = Ca/(Fe²⁺+Mg+Ca+Mn). White mica-rich pseudomorphs after cordierite form lensoidal domains that are aligned with the gneissosity (Fig. 3d).

Figure 3. Phase maps and representative photomicrographs for metamorphic specimen Y5 (a, c, d) and Y3DD (b, e–f). Mineral abbreviations after Whitney and Evans (2010). (c) Garnet wrapped by biotite–sillimanite folia. (d) Pinitized cordierite poikiloblast. (e) Plagioclase clast surrounded by matrix of quartz and muscovite. (f) Radial and fine-grained clots of muscovite. (g) Fractured garnet with biotite overgrowth. (h) Fine-grained clusters of ilmenite along the rim and cleavage plane of low-Ti biotite.

Table 1. Representative EMPA mineral chemistry

Specimen Y3DD is foliated meta-granite with quartz (∼41 %), plagioclase (∼21 %), white mica (∼16 %) and minor K-feldspar, sillimanite, ilmenite and garnet (Fig. 3b). Clasts of plagioclase are wrapped by ribbons of quartz and aluminous folia of biotite, sillimanite and white mica (Fig. 3e). The aluminous folia are overgrown, in part, by radial clots of muscovite (Fig. 3f). There are few garnet grains present in this specimen, and all of them are highly fractured and intergrown with fine-grained white mica and biotite (Fig. 3g). Rims of biotite are decorated by fine-grained clusters of ilmenite (Fig. 3h). Analysis of biotite cores yields Ti = 0.19 apfu and a XMg = 0.42, based on a 22-oxygen formula unit (Table 1).

4.b.2. Plutonic specimens

A total of four specimens were collected for accessory mineral U–Pb geochronology to provide information about the timing of pluton emplacement and deformation in the region. Specimen Y7 was collected from the Karakoram block, while specimens Y9, Y10 and Y12 were collected from the Kohistan arc, south of the Shyok suture (Fig. 1b).

Specimen Y7 was sampled from the northern margin of the Ghamu Bar pluton (Debon et al. 1987), also referred to as Gamughal (Jan et al. 1981) or Ghamubar (Le Fort & Gaetani, 1998), and comprises a K-feldspar + plagioclase + biotite + quartz granodiorite with accessory zircon and titanite (Fig. 4a). It is characterized by K-feldspar megacrysts up to 2 cm in length. The specimen collected and analysed is foliated; however, the same lithology was also observed where it is unfoliated and intruded locally by a garnet-bearing leucogranite.

Figure 4. Representative photomicrographs of plutonic specimens sampled. PPL – plane polarized light; XPL – crossed polars. Mineral abbreviations after Whitney and Evans (2010).

Specimen Y9, collected from the leucocratic phase of the Gindai pluton (Jan et al. 1981), is a K-feldspar + plagioclase + quartz + biotite leucogranite with accessory monazite and zircon. The main minerals are broadly equigranular and medium to coarse in size; no foliation is apparent (Fig. 4b). Specimens Y10 and Y12 were collected from a larger pluton body south of Y9. Both specimens are granodiorite with a mineral assemblage of plagioclase + amphibole + biotite + quartz (Fig. 4c, d); no foliation is apparent. While Y10 yielded accessory zircon after crushing and separation, Y12 was barren.

4.c. Geochronology

4.c.1. Orthogneiss Monazite U(-Th)-Pb

Monazite in specimens Y3DD and Y5 were analysed for U–Th–Pb geochronology to provide information on the timing of metamorphism in the study area. Nine monazite grains in Y3DD were located and mapped for Y, Th, Ca, Si and U distribution to assess elemental zonation and to help guide laser spot placement. With one notable exception, most grains are roughly equidimensional and 30–50 μm across (Fig. S1). Y and Th zonation within monazite cores is generally patchy with some limited sector zoning; thin (≤5 μm) homogeneously high Y or low Y rims are also evident (Fig. S1). All grains are inclusion free except the largest grain analysed, which contained ∼5 μm diameter xenotime inclusions within its core.

One hundred and two monazite analyses across all 9 grains yield a spread in ²³²Th-²⁰⁸Pb dates from 126 ± 11 Ma to 66.1 ± 3.5 Ma (Fig. 5a) with the bulk of the analyses falling between ca. 85 and 70 Ma. Heavy rare earth element (HREE) concentrations vary with date with older analyses generally having flatter HREE slopes (lower Gd/Yb) and younger dates trending towards steeper HREE slopes (higher Gd/Yb) (Fig. 5c). That trend is followed by relative Eu anomalies, which are highest in the oldest spots and progressively decrease with younger dates (Fig. 4c).

Figure 5. Monazite geochronology summary diagrams. (a and b) are ²⁰⁸Pb/²³²Th versus ²⁰⁶Pb^/238U concordia plots coloured by Y concentrations. (c and d) are chondrite normalized trace element spider plots coloured by ²³²Th/²⁰⁸Pb dates. (e and f) are kernel density estimations of ²³²Th/²⁰⁸Pb dates. Plots were created with the ChrontouR package (Larson, 2020) for the open, R scripting platform.

The monazite in Y5 is equant to elongate (up to 2:1) with long axes between 10 and 70 μm (Fig. S2). They are characterized by core and rim internal structures with core regions variably displaying homogeneous, patchy or sector zoning in Y and Th. Grain rims are between ∼1 and 5 μm wide and are commonly low Y domains (Fig. S2). Thirty-two analyses across 10 monazite grains define a spread in ²³²Th-²⁰⁸Pb dates from 107.4 ± 2.6 Ma to 72.3 ± 1.9 Ma that comprise two main populations – an older population between ca. 100 and 90 Ma and a younger population between ca. 80 and 70 Ma (Fig. 5b). HREE slopes in chondrite normalized spider plots of the monazite data define a general trend of increasing slope (Gd/Yb) with younger date (Fig. 5d). This trend is mirrored by normalized Eu concentrations, which are lowest for the oldest dates and broadly increase with younger dates (Fig. 5d).

4.c.2. Apatite U-Pb

Nineteen apatite grains located in specimen Y5 were analysed for U–Pb geochronology (Fig. S3); no suitable grains were identified in specimen Y3DD. Apatite grains in Y5 are generally equant and range in size from ∼50 to ∼200 μm across, with the majority falling on the smaller size of that spectrum. No zoning was evident in back-scattered electron imaging (Fig. S3). Thirty-seven spot analyses across the apatite grains outline a range of analyses that form a single robust regression (Powell et al. 2020) (Fig. 6a) with a lower intercept of 38.3 ± 4.5 Ma (spine [s] = 0.95) in Tera-Wasserburg space.

Figure 6. Top panel shows a Tera-Wasserburg plot of measured apatite U–Pb data in specimen Y5. Bottom panel depicts Rb–Sr geochronological data measured in plagioclase, muscovite and biotite in specimen Y5. The regressions drawn through the data in both panels were calculated using the robust regression method of Powell et al. (2020), where s = spine.

4.c.3. ⁸⁷Rb/⁸⁷Sr

Specimen Y5 was examined for Rb–Sr geochronology to further assess the thermal history of the specimen. Analyses of brown mica, white mica and plagioclase define a wedge-shaped array of data in a ⁸⁷Sr/⁸⁶Sr versus ⁸⁷Rb/⁸⁶Sr plot that opens away from a single initial ⁸⁷Sr/⁸⁶Sr composition (Fig. 6b). The limited spread in the ⁸⁷Sr/⁸⁶Sr data and the uncertainties in the ratios do not permit separation of plagioclase and white mica into multiple populations. The biotite analyses, however, can be grouped to define the sides of the wedge and thereby the minimum and maximum isochrons permissible in the dataset. The older isochron includes all white mica, plagioclase and a subset of brown mica to define a date of 62.5 ± 0.8 Ma (s = 0.45, n = 23/44), whereas the same white mica and plagioclase analyse, but a different subset of brown mica define a younger isochron at 40.1 ± 1.4 Ma (s = 0.39, n = 10/44).

4.c.4. Plutonic zircon and monazite U(-Th)-Pb

Zircon grains separated from specimens Y7, Y9 and Y10 were dated to provide information about the timing of pluton emplacement on both sides of the Shyok suture between rocks of Karakoram (Y7) and Kohistan (Y9, Y10) affinity. Forty spots across the same number of grains (Fig. S4) separated from Y7 define a single age population with a ²⁰⁷Pb corrected (Stacey & Kramers, 1975) ²⁰⁶Pb/²³⁸U weighted mean date of 94.9 ± 0.3 Ma (MSWD = 1.5, n = 40) (Fig. 7). Similarly, 36 spot zircon analyses in Y9 (1 per grain) and 40 spots in 36 zircon from Y10 (Fig. S4) also define single age populations (Fig. 7). The ²⁰⁷Pb corrected (Stacey & Kramers, 1975) ²⁰⁶Pb/²³⁸U weighted mean zircon date for Y9 is 60.8 ± 0.2 Ma (MSWD = 1.7, n = 32/36), while that from Y10 is 60.2 ± 0.2 Ma (MSWD = 1.0, n = 39/40).

Figure 7. Summary of zircon geochronological data for specimens Y7, Y9 and Y10. The top row comprises plots of the measured isotope ratios in Tera-Wasserburg space with ellipses coloured by Y concentration. The middle row of plots depicts weighted mean plots of ²⁰⁷Pb corrected (Stacey & Kramers, 1975) ²⁰⁶Pb/²³⁸U dates for each specimen. The bottom row presents chondrite normalized trace element spider plots for each spot analysis coloured by ²⁰⁷Pb corrected (Stacey & Kramers, 1975) ²⁰⁶Pb/²³⁸U date. All plots were created with the ChrontouR package (Larson, 2020) for the open, R scripting platform.

Zircon trace element concentrations are similar across the three specimens analysed. REE spider plots of data from Y7, Y9 and Y10 show relative depletion in light REE (LREE), with some LREE elements below detection (see supplementary materials), and enrichment in HREE (Fig. 7).

5. Interpretations

The rocks exposed in the Yasin-Ishkuman region record a protracted history of plutonism, metamorphism and post-kinematic metasomatism. The calcic bulk composition and relict igneous microstructures in Y5 and Y3DD attest to both specimens having been derived from silicic plutonic protoliths. Whereas both Y5 and Y3DD preserve a peak metamorphic mineral assemblage with garnet–sillimanite–biotite–plagioclase–quartz–ilmenite (+melt), Y5 also preserves a record of what is likely the subsequent growth of cordierite during decompression from peak-temperature conditions within the stability field of garnet. The replacement of cordierite by white mica (Dunkley et al. 1999) in Y5 and the radial growth of muscovite and growth of ilmenite surrounding lathes of biotite in Y3DD are all consistent with post-kinematic metasomatism (Dharmapriya et al. 2020; Bidgood et al. 2023) and may have resulted from an influx of hydrous fluids during a subsequent tectonothermal event.

5.a. Metamorphic conditions

Figure 8a presents the isochemical phase diagram (pseudosection) calculated for the thin section of specimen Y5. The observed mineral assemblage of garnet–biotite–plagioclase–quartz–ilmenite (+melt) is stable over a broad range of upper-amphibolite facies conditions. However, the overlap of pyrope, grossular and spessartine isopleths refines our estimate of the peak-temperature equilibrium conditions to ∼750 ± 30 °C and 0.66 ± 0.05 GPa. Sillimanite, although observed as intergrowths with biotite (Fig. 8c), is notably absent at the interpreted peak-temperature conditions. Although it is possible that prograde (sub-solidus) sillimanite persisted metastably, as indicated by the dashed prograde arrow in Fig. 8a, we favour the interpretation that the sillimanite present in the specimen was a retrograde phase, stabilized upon exhumation and attendant ductile deformation as the rock’s metamorphic conditions evolved through the garnet–biotite–sillimanite–cordierite–plagioclase–quartz–magnetite stability field (Fig. 8a). Our interpretation is supported, in part, by the overlapping Ti-in-biotite temperatures of 685 ± 24 °C and an avPT results of 709 ± 82 °C, 0.51 ± 0.08 GPa when sillimanite and a stoichiometric cordierite (with acrd = 0.43 ± 0.022 and afcrd = 0.15 ± 0.017) are considered in equilibrium. A portion of the retrograde history of Y5 is recorded by the coexistence of garnet, cordierite and sillimanite, consistent with decompression to ∼0.5 GPa ± 0.1 GPa at ∼675 ± 30 °C. As shown in Fig. 8b, c, modal abundances of both garnet and plagioclase are predicted to increase during prograde metamorphism and decrease along the predicted retrograde path.

Figure 8. Phase diagrams for specimens Y5 (a–c) and Y3DD (d–f), all calculated across a P-T range of 0.4–1.0 GPa and 500–800 °C. (a) Stable phase assemblage diagram overlain by Ti-in-biotite temperature and avPT error ellipses (1 S.D.), yellow polygon denotes stability of garnet core composition (±10% of end-member mole fraction). Dashed arrows mark our interpreted P-T path. (b) Garnet isomodes. (c) Plagioclase isomodes. (d) Stable phase assemblage diagram overlain by Ti-in-biotite temperature and avPT error ellipses (1 S.D.). Dashed arrow marks our interpreted P-T path. (e) Garnet isomodes. (f) Plagioclase isomodes.

As illustrated in the pseudosection shown in Fig. 8d, specimen Y3DD records a similar metamorphic history to Y5, with peak temperatures of ∼790 ± 30 °C based on the coexistence of garnet–K-Feldspar–sillimanite–biotite–plagioclase–quartz. Given the low modal abundance of garnet and its highly fractured and corroded appearance, the pressures of peak-temperature metamorphism were not determined by the position of overlapping isopleths. Instead, we favour the avPT result of 751 ± 42 °C, 0.65 ± 0.13 GPa for our estimate of peak-temperature conditions. The Ti-in-biotite temperature of 683 ± 24 °C likely records the final condition of matrix biotite equilibration as the rock cooled towards the solidus (Dyck et al. 2021). The post-kinematic growth of muscovite could have occurred over a broad temperature range from <500 to 700 °C and need not be related to the peak-metamorphic event recorded in these samples (Fig. 8d). In contrast to Y5, the modal abundance of plagioclase in Y3DD would have decreased during prograde metamorphism once the temperature was in excess of ∼700 °C (Fig. 8f).

5.b. Monazite geochronology

Monazite from the metamorphic specimens analysed record two main growth events; one that peaked between ca. 105 and 95 Ma and a second that peaked ca. 75 Ma (Fig. 5e, f). The older event in both specimen Y5 and specimen Y3DD is characterized by higher Y concentrations, lower Gd/Yb (Fig. 5c, d) and sector zoning in monazite (Figs. S1 and S2). Sector zoning in monazite is typically associated with crystals that have grown from a melt (Stepanov et al. 2012; Catlos, 2013), and as such, the monazite cores are interpreted to reflect the timing of emplacement of the protolith plutonic rocks. The younger population of monazite is spatially associated with patchy or unzoned grain mantles and rims (Figs. S1 and S2). Gd/Yb and Y concentrations generally increase and decrease, respectively, relative to the older population. Assuming that garnet is the main phase controlling HREE availability in these rocks (Foster et al. 2000; Pyle et al. 2001; Foster et al. 2004; e.g. Buick et al. 2006; Larson et al. 2019), decreased availability of Y and HREE for monazite in the younger population would indicate garnet growth/stability at that time. Such an interpretation is consistent with the results of phase equilibria modelling which outline a peak, garnet-stable assemblage, with minimal melt or free water both of which can significantly impact the monazite record (Kohn et al. 2005; Shrestha et al. 2019; Larson et al. 2022). Moreover, the behaviour of Eu, decreasing and increasing towards the younger populations in Y3DD and Y5, respectively (Fig. 5c, d), is also consistent with phase equilibria modelling prediction of plagioclase stability (see Philpotts & Schnetzler, 1968) which increases along the prograde path in Y3DD and decreases in Y5. The ca. 75 Ma monazite population is therefore interpreted to reflect peak metamorphism in the two specimens analysed.

5.c. Rb–Sr and apatite U–Pb geochronology

The wedge-shaped spread in in situ Rb–Sr analyses from specimen Y5 outlines is interpreted to reflect two distinct events. The older, 62.5 ± 0.8 Ma isochron likely reflects post-peak metamorphism cooling through the closure of Sr mobility or Early Paleocene thermal perturbation, perhaps related to metamorphism in the Hunza portion of the Karakoram region at this time (Fraser et al. 2001; Foster et al. 2004). The younger, 41.1 ± 1.4 Ma isochron overlaps with the 38.3 ± 4.5 Ma apatite date from the same rock. A minor, late fluid and/or thermal modification of Y5 may be evidenced by the local chlorite replacement of mica and late, cross-cutting chlorite veins (Fig. 3b). Both apatite U–Pb (Ribeiro et al. 2020) and mica Rb–Sr (Larson et al. 2023) are sensitive to fluid modification and as such it is interpreted that the two geochronological systems may record the same event.

5.d. Plutonic rocks

The zircon from plutonic rocks examined in this study has chemistries that plot across the continental arc-type field in U/Yb vs. Nb/Yb space (Fig. 9; Grimes et al. 2015). An arc affinity is consistent with Y7 (94.9 ± 0.3 Ma) having been sampled from the eastern extent of the Ghamu Bar pluton (Fig. 1b), previously interpreted as part of the broader Karakoram axial batholith generated during mid-Cretaceous subduction beneath the southern margin of Asia (Debon et al. 1987). The Ghamu Bar pluton is part of a linear series of composite continental arc plutons that extend to the southwest away from the study area. To the south-east, these include the ca. 104 Ma Phargam granite (Heuberger et al. 2007) and ca. 110–102 Ma Buni Zom pluton (Faisal et al. 2016), while to the west and north, rocks in the Hunza pluton been dated via whole rock Rb–Sr as 97 ± 17 (Debon et al. 1987), and zircon U–Pb as 95 +4/-6 Ma (lower intercept; Lefort et al. 1983) and 105.7 ± 0.5 Ma (concordant weighted mean; Fraser et al. 2001).

Figure 9. Zircon trace element ratio plot showing Nb/Yb versus U/Yb for each spot analysis. The plot and expected field for continental arc ratios are after Grimes et al. (2015).

Specimens Y9 and Y10 were collected from south of the suture (Fig. 1b). Previous whole rock Rb–Sr geochronology of specimens from the Gindai pluton (Y9) yielded a date of 59 ± 2 Ma (Debon et al. 1987), which overlaps with the zircon U–Pb date of 60.9 ± 0.2 Ma reported herein. The weakly deformed Paleocene plutonic rocks in the Kohistan terrane have been interpreted to reflect subduction under the accreted arc prior to continental collision (Jan et al. 1981; Debon et al. 1987), consistent with the zircon chemistry reported herein.

6. Discussion

Initial mapping of the Hindu Raj area suggested that the plutonic rocks therein may be contiguous with the Karakoram batholith to the east in the Hunza region (Jan et al. 1981; Ebblin, 1988). The Ghamu Bar pluton is not directly connected to the main Hindu Raj – Hunza Karakoram batholith but consists of similar pre-collision granodiorites (Fig. 1b; Gaetani et al. 1996). The new geochronological data presented herein, from both meta-plutonic and plutonic rocks overlap in age with those reported from farther east (Searle et al. 1990 and references therein; Crawford & Searle, 1992) including the interpreted ca. 106 Ma crystallization of the Hunza plutonic unit (Fraser et al. 2001) in the adjacent Hunza region (Fig. 1a). Moreover, zircon trace element signatures are consistent with the plutonic rocks developing as part of a continental arc, again matching that expected for the Karakoram batholith.

The Hindu Raj area records Late Cretaceous age amphibolite facies metamorphism. The timing of peak metamorphism (ca. 80–75 Ma) overlaps with records from across the Hindu Kush to the west (Faisal et al. 2014; Soret et al. 2019) and sparce evidence from the Hunza portion of the Karakoram to the east (Fraser et al. 2001; Foster et al. 2004). Our peak P-T condition estimates of ∼750 ± 30 °C and 0.65 ± 0.05 GPa are higher temperature than those calculated for the same time period in the Hindu Kush (500–600°C and 0.6 GPa; Soret et al. 2019) and the Hunza Karakoram (553 ± 25 °C and <0.6 GPa; Foster et al. 2004).

Previous work in the broader region is divided regarding the potential significance of metamorphism and deformation in Late Cretaceous time. Some researchers suggest that the Late Cretaceous rock record in the south Asian margin reflects the collision of the Kohistan island arc (Petterson & Windley, 1985; Petterson et al. 1991; Weinberg et al. 2000; Fraser et al. 2001; Rehman et al. 2011; Faisal et al. 2014), while others have suggested that docking of Kohistan occurred as recently as 50 Ma (Khan et al. 2009) and that Late Cretaceous deformation and metamorphism may actually be linked to widespread, long-lasting plutonism (Burg, 2011). The plutonic rock dated (Y7) from north of the Shyok suture in the Hindu Raj (Fig. 1b), returned an age of 94.9 ± 0.3 Ma, however, which is significantly older than the timing of peak metamorphism ca. 75 Ma. The lack of widespread and long-lasting Late Cretaceous plutonism in both the current study area and the adjacent Hindu Kush region farther west, which both record metamorphic monazite growth at that time (see Faisal et al. 2014; Soret et al. 2019), is inconsistent with a non-collisional origin. Instead, we conclude that 80–70 Ma monazite growth in the Hindu Raj, and coeval records in adjacent regions, reflects collision of the Kohistan Arc and closure of the Shyok suture.

The record of Late Cretaceous high-temperature metamorphism in the study area appears to have ended by ca. 70 Ma, before the older of the Rb–Sr isochrons that bracket the in situ analyses in Y5 at 62.54 ± 0.79 Ma. That date overlaps an ⁴⁰Ar/³⁹Ar biotite cooling date of 61.6 ± 1.1 Ma (Faisal et al. 2018) from the Buni Zom pluton (110–104 Ma – zircon U-Pb; Faisal et al. 2016) in the Hindu Kush region farther west. The Rb–Sr date also overlaps reported metamorphism in the Hunza portion of the Karakoram region to the east (ca. 63 Ma, Fraser et al. 2001; Foster et al. 2004) indicating a decoupling from the thermal history of the Hindu Raj at that time. The young Rb–Sr isochron that is interpreted to record a fluid or thermal event at 41.1 ± 1.4 Ma is coincident with evidence of widespread continental collision from the southern side of the Kohistan arc at ca. 40 Ma (see summary in Soret et al. 2021).

7. Conclusions

New geochronological and metamorphic data from the Hindu Raj region demonstrate affinity to rocks in the adjacent Hindu Kush and Karakoram regions. The metamorphic record in the Hindu Raj is dominated by the north-directed subduction beneath the Asian continental margin and Late Cretaceous/Paleocene collision of the Kohistan arc with the south Asian margin. A subsequent Paleocene thermal event in the Hindu Raj, however, is coeval with metamorphism in the Karakoram that is not recorded in the Hindu Kush. The geological similarities the Hindu Raj shares with the Karakoram and Hindu Kush areas may not be surprising given its spatial position. The new information presented herein, though, is foundational to linking the geological histories across the south Asian margin prior to continental collision and forming a complete picture of what the antecedent structure of the margin was prior to the arrival of the Indian plate.