Abstract

Hafnium (Hf) isotope composition of zircon has been integrated with U-Pb age to form a long-term (>4 b.y.) record of the evolution of the crust. In contrast, trace element compositions of zircon are most commonly utilized in local- or regional-scale petrological studies, and the most noteworthy applications of trace element studies of detrital zircon have been in “fingerprinting” potential source lithologies. The extent to which zircon trace element compositions varied globally over geological time scales (as, for example, zircon U-Pb age abundance, O isotope composition, and Hf isotope composition seem to have varied) has been little explored, and it is a topic that is well suited to the large data sets produced by detrital zircon studies. In this study we present new detrital zircon U-Pb ages and trace element compositions from a continent-scale basin system in Australia (the Centralian Superbasin) that bear directly on the Proterozoic history of Australia and which may be applicable to broader interpretations of plate-tectonic processes in other regions. U-Pb ages of detrital zircon in the Centralian Superbasin are dominated by populations of ca. 1800, 1600, 1200, and 600 Ma, and secular variations of zircon Hf isotope ratios are correlated with some trace element parameters between these major age populations. In particular, elevated εHf(i) (i.e., radiogenic “juvenile” Hf isotope composition) of detrital zircon in the Centralian Superbasin tends to correspond with relatively high values of Yb/U, Ce anomaly, and Lu/Nd (i.e., depletion of light rare earth elements). These correlations seem to be fundamentally governed by three related factors: elemental compatibility in the continental crust versus mantle, the thickness of continental crust, and the contributions of sediment to magmas. Similar trace element versus εHf(i) patterns among a global zircon data set suggest broad applicability. One particularly intriguing aspect of the global zircon data set is a late Neoproterozoic to Cambrian period during which both zircon εHf(i) and Yb/U reached minima, marking an era of anomalous zircon geochemistry that was related to significant contributions from old continental crust.

INTRODUCTION

The long-term record of zircon Hf isotope composition (176Hf/177Hf, commonly expressed as the Hf isotope ratio at the time of zircon crystallization relative to the chondritic uniform reservoir [CHUR], εHf(i)) is periodic over a variety of time scales (Mitchell et al., 2019), but the clearest fluctuations are related to the supercontinent cycle: periods of relatively non-radiogenic (negative) zircon εHf(i) correspond with supercontinent assembly, whereas radiogenic (positive) εHf(i) was most common during supercontinent breakup. A first-order explanation for this observation is as follows. “Crustal reworking” in continental margin subduction zones, closure of ocean basins, and incipient continental collision during supercontinent assembly tend to produce “evolved” (i.e., non-radiogenic) zircons derived from partial melts of crust that have resided in the crust over a long term (e.g., Roberts, 2012; Gardiner et al., 2016). In contrast, rift- and arc-related magmas generated during supercontinent breakup are strongly influenced by new or nearly new mantle additions in the form of basalt, and they tend to produce “juvenile” zircons because of introduction of new mantle-derived material into magmatic systems, potentially including partial melts of crust with a range of crustal residence times (e.g., Smits et al., 2014). In line with these explanations, two of the largest departures of the zircon Hf isotope record from chondritic values (εHf(i) = 0) are strongly radiogenic compositions at ca. 700 Ma (Rodinia breakup) and particularly non-radiogenic compositions at ca. 600 Ma (Gondwana assembly; e.g., Roberts and Spencer, 2014).

Zircon Hf isotope composition is fundamentally governed by differences between the compatibilities of Lu and Hf in the continental crust versus the mantle (Lu/Hf is about two orders of magnitude greater in the continental crust than in mid-ocean ridge basalt [MORB]; Taylor and McLennan, 1985; Sun and McDonough, 1989; Vervoort, 2014). Because some zircon trace elements also have wide disparities in crust-mantle compatibility, their ratios too may correlate with εHf(i), and they may be potential sources of information about source magma composition that complement zircon Hf isotope ratios. Moreover, they may serve to make depleted mantle (DM) model ages more accurate if they help confine the choices of Lu/Hf selected for the “pre-zircon” history of those elements in the first stage of their ingrowth in the crust. For example, the concentration ratio Lu/Hf would seem to be a logical candidate for a zircon trace element ratio that has secular covariance with zircon εHf(i), although the extent of Hf substitution for Zr in zircon is strongly influenced by fractional crystallization of zircon itself in true granites (Claiborne et al., 2006; Linnen and Keppler, 2002; Grimes et al., 2007), whereas fractionation of these elements in more mafic rocks is relatively unexplored. Among other elements that are commonly measured in zircon, Yb (a heavy and abundant rare earth element [REE]) and U (both of which substitute for Zr in the zircon structure, albeit by different mechanisms; e.g., Grimes et al., 2007) also have particularly large compatibility differences in the continental crust versus mantle, as demonstrated by the fact that Yb/U is ∼25× greater in normal MORB than in bulk continental crust (e.g., Taylor and McLennan, 1985; Hofmann, 1988; Sun and McDonough, 1989; Grimes et al., 2007). Zircon-melt distribution coefficients for U and Yb (Bea et al., 1994) differ by a factor of only ∼1.3, so the bulk rock (i.e., assumed magma) range of Yb/U seems to be closely reflected in the trace element composition of zircon, as evidenced by substantial Yb/U differences between zircons from continental arcs versus those from MORB (Grimes et al., 2007, 2015). Zircon Yb/U is also affected by fractional crystallization (Yb/U decreases as fractional crystallization progresses), although fractional crystallization models suggest that zircon Yb/U is strongly related to the initial Yb/U of the melt (Reimink et al., 2019). In addition, metamorphic zircons commonly have low Yb/U due to competition for heavy REEs (HREEs) with garnet (e.g., Hoskin and Ireland, 2000; Rubatto, 2002; Hoskin and Schaltegger, 2003; Rubatto et al., 2011; Kylander-Clark et al., 2013). The very low Yb/U of some metamorphic zircon (e.g., Kylander-Clark et al., 2013; Takamura et al., 2020) is essentially the result of double distillation: first during formation of continental crust, and secondly during zircon crystallization in the continental crust in the presence of garnet.

For zircon trace element parameters other than Yb/U that may not necessarily be related to crust-mantle compatibility, one avenue of evaluating the underlying causes of their potential secular variability is to simply compare their long-term records with the zircon Hf isotope record. The empirical behaviors of these trace element parameters during the fundamental tectonic processes listed above (arc magmatism, rifting, etc.) could be inferred from a first-order interpretation of zircon Hf isotope variability. To explore these possibilities, we measured U-Pb detrital zircon ages and trace element concentrations from sandstone samples in the Georgina, Amadeus, Warburton, and Officer Basins in central Australia (Fig. 1). These new data form time-integrated zircon trace element records that we compare with the correlative detrital zircon Hf isotope record of central Australia. We divide the following discussion into two parts that reflect two types of commonly encountered zircon geochemistry studies: a regional-scale detrital zircon study (central Australian Neoproterozoic–Paleozoic basins, in this case) and an evaluation of compiled data from zircons that span Earth history.

INTERPRETED TECTONIC HISTORY OF CENTRAL AUSTRALIA

Much of the Australian continent is covered by a series of Neoproterozoic to Paleozoic basins that are collectively referred to as the Centralian Superbasin (Walter et al., 1995). The various components of the Centralian Superbasin are separated by basement-cored uplifts (Fig. 1) in a present-day geographic configuration similar to that of Laramide (Late Cretaceous–early Cenozoic) basins and basement-cored uplifts of the western North American Cordillera. The most areally significant constituent basins of the Centralian Superbasin are the Amadeus Basin in central Australia, the Officer Basin to the southwest, and the Georgina Basin to the northeast (Fig. 1). The Centralian Superbasin comprises a series of latest Mesoproterozoic or early Neoproterozoic to approximately Devonian depocenters, and numerous studies have established broad stratigraphic correlation among and between the individual basins (e.g., Preiss et al., 1978; Walter et al., 1995; Grey and Calver, 2007; Verdel and Campbell, 2017; Normington and Donnellan, 2020). The two major basement-cored uplifts that separate these basins in central Australia are the Musgrave province and Arunta region (Fig. 1), and two of the most notable periods of uplift of these inliers occurred during the Petermann (Ediacaran–Cambrian) and Alice Springs (mid-Paleozoic) orogenies (e.g., Haines et al., 2001; Raimondo et al., 2010; Glorie et al., 2017; Wex et al., 2017). Much, but not all, of the sedimentation in the Centralian Superbasin can be attributed to foreland basin deposition during these two orogenic periods.

The detrital zircon geochronology record of the Centralian Superbasin is one of the more extensively examined in the world. Roughly 7000 detrital zircon U-Pb analyses (essentially all of which were originally collected to constrain maximum depositional ages of various stratigraphic units) have been published from Neoproterozoic to Cambrian strata of the Amadeus, Georgina, and Officer Basins, in addition to detrital zircon data sets from broadly correlative strata of the Adelaide Rift Complex in South Australia (Fig. 1; e.g., Ireland et al., 1998; Rose et al., 2013; Keeman et al., 2020; Lloyd et al., 2020). Some of the previous Amadeus Basin detrital zircon studies (Hollis et al., 2013; Haines et al., 2016; Normington et al., 2019) included Hf isotope compositions, which, as discussed in more detail below, compose a reasonably robust zircon Hf isotope record consisting of >1000 analyses that span >2 b.y. of geologic time and >25 εHf(i) units.

Arunta Region

The detrital zircon chronostratigraphy of the Centralian Superbasin informs the histories of neighboring Archean to Mesoproterozoic basement inliers, principally the Arunta region and Musgrave province (e.g., Collins and Shaw, 1995; Wade et al., 2008; Hollis et al., 2013; Kirkland et al., 2013; Haines et al., 2016). Basement rocks of the Arunta region form the southernmost part of the North Australian craton (Fig. 1), and they are conventionally divided into two provinces: the Aileron province in the north, and a relatively narrow sliver along the southern margin called the Warumpi province (Fig. 1; e.g., Ahmad and Scrimgeour, 2013; Scrimgeour, 2013a, 2013b). A key (although not entirely unequivocal) distinction between these two provinces is age: the Aileron province is generally >1700 Ma, and the Warumpi province is mainly <1700 Ma. Within the Aileron province, one of the most widespread episodes of deformation and magmatism occurred during a late Paleoproterozoic (ca. 1790–1770 Ma) interval that has been variably referred to as the early Strangways orogeny (Collins and Shaw, 1995) and the Yambah event (Fig. 2; Claoué-Long and Hoatson, 2005). The Aileron and Warumpi provinces are separated by the Central Australian suture (Fig. 1). The major period of tectonism along this structure was the ca. 1640 Ma Liebig orogeny (Scrimgeour et al., 2005), which has been interpreted as the accretion of the Aileron and Warumpi provinces (Fig. 2; e.g., Scrimgeour et al., 2005; Betts et al., 2008; Ahmad and Scrimgeour, 2013; Hollis et al., 2013). Both provinces were subsequently affected by a ca. 1580 Ma period of metamorphism (the Chewings orogeny of Teyssier et al., 1988, and Collins and Shaw, 1995) that produced widely documented metamorphic zircon overgrowths (Vry et al., 1996; Williams et al., 1996; Rubatto et al., 2001; Rubatto, 2002; Scrimgeour, 2013a).

In addition to age, a notable difference between the Aileron and Warumpi provinces is zircon Hf isotope composition. Zircons from the older Aileron province have a wide range of εHf(i) from approximately −10 to >+10 (e.g., Hollis et al., 2013; Beyer et al., 2015, 2018; Whelan and Armstrong, 2015), and, in contrast, εHf(i) of zircons from crystalline rocks of the Warumpi province fall into a partially overlapping but different and slightly more restricted range of ∼0 to +15 (Hollis et al., 2013). There is thus an overall transition across the Liebig orogeny toward more radiogenic (i.e., more mantle-like) zircon Hf isotope composition within the Arunta region, and perhaps in other basement provinces of Australia as well. The tectonic model of Hollis et al. (2013) offers a potential explanation for the observed Hf isotope trend. In this model, the substrate of the Warumpi province was rifted from the Aileron province at ca. 1700 Ma, opening an ocean basin that separated the autochthonous Aileron province from the allochthonous continental ribbon of the Warumpi province (Fig. 2). This ocean basin subsequently subducted southward (present coordinates) beneath the Warumpi province until eventual collision of the two provinces during the Liebig orogeny (Fig. 2). According to this model, <1700 Ma magmatism of the Warumpi province characterized by radiogenic zircon is interpreted as arc magmatism.

Musgrave Province

The Musgrave province along the Northern Territory–South Australia border is roughly 200 km south of the Warumpi province and is separated from it by the Amadeus Basin (Fig. 1). It lies along the boundary of the North Australian and South Australian cratons (Fig. 1; e.g., Wade et al., 2008; Smithies et al., 2011). Like the Arunta region, the Musgrave province has a long and complex geological history, which, in the case of the Musgrave province, stretches back to perhaps ca. 1900 Ma (e.g., Kirkland et al., 2013). Much of the Musgrave province consists of calc-alkaline intrusive rocks that have been attributed to a protracted period of subduction from ca. 1600 to ca. 1300 Ma that was terminated during a ca. 1300 Ma phase of contraction (Fig. 2; e.g., Edgoose et al., 2004; Wade et al., 2006, 2008; Smithies et al., 2011; Kirkland et al., 2013; Smits et al., 2014; Yang et al., 2019). Zircons in the Musgrave province are generally younger than those in the Arunta region (e.g., Kirkland et al., 2013), and a particularly noteworthy aspect of the Musgrave province is extensive ca. 1200 Ma granitic magmatism of the Pitjantjatjara Supersuite (Fig. 2; e.g., Smithies et al., 2011), which was coeval with a period of metamorphism and deformation in central Australia known as the Musgrave orogeny (e.g., Wade et al., 2008). The Musgrave province ca. 1200 Ma granitoids have relatively non-radiogenic Nd isotope compositions and have been ascribed to crustal anataxis (Wade et al., 2008).

Subsequent bimodal magmatism occurred in the Musgrave province from ca. 1100 to 1000 Ma (the Giles event) and was followed by the emplacement of mafic dikes at ca. 1000 and 800 Ma (Fig. 2; e.g., Zhao et al., 1994; Wade et al., 2008; Evins et al., 2010; Donnellan et al., 2019). The Musgrave province is commonly considered (in terms of age) the major “Grenvillian” belt of Australia (e.g., Aitken and Betts, 2008; Walsh et al., 2015), although it is doubtful that it was a straightforward along-strike Rodinia continuation of the actual Grenville province of North America, which, instead, likely passes through Antarctica (e.g., Storey et al., 1994; Goodge et al., 2008). Regardless, the Musgrave province is commonly cited as the source for the vast amount of ca. 1200 Ma “Grenvillian-age” detrital zircons across Australia (e.g., Camacho et al., 2002; Haines et al., 2016; Purdy et al., 2016; Armandola et al., 2018). Importantly, ca. 1200 Ma magmatism in the Musgrave province, which is characterized by zircon with a particularly wide range of εHf(i), may have been related to a period of rifting, as was subsequent ca. 1100–1000 Ma bimodal magmatism and the ca. 800 Ma emplacement of mafic dikes (Fig. 2; Zhao et al., 1994; Walter et al., 1995; Close et al., 2003; Wade et al., 2008; Evins et al., 2010; Smithies et al., 2011; Kirkland et al., 2013). Extensional deformation during this time frame has also been described from the Warumpi province (Morrissey et al., 2011). This interpretation suggests an overall rifting framework from ca. 1200 to 800 Ma for central Australia (Fig. 2), the final part of which would be what is commonly referred to as Rodinia breakup.

Like in the Arunta region, the zircon Hf isotope record of the Musgrave province is informative for deciphering its tectonomagmatic history (Kirkland et al., 2013; Smits et al., 2014). Zircons as old as ca. 1600 Ma are common in felsic gneiss of the Musgrave province (e.g., Edgoose et al., 2004), and these zircons are strongly radiogenic, even more so than roughly time-equivalent zircon of the Warumpi province. Based on this difference, Kirkland et al. (2013) inferred that the Warumpi and Musgrave provinces are distinct from each other, i.e., fragments of the Warumpi province do not form the basement of the Musgrave province. At ca. 1600 Ma, the Musgrave province was likely a magmatic arc (e.g., Wade et al., 2006), and subduction continued beneath it until ca. 1300–1200 Ma (Fig. 2). In other words, ca. 1700–1600 Ma magmatism of the Warumpi and Musgrave provinces may have been related to two separate arc systems (perhaps of opposite polarity) that are now positioned on the northern and southern margins of the Amadeus Basin (Fig. 2). In both provinces, ca. 1600 Ma magmatism is marked by radiogenic zircon εHf(i) (Hollis et al., 2013; Kirkland et al., 2013).

DETRITAL ZIRCON Hf ISOTOPE RECORD OF THE AMADEUS BASIN AND NEIGHBORING AREAS

In addition to the major components of the Centralian Superbasin, there are a number of smaller constituent basins (Fig. 1), including the Ngalia Basin (which is situated within the Aileron province itself), the Murraba Basin of Western Australia, and the Irindina province of central Australia, a component of the Arunta region that seems to include metamorphosed equivalents of Amadeus Basin strata (Scrimgeour, 2013c). The Warburton Basin of central and southern Australia (Fig. 1; e.g., Edgoose and Munson, 2013) is entirely within the subsurface but seemingly forms a geographic and stratigraphic link between the Centralian Superbasin and age-equivalent metasediment of eastern Australia (Purdy et al., 2016). In the discussion of compiled zircon Hf isotope results that follows, we focus on the detrital zircon Hf isotope record of the Amadeus Basin (Hollis et al., 2013; Haines et al., 2016; Normington et al., 2019), but we also include a relatively limited detrital zircon Hf isotope data set from the Murraba Basin (Hollis et al., 2013) as well as detrital zircon results from a sample of Ordovician sandstone in the Officer Basin (Haines et al., 2013; Morón et al., 2019). Additionally, some of the new detrital zircon data presented below are from the western part of the Warburton Basin.

The compiled central Australian detrital zircon Hf isotope record outlined above comprises >1400 individual measurements and spans from ca. 3400 to 500 Ma (Fig. 3). The majority of zircon εHf(i) in this data set fall between approximately −15 and +10, and outliers range from −39 to +19. Zircons with U-Pb ages between ca. 3400 and 1900 Ma are relatively uncommon; most of these have εHf(i) between about −5 and +5, although some with U-Pb ages between ca. 2500 and 2400 Ma suggest a departure to strongly non-radiogenic (crustal) εHf(i) of −20. A significant cluster of detrital zircon U-Pb ages occurs between ca. 1900 and 1500 Ma, and εHf(i) of these zircons seem to diverge at ca. 1640 Ma. Detrital zircons between ca. 1900 and 1640 Ma are grouped around ca. 1770 Ma, and they have εHf(i) between about −9 and +9. Between ca. 1640 and 1500 Ma, zircon εHf(i) is more radiogenic that the immediately preceding period, with most values falling between ∼0 and +10. The ca. 1640 Ma divergence thus corresponds with a diminution of non-radiogenic detrital zircon in central Australia. The detrital zircon record is again sparse between ca. 1500 and 1200 Ma, although most zircons in this range have εHf(i) >0. Another major cluster of zircon U-Pb ages occurs between ca. 1200 and 1000 Ma, most of which have εHf(i) between −13 and +10. This late Mesoproterozoic group therefore has a wide range of Hf isotope composition like the preceding ca. 1900–1500 Ma cluster, although the younger group seems to trend toward less radiogenic values. Finally, zircons from ca. 700–500 Ma are characterized by a particularly large range of εHf(i), including values as low as about −30 (Fig. 3).

The central Australian detrital zircon Hf isotope record has been interpreted by previous researchers as follows. Late Paleoproterozoic (ca. 1800 Ma) zircons, which are likely derived primarily from the Aileron province, have an εHf(i) range of ∼20 units (i.e., a “vertical array;” Fig. 3), and thus seem to reflect a combination of melting of Archean crust and a Paleoproterozoic input of juvenile magmatism (basalts; Smits et al., 2014; Beyer et al., 2018). The cluster of detrital zircon U-Pb ages ca. 1770 Ma is coeval with significant granitic and mafic magmatism in the southern part of the Aileron province (i.e., the Yambah event) that seems to have formed in an arc or back-arc environment (Fig. 2; Collins and Shaw, 1995; Claoué-Long and Hoatson, 2005; Maidment et al., 2005; Weisheit et al., 2019). In the broadest terms, these late Paleoproterozoic zircons coincide with the assembly of both the North Australian craton and the Nuna supercontinent (e.g., Betts et al., 2015; Kirscher et al., 2021). The ca. 1640 Ma shift to radiogenic Hf isotope composition seems to coincide with the Liebig orogeny (Fig. 2; Scrimgeour et al., 2005; Smits et al., 2014). The ca. 1200–1000 Ma group of detrital zircons, which (1) is derived primarily from the Musgrave province (Fig. 1), (2) likely formed during a period of crustal extension (Fig. 2), and (3) has a particularly strong Hf isotope vertical array spanning >20 epsilon units (Fig. 3), has been attributed to the combined melting of Archean and Proterozoic crust (Smits et al., 2014). Finally, ca. 700–500 Ma zircons are prevalent in early Paleozoic strata of Australia, but they seem to be largely exotic (e.g., Ireland et al., 1998; Cawood and Nemchin, 2000; Squire et al., 2006; Maidment et al., 2007; Haines et al., 2013; Morón et al., 2019). The extraordinary geochemical compositions of the 700–500 Ma zircons are discussed later in this paper.

With the exception of the youngest group of ca. 700–500 Ma zircons, the central Australian detrital zircon Hf isotope record shares key similarities with combined records from the Aileron, Warumpi, and Musgrave provinces that are reviewed above (e.g., Hollis et al., 2013; Kirkland et al., 2013; Haines et al., 2016; Beyer et al., 2018). The vast majority of zircon U-Pb ages in these records are between ca. 2000 and 1000 Ma, and in this interval εHf(i) values define an “inverted U-shaped pattern” (Smits et al., 2014) that consists of vertical arrays at ca. 1800 and 1200 Ma, separated by an interlude of predominantly radiogenic Hf isotope composition. The global zircon Hf isotope record (Fig. 3; e.g., Belousova et al., 2010; Roberts and Spencer, 2014) follows the same pattern during this time period, which spans much of Earth’s “middle age” (e.g., Cawood and Hawkesworth, 2014).

SAMPLE LOCATIONS, METHODOLOGY, AND RATIONALE

Sample Locations

Detrital zircon data were collected from 13 samples from the Georgina, Amadeus, Officer, and Warburton Basins (Fig. 1; Table S11). Two samples are from Neoproterozoic glacigenic deposits of the southern Georgina Basin. One of these samples is from the early Cryogenian (i.e., Sturtian) Yardida Tillite of the Georgina Basin in drillhole MW005A (Davies et al., 2015), and the other is from an outcrop of the late Cryogenian (i.e., Marinoan) Black Stump Arkose. Six samples are from the Amadeus Basin. The oldest of these is from the Sturtian Areyonga Formation in the far eastern Amadeus Basin. Three samples are from Marinoan glacially influenced strata of the Olympic Formation (e.g., Preiss et al., 1978), including a sample from its type locality (Wells et al., 1967; Kennedy, 1996). One of the Amadeus Basin samples is from a Cryogenian glacial deposit in the northeastern Amadeus Basin that could be either the older Areyonga Formation or the younger Olympic Formation, although the discovery of a clast of what appears to be the mid-Cryogenian Aralka Formation in these sediments suggests that they are part of the late Cryogenian Olympic Formation (Allen et al., 2018). The youngest sample is sand from a modern creek bed at Trephina Gorge in the northeastern Amadeus Basin. Three samples are from the Officer Basin, all of which are from the Marinoan Wahlgu Formation. Finally, two samples are from early Paleozoic sandstone (presumed Hugh River Shale) intersected in a drillhole (McDills 1, drilled by Amerada Petroleum Corporation of Australia in 1965) in the northwestern part of the Warburton Basin (e.g., Edgoose and Munson, 2013).

U-Pb Geochronology and Geochemistry Methodology

Zircons were separated from the samples using standard magnetic and gravimetric techniques. Zircons were analyzed at the Queensland University of Technology (Brisbane, Queensland, Australia) using the methodology outlined in Campbell et al. (2020a). Discordance was calculated by comparing 238U/206Pb and 207Pb/206Pb ages for zircons >1000 Ma, and 238U/206Pb versus 235U/207Pb ages for zircons ≤1000 Ma. Zircons with ages >10% discordant were not considered for further analysis. Yb concentration of zircons was estimated from the concentration of Lu to facilitate comparison of our new data with those of Grimes et al. (2007) and other compilations. To determine an average Yb/Lu relationship, we compiled >9000 zircon analyses from the literature that include both Lu and Yb concentrations. Nearly every zircon in this data set has Yb/Lu between 4 and 10. Linear regression of Yb and Lu suggests an average relationship of [Yb] ≈ 8 × [Lu], which was used to estimate Yb concentration based on Lu concentration. Yb concentrations were estimated for some compiled zircon data sets in the same way. The uncertainty in Yb concentration estimate is mitigated considerably by evaluation of these data on log scales, as described below.

Cerium (Ce) and Eu anomalies were calculated using a version of the lattice strain method of Onuma et al. (1968) and Ballard et al. (2002). Use of this method is driven largely by the fact that [La] is at or near detection limits in many zircons, precluding accurate determination of Ce anomaly by interpolation between La and Pr. For each zircon, concentrations of Ce, Nd, Eu, Dy, and Lu were normalized by their concentrations in average crust (Taylor and McLennan, 1985). A straight line was fit to three pairs of numbers that compose a “trivalent array”: the natural logarithms of the normalized concentrations of Nd, Dy, and Lu (i.e., the natural logarithms of the estimated zircon distribution coefficients for these elements normalized by those of average crust), and the natural logarithms of their corresponding lattice strain parameter (Ballard et al., 2002), thereby establishing a base slope. “Hypothetical” normalized concentrations of Ce and Eu were determined by projection along this slope to the lattice strain parameters of Ce and Eu. These hypothetical trivalent array concentrations are estimates of Ce and Eu concentrations if both were entirely trivalent. Ce and Eu anomalies, as quantified in this paper, refer to the ratios of measured Ce and Eu concentrations to their estimated trivalent array concentrations. Although their magnitudes depart from those of some other methods, Ce and Eu anomalies determined using the lattice strain method have the same sense as other methods: Ce peak becomes taller as Ce anomaly increases (implying increasing Ce4+/Ce3+), and Eu trough deepens as Eu anomaly becomes smaller (implying increasing Eu2+/Eu3+; Fig. S1 [footnote 1]). If zircon Ce and Eu anomalies are controlled entirely by oxidation state, their variations should thus be in phase.

Zircon Trace Element Parameters and Rationale

We focus on four key zircon trace element parameters. The first is Yb/U, a ratio that has been previously used to distinguish detrital zircon derived from oceanic versus continental crust (Grimes et al., 2007, 2015). As described above, the disparity in continental versus oceanic zircon Yb/U seems largely attributable to compatibility differences between Yb and U (Table 1).

The second zircon trace element parameter is light REE (LREE) depletion (Table 1). One of the primary features of zircon REE composition is a pronounced depletion in LREEs relative to HREEs (e.g., Hoskin and Ireland, 2000), a pattern that arises because the ionic radii of HREEs are closer to that of Zr than are those of LREEs (e.g., Hanchar and van Westerenen, 2007). The degree of LREE depletion in zircon can be related to LREE depletion of the melt from which it crystallized (e.g., Sawaki et al., 2017; McKenzie et al., 2018). Numerous factors affect REE compositions of melts, but two observations are particularly relevant to zircon REE geochemistry. The first is that MORBs are strongly depleted in LREEs relative to intermediate and felsic rocks of the continental crust (Taylor and McLennan, 1985; Sun and McDonough, 1989). The second is that in continental arcs, the ratio LREEs/HREEs corresponds with crustal thickness, a correlation that has been ascribed to the influence of garnet and amphibole, both of which are HREE enriched (e.g., Gromet and Silver, 1987; Rubatto, 2002; Whitehouse, 2003; Farner and Lee, 2017; Kotková et al., 2016; McKenzie et al., 2018). A first-order control on zircon REE slope from continental arcs thus seems to be melting and/or fractionation in the presence of garnet- and/or amphibole-rich residues, which require pressures >1 GPa. In other words, zircons from deep crustal melts, where a thickened crust is a necessary condition, are thought to have relatively high LREEs/HREEs due to HREE sequestration by garnet and amphibole. LREEs/HREEs of both rocks and zircons have therefore been used in studies examining secular trends in the thickness of continental crust (Keller and Schoene, 2012; Barth et al., 2013; McKenzie et al., 2018). We use Lu/Nd as a simple measure of the extent of LREE depletion. Based on the factors summarized above, relatively low Lu/Nd (i.e., relatively high LREEs/HREEs) of zircon might be expected to reflect extensive contributions from the continental crust and/or deep melting within thick continental arcs (Table 1). In either case, there is reason to think of zircon Lu/Nd as inversely proportional to continental crust involvement.

The third zircon trace element parameter we evaluate is Ce anomaly (Ce/Ce*; Table 1). “Positive” Ce anomaly (i.e., a peak in Ce on CHUR-normalized REE diagrams; Ce/Ce*>1) is a ubiquitous characteristic of zircon that is attributable to the similarities in charge and ionic radii of Ce4+ and Zr (Ballard et al., 2002). Zircon effectively incorporates all available Ce4+. Zircon Ce anomaly is therefore related to the oxidation state of magmas (larger Ce anomalies correspond with more oxidized magmas; e.g., Ballard et al., 2002; Trail et al., 2012; Smythe and Brenan, 2016), which, in turn, may be influenced by several processes. For example, subduction zones are widely believed to be regions of mantle oxidation (e.g., Ballhaus et al., 1990; Brandon and Draper, 1996; Rowe et al., 2009), suggesting that zircons from magmatic arcs might have greater Ce anomaly than zircons from some other tectonic environments. Another factor that could be an important influence on the Ce anomaly of zircon is the incorporation of sediment in parental melts. The underlying cause of the magnetite- versus ilmenite-bearing granite distinction (which is closely related to the I- and S-type classifications; e.g., Chappell and White, 1974, 2001; Ishihara, 1977, Takagi and Tsukimura, 1997) is thought to be the inclusion of graphite from sedimentary rocks during the formation of “reduced” granites (Flood and Shaw, 1975). Indeed, in the Lachlan orogen of eastern Australia, Ce anomaly of zircon is greater for I-type granites than for S-type (Burnham and Berry, 2017). A third and perhaps related factor that could be important in influencing the Ce anomaly of zircon is crustal thickness. In magmatic arcs, oxygen fugacity of melts decreases as arc thickness increases (Farner and Lee, 2017). This relationship may be fundamentally a reflection of a difference between relatively thin island arcs (relatively oxidized) versus thick continental arcs (relatively reduced). Ce anomaly of zircon also increases as crystallization temperature decreases (Trail et al., 2012). Finally, Ce anomaly of zircon may be affected by monazite Ce (Ce3+) sequestration. For example, deep crustal granites containing both monazite and zircon may produce zircon with relatively low (i.e., relatively flat) Ce anomaly (Gaschnig, 2019).

The final key trace element parameter we evaluate is Eu anomaly (Eu/Eu*; Table 1). Like Ce, Eu occurs in two oxidation states (2+ and 3+ in the case of Eu). Because Eu3+ preferentially substitutes for Zr over Eu2+ in a manner similar to the preferential incorporation of the oxidized form of Ce, Eu anomaly may be, in part, a redox indicator (Ballard et al., 2002; Trail et al., 2012). In this case, strongly positive Ce anomaly (a nominal indicator of prevalent Ce4+ from an oxidized source) should correspond with incorporation of Eu3+ and a minimal “negative” Eu anomaly (i.e., Eu/Eu* <1). The “positive” and “negative” terminology surrounding Ce and Eu anomalies is confusing because the numerical values of both are always >0, so, in accordance with their use as ratios of non-trivalent to trivalent ions, we refer to relatively flat (what might otherwise be called “weakly negative”) Eu anomalies as “high,” and deep (“strongly negative”) Eu anomalies as “low.” Another important consideration in interpreting zircon Eu anomaly is Eu sequestration by plagioclase (e.g., Hoskin and Schaltegger, 2003; Rubatto et al., 2011). Zircons that begin to form after the onset of plagioclase crystallization would experience a plagioclase-fractionated melt composition and thus a lower (i.e., deeper) Eu anomaly. Suppression of plagioclase crystallization due to melting of the deep crust at pressures greater than that of the plagioclase stability field therefore has the opposite effect of producing high (i.e., flat) zircon Eu anomalies (Keller and Schoene, 2012). Large compilations of REE data from both felsic rocks and zircon (Keller and Schoene, 2012; McKenzie et al., 2018), illustrate that, like LREEs/HREEs, Eu anomaly decreased (i.e., Eu troughs became deeper) over much of Earth history, reaching a nadir during the early part of the Neoproterozoic Era, before increasing during the Phanerozoic. Smaller-magnitude variation of Eu anomaly may be related to fractional crystallization (Campbell et al., 2020a), and titanite crystallization may have the effect of increasing Eu anomaly of zircon (Loader et al., 2017).

The four key zircon trace element parameters described above (Yb/U, Lu/Nd, Ce anomaly, and Eu anomaly) form a framework for interpreting zircon REE compositions in terms of tectonic and magmatic processes (Table 1). Based on empirical differences between MORB and continental arc compositions, variation in zircon Yb/U seems to be closely related to inputs from the DM versus the continental crust (basalts versus granites, respectively), as are variations in zircon LREE depletion. For zircons from magmatic arcs, arc thickness (which is essentially synonymous with extent of continental crust involvement) may be an important influence on both the degree of LREE depletion and Ce anomaly. Eu anomaly of zircon is sensitive to redox conditions (like Ce anomaly) and pressure (like LREE depletion). Given that the underlying causes of variation in these parameters seem to be interrelated (Table 1), it might be expected that secular trends of the parameters would be correlated (e.g., Balica et al., 2020). Moreover, given that variation of zircon Hf isotope composition is fundamentally a measure of contributions of contemporaneous mantle melts versus partial melts of continental crust extracted from the mantle in the significant past (e.g., Vervoort, 2014), zircon Hf isotope composition might also be expected to correlate with these trace element parameters. Much of the remainder of this paper is focused on empirical evaluation of the multivariate zircon geochemistry record.

RESULTS

The 13 samples described above produced 1055 concordant zircon U-Pb ages (Fig. 4; Table S2 [footnote 1]). The major peaks in the U-Pb age spectra of individual samples occur at ca. 1800, 1600, 1200, and 600 Ma (Fig. 4). The youngest zircons from Neoproterozoic samples do not constrain depositional age beyond that of existing global Neoproterozoic age models (Rooney et al., 2015) in conjunction with depositional age estimates for central Australia based on carbon isotope chemostratigraphy (Verdel and Campbell, 2017). The maximum depositional age for Cambrian–Ordovician sandstone in the Warburton Basin (the two samples from the McDills 1 drillhole) provides a meaningful constraint, however. The youngest single concordant zircon from these two samples is 496.2 ± 7.6 Ma (2σ). The five youngest concordant zircons in the two samples have U-Pb ages that overlap within 2σ, and their weighted mean age is 501.8 ± 3.9 Ma (2σ). The maximum depositional age of these Warburton Basin sandstones is therefore late Cambrian, a result that is consistent with an inferred depositional age of ca. 503 Ma for the Hugh River Shale (Schmid, 2017).

DISCUSSION

Detrital Zircon Stratigraphy

The overall stratigraphic pattern of detrital zircon U-Pb results (Fig. 4) is similar to that previously described from the Amadeus Basin (Maidment et al., 2007). Samples of Marinoan glacial deposits (Olympic Formation) from the Amadeus Basin have strikingly different detrital zircon age spectra than correlative strata (Wahlgu Formation) in the Officer Basin, reflecting differences in provenance. The combination of new and previous results from the Georgina Basin is intriguing because the emergence of ca. 1000 Ma zircons in the Neoproterozoic–Cambrian detrital zircon stratigraphy of the southern Georgina Basin seems to reflect uplift of the Musgrave province during the Ediacaran to Cambrian Petermann orogeny (Fig. 5). The same pattern does not hold for other parts of the Centralian Superbasin, however, and, in fact, strata in the Amadeus Basin as old as Tonian contain zircons that seem to have been derived from the Musgrave province (e.g., Kositcin et al., 2015). The Georgina Basin detrital zircon stratigraphy may be particularly sensitive to the effects of the Petermann orogeny uplift because it is geographically isolated from the Musgrave province in comparison with the Amadeus and Officer Basins (Fig. 1).

Central Australia Detrital Zircon Trace Element Records, and Comparison with the Central Australia Detrital Zircon Hf Isotope Record

Most of the new detrital zircon ages are ≤1850 Ma (Fig. 4), and we divide these into seven age populations that correspond with the seven age intervals of Figure 2: 1850–1640 Ma (late Paleoproterozoic), 1640–1500 Ma (latest Paleoproterozoic to early Mesoproterozoic), 1500–1250 Ma (mid-Mesoproterozoic), 1250–1050 Ma (late Mesoproterozoic), 1050–850 Ma (latest Mesoproterozoic to Tonian), 850–700 Ma (Tonian to Cryogenian), and 700–500 Ma (Cryogenian to Cambrian). The oldest group of 1850–1640 Ma zircons composes one of the dominant populations in central Australia. Relative to the subsequent age population of ca. 1640–1500 Ma, the older zircons have low Yb/U, low Lu/Nd, and smaller Ce anomaly (Fig. 6). As described above, the ca. 1640 Ma Liebig orogeny of the Arunta region corresponds with a sharp rise of ∼15 epsilon Hf units in zircon from ca. 1700 to 1600 Ma. This transition is mimicked by changes in zircon trace element parameters such as Yb/U, Lu/Nd, and Ce anomaly (Fig. 6). In this interval, younger (1640–1500 Ma) zircons have more radiogenic εHf(i) (suggesting more mantle-like magmatism) and greater Yb/U (more MORB-like). The younger zircons also have slightly greater Lu/Nd (i.e., greater LREE depletion), which is consistent with the patterns of εHf(i) and Yb/U because (1) MORB has greater Lu/Nd than continental crust, and (2) deep crustal melting in thick continental arcs sequester HREEs such as Lu from zircon. Finally, Ce anomaly is also greater in the younger zircon group (Fig. 6), a correspondence that suggests that zircon Ce anomaly can be influenced by contributions of continental crust versus mantle. The apparent inverse relationship between Ce anomaly and continental crust involvement may be attributable to the same processes responsible for ilmenite- versus magnetite-bearing granites, namely that incorporation of sediments results in relatively reduced melts and zircon with relatively low Ce anomaly.

The ca. 1640–1500 Ma zircons as a group have ages that are centered at ca. 1590 Ma. They could therefore potentially be derived from either early Mesoproterozoic metamorphic zircons associated with the Chewings orogeny and/or the Musgrave province, although ca. 1640–1500 Ma zircons seem to be much more prevalent in the Musgrave province than in the Arunta region (e.g., Kirkland et al., 2013; Haines et al., 2016). In the rock record, ca. 1600–1540 Ma calc-alkaline orthogneiss with radiogenic Nd isotope composition is well documented in the Musgrave province (Camacho and Fanning, 1995; Edgoose et al., 2004; Wade et al., 2006). Regardless of the provenance of the ca. 1640–1500 Ma group of detrital zircons, an important observation is that by a variety of geochemical metrics (radiogenic εHf(i), high Yb/U, “oxidized” Ce anomaly, and significantly depleted LREEs), these zircons were less influenced by continental crust than ca. 1850–1640 Ma predecessors presumed to be derived from the Aileron province. Zircon trace element compositions thus support prior inference of a juvenile contribution to ca. 1650–1550 magmatism in the Musgrave province (Kirkland et al., 2013). This latest Paleoproterozoic–earliest Mesoproterozoic period of juvenile magmatism in central Australia could be interpreted as the signal of a new continental arc at the boundary of the North Australian and South Australian cratons (Fig. 2; Wade et al., 2006).

The third age interval (1500–1250 Ma) is not well represented in the central Australia detrital zircon data set (Fig. 6) because regional magmatism was relatively limited (although not absent; Fig. 2) during this period. Zircon parameters during this age interval seem to be broadly similar to those of the preceding 1640–1500 Ma population. The fourth age group (1250–1050 Ma) comprises the “Grenvillian” population of central Australia that largely coincides with the Musgrave orogeny. In terms of both trace element and Hf isotope composition, these late Mesoproterozoic zircons are similar to the 1640–1500 Ma group (Fig. 6). The geochemical compositions of both the 1250–1050 Ma and 1640–1500 Ma groups suggest a juvenile component, although this component could be rift related in one case and nascent arc related in the other (Fig. 2). A potentially important attribute of the 1250–1050 Ma group is that it includes a tail of ca. 1070 Ma zircons with deep (i.e., low) Eu anomalies (Fig. 6), indicative of significant involvement of plagioclase. One interpretation of the deep Eu anomaly zircons is that they signify particularly thin continental crust in central Australia during latest Mesoproterozoic time.

The fifth group consists of 1050–850 Ma zircons that overlap in age with 1100–1000 Ma Giles event bimodal magmatism in the Musgrave province (Fig. 2). Intrusive rocks of this stage have a range of Nd isotope compositions that are generally more radiogenic than those of the preceding Musgrave orogeny (Wade et al., 2008). In contrast, both Yb/U and Ce anomaly of ca. 1050–850 Ma detrital zircons are lower than those of the ca. 1250–1050 Ma group (Fig. 6). This seemingly contradictory relationship may be the product of a fundamental aspect of bimodal volcanism: radiogenic Nd isotope ratios of mafic rocks emplaced during, for example, extensional tectonism may indicate a mantle-derived component, but zircon from coeval felsic rocks produced by anatexis may nevertheless preserve a non-radiogenic composition. Interestingly, the Lu/Nd distribution of the ca. 1050–850 Ma zircons seems to be bimodal (Fig. 6). The sixth detrital zircon group (850–700 Ma) has a relatively limited population. This group overlaps in age with mafic dike swarms in the Musgrave region and mafic volcanic rocks in the Bitter Springs Group of the Amadeus Basin, and has similar trace element compositions as the preceding 1050–850 Ma group (Fig. 6).

A particularly important aspect of the new central Australia detrital zircon trace element record is the composition of the seventh and final group, ca. 700–500 Ma detrital zircons that are principally from Cambrian sandstone samples in the Warburton Basin (Fig. 4). Detrital zircons of this age range in Australia are commonly denoted the “Pacific-Gondwana” population (Ireland et al., 1998; Purdy et al., 2016). They are widespread in early Paleozoic siliciclastic strata across Australia, and, given the absence of clear significant source regions in Australia, they are the most notable “exotic” zircons in the Australian detrital zircon record (e.g., Fergusson et al., 2001; Maidment et al., 2007; Purdy et al., 2016; Kirkland et al., 2020). Although they are not well represented in the current Amadeus Basin detrital zircon Hf isotope record, they are extensively documented in detrital zircon Hf isotope records from other parts of Australia (e.g., Purdy et al., 2016; Martin et al., 2017; Morón et al., 2019; Kirkland et al., 2020). In these regions, the Pacific-Gondwana zircons have a strikingly wide range of εHf(i) from <−40 to +10. In the Warburton Basin samples, the trace element composition of this group is characterized by particularly low Yb/U, low Lu/Nd, low Ce anomaly, and generally greater Eu anomaly than the other groups (Fig. 6). As discussed below, this set of characteristics is arguably the most extraordinary aspect of the global zircon geochemistry record.

Global Zircon Trace Element Record

Motivated by some apparent correspondence of the Amadeus Basin Hf isotope record with the various central Australian zircon trace element records, we added our new Australian results to a global database of zircon trace element composition compiled from 56 published studies that includes >25,000 individual analyses (Fig. 7 and references therein; Table S3 and Fig. S2 [footnote 1]). It is informative that a significant fraction of these compiled data comes from detrital zircon studies in which zircon trace element compositions are reported in supplemental data tables but are nevertheless unmentioned in the actual papers. All REE parameters such as Ce and Eu anomaly were determined uniformly for each zircon analysis from original trace element concentration data tables according to the methods described above for our new central Australia data. Zircon U-Pb ages in the compiled data set range from <0.1 Ma (zircons from Iceland [Carley et al., 2011] and the Juan de Fuca Ridge in the northeastern Pacific [Schmitt et al., 2011]) to >4 Ga (Jack Hills zircons, Western Australia; Peck et al., 2001; Crowley et al., 2005; Bell et al., 2016; Turner et al., 2020). Like in other global zircon compilations, U-Pb ages are not uniformly distributed in the data set. Instead, major U-Pb age peaks occur at roughly 200 Ma, 1700 Ma, and 2700 Ma. In the case of the trace element compilation, the peaks are primarily a result of sampling bias. The major age populations that are overrepresented based on comparison with other zircon compilations (e.g., Campbell and Allen, 2008; Spencer, 2020) are Jurassic zircons from Alaska (Herriott et al., 2019) and <10 Ma zircons from eruptive products of the Yellowstone hotspot (Lu et al., 2016; Rivera et al., 2016), both in the western United States. To illustrate the variability of trace element parameters within the data set, the distributions of some important parameters and REE + U plots are shown in Figures 8 and 9. We divide discussion of this record into three broad time intervals, followed by a final note on three potentially important observations of the global zircon trace element record that span Earth history.

Late Paleoproterozoic and Mesoproterozoic

Owing primarily to the large amount of data from central Australia, South Australia, and western North America, one of the most densely sampled periods of zircon trace element geochemistry is the interval from ca. 2000 to 1000 Ma (Figs. 7 and 8). In this interval, zircon Yb/U, Lu/Nd, and Ce anomaly all follow the same general inverted U-shaped pattern as zircon εHf(i): relatively low values at ca. 1800 Ma (i.e., Aileron province ages in central Australia), a subsequent rise that reaches an apex at ca. 1400 Ma (partially captured by central Australia data, but best represented by detrital zircon from western North America), followed by a decline to low values at the close of the Mesoproterozoic Era (Musgrave province ages in central Australia). Low values of all four parameters in the ca. 2000–1000 Ma interval therefore seem to correspond with periods of supercontinent assembly, whereas elevated values generally correspond with the opening of new ocean basins during the breakup of Nuna. A potential explanation for this covariance is based on the zircon geochemistry framework discussed previously, namely that non-radiogenic εHf(i) values correspond with melting of old continental crust, low values of Yb/U indicate “continental” magmatism, low values of Ce anomaly correspond with relatively reduced magmas, and low values of Lu/Nd (i.e., relatively HREE-depleted compositions) arise from continentally derived magmatism and deep melting (Table 1). Low values of zircon εHf(i), Yb/U, and Lu/Nd are therefore predictable consequences of a significant involvement of continental material in zircon formation because of differences in elemental compatibility between continental crust versus mantle. Low Ce anomaly during supercontinent formation has a counterintuitive explanation, however, given that continental arcs are thought to be regions of mantle oxidation relative to mid-ocean ridges. The global zircon Ce anomaly data set suggests that potential subduction-induced oxidation reflected by zircon trace element composition is overwhelmed by melting of sediments to produce S-type granites that are characteristically reduced and which produce low Ce anomaly zircons (and non-radiogenic Hf isotope composition) during periods of supercontinent assembly. A related observation is that zircon from S-type granites tends to have less radiogenic Hf isotope composition (i.e., more negative εHf(i)) than zircon from I-type granites (e.g., Belousova et al., 2006).

Neoproterozoic to Cambrian

Some of the largest fluctuations in the zircon trace element record occurred during Neoproterozoic to Cambrian time (Figs. 7, 8, 10, and 11). In the late part of the Tonian Period (ca. 750 Ma), median zircon Yb/U rose to elevated values approaching 10, pointing toward juvenile mantle-derived magmatism accompanying Rodinia rifting. With the exception of some earliest Cryogenian zircons from northwestern Canada (Macdonald et al., 2018) that are discussed below, zircons with Cryogenian U-Pb ages are relatively sparse in the compiled data set, but the oldest cluster of Ediacaran zircons (ca. 630 Ma) has Yb/U near the overall median. Similar to global zircon Hf isotope compilations (e.g., Belousova et al., 2010; Roberts and Spencer, 2014), the record of zircon Yb/U (Figs. 7, 8, and 10) reaches a low during late Neoproterozoic to Cambrian time (i.e., coeval with “Pacific-Gondwana” zircons of Australia). The underlying cause of this zircon geochemical anomaly was reworking of continental crust (e.g., Spencer et al., 2014), although different mechanisms and varieties of reworking have been suggested. For example, the nadir of the anomaly is at ca. 580 Ma, ∼55 m.y. after the termination of the global glaciations that marked the Cryogenian Period (e.g., Hoffman et al., 1998). The aftermath of the “snowball Earth” period would have involved massive weathering of the continents and related contribution of continental material to the oceans (e.g., Spencer et al., 2013; Hoffman et al., 2017; Keller et al., 2019; Ward et al., 2019). One explanation for the subsequent relatively short-duration zircon geochemical anomaly is that it signifies subduction of this continentally derived oceanic sediment (Keller et al., 2019). Continental sediment could have also been sourced from weathering and erosion of mountains related to the collision of East and West Gondwana (Squire et al., 2006). Alternatively, the late Neoproterozoic to Cambrian zircon anomaly might be related to extensive passive-margin incorporation in collisional zones during Gondwana assembly (Roberts and Spencer, 2014). A potentially important consideration in evaluating the zircon geochemistry anomaly is that detrital rutile U-Pb cooling ages in North Gondwana are also ca. 580 Ma (Avigad et al., 2017). The overlap of ages suggests that crystallization of the anomalous zircons corresponded with a widespread period of exhumation and cooling, and, in fact, the dominant style of tectonism at ca. 580 Ma in the East African orogen was extensional (e.g., Fritz et al., 2013). In the extant zircon trace element record, the ca. 700–500 Ma zircon Yb/U anomaly is evident in data sets from Antarctica (Takamura et al., 2020), Africa (primarily in detrital zircon data sets from the Congo, Niger, and Nile Rivers; Iizuka et al., 2013; Zhu et al., 2020), Nepal (Neupane et al., 2020), New Caledonia (Campbell et al., 2018, 2020b), New Zealand (Campbell et al., 2020a), South America (primarily detrital zircons from the Rio São Francisco; Zhu et al., 2020), as well as southern (Belousova et al., 2009; Keeman et al., 2020), western (Morón et al., 2019), central (this study), and eastern Australia (Manchee et al., 2019; Zhu et al., 2020).

Over the ∼200 m.y. period spanning from late Tonian to mid-Ediacaran time, there was thus a particularly significant decline in zircon Yb/U, with the median declining by nearly an order of magnitude and outlying values declining by more than two orders of magnitude. Over the same interval, zircon εHf(i) and Lu/Nd also declined (Fig. 10). Not coincidentally, the Nd isotope composition of mudstone decreased over the same interval (Cox et al., 2016). This pattern is probably attributable to the relatively rapid succession of late Tonian mantle-derived magmatism (which was “juvenile” in terms of εHf(i), “oceanic” or “MORB-like” in terms of Yb/U and Lu/Nd, and oxidized or “I-type” in terms of Ce anomaly) and subsequent mid-Ediacaran to Cambrian magmatism during incorporation of old continental crust (“evolved,” “continental,” and reduced or “S-type” by the same metrics) via a variety of potential mechanisms that broadly coincided with Gondwana assembly.

Over the Neoproterozoic to Cambrian interval, the trends of zircon Ce anomaly and Eu anomaly are similar to each other, but different than the trends of εHf(i), Yb/U, and Lu/Nd (Fig. 10). Both Ce and Eu anomalies reach low points at around the Tonian-Cryogenian boundary, and Eu anomaly has an additional low at ca. 850 Ma (Fig. 10). The record of zircon Eu anomaly is similar to that of the long-term trend of Eu anomaly from felsic rocks in that both reach minima in the early to mid-Neoproterozoic (Keller and Schoene, 2012). The Neoproterozoic Eu anomaly low of zircons and felsic rocks may reflect a period when there was relatively little deep crustal melting at pressures greater than that of plagioclase stability, a paucity that could have been related to thinned continental crust during the early stages of Rodinia rifting. In the compiled data set, the major constituent of earliest Cryogenian (ca. 715 Ma) zircons are those from volcanic rocks in northwestern Canada (Macdonald et al., 2018) that may be a distal part of the Franklin large igneous province (Cox et al., 2012). These zircons have relatively high Yb/U and Lu/Nd (i.e., mantle-like), as well as low Ce anomaly and Eu anomaly (Figs. 10 and 11). Their compositions thus seem to point toward a low-oxygen-fugacity mantle source, perhaps a deep mantle plume.

Phanerozoic

In the compiled zircon data set, Yb/U is particularly low at ca. 340 Ma (Mississippian; Figs. 7 and 8). This period is most extensively represented in the data set by Variscan-age zircons from gneiss in central Europe (Kylander-Clark et al., 2013; Kotková et al., 2016; Kozlik et al., 2016) as well as in the European river detrital zircon data set of Zhu et al. (2020). In addition to low Yb/U, the same zircons in the studies of Kotková et al. (2016) and Kozlik et al. (2016) have non-radiogenic εHf(i) of approximately −5 to −10, relatively low Th/U <0.2, and very high U concentration of as much as ∼30,000 ppm, and they seem to have formed during ultrahigh-temperature and ultrahigh-pressure metamorphic conditions. This time period is also the interval of Earth history richest in S-type granites (Zhu et al., 2020). This component of the zircon trace element record could thus be interpreted as metamorphism and incorporation of continental crust during the formation of Pangea. In contrast, zircon Yb/U reaches particularly large values (the greatest Phanerozoic values in the compiled data set) at ca. 240 Ma (Middle Triassic) in a detrital zircon population from Iran (Figs. 7 and 8; Barber et al., 2019). Given that mid-Triassic rifting in Iran was likely related to the fragmentation of Gondwana (Ramezani and Tucker, 2003; Barber et al., 2019), one could view the large-magnitude variation of zircon Yb/U over the ∼100 m.y. period from early Carboniferous to mid-Triassic time (Fig. 8) as fundamentally related to Pangea assembly and subsequent breakup.

General Features

We conclude by noting three intriguing features of the time-integrated global zircon trace element data set not discussed previously. The first is that high (i.e., relatively flat) zircon Eu anomaly seems to be a characteristic of Cu porphyry mineralization (e.g., Ballard et al., 2002; Dilles et al., 2015; Loader et al., 2017; Lu et al., 2016, 2019). The secular record of zircon Eu anomaly (Fig. 8) might thus be viewed as a proxy of Cu porphyry formation, which, in turn, could be seen as an indicator of subduction zones and some form of modern-style plate tectonics (e.g., Barley, 1982; Cawood et al., 2006). The ca. 3.4–3.3 Ga transition and ca. 3.2–3.1 Ga peak in the zircon Eu anomaly record (Fig. 8) may be significant in that context, as well as in delineating subdivisions of the Archean Eon (e.g., Pease et al., 2008; Griffin et al., 2014; Brenner et al., 2020).

The second observation of the overall data set is that zircon crystallization temperatures estimated from Ti-in-zircon thermometry (e.g., Watson et al., 2006) seem to form a mirror image pattern with zircon Eu anomaly after ca. 3 Ga (Fig. 12). In other words, the long-term secular pattern appears to be that zircon crystallization temperature increases as zircon Eu decreases (i.e., as Eu anomalies become deeper). Given that zircon Eu anomaly is largely controlled by plagioclase fractionation and, in turn, crystallization depth and crustal thickness (e.g., McKenzie et al., 2018), the compiled data set suggests that over the last ∼3 b.y., the main factor controlling zircon crystallization temperature is the thickness of the continental crust, such that zircons produced from thin continental crust (deep Eu anomaly) have higher crystallization temperatures than those produced from thick crust (flat Eu anomaly). A logical explanation for this apparent correspondence is that heat advection from the mantle through the continental crust (i.e., the geotherm) has been the main long-term control on zircon crystallization temperature since approximately the Mesoarchean Era.

Finally, in the global zircon trace element data set, fluctuations of Yb/U seem to be larger in amplitude with decreasing age (Fig. 7), the same pattern as for zircon εHf(i) (Fig. 3). Based on this observation, we infer that Yb/U of the DM has changed over Earth history. Yb/U of C1 chondrite is ∼20 (Sun and McDonough, 1989), whereas Yb/U of the present-day DM is estimated at ∼85 (Salters and Stracke, 2004). DM enrichment of Yb relative to U over time seems to be entirely analogous to the DM enrichment in εHf: extraction of U from the mantle to form continental crust has progressively depleted the mantle in U relative to Yb, leading to a present-day DM U/Yb that is <25% of the chondritic ratio.

CONCLUSIONS

The bulk composition of igneous rocks has changed over Earth history (e.g., Keller and Schoene, 2012), and igneous rock geochemistry is reflected in the trace element composition of zircon (Belousova et al., 2002). It therefore stands to reason that bulk zircon trace element composition has changed over Earth history (Balica et al., 2020). Secular variations of some key zircon trace element parameters such as Yb/U, Ce anomaly, and LREE depletion relative to HREEs seem to follow the same general patterns as zircon Hf isotope composition, a correspondence suggesting that the major control on zircon trace element parameters is the balance between inputs of mantle and continental crust. A key time period of Earth history for evaluating these relationships is a late Neoproterozoic to Cambrian interval, when zircon Hf isotope composition became particularly non-radiogenic. Zircons from this age interval span the entirety of εHf(i) space, from that of DM to that of zircons derived from very ancient crust. These zircons also tend to have low Yb/U, low Ce anomaly (a nominal indication of relatively “reduced” magmatic conditions), and low Lu/Nd, a simple metric of HREE enrichment. Based on examination of central Australian and global zircon trace element data sets, we conclude that detrital zircon trace element records are potentially useful in conjunction with isotopic records to investigate tectonomagmatic histories.

ACKNOWLEDGMENTS

We thank Amos Pellas, Nathan Siddle, and Jack Ward for their assistance with this project. Jay Carter, Peter Haines, and Max Heckenberg were instrumental in providing access to drill core in the Northern Territory and Western Australia. Karine Harumi Moromizato was a great help in the lab. The paper benefited from reviews by Brenhin Keller and two anonymous reviewers. This project was supported by Australian Research Council grant DE140100553.

1Supplemental Material. Table S1: Location information for new detrital zircon samples. Table S2: New zircon U-Pb and trace element data. Table S3: Compiled zircon U-Pb and trace element data. Figure S1: Comparison of Ce anomaly (A) and Eu anomaly (B) of zircons calculated using the method of Onuma et al. (1968) and traditional methods. Figure S2: Summary figures from zircon trace element compilation. Please visit https://doi.org/10.1130/GEOS.S.13389110 to access the supplemental material, and contact editing@geosociety.org with any questions. Additional material for Table S2 is located at Github: https://github.com/cverdel/zircon_petrochronology.
Science Editor: Shanaka de Silva
Associate Editor: Christopher J. Spencer
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.