The chemical compositions of magmatic zircon growth zones provide powerful insight into evolving magma compositions due to their ability to record both time and the local chemical environment. In situ U-Pb and Hf isotope analyses of zircon rims from Oligocene–Miocene leucogranites of the Bhutan Himalaya reveal, for the first time, an evolution in melt composition between 32 and 12 Ma. The data indicate a uniform melt source from 32 Ma to 17 Ma, and the progressive addition of an older source component to the melt from at least ca. 17 Ma. Age-corrected ɛHf ratios decrease from between −10 and −15 down to values as low as −23 by 12 Ma. Complementary whole-rock Nd isotope data corroborate the Hf data, with a progressive decrease in ɛNd(t) from ca. 18 to 12 Ma. Published zircon and whole-rock Nd data from different lithotectonic units in the Himalaya suggest a chemical distinction between the younger Greater Himalayan Series (GHS) and the older Lesser Himalayan Series (LHS). The time-dependent isotopic evolution shown in the leucogranites demonstrates a progressive increase in melt contribution from older lithologies, suggestive of increasing LHS involvement in Himalayan melting over time. The time-resolved data are consistent with LHS material being progressively accreted to the base of the GHS from ca. 17 Ma, facilitated by deformation along the Main Central thrust. From 17 Ma, decompression, which had triggered anatexis in the GHS since the Paleogene, enabled melting in older sources from the accreted LHS, now forming the lowermost hanging wall of the thrust.

Crustal melting is a fundamental process both for chemical differentiation and for facilitating ductile deformation of the Earth’s continental crust. In the Himalaya, late Oligocene–Miocene leucogranites provide a well-documented example representing crustal melting induced by continental collision. Oxygen and hafnium isotope data from Himalayan leucogranite zircon rims demonstrate that these granites formed via pure crustal melting with no detectable mantle input (Hopkinson et al., 2017). Melting took place along the Himalaya orogen from at least 25 Ma to 9 Ma (Guo and Wilson, 2012), although in southern Tibet, crustal melting of Eocene age has been reported (Aikman et al., 2012). Despite this extensive period of regional anatexis, no time-dependent change in the magmatic source has been recognized thus far. The evolving thermal regime and/or the continued burial of different material into a stable thermal regime during collisional orogenesis should lead to migration of melt source regions or change in source rocks through time, and to the involvement of sources from different crustal levels (Fiannacca et al., 2017). The recognition of the temporal evolution of melt zones would contribute to more refined models of thermal and structural evolution of the continental crust during orogenesis. The apparent absence of time-dependent trends in the Himalaya may result from the lack of chronological resolution. Here we evaluate the Hf isotope composition of zircon sampled from leucogranites in Bhutan (Eastern Himalaya), along with new whole-rock Nd isotope data. We observe a temporal evolution in the source region, and these data provide, for the first time, direct evidence for a time-dependent change in the mid-crustal material undergoing melting and decompression during Himalayan crustal thickening.

Across the orogen, leucogranites intrude the Greater Himalayan Series (GHS), an upper amphibolite-facies metasedimentary package of primarily Neoproterozoic protolith age (Ahmad et al., 2000; Gehrels et al., 2011). We collected 12 samples from nine localities in Bhutan, where the leucogranites form discrete sheets within the uppermost GHS, close to the tectonic contact with Tethyan sediments. Samples include two-mica (samples 1G03, 3A03, 4D01, 1247, 1251, CWB16, and CWB23), tourmaline-bearing (samples 1G02 and 1215), garnet-bearing (samples 1D01 and 3A02), and pegmatitic two-mica (sample 1G01) types. Samples 1G01, 1G02, and 1G03 are from intersecting sheets from the same exposure, as are samples 3A02 and 3A03 (Fig. 1). We have previously analyzed the isotopic composition of zircon growth zones from all of these samples (Hopkinson et al., 2017).

The base of the GHS is marked by the Main Central thrust, an orogen-parallel thrust zone that has facilitated India-Asia convergence since at least the late Paleogene (Mottram et al., 2014, 2015). Structurally below the Main Central thrust is the Lesser Himalayan Series (LHS), a primarily Paleoproterozoic-sourced stack of metasedimentary rocks (Ahmad et al., 2000; Gehrels et al., 2011). Whereas metamorphism in the GHS commonly reached temperatures >750°C during the Miocene (Harris et al., 2004), peak metamorphic temperatures in the LHS rarely exceeded ∼600°C. Kyanite-grade LHS schists that reached temperatures of ∼640°C have, however, been reported from the upper levels of the Main Central thrust zone in the western Himalaya (Caddick et al., 2007).

Both the GHS and LHS comprise a mix of pelitic, granitic, carbonate, and quartzite compositions. The pelitic assemblages are the most melt fertile, and therefore provide appropriate source materials for anatectic melts at temperatures <∼760°C (at pressures estimated for mid- to lower-orogenic crustal depths; Patiño Douce and Johnston, 1991).

Whole-rock isotope geochemical data from equivalent High Himalayan leucogranites exposed across the Himalayan orogen suggest that they formed largely by partial melting of the pelitic lithologies of the GHS into which the granites intrude (Hopkinson et al., 2017). Small contributions from melting LHS lithologies to some leucogranites have been suggested based on whole-rock data (Aikman et al., 2012; Guo and Wilson, 2012; Hu et al., 2018), but overall, melt production in the Himalaya has generally been ascribed to decompression of the GHS during thrusting along the underlying Main Central thrust. Critically, no temporal control on possible inputs from the LHS has yet been observed.

Uranium-lead (U-Pb) and Hf isotopic analysis were performed at the Natural Environment Research Council (NERC) Isotope Geosciences Laboratory (NIGL; Nottingham, UK) using a Nu Instruments Nu Plasma multicollector–inductively coupled plasma–mass spectrometer (MC-ICP-MS) coupled to a New Wave Research (NWR) 193ss Nd:YAG laser, and a Thermo Neptune Plus MC-ICP-MS coupled to a NWR 193UC excimer laser, respectively; full methods and data are reported in Hopkinson et al. (2017). Common Pb corrections have been reevaluated for the U-Pb data set (Table DR1 in the GSA Data Repository1) to provide accurate age constraints for Figure 2. Whole-rock Sm-Nd isotope analyses (Table DR2) were obtained by thermal ionization mass spectrometry and trace element analyses by laser ablation–ICP–MS (LA-ICP-MS) (Table DR3) at the Open University (Milton Keynes, UK); full methods and data are provided in the Data Repository.

Zircon from all 12 samples yielded rims of Cenozoic age ranging from 34 to 11 Ma (Fig. 2; Table DR1), and the majority have inherited cores of early Paleozoic or Proterozoic age (Table DR2). U-Pb Tera-Wasserburg and concordia diagrams and cathodoluminescence images are provided in the Data Repository (Figs. DR1–DR12).

Zircon rims from multiple grains in individual samples record protracted crystallization that varies between 0.9 and 15.5 m.y. in duration in different samples (Fig. 2; Fig. DR13). This age dispersion in single hand samples is a typical observation in Himalayan granites; previous studies have suggested protracted melt generation, segregation, amalgamation, and crystallization in each individual granite body as causes (e.g., Lederer et al., 2013). Furthermore, our data suggest that most leucogranite bodies are derived from heterogeneous batches of melt accumulated over significant periods of time.

Age-corrected Hf isotope ratios (ɛHf(t)) range from −8.4 to −29.8 across the whole data set. A plot of zircon rim Hf isotopic compositions against corresponding age (Fig. 2) shows that some samples preserve large differences in Hf isotopic composition (e.g., samples 1G03 and 1247). This is likely a result of the lack of homogenization in the evolving magma due to high melt viscosities of low-temperature, H2O-undersaturated siliceous magmas (Harris et al., 2000), which has previously been documented in Himalayan leucogranite whole-rock data (Deniel et al., 1987). Variability in zircon rim Hf isotope compositions at the hand-sample scale can be linked to heterogeneous magma sources (Villaros et al., 2012) or to disequilibrium melting with different mineral phases hosting different reservoirs of radiogenic Hf (Farina et al., 2014; Tang et al., 2014).

Our data show a distinct secular change in εHf(t) values, with increasingly lower values in younger samples occurring after 17 Ma. After this time, the population of all but one sample (1215) decreases to ɛHf(t) <−16. Zircon rims older than 17 Ma have depleted-mantle model ages of 1450–2000 Ma, whereas those of zircon rims younger than 17 Ma range from 2000 to as old as 2620 Ma (Table DR1).

Whole-rock εNd(t) ratios range from −7.4 to −20.9 (Table DR2; Fig. 2, inset). There is a general decreasing trend through time across 11 of the 12 samples. This trend broadly correlates with the in situ zircon εHf data.

An Evolving Melt Source

We attribute the secular change in melt composition indicated in these samples to a changing melt source through time. Significant disequilibrium effects caused by minerals such as garnet in the source would typically lead to uncorrelated Lu-Hf and Sm-Nd isotope compositions of the melt, which we do not observe. We discount disequilibrium melting as a dominant control on the compositions because the trend displayed by the Hf isotopes is directly correlative to that of the Sm-Nd isotope system.

Previous studies implicate the GHS as the dominant source for Himalayan leucogranite melts (Harris and Massey, 1994; Guo and Wilson, 2012; Hopkinson et al., 2017). To address the possible cause of the evolving isotopic signature of the Bhutanese leucogranites in terms of melt source, we compare our new zircon-derived data with existing whole-rock Nd and Hf isotope data from the Himalayan metasedimentary formations. Previous studies describe a clear distinction between GHS and LHS units in terms of their isotopic signature (Richards et al., 2006). For example, whole-rock Nd data taken from across the Main Central thrust in Sikkim (Mottram et al., 2014) indicate εNd values of −12.1 to −18.3 in the GHS (hanging wall), and −23.4 to −27.7 in the LHS (footwall). These overlap directly with the ranges obtained from metasedimentary formations in Bhutan of εNd −12.1 to −17.6 for the GHS and −25.9 to −32.3 for the LHS (Richards et al., 2006). Whole-rock Nd isotope data for the samples in this study (Fig. 2, inset) have values of εNd(t) ranging from −7.4 to −20.9. Bulk rock data from similar leucogranites east of the study region lie on the same array (Aikman et al., 2012).

Our in situ zircon Hf isotope versus age data are plotted against a compilation of available detrital zircon data from the GHS and LHS (Spencer et al., 2018) in Figure 3. For purposes of comparison, all data are age corrected to 17 Ma, a critical age identified with a changing tectonic regime in this and previous studies (Mottram et al., 2014). The compiled data show that the GHS yields a broad range in εHf(t = 17 Ma) values, from −8 to −30, with a dominant population between −8 to −18. The LHS mainly falls between εHf(17 Ma) −21 and −35. As with the Nd isotope data, the leucogranite zircon Hf data are consistent with predominant derivation from a GHS source, but indicate a likely increasing contribution from an older (i.e., LHS) source through time. The lower values of εHf(17 Ma) require discrete melting of more-evolved GHS units or an increasing melt contribution from LHS source rocks. Importantly, isotopic studies of the GHS indicate a homogeneous zircon age and Hf isotopic signature at all structural levels (Spencer et al., 2012); thus, we favor a model of an increasing LHS component through time.

In addition, the Hopkinson et al. (2017) data set shows that samples with 1800–1900 Ma zircon cores are found only in samples with younger zircon rims (Table DR2; Fig. DR14). The well-documented Paleoproterozoic orogenic event of this age affected the LHS (Ahmad et al., 2000; Kohn et al., 2010) but predates the deposition of the GHS during the Neoproterozoic (Richards et al., 2006). Although it appears that some Paleoproterozoic material was reworked into the GHS (Spencer et al., 2012, 2018), together the isotopic evidence and the ages of zircon cores strongly suggest that an increasing component of older and more radiogenic material melted during the period of anatexis from 17 Ma.

Given that the available melt-fertile (i.e., pelitic) source material in the Himalaya is restricted to the LHS and GHS, the results from this study provide clear evidence for a GHS-sourced melt zone throughout the Oligocene, and subsequent tectonic or thermal evolution allowing some melting of the LHS during the early Miocene (since 17 Ma). To quantify this trend, the average pre–17 Ma value of εHf is −12, which is reduced to −22 in several samples by 12 Ma. If we assume that contributions of LHS are characterized by εHf = −27 (peak LHS ratio; Fig. 3), then the data suggest that by 12 Ma, ∼70% of the contributing melts were sourced from the LHS.

The age trend in the source region depends on five samples providing εHf <−17: 1G01, 1G02, 1G03, 1251, and CBW16. We note that where samples were collected from discrete sheets at the same outcrop (samples 1G01, 1G02, and 1G03), all sheets recorded the same trend. Importantly, none of this group of samples contains zircon older than 17 Ma, while some samples with post–17 Ma zircon include εHf values >−17 (samples 1215, 1251, 4D01, 3A03). This suggests that following 17 Ma, the older source regions were sampled by younger magmatic pulses at discrete localities rather than there being a regional shift in magmatic source regions.

Tectonic Implications

The key finding of this study is that increasing components of older material were incorporated into Eastern Himalayan anatectic melts from 17 Ma. Kinematic models, backed up by petrochronological studies, suggest progressive tectonic accretion of LHS material from the footwall to the hanging wall of the Main Central thrust as a consequence of the structure cutting downward over time (Bollinger et al., 2006). Importantly, this progressive structural accommodation appeared to take place from ca. 17 Ma in Sikkim, west of our study area in Bhutan (Mottram et al., 2015). Metamorphism in the Main Central thrust zone developed as heat advected downward from the overriding GHS material. Prior to accretion, the hanging wall of the GHS was entirely comprised of GHS lithologies such that decompression melting led to GHS-sourced magma (Fig. 4A). Following such accretion, exhumation along the hanging wall of the thrust could have led to decompression melting of both GHS and accreted LHS material (Fig. 4B). The introduction of LHS-sourced material into the melt zone may therefore have marked the initiation of accretion of the LHS into the hanging wall of the thrust.

In situ U-Pb and εHf LA-ICP-MS compositions of zircon rims from Himalayan leucogranites provide the first evidence for secular compositional change of melts generated by protracted anatexis during the earliest Oligocene to mid-Miocene. Increasing components from older source regions from 17 Ma to 12 Ma suggest that melts derived from LHS material were incorporated into the aggregating melt bodies during this period. Rapid exhumation of the hanging wall of the Main Central thrust is well documented through the Miocene. The data support a tectonic model in which older LHS material was accreted to the hanging wall of the Main Central thrust, and thus the base of the GHS, during the progressive evolution of the structure, so allowing decompression melting of the LHS from the time of significant accretion.

Hopkinson acknowledges funding from NERC CASE studentship NE/K501074/1. Warren acknowledges funding from NERC Advanced Fellowship NE/HO16279/1. Analytical work was funded by NERC facility grants IMF478/0513 and IP1403/111. Thanks are extended to D. Regis and D. Young for field assistance, A. Wood for assistance with sample preparation at the NERC Isotope Geosciences Laboratory (NIGL), and T. Argles, B. Charlier, D. Regis, and C. Mottram for fruitful discussions. We thank reviewers Stéphane Guillot, Jörg Hermann, and Patrizia Fiannacca for their constructive comments on previous versions of this paper.

1GSA Data Repository item 2020030, detailed methodology and supplemental figures including Terra-Wasserburg and Concordia diagrams (Figures DR1–DR12), ranked age plot (Figure DR13), and Kernel Density Estimate diagram (Figure DR14), and data tables for U-Pb isotopes (Table DR1), Sm-Nd isotopes (Table DR2), and bulk rock elemental analyses (Table DR3), is available online at, or on request from
Gold Open Access: This paper is published under the terms of the CC-BY license.