A new model for zircon saturation in silicate melts

A new model describing zircon saturation in silicate melts is presented that combines the results of 196 data from new experiments with data from previous experimental studies. In the new experiments, the concentration of Zr in melts coexisting with zircon was determined at temperatures between 800 and 1500 °C for 21 compositions (with alumina saturation index, ASI, from 0.20 to 1.15), containing ~ 1 to 16 wt % FeOT and, for a subset of these conditions, at variable pressure (0.0001 to 4.0 GPa) and water content (0 to 15 wt %). The collated dataset contains 626 data, with 430 from 26 literature studies, and covers conditions from 750 to 1620 °C, (including 45 new data and 106 literature data for temperatures < 1000 °C), ASI 0.20 to 2.00, 0.0001 to 4.0 GPa and 0 to 17 wt % H2O. A limitation of previous models of zircon saturation is the choice of parameter used to describe the silicate melt, which may not be appropriate for all compositions and can result in differences in predicted temperatures of over 200 °C for granitic systems. Here we use optical basicity (Λ), which can be easily calculated from the major oxide components of a melt, to parameterise the composition. Using a non-linear least-squares multiple regression, the new zircon saturation model is:logZr=0.96(5)-5790(95)/T-1.28(8)P+12.39(35)Λ+0.83(9)x.H2O+2.06(16)PΛ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log{\mathrm{Zr}} = 0.96(5) - 5790(95)/T - 1.28(8)P + 12.39(35){\varLambda} + 0.83(9)x.\text{H}{_2}\text{O} + 2.06(16)P{\varLambda} $$\end{document}where Zr is in ppm, T is temperature in K, P is pressure in GPa, Λ is the optical basicity of the melt, x.H2O is the mole fraction of water in the melt, and the errors are 1σ. This model confirms that temperature and melt composition are the dominant controls on zircon solubility. In addition, pressure and melt water content exert small but resolvable effects on the solubility and are included, for the first time, in a model. Using this new calibration, 92% of the predicted temperatures are within 10% of the experimental temperatures for the collated dataset (with an average temperature difference of 57 °C), while predicted temperatures for only 78 and 62% of the collated dataset are within 10% of the experimental temperature (with average temperature differences > 80 °C) using the widely cited Watson and Harrison (Earth Planet Sci Lett 64:295–304, 1983) and Boehnke et al. (Chem Geol 351:324–334, 2013) models, respectively. This new model can be extrapolated to temperatures below those included in the calibration with greater accuracy and when applied to melt inclusions from the Bishop Tuff, gives temperatures that are in excellent agreement with independent estimates.


Introduction
Zircon (ZrSiO 4 ) is one of the primary hosts for trace elements such as U, Th, the high field strength elements (HFSE) and the rare earth elements (REE + Y), and can exert a significant control on the trace element budget of a melt during crustal melting and fractionation processes (Hanchar and Watson 2003). It is ubiquitous in a wide range of intermediate-felsic magma compositions (Hoskin and Schaltegger 2003) and is used as a geochronometer, geothermometer and for tracing the evolution of magma chemistry (Hinton and Upton 1991;Ferry and Watson 2007;Hanchar and van Westrenen 2007). As zircon is resistant to both crustal anatexis and chemical and physical weathering, it has a high probability of survival, and the oldest known terrestrial materials are Hadean zircons preserved in younger sedimentary rocks (Mojzsis et al. 2001;Wilde et al. 2001;Cavosie et al. 2006). Therefore, it is important to determine whether zircon is autocrystic or inherited when interpreting its geochemical signature (Siégel et al. 2018). For this reason, zircon has become the focus of saturation thermometry studies to understand its occurrence and stability in igneous rocks (Hanchar and Watson 2003).
The solubility of zircon in silicate melts can be determined using high-temperature, high-pressure experiments, which have found that zircon saturation strongly depends on temperature and melt composition Boehnke et al. 2013). Specifically, for a peraluminous granite composition, the concentration of Zr at which zircon crystallises is low (~ 100 ppm) at 750 °C but increases to ~ 5000 ppm at 1050 °C . The dependence of Zr solubility on melt composition is controlled by the alkali and silica contents. First, Zr saturation strongly depends on the molar proportions of (Na 2 O + K 2 O)/Al 2 O 3 with Zr concentrations increasing from 100 ppm when (Na 2 O + K 2 O)/Al 2 O 3 = 1 to ~ 4 wt % when (Na 2 O + K 2 O)/Al 2 O 3 = 2 at 800 °C (Watson 1979;Watson and Harrison 1983;Boehnke et al. 2013). Furthermore, zircon solubility in peraluminous and metaluminous melts is an order of magnitude lower than zircon solubility in more alkaline compositions (Linnen and Keppler 2002;Gervasoni et al. 2016). Second, zircon solubility increases strongly with decreasing silica content of a melt, with Zr concentrations in lunar basalts five times greater than in a melt with an average peraluminous granite composition at the same temperature (1100 °C, Dickinson and Hess 1982).
Models of zircon saturation in silicate melts have been proposed to describe peraluminous and metaluminous compositions Boehnke et al. 2013;Shao et al. 2020), alkaline to peralkaline compositions (Gervasoni et al. 2016;Shao et al. 2020), and high-temperature basaltic melts (Borisov and Aranovich 2019;Shao et al. 2019). The models of Watson and Harrison (1983) and Boehnke et al. (2013) have been used most extensively in petrogenetic studies of granitic magmatism. These two expressions were calibrated for temperatures between 750 and 1020 °C, and compositions with an alumina saturation index (ASI, the molar ratio of Al 2 O 3 /(CaO + Na 2 O + K 2 O)) between 0.90 and 1.10 (Fig. 1b). The concentration of Zr in the melts containing zircon varied between ~ 150 and 3450 ppm (Fig. 1f). Yet, the typical crystallisation temperature range of granites is 650-800 °C (Scaillet and MacDonald 2001;Schiller and Finger 2019) where Zr concentrations are typically < 150 ppm (Shao et al. 2019). Thus, the application of these models to natural conditions outside of the range for which they were calibrated may be associated with significant uncertainties due to extrapolation, and the predicted zircon saturation temperatures for typical granitic melts (ASI = 1.0) containing 50 ppm Zr are below the granite solidus (Scaillet and MacDonald 2001;Schiller and Finger 2019) using the Boehnke et al. (2013) model (Fig. 2b). Since the model of Watson and Harrison (1983), a large amount of data on zircon saturation has been published and three recent studies have attempted to refine the original  (Zr, ppm) used in each zircon saturation model. For each swarm, the data are offset (jittered) in x based on the kernel density estimate of y for the data in that model, so that the width of the swarm relates to the number of data at each value of y. For the previous studies, the experimental data of that study are shown in grey, and for 'This study' are shown in red. The data from previous studies also included in each model are shown in white. WH83 , B13 (Boehnke et al. 2013), G16 (Gervasoni et al. 2016), BA19 (Borisov and Aranovich 2019), S19 (Shao et al. 2019) and S20 (Shao et al. 2020) Page 3 of 24 71 model (Borisov and Aranovich 2019;Shao et al. 2019Shao et al. , 2020Fig. 1), yet none of these models were calibrated using data for a wide range of compositions and temperatures.
The aim of this study was to produce a model for zircon solubility for an expanded range of silicate melt compositions and temperatures by adding new experimental results to the existing zircon saturation dataset. In addition to data from the previous studies of zircon solubility, the new dataset includes Zr saturation data available indirectly from experimental studies of zircon-melt partitioning. The dataset used here contains 196 new data (see 'Experimental Dataset'), 197 data from the literature not used in previous models and 233 data used in defining the existing models. A new equation was derived, which describes zircon solubility as a function of temperature and melt composition, parameterised as optical basicity (Λ), in addition to water content and pressure, which were not included in previous models.
The previous studies have resulted in seven models, which are given as follows, and include coefficients for temperature and melt composition (variously denoted by M, G or B) in addition to a constant: Watson and Harrison (1983) (1) lnD Zr = 12900∕T − 0.85(M − 1) − 3.80 where Zr is in ppm, D is the zircon-melt partition coefficient for Zr, and T is temperature in K with the exception of Gervasoni et al. (2016) and Shao et al. (2019) where T is in °C.
The relationship between zircon saturation and temperature is well-defined by the experimental data collated from the literature (Fig. 2a excluding unfilled symbols). Nevertheless, the modelled temperature dependencies diverge significantly when extrapolated below 750 °C, which is the lowest temperature for which data was used to calibrate the models, to the conditions most relevant to granitic magmatism (Fig. 2b). For a typical granite composition (ASI = 1.1 with 50 ppm Zr) the predicted zircon saturation temperatures vary by ~ 220 °C between the seven models, with four of the seven models predicting temperatures below the granite solidus. These differences may reflect the temperature range investigated and the number of data < 900 °C used in each study (coloured symbols in Fig. 2a). The latter highlights the challenge of growing zircon experimentally (Hanchar et al. 2001;Rubatto and Hermann 2007;Burnham and Berry 2012) and several studies noted that it was unclear if equilibrium was achieved in experiments < 900 °C (Boehnke et al. 2013;Taylor et al. 2015;Gervasoni et al. 2016;Marxer and Ulmer 2019). The divergence may also reflect the different melt compositions investigated in each study and their effect on zircon saturation.
(2) lnD Zr = 10108(32)∕T − 1.16(15)(M − 1) − 1.48 (9) (3) lnZr = 0.0056(2)T − 1.35 (10) Supplementary Table ST1) and predicted using the seven previous zircon saturation models, as a function of reciprocal temperature and alumina saturation index (ASI); a experimentally determined as a function of 1/T, b predicted as a function of 1/T for a granitic composition with ASI = 1.1, c experimentally determined as a function of ASI, and d predicted as a function of ASI at 1100 °C. The region between the vertical dashed lines in b corresponds to the typical crystallisation temperature range for granites, and the dashed horizontal line represents the average Zr content of granitic melts (50 ppm; Shao et al. 2019). The vertical dashed line in d indicates the average composition of all of the data summarised in Supplementary Table ST1. In a and c the triangles represent exper-imental data from Watson and Harrison (1983) in red (29), Boehnke et al. (2013) in blue (25), Gervasoni et al. (2016) in green (14), Borisov and Aranovich (2019) in black (69), Shao et al. (2019) in yellow (14), and the grey circles (293) represent the data from the other zircon saturation studies in Supplementary Table ST1. The values in parentheses indicate the number of data. Unfilled circles (254) represent data from experiments in Supplementary Table ST1 that were excluded from the model derived in the present study (see 'Model fitting'). The colours of the curves in b and d are the same as those of the data in a and c from which they were obtained, with the addition of Shao et al. (2020) using G (dashed green), and Shao et al. (2020) using M (dashed blue). Note that the Gervasoni et al. (2016) and Shao et al. (2019) models have the same compositional dependency in d Comparing the compositional dependence between the seven models is difficult, as the models use different parameters to describe melt composition. The Watson and Harrison (1983) and Boehnke et al. (2013) models use M, which was specifically designed to describe peraluminous and metaluminous compositions (ASI 0.90-1.10) and is given by the cation molar ratio (Na + K + 2Ca)/(Al × Si) in the melt. To address alkaline melt compositions (ASI < 1.00), Gervasoni et al. (2016) proposed G as the melt parameter, which is defined as the molar ratio (3Al 2 O 3 + SiO 2 )/(Na 2 O + K 2 O + CaO + MgO + FeO T ), and is similar to NBO/T (the ratio of non-bridging oxygens to tetrahedrally coordinated cations; Mysen et al. 1982;Mysen 1990). This parameter was also used in the Shao et al. (2019) and (2020) models. Finally, Borisov and Aranovich (2019) proposed the melt parameter B, given by: 0.14(TiO 2 /SiO 2 ) + 1.3(CaO/ SiO 2 ) + 1.5(Na 2 O/SiO 2 ) − 4.5(K 2 O/SiO 2 ) − 2.7(Al 2 O 3 / SiO 2 ) + (MgO/SiO 2 ) 2 − 3.7(CaO/SiO 2 ) 2 + 75(K 2 O/SiO 2 ) 2 , where all oxides are in moles. The compositions used in the previous zircon saturation models were converted to ASI for enabling comparison in Fig. 2c, where zircon saturation increases with decreasing ASI. Figure 1b shows that the bulk of the data used in previous zircon saturation models cluster near ASI = 1.0, yet the compositional relationships modelled at 1100 °C (where the temperature relationships predicted by the different models agree most closely; Fig. 2b), differ greatly for each model (Fig. 2d). Therefore, it may be possible to identify a different melt parameter that accurately describes zircon solubility for all compositions, and to constrain the temperature dependence over a wider range of conditions, by collating and modelling all the zircon saturation data that is available.
The roles of water, halogens, Fe and oxygen fugacity (fO 2 ) on zircon saturation have also been investigated Keppler 1993;Baker et al. 2002;Linnen 2005;van Lichtervelde et al. 2010;Aseri et al. 2015) but these variables have not been included in any model due to different studies suggesting different effects. Furthermore, zircon-melt partitioning studies, although not directly investigating zircon saturation, provide further data on Zr solubility in silicate melts (Keppler 1993;Linnen and Keppler 2002;Rubatto and Hermann 2007;Luo and Ayers 2009;van Lichtervelde et al. 2011;Burnham and Berry 2012;Stepanov et al. 2012;Trail et al. 2012;Hofmann et al. 2013;Taylor et al. 2015;Marxer and Ulmer 2019;Trail 2019, 2022;Chowdury et al. 2020). These studies are summarised in Supplementary Table ST1 and the experimental conditions cover: 700-1450 °C, ASI of 0.50-1.77, 0.0001-3.0 GPa and 0-17 wt % H 2 O, providing an additional 238 data. A wealth of additional data is available to better constrain zircon saturation models.

Starting materials
Zircon saturation was determined by measuring the Zr concentration of quenched glass after crystallisation of zircon. Twenty-one starting compositions were prepared, corresponding to a range of crustal melt compositions with ASI from 0.20 to 1.15 (Table 1). These include seven compositions that were prepared along a mixing line between an andesitic composition (AB44, Burnham and Berry 2012) and a granitic composition (1690A, Rubatto and Hermann 2007), nine compositions from the CaO-MgO-Al 2 O 3 -SiO 2 (CMAS) system that have low liquidus temperatures (O'Neill and Eggins 2002), three modified CMAS compositions containing Fe, and two compositions corresponding to the average whole-rock compositions of the I-and S-type granites of the Lachlan Fold Belt from Burnham and Berry (2017).
All compositions were prepared by mixing dried reagentgrade powders of SiO 2 , Al 2 O 3 , MgO, CaCO 3 , Ca 3 (PO 4 ) 2 , Na 2 CO 3 , K 2 CO 3 , TiO 2 , and FeO under acetone. The resulting powders were pressed into pellets and fired at 1000 °C for 24 h to decarbonate. Zirconium was added to the oxide mixes in a nitrate solution to produce a slurry that was mixed thoroughly while evaporating to dryness to ensure that Zr was homogenously distributed throughout the sample and thus inhibit pre-nucleation of ZrO 2 . The amount of Zr added was varied according to the melt composition and temperature so not to oversaturate the composition in Zr, (as occurred in Dietrich 1968 andWatson 1979 resulting in the formation of many small, needle-like crystals; Boehnke et al. 2013), and to ensure that zircon was a liquidus phase at the temperatures investigated. The amount of Zr added in each experiment is given in Supplementary Table ST2, and was based on the model of Boehnke et al. (2013), with an additional 15% Zr to ensure that the presence of zircon could be easily confirmed. A large number of trace elements (Li, Sc, Sr, Y, Nb, REE, Hf, Ta, Th, U) were also added as commercial AAS standard nitrate solutions to investigate zircon-melt partitioning. These partitioning results will be reported elsewhere. The powder was once again pressed into a pellet and fired overnight in air at 500 °C to remove the nitrate. The hydrous compositions were initially prepared deficient in SiO 2 or Al 2 O 3 so that water could be added at this stage as silicic acid or Al(OH) 3 .

Sample preparation
The melt compositions given in Table 1 were equilibrated at between 800 and 1500 °C and between 0.0001 and 4.0 GPa. For experiments at atmospheric pressure (0.0001 GPa), the Table 1 Average glass compositions (wt %) determined by EPMA and EDS. AB44 and 1690A are from Burnham and Berry (2012) and Rubatto and Hermann (2007), respectively, and are the two end-members for compositions of a mixing line (named 50:50 etc to represent the percentage of each end-member in the mix). Twelve compositions are from the CMAS(+ Fe) system (O'Neill and Eggins 2002). I-type and S-type are the average whole-rock compositions of the I-and S-type granites from the Lachlan Fold Belt (LFB) reported in Burnham and Berry (2017).
As most of the compositions are Fe-free, fO 2 should have no effect on zircon solubility. Nevertheless, all experiments were conducted at controlled fO 2 using CO-CO 2 gas mixes at atmospheric pressure and mixtures of metal-metal oxides (Pt/PtO 2 , Pd/PdO, Ru/RuO 2 and Re/ReO 2 ) or C/CO 2 at high pressure, to investigate differences in the partitioning of redox variable elements, which will be reported elsewhere. The fO 2 conditions of each experiment are given in Supplementary Table ST2.
The Pt capsules were sectioned along the long axis and the samples equilibrated at 1 atm were lightly crushed to maximise the surface area. Both were then mounted in epoxy resin and polished.

Major oxide composition
Back-scattered electron (BSE) images were acquired using a JEOL 6400 scanning electron microscope (SEM), and the major oxide compositions of any crystals and the quenched melts were determined by energy dispersive spectrometry (EDS). For selected samples the melt composition was also determined by wavelength dispersive spectrometry (WDS) using a Cameca SX-100 electron microprobe (EPMA). The SEM used an accelerating voltage of 15 keV, a beam current of 1 nA and an analysis spot of 3 μm. Major oxide melt concentrations for each sample were averaged from up to nine analyses. Na 2 O and K 2 O counts were monitored to ensure that Na and K were not 'lost' during analysis.
The Cameca SX-100 EPMA was operated using an accelerating voltage of 15 keV, beam current of 100 nA, a 3 × 3 μm focused beam, and TAP (Kα Mg, Al, Si, Fe, P), PET (Kα Ca; Lα Na, Zr) and LLIF (Kα Ti; Lα K) crystals. No interferences occur between the elements present. The acquisition time for each analysis was 360 s and Na and K were measured in the first 30 s to minimise the effect of loss by volatilisation. The standards used were olivine (Mg Kα and Fe Kα), corundum (Al Kα), quartz (Si Kα), apatite (P Kα), wollastonite (Ca Kα), albite (Na Lα), orthoclase (K Lα), rutile (Ti Kα) and zirconia (Zr Lα). Only zircon free areas of melt > 100 μm were analysed to avoid the possibility of secondary Zr fluorescence. The water contents of the hydrous silicate melts were estimated from the difference between the total oxide concentrations and 100 wt % following the method of Devine et al. (1995), which as shown by Mare et al. (2021) gives results that agree well with those determined directly by Fourier Transform Infrared (FTIR) spectroscopy. For samples where the major oxide composition of the glass was determined by both WDS and EDS, the results agree within 7% (Supplementary Table ST2).

Zirconium concentrations
The Zr concentrations of the quenched silicate melts (glasses) were determined by LA-ICP-MS using a Compex ArF Excimer laser (193 nm) coupled to an Agilent 7700S quadrupole ICP-MS with a HelEx ablation cell and He-Ar carrier gas. The laser was operated at a pulse rate of 5 Hz and a constant energy of 80 mJ. The laser beam was apertured to a diameter of 13 μm and the sample moved beneath the beam at a rate of 1 μm/s to assess the homogeneity of Zr in the glass (Fig. 3). The sample was ablated for a maximum of 180 s with the background signal recorded for 30 s preand post-ablation. For comparison, the glasses of selected samples were also analysed by ablating for 40 s using a stationary 37 μm spot, with 20 s of background recorded either side of the ablation.
NIST-610 and -612 glasses were used as the primary and secondary standards, respectively. Glass analyses were obtained in batches of up to 12, bracketed by two analyses of each standard. At the start of each session ThO + /Th + and 44 Ca 2+ / 44 Ca + in analyses of NIST-612 were minimised to < 0.5% to ensure negligible interference from oxides and doubly charged ions. Masses 90 Zr, 91 Zr and 177 Hf were analysed for 0.01 s in addition to masses corresponding to 25 trace elements, which are not reported in this work, giving a total acquisition time in a single sweep of 0.745 s. The SiO 2 concentration of the glass determined by EDS or WDS was used as the internal standard and the data was reduced using Iolite.

Experimental dataset
The new experimental dataset comprises the Zr contents of melts coexisting with zircon for 450 samples prepared between 800 and 1500 °C, 0.0001 and 4.0 GPa, and with compositions with ASI values between 0.20 and 1.15 (Supplementary Table ST2). For all samples, euhedral zircon crystals are present in a quenched silicate melt (glass), which is homogenous in composition (Figs. 3 and 4). Zircon is the only crystalline phase, except for 14 samples prepared at < 1000 °C, which contain minor quartz, feldspar and/or enstatite, and 13 samples prepared at high temperatures, which also contain minor baddeleyite.
The Zr concentrations of each sample were determined from a minimum of three LA-ICP-MS analytical traverses, each covering a minimum distance of 100 μm (Fig. 3b). For each traverse the Zr concentration was constant with distance and the values determined from the traverses agreed within ~ 7.5%, indicating that Zr is homogenous in the glass. The Zr contents of some glasses were determined by LA-ICP-MS using both a stationary spot and the traverse method, and EPMA, and the average difference in the results is less than 9%. The data are given in Supplementary Table  ST2 and shown in Supplementary Fig. A1. Only Zr concentrations determined by the LA-ICP-MS traverse method were used in the modelling that follows.
Although the water content of the melt was not measured directly, all samples were weighed pre-and post-run to check for water loss. The differences between the EPMA totals and 100% are consistent with the amount of water added to the starting compositions.
Of the 450 samples prepared, 406 were Fe-free, and so fO 2 should have no effect on the solubility of zircon in these compositions (see 'Model fitting' and 'FeO content and fO 2 '), as concluded by Burnham and Berry (2012) and Borisov and Aranovich (2019). The difference in Zr concentrations of glasses with the same composition prepared at constant T and P but different fO 2 is ~ 10% (only slightly more than the intra-sample variability), and this difference does not vary systematically with fO 2 . Therefore, Zr concentrations were averaged from samples with the same composition equilibrated at multiple fO 2 s but constant T and P. For example, the Zr contents of the melts of the eight samples with composition AB44 equilibrated at 1400 °C and 1.0 GPa, which were prepared at eight different fO 2 s (Supplementary Table ST2), were averaged to give one value, which is reported in Supplementary Table ST3, where the uncertainty (σZr ppm) is the standard deviation of the average. The resulting 196 data (Supplementary Table ST3), which include all 450 samples, show correlations with temperature, melt composition, pressure and water content (Fig. 5).
The correlation between logZr and reciprocal temperature is in excellent agreement with that determined previously (Fig. 5a). At constant temperature and pressure, logZr is linearly correlated with the melt composition parameterised as optical basicity (Λ; Fig. 5b). Λ is calculated as: Σm i n i Λ i /Σm i n i , where m i is the number of moles of each oxide, n i is the number of oxygens in the oxide (e.g. 3 for Al 2 O 3 ) and Λ i is the optical basicity coefficient of each oxide from Table 4 of Leboutellier and Courtine (1998). The optical basicity for all Fe-bearing melt compositions was calculated Fig. 3 a Back-scattered electron (BSE) image of zircon crystals (light grey/white) in silicate glass (dark grey) (sample C5350b). A LA-ICP-MS traverse is visible to the right of the red line. b Time-resolved LA-ICP-MS profile of the silicate melt obtained using the traverse method, where the sample was moved beneath a 13 μm laser spot at 1 μm/s so that time (s) is equivalent to distance (μm). Counts per second (cps) for selected elements are shown. The start and end of ablation are indicated by the dashed lines assuming Fe 3+ /ΣFe = 0, (where ΣFe = Fe 2+ + Fe 3+ ), as the Fe contents of the literature studies were usually reported as FeO T wt %. This is a valid assumption as the fO 2 of most Fe-bearing experiments was near quartz-fayalite magnetite (QFM), where Fe 2+ is the dominant oxidation state (Kress and Carmichael 1991). The correlation with optical basicity will be compared to those with other melt parameterisations (NBO/T, ASI, M, G and B) in 'Model fitting'. At constant temperature, pressure and composition, the Zr content increases with increasing water content (Fig. 5c), and at constant temperature and composition, the Zr content decreases with increasing pressure (Fig. 5d).

Model fitting
To produce a new zircon saturation model, the data from this experimental study was combined with data from the literature studies summarised in Supplementary Table ST1. The experimental conditions, melt compositions and Zr contents used in the modelling are given in Supplementary Table  ST3. All data from Dietrich (1968), Watson (1979) and Luo and Ayers (2009), and data from experiments ≤ 1000 °C using compositions with ASI < 0.70 from Keppler (1993), Linnen and Keppler (2002), Linnen (2005), van Lichtervelde et al. (2010), Trail et al. (2012) and Aseri et al. (2015), were excluded from the fitting as the Zr contents of these samples differ significantly from those expected for the dependence on temperature defined by the bulk of the experimental results (unfilled symbols in Fig. 2a). Further, the zircon saturation temperatures predicted for these samples by the existing models are all in very poor agreement with the actual temperatures (differences > 500 °C; Supplementary Fig. A2). The resulting collated dataset (new experimental results and experimental data from the literature) comprises 626 data covering 750-1620 °C (with 45 new data and 106 data from previous studies for temperatures < 1000 °C), ASI of 0.20-2.00 or Λ 0.513-0.611, 0.0001-4.0 GPa, and 0-17 wt % H 2 O (Fig. 1). Of these 626  Table ST2) showing euhe-dral zircon (light grey/white) in silicate glass (dark grey). The crystal size decreases with increasing silica content and decreasing temperature. Minor enstatite and/or quartz is visible in g, h and i data, 485 were for hydrous samples but the water contents of only 94 samples were determined directly i.e. by FTIR or SIMS (Baker et al. 2002;Linnen 2005;Hofmann et al. 2013;Zhang and Xu 2016). For the other samples, the water contents were assumed to be those reported for the starting compositions (van Lichtervelde et al. 2010(van Lichtervelde et al. , 2011Trail 2019, 2022;Chowdury et al. 2020) or calculated as the  Table ST3). In b the filled circles represent the new experimental data and open triangles data from the previous studies at 1400 °C (red), 1200 °C (blue), and 1000 °C (black), and the open circles represent the remainder of the new experimental data. In c the filled symbols represent data from experiments conducted at 1300 °C where water content was varied but Λ was kept constant and the open symbols represent the remainder of the 626 data. In d the coloured circles represent data from experiments at 1400 °C using compositions AB44 (red), 50:50 (blue) and 1690A (black), and the open triangles represent data from previous studies conducted at 1200 °C (blue) and 1000 °C (black). Errors (1σ) are smaller than the size of the symbols difference between analysed SEM, EPMA or SIMS oxide totals and 100%.
To model the effect of composition on zircon saturation, it is necessary to select a parameter to describe the melt. Previous models used M, G and B (see 'Previous studies and models of zircon saturation in silicate melts'). M as well as ASI do not include all the compositional components of a melt and thus may be overly simplistic to be used for all compositions. On the other hand, G, Λ and B include all major components of the melt, but G describes the melt structure, Λ describes the melt chemistry (Duffy 1993), and B was found empirically (Borisov and Aranovich 2019). Preliminary fitting of Zr concentration to melt composition only included data from samples in Supplementary Table  ST3 prepared at 1300 °C and 1.0 GPa to isolate the effect of composition on Zr solubility (Fig. 6). The results were compared for melts parameterised by M, G, B, ASI, NBO/T and Λ, and the best correlation between Zr content and melt parameter was obtained for Λ (Table 2). Consequently, Λ was used to parameterise composition in the model that follows.
The 626 data were fit by considering all possible combinations of polynomials and coefficients for all variables using a least-squares, non-linear multiple regression, and ranking the likelihood of each model using Bayes Factors (K). Bayes information criteria (BIC) were calculated for each model using the fit to the data, the weighting of each data point from σZr, the number of variables in the model and the number of observations (Kass and Raftery 1995). The model with the highest BIC is assigned a K of 1 and the ranking of other models relative to this was determined by normalising their BIC to that of the model with K = 1, such that their K describe how many times less likely that model is expected to occur relative to the top ranked model. For example, the second best model using Λ in Table 2 has a normalised BIC of 1/28 and hence is 28 times less likely than the top ranked model using Λ. When the K value for the second ranked model is large, this suggests that the second ranked model is decisively less probable.
A maximum second order polynomial was allowed to model the relationship between Zr content and each variable, and interactions between the variables were included since composition, water content and temperature were codependent in many experiments. The fitting was conducted using log 10 Zr (Zr in ppm) of melt in equilibrium with zircon, 1/T (K), P (GPa), Λ and x.H 2 O (mole fraction of water in the melt). The weighting of each datum in the fitting was 1/σZr (the standard deviation of the Zr concentration). For data where no σZr was reported, a σZr equivalent to 10% was assumed, which is equal to the average σZr of the full dataset. Uncertainties in temperature and pressure were those reported for each study, and uncertainties in Λ and x.H 2 O were determined from the standard deviations of the oxide components in the melt for the new experimental data from this study, and from the uncertainties reported for the oxide components for the literature samples (Supplementary Table  ST3). With four variables, considering linear and second order polynomial functions, and considering all possible cross terms there are 1337 potential models. The best ranked model from the fitting has R 2 = 0.94, is ranked 20 times more likely than the second best model, and is: While all possible cross terms were considered (i.e. TP, , only a term for PΛ was found to be significant. If fO 2 was included as a variable in the model, the value of its coefficient was ~ 0, as expected. The dataset included 57 halogen-bearing compositions (average of ~ 1.71 wt % F and/or Cl), from which a small positive effect on zircon solubility had been determined previously (Baker et al. 2002;Aseri et al. 2015). However, in this study a coefficient of ~ 0 was obtained when the total halogens were included as a variable.
The propagated uncertainty of the calculated Zr concentrations (σlogZr) and calculated temperatures (σT) are given in Supplementary Table ST3. It is algebraically trivial to rearrange Eq. (8) for T, where T is the zircon saturation temperature i.e. the concentration of Zr in equilibrium with zircon for a given composition and pressure. The calculated σT is proportional to the predicted temperature (Supplementary Fig. A3) and for temperatures below 800 °C is ± 48 °C, (~ 6% of the calculated value). The average σT for the total  (Table 2). All models used the same data; (samples prepared at 1300 °C and 1.0 GPa from the collated dataset). The solid horizontal line represents no difference between the actual and predicted temperatures Zr solubility dataset (626 data) is ± 96 °C. Of the total σT, the contribution due to the uncertainty in the Zr melt concentration (σlogZr, measured) is ± 15 °C, irrespective of temperature.
The average difference between the predicted and actual temperature of each experiment (T res ) is 57 °C (43 °C for the new experimental data and 64 °C for the literature data), and predicted temperatures for 92% of the 626 data are within 10% of the experimental temperature. Of the 151 data at temperatures most relevant to natural zircon-crystallising melts (< 1000 °C), 92% of predicted temperatures are within 10% of the experimental temperature. In comparison, predicted temperatures for only 78 and 62% of the 626 data are within 10% of the experimental temperature using the Watson and  and Boehnke et al. (2013) models, respectively, with T res equal to 84 and 114 °C. The results obtained using the new model presented here are compared to those from the previous models, for the same dataset, in Table 3 and Fig. 7. To facilitate use of this model, a spreadsheet is provided in the attached electronic supplement (Supplementary Table ST5) for calculating the temperature at which zircon would crystallise for a particular Zr concentration and composition.

Evidence for equilibrium
In addition to the absence of baddeleyite in the new experimental samples, zircon crystals were well-formed at all temperatures (Fig. 4), suggesting that they were in equilibrium with the melt (Trail et al. 2012;Shao et al. 2019;Borisov and Aranovich, 2019). Further evidence confirming equilibrium, particularly in the new low-temperature samples, is demonstrated by the linear relationship between logZr and 1/T (Fig. 5a). Approximately 60% of the experimental data from the collated dataset is for temperatures > 1200 °C, which is above natural conditions of zircon growth, yet the relationship between logZr and 1/T for our samples < 1200 °C is consistent with the linear trend extrapolated from higher temperatures, and essentially identical results were obtained if only the data below 1200 °C were used in the model Table 2 Statistical results from fitting Zr concentration to melt composition, parameterised in different ways (Λ, NBO/T, G, B, M and ASI), by least-squares non-linear multiple regression. The fitting only used data from samples prepared at 1300 °C and 1.0 GPa (where there were the most data; Supplementary Table ST3) to eliminate any variation in Zr concentrations due to convolution of melt composition with T or P. R 2 is the coefficient of determination; K describes how many times less likely the second ranked model using a particular melt parameter is likely to occur relative to the top ranked model using that same melt parameter, which is assigned a value of 1; Rank shows the K values when the top ranked models for each melt parameter were compared; and T res is the average temperature residual (actual − predicted) when the top ranked model for each melt parameter was applied to the data used in this fitting  Table 3 Comparison of the differences in actual and predicted temperature for the collated dataset (626 data) using different zircon saturation models; n is the number of data used to derive each model and 'Studies' is the number of studies included in each model (see Supplementary Table ST1). σT is the average uncertainty in predicted saturation temperatures from the propagated uncertainties in each model for the collated dataset and σT 800 is the average σT for data below 800 °C. No σT is reported for WH83 as no uncertainties were given for the coefficients. T res is the average temperature residual (actual − predicted), T res < 1000 is the average temperature resid-ual for all samples < 1000 °C, 10% is the percentage of the dataset where the predicted temperature is within 10% of the experimental temperature, 10% < 1000 is the percentage of the data < 1000 °C where the predicted temperature is within 10% of the experimental temperature, and ± 100 °C is the percentage of the dataset where the difference between the actual and predicted zircon saturation temperatures are within this limit. WH83   (2009) and Trail et al. (2012), zircon saturation in the excluded data was determined by the dissolution method, which is unlikely to produce zircon saturation data relevant to natural samples, due to the slow diffusion of Zr in highly silicic melts, particularly at low temperatures Boehnke et al. 2013). Due to the anomalous relationship between logZr and 1/T, and poor predicted temperatures using the existing models, these data are therefore consistent with the experiments not having reached equilibrium. For the samples from Luo andAyers (2009) andTrail et al. (2012), disequilibrium is also indicated by the REE partition coefficients, which show significantly reduced HREE/ LREE fractionation compared to the REE patterns of natural zircon and synthetic zircon from other partitioning studies. Furthermore, a model produced using these data in addition to the data in Supplementary Table ST3 had a poor overall fit (R 2 = 0.58) and only 42% of predicted zircon saturation temperatures for this dataset are within 10% of the experimental temperature (Supplementary Fig. A5).

Temperature
As found by all previous models ( Fig. 2a and b) and shown by the experimental data in this study (Fig. 5a), temperature is the dominant control on zircon solubility in silicate melts. For example, the concentration of Zr in equilibrium with zircon in a typical granitic melt (Λ = 0.525, ~ 75 wt % SiO 2 , ~ 5.5 wt % H 2 O, 0.5 GPa) is calculated using Eq. (8) to increase from 33 ppm at 700 °C to 1790 ppm at 1100 °C.
The difference between the calculated and actual temperatures, or the temperature residual (T res , actual − predicted) for the collated dataset (626 data; Supplementary Table ST3) as a function of temperature, composition, water content, pressure and Zr content are compared for each model in Fig. 8. The accuracy of each model can be assessed from the residuals; if a model perfectly reproduced the data, all points would fall along the zero line. If a parameter was not included or not properly described by a model, this would introduce significant scatter or trends. Each model was calibrated using samples prepared over a particular range of conditions ( Fig. 1; vertical lines in Fig. 8), so a particular model should be accurate within its calibration range. However, these conditions rarely include the temperatures of interest for natural zircon as it is difficult to achieve equilibrium zircon growth at low temperatures, and hence these Fig. 7 Swarm plots showing the distribution of the temperature residuals for the collated dataset (626 data) using the different zircon saturation models expressed as: a the percentage difference between the actual and predicted temperature and b T res (actual temperature − predicted temperature). For each swarm, the data are offset (jittered) in x based on the kernel density estimate of y for that model, so that the width of the swarm relates to the number of data at each value of y. The solid horizontal line in a represents 10% difference and in b represents no difference, and the dashed horizontal lines indicate differences of ± 100 °C. WH83 , B13 (Boehnke et al. 2013), G16 (Gervasoni et al. 2016), BA19 (Borisov and Aranovich 2019), S19 (Shao et al. 2019), S20a (Shao et al. 2020) using G and S20b (Shao et al. 2020)  were not variables in the previous models but the results using these models for samples prepared at various pressures and water contents are shown. The red symbols represent the data from this study (196 data) and the unfilled symbols data from the literature (430 data). In row two, the temperature residuals are plotted against the parameter used to describe the silicate melt composition in each model, and the dashed vertical line corresponds to the average composition of the data used to derive each model, which is close to ASI = 1 in all cases. For all rows, the vertical lines represent the range of each variable in each model, the dashed horizontal lines are the average ± σT of each model and the solid horizontal line indicates T res = 0 (σT is not shown for WH83 as no uncertainties were given for the coefficients). WH83 , B13 (Boehnke et al. 2013), G16 (Gervasoni et al. 2016), BA19 (Borisov and Aranovich 2019), S19 (Shao et al. 2019), S20a (Shao et al. 2020) using G and S20b (Shao et al. 2020) using M models are often extrapolated to outside their calibration range. While no residual trends are observed for the previous models when T res is plotted against temperature in Fig. 8, our model displays the least degree of scatter due to the greater range of temperatures included (Fig. 1a). This is the result of conducting experiments at high temperatures outside of the natural range of zircon growth, which allows more precise values of zircon solubility to be determined, and as a result the trend between logZr and 1/T is better defined. By then collating data over a much larger range, the extrapolations are smaller, and the resulting accuracy at conditions applicable to natural zircon growth is much improved (see 'Application to natural samples').

Melt composition
After temperature, the next most important control on zircon saturation is the melt chemistry (Figs. 2c and d;5b) and across the compositional range included in the model (0.513-0.611 Λ or 40-80 wt % SiO 2 ) the concentration of Zr required for the crystallisation of zircon calculated using Eq. (8) increases from ~ 55 to 1700 ppm at 800 °C and 0.5 GPa. For natural melts the main compositional variable is the silica content, which is negatively correlated with Λ and zircon saturation ( Fig. 9a and b). The amount of silica, along with Al 2 O 3 , is the main factor that determines the degree of polymerisation of a melt, and in mafic melts the low degree of polymerisation results in more non-bridging oxygens (NBO) to coordinate with Zr Shao et al. 2019). As a result, high Zr concentrations, which are unrealistic for natural mafic melts, are required for zircon saturation, as found by Dickinson and Hess (1982) and Borisov and Aranovich (2019). The average Zr content of MORB is 75 ppm (Sun and McDonough 1989), and Eq. (8) suggests that above 900 °C, over 1500 ppm Zr would be required to crystallise zircon. Thus, basaltic melts will tend to dissolve zircon from other sources that it comes into contact with, and have a low preservation potential for inherited zircon. Nevertheless, zircon can occur in mafic rocks (Coogan and Hinton 2006;Fu et al. 2008), and can be mantle derived (Page et al. 2007), but the melts are undersaturated in Zr and hence not in equilibrium with zircon, which has implications for the rate of transport and emplacement. Previous studies have shown that the Zr concentration required for zircon saturation in aluminous to peraluminous granitic melts (ASI 0.9-1.1) increases with alkali content (Na 2 O and K 2 O; Watson and Harrison 1983;Boehnke et al. 2013), and is even higher in alkaline and peralkaline melts (Gervasoni et al. 2016). As network modifiers, the alkalis provide NBOs to coordinate with Zr in the melt, increasing zircon solubility. CaO and MgO also act as network modifiers and the large compositional range included in this study shows that zircon solubility correlates positively with the total NBO content of the melt (Fig. 9b and c), and is greater in mafic melts with a higher overall proportion of NBOs compared to felsic melts, which have higher alkali contents (coloured scale in Fig. 9c).
The network modifiers interpretation is, however, somewhat simplistic as each modifying component affects the melt chemistry differently. For example, CaO and MgO have lower optical basicities than Na 2 O and K 2 O, and so have weaker charge stabilising abilities, meaning CaO and MgO are less likely to stabilise Zr 4+ , which has a high charge, in the melt compared to Na 2 O and K 2 O. The temperature residuals (T res ; actual − predicted) for the models used to rank the suitability of the six different melt parameters (Table 2) are shown in Fig. 6, and show that for this reason, NBO/T (NBO and NBO/T have the same T res ) and G, which include all major oxide components, are not as good at parameterising melts for this application as Λ. In terms of other melt parameters, the temperature residuals in Fig. 8, which are plotted against the melt parameter used in each zircon saturation model, show strong residual trends for M and G, even for compositions within the range for which each model was calibrated. Further, predicted zircon saturation temperatures for samples from the Bishop Tuff (see 'Application to natural samples') differ by 72 °C for the two models from Shao et al. (2020), which used the same dataset but either M or G to describe the melt composition. In contrast, there is no trend and minimal scatter when Λ is used to describe the melt (Figs. 6a and 8). Λ is a measure of the ability of a melt to stabilise the charge on a cation (Duffy 1993). It includes all major melt components unlike ASI and M, includes FeO T unlike B, and differentiates between components such as MgO and CaO, which are treated identically in the NBO/T and G parameterisation. Thus Λ is ranked as the best descriptor of the melt (Table 2; Fig. 6a).

FeO content and fO 2
The proportion of Fe in the silicate melt as Fe 3+ or Fe 2+ can affect its structure, as Fe 3+ is a network former, whereas Fe 2+ is a network modifier (Leboutellier and Courtine 1998). The Fe 3+ /ΣFe ratio of a melt at a given temperature and pressure is controlled by fO 2 . For Fe-bearing samples at constant T, Λ and FeO T content from this study, Dickinson and Hess (1982), Baker et al. (2002), Zhang and Xu (2016) and Borisov and Aranovich (2019), there is no difference in zircon solubility over ten log units in fO 2 ( Supplementary Fig.  A6), for which Fe 3+ /ΣFe varies from 0.01 to 0.99 (Kress and Carmichael 1991). The difference in zircon saturation temperature calculated using Eq. (8) for the Fe-bearing I-type composition in this study (Table 1) with 150 ppm Zr when Fe 3+ /ΣFe = 0 and 1 (Λ differs by 0.002) is only ~ 5 °C. This suggests that neither Fe 2+ nor Fe 3+ affect zircon solubility, in agreement with Borisov and Aranovich (2019). Accordingly, 71 Page 16 of 24 when fO 2 was included as a variable in the model, its coefficient was ~ 0. Baker et al. (2002) proposed that Fe reduces zircon solubility by acting as a network former at high fO 2 s, after comparing their results with those from Harrison and Watson (1983), for samples prepared at 1200 °C with Λ = 0.527 over a range of six log units in fO 2 . However, due to the low Fe content of these samples (< 0.70 wt% FeO T ), Fe has essentially no effect on Λ when Fe 3+ /ΣFe is varied from 0 to 1 (Λ for these samples varies by ~ 0.001 and so T varies by only 3 °C), so fO 2 cannot be the cause of the small variation in zircon solubility that was observed. Both pressure and melt water content differ between these studies, which using Eq. (8) could account for the differences in Zr contents.
In the calculation of Λ, FeO and Fe 2 O 3 are treated like any other oxide components and have no special significance. The dataset used here included data for 44 Fe-bearing samples from the present study, where FeO T ranged from 1.06 to 15.92 wt %, and data for a further 254 Fe-bearing samples from the literature with FeO T contents from 0.10 to 16.34 wt % (Fig. 1c). The average FeO T content of the Fe-bearing samples is 2.04 wt %. The correlation between logZr and both 1/T and Λ are similar for the Fe-bearing and Fe-free samples ( Supplementary Fig. A7), including those where significant Fe-loss from the silicate melt due to alloying of Fe with the Pt or Au 80 Pd 20 metal capsule was reported (Gervasoni et al. 2016;Shao et al. 2019), indicating that zircon saturation for Fe-bearing compositions can be predicted from results obtained using Fe-free compositions, and that Eq. (8) is applicable to Fe-bearing compositions. Fe, and whether it exists as Fe 2+ or Fe 3+ is included in the Λ parameterisation of a melt, unlike the other compositional parameters that have been discussed.

Water content
It is difficult to assess the role of melt water content on zircon solubility, as water is coupled with silica content, as well as temperature in many of the previous studies and the new experimental data. Consequently, the effect of water content on zircon solubility is masked by the effects of temperature and composition when examining the full dataset (unfilled symbols in Fig. 5c). The Λ coefficient of water has not been determined, so although it is a compositional component, the water content was treated separately. Here, we have 130 data for samples with constant composition but varying water contents from this experimental study, Harrison and Watson (1983), Baker et al. (2002), Hofmann et al. (2013), Zhang and Xu (2016) and Shao et al. (2019) (filled symbols in Fig. 5c for data from this subset at 1300 °C), which show that Zr solubility in the melt increases with melt water content. This effect is captured in Eq. (8) and is in agreement with the effect of water content found previously (Harrison Fig. 9 Zr contents of melts at zircon saturation for experimental samples from the collated dataset equilibrated at 1300 °C plotted against a and b Λ, and c NBO (non-bridging oxygens; Mysen et al. 1982;Mysen 1990). Indicative SiO 2 wt % for different Λ values are shown. In a and c, the colour of the symbols corresponds to the mole fraction of total alkalis (Na 2 O + K 2 O) and in b the NBO content. Errors (1σ) are smaller than the size of the symbols Baker et al. 2002;Linnen 2005;Zhang and Xu 2016;Shao et al. 2019). While no temperature residual trends are observed for the previous models when plotted against melt water content (Fig. 8), indicating that the effect of water on Zr solubility is small, Eq. (8) displays the least degree of scatter suggesting that including a coefficient for water improves the model. The average T res for a model using the collated dataset that does not include a term for melt water content is 69 °C, compared to 57 °C for Eq. (8). The effect of water on zircon solubility can be attributed to its effect on the melt structure, as the presence of water as OH − groups decreases the polymerisation of the melt, which increases zircon solubility (Stolper 1982).
For a range in melt water content of 0-10 wt %, the calculated zircon saturation temperatures using Eq. (8) for a granite (Λ = 0.525) at 0.5 GPa with 150 ppm Zr, decrease from 843 to 800 °C, which is similar to the propagated uncertainty (σT ≈ 45 °C). Water contents can be difficult to constrain for natural samples, and are likely to be underestimated by whole-rock compositions, so estimates of zircon saturation temperature using Eq. (8) will represent maximum temperatures if water content is ignored.

Pressure
Previous studies have indicated that the effect of pressure on zircon solubility is insignificant (Boehnke et al. 2013;Shao et al. 2019), although pressure was not investigated directly, and the range was limited (Fig. 1e). Figure 5d shows a negative pressure effect with the Zr concentrations of three compositions from this study all decreasing systematically by ~ 65% between 0.0001 and 4.0 GPa at 1400 °C. Specifically, Zr melt contents decrease from 47,150 ± 4100 ppm at 0.5 GPa to 15,300 ± 1350 ppm at 4.0 GPa for the andesitic end-member composition (AB44), from 29,120 ± 1980 ppm at 0.5 GPa to 11,750 ± 1760 ppm at 4.0 GPa for the intermediate composition (50:50), and from 10,220 ± 948 ppm at 0.5 GPa to 2250 ± 890 ppm at 4.0 GPa for the granitic end-member composition (1690A). In addition, samples prepared at constant T and Λ from Zhang and Xu (2016) and Shao et al. (2019), show a pressure effect equivalent to a decrease in zircon solubility of between 20 and 37% over 1.0 GPa (unfilled triangles in Fig. 5d). Boehnke et al. (2013) conducted experiments at 1.0 and 2.4 GPa but observed no difference in Zr contents between their samples and those from Watson and Harrison (1983) at 0.6 GPa with the same composition and temperature (Λ ~ 0.537 and 930 °C). Zr contents of samples from Rubatto and Hermann (2007) with a similar composition at 900 and 950 °C, and 2.0 GPa were ~ 50% lower than those from both of these studies, but Boehnke et al. (2013) argued that this was due to a matrix mismatch between NIST-612 and hydrous experimental glasses affecting the calibration of the LA-ICP-MS data, rather than a pressure effect. However, a matrix mismatch cannot be argued for the samples with composition AB44 from this study or for those from Zhang and Xu (2016), which were anhydrous. As water content was varied in addition to pressure in the experiments by Boehnke et al. (2013), this may be the reason that a pressure effect was not observed.
The pressure effect shown in Fig. 5d is captured in Eq. (8), and the PΛ term suggests that the effect of pressure is greater for less polymerised melts i.e. high Λ. The temperature residuals in Fig. 8 indicate that including a pressure term in the zircon saturation model improves predicted zircon saturation temperatures for the collated experimental dataset. A systematic residual trend in the temperature residuals with pressure is observed in the collated dataset for the literature models, which increasingly overestimate temperatures for samples > 1.0 GPa. For example, saturation temperatures for the new experimental data from this study > 1.0 GPa and the data from Rubatto and Hermann (2007) are predicted within ± < 50 °C using Eq. (8) compared to ± > 100 °C using all previous models. This indicates a dependence of zircon saturation on pressure, and the average T res for a model using the collated dataset that does not include a coefficient for pressure is 81 °C, compared to 57 °C for Eq. (8).
For a granite (Λ = 0.525) with 150 ppm Zr, the zircon saturation temperature predicted using Eq. (8) changes from 801 °C at 1 km to 842 °C at 30 km (a pressure change of ~ 1.0 GPa, corresponding to a representative thickness of the continental crust), which is similar to the propagated uncertainty (σT ≈ 45 °C). Therefore, any uncertainty in pressure for the common application of the zircon saturation model to granitic systems is unlikely to significantly affect calculated zircon saturation temperatures, although these will be minimum temperatures if pressure is excluded. Ignoring the effect of pressure > 3.0 GPa, however (> 100 km depth; ultrahigh pressures), would result in errors > 130 °C (T = 934 °C at 100 km for a granite with Λ = 0.535 and 150 ppm Zr). This negative pressure effect is interpreted to reflect increasing polymerisation of the melt with increasing pressure as the coordination of Al 3+ increases from four to six-fold and the NBO fraction decreases (Lee et al. 2004;Boehnke et al. 2013;Zhang and Xu 2016), which promotes the formation of zircon. Equation (8) indicates that a melt, which saturated with zircon at depth will become undersaturated on ascent leading to decompression induced zircon dissolution.

Comparison to previous zircon saturation models
By collating all of the available data on zircon saturation in silicate melts (626 data from 27 studies, including the present), Eq. (8) was calibrated using data for a greater range of conditions (temperatures from 750 to 1620 °C, compositions  Supplementary Table ST3 (626 total). Of the 626 data in each panel, the experimental data used to derive that model are shown in black (29,41,55,129,86,168, and 168 data, respectively for each panel, with the exception of h). Of the remaining data, the new experimental data are in red, and the literature data are in white. In h, the new experimental data is shown in red to distinguish it from the literature data, even though all 626 data was used to derive the model. For clarity, a small number of data in b, d and g are not shown, as the predicted temperatures were > 2000 °C. The blue triangles represent data from the Bishop Tuff (Supplementary Table ST4), which were not included in the model fitting. Error bars are 1σ and if not visible are smaller than the size of the symbol. No error bars are shown in a as no uncertainties were reported for the coefficients in this model. The solid line is 1:1 and the shaded area represents the 10% uncertainty range of each model with 40 to 80 wt % SiO 2 , or Λ from 0.513 to 0.611, water contents from 0 to 17 wt % and pressures from 0.0001 to 4.0 GPa) than previous models (Fig. 1), and essentially covers all conditions encountered in crustal settings. In addition, the melt parameter Λ is appropriate for all compositions, including Fe-bearing melts. Each previous model was calibrated for a particular compositional range ( Fig. 1; vertical lines in Fig. 8), over which, that model should be accurate. However, trends in the residual temperatures are observed within the calibrated compositional range of all previous models when applied to the collated dataset (Fig. 8), with the exception of the Shao et al. (2020) model using G as the melt parameter. These trends indicate that the compositional relationships were not well-defined from the limited data used in the previous models. The residuals for Eq. (8) show the least scatter for all variables (Fig. 8) and predicted temperatures for 92% of the collated dataset (626 data) are within 10% of the experimental temperature (Table 3 and Fig. 7), while in comparison, predicted temperatures for < 80% of the collated dataset are within 10% of the experimental temperature using each of the previous zircon saturation models. The average T res is thus smallest for this model (57 °C compared to > 80 °C for all previous models; Table 3). While the propagated uncertainty for Eq. (8) is greater than for the Boehnke et al. (2013) model (96 °C average for total dataset and 48 °C for data < 800 °C compared to 29 °C and 16 °C, respectively), this is not surprising given the additional variables and considerably larger and more diverse dataset used in the construction of Eq. (8).
Collating zircon saturation data over a larger range of conditions allows the relationships with temperature, composition, water content and pressure to be better defined, so extrapolations to lower temperatures beyond the calibration limits of the model are smaller and more accurate. Comparison of the predicted zircon saturation temperatures and the actual experimental temperatures for the collated dataset using each model (Fig. 10), show that in addition to the least scatter in predicted temperatures (92% of samples are predicted within ± 10% of their actual temperature; Table 3), Eq. (8) fits the experimental data < 900 °C well, in contrast to the Boehnke et al. (2013), Gervasoni et al. (2016) and Borisov and Aranovich (2019) models, which underestimate temperatures at these conditions. This likely reflects the limited temperature range of the data included in the Boehnke et al. (2013) and Gervasoni et al. (2016) models, and the emphasis on high-temperature data in the Borisov and Aranovich (2019) model (> 70% is > 1200 °C; Fig. 1a).

Application to natural samples
The Bishop Tuff in California is a zircon-bearing rhyolitic ignimbrite that has been extensively studied. The liquidussolidus temperature interval has been constrained to be between 680 and 750 °C, and the eruption temperature from the early to late units to be between 720 and 800 °C (Fig. 11a) from Fe-Ti oxides (Hildreth and Wilson 2007;Ghiorso and Evans 2008;Jolles and Lange 2019), δ 18 O values (Bindeman and Valley 2002), Ti-in-quartz (Wark et al. 2007; Thomas and Watson 2012) and thermodynamic modelling (Ghiorso and Evans 2008;Gualda et al. 2012;Gualda and Ghiorso 2013;Schiller and Finger 2019).
Zircon saturation temperatures calculated using Eq. (8) for quartz-hosted melt inclusions (Anderson et al. 2000;Schmitt and Simon 2004) are 748 ± 12 °C (n = 16, where n is the number of inclusions), 760 ± 6 °C (n = 6) and 752 ± 15 °C (n = 10), where the uncertainty corresponds to the standard deviation of the n values, for the early, transitional and late eruptive units respectively, which are within the expected range (Fig. 11b). The pressures used in Eq. (8) were those at which the melt would be saturated with the measured CO 2 content and vary between 0.05 and 0.2 GPa (Gualda et al. 2012), over which range the predicted temperature varies by ~ 10 °C. The water contents were determined directly by infrared spectroscopy and SIMS. Zircon saturation temperatures predicted using Eq. (8) from associated pumice clasts (Hildreth and Wilson 2007) and whole-rock samples (Jolles and Lange 2021) are 759 ± 12 °C (n = 17) and 768 ± 25 °C (n = 23), which are also within the expected range.
In contrast, the previous models underestimate the melt inclusion entrapment temperatures (blue circles in Fig. 10 and unfilled circles Fig. 11b), with the exception of the Shao et al. (2020) model using G (average temperature: 764 ± 12 °C). The average predicted temperatures for the melt inclusions are 649 ± 20, 642 ± 28, 678 ± 20, 703 ± 24 and 692 ± 16 °C for the Boehnke et al. (2013), Gervasoni et al. (2016), Borisov and Aranovich (2019), Shao et al. (2019) and Shao et al. (2020) using M models, respectively, with three of these five predictions being below the solidus temperature (Fig. 11b). As discussed previously, these underestimates are attributed to the range of temperature and compositional data used to calibrate the models. The temperature of 712 ± 13 °C predicted by Watson and Harrison (1983) is above the solidus but lower, on average than the estimated eruption temperatures (Fig. 11b). Thus, Eq. (8) improves predictions of zircon saturation temperatures at natural conditions relative to six of the seven previous models, and gives geologically plausible results. The full dataset for this example can be found in Supplementary Table ST4. The zircon saturation model will yield the most accurate results when it is applied to samples where the composition of the melt that was in equilibrium with zircon can be determined. For natural samples, this may be possible using melt inclusions, however, few have been reported that are hosted by zircon, and those preserved in minerals such as quartz may experience post-entrapment re-equilibration (Portnyagin et al. 2008;Moore et al. 2015) or may not have been in equilibrium with zircon (e.g. the hosting mineral crystallised before or after zircon). Unaltered volcanic glasses containing magmatic zircon are potentially ideal samples. However, Hanchar and Watson (2003) noted that almost all the citations of the original Watson and Harrison (1983) model were for application to plutonic rocks rather than volcanic rocks.
For application to plutonic rocks, the whole-rock geochemistry is used as a proxy for the melt. However, this may not be representative of the melt at the time of zircon crystallisation. Figure 12 shows how in a closed tonalitic experimental system the Λ and Zr content of the bulk crystallising component remain constant and are equivalent to the starting melt composition, while Λ and Zr of the residual melt decrease (Marxer and Ulmer 2019). In this example, zircon saturation temperatures predicted using the wholerock composition remain constant (~ 1000 °C i.e. the starting temperature) and are up to 300 °C greater than zircon saturation temperatures predicted using the residual melt from which zircon crystallised. Therefore, zircon saturation temperatures will be overestimated using whole-rock compositions if zircon is a late crystallising phase.
It has been noted that zircon crystallisation temperatures calculated using the Ti-in-zircon thermometer (Watson et al. 2006;Ferry and Watson 2007) are up to 100 °C lower than zircon saturation temperatures (Ickert et al. 2011;Siégel et al. 2018), although this could also be due to overestimates of aTiO 2 (Schiller and Finger 2019;Borisov and Aranovich 2020). Granites are typically intruded at temperatures well above their solidus and have sufficient Zr for zircon to crystallise over a temperature range of ~ 150 °C (Aranovich and Bortnikov 2018;Marxer and Ulmer 2019;Schiller and Finger 2019;Burnham 2020), meaning that zircon saturation will change with the evolution of the melt during fractional crystallisation (Siégel et al. 2018). Therefore, the difference between the calculated zircon saturation based on wholerock Zr contents and Ti-in-zircon temperatures can be used as insight into the timing of zircon crystallisation, and hence whether it crystallised before or after other accessory phases, which is important for interpreting zircon trace element signatures (e.g. Burnham and Berry 2017).
Providing zircon is autocrystic, the zircon saturation model using whole-rock compositions predicts the maximum temperature at which zircon can start to crystallise ( Fig. 12; Schiller and Finger 2019). If zircon is not autocrystic, zircon saturation temperatures of the host may  Shao et al. (2020) using G and M, abbreviated as WH83, B13 G16, BA19, S19, S20a and S20b respectively, and Eq. (8) from this study. Data in Supplementary  Table ST4. In both a and b, the shaded area represents the eruption temperature from the early to late units (Bindeman and Valley 2002) and the dashed lines the approximate liquidus-solidus temperature range bounded by the wet granite solidus at 680 °C (Schiller and Finger 2019). Modified from Schiller and Finger (2019) Author contributions Both authors contributed to the study concept and design. Sample preparation, data collection and analyses were performed by LJC. The first draft of the manuscript was written by LJC, which was revised by both authors.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions.
Availability of data and material (data transparency) All data obtained during this study and collated from previous studies are included in this published article and its supplementary files.

Conflict of interest
The authors have no competing interests to declare that are relevant to the content of this article.

Fig. 12
Melt composition (Λ) versus temperature (°C) for a tonalitic melt (Λ = 0.557; red symbol) undergoing fractional crystallisation at 0.2 GPa (Marxer and Ulmer 2019). As temperature decreases the Λ of the crystallised phases (i.e. the whole-rock composition) remains constant (filled symbols) but Λ decreases for the residual melt (unfilled symbols). The Zr concentrations of the melt at 1000 (also the wholerock), 900 and 800 °C are 1400, 685 and 100 ppm, respectively. Note no Zr was measured in the samples < 800 °C. The modal proportion (wt %) of melt remaining at each temperature is given adjacent to the unfilled symbols 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/.