Zircon petrochronology in large igneous provinces reveals upper crustal contamination processes: new U–Pb ages, Hf and O isotopes, and trace elements from the Central Atlantic magmatic province (CAMP)

Abstract

Zircon occasionally crystallizes in evolved melt pockets in mafic large igneous province (LIP) magmas, and in these cases, it is used to provide high-precision age constraints on LIP events. The precision and accuracy of high-precision ages from LIPs are crucially important, because they may be implicated in mass extinctions. However, why zircon crystallizes in these magmas is not clearly understood, since their mafic compositions should limit zircon saturation. Here, we investigate the occurrence of zircon (and baddeleyite) in intrusive and extrusive mafic rocks from Central Atlantic Magmatic Province (CAMP) using petrography, trace-element analysis, Ti temperatures, Hf and oxygen isotopes, and high-precision U–Pb geochronology, along with petrological and thermal modeling. We provide new ages for CAMP sills that intruded into Paleozoic sediments in Brazil, indicating that the high and low Ti magmatism in this area occurred synchronously over 264 ± 57 ka. We show that upper crustal assimilation, especially of shales, during the emplacement of the CAMP likely led to zircon saturation. Assimilation of upper crustal sediments is also supported by high δ¹⁸O values and some rare negative εHf values in the zircon crystals. The only extrusive sample analyzed was the North Mountain basalt in Nova Scotia, Canada. This sample contains a large age variation in its zircon crystals (up to 4 Ma), and the older crystals have slightly more negative εHf values suggesting the presence of small (micron scale) xenocrystic cores associated with very late-stage sediment assimilation. However, the CAMP dataset as a whole suggests that the presence of xenocrystic cores is rare. Assuming no xenocrystic cores, and considering the zircon undersaturated nature of LIP mafic melts, the oldest zircon age clusters in a population should record the magma emplacement (or time when assimilation occurred), and the younger ages in a population are more likely to reflect Pb loss, especially given the high U concentrations of LIP zircon. Our identification of heterogeneous isotopic and elemental compositions in LIP zircon indicates that zircon in these magmas saturate in isolated minute melt pockets just before the system cools below its solidus.

Introduction

Zircon and other Zr-bearing phases (e.g., baddeleyite and zirconolite) may be found in rocks from mafic large igneous provinces (LIPs) and are frequently used to provide high-precision age control of these events (Svensen et al. 2009; Schoene et al. 2010a; Svensen et al. 2012; 2017; Blackburn et al. 2013; Burgess et al. 2015; Davies et al. 2017; Greber et al. 2020). Precise U–Pb geochronology from mafic LIPs in the Phanerozoic has led to new insights into how these eruptions influence the Earth’s atmosphere, hydrosphere, and biosphere. For example, there is growing evidence that the intrusion of LIP sills into sedimentary basins that contain abundant organic-rich or evaporitic sediments induce degassing of volatiles, a process that can then lead to mass extinctions (Svensen et al. 2007, 2009, 2010, 2017; Ganino and Arndt 2009; Davies et al. 2017; Burgess et al. 2017; Heimdal et al. 2018, 2019, 2020). This model is the preferred explanation for some of the Phanerozoic mass extinction events, in part due to the observation that the ages of sills intruding the sedimentary basins coincide with mass extinction events. It is, therefore, crucial to understand the systematics of zircon (and baddeleyite) formation in LIPs to correctly interpret the U–Pb age models to accurately determine causal effects of events recorded in the stratigraphic record (e.g., Schoene et al. 2010a).

Despite the presence of zircon (and other Zr phases) in these provinces, very little attention has been given to the geochemistry of these minerals, or to the question of why they even exist in LIP magmas. The low Zr concentrations, high crystallization temperatures, and primitive composition of LIP basalts mean that it is very difficult for them to reach zircon saturation (Boehnke et al. 2013). The presence of zircon (and other Zr phases) in these melts could be the result of crustal contamination changing the magmatic composition or due to crystallization in extremely fractionated, tiny residual melt pockets (e.g., Schaltegger and Davies 2017) or a combination of these processes. Zircon and baddeleyite in these types of rocks have remained poorly studied due to their small size (often < 80 μm along the c-axis), limited abundance, and the aforementioned questions regarding the processes governing their formation. However, advances in geochemical methods, for example measuring Hf isotopes from very low concentration solutions (D’Abzac et al. 2016; Bauer and Horstwood 2018), and measuring in situ oxygen isotopes by SIMS in zircon and baddeleyite (see Linge et al. 2017; Davies et al. 2018), now allow for a detailed investigation of the chemistry of these minerals, helping to constrain their origin and LIP petrogenesis.

Determining the composition of the source and degree of crustal contamination of LIP magmas is crucial for understanding their origin and chemical evolution during emplacement. This is typically assessed through a combination of whole-rock elemental and radiogenic isotopic geochemistry (e.g., Carlson 1984; Peng et al. 1998; Ewart 2004; Reichow et al. 2005; Jourdan et al. 2007; Shellnutt 2014; Callegaro et al. 2017; Marzoli et al. 2018). However, these techniques may be compromised by weathered outcrops, late-to-postmagmatic hydrothermal alteration, and by difficulties in distinguishing between small degrees of crystal contamination during emplacement and source heterogeneities. Zircon is relatively robust to secondary alteration and therefore is a reliable recorder of the isotopic and chemical compositions of the melts from which it crystallized. In other magmatic settings, zircon has provided new insights into the processes occurring in felsic magma chambers (e.g., Wotzlaw et al. 2015; Samperton et al. 2015; Farina et al. 2018; Szymanowski et al. 2020) and also in oceanic lithosphere gabbros (Lissenberg et al. 2009; Rioux et al. 2015a). Therefore, investigating the utility of zircon chemistry combined with high-precision U–Pb geochronology from LIP magmas may provide new clues into the petrogenesis of LIP magmas, especially with respect to processes occurring during the later stages of the magmatic evolution and the contribution of upper crustal contamination to these systems. Also, combining chemical information and high-precision ages will allow for a better understanding of why and under what conditions zircon crystallizes in LIP magmas, which will offer more robust interpretations of the ages and more precise correlations between LIPs and events recorded in the stratigraphic record.

Currently, zircon Hf isotope measurements are not commonly used in studying LIPs especially from mafic rocks (baddeleyite Hf isotopes have been used in a few studies). These data are available for the Emeishan LIP, although mostly from the felsic portion of the province (see review by Shellnutt 2014), Central Atlantic magmatic province (CAMP Davies et al. 2017; Callegaro et al. 2017), the Tarim LIP (Zhang et al. 2016), the Paranà-Etendeka (Hartmann et al. 2019; Rocha et al. 2020), Karoo LIP (Greber et al. 2020), and also some mafic dyke swarms, which may have been associated with flood basalt eruptions; e.g., the dykes in the Qiangtang terrane (Zhai et al. 2013). Other LIP zircon geochemical information, e.g., trace elements, have been measured from some CAMP basalts (Schoene et al. 2010b; Davies et al. 2017), and oxygen isotopes have been investigated from mafic zircon in the Tarim LIP (Zhang et al. 2016) and the Scourie dyke swarm (Davies et al. 2015).

For this study, we measured trace elements, hafnium, and oxygen isotopes from zircons extracted from CAMP magmas that have been dated at high precision across the whole province and also provide new U–Pb ages for CAMP sills in the Brazilian Amazonas basin. We also present some oxygen isotope data from CAMP baddeleyite crystals. We show that the geochemical information from zircon can be combined with geochemical models to understand why zircon is present within these melts and under what conditions it crystallizes. The combination of different datasets also enables us to explain the age variations present in some LIP samples, and to determine the correct emplacement age for these complex samples.

CAMP magmatic province overview

Here, we give a very brief overview of the pertinent CAMP geochemistry and of its relationship to the end Triassic mass extinction. Readers interested in the details should refer to the recent reviews (e.g., Marzoli et al. 2018, 2019).

The CAMP is one of the world’s largest LIPs, with an aerial extent of ~ 10⁷ km². The province was emplaced dominantly as a series of tholeiitic dykes and sills at ~ 201 Ma just before the opening of the Central Atlantic Ocean. The province also contains basalt flows which erupted into rift basins in Morocco, Algeria, Portugal, Brazil, Bolivia, and NE North America and two large layered mafic intrusions in Guinea and Sierra Leone. In terms of timing, the province was emplaced in ~ 1 Ma (Marzoli et al. 1999; Blackburn et al. 2013; Davies et al. 2017; Heimdal et al. 2018) which is similar to other well-known Phanerozoic LIPs (Siberian Traps, Karoo-Farrar, Deccan, Columbia river flood basalts Svensen et al. 2012; Burgess et al. 2015; Burgess and Bowring 2015; Greber et al. 2020; Kasbohm and Schoene 2018; Schoene et al. 2019). The CAMP is, however, distinct from these other LIPs, since it has very few preserved flood basalt lava flows and contains gabbroic intrusions that cover an extremely large area over ~ 8000 km from northern France to southern Bolivia (Jourdan et al. 2003; Bertrand et al. 2014) with a volume of at least 3 million km³. The geochemistry of the CAMP rocks suggests that the province is comprised almost entirely of basalt and basaltic andesite compositions (SiO₂ in the range 48–55 wt.%), with no associated alkaline samples found so far, which is also quite distinct from other LIPs. More evolved CAMP samples are very rare and mostly associated with internal differentiation of large mafic sills, dykes, or flows (Shirley 1987; Puffer et al. 2009; Block et al. 2015) rather than felsic intrusive or extrusive magmas. These rare, evolved magma pockets are the target for this study, since in many cases, they contain zircon.

In general, most CAMP samples are moderately evolved mantle melts with MgO ~ 4–10 wt.% (Marzoli et al. 2018), and with only ~ 1% of the analyzed samples having a primitive composition. with MgO > 10 wt.%. The CAMP also contains some high Ti occurrences (TiO₂ > 2 wt.%), which are restricted to a relatively small area in NE South America, Liberia, and Sierra Leone, but the vast majority of CAMP samples are low Ti basalts (Bertrand et al. 1982; Dupuy et al. 1988; Chalokwu 2001; Nomade et al. 2002; De Min et al. 2003; Deckart et al. 2005; Merle et al. 2011; Callegaro et al. 2017). In terms of radiogenic isotopes (Sr, Nd, Pb, Hf, and Os), CAMP magmas with similar compositions are observed across all four of the continents that host these magmas, but basaltic groups with particular compositions have been recognized in some areas (Marzoli et al. 2018). There is minimal overlap between the isotopic compositions of the low and high Ti occurrences, with the low Ti samples mostly having compositions that suggest the involvement of recycled upper and lower crustal material with data trending toward the EMII mantle end-member (Pegram 1990; Puffer 2001; Dorais and Tubrett 2008; Merle et al. 2011, 2014; Callegaro et al. 2013, 2014, 2017; Whalen et al. 2015; Elkins et al. 2020). The high Ti samples, however, are more similar to depleted mantle but with an additional enriched component often interpreted as old subcontinental lithospheric mantle, possibly present in the form of lamproite veins (Deckart et al. 2005; Merle et al. 2011; Callegaro et al. 2017). Most CAMP samples are also considered to have experienced some assimilation of continental crust (up to ~ 10% based on Os isotope data), the composition of which is likely to be specific to each location (Marzoli et al. 2018). However, disentangling mantle source heterogeneity from crustal contamination using radiogenic isotopes can be complicated, especially if the amount of contamination is small. Hafnium isotopes have also been measured from dated zircon and baddeleyite mineral separates from across the CAMP (Davies et al. 2017; Callegaro et al. 2017 and data presented here) which produced Hf isotope values close to chondritic (εHf ~ − 3 to + 4) for the low Ti CAMP samples and εHf values around + 5 to + 8 for the high Ti samples. These data are in agreement with the available whole-rock Hf isotope measurements (Callegaro et al. 2017; Elkins et al. 2020), and therefore reinforce the CAMP source characteristics described above.

The CAMP magmatism is temporally correlated with the End Triassic mass extinction (Marzoli et al. 2004; Schoene et al. 2010a; Blackburn et al. 2013; Davies et al. 2017; Capriolo et al. 2020). This extinction is one of the largest mass extinction events of the Phanerozoic, characterized in the sedimentary record by a loss of Triassic ammonites and severe destruction of scleractinian corals, early Mesozoic vertebrates, and terrestrial megaflora (Benton 1995; McElwain et al. 1999, 2009; Guex et al. 2004; van de Schootbrugge et al. 2009; Alroy 2010; Mander et al. 2010). There is also evidence for large carbon cycle disturbances associated with the paleontological turnovers; these are preserved by multiple negative carbon isotope excursions (Whiteside et al. 2010; Corso et al. 2014). Also, oceanic photic zone euxinia, increased ocean stratification, and mercury anomalies occurred at the same time (Hesselbo et al. 2002; Guex et al. 2004; Ruhl et al. 2009; Deenen et al. 2010; Whiteside et al. 2010; Bartolini et al. 2012; Kasprak et al. 2015; Lindström et al. 2016, 2019; Percival et al. 2017). One of the currently preferred hypotheses for the cause of the mass extinction is climate change caused by metamorphism of volatile-rich sediments and subsequent huge gas release from the Amazonas basin in Brazil in response to the injection of CAMP sills into the sediments (Davies et al. 2017; Heimdal et al. 2018, 2019, 2020). Recently, it has been suggested that the magmas themselves may have also contributed significantly to the volatile budget (Capriolo et al. 2020).

In this study, we focus on samples from across the entire CAMP province, all of which have been dated, either previously (Davies et al. 2017; Callegaro et al. 2017; Marzoli et al. 2018) or during this study and most of which have been previously studied, by whole-rock geochemical techniques, and for trace element and isotope geochemistry. All of the samples apart from the North Mountain basalt are from mafic intrusions, mostly from sills or dykes, although the Kakoulima sample from Guinea and the Freetown layered complex from Sierra Leone are from large mafic intrusions. All new samples from Brazil are from sills which have intruded the Amazonian basin. We use our new data to determine under what conditions zircon crystallized in these rocks.

Analytical techniques

U–Pb geochronology

U–Pb geochronology was undertaken at the University of Geneva using the techniques described in Davies et al. (2017). Extracted zircon grains were annealed in a muffle furnace at 900 °C for 48 h before undergoing chemical abrasion (Mattinson 2005) for 12 h at 210 °C in concentrated HF (see Widmann et al. 2019). The chemically abraded grains were then ultrasonically cleaned in weak HNO₃ before being loaded into 200 μl Savillex microcapsules with ~ 70 μl of HF and ~ 5 mg of EARTHTIME ²⁰²Pb–²⁰⁵Pb–²³³U–²³⁵U tracer solution (calibration 3.0, Condon et al. 2015; McLean et al. 2015) and placed into a Parr digestion vessel and then put into an oven at 210 °C for 48 h. After digestion, the dissolved zircon solutions were processed through two ion-exchange columns, the first column chemically separated Pb and U fractions for dating, and the second column processed the trace-element fraction from the U–Pb separation to purify the Hf for ICP-MS analysis. The Pb and U fractions were loaded onto outgassed Re filaments with a Si-gel emitter and were measured on a Thermo TRITON thermal ionization mass spectrometer with Pb measured in dynamic peak hopping mode using a MasCom discrete-dynode, secondary electron multiplier, and U measured as UO₂ in static mode using Faraday cups equipped with 10¹²Ω resistors. All common Pb was attributed to laboratory blank and was assigned the long-term isotopic composition of the Geneva procedural laboratory blank. All U–Pb data were processed using the Tripoli and Redux software packages following the algorithms of (McLean et al. 2011).

MC-ICP-MS Hf isotope measurements

Hafnium isotopic analysis followed the previously published methodology (see D’Abzac et al. 2016; Farina et al. 2018). Briefly, Hf was isolated from matrix elements through an ion-exchange chromatography (Augland and David 2015), using the washes from U and Pb chemistry. The Hf cut was subsequently dried and re-dissolved in a 0.3 M HNO₃ + 0.005 M HF solution for analysis on a Thermo Scientific Neptune Plus multiple-collector inductively coupled plasma-mass spectrometer (MC-ICP-MS) at the Department of Earth Sciences of the University of Geneva. The cup configuration was adjusted to measure the isotopes ¹⁷²Yb, ¹⁷³Yb, ¹⁷⁵Lu, ¹⁷⁶Hf, ¹⁷⁷Hf, ¹⁷⁸Hf, ¹⁷⁹Hf, ¹⁸⁰Hf, and ¹⁸¹Ta simultaneously in low-resolution mode. During each session, several Plešovice, Temora, and JMC475 Hf standard solutions were analyzed with the samples. Data reduction to obtain the ¹⁷⁶Hf/¹⁷⁷Hf ratio included on-peak zero baseline correction, correction for mass bias, correction of isobaric interferences of ¹⁷⁶Lu and ¹⁷⁶Yb on ¹⁷⁶Hf and an offset correction by adjusting the ¹⁷⁶Hf/¹⁷⁷Hf ratio of the sample to account for the difference between the measured and preferred value of the JMC475 Hf (i.e., 0.282160) (Nowell et al. 2004). The JMC offset correction ranged between 0 and 1 εHf units depending on the session. Data are presented as initial εHf values using the CHUR composition of (Bouvier et al. 2008), the ¹⁷⁶Lu decay constant (λ¹⁷⁶Lu) of 1.86711 × 10^–11 years⁻¹, and an age of 201.5 Ma. The average ¹⁷⁶Hf/¹⁷⁷Hf ratio of all measured Plešovice and Temora solutions is 0.28247 ± 2 and 0.28268 ± 2, respectively. These ratios translate to a Plešovice εHf_337Ma of − 3.7 ± 0.7 and Temora εHf_417Ma of + 5.5 ± 0.6 (2σ), which are identical to published values (Woodhead and Hergt 2005; Sláma et al. 2008).

SIMS oxygen isotope measurements

In samples with abundant separated zircon grains, some crystals were annealed and mounted for oxygen isotope measurements on the Cameca ims 1280 ion microprobe of the SwissSIMS center located at the Institute of Earth Sciences of the University of Lausanne (Switzerland). A ~ 2 nA current Cs⁺ primary beam, accelerated at 10 kV, was used, resulting in a ~ 10 µm beam size. The electron flood gun, with normal incidence, compensated the charges build-up at the surface of the samples. ¹⁶O and ¹⁸O secondary ions, also accelerated at 10 kV, were analyzed at 2400 mass resolving power and collected on faraday cups in multi-collection mode. Faraday cups are calibrated at the beginning of the session, using the calibration routine. Mass calibration was performed every 12 h. Each analysis takes less than 4 min, including pre-sputtering (30 s) and automated centering of secondary ions. Instrumental mass fractionation (IMF) was determined using Plenglai zircon reference material with δ¹⁸O = XA‰ (Li et al. 2010). This setting allowed a reproducibility of 0.26‰ (2σ) on a Plenglai zircon for the whole session, and the internal error for an analysis was usually ~ 0.25‰ [2 standard error (SE)]. One other zircon reference was measured 18 times during the session; Temora (Black et al. 2004) and yielded a δ¹⁸O value of 8.06 ± 0.31 (2SD, n = 18), within error identical to the published bulk value of 8.2 ± 0.1‰ (Black et al. 2004). For the baddeleyite session, synthetic baddeleyite (monoclinic zirconia) crystals with laser fluorination δ¹⁸O values of 6.39 ± 0.5 ‰ 2σ, (see supplementary information) were used to determine the IMF, and UNIL-Q1 (Paine quartz, Seitz et al. 2017) was used to check the instrument stability. Reproducibility was 0.3‰ (2σ) for UNIL-Q1 and 0.89‰ (2σ) for the synthetic baddeleyite standard. The Phalaborwa baddeleyite was measured as a secondary reference material and gave a δ¹⁸O value of 3.98 ± 1.35 ‰ 2σ, which is consistent with the laser fluorination value of 4.6 ± 0.3‰ 2σ (Davies et al. 2018).

Trace elements

After the SIMS analyses, laser ablation-inductively coupled plasma-mass spectrometry (LA-ICP-MS) trace-element measurements were conducted at the University of Bern using a Lambda Physik GeoLas pro 193 nm ArF Excimer laser system coupled with an Perkin Elmer ELAN-DRCe quadrupole ICP-MS. Measurements followed the procedures of (Pettke et al. 2012), employing an energy density of 4–5 J cm² on the sample surface with a 10 Hz laser repetition rate and beam diameters ranging between 20 and 40 μm, the maximum possible to minimize limits of detection. Measurement points were placed next to ion probe spots for maximum correlation of the data. Daily optimisation of ICP-MS performance (using SRM612 from NIST) ensured maximum signal-to-background intensity ratios, robust plasma conditions with equal sensitivities for U and Th, and a ThO production rate below 0.4%. GSD-1G and SRM 612 were employed as primary standards (with concentration data from GeoRem) for data quantification, Temora 2 zircon was measured along with sample zircons and used as a secondary standard for quality control. Data reduction was conducted using SILLS (Guillong et al. 2008), with the improved limit of detection calculations (Pettke et al. 2012).

Ti concentrations in quartz

Quartz Ti concentrations were analyzed following the techniques of (Ackerson et al. 2018) on a JEOL 8530F hyperprobe at the Geophysical Laboratory, Washington DC, at 15 keV and 200 nA with a 10 μm spot diameter. The Ti concentrations were measured using three PET crystals, with peak counting times of 320 s, which yielded detection limits of ~ 6–8 µg/g. A synthetic quartz crystal (AlQ10) with a Ti concentration of 92 ± 10 µg/g was used as a secondary reference material to check for measurement accuracy. Measurements at the Geophysical Laboratory of the AlQ10 crystal yielded within error identical Ti concentrations of 91 ± 14 µg/g.

Images

Cathodoluminescence (CL) images were acquired using a JEOL JSM7001F Thermal Field Emission SEM with a Schottky electron gun at the University of Geneva using an accelerating voltage of 15 kV. CL images were used to characterize the internal textures of the crystal.

Quartz CL images and also the color thin-section CL image were made using the JEOL 8530F hyperprobe at the Geophysical Laboratory, Washington DC.

Results

U–Pb and Hf isotope analyses were made on bulk individual zircon grains. In situ measurements (oxygen and trace-element analysis) were conducted after dating in samples with abundant, leftover zircon grains. For this, the grains were mounted in epoxy. The different in-situ analyses were made on the same grains and are therefore directly comparable. Errors are presented as 2σ unless stated otherwise (Fig. 1).

Fig. 1

modified from Davies et al. (2017) and McHone, (2000). b A modern simplified geological map of Brazil (a more detailed version is provided in the supplementary information). The CAMP samples are located, and also their ages, εHf, and δ¹⁸O values are reported. The asterixis indicates samples that have been previously studied in Davies et al. (2017). The CAMP sills are located in the Paleozoic sediments, and some of the outcrops are not visible at the scale of the map

Geological map of the Brazilian CAMP samples analyzed here. a A map of the circum Atlantic region during the end Triassic time showing the location of the samples analyzed here; this map is

U–Pb ages

Zircons were separated from four new samples from mafic CAMP sills in the Amazonas basin in Brazil. Three samples were from sills with low Ti concentrations (RP134, RP136, and RP144) and the fourth sample was from a sill with high Ti concentrations (RP116). These samples contained prismatic zircon grains up to ~ 150 μm in length with regular (oscillatory or sector) zoning and relatively homogeneous CL emission. There was no evidence for xenocrystic cores or complex alteration. All samples produced a population of grains that overlapped in age, with a few slightly younger data that were attributed to small degrees of residual Pb loss after chemical abrasion (Fig. 2). All of the analyses are concordant within uncertainty, and therefore, weighted means and associated 2 SE errors of the ²⁰⁶Pb/²³⁸U ages are used here to represent the best estimate for the crystallization age. The three low Ti samples produced weighted mean ages of 201.430 ± 0.034 Ma (MSWD = 1.4; n = 5) from sample RP134, 201.348 ± 0.034 Ma (MSWD = 1.2; n = 3) from sample RP136 and sample RP144 gave an age of 201.372 ± 0.062 Ma (MSWD = 0.87; n = 3). Two zircon crystals with the oldest ages from the high Ti sample produced a weighted mean age of 201.495 ± 0.091 Ma (MSWD = 4.1; n = 2). All of these new ages overlap with previous estimates for the ages of the sills in the Amazonas and Solimões basins in Brazil (Fig. 2b). Our new results, therefore, indicate that the high and low Ti magmas that intruded into the Amazonas and Solimões basin were emplaced synchronously over a period of 264 ± 57 ka.

Fig. 2

New U–Pb ages for CAMP sills intruded into the Amazonas sedimentary basin in Brazil. a Vertical bars represent single zircon ²⁰⁶Pb/²³⁸U U–Pb ages, and darker shaded analyses are grouped together and averaged to calculate the age of emplacement for the sample. Lighter shaded analyses were not included in the average calculation due to suspected Pb loss. b Compilation of published high-precision U–Pb CAMP ages, samples are grouped by their geographic occurrence; note that the high Ti and low Ti CAMP magmas intrude simultaneously. References for colors are 1 this study; 2 (Davies et al. 2017), 3 (Marzoli et al. 2019), 4 (Blackburn et al. 2013), 5 (Heimdal et al. 2018), and 6 (Callegaro et al. 2017); this last age is a baddeleyite age presented as a ²⁰⁷Pb/²⁰⁶Pb weight mean age, the reasons for this are discussed in the text. c Two estimates for the age of the onset of the end Triassic mass extinction as well as a probability density function (PDF—blue) and kernel density (KDE—black) curves for the CAMP ages

Twenty-five new analyses were also conducted on zircon separated from the North Mountain basalt (NMB from here on). These analyses are similar to those published in Schoene et al. (2006, 2010a) Blackburn et al. (2013) and Davies et al. (2017), although the new results contain grains with ages up to 204 Ma, which is ~ 2 Ma older than ages found before (unpublished, similarly old ages have also been found by the Princeton U–Pb lab, Blair Schoene, personal communication). Over half of the ages (64%) agree with the eruption age of the basalt published by Davies et al. (2017) of 201.498 ± 0.028 Ma.

Finally, baddeleyite data from the Freetown Layered Complex (FLC) in Sierra Leone reported in Callegaro et al. (2017) are also plotted in Fig. 2; however, the ²⁰⁷Pb/²⁰⁶Pb weighted mean age of 201.19 ± 0.69 Ma is presented and preferred over the ²⁰⁶Pb/²³⁸U weighted mean age of 198.777 ± 0.047 Ma. This older, less precise age accounts for the effects of unidentified Pb loss in the baddeleyite grains and results in a larger uncertainty. The ²⁰⁷Pb/²⁰⁶Pb age overlaps with the rest of the CAMP ages and is in agreement with the Ar–Ar ages from the same intrusion (see Callegaro et al. 2017) for a discussion on the U–Pb age interpretations).

Hf isotopes

Hf isotopes from the dated zircon crystals were averaged for each sample and are plotted in Figs. 3 and 4. The high Ti sample from Brazil (RP116) has an εHf value of 9.49 ± 0.36, and the low Ti Brazilian samples (RP134, RP136, RP144) have consistent εHf values of 2.20 ± 0.36, 2.26 ± 0.34, and 2.42 ± 0.4 respectively. The Hf isotope compositions were also measured from the Amelal sill in Morocco (sample LV34) dated in Marzoli et al. (2019), giving an εHf value of 0.89 ± 0.44. Also, new zircon Hf analysis from the Tarabuco intrusion in Bolivia and the North Mountain basalt (NMB) in Canada gave εHf values similar to those already published in Davies et al. (2017); however, the updated averages for these samples are now εHf of 3.52 ± 2.04 for the Tarabuco sample, and εHf of 0.15 ± 5.24 for the NMB. The 2σ of the εHf values of the Tarabuco and NMB samples are higher than the errors on the secondary reference materials Plešovice and Temora (ca. ± 0.6), which indicates that these zircon grains contain variations in the Hf isotopic compositions outside of the analytical uncertainty. All new Hf isotope values are consistent with previously published values from CAMP samples (Davies et al. 2017; Callegaro et al. 2017; Elkins et al. 2020). Overall, the data indicate that the high Ti CAMP samples have supra-chondritic Hf (εHf between 6 and 10), whereas the low Ti samples are more juvenile with εHf values between 4 and − 4. There is no clear trend relative to the age with the oldest and youngest samples apparently recording the most negative εHf values. The wide εHf range of the NMB zircon is discussed in detail below.

Fig. 3

Age and εHf isotopic compositions of CAMP zircons from this study, and from Davies et al. (2017). Hf isotopic compositions represent average values for the sample with 2σ uncertainties. Double-headed arrow (bottom left) shows the range of εHf values from low Ti samples from eastern USA (Elkins et al. 2020), the double-headed arrow (top right) for the high Ti samples is the range suggested for these samples in this study. Asterisks (*) represent the new U–Pb and Hf samples from this study, capital delta symbols (Δ) indicate samples from (Davies et al. 2017) with updated Hf isotopic compositions

Fig. 4

Oxygen isotope data for the CAMP zircons and also for Freetown Layered Complex FLC baddeleyite. a Box and whisker plots as well as the individual δ¹⁸O measurements; arrows indicate two very anomalous analyses from the Tarabuco and Fouta Djalon samples. b and c Weighted mean δ¹⁸O values for each sample plotted relative to the sample age and εHf values. The FLC baddeleyite has an average δ¹⁸O of 2.74 ± 1.83 ‰ at 2σ, the zircon value calculated for equilibrium with the baddeleyite is indicated by the transparent symbol and arrow, and this value is in agreement with the other CAMP zircon values. See text for further explanation

Oxygen isotopes

The oxygen isotopic composition of zircon from the CAMP samples is relatively homogeneous with most analysis falling in between δ¹⁸O values of 5 and 6‰ (Fig. 4). For most samples, the majority of grains have very consistent δ¹⁸O values; however, most samples also contain 1–3 analyses with higher values, and samples RP134 and RP113 both have one analysis that falls below the main cluster, and the averages reported below and in Fig. 4 are from the main group of data without the outliers. Broadly, the data fall into two groups, the first group contains Kakoulima, Tarabuco, Messejana, Amelal sill, Fouta Dhalon, and the Brazilian high Ti samples, and has a δ¹⁸O value of ~ 5.5–6‰. Whereas, the second group has a δ¹⁸O of ~ 5‰ and contains the NMB and the Brazilian low Ti samples (Fig. 4). When the oxygen isotope data are compared with the sample ages, there is a slight reduction in δ¹⁸O over time from δ¹⁸O of 6.03 ± 0.17 ‰ for the Kakoulima intrusion, to δ¹⁸O of 5.04 ± 0.09 ‰ for RP134 (Fig. 4b). Also, there is no trend between δ¹⁸O and εHf values (Fig. 4c). Similar to the εHf isotope data, the NMB sample has a large range in δ¹⁸O compared with the other samples (Fig. 4a). There are also no correlations between the δ¹⁸O values and zircon trace-element enrichments, which are indicative of secondary alteration (Geisler et al. 2007; Davies et al. 2015; Bell et al. 2019 and supplementary data). Therefore, the δ¹⁸O data are all considered to record magmatic processes, apart from the two outlier measurements, one from Tarabuco which has a δ¹⁸O value of 12.00 ± 0.31 ‰, and one from Fouta Djalon which has a δ¹⁸O value of 9.76 ± 0.31 ‰. These data points are anomalous and not considered further.

Baddeleyite crystals from the FLC were also measured for their oxygen isotope compositions (supplementary data). These measurements have larger uncertainties due to the inhomogeneity of the reference materials used, which is then propagated onto the unknowns (see Davies et al. 2018), and therefore, the individual analyses are not plotted with the zircon data in Fig. 4a. However, the baddeleyite grains have a homogeneous δ¹⁸O value of 2.74 ± 1.83 ‰. Using the equilibrium isotopic fractionation factor between zircon and baddeleyite (Δ_{(zircon-baddeleyite)} = 4.3*10⁶/T²) suggested by Davies et al. (2018), and a temperature of 850 °C (based on Ti in zircon thermometry, see below), the calculated zircon δ¹⁸O value that would be in equilibrium with the FLC baddeleyite is 6.14 ± 1.83 ‰, which is within uncertainty of the other CAMP samples (Fig. 4b, c). Oxygen isotopes have rarely been measured in previous studies on the CAMP magmas, and therefore, there is not much data with which to compare the zircon analysis measured here. One study on the FLC found whole-rock oxygen isotope values ranging between δ¹⁸O of 5.4 and 6.3 ‰ (Chalokwu et al. 1999), which is in agreement with our δ¹⁸O values. Oxygen isotopes were also measured in the NMB from pyroxene mineral separates in silica-rich zones likely related to hydrothermal alteration (Kontak and Dostal 2010). The pyroxenes have extremely high δ¹⁸O values between ~ 10 and 20 ‰, and were interpreted as being altered by, or crystallized from 100 to 200 °C hydrothermal fluids.

Trace elements in zircon and Ti in quartz

The trace-element concentrations of separated CAMP zircon crystals are shown in Fig. 5. In general, the trace-element concentrations are high and variable for all samples. The high U concentrations signify that many of the grains are metamict and are, therefore, more susceptible to fluid overprinting (Geisler et al. 2007). The analyses with elevated Fe, light rare-earth elements (LREE), and Na were considered to have been affected by fluid alteration or contained unidentified fluid inclusions and we omitted all data from the further discussion following the cut-off values for altered zircon compositions of Bell et al. (2016) (supplementary data). None of the previously published trace-element data from CAMP zircon have LREE, Fe, and Na concentrations, since they were measured via electron probe (Davies et al. 2017) or from trace-element washes after U–Pb chemistry (Schoene et al. 2010b). These previously published analyses also have the highest concentrations of trace elements (e.g., > 6000 µg/g U), and therefore, it is possible that they also record secondary trace-element enrichments through fluid alteration, especially due to the metamict nature of some of the high U grains.

Fig. 5

Trace-element concentrations of CAMP zircon from this study (large symbols) and previously published zircon compositions (small symbols) from (Schoene et al. 2010b; Davies et al. 2017). a Zircon compositions plotted relative to the continental zircon and oceanic zircon fields of (Grimes et al. 2007); note the overlap of CAMP zircon with both fields and also that some analysis lies outside the proposed fields. b CAMP zircon plotted relative to the continental arc zircon, S- and I-type granite zircon, TTG zircon, Icelandic zircon, and zircon from the Bishop Tuff. Compilation of literature zircon data is from (Chamberlain et al. 2014; Carley et al. 2014; Reimink et al. 2020) and references therein. c Box and Whisker plots of Ti temperatures in zircon and quartz for the CAMP samples assuming an Si activity of 1 and a Ti activity of 0.6

The unaltered zircon compositions partially overlap with the continental and oceanic zircon fields in U/Yb vs Y representation from (Grimes et al. 2007); however, they also plot outside these fields (Fig. 5a). The grains have very high U concentrations and high Th/U as noted by Schaltegger and Davies (2017). They fall in areas between the fields produced by Icelandic rhyolite zircon (Carley et al. 2014) and the Bishop Tuff zircon dataset from (Chamberlain et al. 2014), both of which represent relatively extreme cases of zircon trace-element enrichment (see the continental arc dataset, S- and I-type granite, and TTG datasets in Fig. 5b). To a first order, the trace-element compositions of LIP type zircon are quite distinct in the amount of trace-element enrichments that they contain. Chondrite normalized REE plots are provided (Supplementary Figure 5) and indicate that all of the zircon have an REE similar pattern apart from the NMB which is both more enriched in LREEs and variably depleted in HREE’s.

The calculated Ti temperatures from these zircon crystals, assuming a Ti activity of 0.6 (after Hayden and Watson 2007), and the lack of associated rutile) and a Si activity of 1 (due to the presence of quartz associated with zircon), are given in Fig. 5c. Notably, almost all of the zircon appears to crystallize at relatively low temperatures between 750 and 850 °C, only the NMB sample has some zircons that crystallized at higher temperatures. The high-temperature NMB crystals appear to correlate with lower δ¹⁸O values (Supplementary data). The low temperatures for most of the CAMP samples are consistent with those found in zircon from mid-ocean ridge gabbros (Rioux et al. 2015b) and gabbro intrusions (Fu et al. 2008). They are also consistent but slightly higher than the Ti temperatures from quartz crystals (calculated assuming the same Ti activity of 0.6) found in the same evolved granophyric pockets as the zircon crystals (Fig. 6). The quartz temperatures range from ~ 600–800 °C for the NMB (Fig. 5c) and are therefore a little lower than the zircon temperatures, suggesting that the bulk of the quartz crystallizes coevally with and just after the zircon. It should also be noted that the NMB zircon temperatures are all above the quartz temperatures, suggesting that a_SiO2 was not 1 during zircon crystallization in this sample (Ferry and Watson 2007). Reducing the activity of silicon in the zircon temperature equation results in lower temperatures, more in line with the quartz, but not overlapping.

Fig. 6

Thin section and representative cathodoluminescence (CL) images of zircon and evolved melt pockets from inside CAMP basalts. Images a and b were taken using a cold-cathode CL detector, whereas the other images (c and f) are made using a panchromatic black and white CL and show zircons extracted from samples, or backscattered electron images from NMB or Foum Zguid thin sections (images d and e)

Discussion

The elemental and isotopic chemistry of zircon from CAMP magmas varies within and between samples. The significance of these variations for understanding the origin and the conditions for zircon crystallization in these magmas are explored below.

When does zircon form in LIP magmas and under what conditions?

Zircon (ZrSiO₄), baddeleyite (ZrO₂), and zirconolite (CaZrTi₂O₇) are found in interstitial melt pockets within CAMP samples (Fig. 6). These melt pockets are typically at the scale of ~ 100 μm, although in the NMB, they can reach the 10 cm scale (Kontak and Dostal 2010). They contain K-feldspar, quartz, apatite, plagioclase, ilmenite, and minor sulfides (that are likely secondary). This mineral assemblage suggests highly fractionated SiO₂-rich melts and is similar to the mineralogy of melt pockets previously identified in CAMP, and other LIPs (Kontak et al. 2002; Kontak and Dostal 2010; Ver Hoeve et al. 2018; Heimdal et al. 2019). The highly variable and enriched trace-element concentrations of the zircons are also consistent with crystallization in highly fractionated, silica-rich melts (Fig. 5). It is clear that none of the Zr-bearing minerals crystallized with the main fractionating assemblage of pyroxene, olivine, and plagioclase, which is not surprising, since zircon saturation conditions in basaltic (tholeiitic) magmas at temperatures > 1000 °C require > 7000 µg/g Zr (Boehnke et al. 2013), which is about 50 times more than the typical content of CAMP basalts. Also, the Ti temperatures of the zircon are mostly around 800 °C suggesting much lower crystallization temperatures (Fig. 4). For a more precise understanding of the zircon saturation conditions, the composition of potential primary magmas can be used (see, e.g., Callegaro et al. 2013; Merle et al. 2014; Whalen et al. 2015; Marzoli et al. 2019). If we use primary CAMP sample AN133 from Marzoli et al. (2019) as an estimate for a primary low Ti magma, with a Zr concentration of 119 µg/g, the zircon saturation conditions can be calculated based on the degree of polymerization of the melt, which can be determined based on the cation ratio (Na + K + 2Ca)/(Al × Si), referred to as M (Watson and Harrison 1983; Boehnke et al. 2013). Higher values of M indicate low degrees of polymerization and enhanced zircon solubility. The AN133 composition gives an M value of 2.77 resulting in an unrealistically low temperature of zircon saturation of 578 °C, which is far below the wet granite solidus ~ 700 °C (e.g., Ebadi and Johannes 1991), clearly indicating that zircon will not crystallize without significant modification through crystallization or assimilation. Conversely, this also means that antecrystic zircon will be readily dissolved in these melts, and therefore, old zircon ages should not normally be attributed to antecrysts.

The zircon saturation conditions in CAMP magmas can be refined through forward modeling using rhyolite-MELTS (version 1.2.x, Gualda et al. 2012; Ghiorso and Gualda 2015), and monitoring the changing composition of the magma as it undergoes fractional crystallization from ~ 1200 °C to ~ 700 °C at 1.5 kbar pressure, which is ~ 4 km depth. The Zr concentration of the melt can be calculated at each step using the mineral phase proportions along with mineral/melt Zr partition coefficients for each phase present (the partition coefficients used are reported in the Supplementary data). The Zr concentration in the melt is then compared to the zircon saturation conditions at every step, calculated using the equations of Boehnke et al. (2013). Where the Zr concentration and the zircon saturation conditions coincide is the likely point of zircon saturation (Fig. 7a). In Fig. 7, the model used sample AN133 as a starting composition (low Ti basalt; Marzoli et al. 2019) and different initial H₂O concentrations to show the effect of initial water on the zircon saturation conditions. Note that the composition of the liquid from which the zircon crystallizes is granitic with ~ 73% SiO₂ for each of the initial H₂O concentrations (see TAS diagram showing the composition of the melt at zircon saturation in the supplementary information). When using the high Ti CAMP sample M13 from Merle et al. (2014) as starting material, zircon saturation was not achieved under any conditions, implying that zircon does not saturate in such a melt, unless they are significantly contaminated (see next paragraph). This is also the case for CS49, which is a low Ti primary CAMP basalt from South-Eastern USA (Callegaro et al. 2013).

Fig. 7

MELTS modeling results of zircon saturation conditions in CAMP magmas. a Plot of modeled Zr concentration in the melt during fractional crystallization vs temperature, combined with the change in zircon saturation conditions at each fractional crystallization step. The point where the Zr concentration curve crosses the zircon saturation conditions is the point at which zircon will saturate in the melt; this point is marked with a star. Models are shown for the primary CAMP sample AN133 (low Ti) with three different initial water concentrations, models are run in 10 °C increments. The colored bar indicates the fraction of melt remaining at each fractional crystallization step, and is only shown on the Zr concentration points for simplicity. Other CAMP primary melts—CS49 (low Ti) and M13 (high Ti)—do not achieve zircon saturation (i.e., the curves do not overlap) with any initial H₂O concentrations. b The variation in zircon saturation point is shown for the AN133 with 0.5% initial H₂O model (the same as in a) after different assimilants are added to the melt, and the assimilants are shale LCC (lower continental crust), granite, and SST (sandstone), the compositions are given in the supplementary data. The zircon saturation point with no assimilation is shown as a white star—at the same location as in a. The assimilants are added after 40% fractional crystallization (at ~ 1100 °C), and the proportion of assimilation is either 10% (closed symbols) or 5% (open symbols). Assimilation after 55 and 60% fractional crystallization is also shown for granite. The insert in the bottom left corner shows the zircon saturation point for SC49 and M13 after 10 and 15% shale assimilation, respectively, zircon saturation is not reached in these melts with the other assimilants, or lower proportions of assimilation. The colored bar indicates the melt fraction remaining at the zircon saturation point. Note that the shale assimilation for AN133 (main part of the figure) causes zircon saturation to be reached at the earliest point during fractional crystallization

It is clear from Fig. 7a that the initial H₂O concentration of the melt is very important in determining when zircon saturates, this is because H₂O suppresses silicate saturation (also see TAS diagram in supplementary information). With an initial H₂O concentration of 0.5 wt. % (similar to the initial water concentrations estimated in Callegaro et al. 2013), the magma becomes zircon saturated at ~ 790 °C, after 85% crystallization. Increasing the initial H₂O concentration of the melt decreases the needed amount of fractional crystallization and also the temperature for zircon saturation. Based on our modeling, zircon saturation is reached after 82% fractional crystallization at ~ 750 °C for an initial H₂O concentration of 0.8 wt. % and at ~ 740 °C after 80% fractional crystallization with an initial H₂O concentration of 1.5 wt. %. Comparing these modeled zircon saturation conditions with the newly determined zircon crystallization (Ti) temperatures of ~ 800 ± 50 °C based on the Ti in zircon thermometer (Fig. 5c) suggests that the 0.5% initial H₂O model provides the most consistent saturation condition estimates. It should be noted that the zircon crystallization temperatures were calculated assuming a silica activity of a_SiO2=1 and a Ti activity of a_TiO2 = 0.6 both of which have an effect on the crystallization temperatures calculated. Reducing the silica activity or increasing the activity of Ti (see discussion on the NMB data above) results in a reduction in calculated zircon temperature, which would be more in agreement with higher initial H₂O. Zircon crystallization under equilibrium conditions is also assumed, which may not be the case in these fast cooling, fractionated melt pockets, if kinetic effects (e.g., Albarède and Bottinga 1972) play a role in the trace-element compositions of the zircon. This has recently been theoretically demonstrated for Zr isotopes in zircon (Chen et al. 2020; Méheut et al. 2021), and hence, the Ti concentrations in our zircon may not only be dependent on temperature. High initial H₂O concentrations have not been previously proposed for the CAMP source rocks, which is in agreement with the typically anhydrous mineralogy of CAMP tholeiites, which contain olivine, pyroxene plagioclase, and oxides without significant amphibole (although amphibole and biotite may be present in some of the evolved melt pockets, Heimdal et al. 2019). In a recent study, however, Capriolo et al. (2020) identified deeply sourced CO₂-rich bubbles trapped in melt inclusions in CAMP minerals, and melt inclusions associated with these gas bubbles were also found to contain ~ 1 wt % H₂O, suggesting that the source for CAMP melts likely had some initial volatile component.

Crustal contamination

Crustal contamination is also thought to play a role in the generation of the CAMP magmas, and this is also likely to affect the zircon saturation conditions, since it will affect the M parameter, and the Zr concentration of the melt. The changes in saturation conditions may even result in xenocrystic zircon being preserved. Examples of CAMP rocks with xenocrystic zircon are the Orange Mountain basalt (one of the CAMP basalts erupted in NE USA; Blackburn et al. 2013; also see NMB discussion below) along with a dyke from Morocco (sample AN733; Davies, unpublished data).

The amount of crustal contamination in the CAMP magmas is thought to be ≤ 10% based on Os isotope modeling (see Merle et al. 2011, 2014; Callegaro et al. 2014, 2017; Marzoli et al. 2018). To assess the impact of this on the zircon saturation conditions, we modeled the effects of assimilation of various crustal rocks, i.e., shale, lower continental crust (LCC), granite, and sandstone (SST; Fig. 7b) (the compositions of the assimilants are given in the supplementary data). The assimilation was simulated using the magma chamber simulator (Bohrson et al. 2014, 2020; Heinonen et al. 2019, 2020), which is a thermodynamic model that calculates the evolution of a composite magmatic system and its fractionally crystallizing minerals. It also contains a number of sub-systems that control the thermodynamics and chemistry of magma recharge, assimilation, stoping, and accumulation of cumulates. The magma chamber simulator uses rhyolite-MELTS to compute phase equilibria, and all of the different sub-systems thermodynamically interact, e.g., the intrusion of magma heats up the wall rock, possibly causing melting which mixes with the melt (see Bohrson et al. 2014, 2020 for a detailed explanation of the model). Here, the magma chamber simulator was used to constrain the effects of assimilating via stoping, different possible contaminants into CAMP primitive magmas, and determine the associated effects on zircon saturation. The bulk composition of the assimilated material was added to the CAMP magma as a liquid after a certain amount of fractionation (the stoping method, see Bohrson et al. 2020), using the magma recharge function of the magma chamber simulator. All assimilation models were run using AN133 as the starting composition to represent a primary CAMP melt, with 0.5% initial H₂O. Also, most CAMP samples are not primary melts, and represent liquids after various degrees of fractional crystallization (10–50%) at low pressure (Marzoli et al. 2018). Therefore, the modeled assimilation proportions were 10 or 5% and assimilation occurred after 40% fractional crystallization. For reference, the effects of assimilating after 55 and 60% fractional crystallization are also shown for granite (Fig. 7b). To model the effects on zircon saturation, the Zr concentration of the melt needs to be known and this was calculated using partition coefficients (given in the supplementary information) between the melt and minerals crystallizing at any time step. It is clear from the modeling results that assimilation has a smaller impact on the zircon saturation conditions than increasing the initial H₂O (note in 7a the position of the zoom shown in 7b), the exception is assimilating shale, which, in our example, had 5% H₂O. Only shale assimilation had a significant effect on the zircon saturation conditions, and resulted in zircon saturation being reached after 79% or 75% fractional crystallization for the 5 and 10% assimilation models respectively. The M value for shale assimilation is also the lowest for all of the models (2 and 1.9 for the 5% and 10% assimilation models, respectively, whereas all other assimilation models have M values ~ 2.2), which explains the drastic change in saturation conditions. Other primary CAMP compositions, for example, CS49 (Callegaro et al. 2013), only reach zircon saturation after 10% assimilation of shale, and even then, it requires 91% fractional crystallization with an M value of 2.7 at saturation, indicating that the melt only just reaches zircon saturation. High Ti CAMP samples (M13, from Merle et al. 2011), which have high Zr, but are geochemically more depleted than the low Ti CAMP samples, reach zircon saturation only after > 15% assimilation of shale. Zircon saturation in these samples is reached after 82% fractional crystallization, with an M value of 2.32 (see insert in Fig. 7b).

Overall, most assimilants do not affect the zircon saturation conditions too much if the melt will saturate zircon on its own without assimilation. The saturation conditions with 0.5% initial H₂O and some degree of assimilation are consistent with the thin-section evidence of zircon in highly fractionated melt pockets (Fig. 6), high trace-element contents in the CAMP zircon, and also the Ti temperatures, suggesting that the modeling may reflect real petrological processes (Fig. 5). Many primary CAMP melt compositions will not saturate zircon without significant contamination, and shale seems to be the contaminant that most easily creates the conditions for zircon saturation. The exact nature of the contaminant for each CAMP sample is likely to be different, since the samples cover the entire province (i.e., over ~ 8000 km distance from the NMB in Canada to the Tarabuco sill in Bolivia). However, it is clear that assimilation (especially of shale) played a role in enabling the magmas to reach zircon saturation. Also, it is crucial (for geochronology) for magmas to reach zircon saturation early, so that zircon has longer time to form, and larger crystals are more likely; only the assimilation of shale causes zircon saturation to occur significantly earlier, therefore, these are the cases where dateable zircon is more likely to be present. Assimilation after > 40% fractional crystallization also causes zircon saturation to occur earlier (see darker blue stars in Fig. 7b) promoting zircon growth.

What can the isotopic compositions of zircon tell us about the contamination processes present in CAMP melts?

To use the O isotopic data to help understand petrologic processes during CAMP magmatism, the zircon δ¹⁸O values need to be corrected for the effects of temperature-dependent fractionation, and converted to basalt values assuming equilibrium isotope fractionation between basalt and zircon. This is a large assumption, since the magma that the zircon crystallized from was not basaltic, and was likely modified by upper crustal contamination (Figs. 5, 7). However, given these caveats, δ¹⁸O values in zircon that are different from the mantle value after the effects of equilibrium isotopic fractionation between basalt and zircon have been accounted for, may be used, in combination with the εHf values and other available data, to understand the magmatic processes that occurred during emplacement. Most of the CAMP samples contain a homogenous population of zircon δ¹⁸O values, with a few outliers, mostly at higher values.

To calculate the equilibrium isotope fractionation between basalt and zircon, the fractionation factor needs to be known. Experiments to determine the equilibrium zircon-WR (whole rock) oxygen isotope factors have been attempted (e.g., Trail et al. 2009), along with empirical observations of a general relationship between WR δ¹⁸O value and that of zircon (e.g., Lackey et al. 2008). However, there currently is no consensus on the correct fractionation factors to use. We applied the 1000 × ln(α_Zrc-Fo) = 0.4 × 10⁶/T² from the first principals’ calculations of (Kieffer 1982) combined with 1000 × ln(α_basalt-Fo) = 1.4 × 10⁶/T² from Eiler (2001), to obtain 1000 × ln(α_basalt-Zrc) = 1.0 × 10⁶/T², with α being the fractionation factor and T being temperature in Kelvin. Using our equilibrium isotope fractionation factor, combined with the average Ti temperature determined for each suite of zircon crystals, we can determine a δ¹⁸O value for the basalt in equilibrium with the zircon. These basaltic values are shown in Fig. 8a along with a range of possible basalt oxygen isotopic compositions that could potentially be produced through fractional crystallization of a tholeiitic basalt at low pressure (after Bucholz et al. 2017). It is clear that some of the calculated basalt δ¹⁸O values are higher than would be expected through fractional crystallization of a primitive basalt alone; therefore, they are likely recording some other processes. Also, as expected, the higher δ¹⁸O values from some of the grains result in extremely high basaltic values of 7.5–8 ‰.

Fig. 8

Relating the oxygen isotope composition of zircon to calculated δ¹⁸O values for basalt assuming equilibrium fractionation and also relating the basalt values to the zircon εHf values. a Zircon δ¹⁸O oxygen isotope data converted to basaltic values assuming equilibrium isotopic fractionation between the zircon and basalt at the Ti temperature recorded by the zircon (see text for further explanation). The zircon δ¹⁸O are separated into the average value for the sample shown by filled circles and also the high values shown by open circles of the same color. The range for the δ¹⁸O of mantle-derived tholeiitic basalt that has undergone high degrees of fractional crystallization is shown by the gray bar (see Bucholz et al. 2017). b Calculated basalt δ¹⁸O values for each sample are compared to the zircon εHf values. Mantle end members are shown along with an example lamproitic composition from (Griffin et al. 2000) and a mixing line between the DMM end-member and lamproite is shown with 2% increments (horizontal mixing line); lamproite values range in εHf between − 10 and − 35 in Griffin et al. (2000); here, we use -20 as an intermediate value. A shaded area representing the lower crust is from (Kempton and Harmon 1992) and (Vervoort et al. 2000). Vertical mixing lines show trajectories for mixing with upper crustal silicic sediments with various compositions (e.g., Pettke et al. 2002; Bindeman 2008; Rickli et al. 2013; Bindeman et al. 2016; Greber et al. 2017). Also shown is the mixing trajectory between CAMP basalt and an example of an upper crustal granite (Valley et al. 1994; Yang et al. 2006)

There are numerous models for the source of the CAMP basalts, and most of these involve mixtures of multiple mantle components and mixing with a small amount of either subducted material or lower continental crust in the case of the low Ti samples, or enriched metasomatic veins in the subcontinental lithospheric mantle in the case of the high Ti samples (Pegram 1990; Puffer 2001; Dorais and Tubrett 2008; Merle et al. 2011, 2014; Callegaro et al. 2013, 2014; Whalen et al. 2015; Marzoli et al. 2018). The different mantle end members are unlikely to explain the δ¹⁸O values seen here (Fig. 8a); also, melts would need to contain up to 100% lower continental crustal material to explain the δ¹⁸O enrichments, since the proportion of oxygen between mafic mantle-derived melts and lower continental crust is approximately equal. The most likely explanation for the high δ¹⁸O values recorded in the CAMP zircon is contamination from high δ¹⁸O upper crustal sediments during emplacement, which is also in agreement with the zircon saturation modeling (Fig. 7) and NMB εHf data (Fig. 9).

Fig. 9

Distribution of ²⁰⁶Pb/²³⁸U dates and εHf values for single zircon grains from the North Mountain basalt (NMB). a Individual zircon dates with their 2σ uncertainties. Analyses are colored by their radiogenic to common Pb ratio (Pb*/Pb_c), and the size of the uncertainties broadly correlates with this ratio, with lower values having higher uncertainties due to less radiogenic Pb. b Comparison of zircon ²⁰⁶Pb/²³⁸U dates with εHf values of the same grains; graph also shows mixing relationships with upper crustal sediments (dashed line) that explain the trends to old ages, and also the potential influence of dissolution precipitation processes which may explain the positive εHf values

A simple mixing relationship between a mantle-derived basalt with an εHf value of + 3.5 (for the low Ti CAMP samples; Elkins et al. 2020), a Hf concentration of 2.88 µg/g (an average of Moroccan CAMP samples, Marzoli et al. 2019), and a δ¹⁸O value of 5.8 ‰ (to reflect a small amount of enriched component in the source), with upper crustal sediments, for example shales, that have εHf of − 10 with 6 µg/g Hf (e.g., Bayon et al. 2006), and a δ¹⁸O value of 16 (Bindeman et al. 2016), can easily explain the variation seen in zircon data. This simple mixing model plotted with 2% mixing increments (steep mixing lines in Fig. 8b) shows that only 10–20% of a sedimentary shale component is required to explain the most positive δ¹⁸O value of ~ 8‰ found in the Amelal sill. However, it is extremely unlikely that a single crustal rock is the contaminant for all of the CAMP magmas. Therefore, the mixing calculations shown in Fig. 8 involve various hypothetical sediment compositions and have a range in Hf concentration and εHf values which covers the entire range of values seen in the zircon data. These hypothetical sediments cover compositional ranges of sediments derived from mafic to felsic rocks, although all have elevated δ¹⁸O values (see compilations in Vervoort et al. 2000; Pettke et al. 2002; Bindeman 2008; Rickli et al. 2013; Bindeman et al. 2014, 2016; Greber et al. 2017). Also shown is a mixing curve with a hypothetical granitic sample with 7 µg/g Hf, εHf of -10, and δ¹⁸O of 8‰, with mixing increments of 2%, and up to a maximum of 50%. It is clear that contaminants with low δ¹⁸O values cannot explain the elevated δ¹⁸O values of zircon. It should be noted that the CAMP Hf and O end members chosen here do not explain all of the zircon data; for example, the NMB data require a CAMP source with lower εHf. The overall conclusion is that some amount of crustal contamination is required to explain the variable δ¹⁸O values in zircon.

This conclusion is corroborated by other independent proxies. For example, samples that show high δ¹⁸O values, notably the Kakoulima and Fouta Djalon intrusions and the Amelal sill, have whole rock and Nd, Sr, and Pb isotopic compositions that suggest significant crustal assimilation (up to 30%) (Deckart et al. 2005). The Amelal sill intrudes lacustrine silts (Marzoli et al. 2019) and the NMB has direct zircon evidence of sediment contamination (Fig. 9). The zircon δ¹⁸O data for the other samples do not require sediment assimilation to explain their isotopic composition, since they can be explained through high degrees of fractional crystallization of basaltic melt. However, it is clear that zircon saturation can be promoted by minor amounts of crustal contamination and this may not change the O-Hf isotopic compositions significantly. The high Ti CAMP samples do not require assimilation to explain their oxygen and hafnium isotope data. The more positive εHf values found for the high Ti sample can be explained using the same mechanisms as those described in Merle et al. (2014) and Callegaro et al. (2017), where a depleted mantle source is mixed with a small amount of highly enriched mantle melt, for example a lamproite (see horizontal mixing line on Fig. 8). However, as discussed above, contamination is likely required to explain the presence of zircon in these samples.

In conclusion, we suggest that late-stage addition of sediments into the melts can create the conditions for zircon saturation to occur and, in some cases, also explains their isotopic compositions. Upper crustal contamination is also consistent with bulk rock radiogenic isotope data. Therefore, zircon isotopic and elemental compositions from LIP magmas can help to understand the conditions in which these crystals grew, but such data are not well suited for determining the source of the melts, or their tectonic setting.

Why does zircon in the North Mountain basalt have such variable ages, trace element, and isotopic compositions?

Zircons from the NMB display more variable ages, Hf and O isotopic compositions, and higher Ti temperatures than zircon from all the other CAMP samples (Figs. 4, 9). The NMB has been dated by high-precision U–Pb techniques in multiple studies (Schoene et al. 2010a; Blackburn et al. 2013; Davies et al. 2017, and the current work), and the recommended eruption age of 201.498 ± 0.028 Ma is an average of these studies. New data obtained here (Fig. 9) including some analyses from Davies et al. (2017) suggest zircon crystallization over ~ 2.5 Ma, which seems highly unlikely given the zircon undersaturated nature of the CAMP basalts. The simplest explanation for the old ages is the presence of xenocrystic cores. Truly inherited grains are unlikely since the ages of the zircon from sediments the NMB intrudes and erupts onto are Carboniferous and older (> 360 Ma, Marzoli et al. 2017). However, it is possible for xenocrystic zircon to survive if they were incorporated just before emplacement and the magma was not able to fully resorb them before cooling below the solidus (e.g., Tang et al. 2014).

To identify potential xenocrystic cores in the NMB zircon, cathodoluminescence (CL) imaging was conducted (Fig. 6; Supplementary data); however, no discernable xenocrystic cores within the grains were detected despite imaging of > 100 grains. Many of the zircons contain melt channels, whereas most have relatively monotonous CL response without significant oscillatory or sector zoning, a typical feature of zircon from mafic magmas (e.g., Rioux et al. 2012; Grimes et al. 2013). Hafnium isotopic compositions can also help to identify xenocrystic zircons. For the NMB, most of the zircon with ages overlapping with the time of the eruption have εHf values of ca. − 0.5 (Fig. 9), similar to other low Ti CAMP magmas (Figs. 3, 4). The older grains have slightly more negative εHf values (around − 1.5 εHf), although this shift is minor and within error of the younger grains. These more negative values at older ages could be explained by the incorporation of tiny (referred to as nano) xenocrystic cores into the grains. The oldest zircon (204.2 ± 0.1 Ma) can be explained by 3% mixing between a xenocrystic core of 360 Ma with εHf − 16 (based on zircon data from similar sediments in the area (Marzoli et al. 2017; Waldron et al. 2019) and NMB magmatic zircon using the age of the eruption and an εHf − 0.5. The mixing trend is very steep with very small changes in the εHf value. Using this mixing proportion, further constraints can be placed on the hypothetical xenocrystic cores based on the 9 pg of radiogenic Pb present within the grain. Only 0.3 pg of (old) Pb is required from the xenocrystic core, suggesting that it was likely extremely small (potentially ≤  ~ 1 μm diameter), assuming that the Hf concentrations are approximately equal between the nano-xenocryst and magmatic zircon. The other grains with ages older than the eruption with similar εHf values require even smaller xenocrystic cores to explain their age and εHf values. There are much older zircon grains in the sediments under the NMB; however, mixing with these older grains cannot explain both the shift toward more negative εHf and the age shift seen in the NMB zircon. Given the fact that the Ti concentrations in the NMB zircon indicate an unusually high zircon crystallization temperature of ~ 900 °C, we suggest a scenario where the sandstones and shales underlying the NMB were entrained into the basaltic melt, quickly dissolved (e.g., Hansen and Nielsen 1999; Heimdal et al. 2019), and locally increased its zircon saturation temperature (see Fig. 7b, where assimilation of sandstone increases the zircon saturation temperature). After dissolution of almost all of the xenocrystic zircon, the last nano-scale zircon cores became the nucleus for magmatic zircon growth. These tiny cores are also unlikely to be present in the CL images due to their small size, resulting in a low chance of exposure during polishing of a grain mount.

There are also two grains with ages close to the eruption age that have εHf values of ~  + 4 and + 9 that cannot be explained by the presence of nano-xenocrystic cores (Fig. 9). It is impossible to add enough juvenile Hf from a xenocrystic core to change the εHf value of a CAMP zircon from -0.5 to an εHf value of + 4 or + 9 without significantly affecting the age. Hafnium isotope variability not associated with xenocrystic zircon is commonly observed in felsic I- and S-type magmas (e.g., Kemp et al. 2007; Villaros et al. 2012; Tang et al. 2014), and multiple potential explanations for this phenomenon have been proposed, including magma mixing (Appleby et al. 2010), antecrystic zircon (Annen et al. 2015), or disequilibrium melting (Tang et al. 2014). However, in mafic systems such as the NMB, the zircon saturation conditions should not allow the preservation of antecrystic or xenocrystic zircon, since they should readily dissolve in the magma (Fig. 7). Isotopic heterogeneities resulting from mixing of magmas would likely be small in low viscosity mafic melts. Therefore, a different justification is needed to explain these data. An alternative model has been suggested by Farina et al. (2014) and Tang et al. (2014) where xenocrystic zircons dissolve in a melt to produce local regions with Hf isotopic compositions that are a mixture between the original melt and the isotopic composition of the dissolved zircon (the composition is skewed toward that of the Hf rich dissolved zircon). The diffusion of Zr in a melt is relatively slow (Zhang 2010), and therefore, newly formed zircon in these local regions use the Zr from the dissolved grain (without retaining any of the radiogenic Pb) and also obtain a Hf isotopic composition close to that of the xenocryst rather than that of the primary melt. In the case of the NMB, a large portion of the zircons in the sediment through which the basalt intrudes had juvenile Hf (see Waldron et al. 2019) associated with island arc magmatism from terrains accreted onto Laurentia during the Appalachian orogen. The dissolution processes described by (Farina et al. 2014) could succinctly explain the variation toward positive εHf values without affecting the zircon age. This process was likely minor in the NMB magma since only two grains apparently record anomalous εHf values.

Magmatic processes such as the presence of isolated, chemically distinct magma batches, and the potential preservation of xenocrystic cores, are often found for zircon from felsic magmatic systems (Wotzlaw et al. 2014; Samperton et al. 2015; Schaltegger and Davies 2017). The data from the NMB suggest that these are also occurring in mafic magmas. However, the scale of these heterogeneous domains is likely much smaller (as indicated by the melt pockets in Fig. 6), especially since zircon is crystallizing so late (see Fig. 7).

Interpretation of LIP zircon U–Pb data

All of the data presented here indicate that the CAMP magmas were generally highly zircon undersaturated, and those that could crystallize zircon without help from assimilation required > 80% fractional crystallization before they reached zircon saturation, usually at a temperature of ~ 800 °C. Other CAMP samples required significant upper crustal assimilation to reach zircon saturation, and even then, this only occurred after similar amounts of fractional crystallization. This suggests that antecrystic zircon is highly unlikely to be present in CAMP (and LIP) samples, and U–Pb ages from LIP zircon should therefore reflect the final stages of crystallization, and be very close to the age of emplacement (or eruption). Following this logic, when a U–Pb dataset from LIP zircon produces an age distribution (not all of the ages overlap—see Fig. 2), the oldest ages in the distribution should be the closest to eruption or emplacement (given the above zircon saturation constraints). The younger portion of the distribution may be more likely to be influenced by Pb loss, especially given the extremely high U contents that can be present in some LIP zircon (Fig. 5b). This interpretation can be tested using the CAMP U–Pb ages produced here and a simple thermal model of a CAMP sill (the thermal model parameters used to solve the one-dimensional heat conduction equation are in the supplementary information). Samples RP136 and RP134 both have zircon age ranges of ~ 160 ka (not including the much younger zircon from RP134). If a cooling CAMP sill stays within the zircon crystallization window for ~ 160 ka (we define the crystallization window as between 900 and 700 °C), then the zircon age range reflects crystallization; if the sill cools much quicker than this, it is likely that the younger ages reflect small amounts of Pb loss. Figure 10 shows two cooling paths for 150 m-thick CAMP sills in the Amazonian basin, one which intrudes at 3 km depth and the other at 0.5 km depth. The thickness and depths were obtained from the maximum and minimum emplacement depths and the maximum thickness for CAMP sills found within drill core from the Amazon basin (Heimdal et al. 2018). Even at 3 km depth, a 150 m-thick sill should cool below the zircon crystallization window within ~ 500 years, which is ~ 300 times shorter than the measured zircon age range. As discussed briefly above, there are various possible reasons for the zircon crystallization temperatures to be inaccurate (unknown activities of Si and Ti, and also possible kinetic effects); however, changing the zircon crystallization window by tens of °C would not be enough to explain the measured zircon age ranges. It is possible that the thermal model significantly underestimates the geothermal gradient (25 °C/km) in the basin, due to the multiple sill intrusions occurring over ~ 300 ka; however, it seems unlikely that the temperature would be elevated enough to maintain the magma above its solidus for ~ 300 times longer than the model currently predicts. Based on our thermal modeling, and zircon crystallization constraints discussed in detail above, we prefer to interpret the older part of the age distribution as the emplacement and crystallization age of the sills. This method of interpreting zircon age distributions is opposite to those used in felsic rocks (intrusive or extrusive) (e.g., Keller et al. 2018), where the younger portion of the distribution is considered to be closest to the eruption or emplacement age.

Fig. 10

Results from one-dimensional thermal models showing the cooling of a 150 m-thick sill intruded at 0.5 and 3 km. Note that, in both cases, the maximum temperature drops below the zircon saturation window in less than ca. 500 years

Unfortunately, the argument above does not seem to be valid for the case of the NMB data, and neither explains the presence of xenocrysts in the Orange Mountain basalt and AN733. The NMB and Orange mountain basalt are both basalt flows, whereas all of the other samples are intrusives, which could indicate a much shorter time between assimilation and freezing of the magma for the aerially exposed basalts, hence the preservation of xenocrysts. We, therefore, argue that when complex age distributions are found in mafic rocks, combinations of elemental, isotopic and U–Pb data, along with modeling, may be required to correctly interpret the data.

Conclusions

Zircon crystallizes in late-stage enriched melt pockets in LIP magmas. It is likely that crustal contamination of the melts is required for zircon to crystallize in many samples of the CAMP and also for LIP magmas in general. Contamination by upper crustal sediments such as shales most easily creates the conditions for crystalizing zircon by increasing the H₂O content and causing zircon saturation to be reached earlier. Zircon ages obtained for LIP rocks, therefore, date evolved portions (pockets) of the magma at high crystallinity. For high-precision U–Pb dating of LIP magmas, this suggests that the older ages of an age distribution should be close to the time of emplacement or eruption, with the younger part of the distribution more likely influenced by Pb loss, especially with the high U concentrations that can be found in LIP zircon crystals. The rare, anomalously old ages found in the NMB are most likely caused by extremely tiny xenocrystic zircon cores from the assimilated sediments. Mass balance considerations suggest that such inherited cores are of the order of a single µm, hence, extremely difficult to image: Hafnium isotope data can help to confirm their presence if the εHf value of the xenocryst is drastically different from the signature of the LIP magma. Complex zircon age and εHf distributions appear to be rare within the CAMP dataset (including Schoene et al. 2010a; Blackburn et al. 2013; Davies et al. 2017; Heimdal et al. 2018), and therefore, the processes occurring in the NMB may be an exception rather than the rule. The NMB zircons are also distinctive in other ways; for example, their high Ti temperatures and extremely high U and trace-element enrichments, which may allow these types of samples to be identified in future studies. Finally, the heterogeneities that we have documented here indicate that in mafic systems, small, isolated, more evolved melt patches may exist, similar to those found in well-known examples from felsic systems.