Mining and exploration companies routinely use four-acid digestion (4AD), inductively coupled plasma, atomic emission spectra/mass spectrometry methods from commercial assay laboratories for analysing drill and rock samples for lithogeochemical assessment and resource reporting. This method is also known to exhibit lower recovery of elements hosted by resistate minerals. To assess the impact of lower recoveries on lithogeochemical interpretation, a suite of commonly used elements for lithogeochemical analysis (high-field-strength elements Zr, Hf, Nb, Ta, Ti and Eu and transition elements V and Sc) was analysed by 4AD and alkali fusion/acid digestion (AFAD). Lower recoveries in the 4AD relative to the AFAD were recorded for Zr, Hf, Nb, Ta, Ti and Eu; Sc and V reported similar concentrations for both decomposition methods. Despite the lower recoveries for Nb, Ta and Ti, element ratios were largely preserved with the 4AD method due to the recoveries covarying at a 1:1 ratio. A plot of Ti/Nb against V/Sc was found to be largely unaffected by decomposition method, producing similar compositional classifications between the two digestion methods. Use of the Eu anomaly (Eu/Eu*) to determine plagioclase fractionation was also found to be unaffected by decomposition method. In contrast, a standard Zr/Ti v. Nb/Y discrimination plot produced incorrect classifications with 4AD producing more mafic and alkaline classifications relative to the AFAD method. Magmatic fertility interpretations utilizing Zr/Hf were also found to be affected in the 4AD results due to the lower recovery of Zr relative to Hf. This resulted in a bias in the 4AD results and produced false-positive anomalism in fertility assessments. Multiple decomposition methods including combinations of acid and fusion methods are recommended for lithogeochemical analysis utilizing large regions of the periodic table. However, if only 4AD data are available, plots such as Ti/Nb v. V/Sc and Nb/Ta, which preserve their ratios, can be used for lithogeochemical classification.

Supplementary material: Wholerock geochemical data and detailed methods are available at https://doi.org/10.6084/m9.figshare.c.6444444

Fractionation trends recorded by igneous rocks during magma evolution can be characterized by assessing changes in ratios of select element pairs (Winchester and Floyd 1977; Deering and Bachmann 2010; Halley 2020). Petrogenic ratios such as Nb/Ta and Zr/Hf can be sensitive to mineralogical changes during magma fractionation that reflect changes to melt composition, temperature, pressure and oxygen fugacity (Linnen and Keppler 2002; Claiborne et al. 2006; Fulmer et al. 2010; Mallmann et al. 2014; Holycross and Cottrell 2020). Various authors have utilized different element pair systematics to infer mineral compositions of melts such as; V–Sc to infer titanomagnetite and magnetite formation (Aeolus Lee et al. 2005; Loucks 2014; Wijbrans et al. 2015), Sr–Y–Eu systematics as a proxy for plagioclase and hornblende crystallization (Drake and Weill 1975; Defant and Drummond 1990; Deering and Bachmann 2010; Richards et al. 2012), Nb–Ta for Ti-oxides and biotite crystallization (Linnen and Keppler 1997; Klemme et al. 2006; John et al. 2011; Stepanov and Hermann 2013) and Zr–Hf for zircon, pyroxene and garnet crystallization (David et al. 2000; Linnen and Keppler 2002; Weyer et al. 2003; Claiborne et al. 2006; Zaraisky et al. 2008).

Recent studies of magma fertility and its links to magmatic–hydrothermal mineralization have focused on trace element ratios in both mineral and whole-rock analysis (Loader et al. 2017; Cheng et al. 2018; Liu and Chen 2018; Wells et al. 2020; Mohammadi et al. 2021). Some of these results are relevant to mineral exploration since they provide a rapid and affordable method of selecting exploration areas using whole-rock geochemistry and can be applied to all scales from regional selection, such as identifying fertile magmatic terranes, to deposit-scale exploration and targeting. Large volumes of whole-rock chemical data are commonly obtained through analytical packages offered at commercial laboratories. The four-acid digestion, inductively coupled plasma mass spectrometry (4AD-ICP-MS) method is frequently used by mineral explorers for its low cost, rapid throughput and low detection limits. This has helped it become the dominant method for trace-element analysis in mineral exploration and resource reporting.

4AD uses combinations of nitric acid (HNO3), perchloric acid (HClO4), hydrofluoric acid (HF) and hydrochloric acid (HCl) to decompose the sample. These digestions are often referred to as ‘near-total’ digests due to their ability to decompose most mineral matrices, including Si-bonds (Chao and Sanzolone 1992). Acid digests produce significantly lower matrix interferences, resulting in improved sensitivity; however, some of the methods can also suffer from partial recovery of resistate minerals such as zircon and rutile (Hoops 1964; Langmyhr and Sveen 1965; Hall and Plant 1992; Totland et al. 1992; Pashkova et al. 2019). This is especially problematic for high-field-strength elements (HFSE) such as Ta (Münker 1998) and sulfates of Ba, Sr and Pb (Hoops 1964). The use of HF also produces fluoride complexes of Si, As, B, Ti, Nb, Ta, Ge and Sb which are lost in open vessel digestion, but also result in improving detection limits by lowering the plasma load during ICP-MS analysis (Chao and Sanzolone 1992). The lower recovery of HFS elements in 4AD is therefore problematic for lithogeochemical analysis which makes extensive use of these elements in interpretation.

Fusion using a lithium–borate flux is often used as an alternative to direct acid digests as it is generally considered to provide a more complete decomposition and produce more accurate results for HFSE when compared to acid digestion (Yu et al. 2001; Potts et al. 2015). However, the higher blank levels and higher salt contents associated with the flux material have the effect of lowering sensitivity and increasing detection limits (Longerich et al. 1990; Totland et al. 1992; Yu et al. 2001; Panteeva et al. 2003). This is particularly important when analysing trace elements with low abundances. Additionally, volatile elements such as Se, Cs, Rb, Pb, Sn, Sb, Tl, some PGE and Zn can potentially be lost due to the high fusion temperatures (Jarvis 1990; Totland et al. 1992; Yu et al. 2001; Senda et al. 2014). Fusion is also known to result in partial decomposition of the sample when highly resistant minerals such as zircon and certain metal oxides, REE phosphates and sulfides are present (Cremer and Schlocker 1976; Potts et al. 1990). Additional treatments are recommended when significant amounts of these minerals are present (Chao and Sanzolone 1992).

Other authors have experimented with combining the fusion and multi-acid digestion methods. Alkali fusion followed by acid digest (AFAD) attempts to combine the benefits of both techniques by producing a total digest without contaminants that may be present in the flux and matrix effects (Senda et al. 2014). These authors found that compared to a standard AD, the AFAD method produced more accurate results, especially for HFSE and HREE where resistate minerals such as zircon were more efficiently digested. However, there were significant losses in volatile elements as discussed above.

In a study of the analysis of HFSE, REE and the impact of digestion type on tectonic interpretation, Hall and Plant (1992) compared REE analysis of LiBO2-ICP-MS with 4AD-ICP-MS. The authors found significant deviations in the MREE–HREE between the two digests. The mixed AD was found to be inefficient at dissolving HREE-bearing resistate minerals such as zircon and garnet, resulting in under-reporting of the elements Gd to Lu. Sholkovitz (1990), in a study using geochemical standards of marine sediments and shales, found that only 0.1 to 0.3% of insoluble zircon (by sample weight) can contribute to between 20 and 100% of the total whole-rock HREE budget. This presents a problem for mineral explorers who might rely on acid digestion data to screen for REE fractionation as a measure of metal fertility (Jamali 2017; Loader et al. 2017; Liu and Chen 2018) and highlights the importance of understanding the analytical requirements prior to method selection as no single method can provide high-quality results across the periodic table. For a detailed summary of the advantages and disadvantages of fusion and ADs, see Totland et al. (1992) and Yu et al. (2001).

We present a case study of Zr–Hf, Nb–Ta-(Ti), V–Sc and Eu–Eu* systematics from the Cambrian Mt Read Volcanics (MRV) metallogenic province of western Tasmania, Australia. The MRV is a 200 × 20 km belt of Middle Cambrian submarine volcanics and associated volcaniclastics with abundant sedimentary and intrusive rocks (Corbett 1992; Fig. 1). The belt formed during a post-collisional phase of extension following the c. 510 Ma Tyennan Orogeny (Berry et al. 2007). It is comprised of mostly calc-alkaline, rhyolitic to basaltic volcanics and volcaniclastics with minor tholeiitic volcanism (Crawford et al. 1992). The MRV is host to several volcanic-associated ore-deposits which have been mined continuously since the late nineteenth century. In this time c. 8.15 Mt Zn, 3.1 Mt Pb, 3.31 Mt Cu, 9.04 kt Ag and 279 t Au is thought to have been extracted, generating $72.4 billion dollars in estimated revenue (Blewett 2012). The MRV samples in the dataset presented here include a range of compositions from basaltic to rhyolitic and include both intrusive and extrusive equivalents and coherent and volcaniclastic facies with varying degrees and types of alteration.

Previous studies of the MRV noted difficulties in interpreting certain geochemical populations with low trace-element abundances not resolvable with the available X-ray fluorescence (XRF) analysis (Crawford et al. 1992). Other studies noted differences between XRF and 4AD-ICP-MS analyses in MRV samples and how these differences translated into variable petrogenetic indicators such as Ti/Zr, making lithogeochemical interpretation difficult (Wu 2014). To understand the fine-scale stratigraphy of the MRV and aid with future mineral exploration, a high-resolution geochemical dataset is required.

This paper focuses on a comparison between fusion methods and direct multi-ADs of whole-rock samples from the MRV and how the different techniques can introduce bias into the subsequent lithogeochemical and magma fertility interpretations. This is expected to have ramifications for exploration geologists who frequently employ AD methods in combination with legacy geochemical datasets.

V–Sc systematics

V–Sc systematics are often used to infer the fO2 evolution of magmatic systems (Anser Li and Aeolus Lee 2004; Aeolus Lee et al. 2005; Mallmann and O'Neill 2009; Wu et al. 2019). Anser Li and Aeolus Lee (2004) and Aeolus Lee et al. (2005) show that V/Sc is robust to mantle modification, including metasomatism, by demonstrating the constancy of V/Sc ratios between Archean basalts and modern basalts and between modern mid-ocean ridge basalt (MORB) and arc environments. In reduced, mafic melts, the V/Sc values are typically c. 7 and are controlled by clinopyroxene and amphibole crystallization where both V and Sc are present in their +3-valence state. In this state, V and Sc have similar partition coefficients and are therefore similarly depleted from the melt during crystallization (Doe 1997). In more oxidized melts, V occurs in V4+ and V5+ oxidation states and will strongly partition into magnetite while Sc remains in its +3-valence state and is much less compatible in magnetite (Wijbrans et al. 2015; Arató and Audétat 2017; Iveson et al. 2018). Additionally, V-Sc systematics can be used to infer the depth of crystal formation in certain settings. High V/Sc ratios in subduction-related magmatic environments can be attributed to accumulations of hornblende due to high-pressure, hydrous conditions. Based on both field observation and experimental data, hornblende crystallization is promoted to the start of the crystallization sequence with magnetite forming late at lower temperature and/or depth (Richards 2011; Loucks 2014). A decreasing V/Sc in an evolving melt can indicate the fractional crystallization of magnetite under oxidized conditions while an increasing V/Sc can indicate accumulations of hornblende under high-pressure, hydrous conditions.

Eu–Eu* Systematics

Unlike the other REE, Eu has long been known to be stable in its bivalent form (Goldschmidt 1937). This ‘anomalous’ behaviour was determined to be a predictable function of temperature and fO2 and the magnitude of the anomaly is related to the abundance of plagioclase where Eu2+ can substitute for Ca2+ (Weill and Drake 1973). The Eu anomaly can be defined as Eu2+/Eu3+ (Weill and Drake 1973); however, it is more commonly defined as Eu/Eu* where Eu* is an interpolant between Sm and Gd. The Eu anomaly is quantified as the enrichment or depletion of measured Eu, relative to the interpolated Eu value (Eu*), normalized to a baseline composition such as C1 carbonaceous chondrites (Taylor and McClennan 1985; Sun and McDonough 1989). In reduced magmas, minerals with strong affinities to Eu2+, such as plagioclase, will exhibit large positive Eu anomalies. Conversely, melts that have evolved from the fractional crystallization of plagioclase will display increasingly negative Eu anomalism with increasing melt fractionation (Deering et al. 2016). In oxidized and/or H2O-rich magmas, the extent of the negative Eu anomaly is expected to be much smaller, either due to Eu being prevalently in the Eu3+, or suppression of plagioclase crystallization, or both (Danyushevsky 2001; Loucks 2014; Tang et al. 2020).

Nb–Ta– (Ti) Systematics

Nb and Ta fractionation occurs in a variety of minerals and was noted to always be associated with Ti-rich phases (Schmidt et al. 2004). John et al. (2011) found that in melts derived from rutile and titanite-bearing eclogites, the Nb–Ta ratio was a function of the modal abundance of rutile and titanite, which itself is a function of pressure. When titanite is much more abundant than rutile, Nb/Ta was very high (>60) and when titanite was less than rutile, Nb/Ta was lower (≤30). For very low ratios (<16), a very high degree of partial melting was required, consuming all Ti-phases. Stepanov and Hermann (2013) compiled Nb–Ta partition coefficients and noted that the Ti-rich phases such as rutile, ilmenite and titanite preferentially incorporate Ta over Nb resulting in a low Nb/Ta restite with corresponding high Nb/Ta partial melt during crustal melting. In contrast, micas such as biotite and high-Ti phengite were found to prefer Nb over Ta, producing a high Nb/Ta restite with corresponding low Nb/Ta partial melt. Amphibole was also noted to preference Nb; however, the overall effect on Nb/Ta is small given the very low compatibility of Nb in amphibole. Low Nb/Ta values with decreasing Ti are therefore inferred to be the signature of the fractionation of mica (Stepanov et al. 2014; Ballouard et al. 2016; Sun et al. 2019; Halley 2020).

Zr–Hf Systematics

Numerous studies have observed wide variations in Zr/Hf across a range of compositions and magmatic settings. These include; MORB and ocean island basalt (Pfänder et al. 2007), abyssal peridotites and harzburgites (Weyer et al. 2003; Niu 2004), carbonatites (Dupuy et al. 1992; Andrade et al. 2002), chondrites and lunar material (Ahrens 1962; Ehmann et al. 1975), alkali granites and pegmatites (Van Lichtervelde et al. 2009; Fulmer et al. 2010; Kynicky et al. 2011) and peraluminous granites and pegmatites (Dostal and Chatterjee 2000; Claiborne et al. 2006; Zaraisky et al. 2008; Van Lichtervelde et al. 2009).

Studies by Bea et al. (2006) and Claiborne et al. (2006) found that zircons formed in late-stage, lower-temperature, high-silica granitic melts are Hf-rich, with the highest concentrations of Hf found in the rims of zircon that formed at the lowest temperature according to Ti-in-zircon geothermometry (Watson et al. 2006). This was interpreted to reflect the fractional crystallization of early, Zr-rich zircons which, after melt segregation, left behind a residual partial melt with a higher proportion of Hf. Subsequent crystallization of lower-temperature zircon enriched in Hf produce lower zircon Zr/Hf values. This effect is exaggerated by water and other complexing ligands such as fluorine which depolymerize the melt to allow for more efficient segregation of early zircon and to lower the temperature at which melt can exist, promoting Hf enrichment in the zircon and prolonging the zircon formation window (Keppler 1993; Linnen 1998; Linnen and Keppler 2002; Van Lichtervelde et al. 2009; Aseri et al. 2015; Shao et al. 2019).

Melts exhibiting decreasing Zr/Hf with increasing silica are interpreted to be due to the fractional crystallization of zircon (Linnen and Keppler 2002; Bea et al. 2006; Claiborne et al. 2006; Zaraisky et al. 2008; Deering and Bachmann 2010; Wang et al. 2010; Breiter et al. 2014; Wu et al. 2017; Yan et al. 2018).

A total of 654 samples from the MRV have been included in the study (Table 1). The samples included a combination of rock chips and drill core from the Mineral Resources Tasmania (MRT) core library and the University of Tasmania, Centre for Ore Deposit and Earth Sciences (CODES) rock collection. Samples were chosen to cover as much of the stratigraphy and the geographic extent of the province as possible and included all rock types (coherent and volcaniclastic facies) and varying degrees of alteration. All samples were analysed at Australian Laboratory Services (ALS), Perth, Western Australia by 4AD-ICP-MS (method name – ME-MS61L™ and MS61L-REE™) while 410 samples were also analysed at ALS (Perth) by AFADICP-MS (method name – ME-MS81s™). Ten samples were analysed at CODES Analytical Laboratories (CAL), at the University of Tasmania using a two-step multi-AD (Yu et al. 2001; Table 1).

At ALS, samples were weighed, crushed and pulverized to >85% passing 75 microns. For the samples that were analysed by 4AD ICP-MS, a 0.25 g aliquot of pulverized material was dissolved in a solution of nitric, perchloric and hydrofluoric acid at 185°C. The residual solution was then leached and diluted in a hydrochloric acid solution. The final solution was then analysed by ICP-MS and ICP-atomic emission spectroscopy.

For the AFAD method, a 0.1 g aliquot was added to a lithium metaborate/lithium tetraborate flux (LiBO2/Li2B4O7) and fused in a furnace at 1025°C. The glass was then cooled and dissolved in an acid mixture of nitric, hydrochloric and hydrofluoric acids and then analysed by ICP-MS.

For analyses at CAL, the existing pulp material as delivered by ALS was used. The pulp was re-pulverized in an agate vessel to ensure homogeneity in the sample. Major element oxides were analysed by XRF while trace elements were analysed by a multi-AD (Digestion Acid System, PicoTrace®; Yu et al. 2001) followed by ICP-MS.

This two-step digestion process is specifically designed to dissolve resistate minerals (Yu et al. 2001). This is achieved by digesting the sample at high temperatures and pressures and extending the digestion times to allow for all the resistate minerals to dissolve. This process is not practicable for high-throughput commercial laboratories and is not a standard method undertaken by ALS. Further method details can be found in the Supplementary material.

Data obtained for the eight elements being assessed (Zr, Hf, Nb, Ta, Ti, V, Sc and Eu) by the three methods are presented in Figures 2, 3 and 4 and in the Supplementary data. Summary statistics for the 4AD and AFAD data are presented in Table 2. Ten samples analysed at CAL were chosen from those analysed by 4AD and AFAD to provide an independent dataset (Table 1). These were chosen to cover a range of values as determined by the different methods. There is generally good consistency between the CAL and AFAD results; however, several outliers are notable.

Sample 143563 had consistently higher AFAD values compared to the CAL value in all eight elements. Sample WSP5_01 had different vanadium values for all three methods with 1.2, 7.12 and 18 ppm for the 4AD, CAL and AFAD methods respectively. Sample 143557 reported >20% variation in the AFAD value and 143544 reported >20% variation in the 4AD value from the CAL value for Sc. The CAL results for Nb, Ta and Ti were all higher in sample HEC6637_01. All outliers have concentrations outside the range of error for the method. This lack of agreement between methods is therefore best explained due to heterogeneities in the pulverized samples. Re-pulverization of the CAL samples may have produced smaller, more digestible grains, improving recovery.

V, Sc and Eu performed the best of all elements as indicated by a 1:1 correlation between the AFAD and 4AD data (Fig. 2). The slope of the line of best fit for V, Sc and Eu is 1.06, 1.04 and 1.02 respectively and all three elements have an R2 value of 0.99 (Table 3). Although these illustrate a general indication of how well the methods performed relative to each other, the gradient values can be misleading as they are more influenced by higher concentrations and can ‘overlook’ variations at low concentration. This is evident for the Nb and Ta diagrams which appear to fit well to a 1:1 correlation with gradients of 1.0 and 0.98 respectively; however, there is significant variation at lower concentrations (Fig. 3a–e). Additionally, gradients <1 (i.e. recovering more from the 4AD method relative to the AFAD method) are not considered to be real effects and are within the method error.

For V, consistency between the methods deteriorates at low concentrations as can be seen in Figure 2b. This is likely a result of the higher detection limits associated with the AFAD method. This effect is also noted in the Sc results, but not in the Eu results owing to sufficiently high sensitivity for Eu in both methods. At higher values (>150 ppm), V starts to show a systematic bias towards the AFAD method. The CAL results are consistent with both methods, falling within 10% of each other. The only exceptions to this were samples 184302 and 137997, which at 3.8 and 2.6 ppm, had the lowest concentrations and are likely approaching the sensitivity limit of the 4AD method.

Sc results are consistent among all three datasets with differences between AFAD and CAL data within 5%. Figure 5 demonstrates that there is very little change in whole-rock V/Sc between the two methods, with the average V/Sc values differing by <5%.

Eu data perform moderately well with 73% of the data falling within 10% of the 1:1 line (Fig. 2g, h). Beyond this range, 13% was biased towards AFAD and 13% was biased towards 4AD. The CAL data are more consistent with the AFAD method with most values fitting within 10% of each other, whereas half of the 4AD samples were >10% lower in concentration than the CAL data. This suggests the bias towards AFAD in the Eu data is a result of under-reporting in the 4AD data.

A value of the amount under-reported or, ‘UR%’ can be calculated as:
(CAFADC4AD)100/CAFAD.
(1)
This is a measure of the percentage of mass that is not recovered in the 4AD analysis relative to the AFAD analysis. Applying this to V and Sc data, the average UR% is 12.7% and 8.4% respectively. This improves to 9.6 and 7.25% respectively when V AFAD values <10 ppm are filtered out to account for the lower sensitivity of the AFAD data.
Calculating the Eu UR% shows that individual samples could have as much as 58% of Eu not recovered from the 4AD method (Fig. 6). However, the average under-reporting overall is only 0.09%. Despite the large deficit in recovery for individual samples, the Eu anomaly, Eu/Eu* (after Taylor and McClennan 1985), shows no bias with recovery with <2% difference between the average AFAD and 4AD Eu/Eu* values (Fig. 6). This might be explained by how Eu* is calculated:
Eu=SmNGdN
(2)
where normalization factors 0.087, 0.231 and 0.306 were used for Eu, Sm and Gd respectively (after Taylor and McClennan 1985). Since Eu, Sm and Gd would be expected to be hosted in the same minerals (i.e. zircon), any loss in recovery would be expected to similarly affect all three elements. This produces a covariation at approximately 1:1 which preserves Eu/Eu* in the 4AD method.

Nb, Ta and Ti all behave in a similar manner (Fig. 3). The data show larger deviations from 1:1 compared to V, Sc and Eu resulting in lower R2 values (Table 3). Both Nb and Ta have c. 70% of the data within 10% of the nominal 1:1 line (Fig. 3a–e). Titanium has only 59% within 10% of the 1:1 line, with most of the values outside of this range strongly biased (98%) towards the AFAD method (Fig. 3g, h). This is also seen in the Nb and Ta plots, but to lesser extents (Fig. 3a–e). As with the Eu data, the consistency between the AFAD and CAL data is suggestive of under-reporting in the 4AD method.

The AFAD values for sample 137997 were c. 16 and 24% below the CAL values in both Nb and Ta respectively. This variation was not seen in the Ti data. The concentrations required to bring the values back within the 10% error range are 2.5 and 0.7 ppm for Nb and Ta respectively. Variations at this scale are well within the range of uncertainty for the method and are not considered true outliers.

Under-report percentage calculations revealed that the maximum UR% for Nb, Ta and Ti was 71% for all three elements, with averages of 3, 3 and 10% respectively. The ratio of the UR% between all three elements is linear at an approximate 1:1 ratio (Fig. 7a). This appears to suggest that a single mineral phase is controlling the recovery of these elements (e.g. rutile or titanite).

Figure 7b demonstrates the relationship between whole-rock Nb/Ta and the Ti UR%. Despite some samples not recovering up to 70% of their Nb, Ta and Ti, the effect on whole-rock Nb/Ta is negligible. This is expected given that the recovery of these elements covary at an approximate one-to-one ratio. This suggests that lithogeochemical techniques utilizing Nb/Ta or Ti/Nb ratios in this dataset will be valid regardless of whether the AFAD or 4AD method is used.

Zr and Hf data have the largest differences between the AFAD and 4AD methods (Fig. 4). This is consistent with these elements being dominantly hosted in highly resistant zircon. A conspicuous bias towards the AFAD values in both elements is noted with more than half of all data falling outside of the 10% error range. Only five samples (1.2% of the data) were biased >10% towards the 4AD method and are considered within error. Comparisons of the AFAD and CAL data show that most samples are consistent and the 4AD data are under-reporting, sometimes significantly, with up to 85% of the Zr in the sample not being recovered. In fact, the summed concentrations from 4AD analysis for Zr and Hf for all 410 samples is 76 and 82% of the AFAD value respectively.

When comparing the UR% between Zr and Hf, there is more Zr under-reporting relative to Hf (Fig. 8a). This has a drastic effect on whole-rock Zr/Hf values which show significant bias towards low Zr/Hf with increasing Zr UR% (Fig. 8b). This renders the 4AD method unfit for purpose when assessing Zr-Hf ratios for lithogeochemical data interpretations.

These results are consistent with those of Pashkova et al. (2019) who found the 4AD method to be ineffective in fully recovering Nb, Ta, V, Zr, Cr and Hf due to resistate accessory minerals in ultramafic meimechites when compared to XRF results. Underestimation of the concentration by 4AD were reported as much as 80% for Nb, 70% for Ta, 50% for V, 40% for Zr and Hf and 20% for Cr. The porphyritic nature of the rocks, as well as the presence of resistate minerals such as Cr-spinel and perovskite were believed to be the main cause in the low recovery of the 4AD (Pashkova et al. 2019).

In each case presented here where a bias was determined, the bias was found to be due to lower recovery in the 4AD method. The elements Sc and V were found to have no significant bias between the methods. While there are several causes for under-estimation in the 4AD method, the most cited cause is the presence of acid-resistant minerals. Minerals rich in HFSE such as zircon, baddeleyite, tantalite, columbite, rutile, ilmenite, titanite, garnet and monazite have long been documented in the literature for their resistance to the 4AD procedure (Ito 1962; Hoops 1964; Jarvis 1990; Chao and Sanzolone 1992; Eggins 2003; Senda et al. 2014; Potts et al. 2015). The AFAD method is therefore the recommended method for the accurate determination of Zr, Hf, Nb, Ta, Ti, Eu and Eu* in these samples.

Having established the improved recovery and therefore increased accuracy of the AFAD over the 4AD method for Zr, Hf, Nb, Ta, Ti and Eu, a brief lithogeochemical interpretation of the MRV dataset utilizing these elements will demonstrate their utility and their shortcomings when used with the 4AD method.

Figure 9 shows how plagioclase fractionation can be identified and relative oxidation states might be inferred between different magma compositions. For the MRV data there is minimal change in the Eu anomaly at higher Sc concentrations. At Sc < 10 ppm, a change in composition is noted where the negative Eu anomaly begins to rapidly increase as the rocks become more depleted in Eu. Within this range, a trend to a positive Eu anomaly at low Sc values is dominated by crystalline samples (e.g. porphyry, granite) whereas a negative Eu anomaly is mostly characteristic of volcanic rocks. This suggests that the most evolved felsic volcanic compositions were derived from magmas that precipitated significant amounts of plagioclase, removing Eu from the melt during crystal separation. This also has implications for the oxidation state of the magma reservoir at this point and suggests a more reduced environment with more pronounced negative Eu anomalies (Drake and Weill 1975; Burnham et al. 2015).

Halley (2020) demonstrated the value of plotting elements such as Sc and Ti against Th, Nb and Zr in discriminating magma compositions as well as the use of Sc as a general indication for maficity, where mafic rocks generally have more Sc than felsic rocks.

A plot of Ti/Nb v. V/Sc is effective in not only discriminating the bulk compositional populations (e.g. mafic v. felsic) but, under certain conditions, can also indicate whether fractional crystallization has occurred and what the temperature and pressure conditions might have been at the time of crystallization. Figure 10 demonstrates how at V/Sc values >5, the Ti/Nb value becomes more variable with over an order of magnitude difference in values for the same V/Sc value. At V/Sc values <5, a gradual decrease in Ti/Nb and V/Sc is noted. To understand what is controlling this relationship, a plot of Sc v. SiO2 is required (Fig. 11).

Only 98 samples in the dataset have both AFAD and 4AD analysis in addition to SiO2 (analysed by XRF at MRT on splits from the same samples used here). This dataset follows the trend observed by Halley (2020) who noted that Sc decreases with increasing SiO2. In Figure 11, the highest Sc values are generally associated with lower silica values, consistent with mafic compositions and lower Sc values are generally consistent with higher silica, felsic compositions. The subset coloured in orange represents all samples with Sc values <10 ppm. These are classified as ‘high silica’ and represent the most evolved rocks in the dataset. This explains the trend towards low Ti/Nb and low V/Sc in Figure 10. Since Nb is more incompatible than Ti, it will enrich relative to Ti in evolved magma compositions. Additionally, Sc is removed during magma evolution, resulting in an overall reduction in Ti/Nb and V/Sc (Fig. 10).

We can now use this 10 ppm Sc cut-off to classify the remainder of the data that do not have SiO2 assays as a proxy for ‘high-silica’ evolved compositions. Figure 12 compares Ti/Nb v. Zr/Hf between the AFAD (Fig. 12a) and 4AD (Fig. 12b) digestion methods. The ‘high-silica’ samples are coloured orange and intrusive samples are represented by open diamonds. Two features are noted in this diagram. Firstly, orange, low-Sc samples are also low in Ti/Nb and Zr/Hf. This can be interpreted to be a result of fractional crystallization. However, this alone does not explain the ‘high-silica’ results at higher Ti/Nb and Zr/Hf and hints at additional controls. The other feature of note is that Ti/Nb is largely unaffected by analysis method. This result was anticipated since it has been shown that Ti and Nb recovery covaries at an approximate 1:1 ratio. The Ti/Nb ratio is therefore preserved, regardless of under reporting of the individual elements. This implies that low recoveries are the result of a single resistate mineral phase affecting both elements equally. It is assumed that rocks with a variety of Ti–Nb rich resistate minerals would not be expected to behave this way. In contrast, Zr/Hf is not preserved between the two digestion methods, with 4AD Zr/Hf having lower values than AFAD Zr/Hf for the same sample (Fig. 12b). Again, this is expected based on the Zr and Hf UR% values (Fig. 8a, b) which demonstrated a bias towards greater Zr under-reporting, resulting in lower Zr/Hf values with intrusive rocks particularly affected by this (Fig. 12b). This has implications for fertility assessments based on Zr/Hf and Nb/Ta which have been used to identify highly fractionated intrusives enriched in incompatible elements (Linnen and Keppler 1997; Linnen 1998; Selway et al. 2005; Zaraisky et al. 2008; Stepanov et al. 2014; Li et al. 2015; Ballouard et al. 2016; Moreno et al. 2016; López-Moro et al. 2017; Simons et al. 2017; Gonçalves et al. 2019; Chandler and Spandler 2020). The bias towards lower Zr/Hf in 4AD analyses from Figure 12b produces false-positive anomalies, making the 4AD method unfit for purpose for the interpretation of Zr/Hf. The corollary to this conclusion is that partial digestion of zircon might actually amplify the whole-rock Zr/Hf signal by only digesting the Hf-rich rims. However, without knowing whether partial digestion has occurred or what extent the zircons were dissolved, Zr/Hf cannot reliably be used in 4AD methods.

Additional indicators for magmatic fertility include La/YbN and Sr/Y fractionation (Richards 2011; Loucks 2014; Jamali 2017; Möller and Williams-Jones 2017; Elliott 2018; Zozulya et al. 2022) whose recovery can be heavily affected by resistate mineralogy such as zircon and garnet (Sholkovitz 1990). Preliminary results indicate that the 4AD method is also prone to producing artificially high Sr/Y ratios, likely due to resistate mineralogy. This has implications for porphyry copper exploration which utilize Sr/Y to identify fertile intrusions (Richards 2011; Loucks 2014).

To investigate the changes to rock classification with analysis method, a volcanic rock classification diagram can be used (Fig. 13; Pearce 1996). This figure shows the 4AD data coloured according to the Pearce classification using the AFAD data. Polygons represent the data extents for each composition from the AFAD data. Black arrows indicate the translation of the classification mean from 4AD to AFAD. This figure shows that, while there is broad agreement between the main compositional classifications, there appears to be a trend in the 4AD data towards more mafic and alkaline classifications compared to the AFAD data. This classification diagram performs poorly with the 4AD data because the UR% for Zr and Ti do not covary (unlike Ti and Nb). This relates to the resistate mineralogy where different minerals are independently controlling recovery such as zircon for Zr and Y and rutile for Ti and Nb.

These conclusions are in agreement with Pashkova et al. (2019) who found that lithogeochemical interpretations of the geodynamic settings of whole-rock samples with known resistate mineralogy were drastically different depending on whether a fusion-based method or an AD method was used. This further highlights the importance of using appropriate analytical methods that are fit for purpose for the intended application.

Figure 13 highlights the potential pitfalls of using trace element discriminant plots produced from 4AD methods that rely on immobile trace elements that tend to concentrate in resistate minerals that are impacted by low recovery. For more accurate composition classifications using HFSE, fusion-based methods such as AFAD are required. If, however, only 4AD data are available, discrimination plots using immobile trace elements that experience similar levels of under-recovery (i.e. Nb and Ti) and elements that are not concentrated in resistate phases (i.e. Sc and V) can be used (Fig. 10).

The comparison between 4AD and AFAD methods in interpreting geochemical data highlights how some of the most informative elements, and their ratios, used in lithogeochemistry are also the most susceptible to under-reporting with the 4AD method. Improper use of discrimination plots using 4AD data can lead to overestimation of maficity and alkalinity as well as potentially produce false-positive anomalism during geochemical fertility assessments. These effects can be attenuated with larger-sample populations as shown by the relatively small translations of the population centroids between the methods (Fig. 13). For individual samples, the bias can be large with as much as 45% difference in Zr/Hf and 81% difference in Zr/Ti values between the two methods. Careful selection of immobile trace elements can also improve the accuracy of the discrimination plots when using 4AD data. Elements such as Nb, Ta and Ti can suffer from incomplete recovery in the 4AD method but can preserve their ratios when controlled by a limited number of resistate minerals. Ratios of these elements will minimize the bias of the incomplete digestion on the discrimination plot, especially when larger datasets are used for interpretation. Caution is recommended when attempting to use immobile trace element plots on the 4AD data where only small sample populations are available.

HFSE such as Nb, Ta, Ti, Zr, Hf and Eu are commonly used in lithogeochemical interpretations to make classifications of rock composition, tectonic setting and magmatic fertility. These same elements are also commonly found in resistate accessory minerals and, as a result, are prone to incomplete recovery in some AD methods such as 4AD. Partial recoveries in the 4AD relative to samples prepared by AFAD were noted for the elements assessed in this study except Sc and V which were largely unaffected by digestion method. On a sample-by-sample basis, these differences were significant with under-reporting up to 85% for Zr, 77% for Hf, 58% for Eu and 71% each for Nb, Ta and Ti, confirming that these elements have incomplete recovery when analysed by the 4AD method. Despite the significant under-reporting for individual elements, certain element ratios are preserved between the analytical methods while others are not. Elements that have similar levels of under-reporting, such as Nb, Ta and Ti, will preserve their element ratios between analytical methods resulting in discrimination plots which can be used in both the 4AD and AFAD data. Plagioclase fractionation modelled using the Eu anomaly behaves similarly since it is calculated as an interpolant of Gd and Sm for which under-reporting levels are similar to Eu. In contrast, Zr is preferentially lost relative to Hf. This results in a bias towards lower Zr/Hf in the 4AD data and can produce false-positive anomalies in magmatic fertility assessments. Other common fertility ratios such as La/YbN and Sr/Y also appear to exhibit a bias in 4AD making it not fit for purpose in utilizing these ratios.

Commonly used discrimination diagrams such as the Pearce (1996) Zr/Ti–Nb/Y plot produce inaccurate classifications when using 4AD data, resulting in samples being classified as more mafic and alkaline compared to AFAD results. These effects are attenuated with larger datasets.

For best practice lithogeochemistry utilizing large portions of the periodic table, multiple digestion methods including combinations of XRF, AD and fusion-based methods are recommended. When only 4AD data are available, plots that utilize elements with similar levels of under-recovery, such as the Ti/Nb v. V/Sc diagram, are recommended since these element ratios are preserved with digestion.

The authors would like Mineral Resources of Tasmania (MRT) for funding the analytical work and to Australian Laboratory Services (ALS) for providing the analytical work. Additionally, we would like to thank the following people for donating their time to provide feedback on draft versions of this manuscript; Ben Cooke, Amanda Stoltze, Mark Pirlo, Michele Ramshaw and Vijay Singh (ALS), Simon Gatehouse (BHP) and Andrew McNeill (MRT).

ZZ: conceptualization (equal), data curation (lead), formal analysis (lead), investigation (lead), methodology (equal), project administration (lead), visualization (equal), writing – original draft (lead); LD: conceptualization (equal), formal analysis (supporting), methodology (equal), supervision (lead), validation (lead), visualization (equal), writing – review & editing (lead); SH: conceptualization (supporting), visualization (supporting), writing – review & editing (supporting); SB: conceptualization (equal), funding acquisition (lead), methodology (supporting), supervision (supporting), writing – review & editing (supporting); MB: supervision (supporting), writing – review & editing (supporting)

This work was funded by the Mineral Resources Tasmania.

Funding for whole-rock analysis was provided by the State Government of Tasmania, Mineral Resources Tasmania (MRT).

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

All data generated or analysed during this study are included in this published article (and if present, its supplementary information files).

This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (http://creativecommons.org/licenses/by/4.0/)