Rare granulitized eclogites exposed in the eastern Himalaya provide insight into conditions and processes deep within the orogen. Sensitive high-resolution ion microprobe (SHRIMP) U-Pb, Ti, and rare earth element (REE) data from zircons in mafic granulitized eclogites located in the upper structural levels of the Greater Himalayan Sequence in Bhutan show that zircon was crystallized under eclogite-facies metamorphic conditions between 15.3 ± 0.3 and 14.4 ± 0.3 Ma, within a couple million years of the later granulite-facies overprint. In conjunction with pressure estimates of the eclogite- and granulite-facies stages of metamorphism, the age data suggest that initial exhumation occurred at plate-tectonic rates (cm yr−1). These extremely rapid synconvergence exhumation rates during the later stages of the India-Asia collision require a revision of theories for the transportation and exhumation of crustal materials during continental collisions. In contrast to western Himalayan examples, the eastern Himalayan eclogites cannot be tectonically related to steep subduction of India beneath Asia. Instead, they more likely represent fragments from the base of the overthickened Tibetan crust. Based on the zircon age and trace-element data, we hypothesize that the protolith of the mafic granulites was middle Miocene mafic intrusions into the lower crust of southern Tibet, linked to Miocene volcanism in the Lhasa block. We suggest that a transient tectonic event—possibly the indenting of a strong Indian crustal ramp into crust under southern Tibet that had been weakened by partial melting—may have promoted exhumation of the eclogitized lower crust under Tibet. The mafic magmatism and volcanism themselves may have been related to the convective thinning of the lithospheric mantle triggered by a reduction in the India-Eurasia convergence rate during the middle Miocene, which in turn could have facilitated the rapid extrusion of the lower crust over the earlier-exhumed middle crust.
In the Himalaya, two contrasting types of high-pressure (HP) metamorphic units have been described (Lombardo and Rolfo, 2000). Close to the suture between the Indian and Asian plates in the western Himalaya (Tso Morari in Ladakh and Kaghan in NW Pakistan), ultrahigh-pressure (UHP) continental rocks were metamorphosed up to >3.9 GPa and >750 °C just prior to, or during, the initial stages of India-Asia collision (ca. 55 Ma) and underwent minor metamorphic overprint during exhumation (Guillot et al., 2008, and references therein). In contrast, evidence for HP metamorphism within the metamorphic core in the central and eastern Himalaya (Fig. 1) is limited to rare mafic rocks hosted by strongly migmatized rocks of Indian affinity hundreds of kilometers from the suture (Chakungal et al., 2010; Guillot et al., 2008; Lombardo and Rolfo, 2000). These mafic rocks show textures compatible with initial eclogite-facies equilibration, followed by pervasive and intensive overprinting to lower-pressure granulite-facies mineralogy during the middle Miocene (Cottle et al., 2009; Groppo et al., 2007; Kali et al., 2010; Lombardo and Rolfo, 2000; Rolfo et al., 2008; Warren et al., 2011b). Their metamorphic evolution and the available geochronological data suggest synconvergent exhumation from lower-crustal levels, but the timing of the precursor eclogite-facies metamorphism is presently unconstrained.
Previously reported numerical geodynamic models of Himalayan tectonics (e.g., Jamieson et al., 2004) predict the exhumation of midcrustal material from under Tibet, driven by a combination of a topographic pressure gradient and focused erosion. These models do not currently predict the exhumation of lower orogenic (>50 km, i.e., >1.4 GPa) crustal material. The geochronological and thermobarometric evolution of the eastern Himalayan HP rocks may therefore provide the key to understanding lower-crustal evolution in continental collision zones in general and the Himalayan orogen in particular.
We analyzed zircons from mafic granulitized eclogites and their host rocks in NW Bhutan (for sample locations, see GSA Data Repository Fig. A11) for U-Pb ages, Ti-in-zircon temperatures, and rare earth element (REE) chemistry to provide constraints on the origin, evolution, and exhumation of eastern Himalayan HP metamorphic rocks. We then combined observations and data from the surface geology of southern Tibet, geophysical data, plate configuration models, geodynamic models, and the geology of the Himalayan metamorphic core into an internally consistent model of the tectonics of the Himalayan orogen.
The studied rocks are exposed within the uppermost Greater Himalayan Sequence in NW Bhutan (Fig. 2). The Greater Himalayan Sequence is bound at its base by the Main Central thrust and at its roof by the South Tibetan detachment system, both of which are north dipping, but which have an opposite shear sense. These crustal-scale coeval shear zones operated between ca. 24 and ca. 12 Ma (Chambers et al., 2011; Daniel et al., 2003; Kellett et al., 2009, 2010). The South Tibetan detachment system is a system of one to several normal-sense brittle faults and/or ductile shear zones. In Bhutan, the South Tibetan detachment system consists of two main structures. The structure located closer to the orogenic front, the outer South Tibetan detachment system (Kellett et al., 2009), is ductile and forms the base of synformal klippen while also separating migmatites and gneisses of the Greater Himalayan Sequence below from the metasedimentary Chekha Group above (Grujic et al., 2002; Kellett et al., 2009). The continuously exposed segment of the South Tibetan detachment system to the north is the inner South Tibetan detachment system (Kellett et al., 2009). Ductile shear on the outer South Tibetan detachment system occurred during the Miocene, from ca. 22 Ma until at least ca. 16 Ma (Chambers et al., 2011; Kellett et al., 2009, 2010), while ductile shear on the inner South Tibetan detachment system ceased by ca. 11 Ma (Edwards et al., 1996, 1999; Kellett et al., 2009; Wu et al., 1998). An out-of-sequence thrust within the Greater Himalayan Sequence, the Kakhtang thrust (Davidson et al., 1997; Gansser, 1983; Grujic et al., 2002), doubles the structural thickness of the Greater Himalayan Sequence, and klippen of Tethyan rocks soled by the South Tibetan detachment system are preserved south of this structure. In northwestern Bhutan, the Kakhtang thrust, which was active between ca. 14 and ca. 10 Ma (Grujic et al., 2002), is probably continuous with the Laya thrust (Fig. 2), which constitutes a sharp contact between the granulite-bearing unit above and the amphibolite-bearing unit below (Davidson et al., 1997; Swapp and Hollister, 1991; Warren et al., 2011b).
The variably retrogressed mafic granulites described in this study and by others (Chakungal et al., 2010; Warren et al., 2011b) are exposed as centimeter- to meter-scale boudins, hosted by migmatized metasediments and felsic orthogneisses within the core of a regional antiform located ∼2–3 km beneath the South Tibetan detachment system, and in the hanging wall of the Laya thrust (Fig. 2). The dikes are oblique to the observed lithological banding in the country rocks and are deformed by all the structures affecting the country rocks. This suggests a primary intrusive association and not a tectonic incorporation of mafic rocks into the felsic rocks. The dominant (gneissic) foliation is affected by a conjugate set of shear bands; north-dipping, top-to-the-NW shear bands are more dominant to the north, whereas south-dipping, top-to-the-SE shear bands are dominant to the south (Fig. 2, inset). The pervasive stretching lineation, often accompanied by sillimanite slicken fibers along foliation planes and shear bands, trends NW-SE (Fig. 2, inset). In the uppermost kilometer of the Greater Himalayan Sequence, a second foliation, parallel to the set of shear bands (Fig. 2, inset), increasingly transposes the dominant foliation and eventually forms the mylonitic foliation of the South Tibetan detachment system. The general regional structure is therefore a recumbent dome with normal fault geometry shearing along the northern flank and a flat thrust along its southern base. Toward the SW, the E-W–trending antiform swings into a NNE-SSW trend, where it is flanked by two opposite-dipping normal faults, the Yadong and Lingshi faults, which shape the Jomolhari Massif into a horst (Fig. 2).
The mafic rocks are medium to coarse grained (porphyroblast size ∼1–20 mm) and contain Grt + Cpx + Pl + Opx + Amp + Qz ± Bt ± Kfs with accessory Ap, Zrn, Rt, Ilm, Ttn, and Mnz (Chakungal et al., 2010; Warren et al., 2011b; mineral abbreviations after Whitney and Evans, 2010). They preserve textural evidence for earlier eclogite-facies metamorphism, including rutile inclusions in garnet, clinopyroxene + plagioclase symplectites, suggesting precursor omphacite, and textural evidence for a plagioclase-absent assemblage at the time of peak pressure metamorphism (Chakungal et al., 2010; Warren et al., 2011b). They are similar in texture to the granulitized eclogites from the Ama Drime Massif in southern Tibet, and from Sikkim in northern India (Cottle et al., 2009; Groppo et al., 2007; Lombardo and Rolfo, 2000; Rolfo et al., 2008). Plagioclase is only found in association with decompression textures: coronas of symplectitic An + Opx ± Amp replaced Grt, and matrix Cpx was replaced by Pl ± Opx ± Amp (Fig. 3) (Chakungal et al., 2010; Warren et al., 2011b). REE concentrations in relict garnet show no negative Eu anomaly, again suggesting a plagioclase-absent assemblage at the time of its growth (Warren et al., 2011b). The pressure-temperature (P-T) path of mafic granulites from the Ama Drime Massif has been suggested to include decompression and heating from HP conditions of >1.5 GPa and >580 °C to high-temperature (HT) conditions of 0.7–1.0 GPa and ∼750 °C, followed by decompression to 0.3 GPa at ∼630 °C (Cottle et al., 2009; Groppo et al., 2007; Lombardo and Rolfo, 2000). Similar metamorphic conditions and a similar path are inferred for the Bhutan granulitized eclogites (Chakungal, 2006; Warren et al., 2011a, 2011b) based on similar mineral compositions and textures (Fig. 3).
Associated metapelites preserve evidence for granulite-facies metamorphism, but not for eclogite-facies metamorphism. The granulite-facies stage is represented by Grt + Opx + Pl + Rt + Qz, which suggest conditions of ∼800 °C and >0.8 GPa (Warren et al., 2011b). Cordierite and spinel may have been introduced during decompression during breakdown of garnet + sillimanite at ∼0.3 GPa and >750 °C (Warren et al., 2011b). The banded migmatitic orthogneiss in the area does not preserve mineralogical evidence for granulite- or eclogite-facies metamorphism and yields peak conditions of ∼700 °C at ∼0.45 GPa (Warren et al., 2011b).
Miocene leucogranite bodies are interspersed throughout the study area. The youngest bodies (∼15.5–11.0 Ma) are Ms + Tur + Crd ± And leucogranites, suggesting crystallization at depths corresponding to ≤0.28 GPa (Kellett et al., 2009). In the uppermost intrusions, a clear strain gradient can be traced from undeformed granite at the base to S/C′ mylonite at the top, caused by ductile top-to-the-NW shear.
Dating of peak pressure conditions in the eastern Himalaya is complicated by a lack of suitable mineral assemblages, and chemical and/or textural overprinting during exhumation. Determination of P-T conditions of the granulite-facies overprint is also complicated by chemical reequilibration at high temperatures (Frost and Chacko, 1989). Furthermore, the elements of use in geobarometry have higher closure temperatures for diffusion than the elements for geothermometry, and one cannot always be certain where on the P-T path the conditions have been “locked in,” at maximum P and T or during cooling (Frost and Chacko, 1989). Pseudosection analysis is similarly limited by uncertainty in the effective bulk composition during the P-T conditions of interest and the equilibrium mineral compositions at that time. Therefore, linking the growth of the dated accessory minerals to specific stages in the growth of the major mineral assemblage is nontrivial, and we acknowledge the uncertainties in our approach and pressure-temperature-time (P-T-t) interpretations.
Using circumstantial evidence and the assumption that the eclogites were related to subduction corresponding to the initial stages of continental collision, Lombardo and Rolfo (2000) concluded that a late Eocene–early Oligocene age (44–33 Ma) was the most likely timing of the eclogitic stage in the Ama Drime Massif of southern Tibet. Sensitive high-resolution ion microprobe (SHRIMP) U-Pb analyses yielded zircon core ages ranging from 110 to 88 Ma, interpreted as protolith ages (Rolfo et al., 2005). A zircon SHRIMP 206Pb/238U age of 29.5 ± 0.4 Ma from a single spot was interpreted as timing of the final closure of the Neotethys (Li et al., 2003). Zircon 206Pb/238U ages of 33 ± 2 Ma from a garnet sillimanite gneiss near Kharta (Liu et al., 2007) were explained as dating the HP stage at 750–800 °C and ∼1.4 GPa. These ages were later reinterpreted as defining the onset of partial melting, most probably during the isobaric temperature increase between 800 °C and 900 °C in the pressure range from 1.4 to 1.8 GPa (Kali et al., 2010).
In the granulitized eclogites of the Ama Drime Massif (Fig. 1), thin zircon rims with low Th/U ratios (0.02–0.03), yielding 206Pb/238U SHRIMP ages of 17.6 ± 0.3 Ma (Li et al., 2003) and 15–12 Ma (Rolfo et al., 2005), were interpreted as dating either the granulite or amphibolite stage, rather than the timing of eclogite-facies metamorphism (Corrie et al., 2009; Groppo et al., 2007). Similarly, ca. 14–12 Ma monazite and xenotime ages from surrounding migmatitic gneisses were interpreted as constraining either the timing of peak granulite metamorphism (Cottle et al., 2009) or of the final crystallization at the onset of cooling after strong adiabatic decompression (Groppo et al., 2007; Kali et al., 2010). Most recently, granulitized eclogites exposed to the south of the Ama Drime range, in the immediate hanging wall of the Main Central thrust (Fig. 1), yielded Lu-Hf garnet dates of 20.7 ± 0.4 Ma (Corrie et al., 2009). This age was interpreted as representing the timing of eclogite-facies conditions.
In NW Bhutan, in the hanging wall of the Laya thrust, zircons from two samples of granulitized eclogite produced U-Pb zircon SIMS (secondary ion mass spectrometry) upper-intercept ages of 1794 ± 58 Ma and 1742 ± 39 Ma, and lower-intercept ages of 256 ± 90 Ma (mean square of weighted deviates [MSWD] = 7.4) and 12.3 ± 4.6 Ma (MSWD = 5.8) (Chakungal et al., 2010). Monazite hosted within the leucosome of neighboring granulite-facies orthopyroxene-bearing metapelite yielded laser ablation–multicollector–inductively coupled plasma–mass spectrometry (LA-MC-ICP-MS) U-Th-Pb ages of 13.9 ± 0.3 Ma, and between 15.4 ± 0.8 Ma and 13.4 ± 0.5 Ma in the host migmatitic orthogneisses (Warren et al., 2011b). Monazite associated with sillimanite-grade metamorphism in gneiss in the footwall of the Laya thrust yielded U-Pb rim ages of 21–17 Ma (Warren et al., 2011). Equally, leucogranites in the granulite-hosting unit yielded zircon SHRIMP ages between 15.5 and 11 Ma (Kellett et al., 2009), while the leucogranites in the lower amphibolite–bearing unit yielded ages ranging between 21 and 16.8 Ma (Carosi et al., 2006) (Fig. 2). These data are consistent with Miocene exhumation of Greater Himalayan Sequence material from a variety of crustal depths at different times; however, they do not constrain the timing of the HP metamorphism, nor the onset of exhumation of the HP rocks.
Linking the absolute timing of deformation and metamorphic events to P-T paths is fundamental to understanding rates of crustal processes. This task is hampered by difficulties in correlating thermobarometric and chronologic data from the same rock sample, and the relatively large uncertainty associated with conventional thermobarometry. The crystallization temperatures of dateable accessory mineral phases can be retrieved using trace-element thermometry: Ti-in-zircon, Zr-in-rutile, and Zr-in-titanite (Ferry and Watson, 2007; Hayden et al., 2008). By linking temperature and U-Pb ages directly, temperature-time points can be immediately established (so long as the closure temperature for diffusion of either species has not been exceeded). Furthermore, the trace-element signature, especially the distinctive REE patterns in accessory minerals, has been demonstrated to represent the metamorphic conditions of their crystallization (Bingen et al., 2004; Corrie and Kohn, 2008; Harley and Kelly, 2007; Rubatto, 2002; Zack et al., 2002). The combination of these two approaches is crucial in order to reliably link U-Pb ages to metamorphic and deformation stages. Our approach for determining the timing of eclogite-facies conditions in Bhutan is to date different accessory minerals from a range of lithologies and connect their growth stage to the corresponding metamorphic stage by collecting trace-element geochemistry at the same analytical spots as those used for age dating (Mazdab and Wooden, 2006; Wooden et al., 2006).
Zircon U-Pb analyses were conducted on the SHRIMP-RG (reverse geometry) ion microprobe co-operated by the U.S. Geological Survey and Stanford University. Zircons were concentrated by standard heavy-mineral separation processes, handpicked for final purity, and then mounted on double-sided sticky tape on glass slides in 1 × 6 mm rows, cast in epoxy, ground, and polished to a 1 μm finish on a 25-mm-diameter by 4-mm-thick disc. All grains were imaged with transmitted and reflected light on a petrographic microscope, and with cathodoluminescence (CL) and backscattered electrons (BSE) as needed on a JEOL 5600 scanning electron microscope (SEM) to identify internal structure, inclusions, and physical defects. The mounted grains were washed with 1 N HCl and distilled water, dried in a vacuum oven, and coated with Au. Mounts were placed into a loading chamber at high pressure (10–7 torr) for several hours before being moved to the source chamber of the SHRIMP-RG. Secondary ions were generated from the target spot with an O2– primary ion beam varying from 4 to 6 nA. The primary ion beam produced spots with a diameter of 20–40 μm and a depth of 1–2 μm for an analysis time of 9–12 min. Nine peaks were measured sequentially on an ETP electron multiplier: 90Zr216O, 204Pb, background 0.050 mass units above 204Pb, 206Pb, 207Pb, 208Pb, 238U, 248Th16O, 254U16O, 40Ca, and 48(Ti,Ca)2. Auto-centering on selected peaks and guide peaks for low or variable abundance peaks (i.e., 96Zr216O 0.165 mass unit below 204Pb) was used to improve the reliability of locating peak centers. The number of scans completed through the mass sequence and counting times on each peak were varied according to sample age and U and Th concentrations to improve counting statistics and age precision. Measurements were made at mass resolutions of 6000–8000 (10% peak height) to eliminate all interfering atomic species.
Concentration data for zircons were standardized against zircon standard R33 (419 Ma, quartz diorite of Braintree complex, Vermont; John Aleinikoff, 2006, personal commun.), which was analyzed approximately every fourth analysis throughout the duration of the session. Errors on standard R33 in individual sessions range from 1% to 2% at 95% confidence. Data reduction followed the methods described by Williams (1997) and Ireland and Williams (2003) and used the Squid and Isoplot programs of Ken Ludwig. Approximately 20 spots from 10 to 15 zircon grains were analyzed for each sample.
Titanium in zircon and REE analyses were performed on the SHRIMP-RG, following the procedure by Barth and Wooden (2010), at the same spots where U-Pb analyses were performed in the previous round. Subsequent to isotopic analysis, the zircon grain mounts were lightly polished to remove the original gold coating and sputtered pits, recoated with gold, and analyzed for a suite of trace elements using an ∼15-μm-diameter, 1–2 nA O2– primary beam. Trace-element analyses were conducted, in so far as possible, in a zircon volume directly beneath or adjacent to that analyzed for isotopic compositions, as well as in additional areas to more fully describe trace-element variations. 254UO, 248ThO, 207Pb, and 206Pb were monitored to allow an estimated Pb/U age (±10%) to be calculated, and, along with reanalysis of Hf concentrations, helped in assuring that trace-element analysis points to be paired with U/Pb ages were sputtered in a compositionally similar zircon volume with respect to luminescence, Hf, U, and Th concentrations, and age. In addition to elements of interest, the elements F, Al, P, Ca, and Fe (and Mg and K in later sessions) were monitored to detect encroachment into the analytic volume of mineral and/or melt inclusions, especially apatite, titanite, and/or oxide minerals, and the possibility of elements being contributed by sources other than zircon, including alteration associated with weathering and metamictization. Detection of phases other than zircon that may contribute to measured Ti concentrations is critical to evaluation of temperature-dependent compositional variations using trace-element thermometry. Trace-element concentrations were standardized against Madagascar Green (MAD) zircons (Mazdab and Wooden, 2006); replicate analyses of fragments of MAD zircons over multiple analytical sessions were used to establish precision of the trace-element analyses. All the trace-element data are listed in Table A1 in the GSA Data Repository (see footnote 1).
Microstructure of Zircons from Mafic Granulites
Three of the dated mafic granulitized eclogite samples have tholeiitic basalt composition (BH219, BH252, BH268), and one has alkaline basalt composition (BH203) (Chakungal et al., 2010). Zircon is present as inclusions in garnet, biotite, and orthopyroxene within the matrix along grain boundaries and adjacent to garnet rims. Up to 30 zircon grains from each of four mafic samples were analyzed. We provide new geochronological data and interpretations from the sample BH219 with respect to those reported by Chakungal et al. (2010).
Zircon grains range in length from 80 to 200 μm along the longest axis (Supplementary Fig. A2 [see footnote 1]). They are subhedral to euhedral in shape, ranging from well-terminated prismatic (3:1) to equant (“soccer-ball”) to rounded/ovoid subhedral grains. The latter two morphologies are typical of zircons that grew under high-grade metamorphic conditions (Hoskin and Black, 2000; Hoskin and Schaltegger, 2003). As evident from scanning electron microscope–cathodoluminescence (SEM-CL) images, less than 10% of grains have xenocrystic cores. These cores are 20–60 μm in length and have variable internal zoning patterns. Some cores have sharp oscillatory to sector zoning, commonly truncated by bright CL and U-poor rims (Fig. 4; Supplementary Fig. A2 [see footnote 1]). Less than 20% of grains contain light, irregular cores with diffuse or blurred zoning. Both faintly zoned rims and fine-scale oscillatory-zoned rims are present (Fig. 4; Supplementary Fig. A2 [see footnote 1]). About 40%–50% of grains show blurred or gently convoluted zoning and overgrowth, which are interpreted as characteristics of metamorphic zircons (Hoskin and Schaltegger, 2003). In addition, faint banded planar and fir-tree sector zoning is present (the latter reflecting strong fluctuations of growth rates). Over 90% of grains are comparatively homogeneous in CL, containing no visible core, but are entirely composed of zircon with zoning texturally similar to the rims present on grains with cores. We call the inner parts of such zircons “centers” in the subsequent text.
Zircon U-Pb SHRIMP-RG Ages
Out of 120 analyses on zircons from mafic rocks, only 12 provide evidence of inheritance in the zircon core (Table 1). These cores have high U contents (883–1678 ppm) and are highly discordant, yielding a range of 206Pb/238U (204Pb-corrected) dates between ca. 2115 and ca. 52 Ma (Table 1). The rest of the zircon grains yield identical core and rim ages within uncertainty. All the concordant 206Pb/238U rim data yield spot ages between 17.8 and 10.1 Ma (Table 1), with the bulk age of each sample lying between 15.3 and 14.4 Ma regardless of the way in which the age data were calculated (intercept or mean age; Fig. 5. In the remaining text, we refer to the intercept ages.
In sample BH219, 31 zircons yield concordant U-Pb data (Fig. 5A). Rim 206Pb/238U dates range from 15.5 ± 1.5 to 13.6 ± 0.5 Ma, with a lower-intercept age of 14.34 ± 0.18 Ma (MSWD = 0.76). Core 206Pb/238U dates range from 15.5 ± 1.1 to 14.1 ± 0.3 Ma, with a lower-intercept age of 14.08 ± 0.52 Ma (MSWD = 0.56).
In sample BH252, 22 zircons yield concordant U-Pb data (Fig. 5B). Rim 206Pb/238U dates range from 16.2 ± 0.8 to 13.6 ± 0.2 to Ma, with a lower-intercept age of 15.18 ± 0.32 Ma (MSWD = 0.32). Dates for three points from centers range from 15.2 ± 0.7 to 14.2 ± 0.5 Ma, with a lower-intercept age of 16 ± 21 Ma (MSWD = 13).
In sample BH268, 28 zircons yield concordant U-Pb data (Fig. 5C). Rim 206Pb/238U coherent dates range from 16.6 ± 0.3 to 14.3 ± 0.4 Ma, with a lower-intercept age of 14.88 ± 0.22 Ma (MSWD = 2.9, resulting from a large spread of ages along the concordia). Dates from centers range from 16.2 ± 1.2 to 14.6 ± 0.2 to Ma, with a lower-intercept age of 15.13 ± 0.52 Ma (MSWD = 2.2).
In contrast, 74 zircons from country rock samples (BH 205, 210, 211, and 244) yield mostly inherited ages. Middle Miocene–aged rims were found on only a few grains; this rim never exceeded 10 μm in width. Zircons from felsic rocks yield ages that mainly fall into three groups: 1900–1800 Ma, ca. 825 Ma, and 445 ± 12 Ma (MSWD = 3.1); different groups are dominant in different samples (Fig. 6; Table 1).
Trace-Element Geochemistry of Zircons from Mafic Granulites
In the zircons from mafic samples, the Th/U ratios for all the rims with concordant ages range between 0.001 and 0.018, with weighted average of 0.0073 ± 0.0012, whereas the inherited cores have Th/U ratios of 0.0037–0.741. The Zr/Hf ratio is constant around 44 ± 1.
The zircons yield two types of REE patterns (Fig. 7). All the zircon cores that yielded discordant Proterozoic ages have a characteristic magmatic REE pattern with a pronounced negative Eu anomaly (Eu/Eu* = 0.04–0.65). The zircon centers and rims with middle Miocene ages yield a different pattern with no statistical difference between the center and rim. In general, the REE content is low, with a weakly negative to slightly positive Eu anomaly (Eu/Eu* = 0.50–1.74; weighted average 1.14 ± 0.06), and there is a distinct depletion of heavy REEs (YbN/GdN = 0.22 – 6.35; weighted average 0.86 ± 0.13). The positive Ce anomaly is present in all analyzed zircon grains (Ce/Ce* = 8.83–39.44; weighted average 21.7 ± 3.5).
Temperatures for zircon crystallization were calculated using the Ti-in-zircon thermometer (Ferry and Watson, 2007), assuming unit aTiO2 and aSiO2, because rutile, ilmenite, and quartz inclusions were found in the garnets together with zircons. Zircon rims yield Ti-in-zircon temperatures between ∼619 °C and 720 °C with a weighted average of 670 ± 6 °C (data from the zircons yielding a statistically consistent age population; Table 2). The uncertainty in the T value was estimated from the Ferry and Watson (2007) error in empirical correlation. There is a negative correlation between temperature and local Hf content (for the rim data), confirming the coherency of Ti-in zircon data (Fig. A3 [see footnote 1]). No correlation between petrographic setting and Ti concentration was noted, nor is intragrain variability apparent. Miocene-aged centers yield the same weighted average temperature of 681 ± 14 °C (Table 2).
Results of Watson et al. (2006), Ferry and Watson (2007), and Tailby et al. (2011) provide experimental evidence of the pressure effect on Ti solubility in zircon. A computational study of this effect (Ferriss et al., 2008) suggests that at 2 GPa, temperatures should be ∼100 °C higher than at the 1 GPa conditions under which the Ferry and Watson (2007) calibration was performed. We thus use the pressure-corrected T because zircon in the mafic samples is considered to have crystallized at 1.6–1.8 GPa (Fig. 7). The “correction” is only approximate because the zircon thermometer is not yet calibrated for these high pressures, and we assume a simple linear relation between pressure and Ti solubility in zircon. Consequently, the temperatures calculated according to the Ferry and Watson (2007) calibration (Table 2) are increased by 60–80 °C in Figure 7A. This correction thus yields the maximum temperature range for zircon crystallization; these are compatible with the P-T conditions of inferred metamorphic reactions (Fig. 8A) and suggest that the zircons in the mafic granulites crystallized during pregranulite prograde metamorphism. Lower temperatures of zircon crystallization would require a more substantial temperature increase during decompression, which would be difficult to explain. In addition, at lower pressures, the Ti-in-zircon data become increasingly mismatched with the thermobarometric data, and the zircon trace-element data are incompatible with the mineral assemblages.
The mafic granulites may represent vestiges of pre-Himalayan dikes and sills emplaced into the Indian continental margin before the India-Asia collision (Chakungal et al., 2010), but they could also represent dikes and sills emplaced into the lower crust of southern Tibet during collision. Previously published mafic rock bulk-rock geochemistry (Chakungal et al., 2010), in combination with the new zircon age and trace-element data and previously published metamorphic P-T-t constraints (Warren et al., 2011b), provides insight into the origin, evolution, and exhumation history of the eastern Himalayan granulitized eclogites, in particular, and evolution of the Himalayan lower crust in general.
Tectonic Affinity of the Granulite-Bearing Units in the Eastern Himalaya
In NW Bhutan, the granulite unit is in the uppermost Greater Himalayan Sequence, situated in the hanging wall of a middle Miocene out-of-sequence ductile thrust (Warren et al., 2011b) and only 1–3 km beneath the South Tibetan detachment system. The ages of orthogneiss protoliths fall into three age groups: ca. 445 Ma, ca. 825 Ma, and 1800–1900 Ma (Fig. 2; Table 1). Similarly, leucogranites in the area have a significant population of inherited cores (Kellett et al., 2009), mostly Cambrian–Ordovician (ca. 580–456 Ma). The Ordovician zircon protolith age is comparable to other gneisses within the Greater Himalayan Sequence, such as the orthogneiss from the Greater Himalayan Sequence west of the Ama Drime Massif, which yielded an age of 473 ± 16 Ma (Cottle et al., 2009), as well as the 465–470 Ma Namche orthogneiss in the Everest region (Viskupic and Hodges, 2001), a 484 ± 9 Ma “Formation III augen gneiss” in the Annapurna region (Godin et al., 2001), and numerous other Cambrian–Ordovician orthogneisses (Cawood et al., 2007). The Cambrian–Ordovician and Neoproterozoic protoliths are considered typical for the Greater Himalayan Sequence metasediments (Richards et al., 2006), while the Mesoproterozoic and Paleoproterozoic protoliths are rare, or absent, in the Greater Himalayan Sequence rocks but are representative for the Lesser Himalayan Sequence rocks (McQuarrie et al., 2008; Richards et al., 2006).
In the Ama Drime Massif, the tectonic affinity of the unit in which the granulitized eclogites are exposed is ambiguous. The presence of Paleoproterozoic intrusive bodies (Cottle et al., 2009; Li et al., 2003) has been used to support the idea that the granulite/eclogite-bearing rocks belong to the Lesser Himalayan Sequence (Cottle et al., 2009; Groppo et al., 2007) and should therefore be situated beneath the Main Central thrust. However, εNd data instead suggest a Greater Himalayan Sequence affinity (Liu et al., 2007). In Bhutan, the unit with granulitized eclogites is in the hanging wall of an out-of-sequence thrust in the Greater Himalayan Sequence and lies structurally well above the Main Central thrust. Numerical models suggest that Greater Himalayan Sequence rocks exhumed in a dome above an out-of-sequence thrust will have a more distal origin and hence be more “Lesser Himalayan Sequence”–like (Jamieson et al., 2006), which may explain the observed mixing of the Greater Himalayan Sequence and Lesser Himalayan Sequence provenance signatures in the uppermost Greater Himalayan Sequence.
Metamorphic Conditions of Zircon Crystallization
The rare inherited cores in some mafic zircons still preserve Th/U ratios typical for magmatic zircon (Hoskin and Schaltegger, 2003; Rubatto, 2002). In contrast, the low Th/U values and distinctive trace-element patterns of the Miocene-aged zircons suggest solid-state crystallization, and could therefore imply metamorphic growth (Bingen et al., 2004; Hoskin and Black, 2000; Rubatto, 2002; Rubatto and Hermann, 2007). These signatures are not necessarily unique, however, nor are they always indicative of metamorphic versus magmatic, but rather they reflect the local assemblage (Rubatto et al., 2009).
The presence of the analyzed (mafic rock) zircon within garnet (rather than associated with garnet breakdown products), the homogeneous age population, the homogeneous zircon crystallization temperature population, and the lack of significantly chronologically, chemically, and/or texturally different zircon rims suggest that the analyzed zircons grew synchronously with garnet during the prograde metamorphic stage. The presence of large amounts (<25% by volume) of modal garnet could have depleted the rock of heavy (H) REEs prior to, or during, zircon growth (Bingen et al., 2004; Harley and Kelly, 2007; Rubatto, 2002; Rubatto and Hermann, 2007). The lack of a negative Eu anomaly in both zircon and garnet, the flat HREE pattern in the zircons, and the enriched HREE pattern in the garnet (Warren et al., 2011a) are suggestive of synchronous growth of zircon and garnet in a plagioclase-absent assemblage (i.e., an assemblage with no sequestering host for Eu) typical of an eclogite-facies assemblage (Rubatto, 2002; Rubatto and Hermann, 2007). The zircon trace-element patterns could also suggest growth throughout an episode of plagioclase reaction (causing Eu release) during, for example, low-degree melting in the presence of garnet (Rubatto et al., 2009); however, the textural association of zircon within the garnet suggests against this.
Protolith of Mafic Granulitized Eclogites
Even though the present structural position of the granulitized eclogite-bearing unit in the Ama Drime range and in NW Bhutan is still debated, their similarity in metamorphic evolution and geochronology suggests a common origin and a common exhumation process. We now explore options for the origin of the magmas from which the granulitized eclogites originally formed.
The Zr/Hf ratio of 44 ± 1 in zircons from mafic samples is similar to the whole-rock Zr/Hf ratio of 40.9 ± 1.7 (Chakungal, 2006), which is not far from Zr/Hf value of 36 ± 3 estimated for Earth's mantle (McDonough and Sun, 1995). Both major- and trace-element bulk-rock data for the mafic rocks suggest a tholeiitic and alkaline basalt chemistry (Chakungal et al., 2010). Negative Nb and Ti anomalies and an enrichment of Th and light REEs accompanied by highly variable Nd isotopes suggest that they were most likely derived from a subcontinental lithospheric mantle source and were later modified by crustal contamination (Chakungal et al., 2010).
In the Ama Drime range, Cretaceous 110–88 Ma zircon ages (Rolfo et al., 2005) or Neoproterozoic zircon ages of 1057–910 Ma (Liu et al., 2007) and 986.6 ± 1.8 Ma (Cottle et al., 2009) have been interpreted as representing the age of the magmatic protolith, while in NW Bhutan, Paleoproterozoic ages of 1794 ± 58 Ma and 1742 ± 39 Ma (Chakungal, 2006) have similarly been interpreted as representing the timing of original magmatic intrusion. While these ages could be interpreted as representing the timing of magma intrusion, zircon crystallization in a mafic magma is relatively rare. During our study in NW Bhutan, the boudinaged mafic dikes were found intruded into Paleoproterozoic, Neoproterozoic, and Cambrian–Ordovician country rock. Field observations indicate that this is a primary association, not a tectonic one. Our preferred explanation for the pre-Miocene ages in the mafic rocks is thus local contamination by zircons from the country rock.
Zircons separated from the mafic granulites generally yield homogeneous middle Miocene ages; less than 10% of grains contain inherited cores. The center and rim ages yielded by most grains suggest a short (1–2 m.y.) period of growth following initial crystallization (Fig. 3). Zircon morphology and CL patterns indicate dominantly continuous growth rather than resorption or dissolution of precursor grains followed by subsequent metamorphic overgrowth. These data therefore suggest one period of zircon growth, implying either that zircon was not initially present in mafic rocks emplaced before the Himalayan orogeny in Indian continental crust or that zircon crystallized during or soon after emplacement of mafic magma into the lower crust beneath Tibet.
Overall, the data therefore do not distinguish between pre- or syn-orogenic scenarios for emplacement of the protolith mafic magma. However, if the mafic rocks were pre-Himalayan and underwent metamorphism during a previous metamorphic event, then a significant component of inherited cores and a greater range of older zircon ages would be expected, as is the case for the zircons from country rocks. Moreover, mafic magmas do not generally crystallize zircon, and most xenocrystic zircon that was present or assimilated in the parental magmas at the time of emplacement would have dissolved in a high-T, mafic magma. Therefore, the zircon age data are compatible with either scenario of protolith for mafic rocks. Conversely, we can state with certainty, based on the presented data, that the zircons formed during a short Miocene-aged metamorphic event.
At crustal depths over 60 km, eclogite-facies minerals must grow in a slowly cooled gabbro. Zircons in such mafic rocks may have crystallized: (1) directly from the magma, (2) during gabbro to eclogite transformation, (3) by replacement of magmatic baddeleyite by zircon due to increased aSiO2 (e.g., Davidson and van Breemen, 1988), or (4) during subsequent heating and the change to granulite-facies mineralogy (by the breakdown of garnet and/or pyroxene). Scenario 1 is unlikely because mafic magmas do not generally preserve or crystallize zircon, there is no evidence of scenario 3, and scenario 4 is also unlikely because the analyzed zircons are found within the garnet, and not associated with its breakdown.
There are lines of evidence to suggest that a synorogenic Miocene-aged magmatic origin for the mafic granulitized eclogites is a plausible hypothesis. Miocene volcanic and magmatic rocks exposed in the Lhasa block of southern Tibet (Nomade et al., 2004) have been explained by subcontinental lithospheric mantle–derived magmatism (Molnar and Stock, 2009). This magmatic event (1) provides a source for the emplacement of mafic magmas into the lower crust under Tibet (comprising underplated Indian crust at this time), and (2) provides a suitably large-magnitude source of heat to the base of the crust to cause the observed granulite-facies metamorphic overprint, and consequently enough weakening of the felsic host rock to allow exhumation. We therefore suggest that the mafic granulites could plausibly have a middle Miocene–aged magmatic origin.
The E-W extension of the Tibetan Plateau is estimated to have started at 18–13.5 Ma, most likely related to processes occurring beneath the plateau (Blisniuk et al., 2001). Presently, the outward growth of the plateau is explained by the slowing of India-Eurasia convergence and by an additional short-term event such as the removal of mantle lithosphere from northern Tibet since 20 Ma (Molnar and Stock, 2009). This process could have allowed intrusion of hot mantle–derived melt into the lower crust. This magma could have solidified as dikes at various levels within the thickened orogenic crust and also reached the surface in southern Tibet. Consequently, we hypothesize that the present location of the Miocene volcanics may indicate the approximate latitude of the granulitized eclogites at the onset of their exhumation (Fig. 9). Located about 350 km north of the orogenic front, this location in the lower crust of Tibet is consistent with the predicted conditions of lower-crustal eclogitization (e.g., Hetényi et al., 2007).
Exhumation and Cooling Rates
The reliability of estimated exhumation rates is dependent on the confidence with which P, T, and t information is linked, and the confidence with which P may be converted to depth. As previously discussed, absolute P-T estimates for the eclogite-facies stage are imprecise. The zircon trace-element data, Ti-in-zircon temperatures, U-Pb ages, and their textural association in the mafic rocks together suggest zircon crystallization at ∼760 °C and >1.5 GPa. Monazites in the host gneisses yield 15.4 ± 0.8–13.5 ± 0.5 Ma U-Th-Pb ages (Warren et al., 2011b), which overlap the 15.3 ± 0.3–14.4 ± 0.3 Ma zircon ages, suggesting continued metamorphism at high temperatures following eclogite-facies metamorphism.
Monazites associated with melting reactions in the felsic granulites yield younger U-Th-Pb ages of 13.9 ± 0.3 Ma (Warren et al., 2011b). These latter monazites are considered to have crystallized during melting under granulite-facies conditions of ∼750 °C–800 °C and 0.8–1.0 GPa (Warren et al., 2011b). Rapid cooling following peak metamorphism is also suggested by rutile U-Pb ages of ca. 10 Ma (Warren et al., 2011a).
Older monazite cores of older than ∼500 Ma (Warren et al., 2011b) preserved in monazites within the host gneisses suggest incomplete monazite dissolution during high-temperature decompression (Kelsey et al., 2008). The monazite ages in metapelites and orthogneisses with partial melt could therefore alternatively be interpreted as representing the onset of vigorous monazite growth as the rocks cooled and the P-T path intersected the 100% dissolution contour (Kelsey et al., 2008), at ∼700–750 °C and ∼0.5–0.4 GPa (Fig. 8A). Together, the zircon and monazite data imply up to ∼1.2 GPa decompression (from ∼1.6 to ∼1.0–0.4 GPa) during a minimum of 1–2 m.y. (Fig. 8). With a mean crustal density of 2800 kg m−3, this corresponds to exhumation of ∼20 to ∼44 km, and thus a mean exhumation rate of up to ∼2–4.4 cm yr−1. These rates are similar to exhumation rates suggested for buoyancy-driven, subduction-related UHP eclogites in the western Himalaya and other orogens (Baldwin et al., 2008; Parrish et al., 2006; Rubatto and Hermann, 2001).
This short surge of rapid exhumation was followed by a period of slower exhumation and faster cooling rates (Fig. 8B). Zircon U-Pb SHRIMP data from leucogranite with magmatic andalusite and cordierite indicate that the rocks in the structurally highest levels of the granulite unit were at ≤0.28 GPa by 11.0 Ma (Kellett et al., 2009). Rutiles from a mafic granulite and mafic amphibolite from the two juxtaposed units in NW Bhutan cooled through 630 ± 50 °C and 580–680 °C at 11–10 Ma (Warren et al., 2011a). Biotite and muscovite 40Ar/39Ar cooling ages of 11.4–10.7 Ma from nearby granites also suggest rapid cooling at this time (Chakungal, 2006; Maluski et al., 1988, respectively). Numerical models suggest that cooling may have been aided by the combination of strongly compressed isotherms above an expulsed dome (e.g., Jamieson et al., 2006, their Fig. 4) and tectonic denudation along the South Tibetan detachment system, which was active until at least 11 Ma (Kellett et al., 2009).
The granulitized eclogites in Bhutan pose a reaction kinetics conundrum: How could such a substantial overprint of an eclogite-facies assemblage be achieved within such a short period of time (see, for instance, Baxter, 2003; O'Brien, 2008)? The complex, nonequilibrium mineral assemblages observed in these rocks (Fig. 3) (Groppo et al., 2007; Warren et al., 2011b) indicate that the mineral-transforming reactions did not have time to go to completion, and they suggest therefore that the time period was short. More importantly, the rocks in the granulite unit are high-temperature tectonites, and it has been shown theoretically and experimentally that crystal plastic deformation enhances the compositional exchange rates, and therefore mineral reaction rates (e.g., Chakraborty, 2008; Grujic et al., 2011; Stünitz, 1998; Yund and Tullis, 1991).
Estimated exhumation rates for the granulitized eclogites require particular geological constraints: (1) low viscosity conditions to allow high strain rates, (2) a heat source to provide the granulite-facies overprint, (3) a suitable driving force to provide rock uplift, and (4) a suitable structural geometry. All of these conditions are provided in a scenario involving removal of the lower lithosphere beneath Tibet, contact between the hot asthenosphere and mantle lithosphere (Molnar and Stock, 2009), partial melting of the lower crust, and insertion of a cold crustal ramp into the weak, partially molten lower crust (Kellett et al., 2010; Warren et al., 2011b). A crustal ramp with an amplitude of ∼35 km and ∼40° dip to the north has been observed seismically underneath southern Tibet (Hauck et al., 1998) (Fig. 9).
The metamorphic transformation from eclogite-facies to granulite-facies assemblages, the widespread melting in the felsic rocks suggested by the zircon age and trace-element data, and the overprinting textural associations in the mafic rocks suggest an increase in temperature during the initial stages of decompression (Fig. 8B). An input of heat to the lower crust under Tibet during the Miocene (Molnar and Stock, 2009) is one plausible mechanism. Despite the high strain rates implied by the rapid exhumation, strain heating (e.g., Kincaid and Silver, 1996; Nabelek et al., 2010) is not a plausible source for the heat source because the rocks were known to be partially molten at this time, and even a small amount of melt significantly reduces the mechanical strength (Rosenberg and Handy, 2005). Removal of mantle lithosphere remains therefore the most plausible source of heat.
The driving force for exhumation could not have included a significant buoyancy component because the crust was relatively homogeneous (most of the rocks in the granulite-bearing unit are felsic), and there is no evidence of lithosphere-scale extension or a subduction channel underneath the central and eastern Himalaya and southern Tibet during the Miocene. A tectonic driver is therefore needed. Numerical models show that exhumation of deeply buried rocks may be achieved by the insertion of a strong (cold) “plunger” into a layer of weakened material (Warren et al., 2008). Rapid initial exhumation and foreland translation of the Greater Himalayan Sequence could have formed a ramp in the incoming Indian crust (Beaumont et al., 2004). In addition, the suggested slowdown of India-Eurasia convergence between 20 and ca. 10 Ma, optimally at 17 Ma (Molnar and Stock, 2009), would favor return flow in the partially molten crustal layer (for details, see Grujic, 2006, and references therein). The rapid onset of melting (Fig. 8B), and thus weakening during exhumation would provide further positive feedback to this process.
This plunger exhumation model also provides an explanation for the difference in age and degree of Miocene metamorphism in NW Bhutan. Structurally lower, 21–17 Ma amphibolite-facies rocks are overthrust by ca. 14 Ma granulite-facies rocks explained by out-of-sequence thrusting along the Kakhtang/Laya thrust system (Warren et al., 2011b). The formation of this thrust, as well as the cessation of ductile shearing along the South Tibetan detachment system beneath the klippen at ca. 16–15 Ma and the northward shift of ductile shearing along the South Tibetan detachment system until ca. 11 Ma (Kellett et al., 2009), may have occurred in response to a pulse of weak, lower-crustal rocks being exhumed over an incoming ramp of colder Indian crust material (Fig. 9).
In the Ama Drime Massif, the granulite-bearing unit forms an approximately N-S–trending antiform flanked by two oppositely dipping normal faults, shaping the unit into a horst (Fig. 1), similar to the SW portion of the granulite-bearing unit in NW Bhutan bounded in the west by the master fault of the Yadong graben (Fig. 2). Normal faulting is thought to have initiated at around 13 Ma (Cottle et al., 2009), prior to ca. 11 Ma (Kali et al., 2010), and in the northern Yadong graben, its maximum age is 11.5 ± 0.4 Ma (Ratschbacher et al., 2011). It has therefore been suggested that these young metamorphic culminations record exhumation of the middle crust, which was accommodated by N-S–striking normal-sense ductile shear zones kinematically linked to orogen-parallel, crustal-scale, E-W extension, following S-directed flow during the early and middle Miocene (Cottle et al., 2009; Jessup et al., 2008). The E-W extension appears to have been responsible for the upper-crustal stage of the exhumation only: In the Ama Drime Massif, it has accounted for ≤0.6 GPa (22 km) of exhumation (Kali et al., 2010), which followed the initial exhumation of the granulitized eclogites from ca. 30 Ma until ca. 13 Ma, above the Main Central thrust and below the South Tibetan detachment system (Kali et al., 2010). Likewise, the Yadong graben fault, (Cogan et al., 1998) is too shallow to have accommodated exhumation of lower-crustal rocks in the footwall. While there is a clear ductile thrust at the base of the granulite-bearing unit in NW Bhutan, the southern border of the Ama Drime Massif, against the Greater Himalayan Sequence, is not known. Nevertheless, the granulitized eclogites in the Ama Drime range and in NW Bhutan might have been exhumed to the upper crust by the same tectonic process, explaining the startling similarity in their geological setting, metamorphic evolution, and geochronology. The different late Miocene–Quaternary tectonics have, however, altered their middle Miocene structures.
The granulitized eclogite-bearing unit of the Greater Himalayan Sequence in NW Bhutan extends the area of known Himalayan Tertiary HP rocks to an over 250-km-long segment of the Himalayan range. This suggests that the formation of HP material is a significant orogenic stage rather than a local tectonic oddity.
In NW Bhutan, the granulitized eclogites are located in the hanging wall of a middle Miocene out-of-sequence thrust and are located in the youngest and structurally highest unit of the Greater Himalayan Sequence. The mafic rocks, interpreted as strongly boudinaged dikes, are possibly the lower-crustal equivalents of middle Miocene mantle-derived volcanism in southern Tibet. The formation and intrusion of this magma may have been related to the deceleration of India-Asia convergence in the middle Miocene, which might have prompted the removal of mantle lithosphere beneath southern Tibet. This process in turn provided the heat for granulite-facies metamorphism, which eventually facilitated the rapid exhumation of the lower-crustal rocks.
Combined SHRIMP U-Pb, Ti, and REE data from zircons in mafic granulitized eclogites suggest that zircon was crystallized in an eclogite-facies metamorphic assemblage between 15.3 ± 0.3 and 14.4 ± 0.3 Ma. The age data, along with pressure estimates of the eclogite- and granulite-facies stages of metamorphism, and the timing of the onset of cooling, suggest that initial exhumation occurred at a rate of ∼2–4.4 cm yr−1 for 1–2 m.y.
The exhumation of eastern Himalayan lower-crustal eclogites was accommodated by insertion of a midcrustal ramp into a weak lower-crustal layer underlying southern Tibet. This model of eclogite exhumation is unusual in that it involves no crustal extension or buoyancy forces.
Field work in the Kingdom of Bhutan was enabled by the invaluable help provided by the people and the Royal Government of Bhutan, and by the Hoch family. We appreciate the motivating discussions with Lincoln Hollister and Daniela Rubatto and thank John Goodge for editorial handling. Critical comments by Micah Jessup, Franco Rolfo, Daniela Rubatto, and Julia de Sigoyer greatly helped to improve our work. The study was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), Blaustein Foundation (Stanford University), and the National Science Foundation (USA). Warren acknowledges support from the Killam Foundation while at Dalhousie University and the support from NERC Fellowship NER/E0114038/1.