Divergent T–ƒO2 paths during crystallisation of H2O-rich and H2O-poor magmas as recorded by Ce and U in zircon, with implications for TitaniQ and TitaniZ geothermometry

During solidification of magma chambers as systems closed to chemical exchange with environs, the residual siliceous melt may follow a trend of rising, constant, or decreasing oxidation state, relative to reference buffers such as nickel + nickel oxide (NNO) or fayalite + magnetite + quartz. Titanomagnetite–hemoilmenite thermometry and oxybarometry on quenched volcanic suites yield temperature versus oxygen fugacity arrays of varied positive and negative slopes, the validity of which has been disputed for several years. We resolve the controversy by introducing a new recorder of magmatic redox evolution employing temperature- and redox-sensitive trace-element abundances in zircon. The zircon/melt partition coefficients of cerium and uranium vary oppositely in response to variation of magma redox state, but vary in tandem as temperature varies. Plots of U/Pr versus Ce4+/Ce3+ in zircon provide a robust test for change in oxidation state of the melt during zircon crystallisation from cooling magma, and the plots discriminate thermally induced from redox-induced variation of Ce4+/Ce3+ in zircon. Temperature-dependent lattice strain causes Ce4+/Ce3+ in zircon to increase strongly as zircon crystallises from cooling magma at constant Ce4+/Ce3+ ratio in the melt. We examine 19 zircon populations from igneous complexes in varied tectonic settings. Variation of zircon Ce4+/Ce3+ due to minor variation in melt oxidation state during crystallisation is resolvable in 11 cases but very subordinate to temperature dependence. In many zircon populations described in published literature, there is no resolvable change in redox state of the melt during tenfold variation of Ce4+/Ce3+ in zircons. Varied magmatic redox trends indicated by different slopes on plots of zircon U/Pr versus Ce4+/Ce3+ are corroborated by Fe–Ti-oxide-based T–ƒO2 trends of correspondingly varied slopes. Zircon and Fe–Ti-oxide compositions agree that exceptionally, H2O-rich arc magmas tend to follow a trend of rising oxidation state of the melt during late stages of fluid-saturated magmatic differentiation at upper-crustal pressures. We suggest that H2 and/or SO3 and/or Fe2+ loss from the melt to segregating fluid is largely responsible. Conversely, zircon and Fe–Ti-oxide compositions agree in indicating that H2O-poor magmas tend to follow a T–ƒO2 trend of decreasing oxidation state of the melt during late stages of magmatic differentiation at upper-crustal pressures, because the precipitating mineral assemblage has higher Fe3+/Fe2+ than coexisting rhyolitic melt. We present new evidence showing that the Fe–Ti-oxide oxybarometer calibration by Ghiorso and Evans (Am J Sci 308(9):957–1039, 2008) retrieves experimentally imposed values of ƒO2 in laboratory syntheses of Fe–Ti-oxide pairs to a precision of ± 0.2 log unit, over a large experimental temperature range, without systematic bias up to at least log ƒO2 ≈ NNO + 4.4. Their titanomagnetite–hemoilmenite geothermometer calibration has large systematic errors in application to Ti-poor oxides that precipitate from very oxidised magmas. A key outcome is validation of Fe–Ti-oxide-based values of melt TiO2 activity for use in Ti-in-zircon thermometry and Ti-in-quartz thermobarometry.


Introduction
Accurate geothermometry is a key that often unlocks access to other geologically useful information from chemically analysed glass and/or mineral assemblages-e.g., rock or magma viscosity, pressure of crystallisation and depth of subsequent erosional exhumation, oxygen or sulphur fugacity, or the dissolved H 2 O content in a magma. These pieces of information may be required by earth scientists investigating problems as diverse as rock mechanics, crustal deformation dynamics, ore-forming processes, and volcanic eruption dynamics. Reliable characterisation of the redox state of magmas is important, because about 20% of the elements in the periodic table vary in ionic charge within the range of oxidation states of terrestrial and lunar silicate melts. In many cases, whether an element is sequestered in early minerals or accumulates in residual melts depends on the oxidation state of the melt, which enhances or thwarts a magma's ability to produce ore deposits of many elements. The oxidation state of sulphur in silicate melts strongly influences its availability for degassing during volcanic eruptions, with climatic implications (Scaillet et al. 1998;Burgisser and Scaillet 2007). Evolution of the oxidation state of silicate melts during generation and crystallisation is normally dominated by the Fe 3+ /Fe 2+ couple, with the S 6+ /S 2− couple usually playing a subordinate but occasionally dominant role (Carmichael and Ghiorso 1986), whereas the other redoxsensitive elements are present in natural silicate magmas in relatively low concentrations and respond passively but promptly to the lead by iron and sulphur.
The widely used Ti-in-zircon geothermometer and Ti-inquartz geothermobarometer are sensitive to errors in estimates of melt TiO 2 activity, which are required as input data (Yakymchuk et al. 2017). Different methods of evaluating a melt TiO 2 give substantially discrepant results. As an example, for the Bishop Tuff melt a melt TiO 2 has been estimated at 0.45-0.65 using the Fe-Ti-oxide method (Ghiorso and Gualda 2013), but Thomas and Watson (2012) obtained a melt TiO 2 = 0.15 using the Rhyolite-MELTS algorithm of Gualda et al. (2012). At an actual crystallisation temperature of, say, 750 °C, these discrepancies in a melt TiO 2 correspond to temperature discrepancies of ~ 60 to 100 °C, and many kilobars discrepancy in pressures calculated by the Ti-in-quartz geobarometer (Kularatne and Audétat 2014). Kularatne  may underlie a dispute in the literature (e.g., Ghiorso and Gualda 2013;Ghiorso et al. 2016;Gualda et al. 2017;Evans and Bachmann 2013;Chamberlain et al. 2015;Evans et al. 2016) as to whether the reducing-cooling trend in the Bishop Tuff ( Fig. 1) and similar T-Δ NNO ƒO 2 trends of Yellowstone caldera tuffs, Cougar Point Tuff, Mt. Katmai, and Toba Tuff are a valid representation of T-ƒO 2 conditions in the melt in their pre-eruption magma chambers.
In this paper, we proceed by examining the reliability of the Fe-Ti-oxide geothermometer-oxybarometer calibration that is the basis of the contentious T-ΔΝΝΟ trends, as plotted in Fig. 1. This is followed by introduction of zircon traceelement indicators of temperature variation and redox variation in magmas during crystallisation of zircon populations from many volcanic and intrusive complexes. Some of the zircon populations are from volcanic suites in which Fe-Ti oxides have also been analysed, providing a test for consistency of results. These tests are followed by interpretation of the causes of variation among temperature-redox trends during crystallisation-differentiation of varied magma types, and by suggestions for improved thermodynamic modelling of TiO 2 activity in silicate melts.

Conflicting opinions on validity of varied Fe-Ti-oxide-based T-ƒO 2 trends
Titanomagnetite-hemoilmenite geothermometry and oxybarometry, as thermodynamically calibrated by Anderson et al. (1993) and by Ghiorso and Evans (2008), are widely applied to volcanic rocks with the aim of recovering information about intra-telluric conditions preceding eruption. The validity of the variation among T-ƒO 2 trends obtained from Fe-Ti oxides in Fig. 1 has been questioned by Ghiorso and Gualda (2013), who note that trends having positive slopes are inconsistent with predictions by the MELTS and Rhyolite-MELTS free-energy-minimisation algorithms (Gualda et al. 2012;Ghiorso and Gualda 2013). Rhyolite-MELTS modelling of Bishop Tuff compositions has been used to propose that the crystallisation temperature for Bishop Tuff phenocrysts lies within a 10 °C range from top to bottom of the eruptive sequence (Gualda et al. 2012(Gualda et al. , 2017, whereas Bindeman and Valley (2002), Hildreth and Wilson (2007), Evans and Bachmann (2013), Chamberlain et al. (2015), and Evans et al. (2016) present mutually consistent results from zircon-saturation thermometry, two-feldspar thermometry, Ti-in-quartz thermometry, 18 O magnetite-quartz thermometry and Fe-Ti-oxide thermometry documenting a vertical gradient of ~ 100 °C in the pre-eruptive Bishop Tuff magma chamber tapped by the eruptive sequence. Bindeman and Valley (2002) also report a vertical gradient of ~ 112 °C indicated by 18 O magnetite-quartz thermometry in pre-eruptive magma chamber of the Lava Creek Tuff (Yellowstone), which corroborates the range, as shown in Fig. 1. Ghiorso and Evans (2008) and Ghiorso and Gualda (2013) and Gardner et al. (2014) show plots of temperature versus oxygen fugacity in the Bishop Tuff derived by applying the Ghiorso and Evans (2008) calibration of the Fe-Ti-oxide oxybarometer-thermometer to compositions of Titanomagnetite-hemoilmenite phenocryst pairs reported by Hildreth and Wilson (2007), as shown in Fig. 1. Ghiorso and Gualda (2013) suggest that the discrepancy between the Rhyolite-MELTS model and the Fe-Ti-oxide-based T-ƒO 2 in the Bishop Tuff eruptive suite implies that the Fe-Tioxide phenocrysts (or ilmenite, in particular) were not in equilibrium with the melt that is represented by enclosing volcanic glass. Ghiorso and Gualda (2013) contend that, in that case, the crystal-melt disequilibrium would propagate to erroneous estimates of the thermodynamic activity of TiO 2 in the melt, which would propagate, in turn, to erroneous estimates of magma temperature based on Ti-in-zircon (TitaniZ) geothermometry (Ferry and Watson 2007) and Ti-in-quartz (TitaniQ) geothermobarometry (Wark and Watson 2006;Huang and Audétat 2012;Kularatne and Audétat 2014), both of which require input of TiO 2 activity for temperature calculations.
In contrast, Evans and Bachmann (2013) present several lines of evidence that eruption of the Bishop Tuff magma preserved equilibrium of exchange components and element concentrations among magnetite, ilmenite, biotite, apatite, zircon, and liquid. Evans et al. (2016) and earlier papers cited therein offer numerous lines of independent thermometric and other evidence that the Fe-Ti-oxide phenocrysts probably also maintained equilibrium with feldspars and quartz in the Bishop Tuff magma. Evans et al. (2016) Gardner et al. (1995); Mt St Helens May 1980 and Shiveluch from Blundy et al. (2006Blundy et al. ( , 2008; South Sister Devil's Hills rhyolitic domes, Stelten and Cooper (2012); Cougar Point Tuff, Cathey and Nash (2004); Toba Tuff, Chesner (1998). Data for Daisen, Fish Canyon Tuff, Katmai, Bishop Tuff, and Yellowstone Lake Creek and Huckleberry Ridge tuffs are from the compilation by Ghiorso and Evans (2008), which cites literature sources. The NNO and FMQ buffer curves are from Pownceby et al. (1994) andfrom O'Neill (1987), respectively contend that the T-ƒO 2 trend recovered from the analytical data of Hildreth and Wilson (2007) and thermodynamic calibration by Ghiorso and Evans (2008) are a valid representation of the T-ƒO 2 conditions in the zoned Bishop Tuff magma chamber. In that case, questions arise regarding the adequacy of the calibration data or thermodynamic procedures used by the Rhyolite-MELTS free-energy-minimization algorithm (Gualda et al. 2012;Ghiorso and Gualda 2013) for calculating a melt TiO 2 , with implications for its reliability in thermodynamic modelling of silicic melts more generally. Trail et al. (2015) present a T-ƒO 2 trend for the Bishop Tuff magma, based on X-ray absorption near-edge structure (XANES) measurements of the Ce 4+ /Ce 3+ ratio in zoned zircons and Ti-in-zircon geothermometry in one pumice sample, showing a rise in Ce 4+ /Ce 3+ ratio toward the lower temperature rims of the zircon crystals. They interpret their rim-ward rise of Ce 4+ /Ce 3+ as supporting the arguments by Ghiorso and Gualda (2013) that the cooling-reducing T-ƒO 2 trend indicated by Fe-Ti oxides in the Bishop Tuff is not valid. In that case, cooling-reducing trends indicated by Fe-Ti oxides in rhyolites at other eruptive centres ( Fig. 1), such as Yellowstone caldera tuffs, Cougar Point Tuff, Katmai/Novarupta tuff, and the Toba Tuff are dubious too, as they also violate melt evolution trends derived by the Rhyolite-MELTS algorithm (Ghiorso and Gualda 2013).

Test of the Ghiorso and Evans (2008) Fe-Ti-oxide thermometer-oxybarometer
To resolve these conflicting propositions, we begin by testing the reliability of the Ghiorso and Evans (2008) thermodynamic calibration of the Fe-Ti-oxide oxybarometer by using it to calculate ƒO 2 values from titanomagnetite-hemoilmenite compositions synthesised at controlled ƒO 2 and temperature in seven experimental studies (listed in Electronic Supplementary Material 1). Most of them are represented in a plot of ΔΝΝΟ vs Τ shown as Fig. 27 of Ghiorso and Evans (2008), but we use here different plotting variables, and we include some additional experimental data (Blatter et al. 2013) that were not included in their calibration dataset.
Excluding only the highly deviant results of the study by Prouteau and Scaillet (2003) at 9.2 kbar in extremely hydrous felsic melts, we show in Fig. 2a that the Ghiorso and Evans (2008) oxybarometer retrieves the reported experimental ƒO 2 values of all 45 experiments in the other seven studies with a standard error of just ± 0.2 log unit ƒO 2 , which does not exceed the 0.2-0.3 log unit uncertainty in experimental control of ƒO 2 (Scaillet and Evans 1999;Pichavant et al. 2002). Figure 2b shows large deviations of calculated temperatures from the reported experimental temperatures for Ti-poor oxide pairs synthesised at high ƒO 2 . Under those conditions, experimental calibration data available to Ghiorso and Evans were scant; in a Roozeboom plot of X ilm versus X ulv (Ghiorso and Evans 2008, their Fig. 14),  Ghiorso and Evans (2008) thermodynamic calibration of the Fe-Ti-oxide oxybarometer is tested by using it to calculate ƒO 2 values (at experiment temperatures) from titanomagnetite-hemoilmenite compositions reported in 45 experiments by seven investigative groups over an experimentally controlled (buffered) log ƒO 2 range from NNO to NNO + 4.4, 0.5 to 7 kbar (mostly 2-3 kbar), and 750-1050 °C in hydrous, metaluminous bulk compositions. Colours represent experimental temperature intervals; symbols of different shape identify different investigations. The abscissa variable is the mole fraction of FeTiO 3 in Mg-Fe-Mn-Al-Ti hemoilmenite projected onto the FeTiO 3 -Fe 2 O 3 binary, calculated according to the procedure of Ghiorso and Evans (2008). b There are large systematic deviations of temperatures calculated by the Ghiorso and Evans (2008) geothermometer from Fe-Ti-oxide pairs synthesised at high ƒO 2 (low mole fraction FeTiO 3 ), evident in this plot of the same experiments as in a. But over the composition range of Fe-Ti-oxide pairs in the Bishop Tuff (grey bar ;Hildreth 1977;Hildreth and Wilson 2007), the geothermometer is accurate. The plotted data are listed in the Electronic Supplementary Material 1 isotherms converge tightly, causing poor temperature resolution in Ti-poor compositions. The data trend at low X' ilmenite in Fig. 2b may be used to estimate empirical corrections to calculated temperatures of Ti-poor oxide pairs until an improved thermodynamic model becomes available. However, Fig. 2b also indicates that the Ghiorso and Evans (2008) Fe-Ti-oxide geothermometer is accurate in the composition range of Fe-Ti oxides in the Bishop Tuff and other suites plotted in Fig. 1 (X' ilmenite > 0.45, ΔΝΝΟ ≤ 2). Therefore, we conclude that mis-calibration of the titanomagnetite-hemoilmenite oxybarometer-thermometer is not responsible for the disagreement of the positive T-ƒO 2 trend of the Bishop Tuff Fe-Ti-oxide data in Fig. 1 with the negative slope of the magmatic T-ƒO 2 trend predicted by Rhyolite-MELTS or negative T-ƒO 2 slope inferred by Trail et al. (2015) from zircon Ce 4+ /Ce 3+ measurements and Tiin-zircon temperatures.

Independent test of magmatic T-ƒO 2 trends using Ce and U redox couples
Zircon takes in cerium and uranium in quantities that are easy to measure quantitatively. Zircon is very resistant to adjustment of its initial composition by ion diffusion within the crystal lattice (Cherniak and Watson 2003). Ce and U vary in ionic charge within the ƒO 2 range of terrestrial silicate melts; zircon/melt partition coefficients of Ce and U vary sensitively with ƒO 2 and can, in principle, be used as magmatic oxybarometers. We examine first the variation in ionic charge of Ce and U in the melt as its T and ƒO 2 vary and follow with an examination of the role of temperaturesensitive lattice strain on zircon/melt partition coefficients of U 4+ , U 5+ , Ce 4+ , Ce 3+ , and trivalent lanthanides. That is followed by introduction of plots of U/Pr versus Ce 4+ /Ce 3+ in zircon as robust indicators of T-ƒO 2 trends as zircon crystallises during magmatic differentiation. Those zircon-based magmatic redox trends are to be compared with Fe-Ti-oxide trends in Fig. 1 to test for consistency.

Relative sensitivity of Ce, U, and Fe redox couples in the melt to variation of T and ƒO 2
The effect of varying temperature on Ce 4+ /Ce 3+ and U 5+ / U 4+ in silicate melts is illustrated in Fig. 3. These show that the T-ƒO 2 slopes of contours of Ce 4+ /Ce 3+ and U 5+ /U 4+ in the melt are nearly indistinguishable from the slope of Fe 3+ / Fe 2+ contours in the melt and from the NNO and FMQ reference buffers. Schreiber (1987) and Schreiber et al. (1987) use redox micro-titrations and optical spectroscopy to measure redox ratios in silicate glasses quenched at varied ƒO 2 . They report Ce 4+ /Ce 3+ = 1 at log ƒO 2 = − 0.25 at 1150 °C and at log ƒO 2 = + 3.13 at 1500 °C (uncertainties are not given). They report that U 5+ /U 4+ = 1 at log ƒO 2 = − 8.8 at 1150 °C and log ƒO 2 = − 5.1 at 1500 °C. Those results are in satisfactory accord with ƒO 2 dependence of mineral/melt partition coefficients of Ce and U, as experimentally constrained by Burnham and Berry (2012) and of U by Fonseca et al. (2014), and in accord with range brackets by Farges et al. (1992) on U 5+ /U 4+ measured by optical spectroscopy in aluminosilicate glasses synthesised at varied ƒO 2 .
We can calculate the average enthalpy of the Ce and U redox reactions in the melt over that temperature interval using the Van't Hoff relation (T in kelvins): We obtain ΔΗ = − 466.55 kJ for the homogenous reaction 2 Ce 2 O 3 + O 2 = 4 CeO 2 in the melt. That value compares closely with − 463.44 kJ for the NNO buffer reaction 2Ni + O 2 = 2 NiO (Pownceby and O'Neill (1994). For the uranium oxidation reaction 4 UO 2 + O 2 = 4 UO 2.5 in the melt, we obtain Δ H = − 510.72 kJ. That result is close to ∆ H = − 499.18 kJ for the FMQ buffer reaction, 3 Fe 2 SiO 4 + O 2 = 2 Fe 3 O 4 + 3 SiO 2 in the same temperature interval (O'Neill 1987). Our result for the cerium oxidation reaction is in accord with the value − 436.8 ± 38.4 kJ (2σ) reported by Smythe and Brenan (2015) in experiments at 1200-1500 °C (1)

Fig. 3
Position and slope of the fayalite-magnetite-quartz equilibrium and nickel metal-nickel oxide equilibrium, as calibrated by O'Neill (1987) and Pownceby and O'Neill (1994), are compared with experimentally calibrated contours along which Ce 4+ /Ce 3+ and U 5+ / U 4+ and Fe 3+ /Fe 2+ in the silicate melt are unity. The cerium and uranium redox ratios in synthetic aluminosilicate glasses at 1150 and 1500 °C are from Schreiber (1987) and Schreiber et al. (1987). The curve for Fe 3+ /Fe 2+ is according to the experimental calibration by Kress and Carmichael (1991). Near parallelism of the curves implies very similar enthalpies for these redox reactions and 1 bar ƒΟ 2 on anhydrous, Fe-free, metaluminous haplobasalt, haplo-andesite, and haplo-rhyolite compositions in which Ce 4+ /Ce 3+ was determined potentiometrically and by XANES spectroscopy. Our result for Ce is also compatible with a value of − 397 ± 292 kJ (2σ) reported by Burnham and Berry (2014) in experiments at 1300-1500 °C in melt of anorthite-diopside eutectic composition in which they measured Ce4+/Ce3+ by XANES in the quenched glass. We are unaware of any independent constraints on the enthalpy of the U 4+ → U 5+ reaction. The results calculated above and shown in Fig. 3 indicate that T-ƒO 2 slopes of the Ce and U redox reactions in aluminosilicate melts are statistically indistinguishable from the T-ƒO 2 slopes of the NNO and FMQ reference buffers. We expect that magmatic differentiation trends at essentially constant Ce 4+ /Ce 3+ and U 5+ /U 4+ in the melt are the usual case over most of the course of magmatic differentiation in upper-crustal magma chambers. Carmichael (1991) reports that experimentally constrained T-ƒO 2 slopes of contours of Fe 3+ /Fe 2+ in metaluminous aluminosilicate melts are also indistinguishable from those of the NNO and FMQ buffers, as we illustrate in Fig. 3. As he shows, natural metaluminous magmatic differentiation series that are dominated by fractional crystallisation in convergent and divergent platemargin and intra-plate tectonic settings typically evolve from basaltic-andesitic through dacitic compositions along T-ƒO 2 trends subparallel to NNO or FMQ: • Hildreth (1983) uses Fe-Ti oxide phenocryst pairs to document a T-ƒO 2 trend from andesitic to felsic melts subparallel to the FMQ buffer at Mt. Katmai, Alaska. • Blundy et al. (2006Blundy et al. ( , 2008 demonstrate Fe-Ti-oxidebased T-ƒO 2 trends parallel to NNO in arc andesite-todacite differentiation series at Shiveluch (Kamchatka) and Mt St Helens (Washington). • Rohrlach (2002) and Rohrlach and Loucks (2005) show Fe-Ti-oxide-based T-ƒO 2 trends parallel to NNO in arc andesite-to-dacite differentiation series at Tampakan volcano (Mindanao, Philippines). • Perfit and Fornari (1983a, b) report that along the Eastern Galapagos Rift spreading ridge and Inca Transform, the magmatic differentiation series from olivine tholeiitic melts to Fe-Ti basalt to icelandite and ferro-dacite (glass compositions 51-64 wt% SiO 2 and 7.5 to 1.0 wt% MgO) evolved subparallel to the FMQ buffer along the liquid line of descent from ~ 1200 °C to < 910 °C. • Righter et al. (1998) demonstrate a T-ƒO 2 trajectory parallel to FMQ over the course of magmatic differentiation from olivine basalt to rhyolite at Volcán Alcedo, Galápagos. • Carmichael (1967) demonstrates the same trajectory for the basalt-to-rhyodacite series at Thingmuli, Iceland.
• Delano (2001) examines a large database for mid-oceanridge igneous suites and demonstrates that a T-ƒO 2 trajectory parallel to FMQ is typical.
According to Carmichael (1991), this behavior is because (1) the magma chambers are effectively closed to O 2 exchange with environs and (2) during cooling and crystallisation, the poly-phase cotectic mineral assemblage extracts Fe 3+ and Fe 2+ in essentially the same proportions as were present in the melt at the onset of Fe-Ti-oxide precipitation. To the degree that a magmatic differentiation series evolves along such a trend of constant Fe 3+ /Fe 2+ in the melt, it also evolves along a trend of essentially constant Ce 4+ /Ce 3+ and U 5+ /U 4+ in the melt, as shown in Fig. 3. In a following section, we introduce an independent, robust demonstration that igneous intrusions usually follow a liquid line of descent from basaltic-andesitic through dacitic compositions controlled by fractional crystallisation in a system closed to oxygen and hydrogen. Blundy and Wood (1994) introduced a method for estimating mineral-melt trace-element partition coefficients, based on the elasticity of the crystal lattice and the degree to which substituent ions of various sizes misfit the lattice sites of the host mineral. Their method extends to crystal-melt ionexchange reactions an expression derived by Brice (1975) that relates the mechanical strain energy, ΔG strain , around a size-mis-fitting ion to the elasticity (Young's Modulus, E) of the host crystal:

Temperature-dependent lattice-strain influence on M 4+ and M 3+ partition coefficients
wherein ion radius r i and partition coefficient D i refer to the substituent element i in zircon (z) and melt (m), and r o and D o refer to an element that is an ideal, strain-free fit in the lattice site. N A is Avogadro's number, R is the gas constant, and T is in kelvins. This relation is of the same form as the familiar relation between the Gibbs energy and equilibrium constant for an ion-exchange reaction, ΔG° = − RΤ In Κ exchg , but Nernst-type partition coefficients, D i (in units of concentration by weight) can be used, because conversion factors for ppmw to mole fractions of components in zircon and melt cancel in the ratio. Division of Eq. (2) by RT yields a relation that illustrates the reciprocal variation of partition coefficients with temperature.
We apply Eq.
(2) to zircon/melt partition coefficients reported for a series of tetravalent ions and trivalent ions, to evaluate r o and D o and E by least-squares regression on each valence series. Ionic radii, r i , are listed in Table 1 (from Shannon 1976). Figure 4a illustrates the results for a series of tetravalent and trivalent ions in the Zr-dominated sublattice of zircon, as measured by ion probe on zircons and haplo-andesitic melt in experiments at 1300 °C and 1 atm and log ƒO 2 ≤ FMQ by Burnham and Berry (2012). Figure 4a also shows a series of partition coefficients, measured by ion probe on unpolished surfaces of natural zircon crystals and on enclosing rhyolitic volcanic glass, reported by Claiborne et al. (2018). Figure 4a shows that parabolas fitted to the data are tighter for tetravalent ions than for trivalent ions at each temperature, and the parabolas are narrower at lower temperature, reflecting decrease in lattice-site elasticity with increasing ionic charge and with decreasing temperature.
Least-squares fits of Eq.
(2) to those and other volcanic samples described by Claiborne et al. (2018) yield the parameters listed in Table 2. The derived values in Table 2  to increase with decreasing temperature is suggested.
In the case of mis-fitting trivalent ions, there is a competition between the tendency of crystal/melt partition coefficients to increase with decreasing temperature, and the increasing inability of the lattice to deform to accommodate them as temperature decreases. Figure 4b shows a plot of zircon/melt partition coefficients of Ce 4+ and Ce 3+ (Table 2) at the experimental temperature 1300 °C or at the zircon-saturation temperature reported by Claiborne et al. for natural zircon-glass pairs. Dispersion of the data for volcanic suites is due mainly to uncertainty in experimental calibration of the zircon-saturation geothermometer (around ± 15 °C, 1σ; Hanchar and Watson 2003), and due to micrometer-scale compositional zoning in analysis spots 15 µm diameter on "rims" of natural crystals in samples IHB, ITHn, and HRL21. Natural surfaces of zircon crystals were analysed in SHL samples. Magma chamber replenishments and mixing preceding and accompanying volcanic eruptions may, in some cases, cause mis-match of zircon crystals with the melt from which they grew. Least-squares regression of the array in In Fig. 4a, the measured values of D z/m U plot below the parabolas fitted to other tetravalent ions. Blundy and Wood (2003), Ballard et al. (2002), and others have noted the same sense of deviation of measured of D z/m U from the value predicted for U 4+ in natural zircons from peralkaline and metaluminous rocks. The deviation may be due to incorporation of a minor proportion of the uranium as U 5+ in zircon via the charge-coupled substitution 2 Zr 4+ ↔ REE 3+ + U 5+ . Burnham and Berry (2017) report that zircons from peraluminous S-type granites have phosphorus contents near the amounts expected for charge balance of the measured abundances of Y + REE 3+ according to the charge-coupled substitution Zr 4+ + Si 4+ ↔ REE 3+ + P 5+ , but they report that zircons from metaluminous I-type granitoids have a systematic deficit of charge-balancing phosphorus. This is understandable in light of experiments by Pichavant et al. (1992) showing that, at a given temperature and pressure and wt% H 2 O in the melt, apatite solubility and phosphorus content of silicate melts increase dramatically with increasing Al/(Ca + Na + K), due Table 1 Effective ionic radii (Å) of cations in (A) eightfold coordination with oxygen that may substitute for Zr in zircon (Shannon 1976) and (B) fourfold coordination with oxygen that may substitute for Si in zircon (Shannon 1976 to increase of Ca-Al complexing in the melt with increasing A/CNK. Therefore, in metaluminous melts with higher lime activity and lower phosphate activity at apatite saturation, other charge-balancing mechanisms are implicated for incorporation of REE 3+ into zircon. Ionic radii in Table 1 suggest that U 5+ should be a lower strain candidate for that charge-coupled substitution role than, say, Nb 5+ , Ta 5+ , or V 5+ . Variation of the charge-coupled substitution mechanism with melt A/CNK is expected to diminish a correlation of E 3+ with temperature in Table 2.

Cooling-induced variation of Ce 4+ /Ce 3+ and U 5+ /U 4+ in zircon at constant ratios in the melt
In Fig. 5, we present new measurements illustrating strongly temperature-dependent variation of Ce 4+ /Ce 3+ in a natural zircon assemblage from a calc-alkalic volcanic arc suite. A few drill-core and outcrop samples of eruptive andesites and dacites from the Tampakan volcanic complex are sufficiently free of hydrothermal alteration to permit application of the titanomagnetite-ilmenite geothermometer (calibration of Ghiorso and Evans 2008) and the plagioclase-hornblende geothermometer (calibration of Holland and Blundy 1994). New data for Tampakan phenocryst compositions that are  rison 1983). Red symbols represent ion-probe analyses of zircon rims and enclosing haplo-andesitic glass produced in experiments at 1300 °C and 1 atm by Burnham and Berry (2012). Ionic radii are from Shannon (1976). b Temperature dependence of Ce 4+ /Ce 3+ fractiona-tion between zircon and melt due to lattice strain, as represented by variation of the zircon/melt partition coefficients of Ce 4+ and Ce 3+ with temperature, determined by least-squares fitting (as in a) of Eq.
(2) to partition coefficients and zircon-saturation temperatures reported for natural phenocryst-glass pairs by Claiborne et al. (2018) (black and blue symbols; standard errors not reported) and fitted to experimental (1300 °C) zircon/melt partition coefficients analysed by SIMS and by LA-ICPMS by Burnham and Berry (2012 Wood (1994, 2003). All available paired T-Ce 4+ /Ce 3+ data are plotted in Fig. 5, which shows strong lattice-strain discrimination, as the zircon lattice thermally contracts during growth from cooling magma, in favour of smaller, better-fitting Ce 4+ relative to larger, size-and charge-mis-fitting Ce 3+ ions in Zr 4+ lattice sites. The T-ƒO 2 trend of Tampakan Fe-Ti oxide data in Fig. 1 and trace-element compositions of Tampakan zircons presented ahead (Fig. 9) demonstrate that the Tampakan andesite-dacite series developed at essentially constant oxidation state of Fe, Ce, and U in the melt. The effect of thermally modulated lattice strain on incorporation of U 4+ , U 5+ , Ce 4+ , Ce 3+ , and other trivalent lanthanides is illustrated in Fig. 6, wherein each panel shows a compilation of experimental results at varying temperature but constant pressure, constant bulk composition of the charge, and constant oxidation state of Ce and U in the melt. In each suite of experiments, those at the three highest temperatures are illustrated, as they are most likely to approximate equilibrium partitioning. Each panel shows experimentally determined zircon/melt partition coefficients of an ion, U 4+ or Ce 4+ , that fits relatively well in the Zr 4+ sublattice of zircon, and also shows partition coefficients of a larger ion that fits relatively poorly in the Zr 4+ sublattice (see ionic radii in Table 1). In each panel of Fig. 6, the ratio of the better-fitting, smaller ion to poorer fitting, larger ion increases in zircon with decreasing temperature at constant redox state of the melt-including the ratio Ce 4+ /Ce 3+ in zircon (Fig. 6d). All experiments were buffered at NNO, which is shown in Fig. 3 to correspond to effectively constant oxidation state of U and Ce in the melt. These experiments corroborate calculated results in Fig. 4 for temperature dependence of lattice-strain influence on zircon/melt partition coefficients of ions of various sizes.
In thoughtfully designed experimental studies aiming to calibrate Ce 4+ /Ce 3+ in zircon as a magmatic oxybarometer, Trail et al. (2011Trail et al. ( , 2012 recognise substantial temperature dependence of zircon/melt partition coefficients of Ce 3+ and Ce 4+ . However, zircon/melt disequilibrium partitioning is problematic in laboratory experiments (and in nature), as discussed by Luo andAyers (2009), Hofmann et al. (2014), and Trail et al. (2015), so the temperature dependence has not been well characterised by laboratory experiments on granitoid magmas. However, Claiborne et al. (2018) address much of that deficiency by determination of zircon/glass partition coefficients in natural volcanic suites, coupled with zircon-saturation geothermometry to characterise the varying temperature dependence of zircon/melt partition coefficients for ions of charge 5 +, 4 +, and 3 +. Figure 7 illustrates the effect of varying magmatic ƒO 2 on zircon/melt partition coefficients of Ce and U at constant pressure (1 atm) and temperature (1300 °C) and constant elemental composition of the melt (Fe-free haplo-andesite). The experiments show that zircon/melt partition coefficients of U and Ce have opposite senses of variation in response to isothermally varying ƒO 2 . Zircon/melt partition coefficients of the charge-invariant elements La, Pr, Nd, Sm, Gd and heavier REEs do not vary over that range of ƒO 2 . Therefore, within the log ƒO 2 range FMQ−2 to FMQ+5 that is relevant to most terrestrial magmas, nearly all Ce in zircon is Ce 4+ , and the element ratios U/Pr and Ce/Nd show opposite variation trends. Ce/Nd is a proxy for Ce 4+ /Ce 3+ in these experiments, in which Ce valence was not measured (Fig. 7). Although the bulk composition of experimental charges includes an ample supply of Sc 3+ , Y 3+ and REE 3+ and of P 5+ , Nb 5+ and Ta 5+ ions that could team up for coupled substitutions of the type 2 Zr 4+ ↔ M 3+ + M 5+ or Zr 4+ + Si 4+ ↔ M 3+ + M 5+ , the trends of zircon/melt partition coefficients of Ce and U demonstrate that coupled substitutions of For zircon-bearing silicic andesites and dacites from the Tampakan volcano, Mindanao, that are sufficiently free of hydrothermal alteration to permit geothermometry on associated phenocrysts, we combine whole-rock ICPMS analyses of trace elements with analyses of zircons by ArF ELA-ICPMS (Rohrlach 2002) to calculate zircon Ce 4+ /Ce 3+ ratios by the lattice-strain method of Blundy and Wood (1994). Those Ce 4+ /Ce 3+ ratios are plotted against temperatures obtained for the same samples on touching titanomagnetite-ilmenite phenocryst pairs (calibration of Ghiorso and Evans 2008) and on hornblende-plagioclase phenocryst pairs (calibration by Holland and Blundy 1994). Estimated uncertainty bars are ± 1 SE charge-and size-misfitting ions are strongly disfavoured by zircon. The variation trends can be described by a sigmoid function, as derived by Burnham and Berry (2012). Burnham and Berry (2012) question the validity of the U trend in Fig. 7 because of suspected (but unsubstantiated) volatilisation of UO 3 from the melt in the 1-atm gas-flow furnace used for these experiments. However, Luo and Ayers (2009)  That 10-fold range in D zirc/melt U is the same as in the equivalent ∆FMQ range in Fig. 7. Various spectroscopic methods have been employed by Schreiber and Andrews (1980), Farges et al. (1992) and Berry et al. (2008) to demonstrate that the oxidation of U 4+ to U 6+ in the melt proceeds via a U 5+ intermediate stage which has a predominance range a bit above FMQ at magmatic temperatures.

Introduction of zircon U/Pr versus Ce 4+ /Ce 3+ plots
The experiments reported in Fig. 7 demonstrate that, at oxidation states of common terrestrial magmas, essentially all Ce in zircon occurs as Ce 4+ and nearly all U as U 4+ . We introduce plots of U/Pr versus Ce 4+ /Ce 3+ in zircon to discriminate whether Ce 4+ /Ce 3+ variations in zircon are due to redox variations in the melt during zircon crystallisation and/ or due to temperature-regulated fractionation of Ce 4+ from ) of U, Ce, La, Pr, Gd, and Dy in experiments in which ƒO 2 was buffered by Ni + NiO, at constant pressure, and nearly constant major-element composition of the hydrous granitoid melts. In each panel, the element that is most compatible in zircon is red; the element least compatible in zircon is blue. In each panel the ratio of the better-fitting ion to the poorer fitting ion rises with decreasing temperature. a, b Experimen-tal results by Rubatto and Hermann (2007) in hydrous peraluminous tonalite, 66.1-69.8 wt% SiO 2 at 20 kbar and ƒO 2 buffered at NNO. c Experimental results by Luo and Ayers (2009)  Ce 3+ by lattice strain. Table 1 shows that the ionic radius of U 4+ in eightfold coordination with oxygen is 1.00 Å, very similar to that of Ce 4+ , 0.97 Å, and shows that the U 4+ /Pr 3+ ionic radius ratio, 0.888, is very similar to the Ce 4+ /Ce 3+ ionic radius ratio, 0.849. Therefore, the effect of temperature variation on the relative zircon/melt partition coefficients of these ratios must be nearly identical. The close similarity of ionic radii of Pr 3+ and Ce 3+ is the rationale for choosing Pr as a denominator element. However, those ratios in zircon vary oppositely in response to ƒO 2 variation, as indicated in Fig. 7b. This test of T-ƒO 2 trends during magmatic differentiation is applicable to slowly cooled intrusive rocks, and it can test the reliability of T-ƒO 2 trends recovered from titanomagnetite-hemoilmenite phenocryst pairs in volcanic rocks. The Fe-Ti-oxide oxybarometer-thermometer is reliably applicable only to quickly quenched volcanic rocks, inasmuch as Fe-Ti oxides re-equilibrate on the time scale of days to weeks at near-solidus temperatures (Venezky and Rutherford 1997).

Analytical and computational methods and sources of data dispersion
Electronic Supplementary Material 3 lists previously unpublished analyses of zircons from the Tampakan volcanic complex, Mindanao, done by ArF ELA-ICPMS at the Research School of Earth Sciences, Australian National University. In following paragraphs, we show zircon U/Pr versus Ce 4+ / Ce 3+ derived from ion-probe and LA-ICPMS analyses that were compiled from the literature and combined with our data for the Tampakan complex. Calculation of Ce 4+ /Ce 3+ by the lattice-strain method of Wood (1994, 2003) utilises zircon/melt partition coefficients obtained by combining analyses of zircon crystal rims or surfaces with analyses of associated glass (in volcanic suites) or with analyses of host whole rocks (in intrusive suites).
Zircons in all suites in the following data plots have complex oscillatory growth bands and sector zoning similar to that of Bishop Tuff zircon in Fig. 8. There is evidence from synthetic and natural zircons that very slow diffusion of high-electric-field-strength ions of charge 3+, 4 + and 5 + through a microns-thick melt boundary layer adjacent to a growing zircon crystal surface results in disequilibrium between the boundary-layer melt film and the bulk melt, and hence between zircon and bulk melt; that diffusional disequilibrium causes oscillatory compositional zoning of zircon at sub-micron to tens-of-microns scale (Luo and Ayers 2009). Such disequilibrium oscillations in composition are likely to be especially pronounced in zircons that grew from very viscous, H 2 O-poor, high-silica rhyolitic melts, as documented in the Younger Toba Tuff and Bishop Tuff by Hofmann et al. (2014). There is also persuasive evidence that the zircon/ melt partition coefficients of high-field-strength elements vary between prism and various pyramidal crystal faces along the same growth band in zircon, as demonstrated by Fig. 7 Nernst partition coefficients of Ce and U and other trace elements between zircon and synthetic andesitic melt were measured by secondary-ion mass spectrometry (SIMS, ion probe) and by LA-ICPMS in Burnham and Berry's (2012) experiments at 1300 °C and 14 log units variation of oxygen fugacity relative to the fayalite-magnetite-quartz reference buffer (FMQ). a Oxidation state of cerium varies from entirely Ce 3+ in melt at the lowest ƒO 2 to increased (but still small) amounts of Ce 4+ in the melt at the highest ƒO 2 , over which range the Ce content of zircon rises ~ 100-fold as the proportion of Ce 4+ in the melt increases. Uranium shows the opposite trend, being entirely U 4+ in the melt at the lowest ƒO 2 , but with increasing amounts of U 5+ the melt as ƒO 2 increases. The data show a very strong preference for U 4+ in Zr 4+ lattice sites, and the zircon/melt partition coefficient of U declines approximately 20-fold as the melt's supply of U 4+ diminishes with rising ƒO 2 . Burnham and Berry report standard errors of D z/m Ce that average ± 74% of the measured value; average standard errors of D z/m U are ± 34% of the measured values. b Ce/Nd ratio in zircon (a proxy for ΣCe/Ce 3+ in this experimental dataset in which Ce valence was not reported) rises approximately 100-fold over the experimental ƒO 2 range, while U/Pr in zircon declines by a factor of about 20 Tailby et al. (2011) and in Bishop Tuff zircons by Chamberlain et al. (2014). Although all utilised zircon analyses in Figs. 9, 10, and 11 are described by the respective data sources as "rims", the compositional heterogeneity within outer growth bands can be substantial at the 15-30 μm scale of ion-probe and LA-ICPMS beam spots. This contributes scatter to the zircon composition arrays that exceeds instrumental analytical errors on homogeneous material. However, the data dispersion within zircon arrays typically does not substantially exceed the dispersion within Fe-Ti-oxide arrays for volcanic suites in Fig. 1. Dispersion of zircon and Fe-Ti-oxide compositions within volcanic suites is in many cases due in part to chamber replenishments and magma mixing immediately preceding the eruptions, as documented for the Bishop Tuff, Younger Toba Tuff, and Yellowstone supereruptions, among others (Evans and Bachmann 2013;Chesner and Luhr 2010;Wotzlaw et al. 2015). Eruption of inhomogeneous melt contributes some uncertainty in identifying equilibrium pairs of zircon rims and host glass in those volcanic suites.

Suites showing no change in oxidation state during cooling
Electronic Supplementary Material 3 lists Tampakan zircon trace-element analyses by ArF ELA-ICPMS at the Research School of Earth Sciences, Australian National University. Our analyses of Tampakan zircons have been combined with our compilation from published literature of zircon compositions in all igneous suites for which Ce 4+ /Ce 3+ in zircon was calculated by the respective authors using the lattice-strain method, and which have enough data points to perform a test of trend. No examples were found that show a negative correlation of U/Pr with Ce 4+ /Ce 3+ , as would represent variation of Ce 4+ /Ce 3+ dominated by changing oxidation state of the melt. Figure 9 shows strong positive correlations of zircon Ce 4+ /Ce 3+ with U/Pr in eight igneous complexes. In each of them, the regression slope is near 1.0, which means that the variations in both ratios are attributable entirely to cooling during crystallisation. Accordingly, there is no resolvable variation in redox state of the melt (as represented by redox indicators in Fig. 3) within any zircon suite shown in Fig. 9, including the Tampakan suite that is also illustrated in Fig. 5. Of the igneous suites represented in Fig. 9, the Tampakan volcanic suite is the only one for which analyses of Fe-Ti-oxide phenocrysts are available for independent corroboration of the T-ƒO 2 trend (Electronic Supplementary Material 2). The Tampakan Fe-Ti oxide T-ƒO 2 array in Fig. 1 shows no resolvable change in oxidation state of the melt over the 920-705 °C temperature interval, which corroborates the redox independence of the Tampakan zircon trend in Fig. 9 showing ~ 1:1 co-variation of U/Pr with Ce 4+ /Ce 3+ .

Suites showing oxidation during cooling
In Fig. 10, we show a collection of zircon analyses from rhyolitic flow domes at South Sister volcano (Cascade arc) and from the 3.5 ka Yn tephra erupted from Mt St Helens (Cascade arc) and from six subduction-related calc-alkalic intrusions. All except South Sister and Montecristo are "adakitic" (high Sr/Y and La/Yb, modest or no Eu anomaly in the whole-rock REE pattern). Montecristo is a pre-ore intrusion in the Chuquicamata district. The chondrite-normalised REE patterns of the South Sister's Devil's Hills rhyolite lava domes are "adakite-like", steeply dipping, with a minimum at Ho and La N /Yb N = 14 and Eu/Eu* = 0.96 to 0.97 in the two rock analyses reported by Stelten and Cooper (2012).
All the Ce 4+ /Ce 3+ values plotted in Fig. 10 are those reported in the cited publications (following the calculation  Hofmann et al (2014) and Chamberlain et al. (2014) to correspond to significant variations of traceelement composition from band to band within a sector, and from sector to sector along individual growth bands. The outer oscillatory growth bands are narrower than the usual 10-30 µm diameter of ionprobe or LA-ICPMS beam spot analyses on zircon "rims", so such heterogeneity contributes dispersion to arrays of precise analyses in Figs. 9, 10, and 11 procedure of Ballard et al. 2002, using whole-rock compositions), except South Sister and Mt St Helens Yn tephra, which we calculated using glass analyses reported by the respective authors (Electronic Supplementary Material 4). All data arrays in Fig. 10 have regression slopes significantly < 1. As indicated by the inset vector-decompositions of the data trends, we interpret these regression trends as signifying magma oxidation during crystallisation of the analysed zircon suites. The Devil's Hills rhyolitic lava domes in the South Sister volcanic complex and the Mt St Helens Yn tephra are the only suites in Fig. 10 for which Fe-Ti oxide data are available for an independent test of redox trend. The negative slopes of the Mt St Helens Yn tephra and South Sister Fe-Ti oxide array in Fig. 1 corroborate the zircon-based U/Pr versus Ce 4+ /Ce 3+ trends in Fig. 10 that indicates melt oxidation during cooling and crystallisation. Prouteau and Scaillet (2003), Rohrlach and Loucks (2005), Loucks (2014) and Lu et al (2015), among others, Fig. 9 Ce 4+ /Ce 3+ ratios in zircons (calculated by the lattice-strain method of Wood 1994, 2003) from eight igneous units are plotted against U/Pr of the zircons, measured by ArF ELA-ICPMS. As discussed in the text, these ratios must have nearly identical temperature sensitivity of their zircon/melt partition coefficients. Note that the slopes of the arrays do not differ significantly from 1.0 in these eight examples. Standard error of the slope and intercept are ± 1σ. After considering the presence of non-negligible (but unquanti-fiable) analytical uncertainty in these datasets, it is apparent from the slopes and R 2 values that there is no resolvable component of the variance in Ce 4+ /Ce 3+ within each array that is attributable to variation of magmatic oxidation state. Data sources: Koksai, Erdenet, Nurkazghan, and Bozshakol porphyry copper deposits: Shen et al. (2015); Shiyaogou porphyry molybdenum deposit: Han et al. (2013); El Abra porphyry copper igneous complex: Ballard et al. (2002); Tampakan volcanic complex: Electronic Supplementary Material 3 demonstrate that adakitic magmas are exceptionally hydrous, commonly with 8-14 wt% H 2 O dissolved in the silicate melt over most of the course of high-pressure magmatic differentiation. As silicate melts cannot dissolve that much H 2 O at upper-crustal pressures, it is implied that fluid exsolved during ascent of andesitic and dacitic melts to subvolcanic magma chambers, with crystallisation at fluid-saturated conditions upon epizonal emplacement. Oxidation of iron and other redox couples in the melt accompanies dehydration of the melt (Czamanske and Wones 1973;Candela 1986;Gaillard et al. 2002;Rowe et al. 2006). Figure 1 shows that the adakitic suites (Pinatubo, Daisen, Fish Canyon Tuff, and Mt St Helens Yn tephra) have systematically higher ƒO 2 and negative T-ƒO 2 slopes in the dacitic to rhyolitic composition range. We suggest that exsolution of hydrothermal fluid during magma ascent and epizonal crystallisation is accompanied by partitioning of oxidised sulphur from melt into hydrothermal fluid according to the heterogeneous reaction: wherein exsolution of one mole of S 6+ in the melt as S 4+ in the fluid is accompanied by oxidation of two moles of Fe 2+ in the melt. The effect on Fe 3+ /Fe 2+ ratio in the residual melt would be greatest in oxidised, Fe-poor felsic melts in which the S:Fe ratio is relatively large at the onset of fluid saturation. Scaillet et al. (1998) demonstrate experimentally that oxidised sulphur is much more volatile as SO 2 than as SO 3 , and that at log ƒO 2 > NNO + 1.3, the fluid/melt partition coefficient of sulphur is ≥ 1000 in dacitic and more silicic magmas, so this provides a plausible explanation of the negative T-ƒO 2 slopes at high ƒO 2 in Fig. 1 and slopes < 1 in Fig. 10. Experiments by Scaillet et al. (1998) indicate that the fluid/melt partition coefficient of sulphur is about 300 times smaller at log ƒO 2 ≤ NNO + 0.5, in the stability field of S 2− predominance in the melt. Therefore, gravitational segregation of buoyant hydrothermal fluid as it exsolves should induce some oxidation of the Ce 4+ /Ce 3+ and U 5+ /U 4+ and Fe 3+ /Fe 2+ redox couples in the residual silicate melt and in zircons that grow during decompression and loss of SO 3 to exsolving fluid. Baker and Rutherford (1996) find in experiments at controlled ƒO 2 that change of melt structure with variation of dissolved H 2 O content of metaluminous rhyolite has no effect on the Fe 3+ /Fe 2+ ratio in the melt at temperatures below 850 °C. However, at magmatic temperatures molecular H 2 O in the melt and in exsolved fluid partially dissociates according to the reaction H 2 O → H 2 + ½ O 2 , and loss of H 2 to segregating fluid oxidises iron in the residual melt by the reaction 2FeO + H 2 O → Fe 2 O 3 + H 2 . As pointed out by Carmichael (1991), in relatively H 2 O-poor silicic melts, such as the Toba Tuffs (Chesner and Luhr 2010) and Bishop Tuff (Wallace et al. 1999;Anderson et al. 2000;Roberge et al. 2013) and Yellowstone rhyolites (Almeev et al. 2012;Befus and Gardner 2016), most of the dissolved H 2 O occurs as hydroxyl. Therefore, the H 2 O molecular dissociation reaction does not proceed to as significant a degree as in extremely hydrous magmas.
H 2 transfer from silicate melt to exsolving hydrothermal fluid and gravitational segregation of buoyant bubbles from the conjugate melt (forming a fluid cap at the top of the magma column) may provide a supplementary mechanism for inducing the oxidising-cooling trend revealed by the zircon data in Fig. 10 and by the Fe-Ti-oxide T-ƒO 2 trends ( Fig. 1) of the low-temperature (< 800 °C) dacitic/rhyolitic segments of the Pinatubo and St Helens Yn arrays and similar cooling-oxidising trends of other adakitic suites (e.g., Daisen and Fish Canyon). Czamanske and Wones (1973) attribute to segregation of exsolving magmatic hydrothermal fluid and H 2 loss from the melt an oxidising-cooling trend during crystallisation of the Finnmarka granitoid complex, Norway. Cornejo and Mahood (1997) attribute to H 2 loss the oxidising-cooling T-ƒO 2 trend in the adakitic La Gloria pluton of Miocene age in central Chile. Bell and Simon (2011) have identified an additional mechanism for increasing the Fe 3+ /Fe 2+ ratio in the residual melt as hydrothermal fluid exsolves and physically segregates. Their experiments and thermodynamic data and modelling studies show that over a large range of ƒO 2 relevant to natural magmas, Fe 2+ partitions into a chloride-bearing magmatic hydrothermal fluid much more strongly than Fe 3+ , with the result that fractional segregation of the brine raises the oxidation state of the residual melt. The effect is especially pronounced during exsolution of a relatively large mass of brine from relatively Fe-poor felsic melts.

Suites showing reduction during cooling
In contrast to the differentiation trends of wet magmas that have slopes < 1 in Fig. 10 and negative slopes in Fig. 1, we   Fig. 10 Slopes of zircon composition arrays test whether the melt was undergoing reduction or oxidation as the zircon suite crystallised. Zircon Ce 4+ /Ce 3+ values reported by authors of the source publications (evaluated by the lattice-strain method of Blundy and Wood 1994;Ballard et al. 2002) are plotted against the U/Pr ratios in those zircons, except for South Sister and Mt St Helens suites, for which we calculated zircon Ce 4+ /Ce 3+ by the same procedure. All arrays have regression slopes < 1, which implies oxidation of the melt during zircon crystallisation, as indicated in three panels on the right, in which the data trends are decomposed into vector components representing cooling (slope + 1.0) and oxidation (slope − 1.0). Standard error of the slope and intercept are ± 1σ. Similar vector diagrams apply to panels on the left. South Sister data are from Stelten and Cooper (2012); Mt St Helens data from Claiborne (2011) and Claiborne et al. (2018); El Teniente porphyry copper igneous complex (Chile) data from Muñoz et al. (2012); Lower Yangtze Valley data from Wang et al. (2013); Chuquicamata porphyry copper complex and Montecristo (Chile) data are from Ballard et al. (2002) ◂ show in Fig. 11 log U/Pr versus log Ce 4+ /Ce 3+ trends for zircons in H 2 O-poor, fayalite-and pyroxene-phyric, high-silica rhyolites that have positive slopes of the T-ƒO 2 trends in Fig. 1. The data and procedures for evaluating zircon Ce 4+ / Ce 3+ in the suites plotted in Fig. 11 are shown in Electronic Supplementary Material 4. Those trends having positive slopes are a subject of dispute, as described in the introduction of this paper. In Fig. 1, the T-ƒO 2 trend of the Yellowstone Lake Creek fayalite rhyolite and Toba Tuff pyroxene dacite and Toba Tuff high-silica fayalite rhyolite show positive T-ƒO 2 slopes like the Bishop Tuff pyroxene rhyolite. The zircon composition trends in all three panels of Fig. 11 have slopes substantially > 1, which corroborate the cooling-reducing trend from Fe-Ti oxides for those three suites shown in Fig. 1.
In each of these volcanic suites in Fig. 11, and in the Tampakan and South Sister and Mt St Helens Yn volcanic suites in Figs. 9 and 10, the data dispersion in each array is similar to the degree of dispersion within T-ƒO 2 arrays of Fe-Ti oxide phenocrysts from the same volcanic suites in Fig. 1. The Bishop Tuff and Yellowstone zircon Ce 4+ /Ce 3+ values in Fig. 11 were calculated by pairing the average zircon "rim" composition (outer 10-15 µm) in each pumice sample with glass analyses from the same eruptive unit and sampling area, reported by other analysts, to derive apparent zircon/ melt partition coefficients of REEs, Th, Zr and Hf, which were then used to estimate Ce 4+ /Ce 3+ in zircon according to the lattice-strain method of Wood (1994, 2003) (see Electronic Supplementary Material 4). The assumption that zircon and glass represent equilibrium pairs may not be perfect in every case, because chamber replenishment and magma mixing that immediately preceded the eruptions have been well documented in the Bishop Tuff (Evans and Bachmann 2013), Toba Tuff (Chesner and Luhr 2010), and Fig. 11 Slopes of zircon composition arrays test whether the melt was undergoing reduction or oxidation as the zircon suite crystallised. The least-squares regression in each panel has a slope > 1.0, which is decomposed into vector components representing cooling (slope + 1) and decrease in oxidation state of the melt (slope − 1). Standard error of the slope and intercept are ± 1σ. In each panel the data trend shows that the redox state of the melt decreased as zircon crystallised from cooling melt. The trends of these zircon arrays corroborate the Fe-Ti-oxide T-ƒO 2 trends for these volcanic suites in Fig. 1. a Bishop Tuff zircon crystal rims analysed by Smythe and Brenan (2016, sample BTZ) and by Reid et al. (2011;other sample numbers). The label following the comma is the eruptive stratigraphic unit from which the sample was collected (nomenclature of Hildreth and Wilson 2007). Calculation of zircon Ce 4+ /Ce 3+ utilises analyses of glass inclusions in Bishop Tuff phenocrysts by Wallace et al. (1999), Anderson et al. (2000), and Roberge et al. (2013). b Zircon crystal-rim compositions in eruptive Unit B of the Yellowstone caldera-forming Lava Creek Tuff supereruption (0.62 Ma), as reported by Wotzlaw et al. (2015) are combined with Lava Creek Tuff Unit B pumice composition as reported by Christiansen (2001). c Toba Tuff zircons and host glass analyses from two pumice samples reported by Gaither (2011) ▸ Yellowstone Lake Creek Tuff (Wotzlaw et al. 2015). Such possible heterogeneity in composition of the host glass may contribute analytical dispersion in addition to the compositional heterogeneity within zoned crystals (Fig. 8).
The Bishop Tuff data trend in Fig. 11 implies that successive eruptive units tapped successively deeper levels of a magma chamber that was zoned in oxidation state, with early-erupted, more reduced, cooler magma overlying hotter, less reduced magma that was tapped later in the eruption sequence, corroborating interpretations by Hildreth and Wilson (2007) based on independent evidence from Fe-Ti oxides (Fig. 1) and other phenocrysts. Their interpretation that the Bishop Tuff magma chamber had strong vertical gradients in temperature and oxidation state was disputed by Gualda et al. (2012Gualda et al. ( , 2017 who claimed, on the basis of Rhyolite-MELTS modelling, that the crystallisation temperature for Bishop Tuff phenocrysts lies within a 10 °C range from top to bottom of the eruptive sequence. The cause of the reducing trend during cooling of the Bishop Tuff magma is illuminated by consideration of Fe 3+ / Fe 2+ ratios in crystals and coexisting melt. Hildreth and Wilson (2007) report that ilmenite phenocrysts are exceedingly sparse in all eruptive units of the Bishop Tuff, and very subordinate to titanomagnetite in modal abundance. Hildreth (1977) carries out wet chemical measurements of Fe 3+ /Fe 2+ ratios in biotite in several eruptive units of the Bishop Tuff and reports that 40-80% of the iron in biotite is Fe 3+ . Multiple lines of evidence summarized by Evans et al. (2016) and references therein show that the pyroxene phenocrysts found in late-erupted units that tapped deeper, hotter parts of the magma chamber are xenocrysts out of chemical equilibrium with all other phenocrysts and glass in the samples. Therefore, the only modally significant Febearing phenocrysts precipitating from the melts represented by erupted pumice are titanomagnetite with Fe 3+ /Fe 2+ ≈ 1.2 on average (Hildreth and Wilson 2007) and biotite with a similar Fe 3+ /Fe 2+ ratio. At 800-850 °C, the Bishop melt was at log ƒO 2 ≈ NNO + 0.5 (Fig. 1), so Fe 3+ /Fe 2+ in the melt was ~ 0.06 to 0.07 (Gaillard et al. 2001;Wilke et al. 2002). Precipitation of a phenocryst assemblage having a higher Fe 3+ /Fe 2+ than the high-silica, Fe-Mg-poor melt can account for the evolution of the melt toward lower Fe 3+ /Fe 2 with decreasing temperature, as illustrated in Figs. 1 and 11a. Zircon arrays from slowly cooled intrusive suites plotted in Figs. 9 and 10 tend to show less data dispersion.

Implications for formulation of TiO 2 activity by the MELTS and Rhyolite-MELTS thermodynamic models
Rutile solubility experiments in hydrous granitoid melts performed and compiled by Kularatne and Audétat (2014) underpin their thermodynamic model for TiO 2 activity in rhyolitic melts. They apply their TiO 2 activity model to quartz-hosted glass inclusions in metaluminous and peraluminous rhyolitic volcanic rocks to calculate values of TiO 2 activity that are in good agreement with TiO 2 activities calculated from Fe-Ti-oxide phenocrysts in those samples using the method of Ghiorso and Gualda (2013). However, Kularatne and Audétat (2014) report that titania solubilities predicted by MELTS (Ghiorso and Sack 1995) and Rhyolite-MELTS (Gualda et al. 2012) free-energy-minimisation programs are 1.5-2.9 times higher than actually observed TiO 2 solubilities in all subaluminous and peraluminous melt compositions investigated in their experimental dataset. Kularatne and Audétat's experiments indicate that rutile solubility is insensitive to alumina saturation index in metaluminous melts having Al/Na + K ≥ 1, but rutile solubility in peralkaline melts increases strongly with decreasing Al/ Na + K. The sensitivity of rutile solubility to alumina saturation index in peralkaline melts is consistent with the formation of alkali-titanate complexes as suggested by Dickinson and Hess (1985) and later confirmed by X-ray absorption spectroscopy (e.g., Farges et al. 1996;Farges and Brown 1997).
We infer that the errors in the formulation for a melt TiO 2 by the MELTS and Rhyolite-MELTS thermodynamic models, which underlie the dispute over T-ƒO 2 trends described in the introduction, may be due to failure to adequately account for the varying structural role of Ti with varying alumina saturation index and hydration state of the melt, with Ti being largely a "network forming" component (tetrahedrally coordinated with oxygen) in low-pressure peralkaline melts and predominantly a "network-modifying" component (octahedrally coordinated, as in rutile) in peraluminous melts and high-pressure metaluminous melts. If so, then the enthalpy of fusion of rutile must be combined with variable enthalpy of solvation (mixing) of molten TiO 2 in aluminosilicate melts of varied A/CNK and hydration state (taking Na-OH complexing into account; Kohn et al. 1998), to thermodynamically model a melt TiO 2 in thermally and compositionally varied aluminosilicate melts. Our literature survey indicates that rutile saturation at mid-to upper-crustal pressures in nature is largely limited to peraluminous granitoids. We expect that the configurational behavior of Ti in silicate melts resembles that of Al, which is entirely a tetrahedrally coordinated "network former" in low-pressure peralkaline melts (in which liquidus minerals rich in octahedral Al are unstable), but largely an octahedrally coordinated "network modifier" in peraluminous melts and high-pressure melts that stabilise liquidus minerals rich in octahedrally coordinated Al. By extension of the Burnham (1981) "quasi-crystalline model" for configurational properties of aluminosilicate liquids, MELTS-type thermodynamic models may benefit by transitioning among pressure-and bulk-composition-dependent alternative sets of thermodynamic components of the melt (e.g., using the compositions and fusion enthalpies of CaTS or grossular rather than anorthite, and of jadeite rather than albite) to simplify treatment of configurational entropy and enthalpic mixing properties of aluminosilicate melts over a range of mantle to crustal pressures and varied alumina saturation index. Thermodynamic properties of biotite solid solutions, with particular attention to effective biotite/melt partition coefficients of Fe 3+ and Fe 2+ , may also warrant re-examination in the MELTS and Rhyolite-MELTS programs.

Conclusions
This work presents new evidence that shows that the Fe-Ti-oxide oxybarometer calibration by Ghiorso and Evans (2008) retrieves experimentally imposed values of ƒO 2 in laboratory syntheses of Fe-Ti-oxide pairs to a precision of ± 0.2 log unit, over a large experimental temperature range, without systematic bias up to at least log ƒO 2 ≈ NNO + 4.4. Proceeding from that result, this study addresses a recent dispute (Ghiorso and Gualda 2013;Ghiorso et al 2016;Evans and Bachmann 2013;Evans et al. 2016) regarding the reliability of T-fO 2 trends in the Bishop Tuff magma chamber recorded by titanomagnetite-hemoilmenite phenocryst pairs and trends recovered by Trail et al. (2015) from zircon Ce 4+ /Ce 3+ ratios and Ti-in-zircon temperatures in the Bishop Tuff. Results from our study show that lattice-strain effects cause Ce 4+ /Ce 3+ in zircon to increase strongly as zircon crystallises from cooling magma at constant Ce 4+ /Ce 3+ ratio in the melt. In many natural assemblages described in published literature and interpreted as redox indicators, there is no resolvable change in redox state of the melt over a large range of Ce 4+ /Ce 3+ variation in zircons. We indicate that those magma chambers have crystallised at constant redox state of the melt in systems closed to oxygen (or hydrogen) exchange with the environs. Exceptionally H 2 O-rich arc magmas tend to follow a trend of rising oxidation state of the melt during late stages of fluid-saturated magmatic differentiation at upper-crustal pressures. We suggest that SO 3 and/or H 2 and/or Fe 2+ loss from the melt are responsible. H 2 O-poor magmas tend to follow a T-ƒO 2 trend of decreasing oxidation state of the melt during late stages of magmatic differentiation at upper-crustal pressures, because the precipitating mineral assemblage has higher Fe 3+ /Fe 2+ than coexisting rhyolitic melt. In conclusion, this study corroborates the validity of the Ghiorso and Evans (2008) Fe-Ti-oxide oxybarometer, and supports the methodology of Ghiorso and Gualda (2013) for evaluating TiO 2 activity from quenched Fe-Ti-oxide phenocrysts, and thereby supports use of those TiO 2 activity values in application to Tiin-zircon and Ti-in-quartz geothermometry to address a wide range of fundamental geological problems. We offer suggestions for improving the Rhyolite-MELTS program of Gualda et al. (2012).