Do Supercontinent-Superplume Cycles Control the Growth and Evolution of Continental Crust?

Characterizing the secular changes in the pattern of crustal growth and its reworking is essential for understanding how the Earth system works (Arndt and Davaille, 2013; Dhuime et al., 2012; Rolf et al., 2012; Wang et al., 2012, 2009; Lancaster et al., 2011; Belousova et al., 2010; Li and Zhong, 2009). However, models for continental growth are still highly debated and commonly contradictory (Dhuime et al., 2012; Voice et al., 2011; Belousova et al., 2010; Condie et al., 2009; Parman, 2007; Kemp et al., 2006; Condie, 1998; Taylor and McLennan, 1995; McCulloch and Bennett, 1994; Allègre and Rousseau, 1984; Armstrong, 1981). Estimating the proportion of new crustal addition versus reworking of older crust through time is crucial for addressing this issue. However, it is still a challenge to distinguish new continental crustal addition from reworking of older crust using geochemical and isotopic indices (e.g., Iizuka et al., 2013; Dhuime et al., 2011). In this paper, we critically assess all parameters available for study and provide a revised approach for defining juvenile continental crust. We place a strong emphasis on the mechanisms involved in the generation of juvenile continental crust and how to use zircon data to constrain crustal growth. Our analysis demonstrates that the evolution of continental crust can be directly linked to whole of Earth geodynamics that control the supercontinent-superplume cycles.

The paradox of the Earth's continental crust is that although it is ultimately derived from the mantle (Hofmann, 1988), its andesitic (SiO₂=60.6 wt.%; Rudnick and Gao, 2003) to dacitic (SiO₂=65.2 wt.%; Hacker et al., 2011) bulk composition contrasts with mantle-derived melts that are dominantly basaltic in composition. There are three main hypotheses to account for this paradox. The first (the basalt-input model) postulates that juvenile continental crust originally had a basaltic composition formed by melting of peridotitic mantle, which then evolved to an andesitic composition through stripping off the mafic-ultramafic high- density components from the lower continental crust into the convective mantle by foundering processes (also called delamination; Rudnick, 1995; Kay and Kay, 1993; Arndt and Goldstein, 1989). This was accompanied by uprise of hot asthenospheric mantle, resulting in remelting of the thinned juvenile crust to produce rocks of more felsic compositions. The second is the island-arc model which argues that juvenile continental crust has an andesitic composition because it was generated through the accretion of island arcs (Castro et al., 2013; Kelemen et al., 2007, 2003; Kelemen, 1995; Taylor, 1967). The third model proposes that juvenile continental crust is mainly formed by partial melting of a mixture of subducted oceanic basalts and sediments (Castro et al., 2013, 2010; Niu et al., 2013).

Because the major element composition of the bulk continental crust is more felsic than that of mantle melts (Rudnick and Gao, 2003; Taylor, 1967), the basalt-input model requires a mechanism to convert an initial basaltic composition into an andesitic composition via the elimination of an unseen ultramafic residue, which is commonly considered to have sunk into the underlying mantle (Gao et al., 2004; Rudnick, 1995; Kay and Kay, 1993; Arndt and Goldstein, 1989). However, this is inconsistent with thermo-mechanical modeling, the composition and mineral assemblage of lower crustal xenoliths, and laboratory experimental results (Castro et al., 2013, 2010; Hacker et al., 2011). The island-arc model also has many shortcomings (see summary in Niu et al., 2013), including: (1) the dominant primitive mantle input into an arc is basalt, not andesite (Niu et al., 2013; Gao et al., 2004; Rudnick, 1995; Kay and Kay, 1993; Arndt and Goldstein, 1989); (2) arc magma has variable but large Sr excesses, whereas the bulk continental crust is slightly Sr deficient (Niu et al., 2013); and (3) arc crust production is mass-balanced by subduction erosion and sediment recycling, thus contributing no net gain to the continental crust, at least in the Phanerozoic (see Niu et al., 2013; Stern, 2011 for reviews). The third model, based on partial melting of subducted oceanic crust, emphasizes the role played by recycling of mafic crust.

Before discussing the growth and evolution of continental crust, it is necessary to define what is meant by 'juvenile crust'. Essentially, this is young crust that has a short crustal residence time. Because the crustal residence time is mainly constrained by radiogenic isotopes, in most cases 'juvenile continental crust' refers to its characteristic juvenile radiogenic isotopes. Because the formation of continental crust is ultimately dependent on mantle input, the primary crust should be in equilibrium with peridotitic mantle in terms of both major element and radiogenic isotope compositions.

The earliest preserved continental crust (3.8–2.5 Ga) is mainly composed of tonalite-trondhjemite-granodiorite (TTG) (e.g., Moyen and Martin, 2012; Martin et al., 2005; Rapp et al., 2003, 1991; Foley et al., 2002; Rapp and Watson, 1995; Martin, 1986). The origin of TTG has been extensively studied and discussed over several decades, with ideas ranging from formation by shallow mantle melting and fractional crystallization, partial melting of eclogite/amphibolite, melting of subducted oceanic plateaus, and incipient melting of hydrous and fertile upper mantle peridotite (Martin et al., 2014; Moyen and Martin, 2012; Ernst, 2009; Rapp et al., 2003; Foley et al., 2002; Condie and Abbott, 1999). However, the composition of Archean TTGs is very similar to that of typical Phanerozoic subduction-related igneous rocks, including adakites (Castillo, 2012; Defant and Drummond, 1990), sanukitoids (Tatsumi, 2008), and high-Mg andesites (Kawabata and Shuto, 2005; Kelemen, 1995). Such a compositional similarity has led to the popular view that TTGs must have formed in subduction settings (Martin et al., 2014, 2005; Rapp et al., 2003, 1991; Foley et al., 2002; Rapp and Watson, 1995; Martin, 1986). Furthermore, the types of rocks present in Archean crust, their patterns of deformation, their magmatic and metamorphic features, and their trace element geochemical signatures suggest that geodynamic processes similar to those operative in Phanerozoic subduction zones played a key role in the growth of continental crust through time (Polat et al., 2011, 2002; Furnes et al., 2009, 2007; Jenner et al., 2009; Nutman and Friend, 2009; Polat and Münker, 2004). In addition, high ambient mantle temperatures during Earth's early history would more readily facilitate partial melting of oceanic crust. Thus, Archean TTGs can be produced by partial melting of oceanic crust (Martin et al., 2005; Foley et al., 2002; Rapp and Watson, 1995; Drummond and Defant, 1990; Martin, 1986) or mafic lower crust in tectonically or magmatically thickened arcs or oceanic plateaus (Martin et al., 2014; Smithies et al., 2009, 2003; Bédard, 2006; Condie, 2005; Rapp et al., 2003; Smithies, 2000).

Proterozoic continental crust was generated mainly through addition of vast granitic batholiths and their eruptive equivalents at convergent plate boundaries (e.g., Zhang et al., 2018; Castro et al., 2013; Niu et al., 2013; Spandler and Pirard, 2013; Hacker et al., 2011; Jagoutz et al., 2011; Jagoutz, 2010; Kröner et al., 2007; Jahn, 2004; Jahn et al., 2004; Rudnick, 1995; Rudnick and Fountain, 1995; Taylor and McLennan, 1995). In the Phanerozoic, the Central Asian Orogenic Belt (CAOB) (Tang et al., 2012; Kröner et al., 2007; Jahn, 2004; Jahn et al., 2000; Wu et al., 2000) and the Greater Tibetan Plateau (Niu et al., 2013 and references therein) represent the two largest belts of juvenile continental crust on Earth, and likewise contain abundant granitoids. The generally felsic chemical composition of the crust precludes direct peridotitic mantle input in the formation of such huge volumes of granitoids (Niu et al., 2013; Chappell, 1996a, b; Rudnick, 1995; Hamilton et al., 1980). However, the radiogenic isotopic (Sr-Nd- Pb) features of rocks in the CAOB and Tibetan Plateau suggest a large depleted mantle-like input in their sources (Niu et al., 2013; Bonin, 2007; Wu et al., 2000). For example, recently-compiled whole-rock Nd isotope data (Yang et al., 2017) suggest that about 64% of the continental crust in the southeastern CAOB is juvenile, but more than 80% of these juvenile igneous rocks have felsic major element compositions. This implies that the contrast between felsic major element composition (not in equilibrium with the mantle) and strongly depleted mantle-like radiogenic isotope signatures is a primary feature of juvenile continental crust. Although experimental evidence shows that water-saturated melting of peridotite (H₂O≥6 wt.%) with residual refractory harzburgite can produce andesitic melts (e.g., Hirose, 1997), andesite production by such hydrous peridotite melting is exceptionally rare (see Straub et al., 2013 for a review).

In-situ analysis of Hf and O isotopes in single zircon grains from the 425–390 Ma Lachlan Ⅰ-type suite of Eastern Australia indicates 50% to 85% of juvenile contribution in their generation (Kemp et al., 2007). However, the felsic whole-rock chemical composition of Ⅰ-type granitoids (Chappell, 1996a, b; Hamilton et al., 1980) implies that such a high proportion of juvenile components cannot be derived directly from mantle melting. Moreover, recent geochemical studies show that juvenile continental crust in the vast juvenile granitoid provinces of the CAOB (Xu et al., 2013; Tang et al., 2012) and the Greater Tibetan Plateau (Niu et al., 2013) may have been formed by partial melting of subducted oceanic basalts, plus variable amounts of sediment (Wei et al., 2002). In addition, lower crustal xenoliths, including granulites and pyroxenites from the Cenozoic Hannuoba basalts (North China Craton), also provide evidence for continental crustal growth by melting of recycled oceanic slabs (Liu et al., 2005; Xu, 2002).

Current knowledge about global-scale differentiation of the bulk silicate Earth and Hadean–Eoarchean crust is mainly derived from ancient zircon grains and short-lived ¹⁴⁶Sm-¹⁴²Nd isotope data. The composition and origin of Earth's early crust remains hotly debated. Mafic Hadean proto-crust was inferred to be present at ≥4.3 Ga based on ¹⁴²Nd isotope data (Rizo et al., 2013, 2012, 2011; O'Neil et al., 2008). Likewise, the Jack Hills detrital zircon suite provides evidence that 4.4–4.0 Ga felsic continental crust was possibly derived from a recycled and enriched basaltic progenitor (Blichert-Toft and Albarède, 2008) or consolidated magma ocean of mafic composition (Kemp et al., 2010). Howerver, the inferred equilibrium melts based on partition coefficients from 4.3–3.3 Ga Jack Hills zircon grains are similar to andesites formed in modern subduction settings (Turner et al., 2020), and the Hf data support the long-term recycling of original Hadean crust into the Neoarchean (Kemp et al., 2010).

Numerical modeling (Zhu et al., 2013; Vogt et al., 2012), geochemical (Castro et al., 2013, 2010; Goss and Kay, 2009; Plank et al., 2009; Wood and Turner, 2009; Gómez-Tuena et al., 2008; Klimm et al., 2008; Hermann and Spandler, 2007; Kawabata and Shuto, 2005; Tsuchiya et al., 2005; Beate et al., 2001; Kelemen, 1995), and experimental (Behn et al., 2011; Rapp and Watson, 1995) studies suggest that melting of a mixture of subducted oceanic basalt and sediment may have played a key role in continental crustal growth throughout Earth's history. Basaltic oceanic crust is the major product of mantle melting and would be recycled into the mantle convection system or be involved in building juvenile continental crust shortly after extraction from its mantle source (≤200 Ma). Therefore, the short residence time and basaltic major element composition of oceanic crust ensures that its derivative melts acquire a depleted mantle-like radiogenic isotopic signature. However, its oxygen isotopic signature will vary from normal mantle-like to higher or lower values, depending on the degree and nature of water/rock interaction as the basalts interact with seawater. Subsequent melting will lead to the production of melts with andesitic to felsic major element compositions.

Recognition of primitive end-member sources that contributed to build continental crust (Castro et al., 2013, 2010; Niu et al., 2013; Yang et al., 2007a, b, 2004; Gray, 1984; DePaolo, 1981a) is crucial for understanding the growth and evolution of continental crust. However, homogenization is usually very common and effective in obscuring the discrete origin of different magma batches (Barbarin, 2005; Poli et al., 1996; DePaolo, 1981b). This hampers the usage of whole-rock composition to decipher the extent of juvenile addition versus crustal reworking in the generation of granitoids. Furthermore, although source mixing may play an important role in producing continental igneous rocks, the highly complex nature of the processes involved usually make it impossible to identify the end-member sources of most granitic rocks. Nonetheless, magmatic zircon can preserve details of Hf-O isotopic variations in the end-member magma sources (Yang et al., 2007a). In-situ zircon U-Pb dating and Hf-O isotope analyses together provide a way for tracing continental crustal growth (Li et al., 2009b; Kemp et al., 2007; Yang et al., 2007a).

Zircon grains that retain precise and concordant U-Pb ages can provide not only a record of magmatic and thermal events (Hawkesworth et al., 2010), but also insights into continental crustal growth and its evolution (Dhuime et al., 2012; Wang X-C et al., 2012, 2009; Belousova et al., 2010; Hawkesworth et al., 2010; Condie et al., 2009; Kemp et al., 2007, 2006; Iizuka et al., 2005; Condie, 1998). They can provide at least three important pieces of information: U-Pb (or Pb-Pb) ages, Hf-O isotopic signatures, and Hf model ages. Zircon oxygen isotopes have been widely utilized as a diagnostic index for recognizing juvenile continental crust (Cawood et al., 2013; Iizuka et al., 2013; Wang F Y et al., 2013; Dhuime et al., 2012; Wang C Y et al., 2011, 2009; Hawkesworth et al., 2010; Kemp et al., 2007; Hawkesworth and Kemp, 2006a, b) and both zircon U-Pb crystallization ages and Hf model ages have also been applied to investigate the growth and reworking of continental crust. However, care is needed in interpreting these data and below we provide a new approach to evaluate how and when continental crust was generated.

With rare exceptions, the oxygen isotopic ratio of ¹⁸O/¹⁶O (usually expressed as δ¹⁸O values, reported relative to Vienna Standard Mean Ocean Water, VSMOW) is remarkably homogeneous in the upper mantle with a mean value of 5.3‰±0.3‰ (1σ) and has remained so throughout Earth's history (Valley et al., 1998; Mattey et al., 1994). In contrast, supracrustal materials or their recycled components usually have highly variable δ¹⁸O values (to +30‰) (e.g., Bindeman, 2008; Simon and Lécuyer, 2005; Valley et al., 2003; Eiler, 2001). Low-temperature alteration produces high δ¹⁸O values, whereas high-temperature alteration produces low δ¹⁸O values, with negative δ¹⁸O values usually related to glacial events (e.g., Bindeman, 2008; Eiler, 2001). Magmas that have incorporated these materials through crustal assimilation, magma or source mixing will therefore show a range of values reflecting the mixing proportions when zircon crystallises (e.g., Bindeman, 2011, 2008; Wang X-C et al., 2011; Eiler, 2001). For instance, hydrothermal alteration is an important feature of the upper part of the oceanic lithosphere (pillow lavas, sheeted dykes, and gabbros) and can result in δ¹⁸O values ranging from about 2‰ to 16‰, with peak values of 9‰–12‰ for pillow lavas, 6‰–11‰ for sheeted dykes, and 4‰–6‰ for gabbros (Yamaoka et al., 2012; Miller et al., 2001). Thus, juvenile continental crust formed through partial melting of subducted upper oceanic crust can possess a large range of δ¹⁸O values that bear no resemblance to the average mantle value. Juvenile crustal materials that grew from sub-arc mantle may also possess a large range of oxygen isotopic compositions that plot outside the mantle zircon δ¹⁸O range of 4.5‰–6.5‰ (Valley et al., 1998) because sub-arc mantle can be significantly modified by the injection of small proportions of slab-derived fluids (Prouteau et al., 2001; Eiler et al., 2000, 1998). Therefore, the assumption that normal mantle-like oxygen isotopic signatures are both diagnostic of, and necessary for the identification of juvenile continental crust is highly questionable. A normal mantle-like zircon oxygen isotope signature may only provide an additional constraint on identifying juvenile continental crust when the targeted igneous rocks did not experience water-rock reaction at either low or high temperatures.

This conclusion is confirmed by a global compilation of chemical and isotopic data from Late Mesozoic to modern igneous rocks generated at convergent plate boundaries. As shown in Fig. 1, the majority (> 50%) of arc rocks with < 57 wt.% SiO₂ have zircon δ¹⁸O values that plot away from the mantle field. Because metasomatized sub-arc lithospheric mantle has an enriched δ¹⁸O value of δ¹⁸O=8.3‰±0.28‰ (1σ) (Liu et al., 2014), generation from this source can account for the high δ¹⁸O signatures of arc basalts. Although crustal contamination could also result in a high δ¹⁸O in evolved magmas, it would require more than 25% crustal input if supracrustal end-member and parental magma have δ¹⁸O values of 16‰ and 5.3‰, respectively. Input of such a high proportion of crustal materials should result in high silica contents and unradiogenic isotopes, but this is inconsistent with the basaltic to andesitic major element compositions and juvenile radiogenic isotopes of arc basalt and andesite as shown in Fig. 2. Although it is not possible to strictly exclude a contribution from crustal contamination for any given example, this observation demonstrates that the large range of oxygen isotope values is most likely a primary feature of the sub-arc mantle. Therefore, non-mantle-like zircon oxygen isotopic signatures can be present in juvenile crust that grew from sub-arc mantle. Consequently, mantle-like zircon oxygen isotopic signature is not a diagnostic index for juvenile crust that grew from sub-arc mantle. Furthermore, the similarities in high δ¹⁸O values and other geochemical features between modern adakite (Bindeman et al., 2005) and Archean TTG (Hiess et al., 2009; Valley et al., 2005) support the hypothesis that recycled oceanic crust played an important role in the generation of Archean TTG (Martin et al., 2005; Defant and Drummond, 1990; Martin, 1986).

Figure 1.  Plot of re-calculated zircon δ¹⁸O values against SiO₂ for modern to Late-Mesozoic arc magmas from the Banda arc (Nebel et al., 2011; Downes et al., 2001), the Lesser Antilles arc (Davidson and Wilson, 2011; Bouvier et al., 2008a; van Soest et al., 2002), the Izu-Bonin intra-oceanic arc (Macpherson and Mattey, 1998; Harmon and Gerbe, 1992), the Aleutian arc (Finney et al., 2008), the Japan arc (Kimura et al., 2002), the Central Mexican volcanic belt (Johnson et al., 2009), the Sunda arc (Gertisser and Keller, 2003), the intra-oceanic Kermadec arc in the Southwest Pacific (Barker et al., 2012), the Lassen continental arc (Feeley et al., 2008), the Kamchatka arc (Bindeman et al., 2004; Dorendorf et al., 2000; Pineau et al., 1999), and the Andes arc (Hidalgo et al., 2012). With the exception of oxygen isotope data from Tamura et al. (2009) obtained on fresh glasses, all other data were obtained from phenocrysts, including olivine, pyroxene, plagioclase, and amphibole. Samples that underwent alteration were stripped off according to the original case studies. The melt oxygen isotope values (δ¹⁸O) were calculated from phenocryst δ¹⁸O data using theoretical and empirical mineral-melt fractionation factors. The calculated melt δ¹⁸O values and whole-rock silica contents were then used to calculate zircon δ¹⁸O values according to δ¹⁸O (zircon)–δ¹⁸O (whole rock)≈ -0.061 2×(wt.%SiO₂)+2.5 (Valley et al., 2005). The mantle-like zircon δ¹⁸O range of 4.5‰–6.5‰ (Valley et al., 1998) is highlighted by the grey band. The weighted mean value of mantle zircon δ¹⁸O=5.3‰±0.3‰ (Valley et al., 1998) is shown by the blue dashed line. The vertical red line separates basalt and basaltic andesite (SiO₂ < 57 wt.%) that were derived mainly from a mantle source, from andesite and rhyolite (SiO₂≥57 wt.%) that were mainly melts of a crustal source (oceanic or continental). Juvenile continental crust formed by arc magmas should display mantle-derived low-silica and high-ε_Nd signatures. In this study, magmas with SiO₂ < 57 wt.% and ε_Nd > 0 are regarded as representing juvenile continental crust.

Figure 2.  Whole-rock Hf isotopic composition of juvenile crust versus magma formation ages. The Hf isotopic composition of juvenile continental crust is defined by juvenile arc magmas (Kröner et al., 2013; Li J X et al., 2013; Mohan et al., 2013; Zeh et al., 2013; Choulet et al., 2012; Muñoz et al., 2012; Bouilhol et al., 2011; Yang et al., 2011; Matteini et al., 2010; Li X H et al., 2009a; Polat and Münker, 2004; Samson et al., 2003), plume-related basalts, including the Hawaiian OIB (Frey et al., 2005; Stracke et al., 2000), the Siberian continental flood basalts (Carlson et al., 2006), the central Kerguelen Archipelago flood basalts (Xu et al., 2007), the Iceland flood basalts (Søager and Holm, 2011), the Bikou Neoproterozoic continental flood basalts (Wang et al., 2008), the Deccan continental flood basalts (Sen et al., 2009), oceanic plateau lavas (Timm et al., 2011; Hastie et al., 2008; Thompson et al., 2004), and komatiites (Blichert-Toft and Puchtel, 2010; Blichert-Toft et al., 2004, 1999); other juvenile continental crustal materials include plume-related syenite (Zhong et al., 2009), Mesoarchean to Mesozoic mantle-derived rocks and Eoarchean gneiss from West Greenland (Vervoort and Blichert-Toft, 1999). The Hf isotopic evolution curves for the depleted mantle (Griffin et al., 2000), new continental crust (NCC; Dhuime et al., 2011), and sub-arc mantle (grey zone; Iizuka et al., 2013) are shown. The range of modern sub-arc mantle is constrained by Neogene island arc basalts (Iizuka et al., 2013). The Hf isotopes defining the upper limit for juvenile continental crust can be divided into two stages. The first stage is 4.56–3.5 Ga derived from a super-chondritic bulk silicate Earth with ¹⁷⁶Lu/¹⁷⁷Hf=0.037 5 (Caro and Bourdon, 2010) and evolving along a curve with ¹⁷⁶Lu/¹⁷⁷Hf=0.037 5+0.007 6t, where t=0 indicates the time of Earth formation. The second stage is 3.5–0 Ga that evolved along a curve with constant ¹⁷⁶Lu/¹⁷⁷Hf=0.038 5 and ε_Hf(3.5 Ga)= +7.3. The lower limit for juvenile crust is defined by ε_Hf(0)= +3.5 and ε_Hf(4.56 Ga)=0. The value of Δε_Hf(t)≥0, where Δε_Hf(t)=ε_Hf(t)_sample–ε_Hf(t)_LJCC, ε_Hf(t)_LJCC, indicates the lower limit of Hf isotopic composition of juvenile continental crust at time t. The present-day chondritic parameters of ¹⁷⁶Lu/¹⁷⁷Hf=0.033 6 and ¹⁷⁶Hf/¹⁷⁷Hf=0.282 785 (Bouvier et al., 2008b) were used to calculate ε_Hf(t). The Hf isotope evolution of protoliths extracted from a mantle source with ¹⁷⁶Lu/¹⁷⁷Hf of 0.022 (mafic lower crust; Amelin et al., 1999), 0.015 (average continental crust; Griffin et al., 2000), and 0.009 3 (felsic upper continental crust; Amelin et al., 1999) are also shown. The inset histogram shows the Hf isotope range of MORB with a weighted mean present-day ε_Hf value of 13.6±6.6, based on MORB data compiled by Meyzen et al. (2007).

The above discussion demonstrates that juvenile continental crust can display a large range of oxygen isotope compositions. Zircon grains formed from melts generated by partial melting of subducted oceanic crust or metasomatized sub-arc mantle can therefore possess both mantle-like and non-mantle-like oxygen signatures. On the other hand, because oxygen isotopes of igneous rocks are ultimately related to water-rock interaction, felsic melts derived from partial melting of ancient pristine middle to lower crustal rocks that have not undergone interaction with surface water can still possess mantle-like oxygen signatures. In summary, although oxygen isotopes are sensitive to recycling of materials that were once modified by water-rock interaction at shallow depth, they cannot be used as a diagnostic index for distinguishing crustal addition from the reworking of pre-existing continental crust. This shows that oxygen isotopes must be integrated with other parameters (e.g., Hf isotopes) to determine if the continental crust is juvenile.

Determining the Hf characteristics that can be applied to define juvenile continental crust is crucial for distinguishing between crustal growth and crustal reworking. The Hf isotopic evolution of the depleted mantle, based on mid-ocean ridge basalts (MORB), has generally been accepted for depicting the Hf characteristic of juvenile crust (e.g., Griffin et al., 2000; DePaolo, 1981b). However, modern continental crust is dominantly created at convergent margins and not spreading centers, so the isotopic composition of juvenile island arc crust is different from that of the depleted mantle (e.g., Iizuka et al., 2013; Dhuime et al., 2011). Hf isotope compositions of modern juvenile continental crust have also been constrained by Neogene island arc basalts, which have a large range of ¹⁷⁶Hf/¹¹⁷Hf isotope ratios, varying from 0.283 271 (ε_Hf= +17.2) to 0.282 913 (ε_Hf= +4.5) (Iizuka et al., 2013). This emphasises that Hf isotope systematics of juvenile continental crust have a range of values rather than a specific value.

Combined whole-rock and in-situ zircon Lu-Hf isotope data for a large number of Archean TTG suites shows that early continental crust most likely grew from nearly primordial unfractionated materials extracted from the deep mantle via rising plumes that left a depleted melt residue in the upper mantle during the Archean (Guitreau et al., 2012). Therefore, the Hf isotopic composition of early continental crust was different from that of modern depleted mantle (Griffin et al., 2000) or modern arc mantle (Dhuime et al., 2012). This implies that applying present-day average values back into Earth's early history may not adequately characterise the isotopic composition of early juvenile continental crust. Knowledge about the Hf isotopic composition of juvenile continental crust and its change through time is thus a fundamental requirement for distinguishing crustal growth from reworking.

Because juvenile continental crust is ultimately derived from the mantle (Hofmann, 1988), the radiogenic isotope composition and evolution of mantle melts provide the best candidates for characterizing juvenile continental crustal throughout Earth's history. A global compilation through time of whole-rock ¹⁷⁶Hf/¹⁷⁷Hf isotope data from juvenile arc magmas, plume-related flood basalts and Archean komatiites and plume-induced syenites that were taken to represent juvenile continental crustal components generated by plume activities, and Mesoarchean to Mesozoic mantle-derived rocks (Fig. 2) therefore provides a characterisation of the Hf isotopic signature of juvenile continental crust. Figure 2 reveals two striking features. One is that juvenile continental crust displays a large range of ε_Hf(t) values throughout Earth's history. This again demonstrates that the widely-used Hf isotopic evolution curves based on either MORB (red dashed line) or the average of modern arc basalts (blue line) do not accurately represent the Hf isotopic characteristics of juvenile continental crust. The other prominent feature is that the highest ε_Hf(t) values define a two-stage Hf isotopic evolution curve (solid green line in Fig. 2). Here, the ε_Hf(t) values appear to increase rapidly from +2.0 at about 4.0–3.9 Ga to +7.5 at about 3.5 Ga. After that, the increase in ε_Hf(t) values is more gradual and is consistent with the depleted mantle growth curve, with constant ¹⁷⁶Lu/¹⁷⁷Hf ratios from about 3.5 Ga to the present-day (Fig. 2). The lower limit of ε_Hf(t) values for juvenile continental crust (dashed green line in Fig. 2) is defined by the lowest values of less depleted mantle-derived basaltic rocks, including the 25.7 Ma Mt. Capitole flood basalts in the central Kerguelen Archipelago (Xu et al., 2007), the 880 Ma Caribbean Plateau basalts (Thompson et al., 2004), the 252 Ma Siberia continental flood basalts (Carlson et al., 2006), the 820–810 Ma Bikou continental flood basalts in South China (Wang et al., 2008), the 3.5 Ga Barberton komatiites (Puchtel et al., 2013; Blichert-Toft and Arndt, 1999), juvenile arc magmas from post-collisional ore-bearing adakitic porphyries in southern Tibet (Li et al., 2011), the Vijiayan arc complex in Sri Lanka (Kröner et al., 2013), and the Goiás arc in Central Brazil (Matteini et al., 2010). Such a lower limit is described by the line of ε_Hf(t)=3.5–0.766 5T in Fig. 2, where T is in Ga and starts from the present-day. This lower limit of ε_Hf values for juvenile continental crust matches well with the lower limit values of MORB (present-day ε_Hf=+3.6; Fig. 2). Because this approach integrates both mantle heterogeneity and isotopic evolution of direct mantle-derived basaltic material, it can provide the first-order constraint on Hf isotope characteristics of strictly-defined juvenile continental crust.

In this study, zircon grains that have Hf isotope compositions plotting within the field defined by the dashed and solid green lines in Fig. 2 are regarded as being derived from juvenile continental crust. These grains are characterized by Δε_Hf(t)≥0, where the Δε_Hf(t) value was calculated by Δε_Hf(t)=ε_Hf(t)_sample– ε_Hf(t)_LJCC, with LJCC indicating the lower limit of juvenile continental crust defined by ε_Hf(0)= +3.5 and ε_Hf(4.56 Ga)=0 (dashed green line in Fig. 2). It should be noted that this method cannot completely exclude rapid juvenile crustal reworking because the source rock(s) of the bulk of these zircon grains is unknown. Fortunately, this weakness can be overcome by combining geological mapping, whole-rock (element and isotope) analyses, and tectonic constrains to evaluate such a scenario. Therefore, zircon age peaks defined by grains with Δε_Hf(t)≥0 can be used to constrain the timing of formation of juvenile continental crust.

The uncritical use of detrital zircon Hf model ages is questioned here because of the large uncertainties involved in their interpretation, as outlined below.

Firstly, the Lu/Hf ratio is the key parameter for calculating zircon Hf model ages. This calculation requires knowledge of the Lu/Hf ratios of both the host rock and the mantle reservoir from which this material was derived. Juvenile radiogenic isotopes can be recorded by both mafic and felsic rocks (Fig. 2), but they are significantly different. For example, the mafic lower continental crust and felsic upper continental crust have Lu/Hf=0.022 and 0.009 3, respectively (Amelin et al., 1999), whereas the bulk continental crust is considered to have Lu/Hf=0.015 (Griffin et al., 2000), based on an almost 50 : 50 mixture of lower and upper crust. Furthermore, the Lu/Hf ratios of mantle reservoirs are also a deciding factor in the calculation of zircon Hf model ages. This is because the ¹⁷⁶Lu/¹⁷⁷Hf ratios determine the path of Hf isotope evolution lines for both the host rocks and their sources, and the calculated zircon Hf model ages are projected back to their sources using evolution curves based on the ¹⁷⁶Lu/¹⁷⁷Hf ratios assigned to the host rocks. As shown in Fig. 2, the Lu/Hf ratios are a key factor in determining the Hf isotope evolution lines for new continental crust (blue dashed line), depleted mantle (red dashed line), and ultra- depleted mantle (solid green line) and the trajectories of the host rock (solid lines with arrows). In addition, it has been well documented that: (1) the chemical composition of Earth's mantle is highly heterogeneous (Hofmann, 1997; Zindler and Hart, 1986); and (2) source- and magma-mixing are commonly involved in formation of continental crust (Li et al., 2009b; Kemp et al., 2007). Thus, unless the Lu/Hf ratios of both the host rocks and their mantle sources can be determined for each individual zircon grain, the zircon model ages are meaningless, regardless of the specific geochemical model involved in the Hf model age calculations. As shown in Fig. 2, adjusting the Lu/Hf ratios of the melt/source will result in large changes in the Hf model ages. In the case of detrital zircon, each individual grain could potentially have crystallised from a different igneous rock. Hence, it is impossible to constrain the Lu/Hf ratio of the host rock from which an individual detrital zircon was derived.

Following on from the above, model ages are always going to be suspect because the source rocks of individual detrital zircons are unknown. A model age is valid only if the host rock of a zircon was derived from a single source (e.g., Arndt and Goldstein, 1987). However, magma- or source-mixing is inevitable in the formation of crustal rocks. Experimental data (Castro et al., 2010), partial melting models (Niu et al., 2013), and geochemical evidence (Tang et al., 2012) show that partial melting of a mixture of subducted oceanic basalts and sediments is an important way of producing juvenile continental crust. Furthermore, we have demonstrated above that zircon oxygen isotope signatures cannot alone distinguish new crustal additions from reworking of pre-existing continental crust. As illustrated in Fig. 2, igneous rocks from the Vijiayan Complex show strong mantle-like Hf isotopic signatures, suggesting that juvenile crustal addition occurred around 1.1–1.0 Ga. However, the calculated zircon Hf model ages based on depleted mantle (T_DM) or arc mantle (T_NCC) vary from about 1.2 to 2.0 Ga. If detrital zircon grains were later derived from this complex, the resultant zircon Hf model ages would not constrain the timing of the crustal growth event. In contrast, U-Pb crystallisation ages of zircon grains that plot within the juvenile crustal field (Fig. 2) can preserve a record of such a crustal growth event. It is therefore inappropriate to use their Hf model ages, which may have no real geological meaning, for estimating the time of continental crustal growth. The use of detrital zircon Hf model ages is therefore not as straightforward as some studies would suggest (e.g., Dhuime et al., 2012).

Finally, if there is a large range of Hf isotope compositions in the source region from which juvenile continental crust was derived, this will also result in large uncertainties in the calculation of zircon Hf model ages. We have demonstrated in the previous section that juvenile continental crust can be derived from a range of sources, including recycled igneous oceanic crust and basalt derived from other mantle sources, including metasomatised sub-arc mantle, typical depleted mantle, or less depleted mantle (lower mantle reservoirs). Furthermore, each individual reservoir can potentially have a large range of Hf isotope compositions. For example, sub-arc mantle has ¹⁷⁶Hf/¹¹⁷Hf isotope ratios varying from 0.283 271 (ε_Hf= +17.2) to 0.282 913 (ε_Hf=+4.5) (Iizuka et al., 2013), whereas MORB source mantle ranges from 0.292 874 (ε_Hf=+3.5) to 0.283 454 (ε_Hf= +24.1) (Meyzen et al., 2007). Because juvenile continental crust will have a Hf isotopic component that reflects the composition of the parental basaltic melt derived from such reservoirs, it should also possess a large range of Hf isotope values. Thus, the application of a uniform value, such as an average for depleted mantle (red dashed line in Fig. 2) or modern arc basalts (NCC, blue dashed line in Fig. 2), to calculate zircon Hf model ages is also highly questionable. As illustrated in Fig. 2, in the case of the 1.1–1.0 Ga Vijiayan Complex, adjusting Hf isotopes for the melt source from the lower limit value of MORB (ε_Hf=+3.5) to the average value of sub-arc mantle (NCC, blue dashed line) or the depleted mantle (red dashed line) will result in large changes in the Hf model ages, as shown earlier. This implies that zircon Hf model ages are also strongly dependant on the nature of the magma reservoir.

In summary, without information about the source reservoir, extent of partial melting and magmatic evolution of an igneous rock from which zircon crystallises, the zircon Hf model age is not only difficult to calculate but its meaning is hard to evaluate and thus highly questionable. Indeed, it would not be possible to acquire all the information required for a correct interpretation of Hf model ages for detrital zircon grains. Thus, we suggest that detrital zircon Hf model ages on their own are inappropriate for constraining such a critical issue as the growth of continental crust (see Dhuime et al., 2012 for an alternative view). In contrast, regardless of whether the host rock involved magma mixing or was derived from a chemically heterogeneous source, U-Pb crystallisation ages of juvenile zircon grains can record the timing of addition of juvenile material to the crust.

Because O isotopes are sensitive to the contribution of recycled hydrated components in the source region of igneous rocks and Hf isotopes can potentially determine the timing of melt extraction from the mantle, by combining igneous zircon Hf and O isotopes we can investigate the potential conversion of mafic components into andesitic felsic crust. Intense hydrothermal alteration can significantly affect the O isotopic composition of upper oceanic crust, resulting in supracrustal O isotope signatures, whereas water-rock interaction has little or no effect on the Hf isotopic composition. Thus, the most prominent characteristic of upper oceanic crust is the decoupling of hafnium and oxygen isotopes, resulting in depleted mantle-like hafnium isotopes and supracrustal oxygen isotopes. Because partial melting of subducted oceanic crust played an important role in the generation of juvenile continental crust, the co-existence of depleted mantle-like Hf isotopes and supracrustal oxygen isotopes in a single zircon grain is likely to be a widespread phenomenon.

A global compilation of ¹⁷⁶Hf/¹⁷⁷Hf data from both detrital and magmatic zircon grains with mantle-like and supracrustal oxygen isotopes are plotted as a function of their crystallization ages (Fig. 3a). Zircon grains with δ¹⁸O varying from +4.5‰ to +6.5‰ display large ranges in initial Hf isotopes at any given time in geological history (yellow dots for detrital and green triangles for magmatic zircon grains in Fig. 3a) with a marked increase in the spread after 3.0 Ga. As shown in Fig. 3b, zircon grains with Hf isotope compositions plotting within the field of juvenile continental crust, as defined in Fig. 2, display a large range in O isotopic ratios, varying from about +2‰ to > +10‰ with a large proportion of these grains (nearly half of the total analyses) displaying a non-mantle-like oxygen isotopic signature (δ¹⁸O > +6.5‰ or < +4.5‰). The mantle-like Hf isotopic compositions imply that their host rocks were derived directly from the mantle or by remelting of basaltic rocks shortly after extraction from the mantle. Conversely, the non-mantle-like oxygen isotopic signatures suggest that the source of their host rocks was exposed to surface conditions and underwent hydrothermal alteration, or else entrained recycled hydrothermally- altered material (e.g., Bindeman et al., 2005; Eiler, 2001; Eiler et al., 2000, 1998). Therefore, the co-existence of mantle- derived radiogenic and non-mantle-like O isotopic signatures in juvenile rocks is likely attributable to partial melting of hydrothermally-altered juvenile crustal materials (subducted oceanic basalts, plus variable amounts of recycled sediments) or of the metasomatised mantle wedge.

Figure 3.  (a) Plot of ε_Hf(t) versus U-Pb age for magmatic and detrital zircon grains with different δ¹⁸O values. (b) Global compilation of δ¹⁸O versus U-Pb age for detrital zircon grains with Δε_Hf(t)≥0. Data sources: detrital zircon from sediments (Hervé et al., 2013; Kirkland et al., 2013; Dhuime et al., 2012; Li et al., 2012; Wang X-C et al., 2012, 2011; Yin et al., 2012; Lancaster et al., 2011; Rapela et al., 2011; Pietranik et al., 2008; Kemp et al., 2006); modern river sands (Iizuka et al., 2013; Wang C Y et al., 2011, 2009); and magmatic zircon grains (Chen et al., 2013; Li X H et al., 2013, 2010, 2009b; Lu et al., 2013; Pietranik et al., 2013; Wang F Y et al., 2013; Wang M X et al., 2013; Wang X-C et al., 2013b; Dan et al., 2012; Guo et al., 2012; Heinonen et al., 2012; Jiang et al., 2012; Liu X et al., 2012; Muñoz et al., 2012; Szilas et al., 2012; Zheng Y C et al., 2012; Dai et al., 2011; Su et al., 2011; Zhu et al., 2011; Marschall et al., 2010; Sun et al., 2010; van Dongen et al., 2010; Liu D Y et al., 2009; Bolhar et al., 2008; Zheng Y F et al., 2007). Only zircon grains with concordant U-Pb ages (percentage discordance D ranging from +10 to -10) were used for constraining the zircon Hf isotope compositions. Percentage discordance D is defined as if (e^235×T7/6–1) > ²⁰⁷Pb/²³⁵U_m, D=((e^235×T7/6–1)–²⁰⁷Pb/²³⁵U_m)²+((e^238×T7/6–1)–²⁰⁶Pb/²³⁸U_m)^0.5/((²⁰⁷Pb/²³⁵U_m)²+ (²⁰⁶Pb/²³⁸U_m)²)^0.5, or D= -((e^235×T7/6–1)–²⁰⁷Pb/²³⁵U_m)²+((e^238×T7/6–1)–²⁰⁶Pb/²³⁸U_m)^0.5/((²⁰⁷Pb/²³⁵U_m)²+(²⁰⁶Pb/²³⁸U_m)²)^0.5, where m and T7/6 indicate measured values and ²⁰⁷Pb/²⁰⁶Pb age, respectively. The decay constants used were: λ²³⁵=9.848 50E–10, λ²³⁸=1.551 25E–10. The grey fields in (a) and (b) indicate the Hf isotope composition of juvenile crust (as defined in Fig. 2) and mantle-like zircon δ¹⁸O values, respectively.

Although partial melting of metasomatised mantle wedge by fluids released from subducted slabs can generate depleted mantle-like Hf isotopes with non-mantle-like oxygen isotope signatures, the resultant melts are dominantly mafic (Auer et al., 2009; Dorendorf et al., 2000) rather than the more silica-rich andesites that characterise the bulk continental crust. This process alone therefore cannot explain the compositional paradox that juvenile continental crust displays depleted mantle-like radiogenic isotopes yet has a non-mafic major element composition. In contrast, partial melting of juvenile crustal materials (oceanic crust or juvenile mafic arc continental crust) can explain the reason why there is a distinction between depleted mantle-like radiogenic isotopes and the non-mantle signature of the major elements. Numerical simulations show that melts/fluids released from subducted oceanic basalts and sediments are controlling factors in crustal growth at continental margins (e.g., Zhu et al., 2013; Vogt et al., 2012). Experimental data suggest that melting of subducted slabs and diapirism may be the rule rather than the exception beneath arc fronts in most subduction zones (Behn et al., 2011; Plank et al., 2009; Klimm et al., 2008; Hermann and Spandler, 2007), and can produce melts that strongly resemble andesitic arc magmas (Castro et al., 2010). This suggests that andesitic continental crust can be constructed directly from mantle- modified slab melts (Castro et al., 2013, 2010; Behn et al., 2011; Gómez-Tuena et al., 2008), or partial melting of mantle metasomatized by slab-derived melts/fluids (Chen et al., 2014). This is consistent with the previous conclusions that non-mantle δ¹⁸O values in primary arc magmas are ultimately attributable to fluids/melts from subducted oceanic basalts and sediments (Liu et al., 2014; Auer et al., 2009; Dorendorf et al., 2000; Eiler et al., 2000, 1998). Thus, the decoupling of hafnium isotopes from oxygen isotopes in single zircon grains provides convincing evidence that partial melting of subducted oceanic basalts and sediments plays a key role in the generation of juvenile continental crust.

As is evident from the above discussion, both Hf isotope systematics and U-Pb crystallisation age of magmatic zircon are crucial for constraining the timing and rate of crustal growth. Based on this, we define a new approach whereby whole-rock Hf isotope data (not zircon data) obtained from depleted and less depleted mantle-derived rocks is integrated with data from juvenile crustal igneous rocks, and used to determine the Hf compositions of juvenile continental crust throughout Earth's history (Fig. 2). The newly-defined range of Hf isotopes for juvenile continental crust establishes a new criterion for the recognition of zircon grains that were derived from juvenile or dominately juvenile crustal rocks.

Based on the above evaluation, U-Pb crystallisation ages of juvenile zircon grains can now be used to constrain patterns of crustal growth, to test whether or not the growth of continental crust was episodic, and to unravel the possible factors that control/affect the growth and evolution of continental crust.

Using the Hf isotopic range of juvenile continental crust defined by this study in Fig. 2, we investigate the potential secular changes in the formation of new continental crust. As shown in Fig. 4a, the proportion of zircon grains from juvenile continental crust, when compared to the total of dated zircon at any given time, drops dramatically from ~72% at ca. 3.5 Ga to about 30% at ca. 3.4 Ga. From 3.4 Ga to the present day, the proportion periodically fluctuates, with four cycles recognized (3.4–2.1, 2.1–1.3, 1.3–0.7, and 0.7–0 Ga). This cyclical nature and its timing is almost identical to the cycle of assembly and breakup of supercontinents (Li and Zhong, 2009). Within each cycle, the proportion of juvenile continental crust formation ages decreases during the assembly phase of the supercontinent and then increases toward its breakup (see 'supercontinent cycles' at the top of Fig. 4a). Overall, the currently-available data imply that the most dramatic shift in the proportion of juvenile continental crust formation ages occurred at ca. 3.5 Ga. This is different from the 'U'-shaped curve of the proportion of juvenile continental crust formation ages through Earth's history (dashed black line in Fig. 4a) based on zircon oxygen isotope data and Hf model ages (Dhuime et al., 2012). Dhuime et al. (2012) implied a dramatic change in the proportion at ca. 3.0 Ga, which they interpreted as the start of subduction-driven plate tectonics. However, our results indicate that this dramatic change occurred ca. 500 Ma earlier at ca. 3.5 Ga (Fig. 4a). Because the number of samples older than 3.4 Ga is small, we cannot entirely rule out sampling bias, but the available data clearly indicate a sharp change at this time.

Figure 4.  (a) Juvenile continental crust formation ages and crustal reworking ages throughout Earth's history. Zircon grains crystallized in juvenile magmas have Δε_Hf(t)≥0, whereas those formed during reworking of crustal materials have Δε_Hf(t) < 0. The data points represent the proportion of U-Pb ages associated with new crust generation in the worldwide zircon record, calculated using a 100 Ma sliding age window with 100 Ma steps, only when there are n≥2 analyses for each time window. P. Pangea; G. Gondwana; R. Rodinia; C/N. Nuna/Columbia; K. Kenorland. Data sources are the same as in Fig. 3. Age data for assembly of supercontinents: Pangea-Gondwana 0.60–0.32 Ga, Rodinia 1.3–0.9 Ga, Nuna/Columbia 2.1–1.8/1.6 Ga; and for rifting-breakup: Pangea 0.25–0.90 Ga; Rodinia 0.85–0.60 Ga; Columbia/Nuna 1.4–1.2 Ga. The orange color bands highlight the lifespan of stable supercontinents. Age data for supercontinent cycles are from Li et al. (2008), Zhang et al. (2012), and Stampfli et al. (2013). The possible pre-Nuna supercontinent (Kenorland) is poorly constrained. (b) Continental growth curves. The cumulative volume of crust is calculated for 100-million-year time intervals. Curve 1 (green) is based on zircon U-Pb ages from grains with Δε_Hf(t)≥0. Curve 2 (red) is based on zircon U-Pb ages with Δε_Hf(t)≥0 (curve 1) and the Global Lithospheric Architecture Mapping (GLAM) model (Belousova et al., 2010). Curve 3 is the continental crustal growth curve from geochemical modeling of REEs in conjunction with isotopic data (Taylor and McLennan, 1995). Curve 4 is based on zircon Hf model ages with mantle-like (MORB) oxygen isotopes (Dhuime et al., 2012). Curve 5 (grey shaded zone) is based on the evolution of the atmospheric ⁴⁰Ar/³⁶Ar ratio as a function of time (Pujol et al., 2013). Curve 6 is based on the mass balance between generation and recycling of continental crust (Armstrong and Harmon, 1981). Curve 7 is based on age distribution peaks in orogenic granitoids (Condie and Aster, 2010). Curve 8 is based on the geographic distribution of Rb-Sr and K-Ar isotope ages (Hurley and Rand, 1969). Curve 9 is based on zircon Hf model ages and the GLAM model (Belousova et al., 2010). Inset in (b) shows the details of curve 2 in the post-Archean.

Figure 4b shows various estimates of continental growth through time. The crustal growth curve of Taylor and McLennan (1995) (curve 3 in Fig. 4b) was constrained by geochemical modeling of REEs in conjunction with isotopic data, based on the island-arc model; this still remains the most popular model. Growth curves 4 and 9 are based on age distribution peaks in detrital zircon, whereas curve 7 is based on orogenic granitoids. Growth curve 5 (grey band in Fig. 4b) is based on the argon isotopic composition of the ancient atmosphere (Pujol et al., 2013), whereas growth curves 6 (Armstrong and Harmon, 1981) and 8 (Hurley and Rand, 1969) are representative of two contrasting end-member estimates, based on early data and assumptions that are not currently held. Because detrital zircon Hf model ages are inappropriate for constraining the growth of continental crust, growth curves 4 and 9 are highly questionable. Bradley (2011) pointed out that the approach for estimating growth curve 7 of Condie and Aster (2010) has a number of problems: (1) global data coverage is uneven, being less in remote areas; (2) countless plutons have been lost from the geologic record, either buried beneath sedimentary cover or totally destroyed; (3) classification of some granites as "orogenic" is equivocal; and (4) ⁴⁰Ar/³⁹Ar ages were not included, resulting in an underestimate of Mesozoic and Cenozoic igneous rocks. Although growth curve 3 is widely-accepted, the island-arc model may not account for crustal addition in terms of mass balance (Niu et al., 2013). More importantly, with the exception of growth curve 3, all the other curves do not consider the generation mechanism and isotopic charactersitcs of juvenile continental crust. Thus, new estimates incoporating the continuously-updated global database and defining the characterstics of juvenile crust are required.

Our initial crustal growth curve 1 (green curve in Fig. 4b) is based on the newly-defined whole-rock Hf isotopic characteristics of juvenile continental crust, as outlined above. Our crustal growth curve 2 combines growth curve 1 with the Global Lithospheric Architecture Mapping (GLAM) model of Belousova et al. (2010). The GLAM model integrates geophysical, geological, geochronological, and geochemical data for the crust and lithospheric mantle to produce maps of lithospheric composition and architecture. According to the GLAM model, 70.5% of existing continental upper lithosphere was formed in the Archean (Archon), 19.3% in the Proterozoic (Proton) and 10.2% in the Phanerozoic (Tecton). For this study, zircon U-Pb ages from juvenile continental crust were classified into the same three domains: Archon- origin (> 2.5 Ga, 8.8%), Proton-origin (2.5–1 Ga, 46.7%) and Tecton-origin (< 1 Ga, 44.5%). The factors, (8.01 for Archon- origin, 0.41 for Proton-origin and 0.23 for Tecton-origin), as defined by ratios between the GLAM estimates and the proportions of zircon data for different time intervals, were applied to calculate the cumulative growth curve 2 in Fig. 4b (see details in Belousova et al., 2010). Our favoured crustal growth model (curve 2) shows that continental crust underwent two stages of rapid growth (3.7–3.5 and 3.0–2.5 Ga), with long-term slow progressive growth (post-Archean), punctuated by episodes of more rapid crust formation. Our estimate suggests that about 70% of the continental crust had been generated by 2.5 Ga. This is consistent with the constraints from Nd isotope data for Australian shales (e.g., Cawood et al., 2013): it also includes the major period of formation of Archean komatiites. For example, some 70% of all known komatiites were erupted in the interval 2.73–2.70 Ga (e.g., Arndt and Davaille, 2013). Overall, our estimate, though still based on a moderate number of juvenile zircon datasets, is very similar to the Taylor and McLennan (1995) crustal growth (curve 3 in Fig. 4b), but differs significantly from the Belousova et al. (2010) model (curve 9 in Fig. 4b) based solely on zircon Hf model ages.

The currently-accepted Earth geodynamic hypothesis postulates that mantle upwellings (plume) and plate tectonics- driven deep subduction control global material circulation and energy exchange between Earth's lithosphere and its interior. These are two of the first-order geotectonic phenomena that potentially operated throughout much of Earth history (e.g., Li and Zhong, 2009; Zindler and Hart, 1986; Hofmann and White, 1982). Multiple subduction zones promote the rapid amalgamation of continental fragments into supercontinents, and also act as major zones of material circulation in the Earth system (Santosh, 2010). The subducted material accumulates at the mantle transition zone and either sinks down to the core-mantle boundary or is returned back to the upper mantle, where it triggers/controls the generation of plumes and superplumes that ultimately fragment the supercontinents (Yoshida and Santosh, 2011; Li and Zhong, 2009). Geological evidence related to the fragmentation and dispersal of the three major supercontinents that are generally accepted (Nuna/Columbia, Rodinia, and Gondwana) attests to the involvement of plumes as well as plate tectonics (Yoshida and Santosh, 2011; Li et al., 2008). Thus, although it is still debated, whole-of-Earth geodynamic theory suggests that global-scale material circulation between Earth's crust and its interior is most likely controlled by the first-order cycle of plumes, plate tectonics and their consequence-supercontinents and superplumes. The growth of continental crust, as a part of this global-scale material circulation, is therefore also controlled by these first-order geodynamic cycles, and we examine this in detail below.

A compilation of published global U-Pb data comprising 220 000 zircon grains, plus 7 867 magmatic orogenic granitoid rocks reveal seven prominent and statistically significant age populations: 3.5–3.3, 2.9–2.4, 2.2–1.7, 1.3–1.0, 0.9–0.6, 0.6–0.4, and 0.3–0.1 Ga (Fig. 5a). Similar zircon crystallisation age peaks have been documented on all seven continents (Voice et al., 2011). The peaks at 2.2–1.7, 1.3–1.0, and 0.6–0.4 Ga correlate well with the assembly of the supercontinents Gondwana-Pangea (0.60–0.32 Ga), Rodinia (1.3–0.9 Ga), and Nuna/Columbia (2.1–1.8/1.6 Ga). The peak at 0.9–0.6 Ga is correlated with the rifting and breakup of Rodinia (0.85–0.60 Ga), whereas the peak at 0.3–0.1 Ga coincides with final assembly (ca. 0.32 Ga) and rifting and breakup of Pangea (0.25–0.09 Ga). This has previously been linked to episodic generation of continental crust (Voice et al., 2011; Condie and Aster, 2010; Condie et al., 2009; Li and Zhong, 2009; Campbell and Allen, 2008; Condie, 2000, 1998). Zircon U-Pb ages in grains with Δε_Hf(t)≥0 reveal six statistically significant peaks: 2.8–2.4, 2.2–1.8, 1.7–1.5, 1.3–1.0, 0.9–0.6, and 0.4–0.1 Ga (Fig. 5b). The peak at 0.9–0.6 Ga also correlates well with the breakup of Rodinia (0.85–0.60 Ga). The peak at 1.3–1.0 Ga spans the late-stage breakup of Nuna/ Columbia (1.4–1.2 Ga) and the early-stage assembly of Rodinia, whereas the peak at 0.4–0.1 Ga covers the late-stage assembly of Pangea (final assembly at ca. 0.32 Ga) to its rifting and breakup (0.25–0.09 Ga) (Fig. 5b). The zircon U-Pb ages of grains with Δε_Hf(t) < 0 define four statistically significant peaks: 2.8–2.4, 2.1–1.5, 1.2–0.9, and 0.7–0.4 Ga (Fig. 5c), and also broadly correlate with the timing of supercontinent assembly.

Figure 5.  Episodic evolution and differentiation of continental crust. (a) Histogram of crystallization age defined from a global compilation of detrital zircon U-Pb ages with concordance > 85% (data source is the same as that in Fig. 3) and orogenic granitoid ages (Condie and Aster, 2010). (b) Histogram and Gaussian distribution (red curve) of juvenile crust formation ages defined by zircon grains with Δε_Hf(t)≥0. (c) Histogram and Gaussian distribution (red curve) of crustal reworking ages defined by zircon grains with Δε_Hf(t) < 0. (d) Intensity of LIP events (Prokoph et al., 2004). The Gaussian distribution for each event is based on age and uncertainty (applying a uniform uncertainty of 1.5% SD to all databases). The best-fit sine wavelengths are (a) 720 Ma for the total zircon database; (b) 650 Ma for the juvenile crustal age; and (c) 750 Ma for the reworking crustal age. Dashed blue curve in (d) is drawn from Li and Zhong (2009) according to the cyclical nature of LIP age records with the ca. 750–550 Ma wavelengths obtained by Prokoph et al. (2004) using an integrated wavelet, a spectral, and a cross-spectral approach. It is noted that peaks at 0.9–0.7 and 0.3–0.1 Ga defined by the total zircon database in (a) and at 1.7–1.5 Ga defined by juvenile zircon grains in (b) do not match the cyclicity of the whole database; see text for discussion. The algorithms for wavelet analysis in (a)–(c) are the same as in (d) by Prokoph et al. (2004).

Knowledge of the history of melt extraction from the mantle is essential for an understanding of continental dynamics and the long-term stability of ancient continental landmasses (Griffin et al., 2004). Such information can be gleaned by examining large igneous provinces (LIPs, volcanic provinces characterized by anomalously high rates of mantle melting that represent the largest volcanic events in Earth history). The formation of LIPs over the past 3.5 Ga cluster in five peaks at 2.8–2.4, 2.2–1.8, 1.3–1.1, 0.9–0.6, and 0.3–0 Ga (Fig. 5d), also correlating well with the assembly and breakup of supercontinents (Doucet et al., 2020; Gamal El Dien et al., 2019). The history of melt extraction from the mantle can also be revealed by examining the Re-Os isotopic system (Shirey and Walker, 1998; Walker et al., 1989). Although osmium model ages are model-dependent with large uncertainties, they are widely considered the best parameter to constrain the history of melt extraction from the mantle (Pearson et al., 2007; Brandon et al., 2006; Aulbach et al., 2004; Griffin et al., 2004, 2002; Alard et al., 2002). For instance, recent Os isotope results for PGE-alloy grains from ophiolites of various ages suggest that the Os isotope composition of the mantle records "pulses" of mantle depletion attributable to discrete continental crust extraction events (Pearson et al., 2007). A compilation of published Re-Os isotope data from osmium-rich platinum-group alloy grains from mantle xenoliths and ophiolites (n=1 290) shows that Os model ages (T_RD) also cluster at four distinct periods in Earth's history: 2.9–2.6, 1.6–1.1, 0.9–0.6, and 0.3–0.2 Ga (Fig. 6a) and Os T_RD from whole-rock samples cluster at four similar periods: 2.9–2.4, 1.6–1.1, 0.9–0.6, and 0.3–0 Ga (Fig. 6b). Osmium model age peaks at 2.9–2.4, 0.9–0.6, and 0.3–0 Ga (Figs. 6a and 6b) closely match juvenile zircon U-Pb age peaks at 2.8–2.4, 0.9–0.6, and 0.4–0.1 Ga (Fig. 5b) and LIP formation age peaks at 2.8–2.4, 0.9–0.6, and 0.3–0 Ga (Fig. 5d). Considering the large uncertainty associated with Os model ages, the peak at 1.6–1.1 Ga can be viewed as comparable with the juvenile zircon U-Pb age peak at 1.3–1.0 Ga. The coincidence of the age peaks between LIPs, mantle melting events and juvenile continental crust suggests that there is a common driver controlling the episodic nature of large-scale melt extraction from the mantle.

Figure 6.  (a) Mantle Os model ages (T_RD, Ga) defined by osmium-rich platinum- group alloy grains (Shi et al., 2012, 2007; Wang et al., 2010; Xu et al., 2008; Pearson et al., 2007; Brandon et al., 2006; Aulbach et al., 2004; Griffin et al., 2004, 2002; Alard et al., 2002). The mantle Os model ages were calculated using the Enstatite Chondritic Reservoir (ECR) that has a present-day ¹⁸⁷Os/¹⁸⁸Os value of 0.128 1±0.000 4 and ¹⁸⁷Re/¹⁸⁸Os=0.421±0.013 (Walker et al., 2002). (b) Mantle Os model ages defined by peridotite (Wittig et al., 2010; Chu et al., 2009; Schilling et al., 2008; Xu et al., 2008; Handler et al., 2003). (c) Correspondence of OIB ⁴He/³He peaks with juvenile continental crustal zircon age peaks. A probability distribution function of OIB (Parman, 2007) is shown along the right-hand side of Fig. 6c. The density plot of juvenile continental crust formation ages is from Fig. 5b. Six juvenile zircon age peaks (2.5, 1.6, 1.2, 1.1, 0.75, and 0.3 Ga) correlate with ⁴He/³He peaks in the OIB distributions (connected by solid blue lines). A regression (red dashed line) through these points (green filled circles) yields ⁴He/³He (10³)=88.89e^-0.593t, R²=0.993 and t in Ga. Such a regression matches well with the He isotope evolution for depleted mantle that is a residue of extraction of continental crust (Lee et al., 2010). The dominant ⁴He/³He peak in the mid-ocean-ridge basalts (91.0±1.5)×10³ (Parman, 2007) and the initial ⁴He/³He of the Earth (6.0±0.2)×10³ (Mahaffy et al., 1998) are highlighted by filled purple circles.

The Morlet Wavelet transform is a powerful tool in time-frequency analysis (Torrence and Compo, 1998). It is applied to delineate temporal variations of cycle amplitude and phase over a 0–3 500 Ma spectrum for all zircon datasets. The Morlet wavelet analyses show that distribution patterns of zircon ages have 650–750 Ma best-fit sine wavelets (Figs. 5a–5c).

As shown in Fig. 5a, the Morlet Wavelet transform defines the cyclicity of five age peaks at 3.5–3.3, 2.9–2.4, 2.1–1.7, 1.3–1.0, and 0.6–0.4 Ga, but two age peaks at 0.9−0.6 and 0.3−0.1 Ga do not fit the pattern. Supercontinents can form in two tectonic scenarios: extroversion and introversion (Murphy and Nance, 2003). Extroversion-type supercontinents are characterized by juvenile crustal growth, whereas introversion-type supercontinents are characterized by crustal reworking (Murphy and Nance, 2003). The Pangea supercontinent is an extroversion- type with an increase of juvenile Nd isotopes in crustal rocks and so the age peak at 0.3–0.1 Ga is consistent with constraints from crustal Nd isotopes (Murphy and Nance, 2003). This implies that whether a supercontinent is of extroversion or introversion type has implications for continental crustal evolution. The age peak at 0.9–0.6 Ga matches well with the widespread 850–630 Ma magmatic events in South China (Zhao and Asimow, 2018; Lyu et al., 2017; Yang et al., 2016; Wang et al., 2009, 2008, 2007), which were interpreted to record the initiation of the breakup of Rodinia (Li et al., 2008).

When zircon grains are classified into juvenile and reworked types by their Δε_Hf values, the Morlet wavelet analyses generate 650 (Fig. 5b) and 750 Ma (Fig. 5c) best-fit sine wavelets, respectively. The whole juvenile zircon data set is cyclical and defines five age peaks at 2.8–2.4, 2.2–1.8, 1.3–1.0, 0.9–0.6, and 0.4–0.1 Ga (Fig. 5b). This distribution pattern is similar to the 550–750 Ma periodicity of LIP events through time (Fig. 5d; Prokoph et al., 2004), implying that heat, and possibly materials, derived from uprising hot mantle (plume or superplume) is important for continental crustal growth. It should be pointed out that the age peak at 1.7–1.5 Ga defined by juvenile zircon grains does not match the peaks defined by the full data set (Fig. 5b) and the LIP record (Fig. 5d).

The intensified crustal reworking events at 0.6–0.5 Ga (Fig. 5c) strongly affected the wavelet analysis and its intensity was artificially reduced by about 40% to get a best-fit sine wavelength for the whole dataset. The full data set defines four age peaks at 2.8–2.4, 2.1–1.5, 1.2–0.9, and 0.7–0.4 Ga (Fig. 5c).

Taken together, zircon age distribution patterns have 650–750 Ma cycles (Figs. 5a–5c), which is similar to the duration of supercontinent-superplume cycles as proposed by Li and Zhong (2009). Apparent frequency peaks in the ³He/⁴He global OIB dataset may also record global-scale mantle depletion events (Parman, 2007; Parman et al., 2005). Thus, the correlation between peaks in the distribution of ⁴He/³He in ocean island basalt (OIB) and peaks in crustal zircon ages was attributed to coupling of global-scale mantle depletion and juvenile continental crust formation (Parman, 2007). However, Rudge (2008) proposed that the correlation of Parman (2007) was likely to be statistically insignificant. Furthermore, the Parman (2007) correlation generated a linear trend for helium evolution for the convective mantle (purple dashed line in Fig. 6c). Such a linear evolution trend is inconsistent with estimates for the convective mantle (thick grey line in Fig. 6c). We re-investigated the correlation of age peaks of juvenile continental crust formation and peaks in the OIB ⁴He/³He distribution and found that the helium isotopes of the convective mantle evolved along a curve of ⁴He/³He (10³)=88.89e^-0.593T(Ga). Such a helium isotope evolution trend is similar to the evolution curve constrained by top-down differentiation and mass flow between mantle and crust (grey thick line in Fig. 6c; Lee et al., 2010). Thus, global mantle depletion events do, indeed, correspond well with the age peaks of juvenile continental crust formation.

The episodic nature of mantle melting, as evidenced by LIPs and mantle Re-Os and ⁴He/³He isotope records, most likely reflects the thermochemical state of the mantle and its convection style. The association between the timing of komatiite eruption and granite ages provides a direct link between plume activity and the U-Pb age peaks (Arndt and Davaille, 2013 and references therein). The ages of Late Archean and Early Proterozoic greenstone belts also coincide with these peaks (Arndt and Davaille, 2013 and references therein). As noted above, Arndt and Davaille (2013) suggested that the growth age peaks of continental crust in the Archean and Proterozoic record periods of enhanced mantle convection, leading to more rapid subduction and consequently to periods of accelerated growth of continental crust. They considered plumes in the upper mantle to be the driving force of episodic mantle convection and regarded enhanced subduction as a consequence of the arrival of a hot plume in the upper mantle.

LIPs are one of the most prominent surface expressions of mantle melting events. A whole variety of models have been proposed to explain the origin of LIPs (see summaries in Bryan and Ernst, 2008; Ernst et al., 2005; Saunders, 2005), including mantle plumes emanating from the core-mantle boundary (e.g., Campbell, 2007; Coffin and Eldholm, 1992; Campbell and Griffiths, 1990; Richards et al., 1989), impact-induced decompression melting (e.g., Ingle and Coffin, 2004; Jones et al., 2002), lithospheric delamination (Hales et al., 2005; Tanton and Hager, 2000); decompression melting during rifting (White and McKenzie, 1989) or following mantle heating beneath supercontinents (Coltice et al., 2009, 2007); edge driven convection (King and Anderson, 1998); melting of fertile mantle without excess heat (Anderson, 2005) or shallow melting anomalies generated by plate tectonic-related processes (Foulger, 2007); horizontally extended slabs triggering partial melting (Cox, 1978), stress- induced lithospheric fracturing and drainage of a relatively slowly-accumulating sub-lithospheric basaltic magma reservoir (Silver et al., 2006); and back-arc rifting (Smith, 1992). The debate of a plume versus a non-plume origin for LIPs is less relevant here, but it is necessary to discuss it briefly.

Numerical modeling supports the hypothesis that mantle plumes and plate tectonics are both essential parts of the Earth's self-organization geodynamic system (see Wang X-C et al., 2013a; Li and Zhong, 2009 for reviews). Recent studies show that an ancient mantle reservoir, generally considered to be a dense chemical layer located at the core-mantle boundary, played an important role in the generation of LIPs (Wang X-C et al., 2013a, b; Campbell and O'Neill, 2012; Caro, 2011; Jackson and Carlson, 2011; Nomura et al., 2011; Jackson et al., 2010). Furthermore, such ancient reservoirs may also have co-existed with young recycled components in the mantle source of Late Cenozoic basalts (Wang X-C et al., 2013a). Thus, mantle plumes are intimately linked to plate tectonics (Wang X-C et al., 2013a, b; Steinberger and Torsvik, 2012; Li and Zhong, 2009) and both contribute to the formation of LIPs.

The co-existence of OIB-like and arc-like basalts within a single intra-continental basaltic magma province (Wang X-C et al., 2014, 2009, 2008; Ivanov and Litasov, 2014; Jourdan et al., 2007; Puffer, 2001) is generally attributed to plume-lithosphere interaction (Wang X-C et al., 2014, 2009, 2008; Puffer, 2001; Hawkesworth et al., 2000, 1995; Turner and Hawkesworth, 1995). However, geochemical evidence from LIPs formed during breakup of the Pangea supercontinent (Wang et al., 2016, 2015; Heinonen et al., 2014, 2013, 2010; Ivanov and Litasov, 2014; Merle et al., 2013; Rooney et al., 2012; Sobolev et al., 2011; Jourdan et al., 2007; Ewart et al., 2004) show that such a feature may also be explained by plate tectonic-driven deep-Earth volatile cycling, as proposed by Cox (1978), Ivanov and Litasov (2014) and Wang et al.(2016, 2015). The role of mantle plumes and deep-volatile cycling in crustal growth is briefly outlined below.

The following lines of evidence imply that the integrated effect of ponding hot material beneath the transition zone, and consequent wet upwelling, was responsible for destabilisation of stagnant slabs and that it enhances mantle convection and thereby exerts a control on continental crustal growth. Firstly, the initial rifting of supercontinents seems to be earlier than the onset of plume-induced magmatic events (Wang et al., 2016). This is at odds with the classical plume hypothesis that proposes that mantle plumes trigger the rifting and breakup of supercontinents (Lyu et al., 2017). Slab avalanches within the hydrated transition zone would generate large-scale wet upwellings that hydrate the sub-continental lithospheric mantle, resulting in instability of the lowest part of the continental lithosphere (Peslier et al., 2010; Windley et al., 2010). Hydrous basaltic melts would prefer to pond at the boundary between the lithosphere and asthenosphere (Crépisson et al., 2014; Sakamaki et al., 2013), providing a lubricating mechanism for speeding-up plate motion (Lyu et al., 2017; Wang et al., 2016). Slab avalanche-driven upwelling and mobility of hydrous basaltic melts would hence promote the breakup of a supercontinent. It is therefore possible that slab avalanche-driven upwelling triggers the initial weakening of supercontinents. The evidence from the Rodinia supercontinent shows that continental extension, anorogenic magmatism, and rifting occurred ca. 30 to 50 Ma earlier than the first well- documented mantle plume event at 825–800 Ma (Lyu et al., 2017), which included large-scale pre-LIP uplifting, emplacement of high-temperature magmas and 825–810 Ma continental flood basalts (Wang et al., 2009, 2008, 2007; Li et al., 2008). However, it should be noted that the reconstruction model for the Rodinia supercontinent (Li et al., 2008) is still debated, but this is beyond the scope of this study.

Importantly, the breakup and assembly of supercontinents are likely to overlap. For example, Rodinia broke up in stages from ca. 850 to ca. 600 Ma. Before Rodinia had completely broken apart, some of its pieces had already been reassembled into a new configuration (Gondwanaland; Bradley, 2011). The peak of juvenile continental growth at 0.4–0.1 Ga covers both the late-stage assembly of Pangea and its rifting and breakup; likewise, the peak at 0.9–0.6 Ga coincides with the final assembly of Rodinia and its subsequent rifting and breakup (Li X H et al., 2010; Li W X et al., 2010; Li and Zhong, 2009; Wang X-C et al., 2009; Li Z-X et al., 2008) (Fig. 5b). The 1.3–1.0 Ga peak, also coincides with the final breakup of Nuna/Columbia and the early assembly of Rodinia. Therefore, if breakup of supercontinents is caused by superplume activity (Li and Zhong, 2009) and ponding at the base of the mantle transition zone is a primary feature of superplumes (Wang et al., 2016, 2015; Cao et al., 2011; Boschi et al., 2007), we speculate that periods of enhanced mantle convection may be initiated by ponding of a hot mantle plume beneath the mantle transition zone, accompanied by destabilisation of the stagnant slabs residing there (Fukao et al., 2009).

Although the exact timing of global-scale mantle partial melting and enhanced formation of juvenile continental crust is still debated, due to both the small size of the global database and to sampling bias, they do seem to be well correlated when considering the Earth as a self-consistent geodynamic system. Cycles of supercontinent and LIP events (plume/superplume activities) are surface expressions of Earth's self-consistent geodynamic system, and they appear to have controlled the first-order pattern of crustal growth. Although the peak at 1.7–1.5 Ga does not seem to be related to the LIP cycles (Fig. 5b), it does coincide with local rifting events within the Nuna supercontinent that began as early as ca. 1.8–1.6 Ga (e.g., Zhai et al., 2015; Evans and Mitchell, 2011), as evidenced by the Dubawnt granitoids in central Laurentia (Rainbird and Davis, 2007; Rainbird et al., 2006), the coeval Hekla Sund volcanic rocks in northern Greenland (Pedersen et al., 2002), the Ulkan and Urik-Iya grabens in southern Siberia (Pisarevsky et al., 2008), the Cleaver dikes (Irving et al., 2004) plus Bonnet Plume River intrusions (Thorkelson et al., 2001) in Northwest Laurentia, and rift-related mafic dykes in the North China Craton (e.g., Li et al., 2015; Peng, 2015; Zhai et al., 2015). Although these early local intracontinental rifting events predated the arrival of the actual plume that resulted in the breakup of the Nuna/Columbia supercontinent, we propose that wet mantle upwelling driven by deep-volatile cycling (Lyu et al., 2017; Wang et al., 2016, 2015) most likely played a key role in triggering the rifting.

The coincidence between the earliest phases of juvenile continental crust formation at 2.8–2.4 and 2.2–1.8 Ga with the two major episodes of komatiite eruption at 2.73–2.70 Ga (70% of all known komatiites) and at ca. 1.9 Ga (Arndt and Davaille, 2013 and references therein) and the LIP records (Fig. 5d) suggests that mantle thermal events were more important than supercontinent cycles for promoting crustal growth in the early Earth because of higher mantle potential temperatures and possibly more intensive mantle plume events. However, there is a poorly-constrained supercontinent, Kenorland, speculated to have formed at 2.8–2.5 Ga (Yakubchuk, 2019; Strand and Köykkä, 2012; Lubnina and Slabunov, 2011; Aspler and Chiarenzelli, 1998; Williams et al., 1991). The breakup of this supercontinent may have extended from ca. 2.5 to 2.1 Ga (Strand and Köykkä, 2012; Aspler and Chiarenzelli, 1998). The earliest phases of juvenile continental crust formation at 2.8–2.4 Ga may therefore relate to both intensive plume activity and assembly of the Kenorland supercontinent. The other peak at 2.2–1.8 Ga extends from the breakup of a potential Kenorland supercontinent to the assembly of Nuna. It should be pointed out that the assembly-breakup of Kenorland(?) and Nuna/Columbia remains poorly constrained. This will partially affect our examination of the relationship between the two earliest phases of juvenile continental crust formation.

In summary, the cyclical nature of growth and evolution of continental crust, including reworking, reflect the integrated effect of superplume-supercontinent cycles, whole mantle convection style, large-scale mantle partial melting events, and global- scale material circulation and energy exchange between Earth's lithosphere and its interior at each episode of supercontinent- superplume activity, and the driving force is the integrated effect of these fundamental processes. For example, with decreasing mantle potential temperature over time, thermal anomalies (plume events) likely play a much less important role than plate tectonic processes in the growth and evolution of continental crust. This may explain why the early age peaks at 2.8–2.5 and 2.2–1.8 Ga correlate with intensive plume activity (komatiite eruption), whereas the young age peaks correlate with the ultimate expression of plate tectonics, which is the formation and destruction of supercontinents (Fig. 5).

The detrital zircon U-Pb age distribution patterns presented in Figs. 5a–5c reflect the integrating of magmatism and zircon generation rate, but possibly also reflect selective preservation. Zirconium saturation is controlled by the chemical composition of the host rock (Watson and Harrison, 1983; Watson, 1979), and zircon is most common in igneous rocks of intermediate to felsic composition and less common in mafic igneous rocks. Magmas associated with mountain building during supercontinent assembly are dominated by andesitic to more felsic types. Generally, they have high zircon generation rates. In contrast, the epochs of supercontinent breakup are characterized by mantle upwellings (e.g., Arndt and Davaille, 2013; Rolf et al., 2012; Yoshida and Santosh, 2011; Li and Zhong, 2009), resulting in the addition of a large volume of mantle-derived materials into the crust, dominated by basaltic magmatism (exemplified by flood basalts), and rocks from this phase are relatively depleted in zircon.

Alternatively, recent studies have proposed that the U-Pb zircon age peaks do not correspond to pulses of true crustal growth, but to periods of enhanced preservation of continental crust (Cawood et al., 2013; Condie and Kröner, 2013; Roberts, 2012; Condie et al., 2011; Lancaster et al., 2011; Hawkesworth et al., 2010, 2009). The preservation model postulates that both magma volume and zircon generation rates are high in the subduction stage but low in the collision and breakup stages, and that the peaks in zircon U-Pb ages correspond with the time of maximum supercontinent aggregation (Cawood et al., 2013; Hawkesworth et al., 2010, 2009). However, zircon is not ubiquitous in all crustal rocks; as just stated, it is most common in igneous rocks of intermediate to felsic composition and least common in less saturated rocks, such as primary hydrous arc magmas. Experimental studies have demonstrated that the saturation behavior of zirconium in crustal melts is a function of both temperature and melt composition, including water concentration (e.g., Boehnke et al., 2013; Watson and Harrison, 1983; Watson, 1979). One problem with the preservation model is that it does not consider the geochemical behavior of Zr in hydrous arc magmas. Experimental studies on andesitic to basaltic compositions show conclusively that these liquids require an unrealistically high abundance of > 5 000 ppm Zr to directly crystallize zircon, and thus zircons found in mafic environments must have crystallized from late stage evolved melts (Boehnke et al., 2013). The low solubility of Zr in slab-derived fluids (Bernini et al., 2013) would further suppress zircon crystallization from primary hydrous arc magma. Zircon is therefore generally absent in primary arc volcanic rocks, and only gradually starts to appear when these rocks are reworked.

The preservation model is also inconsistent with some other observations. The fundamental assumption of this model is that new continental crust is mainly created by arc magmatism at convergent plate margins. However, several studies indicate that large amounts of continental crust were created in the geological past during short "super events" at rates that are difficult to explain by "typical" subduction activity (Arndt and Davaille, 2013; Niu et al., 2013; O'Neill et al., 2007; Stein and Ben-Avraham, 2007; Condie, 1998; Stein and Goldstein, 1996; Stein and Hofmann, 1994; Schubert and Sandwell, 1989; Reymer and Schubert, 1986, 1984). This implies that hot spot or mantle plume magmatism (e.g., oceanic plateau and basaltic underplating; see summary in Stein and Ben-Avraham, 2007), or syn-collision magmatism (see Niu et al., 2013 for a review) also played an important role in continental crustal growth. For instance, it is considered that the formation of oceanic plateaus and their accretion to continents was an essential part of continental crust production throughout the past 2.7–2.9 Ga of Earth's history (Stein and Ben-Avraham, 2007; Stein and Goldstein, 1996; Stein and Hofmann, 1994; Schubert and Sandwell, 1989). Arndt and Davaille (2013) proposed that mantle plumes may play an important role in the formation of juvenile continental crust through enhancing the rates of subduction. The calculation made by Schubert and Sandwell (1989) suggests that accretion of all oceanic plateaus to the continents on a time scale of 100 Ma would result in a rate of continental growth of 3.7 km³/year. This rate is much higher than the continental growth rate based solely on accretion of island arcs, which is 1.1 km³/year (Reymer and Schubert, 1986). Geochemical and petrological evidence from the CAOB (Tang et al., 2012; Kröner et al., 2007; Jahn, 2004; Jahn et al., 2000; Wu et al., 2000) and the Greater Tibetan Plateau (Niu et al., 2013 and references therein) demonstrate that partial melting of recycled mafic crustal components has also played an important role in Phanerozoic continental crustal growth. Furthermore, the rates of growth of major continental crustal domains (e.g., west-central United States, Canadian Shield, the Svecokarelian Province of northern Europe, and the Arabian-Nubian Shield) significantly exceed Mesozoic–Cenozoic arc-addition rates (Reymer and Schubert, 1986, 1984), implying island-arc accretion alone is insufficient to account for the amount of crust that was produced in each of these terranes. However, the preservation model totally ignores non-arc contributions that appear to be episodic and also play an important role in juvenile addition to continental crust.

It is widely considered that the rate of continental crustal growth is controlled by mantle convection and global-scale material circulation (Arndt and Davaille, 2013; O'Neill et al., 2007; Condie, 1998; Davies, 1995; Stein and Hofmann, 1994; Reymer and Schubert, 1986). Geological evidence related to supercontinent reconstruction (e.g., Li and Zhong, 2009) shows that both the location and formation of superplumes were dominantly controlled by the first order geometry of global subduction zones. The cyclical nature of supercontinents (Li and Zhong, 2009) and LIP events (Gamal El Dien et al., 2019; Prokoph et al., 2004) thus provide strong supporting evidence for episodic convection of the mantle. Stein and Hofmann (1994) proposed that such a process can account for many hitherto puzzling features of the isotope and trace-element geochemistry of the Earth's mantle and crust. Numerical modeling (Arndt and Davaille, 2013; Steinberger and Torsvik, 2012), geochemical evidence (Wang X-C et al., 2013a, b; Jackson and Carlson, 2011; Jackson et al., 2010), and global supercontinent reconstructions (Li and Zhong, 2009) suggest a cyclical nature for global-scale material circulation between Earth's surface and its interior, leading to pulses of continental crustal growth. Additionally, recent numerical models indicate the rate of continental crustal growth at active continental margins is highly variable and strongly dependent on the evolution of individual arc systems (Zhu et al., 2013) and the type of continental arc (Vogt et al., 2012). Because the evolution of continental arc systems is directly associated with continent configuration and oceanic plate tectonics, continental crustal growth at continental margins is not a uniform process. This is consistent with paleomagnetic data, which present strong evidence for plate- driven episodicity in the Precambrian (Li and Zhong, 2009; Li et al., 2008; O'Neill et al., 2007). The periods of rapid plate motions coincide with the observed peaks in crustal age distribution (O'Neill et al., 2007). O'Neill et al. (2007) and Li and Zhong (2009) considered subduction as the driving force of episodic mantle convection, whereas Arndt and Davaille (2013) proposed that enhanced subduction is a consequence of the arrival in the upper mantle of major mantle plumes. Regardless of the precise interaction of these driving forces, it appears evident that continental crustal growth was most likely an episodic rather than a continuous process in the Precambrian, contrasting with what has been proposed by the preservation model (e.g., Dhuime et al., 2012).

The preservation model also proposes that the age peaks coincide with the assembly of supercontinents (Cawood et al., 2013; Hawkesworth et al., 2010, 2009). However, as shown in Figs. 5a–5c, the zircon age peaks and cyclical nature of the complete zircon dataset are strongly dependent on the nature of their host rocks (mixed in 5a, juvenile in 5b, and reworked in 5c). For example, the age peak at 0.9–0.6 Ga defined by juvenile zircon grains (Fig. 5b) and the age peaks at 0.9–0.7 and 0.3–0.1 Ga defined by the total zircon database (Fig. 5a) broadly correlate with the rifting and breakup of Rodinia (0.85–0.60 Ga) and Pangea (0.32–0.09 Ga), respectively, rather than the assembly stage as proposed by the preservation model. In contrast, the other two age peaks at 2.1–1.7 and 1.3–1.0 Ga defined by the total zircon database (Fig. 5a) and the three major age peaks at 2.1–1.5, 1.3–0.9, 0.7–0.4 Ga (Fig. 5c), defined by zircon grains from reworked crustal materials, broadly correlate with the assembly of Columbia/Nuna (2.1–1.8/1.6 Ga), Rodinia (1.3–0.9 Ga), and Gondwana-Pangea (0.6–0.32 Ga). This is contrary to the main prediction of the preservation model.

The preservation model also proposes that "Globally, subduction is continuous, which in turn suggests that the processes of crust generation should result in a continuum of ages" (Hawkesworth et al., 2009) and "Growth of continental crust appears to have been a continuous process" (Dhuime et al., 2012). However, this is contrary to the evidence presented here that establishes that juvenile zircon age distribution patterns display a cyclical nature (Fig. 5b). The juvenile zircon age peaks at 2.8–2.4, 2.2–1.8, 1.3–1.0, and 0.9–0.6 Ga broadly correlate with rapid crustal growth pulses at 2.9–2.7 Ga (2.7 Ga Superior Province and 2.9–2.8 Ga Sumozero-Kenozero and Kostomuksha orogens), 2.2–2.1 Ga (Birimian Orogen), 1.3–1.0 Ga (Grenville Orogen) and 0.9–0.6 Ga (Pan-African Orogen), which have been attributed to episodes of enhanced exchange between the lower and upper mantle (Stein and Hofmann, 1994). Importantly, this is supported by juvenile zircon age peaks correlating with LIP cycles (Fig. 5d) and OIB ³He/⁴He isotope peaks (mantle depletion events, Fig. 6c).

Furthermore, Bradley (2011) pointed out that the preservation model cannot readily account for the broad similarities between the abundance of passive margins and the age distribution patterns of detrital zircons, or fluctuations in seawater ⁸⁷Sr/⁸⁶Sr.

The weight of evidence thus leads to a rejection of the hypothesis that the U-Pb zircon age peaks correspond to times of enhanced crustal preservation; instead, it indicates episodic convection of the mantle and the cyclical nature of global-scale material circulation, most likely linked to the coupled supercontinent-superplume cycles.

This study has evaluated the role of subducted oceanic crust in the generation of juvenile continental crust and the validity of various isotopic parameters (U-Pb crystallization age, whole-rock Hf model age, and Hf-O isotopes) in constraining continental crustal growth and its reworking. We demonstrate that a mantle- like oxygen isotopic signature is not a diagnostic feature of juvenile continental crust, and that it is inappropriate to use a combination of Hf model age and mantle-like oxygen isotopic signatures of zircon grains to define the growth and evolution of continental crust. A global compilation of whole-rock Hf isotope data demonstrates that juvenile continental crust has a large range in Hf isotope compositions throughout Earth's history. This implies that applying any specific value (e.g., present-day average Hf isotope values of MORB) to evaluate juvenile continental crustal growth in Earth's early history may not be correct. The co-existence of mantle-derived Hf and non-mantle-like oxygen isotopic signatures in single zircon grains provides solid evidence for significant contributions from partially-melted subducted oceanic crust in the generation of juvenile continental crust. Utilizing a global compilation of whole-rock Hf data from mantle-derived rocks, integrated with zircon U-Pb age and O-Hf isotope data from more precisely-defined juvenile crustal igneous rocks, leads us to propose that the U-Pb crystallisation ages of juvenile zircon grains can be used to unravel crustal growth through time and to test its relationship with supercontinent-superplume cycles.

Our new analysis shows that continental crust underwent two stages of rapid growth (3.7–3.5 and 3.0–2.5 Ga), and a long-term, slow but progressive growth stage (post-Archean), punctuated by periods of more rapid growth. What happened prior to this in the Hadean–Eoarchean is difficult to decipher, since it is based on a small zircon database and sparse short-lived ¹⁴⁶Sm-¹⁴²Nd isotope data, hence, we are currently unable to incorporate them in our model. Importantly, the age peaks of juvenile continental crust formation correspond with the timing of global mantle melting (depletion), as recorded by Os and He isotopes, and the intensity of LIP events. Conversely, the zircon U-Pb age peaks for magmas dominantly generated by crustal reworking correlate with the timing of supercontinent assembly. It thus appears that the cyclical nature of Earth's dynamic system involving supercontinent-superplume cycles, provides a first-order control on the episodic evolution of continental crust. Because of the extremely high mantle potential temperature in the Archean to Early Paleoproterozoic, mantle thermal anomalies likely played a much more important role than plate tectonic processes in the growth and evolution of continental crust at that time. This may explain why the early age peaks correlate with both intensive plume activities (komatiite eruption) and assembly of supercontinents, whereas the younger age peaks correlate with the breakup of supercontinents.

Guo-Qiang Tang from the Institute of Geology and Geophysics, Chinese Academy of Sciences is thanked for help with conducting wavelet analysis. This study was supported by the National Key R & D Program of China (No. 2017YFC0601302), the Research Start-up Project for Introduced Talent of Yunnan University (No. 20190043) and the Australian Research Council to Zheng-Xiang Li (Nos. DP0770228, FL150100133). This is a contribution to IGCP 648. The final publication is available at Springer via https://doi.org/10.1007/s12583-020-1077-4.