Simultaneous acquisition of detrital zircon Pb-Pb ages and trace element abundances from grains collected across the Indian craton, spanning ∼3 b.y., reveals prominent shifts in Eu/Eu* and light and middle to heavy rare earth element ratios. These shifts correspond to a ca. 3.0–2.2 Ga interval of crustal thickening during Indian craton formation, followed by a period wherein arc magmatism occurred along thinner craton margins from ca. 1.9 to 1.0 Ga, with arc magmatism concentrated along attenuated continental margins after ca. 1.0 Ga. Similar temporal shifts in trace element concentrations are recognized in global whole-rock compilations. We propose that the post–1.0 Ga increase in juvenile magmatism reflects a switch to lateral arc terrane accretion as the primary style of continental growth over the past billion years.

In contrast to thin dense mafic oceanic crust, thick continental crust is composed of buoyant, intermediate to felsic rock. These compositional distinctions govern how modern plate tectonics operate: dense oceanic crust is consumed along subduction zones, whereas long-lived continents resist subduction (Rudnick and Gao, 2003; Korenaga, 2013). Because tectonic processes regulate conditions on Earth’s surface environment, primarily through volcanism and chemical weathering, and the presence of exposed continental crust is critical for maintaining habitability, investigating the chemical and structural evolution of continental lithosphere is essential for understanding Earth system evolution.

Various geochemical proxies have been used to track crustal growth through time (e.g., Taylor and McLennan, 1995; Rudnick and Gao, 2003), including zircon U-Pb age compilations (Campbell and Allen, 2008; Condie and Aster, 2010). Potential preservational biases in the U-Pb zircon record directed attention to zircon Hf isotopic compositions under the assumptions that these data record extraction of melts from the upper mantle (Kemp et al., 2006; Belousova et al., 2010; Voice et al., 2011). While zircons with mantle-like δ18O values avoid signals from crustal recycling (Dhuime et al., 2012; Kemp et al., 2006), thus serving as a viable means to track volumetric extraction of continental crust from the mantle, these data lack information on the compositional evolution of continents via contributions from mixed, highly fractionated melts.

Here we explore the potential for detrital zircon trace element (ZrTE) compositions as a proxy for crustal evolution. Detrital ZrTE data can be used as a provenance tool to distinguish between igneous source rocks (e.g., Belousova et al., 2002; Barth et al., 2013; Grimes et al., 2015); however, these data have not been systematically evaluated with respect to crystallization age over time scales relevant to continent formation. Zircon saturation is elevated in silicic melts, which are often generated along continental arcs (Lee and Bachmann, 2014). The detrital zircon record is therefore strongly influenced by regional arc magmatism (Lee et al., 2016; McKenzie et al., 2016), an important mechanism for crustal addition via water-induced melting of the upper mantle (Rudnick and Gao, 2003). Before the onset of lateral plate tectonics, tonalite-trondhjemite-granodiorites (TTGs) may have been prominent sources of detrital zircon (Moyen and Martin, 2012). However, the detrital zircon record is dominated by zircons younger than 3 Ga (Lee et al., 2016), and so largely reflects secular changes in arc magmatism since plate tectonic initiation (Dhuime et al., 2012, 2015).

Three factors suggest that igneous zircon rare earth element (REE) abundances faithfully preserve a record of the composition of silicic parental melts.

  • (1) REE diffusion in zircon is exceedingly slow (Cherniak et al., 1997), precluding diffusive equilibration after crystallization.

  • (2) Zr is an incompatible trace element of moderately high abundance (∼101–102 ppm) in silica-rich and intermediate composition magmas, meaning that fractional crystallization elevates melt Zr concentrations to the point of zircon saturation for a large range of melt compositions, H2O contents, and, critically, temperatures (Watson and Harrison, 1983); typical arc magmas are expected to reach zircon saturation at temperatures in excess of ∼750 °C (Lee and Bachmann, 2014).

  • (3) Zircon has a propensity to sequester heavy (H) REEs relative to light (L) REEs and middle (M) REEs from the host melt. This is likely driven by the xenotime coupled substitution mechanism (Y + REE)3+ + P5+ = Zr4+ + Si4+ in which the capacity of zircon to accommodate LREE3+ is limited by the large ionic radii of LREEs and attendant lattice strain at the Zr site (Speer and Cooper, 1982). Conversely, the smaller mismatch between ionic radii of Zr4+ and the HREEs accounts for partition coefficients >101 for the REEs Gd and Lu (Hanchar and van Westrenen, 2007). These crystal-chemical controls on partitioning mean that REE abundances in zircon are sensitive to the presence of cogenetic HREE-compatible minerals phases, notably garnet or amphibole, in the parental melt.

Application of combined detrital zircon U-Pb geochronology and TE abundances has been limited by collection of ages and elemental data from different portions of the same grain (Hoskin and Ireland, 2000). Laser ablation split stream–inductively coupled plasma–mass spectrometry (LASS-ICP-MS) circumvents this shortfall by simultaneous collection of U-Pb isotopic and TE abundance data from the same analytical volume, enabling zircon crystallization ages to be linked to melt composition. Provided that ZrTE abundances reflect the primary composition of the parental melt, this technique is ideally suited to assess secular trends in TE chemistry of detrital accessory phases.

Our study focuses on a detrital zircon data set derived exclusively from the Indian subcontinent. Samples are from southern India (Kaldagi Basin), central India (Vindhyan, Aravalli-Delhi, and Marwar sectors), and northern India (Himalaya) (Fig. 1) (see the GSA Data Repository1). LASS analyses, following the methodology presented by Kylander-Clark et al. (2013), were undertaken at the University of Texas (Austin, Texas, USA) using two ThermoFisher Element 2 high-resolution ICP-MS instruments, coupled to a Photon Machines Analyte G.2 ArF 192 nm excimer laser ablation system (www.teledynecetac.com/; see the Data Repository for analytical details). Our ZrTE and U-Pb data set comprises 574 single-grain analyses with Pb-Pb ages between ca. 0.4 and 3.4 G—an ∼3 b.y. record. To assess variations of mean ZrTE concentration with time, we used Monte Carlo bootstrap resampling (see the Data Repository) (Fig. 2). The bootstrap analysis yields an estimate of the average Indian ZrTE composition through time. To avoid inclusions and metamict grains, ZrTE analyses with Ti > 50 ppm and REE + Y > 1 wt% (Hoskin and Schaltegger, 2003) were discarded (n = 82) from the bootstrap analysis.

The ZrTE data set shows the following trends. Eu anomalies (Eu/Eu*), calculated as

forumla decrease in magnitude from 3 to 1 Ga, after which they increase. The effect of oxidation on Eu/Eu* in zircon is subordinate to plagioclase fractionation, which consumes Eu2+ causing a negative anomaly in the residual melt (Trail et al., 2012). The observed Eu/Eu* trend is therefore consistent with suppression of plagioclase crystallization during deep crustal differentiation through time; the increase in Eu/Eu* after ca. 1.0 Ga suggests that there has been a recent increase in the depth of differentiation. We note that the Eu/Eu* trend will also be affected by the H2O content of parental melts, which delays zircon saturation relative to plagioclase crystallization (Keller et al., 2015). Broad increases in mean LREE/HREE and MREE/HREE ratios and Y concentrations occur between ca. 3 and 2.2 Ga, after which values decrease (Fig. 2; Fig. DR1 in the Data Repository). The concomitant increase in LREE/HREE and MREE/HREE ratios and Y concentrations prior to ca. 2.2 Ga is consistent with early deep melting and associated competition for HREEs between zircon and garnet, or amphibole, in parental melts. Mean Th/U decreases monotonically from values between 0.6 and 0.8 before 3 Ga to between 0.3 and 0.5 ca. 500 Ma. This serves as an important check against potential contamination by metamorphic zircons, expected to have low Th/U values (∼<0.1) (Ahrens et al., 1967; Rubatto, 2002; Hoskin and Schaltegger, 2003). Because Th/U in granitoids is strongly influenced by fractionation of accessory phases (Bea, 1996), a decrease in mean zircon Th/U since the Archean could be explained by increasing importance of removal of Th-rich accessory phases in minerals that saturate prior to zircon. Low concentrations of Ce preclude accurate determination of Ce anomalies (Fig. DR1).

Geochemical compilations of intermediate rocks from Phanerozoic continental arcs show that LREE/HREE to MREE/HREE ratios correlate with regional-scale trends in the depth of magmatic diversification (Chapman et al., 2015; Profeta et al., 2015; Farner and Lee, 2017). Elevated LREE/HREE to MREE/HREE ratios are predicted to occur in rocks that have undergone fractionation in the presence of either amphibole- or garnet-rich residues, both of which require pressures >∼1 GPa (Rapp and Watson, 1995). Separating out the individual effects of amphibole and garnet fractionation using geochemical indices alone has proven difficult. Furthermore, melting of subducted oceanic crust can also generate intermediate magmas with low LREE/HREE ratios, due to the LREE-depleted nature of mid-oceanic ridge basalts (MORBs) (Defant and Drummond, 1990). However, covariance between Sr/Y and La/Yb in a large number of intermediate arcs, combined with predictions from thermal models, strongly supports the notion that arc magmatic rocks record geochemical signatures imparted by crustal magmatic processes (Ducea and Barton, 2007; Profeta et al., 2015). Covariance between crustal thickness and REE fractionation of Archean magmatic products is complicated because TTG suites exhibit strongly fractionated LREE/HREE to MREE/HREE ratios (Condie, 1994). Petrogenesis of TTGs is considered to occur in response to either slab melting or by partial melting of a thickened tholeiitic crust (Condie, 1994); the former setting provides an additional mechanism to fractionate REE ratios other than crustal thickness. This caveat aside, trends in detrital zircon LREE/HREE to MREE/HREE ratios potentially represent a spatially and temporally averaged record of crustal thickness through time. We now discuss the ZrTE trends in context of the geologic history of India.

The 3.5–2.4 Ga zircons were likely sourced from a mix of TTGs (mostly older than 3.0 Ga) and subduction-related silicic rocks, and their TE data track formation of the thick stable craton interior coincident with global cratonization (Taylor and McLennan, 1995). Coeval Andean-type continental arc systems formed along the Aravalli belt in west-central India (Buick et al., 2006) and the north Indian margin (Kohn et al., 2010) from ca. 1.9 to 1.6 Ga. Magmatic arc systems along the Eastern Ghats belt (Dasgupta et al., 2013) were probable sources of younger ca. 1.4–1.0 Ga zircons. The 1.0–0.8 detrital zircons, which are abundant in Marwar and Himalaya strata (McKenzie et al., 2011), were likely derived from magmatic arc systems that spanned the Aravalli-Delhi (Just et al., 2011) and the northern margins. Early Paleozoic (0.6–0.4 Ga) Himalayan zircons were likely sourced from arc systems along the north Indian margin (Gehrels et al., 2011). It is possible that some grains were derived from non-Indian terranes during supercontinent amalgamation, but given that source rocks for each age population are known within India proximal to the various basins, grains were probably locally derived. The 1.9–1.0 Ga plutonic rocks are adjacent to cratonic blocks and are known to intrude local basement (Deb and Thorpe, 2004; Buick et al., 2006), whereas the various 0.9–0.4 Ga arc-associated rocks are generally not known to intrude basement or the older Paleoproterozoic magmatic sequences. These Neoproterozoic–Paleozoic arcs were interpreted as outboard arc systems (Gehrels et al., 2011) similar to modern East Asia margins. Accordingly, those magmas were emplaced in transitional or extended continental crust, consistent with the corresponding Eu/Eu* and LREE/HREE trends (Fig. 2).

Our data are similar to global geochemical data sets that indicate a shift to thicker, differentiated continents from ca. 3.0 to 2.4 Ga (Keller and Schoene, 2012; Dhuime et al., 2015; Lee et al., 2016; Tang et al., 2016) (Fig. 3). The LREE/HREE to MREE/HREE reduction after ca. 2.2 Ga could have been caused by decreased crystal-liquid fractionation with melts penetrating relatively thinner crust (Keller and Schoene, 2012). Dhuime et al. (2015) attributed decreased Rb/Sr whole-rock values from 1.0 to 0.2 Ga (Fig. 3) to increased crustal erosion such that continental destruction outpaced crustal generation. We find that similar coeval shifts in Indian ZrTE Eu/Eu* and LREE/HREE to MREE/HREE ratios correspond with magmatism in outboard arc systems, rather than direct emplacement within the thickened craton. Incompatible/compatible element ratios (Rb/Sr) should track magmatic processes, such as liquid-crystal fractionation, rather than net volumetric changes in continental crust (Lee and McKenzie, 2015); therefore, we postulate that the post ca. 1.0 Ga geochemical trends represent an overall increase in juvenile magmatism and lateral accretionary continental growth. The western margins of North America and South America are largely composed of accreted terranes with thin juvenile basement (e.g., Coney et al., 1980; Ramos, 2009). By 2.5–2.2 Ga, cratons may have reached some critical thickness that hindered magmatic addition (Dhuime et al., 2015). Subsequent changes in Earth’s thermal boundary conditions may have influenced tectonic processes (Korenaga, 2013), promoting slab rollback and magmatism in thin attenuated continental margins, increasing lateral aggradation via arc-terrane accretion. This serves as an alternative explanation to the hypothesized rate increase of crustal destruction over the past ∼1 b.y., which is further challenged by the difficulty of foundering buoyant cratonic material into the mantle.

Despite observations that show that ZrTEs exhibit intrasample variability (e.g., Grimes et al., 2015), general agreement between Indian data and whole-rock compilations (Fig. 3) suggests that detrital ZrTE data have the potential to preserve a signal of parental melt compositions. Our study demonstrates how U-Pb and TE detrital zircon data acquired by LASS-ICP-MS provide a tool to rapidly generate records of crustal processes. By filling in petrogenetic gaps in the whole-rock record of crustal composition, this approach will enhance our understanding of tectonic and crustal evolution throughout Earth history.

We thank T. Kemp and anonymous reviewers for constructive comments. Support was provided by U.S. National Science Foundation Tectonics grant EAR-1348031 to McKenzie and Stockli, and Hong Kong RGC-ECS (Research Grants Council–Early Career Scheme) grant 27305417 to McKenzie. Smye and McKenzie acknowledge support from the Jackson Postdoctoral Fellowship Program (University of Texas at Austin).

1GSA Data Repository item 2018077, methods and data tables, is available online at http://www.geosociety.org/datarepository/2018/ or on request from [email protected].