Spot the difference: Zircon disparity tracks crustal evolution

Despite application of an increasing array of geochemical tools, facets of the evolution of Earth’s continental crust remain enigmatic, most significantly because of denudation and subduction recycling of older crust (Condie, 2018; Hawkesworth et al., 2018). Zircon crystals grow in a range of magmatic rocks, are physicochemically robust, and accumulate as common accessory phases in sedimentary detritus. The stability of the ZrSiO₄ lattice and slow element diffusion enable primary age and geochemical data to be preserved in zircon grains over geological time scales, through weathering, erosion, diagenesis, and metamorphism, permitting the tracking of crustal evolution from at least 4.4 Ga (Hawkesworth and Kemp, 2006; Spencer et al., 2014). Despite source rock fertility influences and transport modification, investigation of detrital zircon populations captures snapshots of a drainage system’s geological history and reduces preservation and sampling bias inherent in the available crystalline rock record (Hawkesworth et al., 2018).

While the age of “big data” is shedding new light on the development of continental crust (e.g., Voice et al., 2011), statistical tools are needed to resolve meaningful trends in data sets. Here, we demonstrate the potential of a new detrital zircon–based technique—temporally binned Kolmogorov-Smirnov distance-based multidimensional scaling (Vermeesch, 2013)—to assess geological relationships between discrete crustal regions across Earth history. We propose that this approach to global data provides novel insight into supercontinent development and fundamental tectonic processes that couple deep (magmatic-generational) and surface (sedimentary-preservational) systems.

Detrital zircon analyses from a recently created Greenland geochronology database (Kirkland et al., 2016) were integrated with those extracted from the study of Voice et al. (2011) to yield 149,767 detrital zircon analyses spanning 4–0 Ga (Table DR1 in the GSA Data Repository¹). Analyses were filtered to include only those within 10% of concordia. Analyses older than 4 Ga (298) are rare and geographically biased to Western Australia, and so they were omitted from global comparison. Data were divided geographically according to their sampling site and plotted as kernel density estimation age spectra (Fig. 1; Vermeesch, 2012). For grains older than 1500 Ma, ²⁰⁷Pb/²⁰⁶Pb ages were used, and ²³⁸U/²⁰⁶Pb ages were used for younger grains. Data from each geographic region were binned at 200 m.y. intervals to increase data per bin and thereby reduce the impact of local data significance, i.e., noise, while also resolving trends on a sub–500–300 m.y. supercontinent periodicity. Bins for any region with less than 30 analyses were deemed insufficient to meaningfully capture major peak age characteristics (Vermeesch, 2004).

We measured the similarity of detrital zircon spectra between different regions for each 200 m.y. time slice using the Kolmogorov-Smirnov distance statistic (Massey, 1951) expressed in multidimensional scaling (MDS) space (Vermeesch, 2013, 2018; see the Data Repository for a more detailed account of the methods). MDS diagrams graphically represent the similarity of samples by their proximity on a plot (Fig. 1B). We applied cluster analysis within MDS space to determine changes in the degree of similarity in detrital zircon populations through time. To avoid bias in measures of clustering, both average distance to a centroid (CD) and average nearest neighbor (NN) values were determined in MS Excel^® and PAST (https://folk.uio.no/ohammer/past/), respectively (Fig. 1B; Hammer et al., 2001). The CD metric focuses on the proximity (similarity) of all regions relative to each other, while the NN metric is more sensitive to subclusters of regions (connections between some regions regardless of similarity to the total data set).

The detrital zircon time-disparity pattern was assessed for correlation with the supercontinent cycle via a point-biserial test in SPSS, using disparity measures as a continuous record against a simplified, binary supercontinent/no-supercontinent record (Kornbrot, 2014). There are no qualitative deep-time measures of the degree of supercontinentality. Instead, we adopted the definition of a “supercontinent” as a substantial grouping of formerly dispersed continents (Bradley, 2011). The most recent constraints on the timings of supercontinents were compiled for the Phanerozoic and Proterozoic (Table DR2), including relevant constructive phases of the supercontinents. We calculated r and p values to establish the sense and statistical significance of any relationship between the variables (see the Data Repositoryfor further explanation).

Detrital zircon age spectra exhibit broad similarities in age components and define global age peaks at ca. 0.5 Ga, ca. 1.0 Ga, ca. 1.8 Ga, ca. 2.5 Ga, and ca. 2.7 Ga, which correspond to intervals of supercontinentality (Fig. 1A; Hawkesworth et al., 2018). A distinct contribution of detrital zircon ages outside of a supercontinent episode from the North and Central American region provides a “Western Cordillera” (<0.25 Ga) sampling bias in the global data set from ongoing Pacific subduction. Disparity data vary by close to one order of magnitude, with values reflecting a relative metric within the time series being considered (Fig. 2; Table DR3). Overall, global detrital zircon disparity values are low through the Eoarchean to Mesoarchean, rise through the Neoarchean, and fall again through the first half of the Paleoproterozoic. Global detrital zircon disparity measures oscillate through the rest of the Proterozoic and Phanerozoic, with peaks in disparity centered at ca. 1.5 Ga, ca. 0.9 Ga, ca. 0.5 Ga, and 0.1 Ga (Fig. 2). Temporal variations in global detrital zircon disparity measure are broadly consistent regardless of the binning interval used (400, 200, or 100 m.y.), with reduced data per bin more likely to express noise from local signals in shorter intervals, while longer bins lack sufficient resolution for geologically significant processes such as the supercontinent cycle (Nance et al., 2014).

Ideally, when interpreting global detrital zircon age data, consideration should be given to the potential of multicycle detritus, as well as renewed/enhanced erosion of preexisting, regionally specific, age-diagnostic crystalline rocks during later continent associations. Theoretically, interregional detrital zircon homogenization via exchange of preexisting diagnostic detrital zircon age signatures during supercontinent assembly has the potential to falsely indicate older, shared geological histories. However, Parman (2015) demonstrated that although distinctive age peaks in the global detrital zircon record are persistent through time, in general, diagnostic peaks are most prominent in preserved sediments laid down soonest after source crystallization. Furthermore, areas currently juxtaposed, or areas that were adjacent during recent supercontinent events (e.g., Antarctica, Australia, and India) show distinct total population detrital zircon spectra (Fig. 1). These observations support the dominantly intracrustal block erosion-preservation coupling that accounts for characteristic detrital zircon age peaks within geographically pooled data sets, rather than secondary exchange during later collisions. This, in turn, implies that the global detrital zircon database predominantly reflects fundamental tectono-magmatic production and preservation associations, and that detrital zircon global disparity patterns provide a tool that can track first-order crustal assemblages and processes through time.

Detrital zircon data were grouped according to present-day geographical sampling locations, reflecting one approach to areal pooling that mitigates limitations of insufficient data from smaller regions (impact of local/sampling bias), as well as complexities of tracking modern sampling back through complex geological histories for dissected or amalgamated crustal elements. While the geographic regions utilized do not necessarily reflect defined crustal blocks, many do contain discrete cratonic cores. Therefore, the approach provides a spatially unitized reference that can provide meaningful insights regarding geologically influenced disparity variations within the global detrital zircon population.

Disparity measures of geographically distinct detrital zircon age populations, based on both NN and CD metrics, largely trace the same patterns of similarity through geological time (r value = 0.8, p value << 0.000), except during the 2.6–2.4 Ga interval. During this period, the NN metric tracks two relatively distinct subclusters of regional detrital zircon similarity and consequently indicates lower disparity, while the CD value is insensitive to subclusters and indicates higher disparity (Fig. 1B). The components of the two clusters identified by NN analysis within the 2.6–2.4 Ga interval broadly correspond with proposals for Kenorland and Nunavutia (Pehrsson et al., 2013), during one of the earliest known intervals of craton assembly (Evans, 2013). These observations, and the coherence of detrital zircon disparity patterns through time across different temporal bin widths, strongly support the interpretation of a real geological signal within global detrital zircon disparity measures.

Myriad studies have previously identified patterns in the relative abundance of zircon crystallization ages, as well as the Hf- and O-isotope records preserved in the same grains, to help constrain crustal evolution during supercontinent assembly and breakup (Fig. 3). Broad temporal agreement across independent isotopic systems (e.g., O, Hf, Sm-Nd; Condie, 2014; Gardiner et al., 2016), crustal geochronology, and paleomagnetism (Veikkolainen et al., 2017) has improved confidence on the timing of Earth’s supercontinent cycle. Whether some of these proxies track fundamental changes in crust generation and destruction (Condie, 2005) or preservational bias (Hawkesworth et al., 2009) associated with continental assembly and breakup continues to be debated. Regardless of which model is correct, here we demonstrate the excellent correlation between detrital zircon disparity measures and the perceived supercontinent cycle (p values 0.009–0.003; Table DR2). Sensitivity analysis of the disparity measure indicates the metric is robust to spurious signals (see the Data Repository).

Supercontinents involve shared orogenic belts and enhanced sediment dispersal and ultimate preservation from these same uplifted crystalline source regions. Conversely, separated continental fragments are characterized by unique magmatic events and independent sediment systems. Thus, overall, global detrital zircon age spectra are more similar (lower disparity) during supercontinent assembly and less similar (higher disparity), and regionally unique, during separation (Figs. 3 and 4).

An increase in global detrital zircon disparity occurs in the latest Neoarchean, before returning to low values during the 2.0–1.6 Ga interval, correlated with Nuna formation (Figs. 2 and 3). This pronounced excursion is interpreted to represent significant changes in crustal preservation and recycling during the transition from the early Earth’s relatively hot crust and mantle conditions (Fig. 4).

Increases in measured detrital zircon disparity from ca. 2.8 Ga are interpreted to reflect increasingly emergent, isolated crustal components being progressively incorporated into the preserved rock record (Figs. 2–4; Flament et al., 2008). Increased zircon abundance is recognized globally, along with an interpreted change from dominantly destructive crustal processes to production processes (Parman, 2015). These interpretations are supported by the appearance of passive continental margins in the geological record (Bradley, 2011), and geochemical evidence for significant pulses of juvenile continental crust growth between 3.0 and 2.6 Ga (Cawood et al., 2009), as well as increasing preservation and abundance of discrete continental crust fragments that formed in the Neoarchean.

Although geochemical evidence for juvenile crust addition and subduction processes has been identified in rocks from ca. 3.2 Ga (Van Kranendonk and Kirkland, 2016), it is likely that the transition to a world dominated by plate tectonics was gradual (Condie, 2018). The establishment of widespread subduction was associated with stiffening, thickening, and strengthening of the continental lithosphere from 2.6 to 2.4 Ga (Van Kranendonk and Kirkland, 2016), as also evidenced by the geochemistry of mafic rocks that reflect cooler mantle from the introduction of cold slabs (Condie, 2018). Coupling of zircon O- and Hf-isotope curves becomes apparent from approximately the Neoarchean, when isotope records become synchronized and more variable in response to increasing crust-mantle differentiation and later crustal reworking during the supercontinent cycle (Fig. 3). Zircon O-isotopes exhibit more positive values (with a distinctly different mean oscillation level after the Archean), corresponding to an increasing degree of crustal reworking (Fig. 3; Valley et al., 2005). This transition to widespread subduction and crustal recycling (Spencer et al., 2014; Roberts and Spencer, 2015), and consequent assembly of crustal fragments, is supported by a number of independent indicators of orogenesis, such as metamorphism (Brown, 2014) and orogenic gold mineralization (Goldfarb et al., 2010). Thus, reduced global detrital zircon disparity in the early Paleoproterozoic tracks increasing crustal recycling and orogenesis, and reflects the fundamental shift to a world with supercontinents, coincident with an increasing influence of collisional tectonics (Fig. 4).

Time-series disparity measures of global detrital zircon age spectra provide a new technique with which to investigate a number of geologically significant lithospheric to surface processes. Fundamentally, higher global detrital zircon disparity correlates with regionally unique age records driven by (1) more widespread crustal preservation of discrete crustal fragments (Neoarchean), or (2) separate evolution of disconnected continental areas (Proterozoic–Phanerozoic). Conversely, lower global detrital zircon disparity correlates with matching age records dictated by (1) global signatures on limited crustal blocks (Archean), and (2) increasing amalgamation of continental areas via subduction, and the similarity of preserved detritus from shared tectono-magmatic events during continent assembly (Proterozoic–Phanerozoic). Zircon disparity data are sensitive to crust generation and destruction both prior to plate tectonics and through the supercontinent cycle, and they illustrate a fundamental coupling of magmatic crust-formation and sedimentation processes. Furthermore, disparity time-series analysis holds significant promise for more spatially or temporally focused applications to assess evolving catchment-basin-craton relationships.

We thank D. Brown, D. Bradley, F. Corfu, C. Hawkesworth, and two anonymous reviewers for comments that improved this manuscript.