Characterization of the zircon Ce anomaly for estimation of oxidation state of magmas: a revised Ce/Ce* method

Abstract

As a proxy of magmatic oxidation state, the accurate characterization of the Ce anomaly of zircon is of great significance since it can give important information for provenance studies of rocks as well as for exploration of intrusion-related mineral deposits. The magnitude of the zircon Ce anomaly has been traditionally described by Ce/Ce*, where Ce* is the theoretical Ce value derived from a chondrite-normalized rare earth element (REE) pattern. More recently, the Ce⁴⁺/Ce³⁺ method based on the lattice strain model has been proposed, since the latter method does not need La and Pr contents for zircon, both of which are commonly below the limit of detection and susceptible to contamination from melt/mineral inclusions. In this contribution we show that the Ce⁴⁺/Ce³⁺ method is confronted with some problems in practice and should be further improved. In contrast, by re-examining chondrite-normalized REE patterns of zircon, we find that Ce* can be estimated according to a logarithmic function curve without involvement of La and Pr contents. Application of this new method to zircon data from 11 giant to supergiant porphyry Cu deposits suggests this revised method as a more valid measure in evaluating magmatic oxidation state. The revised Ce/Ce* method is of particular importance for analyses where the provenance of the analyzed zircon is unknown or in question, since the method does not require knowledge of the melt composition.

Introduction

As one of the most important physico-chemical parameters of magmas (Barberi et al. 1971), there is a consensus that the oxidation state can act as a valuable archive of geochemical information. For instance, I-type granites are normally characterized by higher oxidation state than S-type granites (Blevin and Chappell 1992); therefore, geologists can obtain provenance information of magmatic rocks from the oxidation state of the magmas. Also, studies have demonstrated that oxygen fugacity exerts a fundamental control on magma metal fertility for many intrusion-related mineral deposits (Zhong et al. 2017), which are currently the main sources of a number of economically important metals on earth (Gardiner et al. 2017).

Zircon (ZrSiO₄) is an accessory mineral enriched in REEs that is common in igneous rocks of intermediate to felsic composition. Experimental studies have established that Ce enters zircon in both 3+ and 4+ oxidation states (Trail et al. 2015) and that Ce (as well as other REEs) will only be incorporated into the Zr site in the zircon structure (Finch et al. 2001). Owing to these reasons, in chondrite-normalized REE patterns, zircon is commonly characterized by a positive anomaly relative to neighboring La and Pr which exclusively form trivalent ions. Compared to Ce³⁺, Ce⁴⁺ is preferentially incorporated into the zircon structure due to its identical charge and similar size to Zr. Thus, the degree of the Ce anomaly is largely controlled by the oxidation state of the magma, although the water content of the melt, temperature and melt composition also have an influence (Smythe and Brenan 2015). High oxygen fugacity will oxidize Ce³⁺ to Ce⁴⁺ and permits more Ce to be incorporated into the zircon crystal structure, thus resulting in greater Ce anomalies in chondrite-normalized patterns. This relationship between zircon Ce⁴⁺/Ce³⁺ and magmatic oxygen fugacity has now been firmly demonstrated by experimental studies (Burnham and Berry 2012; Trail et al. 2012; Trail et al. 2011; Smythe and Brenan 2016).

As a proxy for the oxidation state of magmas, a reliable characterization of the zircon Ce anomaly appears merited and has now attracted the interest of many geologists. The Ce anomaly was generally estimated using La and Pr concentrations (e.g., Claiborne et al. 2010; Munoz et al. 2012; Gardiner et al. 2017). This method was recently revised by Ballard et al. (2002) and Loader et al. (2017). In the study of Loader et al. (2017), Nd and Sm replacing La and Pr were used in the calculation procedure for the zircon Ce anomaly. They stated that their revised method was more robust than that conventionally used as it did not need precise analyses of La and Pr contents for zircons, both of which are commonly below the limit of detection and susceptible to contamination from mineral inclusions. Alternatively, zircon Ce anomaly can be described using Ce⁴⁺/Ce³⁺. Trail et al. (2015) argued that Ce valence of zircon could be directly measured by X-ray Absorption Near Edge Structure (XANES). Nevertheless, direct determination of zircon Ce⁴⁺/Ce³⁺ by spectroscopic methods is now less common in literature because γ-irradiation experiments suggest that decrease of Ce⁴⁺ will occur due to the reductive reaction of Ce⁴⁺ to Ce³⁺ by the radiation effect of U and Th in the zircon (Takahashi et al. 2003). Ballard et al. (2002) proposed an indirect method to describe zircon Ce⁴⁺/Ce³⁺ ratios according to the lattice strain model. Compared to the Ce/Ce* method, Ballard et al.’ method may be less sensitive to contamination from light REE (LREE)-enriched inclusions since it no longer depends on just a pair of LREEs (La-Pr or Nd-Sm), but it needs REE compositions of not only zircons but also melts from which these zircons crystallized.

In this study, using published databases for zircon, we compare these different methods in characterizing zircon Ce anomalies for estimation of the oxidation state of magmas. It is shown that the above methods are yet neither satisfactory nor convincing, and in some cases may provide misleading information for oxidation state of magmas. By re-examining the chondrite-normalized REE patterns of zircon, a revised method to describe zircon Ce (and also Eu) anomalies is proposed. We demonstrate that the revised method in this study is valid and sensitive in examining magmatic oxidation state.

The Ce⁴⁺/Ce³⁺ method

The crystal lattice strain model

According to the method of Ballard et al. (2002), the total concentrations of Ce in zircon and co-existing melt can be considered as sums of Ce³⁺ and Ce⁴⁺ in each phase:

$$ {\mathrm{Ce}}_{\mathrm{zircon}}={\mathrm{Ce}}_{\mathrm{zircon}}^{3+}+{\mathrm{Ce}}_{\mathrm{zircon}}^{4+} $$

(1)

$$ {\mathrm{Ce}}_{\mathrm{melt}}={\mathrm{Ce}}_{\mathrm{melt}}^{3+}+{\mathrm{Ce}}_{\mathrm{melt}}^{4+} $$

(2)

By introducing partition coefficients (${\mathrm{D}}_{{\mathrm{Ce}}^{3+}}^{\mathrm{zircon}/\mathrm{melt}}$ and ${\mathrm{D}}_{{\mathrm{Ce}}^{4+}}^{\mathrm{zircon}/\mathrm{melt}}$) for Ce³⁺ and Ce⁴⁺, respectively, Eqs. 1 and 2 can be then substituted, combined and rearranged to produce the following expression:

$$ {\left(\frac{{\mathrm{Ce}}^{4+}}{{\mathrm{Ce}}^{3+}}\right)}_{\mathrm{zircon}}=\frac{{\mathrm{Ce}}_{\mathrm{melt}}-\frac{{\mathrm{Ce}}_{\mathrm{zircon}}}{{\mathrm{D}}_{{\mathrm{Ce}}^{3+}}^{\mathrm{zircon}/\mathrm{melt}}}}{\frac{{\mathrm{Ce}}_{\mathrm{zircon}}}{{\mathrm{D}}_{{\mathrm{Ce}}^{4+}}^{\mathrm{zircon}/\mathrm{melt}}}-{\mathrm{Ce}}_{\mathrm{melt}}} $$

(3)

where Ce_zircon can be measured directly by laser ablation–inductively coupled plasma–mass spectrometry (LA–ICP–MS) or secondary ion mass spectrometry (SIMS) methods; whole rock Ce is suggested as an equivalent of Ce_melt; and ${\mathrm{D}}_{{\mathrm{Ce}}^{3+}}^{\mathrm{zircon}/\mathrm{melt}}$ and ${\mathrm{D}}_{{\mathrm{Ce}}^{4+}}^{\mathrm{zircon}/\mathrm{melt}}$ can be calculated using the lattice strain model proposed by Blundy and Wood (1994).

According to the lattice strain model, for an isovalent series of ions with charge n + and radius r_i entering the crystal lattice site, the partition coefficient, D_i, can be described in terms of three parameters (Eq. 3): ${r}_0^{\mathrm{n}+}$, the optimum radius of that site; Eⁿ⁺, the Young’s Modulus of that site caused by ions that are bigger or smaller than ${r}_0^{\mathrm{n}+}$; and ${D}_0^{\mathrm{n}+}$, the “strain-compensated” partition coefficient for a (fictive) ion with radius ${r}_0^{\mathrm{n}+}$:

$$ \ln \kern.2em {\mathrm{D}}_{\mathrm{i}}=\ln \kern.2em {\mathrm{D}}_0^{\mathrm{n}+}-\frac{4\uppi {\mathrm{N}}_{\mathrm{A}}{\mathrm{E}}^{\mathrm{n}+}}{\mathrm{RT}}\left(\frac{{\mathrm{r}}_{\mathrm{i}}}{3}+\frac{{\mathrm{r}}_0^{\mathrm{n}+}}{6}\right){\left({\mathrm{r}}_{\mathrm{i}}-{\mathrm{r}}_0^{\mathrm{n}+}\right)}^2 $$

(4)

where N_A is Avogadro’s Number, R is the universal gas constant and T is in K.

When plotted according to Eq. 4, linear arrays can be found on a lnD_i versus $\left(\frac{{\mathrm{r}}_{\mathrm{i}}}{3}+\frac{{\mathrm{r}}_0^{\mathrm{n}+}}{6}\right){\left({\mathrm{r}}_{\mathrm{i}}-{\mathrm{r}}_0^{\mathrm{n}+}\right)}^2$ plot for an isovalent series of ions. By fitting these lines to trivalent REE and tetravalent Zr-Hf-U-Th arrays, values for ${\mathrm{D}}_{{\mathrm{Ce}}^{3+}}^{\mathrm{zircon}/\mathrm{melt}}$ and ${\mathrm{D}}_{{\mathrm{Ce}}^{4+}}^{\mathrm{zircon}/\mathrm{melt}}$ can be estimated. The resulting partition coefficients for Ce³⁺ and Ce⁴⁺ can be substituted into Eq. 3 and used in the calculation of the Ce⁴⁺/Ce³⁺ values. A similar procedure can be applied to Eu, which is also an important REE with variable valences (Eu²⁺ and Eu³⁺).

Assessment of assumptions for Ce⁴⁺/Ce³⁺ calculation

According to the procedure of Ballard et al. (2002), two assumptions have to be made to obtain Ce⁴⁺/Ce³⁺ values by the lattice strain model. Firstly, whole rock geochemistry is assumed to be identical to melts from which zircon crystalizes, so whole rock REE (as well as Hf, Th, U, Zr) data can be used as a proxy for melt chemistry when calculating D_i. Secondly, REEs and tetravalent cations (including Hf, Th and U) will only substitute the Zr site in zircon crystals, thus in previous studies it was assumed that both ${\mathrm{r}}_0^{3+}$ and ${\mathrm{r}}_0^{4+}$ were identical to the ionic radius of Zr in eight-fold coordination (0.84 Å).

Unfortunately, the first assumption has been demonstrated untenable. By studying zircon samples from the world class Oyu Tolgoi porphyry Cu-Au deposit, Loader et al. (2017) demonstrated that if titanite crystallization precedes or accompanies zircon crystallization, REE_melt would be lower than REE_WR, causing estimates of Ce⁴⁺/Ce³⁺ to be erroneously high. Crystallization of other accessory minerals, e.g., apatite and monazite, which are also common REE-enriched minerals in magmatic rocks, may play a similar role in affecting melt compositions from which zircon crystallized. Besides, zircons from many intrusions record a protracted and complex history of fractional crystallization and crustal assimilation, and thus even zircons from the same hand specimens may be characterized by distinct melt compositions (Claiborne et al. 2010). Therefore, using whole rock data as a proxy for melt compositions is inherently problematic under certain circumstances.

The second assumption is also problematic. According to Eq. 4, an inverse parabolic distribution can be expected in an ionic radius versus lnDi diagram (also called Onuma diagram) for isovalent ions entering the Zr site, and that the “best-fit” ion should have a radius ${\mathrm{r}}_0^{\mathrm{n}+}$ and a corresponding partition coefficient ${\mathrm{D}}_0^{\mathrm{n}+}$ at the peak of the parabola (Fig. 1). However, Blundy and Wood (1994) argued that values for the optimum site radius, ${\mathrm{r}}_0^{\mathrm{n}+}$, were not necessarily the same as the radius of the cation normally resident at the site of interest. Wood and Blundy (1997) further pointed out that ${\mathrm{r}}_0^{\mathrm{n}+}$ was actually sensitive to variations of crystal composition, as well as temperature and pressure. To find the real ${\mathrm{r}}_0^{3+}$ for zircon is not a trivial matter since changing ${r}_0^{3+}$ will result in variations of zircon Ce⁴⁺/Ce³⁺ values and determine to what degree rocks with different oxidation state of magmas can be distinguished using Ce⁴⁺/Ce³⁺ (Fig. 2). In the following section, we demonstrate that using the radius of Zr (0.84 Å) to replace ${\mathrm{r}}_0^{3+}$ in the procedure of Ballard et al. (2002) is not accurate.

Fig. 1

Cartoon illustrating the lattice strain model for isovalent ions entering into a zircon crystal (modified after Blundy and Wood 2003)

Fig. 2

Comparison of zircon Ce⁴⁺/Ce³⁺ values calculated by assuming ${\mathrm{r}}_0^{3+}$ = 0.84 and 0.953, respectively. ${\mathrm{r}}_0^{3+}$ = 0.953 was chosen for comparison according to the second section. The low coefficient of determination (R² = 0.1888) indicates that the change of ${\mathrm{r}}_0^{3+}$ will result in unpredictable variations of zircon Ce⁴⁺/Ce³⁺ values. Zircon data are from porphyry copper deposits of northern Chile (Ballard et al. 2002), with fertile and barren intrusions representing magmas of high and low oxidation state, respectively

Revised ${\mathrm{r}}_0^{3+}$ for zircon

In Online Resource 1, trace element data of a total of 35 experimental and natural zircon samples were considered to examine the optimum radius (${\mathrm{r}}_0^{\mathrm{n}+}$) for ions substituting the Zr site of zircon lattice. They were complied from publications of Burnham and Berry (2012); Colombini et al. (2011); Luo and Ayers (2009); Mahood and Hildreth (1983); Marshall et al. (2009); Murali et al. (1983); Sano et al. (2002) and Thomas et al. (2002). This database was chosen because not only zircon, but also melt compositions (by analyzing co-existing groundmass in volcanic rocks or melt inclusions) were measured directly, thus giving exact estimation for partition coefficients of trace elements. In this table, ${\mathrm{r}}_0^{3+}$ values were calculated by fitting the partition coefficients of trivalent REE (and Sc where possible) ions to a parabolic curve in the Onuma diagram.

Variations of ${\mathrm{r}}_0^{3+}$ are observed between samples and between different zircon crystals from the same samples. Nevertheless, the majority of ${\mathrm{r}}_0^{3+}$ values for zircon crystals range from 0.914 to 1.010 Å, which are significantly greater than the radius for Zr (0.84 Å) (Fig. 3). Only data from Thomas et al. (2002) show ${\mathrm{r}}_0^{3+}$ values smaller than 0.84 Å, which however may not be convincing because only a few REEs were reported in their study and thus cannot constrain the parabola well. Figures 1 and 3 manifests that if there is no element with ${\mathrm{r}}_{\mathrm{i}}^{3+}$<${\mathrm{r}}_0^{3+}$, the predicted ${\mathrm{r}}_0^{3+}$ by the lattice strain model will also be inaccurate (Blundy and Wood 2003). In Online Resource 1, only five zircon grains from Colombini et al. (2011) and Burnham and Berry (2012) reported D_Sc (r_Sc=0.87 Å) and the full REEs simultaneously, giving a narrow range from 0.948 to 0.959 Å with an average of 0.953 Å. This value is currently the best estimation for zircon ${\mathrm{r}}_0^{3+}$. A recent study for zircon by Smythe and Brenan (2016) also found the Zr radius could not represent the true ${\mathrm{r}}_0^{3+}$. They used ${\mathrm{r}}_0^{3+}$ = 0.93 Å in their zircon Ce⁴⁺/Ce³⁺ calculation, which is almost identical to the result from this study. All these indicate that “real”${\mathrm{r}}_0^{3+}$ for zircon should be at least greater than 0.84 Å.

Fig. 3

Onuma diagrams showing the parabola was obtained by fitting the reported zircon/melt partition coefficients of REEs (plus Y and Sc) to the lattice strain model of Blundy and Wood (1994). Information in legends corresponds to data sources and sample names (see Online Resource 1 for more details). The asterisk indicates data were from experimental samples, whereas others from natural samples. Ce and Eu were not included because they are susceptible to oxidation condition and/or feldspar fractionation, respectively. The inserted figure is a cartoon illustrating the lattice strain model of trace element partitioning for a crystal

The ${\mathrm{r}}_0^{4+}$ values, however, cannot be estimated in the same way as for ${\mathrm{r}}_0^{3+}$, because there are too few elements (normally Hf, Th, U and Zr) to constrain the parabolas for tetravalent ions. Nonetheless, unlike ${\mathrm{r}}_0^{3+}$, the assumption that ${\mathrm{r}}_0^{4+}$ = 0.84 Å may be reasonable. Since r_Hf < r_Zr > r_{Th, U} and for most terrestrial zircons D_Zr > D_Hf,Th,U, it is not hard to see from Fig. 1 that the “best-fit” tetravalent ion entering into zircon may have a radius similar to Zr.

It should be noted that D_REE was found to be related to temperature (e.g., Rubatto and Hermann 2007), and therefore the calculated ${\mathrm{r}}_0^{3+}$ values, according to the lattice strain model, may vary according to the melt temperatures. However, a correlation between the calculated ${\mathrm{r}}_0^{3+}$ values and melt temperatures is not observed (see Online Resource 1), ruling out that the observed high ${\mathrm{r}}_0^{3+}$ values reflect discrepancies due to melt temperatures. Nevertheless, considering the variations of ${\mathrm{r}}_0^{3+}$ observed in Online Resource 1, we recommend using ${\mathrm{r}}_0^{3+}$ = 0.954 Å in calculation of zircon Ce⁴⁺/Ce³⁺ only after it is further verified by more experimental studies, or we suggest that ${\mathrm{r}}_0^{3+}$ must be determined by regression in each case. In either case, using the method of Ballard et al. (2002) to calculate zircon Ce⁴⁺/Ce³⁺ may be problematic in practice, and caution must be taken when using it to evaluate oxidation state of magmas.

Alternative method to calculate Ce⁴⁺/Ce³⁺

Rearrangement of Eq. 3 can give the following equation:

$$ \ln \kern.2em {\mathrm{D}}_{\mathrm{i}}=\frac{\mathrm{A}}{6}\left[-2{\mathrm{r}}_{\mathrm{i}}^3+3{\mathrm{r}}_0^{\mathrm{n}+}{\mathrm{r}}_{\mathrm{i}}^2+\left(\frac{6}{\mathrm{A}}\ln \kern.2em {\mathrm{D}}_0^{\mathrm{n}+}-{\left({\mathrm{r}}_0^{\mathrm{n}+}\right)}^3\right)\right] $$

(5)

where $\mathrm{A}=\frac{4\uppi {\mathrm{N}}_{\mathrm{A}}{\mathrm{E}}^{\mathrm{n}+}}{\mathrm{RT}}$; and $\ln {\mathrm{D}}_0^{\mathrm{n}+}$ and ${\mathrm{r}}_0^{\mathrm{n}+}$ are both constants for a given zircon crystal.

Thus, Eq. 5 can be regarded as a cubic function in the form f(x) = ax³ + cx² + d, where f(x) = lnD_i and x = r_i. When lnD_i and r_i are plotted according to Eq. 5, a parabola can be obtained in the Onuma diagram, by which ${\mathrm{D}}_{{\mathrm{Ce}}^{3+}}^{\mathrm{zircon}/\mathrm{melt}}$ can be estimated directly without assigning a value to ${\mathrm{r}}_0^{3+}$. The same method was used by Burnham and Berry (2017). ${\mathrm{D}}_{{\mathrm{Ce}}^{4+}}^{\mathrm{zircon}/\mathrm{melt}}$ can be estimated in a similar way. By substituting the partition coefficients for Ce⁴⁺ and Ce³⁺ into Eq. 3, the Ce⁴⁺/Ce³⁺ ratio can be calculated.

However, this method should still be used with caution, since the exact melt composition from which zircon crystalized is still needed in the calculation process. Also, Sc contents of zircon and melts are needed to constrain the ${\mathrm{r}}_0^{3+}$ (and thus the parabola for trivalent ions), but unfortunately most studies of zircon do not report Sc (Burnham and Berry 2012). This is partly due to the analysis of Sc concentrations in zircon being complicated because, in addition to ²⁹Si¹⁶O, ⁹⁰Zr²⁺ interferes with Sc at mass 45 (Wopenka et al. 1996).

The Ce/Ce* method

The existing methods

The most classic method to quantify Ce anomaly is Ce/Ce*, where Ce* = $\sqrt{\mathrm{L}{\mathrm{a}}_{\mathrm{N}}\times \mathrm{P}{\mathrm{r}}_{\mathrm{N}}}$ (the subscript indicates chondrite normalization). However, La and Pr are very low in abundance and frequently below the limit of detection. Therefore, this method often incurs criticism where used, since it is susceptible to contamination by small LREE-enriched melt/mineral inclusions that are common in zircon (Bindeman et al. 2014). Loader et al. (2017) modified the traditional method and proposed an alternative formulation for Ce*, based on a similar geometric relationship previously used:

$$ {\mathrm{N}\mathrm{d}}_{\mathrm{N}}=\sqrt{\mathrm{C}{\mathrm{e}}^{\ast}\times \mathrm{S}{\mathrm{m}}_{\mathrm{N}}} $$

(6)

which, when rearranged, yields an expression for Ce*,:

$$ \mathrm{C}{\mathrm{e}}^{\ast }=\frac{\mathrm{N}{{\mathrm{d}}_{\mathrm{N}}}^2}{\mathrm{S}{\mathrm{m}}_{\mathrm{N}}} $$

(7)

This formulation does not require the accurate determination of either La or Pr, thus it is believed to be more robust than the conventional method. However, many studies have argued that in chondrite-normalized REE patterns, zircons REEs show a concave-downwards shape, rather than a linear shape. Thus, both $\sqrt{\mathrm{L}{\mathrm{a}}_{\mathrm{N}}\times \mathrm{P}{\mathrm{r}}_{\mathrm{N}}}$ and $\frac{\mathrm{N}{{\mathrm{d}}_{\mathrm{N}}}^2}{\mathrm{S}{\mathrm{m}}_{\mathrm{N}}}$ are not the best estimation for Ce*, and will result in underestimation (the conventional method) or overestimation (Loader et al.’s (2017) method) of the true Ce*. Besides, Loader et al. (2017) showed that Ce/Ce* calculated by their method showed a great overlap between fertile and barren rocks, indicating that it is not a robust tool for estimation of the oxidation state of magmas.

Revised Ce/Ce* method by this study

To determine Ce*, the key is to find a function which can perfectly match the position of LREEs and middle REEs (MREEs) with 3+ cations in the chondrite-normalized REE pattern for zircon. We therefore examined ca. 1000 zircon trace element data from intrusions in 11 giant to supergiant porphyry Cu(-Au) deposits worldwide. Complied descriptions about these deposits and literature sources were listed in Online Resource 2. In this study, La to Lu (15 elements) were assigned numbers 1 to 15 (Pm was 5 although it does not exist in nature) in the chondrite-normalized REE pattern. Regression calculation clearly shows that REE patterns of zircons from Online Resource 2 can be matched by the logarithmic function (y = alog(x) + b, where y is the logarithm (base 10) of chondrite-normalized REEs) (Fig. 4), with more than 99% of 847 fitting curves characterized by the coefficient of determination (R²) of >0.95 (Online Resource 2). It should be noted that during the regression calculation La and Pr were not included due to their uncertain concentrations. With these fitting curves, Ce* as well as Eu* (previously constrained by $\sqrt{\mathrm{S}{\mathrm{m}}_{\mathrm{N}}\times \mathrm{G}{\mathrm{d}}_{\mathrm{N}}}$) can be precisely determined.

Fig. 4

Chondrite-normalized REE patterns for zircon from the Bozshakol porphyry Cu deposit, Kazakhstan (Shen et al. 2015). All these zircon REE compositions can be well described by logarithmic function curves with high coefficient of determination (R² > 0.95). Note that La, Ce, Pr and Eu were not involved during the regression calculation. The fitting curves seem to support that the La concentration of natural zircon should be far below that of chondrite (0.237 ppm). CI chondrite data are from McDonough and Sun (1995)

Using this revised method, we re-calculated zircon Ce/Ce* and Eu/Eu* values from 11 porphyry Cu(-Au) deposits (see Online Resource 2 for more details). Compared to the existing method, both Ce/Ce* and Eu/Eu* values calculated by the revised method show a much narrower range (Fig. 5). Moreover, although zircon Ce/Ce* and Eu/Eu* values calculated by either method exhibit some overlaps between barren and fertile intrusions, those calculated by the revised method behave more regularly and can more clearly distinguish fertile from barren intrusions (Fig. 5). Since the high oxygen fugacities of the fertile intrusions in these porphyry deposits have been demonstrated previously, all these data manifest the advantages of the revised method in determining magmatic oxidation state and distinguishing between fertile and barren intrusions for porphyry Cu systems.

Fig. 5

Comparison of different methods in calculating zircon Ce and Eu anomalies. Zircon data were from 11 giant to supergiant porphyry Cu (-Au) deposits, where fertile and barren intrusions had been demonstrated related to magmas of high and low oxidation state, respectively. (a) Ce/Ce* and Eu/Eu* were calculated by the conventional method following the formulas: Ce* =$\sqrt{\mathrm{L}{\mathrm{a}}_{\mathrm{N}}\times \mathrm{P}{\mathrm{r}}_{\mathrm{N}}}$ and Eu* =$\sqrt{\mathrm{S}{\mathrm{m}}_{\mathrm{N}}\times \mathrm{G}{\mathrm{d}}_{\mathrm{N}}}$. (b) Eu* was calculated using the conventional method, whereas Ce* follows the procedure proposed by Loader et al. (2017). (c) Ce* and Eu* were calculated by the revised method in this study by fitting zircon REE data to logarithmic function curves. The insert represents the full results which show a large range of Eu/Eu*

Implication for La entering into zircon lattice

Figure 4 shows that the predicted La values (<0.001 ppm) by the logarithmic function curve for zircon from the Bozshakol porphyry Cu deposit are far below the measured values (0.001–0.20 ppm) (Online Resource 2; Shen et al. 2015). In Online Resource 2, predicted zircon La values for ten other deposits display a similar behavior and are all far below 0.03 ppm (~0.1 times that for chondrite). A similar case is found for Pr.

The following scenario may account for the deviation of the measured La (as well as Pr) concentrations from the fitting lines: (1) inapplicability of the logarithmic function in describing the zircon REE behaviors; (2) contamination from LREE-enriched minerals in the measured La values; and (3) modification by hydrothermal alteration of the measured zircon. The fitting logarithmic function lines obtained in our study are characterized by very high R² (>0.95), which indicates that they can describe the behaviors of REEs. Nevertheless, it is quite possible that the logarithmic function line found in our study is not applicable to hydrothermal zircon, considering that zircon-fluid partition coefficients for REEs may be quite different from zircon-melt partition coefficients. Theoretically, it is possible that hydrothermal alteration can explain the abnormally elevated La and Pr contents measured in some zircons if zircon-fluid partition coefficients for La and Pr are higher than their zircon-melt partition coefficients. However, direct experimental data are absent to demonstrate that zircon crystallized from fluid can have higher La and Pr partition coefficients. Besides, based on petrographic observation, trace element compositions and ages of those zircons, hydrothermal alteration on these studied zircons had been excluded (Ballard et al. 2002; Loader et al. 2017; Munoz et al. 2012; Shen et al. 2015; Zhu et al. 2018).

Therefore, the most likely explanation for the high La and Pr content in some natural zircons versus the modelled values is contamination by mineral, fluid and/or melt inclusions. Zhong et al. (2018) showed at least one case where abundant mineral inclusions are present in zircon, some of which were supposedly sub-micrometer in size and hard to monitor during the zircon trace element analysis. They demonstrated that only 0.05 vol.% contamination from those small LREE-enriched mineral inclusions (e.g., apatite and titanite) would result in evident elevation of La and Pr for zircon. Thus, it is quite possible that the deviation of predicted La contents by the logarithmic function from measured La contents just results from mineral inclusion contamination in zircon during in situ microanalysis (Fig. 4 and Online Resource 2). If this conclusion is correct and La = 0.03 ppm is taken as a critical value for clean zircon, approximately 40% of zircons in Online Resource 2 suffered contamination from LREE-enriched minerals during the analysis process. We thus suggest that La entering into zircon is much lower than what commonly measured, and that accidental sampling of melt/mineral inclusions during the zircon trace element analysis is very common due to their widespread occurrence and small sizes.

It should be noted that a recent study by Zou et al. (2019) also proposed a cut-off La value (< 0.1 ppm) for inclusion-free zircon based on a compilation of zircon data from the GEOROC database (which shows that most zircons are characterized by La contents as low as 0.1 ppm). According to our study, La entering into zircon should be much lower than this value and most zircon trace element data (at least for LREEs) in the GEOROC database appear to be contaminated by LREE-enriched melt/mineral inclusions.

Advantages and applications

The revised Ce/Ce* method introduced by our study shows some advantage in calculation of the Ce anomaly and evaluation of the oxidation state of magmas. Compared to the Ce⁴⁺/Ce³⁺ methods, the revised Ce/Ce* method does not involve melt composition which is generally inaccessible for intrusive rocks and detrital zircons. Also, our method only requires REE concentrations in zircon when compared to Ballard et al. (2002) that requires some less constrained parameters (e.g., ${\mathrm{r}}_0^{\mathrm{n}+}$), which can result in high errors of the calculated Ce⁴⁺/Ce³⁺. Furthermore, using published databases for zircon, our study shows that the new method is more valid in discriminating the oxidation state of different magmas.

However, the presented method still should be used carefully when evaluating the magmatic oxidation state. Points to keep in mind are other possible controls on the Ce anomaly apart from fO₂ of magmas and the behaviors of LREE entering into zircon. Current studies have illustrated that the Ce anomaly is not just a function of fO₂, but also of the water content of the melt, temperature and melt composition (Smythe and Brenan 2015). This may explain why in Fig. 5c there is still an overlapping area between barren and fertile samples. Nevertheless, the water content of melt, temperature and melt composition seem to play a less predominant role on the Ce anomaly compared to fO₂. Figure 6 shows the variations of Ce⁴⁺/Ce³⁺ of Ce doped glasses with fO₂ from Smythe and Brenan (2015). These glasses varied in composition (basalt to rhyolite) and were synthesized under quite different temperatures (1200 to 1500 °C), water contents (<2.8 to 9.1 wt%) and oxidation states (FMQ −4.0 to FMQ +8.4). It shows that although the water content, temperature and melt composition are quite different, the Ce anomaly and fO₂ of synthesized glasses still show a strong correlation. Same information can also be gained by examining the coefficients of Eq. 15 in Smythe and Brenan (2015). This means that despite probably being defective, using the Ce anomaly to reflect the degree of fO₂ is still feasible in practice.

Fig. 6

log(fO₂) vs. log(Ce⁴⁺/Ce³⁺) for a series of synthesized glasses ranging in composition from basalt to rhyolite, under fO₂ conditions varying from FMQ −4.0 to FMQ +8.4, water contents from <2.8 to 9.1 wt% and temperatures from 1200 to 1500 °C. Data are from Smythe and Brenan (2016)

The lack of the exact knowledge on the behaviors of LREE entering into zircon is a limiting factor for the use of the method (and also for all other existing methods using the zircon Ce anomaly to evaluate fO₂). Fundamental to our new method is the assumption of the consistent mechanism of LREE with MREE and heavy REE (HREE) entering from melt into zircon. This is hardly to be verified since La and Pr contents in measured zircons are commonly below the limits of detection of most analytical technologies. Alternatively, using the lattice-strain model of Blundy and Wood (1994), the zircon-melt partition coefficients for the LREE can be predicted, then the LREE values in the zircon can be calculated theoretically if melt compositions are known, which then can be compared to the values determined by our method. However, the exact composition of melts from which zircon crystallized is nearly inaccessible as discussed in the second section, even for experimental and/or natural phenocryst (zircon)-glass pairs.

Overall, although on-going works are still needed to demonstrate its effectiveness, the presented method provides a simple way to predict the magmatic oxidation state. In practice, the revised Ce/Ce* method is particularly important for analyses where the provenance of the analyzed zircon is unknown or in question (e.g., detrital zircons). It may also be important for analyses of rock samples from hydrothermal mineral deposits, which are generally strongly modified by hydrothermal alteration. In this case, zircons collected from rocks may be still of magmatic origin but whole rock geochemistry cannot represent the real melt from which the rock formed, let alone represent the melt where zircon crystallized. Moreover, for zircon studies where REEs are not fully provided, we can also use the logarithmic curve on chondrite-normalized diagrams to give the missing REE concentrations, especially for the high atomic number REEs.

Conclusions

Although feasible in theory, both the existing Ce/Ce* methods and the Ce⁴⁺/Ce³⁺ method proposed by Ballard et al. (2002) are problematic in practice and thus may not be a good proxy for the oxidation state of magmas. The correct ${\mathrm{r}}_0^{3+}$ for zircon lattice strain modelling is not 0.84 Å and may approach 0.953 Å. A revised Ce/Ce* method is proposed in our study, which is based on the new finding that zircon REEs follow a logarithmic function in chondrite-normalized patterns. Application of this new method to published databases for zircons shows that it is suitable to estimate oxygen fugacity and more effective in discriminating fertile from barren rocks in porphyry Cu systems. The study also manifests that La contents in magmatic zircon should be generally very low (perhaps <0.03 ppm) and many reported magmatic zircons with high La contents were probably just contaminated by small melt/mineral inclusions during analyses.