Application of Onuma and lattice strain derived methods to calculate missing REE and Ce and Eu anomalies in magmatic zircons

The shape of chondrite normalized zircon rare earth element (REE) patterns, including the magnitude of the Cerium (Ce) and Europium (Eu) anomalies, provide valuable insights into the magmatic conditions under which a zircon formed. However, lanthanum (La) and praseodymium (Pr), which are essential for the calculation of Ce anomalies, are often present at concentrations close to or below the detection limit of most analytical methods. We propose two new methods to calculate missing REE, based on Onuma diagrams (Chondrite-Onuma) and the lattice strain theory (Chondrite-Lattice), but using chondrite normalized values instead of partition coefficients. We compiled a dataset of ~ 1500 zircons with known REE + Y concentrations and used it to test and calibrate these methods and demonstrate that they are more accurate than other previously published models, with the Chondrite-Onuma method performing better than the Chondrite-Lattice method. These methods require analyses of as few as five REEs to impute the missing REE data or to estimate La and Pr concentrations or Ce anomalies in magmatic zircons, which allows a reduction in the number of REE analysed, where desirable, or to impute missing REEs in legacy data. The imputeREE package for the R programming language was written with a set of tools to apply these methods. A companion app is available to calculate missing REE and Ce and Eu anomalies.


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. 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. 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 (ZrSiO 4 ), 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 eightfold 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 REE 3+ . Over 99.99% of the variance in the data is accounted for by both the partition coefficients and the chondrite-normalized values (R 2 > 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.

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 (D i ) 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: where D 0 is the hypothetical strain-free partition coefficient for the crystal lattice site (e.g., D 0 = Zr 4+ and D i = Ce 4+ 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 r i and the strain-free site r 0 through Young's modulus (E), which describes the linear elasticity between the stress and strain of the host crystal as: where N A is Avogadro's number (to transform atoms into moles). Replacing ΔG from Eq. 1 in Eq. 2 gives: The terms 4πE/RT and N A are constants and, for any D i and r i of a lattice site, D 0 and the optimal size r 0 will be constant. Therefore, ln(D i ) 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,  Fig. 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) 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 (r i /3 + r 0 /6)(r i − r 0 ) 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 R 2 values over 0.9999 for both the partition coefficients and chondrite-normalized values.
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 r i in the lattice strain model. If the Ce 3+ and Eu 3+ 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.

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.
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 y = 2.69 163 x R 2 = 1.00 y = 2.65 177 x R 2 = 1.00  Fig. 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 = log 10 (x). The X-axis is reversed. The shading indicates the 95% confidence interval for the linear regression. A Partition coefficients vs (r i /3 + r 0 /6)(r i -r 0 ) 2 , the best fit is obtained from an average r 0 = 0.91 Å. B Chondrite normalised values vs (r i /3 + r 0 /6)(r i -r 0 ). 2 , the best fit is obtained from an average r 0 = 0.84 Å. Chondrite values are from Palme and O'Neill (2014) 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 log 10 (REE chondrite ) = log 10 (REE glass ) * 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 R 2 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 R 2 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).
The correlation observed between the chondrite and glass REE composition shows why the lattice strain theory and  (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. 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 LREErich 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).
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 r 0 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 r 0 and an analysis of misfitting models, and a discussion of its applications to other minerals can be found in the ESM 4.

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 (R 2 < 0.9 for Chondrite-Lattice, R 2 < 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. 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 r 0 in each zircon grain following Loader et al. (2022). B Chondrite-Onuma method

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. 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.

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. 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. 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 dif-ference to the accuracy when compared with the results based on 3 REE.
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.

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  (Table 1). Figure 8A compares the coefficient of determination (R 2 ) 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 R 2 for most REE are obtained using the Chondrite-Onuma regression, except for Sm where Zhong et al. (2019) method gives slightly better results.
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 R 2 , 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  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 r i 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 r 0 (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.
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 [Eu N /Eu*, where Eu * = (Gd N × Sm N ) 0.5 ], are highly correlated with those calculated using the Chondrite-Onuma method (Fig. A7, ESM 1) with an R 2 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, 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) 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.

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. editing: Funding Open Access funding enabled and organized by CAUL and its Member Institutions.

Data availability
The data that support the findings of this study were derived from the resources available in the public domain which are listed in Table 1 and available in the Electronic Supplementary Information 2.

Code availability
The code used to calculate REE patterns has been documented in the imputeREE package published in the comprehensive R archive network (CRAN) at https:// CRAN.R-proje ct. org/ packa ge= imput eREE. The package includes functions to apply the Chondrite-Lattice, Chondrite-Onuma and the Zhong et al. (2019) methods to calculate REE and Ce*. A companion application, based on this package to calculate REE was written in R shiny that allows to apply this three methods with an user friendly interface. A link to the source code and a link to the app are available at the development website of the imputeREE package at https:// github. com/ cicar rascog/ imput eREE.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. 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)