Introduction

Zircon is present in a wide range of geological environments and resists both chemical weathering and physical abrasion. Its trace element and isotopic data can provide information critical to understanding Earth’s processes (Finch and Hanchar 2003). Igneous zircons are important reservoirs for U, Th, Hf, and the rare earth elements (REE). These elements allow the identification of fundamental petrogenetic processes or are the key components in relevant isotopic systems (Hoskin and Schaltegger 2003 and references therein). Moreover, the incompatibility of Pb in zircon structure, and the relatively high abundance of U and Th, make it the mineral of choice for U-Th-Pb geochronology (Hoskin and Schaltegger 2003; Ireland and Williams 2003). Finally, the concentrations of Ce, Eu and Ti in zircon can provide valuable insights into petrogenetic processes such as the Watson et al. (2006) and Ferry and Watson (2007) Ti-in-Zr-geothermometer or the oxybarometer of Loucks et al. (2020).

One of the challenges encountered during the analyses of REE in zircon is the ability to accurately measure the light elements in the REE spectrum, which are incompatible with the mineral structure and are present at concentrations close to or below instrumental detection limits (e.g., Burnham 2020). Also, elements of lesser interest, for instance, Tb, Tm, and Ho, may be omitted during analysis by LA-ICP-MS, to allow more counting time for trace elements of interest, such as Ti for thermometry and/or isotopes, such as U and Pb isotopes for zircon dating (e.g., Zhu et al. 2020). This is especially the case where the Laser Ablation Split Stream (LASS) technique is used to simultaneously date the zircons by the U–Pb method, and measure selected trace elements and Hf isotopes. Here, keeping the number of trace elements to a minimum is essential, if only half or less than half of the gas stream is directed through the LA-ICP-MS, so as not to compromise the precision of the most important trace elements and isotopes.

When REEs are measured in zircons it is possible to calculate Ce and Eu anomalies, which require accurate measurement of the REE adjacent to these elements. They are traditionally calculated using the ratio of the element (Ce or Eu) and the geometric mean of the two adjacent elements normalized to their chondritic value. Lanthanum and Pr, for example, are essential for the calculation of Ce anomalies, but their concentrations in zircon are low, usually close to or below the LA-ICP-MS detection limit and, where they are measured, the uncertainty is usually high (Claiborne et al. 2018; Zou et al. 2019; Zhong et al. 2019). Furthermore, their measured concentrations can be compromised by minute inclusions of apatite and other light REE (LREE) rich minerals, which are often overlooked (Burnham 2020; Dilles et al. 2015; Loader et al. 2017; Zhong et al. 2018, 2019; Zou et al. 2019). REEs, Ce and Eu anomalies in zircon are commonly used in geochemistry to monitor changes in redox state and/or the influence of co-crystallizing phases during magma evolution (Ballard et al. 2002; Cocker et al. 2015; Loader et al. 2017, 2022; Loucks et al. 2018).

The calculation of Ce/Ce* ratios is often avoided due to the analytical challenges mentioned above and several methods have been proposed to overcome the analytical problems. The traditional method considers the geometric mean between LaN and PrN to obtain Ce* [Ce* = (LaN × PrN)0.5 where N stands for chondrite normalized values]. Loader et al. (2017) suggested that the geometric mean of NdN = (Ce* × SmN)0.5 can be used to calculate Ce* without requiring the measurement of Pr and La and argue that it can provide a more robust estimation of Ce* because Nd and Sm concentrations can be accurately measured in zircon. Both the La-Pr or Nd-Sm method assumes that the REE pattern of zircon is linear, which can over or underestimate the value of Ce*, due to the curved shape of the zircon REE pattern. Lu et al. (2016) suggested Ce/Nd ratios as a proxy for Ce anomalies, but Ce/Nd and Ce/Ce* are poorly correlated with an R2 of only 0.46. Zhong et al. (2019) suggested a logarithmic regression between the REE atomic number and chondrite-normalized REE values in zircon, which might provide a better approximation. Burnham (2020) reviewed the challenges associated with analysing REE in zircon and suggested that La, Ce* and Pr could be obtained from a lattice strain theory regression. Loader et al. (2022) were the first to use lattice strain theory to calculate Ce* from zircon-glass pairs and compared the results with the methods mentioned above but used Ce/Nd ratios to facilitate comparisons with other published data, and showed a higher correlation between the Ce/Ce* and Ce/Nd than Lu et al. (2016). The multiple methods that have been suggested lack consistency in the literature which makes a comparison of Ce anomalies, calculated by the various methods, difficult.

Most REEs have a 3 + charge in geological processes, which allows them to be fitted with the lattice strain theory (Brice 1975; Blundy and Wood 1994) or to Onuma diagrams (Onuma et al. 1968). The lattice strain theory describes a linear relationship between the logarithm of the partition coefficient of an element, its size, and the elasticity of its preferred lattice site in a crystal, whereas Onuma diagrams describe a quadratic relationship between partition coefficients and ionic radius. In this paper, we present two new methods to calculate REE + Y, Ce* and Eu*, based on the lattice strain theory and Onuma diagrams, but taking advantage of an empirical observation that chondrite-normalized REE values also follow a linear and quadratic fit, respectively. This allows us to calculate REE + Y, Ce* and Eu* for detrital zircons, and for zircons whose host rock composition is unknown, whereas conventional lattice strain theory and Onuma diagrams require prior knowledge of the composition of the melt from which the zircon crystallized to be able to obtain the required partition coefficients. Furthermore, Zhong et al. (2021) have shown that there are up to three orders of magnitude of disagreement between natural and experimental REE partitioning. They suggest that partition coefficients from natural rocks are more reliable but also highlight the difficulty of obtaining true melt compositions in igneous rocks. In the following sections, we will compare both methods, the lattice strain theory and Onuma diagrams, using partition coefficients and chondrite normalized values and compare them with the methods proposed in previous studies.

Ionic radius, Onuma diagrams and parabolas

Onuma et al. (1968) showed that for a set of isovalent cations, the logarithm of their partition coefficients has a quadratic relationship with their ionic radius. The vertex (and highest partition coefficient) of the parabola is the optimal size for that site in the crystal lattice for an element with a given charge. For example, in zircon (ZrSiO4), the zirconium site in the lattice has eight-fold coordination. The optimal size for an ion with a 4 + charge is 0.84 Å (Shannon 1976) but most estimates for cations with a 3 + are close to 0.92 Å (Burnham 2020; Colombini et al. 2011; Loader et al. 2022; Zhong et al. 2019). For the REE with a 3 + charge, partitioning into the zircon in the eight-fold coordination site (and excluding Ce and Eu that have multiple oxidation states), the partition coefficient describes a quadratic function that decreases as the ionic radius of the element becomes higher or lower than the optimal ionic radii for 3 + cations. Elements with a 4 + charge in zircon (e.g., Ti, Th, U, Hf) also describe a quadratic function with a higher aperture (a higher partition coefficients) but with the inflexion point at the Zr partition coefficient and ionic radius.

Figure 1 shows the comparison of the ratio of zircon to whole-rock REE compositions as a proxy of partition coefficients (Fig. 1A) and chondrite-normalized values compared with their ionic radii (Fig. 1B) for three zircons from the Chuquicamata-El Abra District, Northern Chile (Ballard et al. 2002). A quadratic model was fitted for each grain from Nd to Lu. Cerium and Eu were excluded because they can be present in more than one oxidation state, and their position will lie between the 3 + and 4 + parabolas and the 2 + and 3 + parabolas, respectively. Lanthanum and Pr concentrations in zircon usually have a high uncertainty or can be over-measured due to LREE-bearing inclusions in zircon (e.g., apatite, allanite; Zhong et al. 2018). Yttrium was also excluded because its measured values do not fit the quadratic regression as well as the REE3+. Over 99.99% of the variance in the data is accounted for by both the partition coefficients and the chondrite-normalized values (R2 > 0.9999). Because the observed quadratic relationship is with the logarithm of the partition coefficients, minor deviations, lack of data, or a misfitted observation, can heavily affect the interpolations, especially towards the borders of the parabola and can lead to unrealistic high or low values, especially for La and Lu.

Fig. 1
figure 1

Onuma diagrams for three zircons from Ballard et al. (2002), two from the El Abra and Chuquicamata porphyry copper deposits and one from a Los Picos barren intrusion. We have used the whole rock composition as a proxy of the melt composition to estimate the partition coefficient. Equations for each parabola are colour-coded, where y = log10(x). The X-axis is reversed (ionic radius decreases to the right). The shading indicates a 95% confidence interval for the quadratic fit. Lanthanum, Ce, Eu and Y are excluded from the regression. A Partition coefficient vs ionic radii. B Chondrite normalised values vs ionic radii. Chondrite values were taken from Palme and O’Neill (2014)

Full size image

The lattice strain theory

The lattice strain theory was proposed by Brice (1975) to explain the incorporation of misfit cations during the growth of doped synthetic crystals. Blundy and Wood (1994) applied this theory to geological processes. According to the lattice strain theory, the mineral-melt partition coefficient (Di) for a cation, which misfits into a given crystal lattice, is related to the Gibbs free energy (ΔG) required to introduce it into the lattice as:

$${D}_{i}={D}_{0}e^{-\Delta G/RT}$$
(1)

where D0 is the hypothetical strain-free partition coefficient for the crystal lattice site (e.g., D0 = Zr4+ and Di = Ce4+ in a zircon crystal), R is the gas constant, and T is the temperature in Kelvin. The ΔG is controlled by the difference between the ionic radius of the cation ri and the strain-free site r0 through Young’s modulus (E), which describes the linear elasticity between the stress and strain of the host crystal as:

$$\Delta G=4\pi E{N}_{A}\left[\frac{{r}_{0}}{2}{\left({r}_{i}-{r}_{0}\right)}^{2}+\frac{1}{3}{\left({r}_{i}-{r}_{0}\right)}^{3}\right]$$
(2)

where NA is Avogadro’s number (to transform atoms into moles). Replacing ΔG from Eq. 1 in Eq. 2 gives:

$${\text{ln}}\left({D}_{i}\right)={\text{ln}}\left({D}_{0}\right)-{N}_{A}\frac{4\pi {\text{E}}}{\text{RT}}\left(\frac{{r}_{i}}{3}+\frac{{r}_{0}}{6}\right){\left({r}_{i}-{r}_{0}\right)}^{2}$$
(3)

The terms 4πE/RT and NA are constants and, for any Di and ri of a lattice site, D0 and the optimal size r0 will be constant. Therefore, ln(Di) will depend on the difference between the size of the lattice site and the ionic radii of the cations, which have a linear relationship in which the charge of isovalent cations defines the slope. In zircons, the slope for 4 + cations will be different from the slope for 3 + cations. Elements that can be present in more than one oxidation state, like Ce and Eu, will fall between the lines for their oxidation states in the proportion they enter into the crystal lattice.

Figure 2 shows the comparison of the ratio zircon to whole-rock REE as a proxy for partition coefficients (Fig. 2A) and chondrite-normalized REE (Fig. 2B) plotted against (ri/3 + r0/6)(ri − r0)2 for the same zircons as in Fig. 1. The excluded elements in each model are the same as in Fig. 1. In both cases, the relationship is log-linear, with R2 values over 0.9999 for both the partition coefficients and chondrite-normalized values.

Fig. 2
figure 2

Lattice strain model for three zircons from Ballard et al. (2002), two from the El Abra and Chuquicamata porphyry copper deposits and one from a Los Picos barren intrusion, represented by different colours as in Fig. 1. We have used whole rock composition as a proxy of the melt composition to estimate the partition coefficients. Equations for each linear regression are colour-coded, where y = log10(x). The X-axis is reversed. The shading indicates the 95% confidence interval for the linear regression. A Partition coefficients vs (ri/3 + r0/6)(ri-r0)2, the best fit is obtained from an average r0 = 0.91 Å. B Chondrite normalised values vs (ri/3 + r0/6)(ri-r0).2, the best fit is obtained from an average r0 = 0.84 Å. Chondrite values are from Palme and O’Neill (2014)

Full size image

The equations illustrated in Figs. 1 and 2 can be used to determine the REE and Y concentrations by taking the ionic radius of each element as the variable x in the Onuma diagrams or ri in the lattice strain model. If the Ce3+ and Eu3+ in eight-fold coordination ionic radius are used in the equation (Ce = 1.143 Å, Eu 1.066 Å; Shannon 1976), the values obtained will correspond to Ce* and Eu*, respectively, which can be used to calculate Ce and Eu anomalies.

We compiled a zircon database with no missing data for the range of REE from La-Lu and Y. The database has been used to build regression models to calculate the REE using the Onuma diagram parabola (Chondrite-Onuma method) and the lattice strain theory (Chondrite-Lattice method), both using chondrite-normalized values instead of partition coefficients. The calculated values were then compared to the measured values, and other methods, for determining REE concentrations and Ce* and Eu* values.

Methods

Data compilation

Nearly 1500 zircon grains, which were analysed for the complete range of REE + Y, were compiled from the literature. Only datasets from publications that used a screening method, such as cathodoluminescence (CL) imaging or a geochemical technique, to exclude zircons with mineral inclusions (e.g., apatite, titanite), were considered. The data sources are summarized in Table 1. We will refer to this dataset as the main dataset.

Table 1 References for the data used in this study
Full size table

A second dataset containing analysis of experimental (n = 6; Taylor et al. 2015) and natural zircon-matrix glass pairs (n = 104; Colombini et al. 2011; Claiborne et al. 2018) was used to compare our methods with lattice strain estimates following Loader et al. (2022). In order to prevent any confusion with the main dataset, we will refer to this particular collection as the partition coefficients dataset.

The relationship between REE rock and chondrite compositions

The observation that chondrite-normalized values and partition coefficients for REE produce almost identical lines when used with Onuma the lattice strain theory suggests that they are correlated (Figs. 1 and 2). However, the trace element concentration in zircon is the numerator for both variables, which results in a spurious correlation (Fig. A1, Electronic Supplementary Material (ESM) 1; Rollinson and Pease 2021). Therefore, the extent to which using lattice strain or Onuma diagrams using partition coefficients differs from using these methods with chondrite-normalized values, depends on the difference between the melt and chondrite compositions. To test whether this is a significant factor, we evaluated the correlation between chondrite compositions and glass compositions using the zircon partition coefficients dataset. Figures A2 and A3 (ESM 1) show the REE concentrations from the partition coefficient dataset plotted against the REE chondrite composition of Palme and O’Neill (2014). For each of the glass-chondrite REE pairs is possible to fit a linear regression in the form of log10(REEchondrite) = log10 (REEglass) * m + b, where m and b are constants that represent a slope and intercept that is specific for each sample. Because the chondrite REE composition is constant, the regression parameters depend solely on the composition of the melt. Figure A2 (ESM 1) shows that all of the samples can be fitted with a linear regression with an R2 between 0.68 and 0.93. Figure A2 also shows that La and Pr fall off the regression trend and could be excluded bearing in mind that the analysis of these elements in zircon is usually compromised by contamination (e.g., Zhong et al. 2018). After exclusion of La and Pr, the R2 of the linear regression between melt and the chondrite composition increases from 0.68–0.93 to 0.87–0.98 (Fig. A3, ESM 1).

Figure 3 summarizes this relationship as the distribution of the Pearson correlation coefficient between the natural logarithm of the glass-chondrite concentration pairs for the samples in the partition coefficient dataset. All of the pairs have a positive correlation within the range of 0.82 to 0.99 for the full range of trivalent REE (Fig. 3A), which narrows to 0.92–0.99 when La and Pr are excluded (Fig. 3B).

Fig. 3
figure 3

Kernel density estimate of Pearson correlation coefficient between the natural logarithms of the trivalent REE concentrations in the glasses from Colombini et al. (2011), Taylor et al. (2015), Claiborne et al. (2018), using the chondrite of Palme and O’Neill (2014) from A La to Lu and B Nd to Lu. Ce and Eu are excluded in both cases

Full size image

The correlation observed between the chondrite and glass REE composition shows why the lattice strain theory and Onuma diagrams can be successfully applied to chondrite normalized values for the Nd to Lu range. However, the data presented in Fig. 3 only covers a limited range of compositions and may not be generalizable to a broader range of melt compositions.

To test if this empirical observation is present in other melt compositions, we have used precompiled files from the GEOROC database (DIGIS Team 2022a, b, c, d, e, f, g, h, i; Lehnert et al. 2000) for different geological settings, including the Andean and Tonga Arcs (Fig. 4A and B), the Karoo-Ferrar large igneous province (Fig. 4C), the East African Rift (Fig. 4D), the Western Australian Craton (Fig. 4E), and the Hawaiian ocean island basalts (Fig. 4F). Each dataset was filtered for igneous rocks with measured concentrations for the whole range of REE. It is important to recognise that GEOROC data do not necessarily represent a melt composition.

Fig. 4
figure 4

Person correlation coefficient between the natural logarithms of the chondrite values and rock composition for the pre-compiled dataset from the GEOROC database, filtered for samples with no missing REE data. A reference list of the original publications is in ESM 3. N = number of observations, r > 0.8 = X% indicates the proportion of the data with a Pearson r coefficient higher than 0.8. see text for details. A Andean arc. B Tonga Arc. C Karoo-Ferrar Large igneous province. D East African Rift. E West Australian Craton. F Hawaiian ocean island basalts

Full size image

Figure 4 shows the distribution of the Pearson correlation coefficients between chondrite and whole-rock composition for each of these geological environments. Most of the whole rock–chondrite pairs have a positive correlation > 0.8 for the Nd-Lu REE range, which suggests a relationship between the REE concentrations in rocks and the chondritic values. This correlation could be the result of the REE behaving as an incompatible group, which leads to smooth chondrite normalized patterns, especially towards the heavy REE end of the spectrum where the ionic radius differences are smaller, which can be observed in typical REE patterns of continental or oceanic crust (e.g., Rudnick and Gao 2014; White and Klein 2014).

A possible cause of concern, when fitting linear REE models to whole rock data, where the REE concentrations of the rock are unknown, is the assumption that the correlation between rock and chondrite compositions is valid. However, a poor correlation would result in poor regression fit and the analysis being discarded.

Disagreement between predicted and measured concentrations

We evaluated the accuracy of the Chondrite-Lattice and the Chondrite-Onuma methods for calculating missing REE and Y data by comparing the calculated values with the measured values for each trivalent REE from Nd to Lu and Y on the main dataset.

Cerium, Eu, Y, La and Pr were excluded because they do not fit the regression obtained from either model. Ce and Eu can occur in more than one oxidation state in magmas, therefore they would not fit the 3 + curves for either method.

Yttrium is trivalent, with an ionic radius similar to Ho and is it therefore expected to be Ho’s elemental twin. However, Y’s chondrite-normalized values deviate from both lattice strain and Onuma fitting models, as seen in Figs. 1 and 2 where Y concentrations fall slightly above the regression line with chondrite normalized values closer to Er. The reason for this difference is unclear but it reflects preferential partitioning of Ho over Y in HREE-compatible minerals. Magmatic zircons, hornblende and titanite, for example, have Ho partition coefficients that are up to 19%, 14% and 47% higher than for Y, respectively (Colombini et al. 2011; Claiborne et al. 2018; Nandedkar et al. 2016).

Lanthanum is highly incompatible with the zircon lattice. Modelling of partition coefficient and REE concentrations indicates that the zircon lattice can accommodate up to 2 to 100 ppb La (Burnham 2020; Claiborne et al. 2018; Zhong et al. 2019; Zou et al. 2019). This makes zircon crystals very sensitive to micro inclusions of LREE-rich minerals. For instance, Burnham (2020) found that 0.005% allanite contamination is enough to raise zircon La concentrations to 1.5–2.5 ppm. Therefore, La was also excluded from the modelling. Similarly, Pr which is also present in low concentrations in zircon and can be susceptible to LREE-rich inclusions and was also excluded (Burnham 2020; Zhong et al. 2018; Zou et al. 2019).

Figure 5 shows boxplots of the ratio between the calculated and measured concentrations in zircon for the REE + Y from Nd to Lu, ordered by ionic radius. Overall, both methods show a small disagreement between measured and calculated values. This effect is more prominent for the Chondrite-Lattice model (Fig. 5A) for which the boxplots are wider and are farther from the ratio = 1 line, especially towards the extremes of the REE spectrum. In contrast, the results from the Chondrite-Onuma (Fig. 5B) are less spread, as shown by narrower boxplots, and closer to the calculated/measured ratio of 1 for the whole REE spectrum. Both methods underestimate Y. This disagreement can be corrected if it is assumed that the disagreement is constant. This process is known as model calibration, where the predicted values are adjusted to the actual measurements (Kuhn and Johnson 2013). The Chondrite-Lattice strain method underestimates Y, Lu, Yb, Tb and Nd, in median, by 29%, 22%, 14%, 11% and 11%, respectively. Tm and Er are underestimated by nearly 1% whereas Sm to Ho are overestimated by 1–7%. A correction factor using the reciprocal of the median value of the difference can be applied, which is equivalent to centring each boxplot to the calculated/measured ratio = 1 line in Fig. 5 (Fig. A4 ESM 1). The deviations for the Chondrite-Onuma methods are under 10% for all REE, except Y which gives calculated values that are 44% lower than the measured concentrations. The correction factors for both methods are listed in Tables A1 and A2 (ESM 1).

Fig. 5
figure 5

Boxplots for the ratio between calculated and measured concentrations for each REE, excluding La, Pr, Ce and Eu. Calculated values will plot along the y = 1 line if they are equal to the measured values. The blue segmented lines, blue continuous lines and red lines are reference lines for a discrepancy between measured and calculated values of 10%, 50% and 100%, respectively. Each box contains 50% of the data. Each box and whiskers represent 99.3% of the data. The grey dots are outliers and correspond to 0.7% of the data. A Chondrite-Lattice strain method iterating for an optimum r0 in each zircon grain following Loader et al. (2022). B Chondrite-Onuma method

Full size image

After corrections, more than 50% of the calculated values from Nd to Tm and Y have discrepancies smaller than 10% when compared to the measured concentrations when using the Chondrite-Lattice method. Furthermore, for most REE, nearly all the data show a disagreement of less than 50% when compared to the calculated values, except for Pr and Lu (Fig. A4). In contrast, nearly all the calculated values have a disagreement of less than 10% after correction with the Chondrite-Onuma method (Fig. A5).

The source of these differences is uncertain but could be due to (i) poor fitting between the LREE and HREE in the model as would be characteristic for a robust model like a linear regression, where overestimated values fall below the regression line and underestimated values above; (ii) uncertainty in the chondrite values used to normalize the REE concentrations. For instance, if chondritic values of McDonough and Sun (1995) are used, Y, Yb and Lu values are underestimated by 21%, 15% and 20%, respectively, for the Chondrite-Lattice method, which shows that the model results vary with the choice of chondritic values; and (iii) for the Chondrite-Lattice method the r0 value might not be well constrained due to the differences between chondritic values and true melt composition.

We will present the results for the chondrite Onuma method, as it is a simpler alternative to the Chondrite-Lattice method and because it provides more accurate results without requiring corrections. The results on the Chondrite-Lattice method, including the selection of the best fitting r0 and an analysis of misfitting models, and a discussion of its applications to other minerals can be found in the ESM 4.

Results

The La and Pr concentrations

As previously stated, La and Pr can be subject to contamination. We have compared the calculated values vs the measured values for both elements. The results for the Chondrite-Onuma method (Fig. 6A and B), give La values that are nearly 35% lower than the measured concentrations, with a median of 19 ppb, which is only 1.5 times lower than the measured average, whereas Pr calculated concentrations are similar to the measured concentrations (Fig. 6B). There are zircons in the dataset with La and Pr contents as high as 100 ppm which have a poor model fitting (R2 < 0.9 for Chondrite-Lattice, R2 < 0.98 for Chondrite-Onuma). These high LREE zircons are from the Bajo de la Alumbrera (Buret et al. 2016) and the high concentration may be due to unidentified or missed inclusions during data screening. The discrepancy between measured and calculated values for La and Pr, in addition to observations by other authors (e.g., Claiborne et al. 2018; Zhong et al. 2019; Zou et al. 2019) suggests, in most cases, it is better to calculate the La concentration than to measure it.

Fig. 6
figure 6

A Measured La vs calculated La and B measured Pr vs calculated Pr using the Chondrite-Onuma method. The colour represents the coefficient of determination (R2). The continuous black line is the 1:1 ratio between the calculated and measured values. Each segmented line indicates a 0.5 order of magnitude step

Full size image

Model testing

The objectives of the method are to (i) calculate REE and Y data where data are missing, and (ii) calculate Ce* and Eu*. Therefore, it is necessary to evaluate how well these methods predict REE concentrations when data are missing. Six different scenarios were considered in which one or more REEs were excluded during modelling to represent missing data. The results are summarized in Fig. 7A–F for the Chondrite-Onuma method. Each scenario is described below:

  1. 1.

    Calculations based on three light REEs, Nd, Sm and Gd (Fig. 7A): This scenario yields good results for the light REEs but becomes increasingly less reliable as the REE becomes heavier. We have modelled this scenario for illustration purposes only to show that an interpolation using only an extreme from the REE pattern should be avoided.

  2. 2.

    Calculations based on three heavy REEs Ho, Tm and Lu (Fig. 7B): As with the previous case and it produces unrealistic values at the other extreme of the REE pattern and is not recommended.

  3. 3.

    Calculations based on three evenly distributed REEs, Nd, Dy and Tm (Fig. 7C): These models show that for Nd to Lu, the differences between calculated and measured data, after correction, are minimal.

  4. 4.

    Calculations based on four evenly distributed REEs, Nd, Sm, Dy, and Lu, as in the river zircon database from Zhu et al. (2020) (Fig. 7D): The addition of an additional REE into the calculation (Fig. 7E, F), makes little difference to the accuracy when compared with the results based on 3 REE.

Fig. 7
figure 7

Figures A–F show boxplots of the ratio between the calculated and measured values using the Chondrite-Onuma method and applying the reciprocal of the median as a correction. The REE used in each regression is at the top of each plot. The boxes in a perfect regression would lie along the central line, with the size of the box representing the spread of calculated concentrations. The blue dashed line, blue continuous line and red line indicate a difference of 10%, 50% and 500% between calculated and measured values, respectively. Each box contains 50% of the data. Each box and whisker include 99.3% of the data. The dots indicate outliers and account for 0.7% of the data

Full size image

Our analysis shows that 3 or more REEs distributed over the full REE spectrum can be used to obtain reliable calculated values for REEs that were not analysed. We suggest the analysis of at least five REEs: Ce, Eu, Nd, Dy, and a heavy REE (Tm, Yb or Lu) when it is desirable to reduce the number of REEs during ICP-MS analyses. If more REEs are considered for analysis, Lu or Sm should be included since the model, shows a slightly higher uncertainty for these elements (Fig. 5B). However, the use of this method, as an imputing technique, should only be used if the element of interest was not analysed, or if the measured value is judged to be unreliable as is usually the case for La and Pr.

Discussions

Comparison to other methods to calculate REE

There are only two other methods that do not require knowledge of the rock or melt compositions to calculate REE in zircon. One approach is to use the geometric mean from adjacent REE, as Ce* and Eu*, are traditionally calculated. However, this method is only an approximation if the pattern is curved, as it is for zircon, and cannot be applied when data are missing for more than two consecutive REEs. The other method is from Zhong et al. (2019), which uses a logarithmic regression between the logarithm of the chondrite normalized values and the atomic number Z. We have compared Zhong et al. (2019) regression to our method using the main dataset (Table 1).

Figure 8A compares the coefficient of determination (R2) between the calculated and measured values for each REE using the method in this study, with and without corrections, with the logarithmic regression of Zhong et al. (2019). The highest R2 for most REE are obtained using the Chondrite-Onuma regression, except for Sm where Zhong et al. (2019) method gives slightly better results.

Fig. 8
figure 8

A Coefficient of determination for each pair of measured vs calculated REE in the main dataset. Because the median correction shifts each element by a constant amount it does not change the R2 value so the lines of median-corrected models fall exactly on those that are not corrected. B Normalized root mean square error for each pair of measured vs calculated REE, normalized by the mean of each element. Values of R2 closer to 1, and lower values of RSME, indicate a closer match between measured and calculated REEs. The legend gives the reference, the method and corrections if applied

Full size image

The coefficient of determination is a measure of how well the regression line fits the data. We have also calculated the root mean square error (RMSE) normalized to the mean value of each element, which is a measure of the average deviation of the predictions from the actual values. Lower values of RMSE indicate that the calculated values are closer to the measured values. The RMSE obtained between the calculated and measured chondrite normalized values for different models are shown in Fig. 8B. As with R2, the best models use the Chondrite-Onuma method, and the poorest results are from the Zhong et al. (2019) model, except for Sm. There is not much difference between corrected and uncorrected models, as shown in Fig. 8A and B, except for Y and Tm, for which their RMSE is reduced by half after corrections.

The method of Zhong et al. (2019) is also subject to the deviations seen in Fig. 5 for the Chondrite-Onuma and Chondrite-Lattice methods, as illustrated in Fig. A6 (ESM 1). Therefore, similar corrections could be applied. However, the method proposed by Zhong et al. (2019) cannot impute Y concentrations and gives less consistent results than the method we propose. We recommend using the Chondrite-Onuma method to impute missing trivalent REE in the range Nd-Lu and for Y concentrations.

La and Pr are unique in that they are susceptible to contamination, as previously discussed. All the methods predict similar Pr concentrations that are lower than the measured concentrations.

The La results are more strongly method dependent than Pr. All methods predict La concentrations that are lower than the measured concentrations. However, the method of Zhong et al. (2019) predicts values that are 1 to 2 orders of magnitude lower than the Chondrite-Onuma method. The median La calculated by the Zhong et al. (2019) method is 0.34 ppb, more than 100 times lower than the median measured concentration. However, it is challenging to identify which method would provide the most accurate estimates for La, Pr, and Ce* in zircons, as there are no reliable concentration measurements to compare them with. Different authors have proposed different thresholds for determining whether a zircon is “clean”, or inclusion free. For example, (Zou et al. (2019) suggest that a zircon should contain less than 100 ppb of La to be considered clean, while Burnham (2020) suggests that even zircons with < 10 ppb of La may still have “excess” La. Claiborne et al. (2018) derived partition coefficient data using the lattice strain theory and concluded that zircon can accommodate no more than 2 ppb La, which agrees with Zhong et al. (2019), who propose a threshold of less than 3 ppb. These authors suggest that the concentration of La in zircon should be lower than most measured values.

Therefore, because the concentration of La and Pr in zircon cannot be accurately measured, we have used the best estimate based on lattice strain theory. The values for La and Pr in this case can be obtained by using their ionic radii as ri in the equations from Fig. 2A to obtain the partition coefficient which is then used to estimate the concentrations of these elements in the zircon. Figure 9A and B show the concentrations of La and Pr, respectively, obtained by using the lattice strain theory (from partition coefficients) and iterating for the optimal r0 (following Loader et al. 2022), against the Chondrite-Onuma and Zhong et al. (2019) methods, using the partition coefficient dataset (glass-zircon pairs). It is clear from Fig. 9B that the two methods give Pr calculated values similar to the lattice strain estimates, however, they are slightly better for the Chondrite-Onuma method (lower RMSE), particularly at lower concentrations. In contrast, the method proposed by Zhong et al. (2019) has been found to significantly underestimate the concentration of La in zircon with discrepancies approaching two orders of magnitude when compared to estimates obtained using the lattice strain theory. The magnitude of these differences increases as the value of the La lattice strain estimate decrease. On the other hand, the Chondrite-Onuma method tends to overestimate the La concentrations by about 0.5 order of magnitude, but it is consistent in that it keeps the same distance from the identity line in Fig. 9A. The RSME for the Chondrite-Onuma method is almost three times lower than for the method of Zhong et al. (2019), which makes the Chondrite-Onuma method the most appropriate for estimating La concentrations in zircons.

Fig. 9
figure 9

A La and B Pr concentrations obtained by the lattice strain theory vs the calculated values using the Chondrite-Onuma and the Zhong et al. (2019) methods. The blue line indicates the identity line. The segmented line represents 0.5 orders of magnitude steps from the identity line. RSME are colour coded according to the method. Zircon-glass pairs are from Colombini et al. (2011), Taylor et al. (2015) and Claiborne et al. (2018)

Full size image

One way to calibrate the La and Pr concentrations for each method would be to project them to the identity line (e.g., Fig. 9), using the reciprocal of the median correction that was applied to the other REE. However, we advise against this. The La and Pr concentrations calculated by the lattice strain theory, despite being based on physical principles rather than empirical observations, are still estimates derived from a linear regression, and there are no measured values to use as a benchmark, as is the case for the other methods presented in this study.

A comparison of the different methods to calculate Ce and Eu anomalies in zircon

The Chondrite-Onuma method can be used to calculate Ce* and Eu* for each zircon by employing the regression equations, based on their trivalent ionic radius, as in Fig. 1. This allows for the calculation of Eu and Ce anomalies, which we now compare with other methods in the literature.

The results for Eu/Eu*, calculated using the traditional geometric mean [EuN/Eu*, where Eu* = (GdN × SmN)0.5], are highly correlated with those calculated using the Chondrite-Onuma method (Fig. A7, ESM 1) with an R2 of > 0.99 for the 1486 zircons grains compiled for this study. The regression parameters have an intercept close to 0 and a slope of ~ 1. This indicates that there is no advantage in using either method to calculate Eu anomalies, except when Sm or Gd is not available.

The problems encountered in measuring La and Pr in zircon make Ce anomalies, calculated by the geometric mean, unreliable. As with La and Pr, the best estimates for Ce* are derived from the Lattice strain theory. Figure 10 compares the Ce*, calculated by the different methods, with values calculated by the lattice strain theory for the partition coefficient dataset. The Chondrite-Onuma gives the best results and they are almost identical to those from the lattice strain estimates. In contrast, the method of Zhong et al. (2019), as it does with La, tends to increasingly underestimate Ce* as the lattice strain estimate for Ce* is lower, which can potentially lead to overestimated Ce/Ce* ratios.

Fig. 10
figure 10

Ce* obtained by the lattice strain theory vs the calculated values using the Chondrite-Onuma and the Zhong et al. (2019) methods. The blue line shows the identity line. Each segmented line represents 0.5 orders of magnitude. RSME are colour coded according to the method. Zircon-glass pairs are from Colombini et al. (2011), Taylor et al. (2015) and Claiborne et al. (2018)

Full size image

Conclusions

We have tested two new methods to calculate REE and Ce and Eu anomalies in magmatic zircon based on the empirical observation that chondrite-normalized values can be used, instead of partition coefficients, when applying Onuma diagrams and Lattice strain theory to REE in the range Nd to Lu. This is possible because most rocks preserve an REE signature similar to the CI chondrites between Nd and Lu, as shown by the strong correlation between chondritic values and rock compositions.

Onuma diagrams and Lattice strain theory, using chondrite normalized values instead of partition coefficients, can be used to robustly calculate missing REE and Ce* and Eu* from as few as five REE (Ce, Eu, Nd, Dy and Lu). The principal application of these methods is to calculate the missing REE + Y from legacy data and to estimate La and Pr concentrations and Ce* in zircons. The methods allow more than half of the REE to be omitted during analyses, which gives additional counting time for elements and isotopes of greater interest when analytical capabilities are limited. Based on the results of this study, we recommend the use of the Chondrite-Onuma method, in preference to the Chondrite-Lattice method, because of the simplicity and lower uncertainty of the former.