Plumbing the depths of magma crystallization using ¹⁷⁶Lu/¹⁷⁷Hf in zircon as a pressure proxy

Magmas rise within the continental crust until they stall and solidify (Brown, 1994; Petford et al., 2000; Annen et al., 2006). Their crystallization tends to be shallower in extensional settings than in compressional settings due to contrasting normal stress gradients and conditions of magma propagation (Hogan and Gilbert, 1995; Vigneresse et al., 1999; Loucks, 2021). In geologically young terrains, the crystallization depth of magmas can be measured through geophysical approaches (e.g., Thybo and Nielsen, 2009). However, as a result of erosion and reworking of the crust through time, much of the ancient geologic record has been destroyed. As a consequence, there is little information with which to constrain the depths of crystallization and the periods of crustal extension and compression in the deep geologic past.

Zircon (ZrSiO₄) crystallizes from a range of SiO₂-rich magmas, and it is resistant to weathering and erosion processes, such that it is widely preserved in continental sediments across geologic time (Belousova et al., 2010; Voice et al., 2011). Many studies have explored the zircon archive to draw insights into Earth’s geodynamic regimes over time (e.g., Belousova et al., 2010; Voice et al., 2011; Dhuime et al., 2012; McKenzie et al., 2018; Balica et al., 2020). These include the use of U-Pb ages to infer the relative volume of magmatism (e.g., Schoene et al., 2012), Hf isotopes as a measure of mantle input and crustal reworking/recycling (e.g., Gardiner et al., 2016), trace-element ratios to infer mean crustal thickness (e.g., Tang et al., 2021), and oxygen isotopes as a proxy for crustal reworking (e.g., Dhuime et al., 2012; Spencer et al., 2014).

Variations in zircon ¹⁷⁶Lu/¹⁷⁷Hf have been documented to describe changes in magma conditions and metamorphic reactions at local scales (e.g., O’Brien and Miller, 2014; Rubatto, 2017), but global data sets of ¹⁷⁶Lu/¹⁷⁷Hf in zircon remain unexplored. Considering that variations in zircon composition may reflect the minerals crystallizing from the host magmas and residues during partial melting, we evaluated the relationship between zircon ¹⁷⁶Lu/¹⁷⁷Hf and pressure. We show that the ¹⁷⁶Lu/¹⁷⁷Hf ratio can be used to ascertain the depth of crystallization of zircon-bearing magmas, and we explore the implications for large-scale tectonic variations through time.

Figure 1 shows the variation of the geometric mean of ¹⁷⁶Lu/¹⁷⁷Hf in zircons from a worldwide data set of >120,000 analyses, subdivided into 0.1 b.y. time slices (Data Set S1 in the Supplemental Material¹). The variation of the zircon ¹⁸O/¹⁶O ratios (expressed as δ¹⁸O relative to Vienna standard mean ocean water [VSMOW]) from an independent global database (Spencer et al., 2014) is plotted for comparison. It is striking that major troughs in ¹⁷⁶Lu/¹⁷⁷Hf correlate with peaks in δ¹⁸O, and these in turn broadly coincide with periods of supercontinent assembly. This observation was confirmed by the calculated anomaly in the running means of 0.1 b.y. intervals (inset Fig. 1), with a 0.5–0.6 b.y. cyclicity in δ¹⁸O peaks and ¹⁷⁶Lu/¹⁷⁷Hf troughs at 2.75–2.45 Ga, 2.1–1.9 Ga, 1.3–1.05 Ga, and 0.75–0.55 Ga (as indicated by coherent statistical variations between the two data sets; Tables S1–S5).

Higher values of δ¹⁸O in zircon can be linked to increased assimilation/reworking of supracrustal components and crustal thickening at times of major continental collisions (Spencer et al., 2014). The concomitant lower mean values of ¹⁷⁶Lu/¹⁷⁷Hf during periods of higher δ¹⁸O may similarly reflect crustal thickening and collisional tectonics.

Horizontal compression inhibits the ascent of magmas within the crust (Loucks, 2021), and crustal thickening expands the garnet stability field (e.g., Taylor et al., 2016; Rubatto, 2017). Given that Lu partitions preferentially into garnet, crustal thickening and an expanded garnet inventory in the crust may be influential factors. Zircon crystallization increases ¹⁷⁶Lu/¹⁷⁷Hf in the coexisting melts, and more fractionated magmas tend to have higher ¹⁷⁶Lu/¹⁷⁷Hf values. Thus, the ¹⁷⁶Lu/¹⁷⁷Hf of magmas and their zircons may be influenced by greater depths/pressures of crystallization in orogenic zones, and by the degree of fractionation as magmas ascend within the crust (Fig. 2A).

To assess the controls on ¹⁷⁶Lu/¹⁷⁷Hf in zircon from igneous magmas, we compiled studies with pressure estimates of zircon-bearing igneous rocks from which zircons were measured for ¹⁷⁶Lu/¹⁷⁷Hf isotopes. These studies generally relied on Al-in-hornblende to determine pressures of crystallization, while other petrologic/geochemical approaches were used to define pressures of magma generation (e.g., experimental petrology, thermodynamic modeling).

Zircon saturation-crystallization may occur from depths of magma generation to depths of final emplacement (Schoene et al., 2012), depending on the composition, temperature, and alkali content of the host crystal-rich magma reservoir (Boehnke et al., 2013; Laurent et al., 2020). At present, it is not possible to directly measure pressures of crystallization in zircon, but these can be estimated by analyzing mineral phases interpreted to be cogenetic with zircon crystallization (e.g., Barnes et al., 2019). Two or more minerals from a single paragenesis can show contrasting pressures of crystallization, and the reasons are often debated (e.g., Walker et al., 2013). To overcome this issue, our data set included a range of pressures from partial melting to emplacement depth estimates to better represent the complexity of transcrustal magmatic systems (Fig. 2A; Cashman et al., 2017).

The data spanned near-surface depths to ~2.0 GPa pressures and ranged in ¹⁷⁶Lu/¹⁷⁷Hf from 0.0005 (ultrahigh-pressure tonalites) to >0.002 (volcanic zircons). An exponential relationship of decreasing ¹⁷⁶Lu/¹⁷⁷Hf with increasing pressure is described in Equation 1:

where GPa_(sample) is the pressure determined for crystallization, and ¹⁷⁶Lu/¹⁷⁷Hf_{(igneous zircon)} is the geometric mean of many determinations for a sample or group of related samples from the same geologic unit. For an individual sample, the ¹⁷⁶Lu/¹⁷⁷Hf ratios can vary by ±25% in cogenetic zircons, whether volcanic or plutonic. As a result, ~10 or more analyses are required to estimate the mean ¹⁷⁶Lu/¹⁷⁷Hf and associated pressure for a given data point (Fig. 2B). One important caveat is therefore that Equation 1 is not appropriate for single-zircon pressure determinations.

Metamorphic zircons represent a modest component in the global geologic record (e.g., ~8% estimated by Balica et al., 2020). Similar to the approach for igneous zircons, we compiled data for metamorphic zircons using their age, ¹⁷⁶Lu/¹⁷⁷Hf, and phase equilibria–determined pressure-temperature conditions (Fig. 2B). Igneous zircons with ¹⁷⁶Lu/¹⁷⁷Hf < 0.0003 are scarce, and a large majority of metamorphic zircons have ¹⁷⁶Lu/¹⁷⁷Hf < 0.0005 (Fig. 2C). The overlap between igneous and metamorphic zircons is therefore small, and the ¹⁷⁶Lu/¹⁷⁷Hf ratio can be used for screening zircons of unknown origin when internal textures are ambiguous or unavailable (see the Supplemental Material).

The main sink of Hf in igneous rocks is zircon, whereas Lu is incorporated in a number of other minerals aside from garnet (e.g., amphibole, pyroxene, apatite, titanite). Lu is less compatible in zircon than Hf, and so the ¹⁷⁶Lu/¹⁷⁷Hf ratio is substantially lower in zircon than in the melt from which it crystallized, approximately by a factor of 5–30 (Data Set S2). Garnet has a much greater affinity for Lu, it can be abundant in the source of magmas in the middle and deep crust, and any melt generated in equilibrium with garnet is likely to have a lower ¹⁷⁶Lu/¹⁷⁷Hf ratio than the residue (Taylor et al., 2016; Rubatto, 2017; Gardiner et al., 2018). Other minerals such as amphibole, pyroxene, apatite, titanite, and ilmenite can influence ¹⁷⁶Lu/¹⁷⁷Hf in the melt, but in the presence of zircon and/or garnet, they have a minor effect on ¹⁷⁶Lu/¹⁷⁷Hf in the evolving magma or in zircon (see the Supplemental Material).

The ¹⁷⁶Lu/¹⁷⁷Hf ratio of a melt is sensitive to the melt fraction and to the proportion of garnet in the residue (Fig. 3A). Once this melt forms and starts to ascend within the crust, zircon saturation and fractionation become the dominant controls on ¹⁷⁶Lu/¹⁷⁷Hf in the evolving melt (assuming the lack of igneous garnet influence on the liquidus), by subtracting Hf from the system (Fig. 3B). The variation of ¹⁷⁶Lu/¹⁷⁷Hf with pressure is therefore thought to reflect higher amounts of garnet in the source at greater depths, combined with the increase in magma fractionation with decreasing depth (Dhuime et al., 2015). Other possible factors controlling Lu/Hf variation are temperature (Claiborne et al., 2017) and potential changes in magma types through time (Laurent et al., 2014). The impact of these two factors on the ¹⁷⁶Lu/¹⁷⁷Hf ratio in zircon is minor (see the Supplemental Material), and we conclude that garnet as a residual phase and zircon fractionation and increasing magma fractionation with decreasing depth are the major factors that control ¹⁷⁶Lu/¹⁷⁷Hf variations with depth (Fig. 2). Crystallization is the dominant control on ¹⁷⁶Lu/¹⁷⁷Hf through most of the depth range considered (Fig. 2B), and, as such, this is the focus of the following discussion.

We applied Equation 1 to geometric means of ¹⁷⁶Lu/¹⁷⁷Hf from the global zircon database to calculate the apparent mean depths of crystallization of zircon-bearing igneous rocks through time (Fig. 4, red curve). The deepest mean depths of crystallization (~30–35 km) occurred at ca. 2.6 and ca. 1.9 Ga. Shallower depths (~10–15 km) are observed for the Hadean and Paleoarchean/Mesoarchean, suggesting that changes in the depths of crystallization are not simply a function of time.

Strikingly, a marked cyclicity was observed over the last 3.0 b.y., but it is absent in the earlier record (Fig. 4). The broad coincidence of troughs in the crystallization depth curve (yellow stars, associated with the lowest ¹⁷⁶Lu/¹⁷⁷Hf ratios) with peaks in δ¹⁸O and periods of supercontinent assembly (Fig. 1) suggests a firstorder link with collisional processes, such as thickening and reworking of continental crust. This pattern implies either deeper exhumation through erosional and tectonic processes associated with continental collision and thickening (Blackburn et al., 2018), more efficient entrapment of magmas at deeper levels during global compressional episodes (e.g., Loucks, 2021), or a combination of both processes. Reciprocally, the peaks in the crystallization depth curve (green stars, associated with the highest ¹⁷⁶Lu/¹⁷⁷Hf ratios) reflect periods of maximum shallowing during tectonic extension, as in supercontinent breakup.

An unexpected feature of Figure 4 is that both the depths of deepest crystallization (yellow stars) and the peaks of shallowest crystallization (green stars) have progressively decreased since the Neoarchean, with slopes of ~0.015 and ~0.010 km/m.y., respectively. This may be interpreted as follows: (1) The red curve is a primary signal for mean crystallization depths, with deeper magmas being more abundant during older collisional events; or (2) the progressive apparent decrease of crystallization depths is due to preferential removal of uppercrustal zircons from the global zircon archive via erosion and physical destruction. The latter hypothesis does not necessarily involve changes in the primary mean depths of crystallization.

The primary mean depth of crystallization has decreased with time since 3 Ga (Fig. 4), and this may be linked to the secular cooling of the mantle (Herzberg et al., 2010) or to changes in average crustal thickness (Tang et al., 2021), although these are complex issues to address with little consensus (e.g., Balica et al., 2020; Tamblyn et al., 2022). In contrast, cumulative erosion is known to increase with time, and in young orogens in which the upper crust remains well preserved, deeper orogenic rocks are present in smaller volumes than in older orogens (e.g., Blackburn, et al., 2018). For example, the comparison between Tibet and older orogens with more deeply eroded roots (e.g., ca. 0.6 Ga Braziliano, ca. 1.1 Ga Grenvillian, and ca. 1.9 Ga Trans-Hudson collisional orogens) revealed a much greater exposure of middle to deep crust in the older orogens (e.g., Weller et al., 2021). In addition, continental collision involves significant reworking and exhumation of deeper crustal levels that may contain ancient reworked basement and continental margin sediments. This results in a progressive removal of upper-crust zircons and exposure of the grains that crystallized at greater depths.

It is likely that the ¹⁷⁶Lu/¹⁷⁷Hf global zircon record has been affected by both erosion and geothermal lithospheric conditions over billions of years, and a better understanding of the observed shallowing of maximum magma crystallization depths through Earth history may involve more detailed modeling of the interaction between erosion and thermal history through the cyclical generation and destruction of supercontinents. Nonetheless, it is striking that the cycles of magma crystallization depths started in response to global supercontinent activity between 3.0 and 2.5 Ga, which may be further evidence that plate tectonics became global at that time.

This project was funded by the European Research Council under the European Union’s Horizon 2020 research and innovation program (no. 817934). We thank the anonymous reviewers, Steven L. Goldstein, and Oscar Laurent for constructive and insightful comments that allowed us to clarify the manuscript. We also thank Urs Schaltegger for very helpful editorial guidance.