Melting of Peridotites through to Granites: A Simple Thermodynamic Model in the System KNCFMASHTOCr

ABSTRACT

A new set of thermodynamic models is presented for calculating phase relations in bulk compositions extending from peridotite to granite, from 0·001 to 70 kbar and from 650°C to peridotite liquidus temperatures, in the system K₂O–Na₂O–CaO–FeO–MgO–Al₂O₃–SiO₂–H₂O–TiO₂–Fe₂O₃–Cr₂O₃ (KNCFMASHTOCr). The models may be used to calculate phase equilibria in partial melting of a large range of mantle and crustal compositions. They provide a good fit to experimental phase relation topologies and melt compositions across the compositional range of the model. Compared with the preliminary model of Jennings, E. S. & Holland, T. J. B. (2015) (A simple thermodynamic model for melting of peridotite in the system NCFMASOCr. Journal of Petrology56, 869–892) for peridotite–basalt melting relations, the inclusion of K₂O and TiO₂ allows for better modelling of small melt fractions in peridotite melting, and in reproducing rutile-bearing eclogite melting at high pressures. An improved order–disorder model for spinel is now incorporated. Above 10 kbar pressure, wet partial melting relations may be significantly affected by the dissolution of silicates in aqueous fluid, so the set of models includes an aqueous low-density silicate-bearing fluid in addition to a high-density H₂O-bearing silicate melt. Oxygen fugacity may be readily calculated for the whole range of bulk compositions investigated, and the effect of water content on melt fO₂ is assessed.

basalt, mantle, partial melting, peridotite, thermodynamics

INTRODUCTION

The modelling of melting relations in rock bulk compositions is an important goal in petrology. Modelling makes it possible to predict mineral + melt assemblages at pressures, temperatures and compositions where existing experimental data must be interpolated. Model predictions may even, with caution, be extrapolated to regimes where experimental constraints would be challenging or time-consuming to obtain. Partial melting plays a major role in both modulating and probing Earth processes, the former by advecting heat and matter, and the latter by influencing seismic velocities. Consequently, the capacity to make accurate predictions of partial melting relations has great significance for Earth system science.

In order to model partial melting, a thermodynamic description is needed for each of the phases involved. The silicate melt is perhaps the most difficult. This is fundamentally because there is no obvious way to write the entropy of the melt, unlike in a crystalline solid, for which the configurational entropy of cation mixing on lattice sites is well defined. However, the problem is exacerbated by the enormous extent of composition space that melts are able to occupy, relative to any solid solution. Thermodynamic modellers have previously handled this problem by modelling subsets of this potential composition space, relevant to melt in a limited range of settings.

Using this approach, partial melting can now be handled in a wide range of contexts using different thermodynamic models. Basaltic melting can be addressed with the MELTS models of Ghiorso & Sack (1995), extended by Ghiorso et al. (2002), or with those of Jennings & Holland (2015), on which the work in this paper is based. Models generating melt of wet granitic composition include those of Holland & Powell (2001) and White et al. (2008), as well as those of Gualda et al. (2012) who extended the MELTS package, as rhyolite-MELTS, in an attempt to deal more effectively with quartz- and K-feldspar-bearing liquids, albeit at the expense of problematic adjustments to the enthalpies of solid quartz and K-feldspar. The models of Green et al. (2016) represent the melting of metabasite rocks, generating a melt more calcic than a typical granitic melt.

Between them, the models listed above can be used to represent partial melting in most geological contexts. The present study describes a simple new thermodynamic model that extends the scope and application of the Holland & Powell (2011) dataset to melting in a wide variety of rocks from peridotite to granite. It modifies and extends the earlier model of Jennings & Holland (2015). It offers, therefore, an alternative calibration to the existing rhyolite-MELTS package, which currently cannot be used to model partial melts in equilibrium with biotite- or hornblende-dominated assemblages, or near-solidus phase equilibria, as explained on the MELTS web page (http://melts.ofm-research.org/).

THERMODYNAMIC MODELS

Full details of all models are given in the Appendix. Thermodynamic data for all of the end-members used here are presented in file tc-ds633.txt on the website listed in the Appendix, as are activity–composition (a–X) relations for the phases involved.

A Monte Carlo approach was adopted for calibrating the models, as described in Green et al. (2016). In the Monte Carlo simulations, filters were applied, both as temperature brackets for phase equilibrium features such as the solidus and for melt compositions expressed in terms of the end-member proportions. Initial estimates for values of parameters (with an associated variation range) were taken, and runs of THERMOCALC were made with relatively wide bracket constraints on phase boundaries and melt compositions. These were progressively refined by narrowing the constraint brackets until all the observations could be satisfied within the estimated uncertainties. Rather than fitting with respect to a large database of individual experimental runs, the Monte Carlo procedure was constrained using a few key runs from each of a set of experimental studies, covering the required range of bulk composition.

This step involves review and interpretation of the available experimental data in order to draw out the salient features of phase relations. When a database of petrological experiments is used to characterize trends in thermodynamic properties, the data typically contain inconsistencies as well as large amounts of redundancy with respect to the variation in properties along some of the relevant axes of pressure, temperature and composition. Meanwhile, along other axes, the variations are highly correlated and are individually under-constrained. Further, such databases generally contain biased weightings, particularly in favour of experiments at 1 atm. Though the interpretation of the phase relations that we adopt is necessarily non-unique, we consider this manual pre-processing to be a valuable step in understanding the data and extracting thermodynamic meaning from them. The constraints interpreted on the basis of these key experiments were either temperature brackets for phase equilibrium boundaries or melt compositions. Constraints were chosen with a view to being able to calculate pseudosections at the bulk compositions of each experimental study that would closely resemble the experimentally-determined phase relations. Monte Carlo runs were made by successively imposing the constraints above until 100 successes at each stage were achieved through varying the set of parameters. The resulting model, while satisfying the constraints, may still be non-unique, but is sufficient within the scope of the objectives in this system. Parameters varied in the fitting are the pairwise interaction energies (W) and modifications to the enthalpies of end-members from the Holland & Powell (2011) dataset (H^mod). Below, the interaction energy for mixing between end-members mm and nn is denoted as W_mm–nn.

Solid phases in KNCFMASHTOCr

Thermodynamic data on mixing models for pyroxenes, olivine and garnet are based on those used in Jennings & Holland (2015). They are extended by inclusion of K and Ti in cpx, Ti in garnet, and Na₂O in opx to allow for better expressing phase relations at small melt fractions. Spinel in the Jennings & Holland (2015) model used a simplified random mixing model based on Bryndzia & Wood (1990). A new model for spinel is presented here (see Appendix) which incorporates order–disorder in a manner similar to that of Hill & Sack (1987) and Sack & Ghiorso (1991a, 1991b). It involves the end-members spinel (MgAl₂O₄), inverse spinel (AlMgAlO₄), hercynite (FeAl₂O₄), inverse hercynite (AlFeAlO₄), inverse magnetite (Fe³⁺FeFe³⁺O₄), normal magnetite (FeFe³⁺Fe³⁺O₄), picrochromite (MgCr₂O₄) and qandilite (MgTiMgO₄). The ternary plagioclase mixing model from Holland & Powell (2003) is used here, allowing incorporation of K₂O in plagioclase and CaO in alkali feldspars. Olivines in the new model incorporate Ca via a monticellite end-member (CaMgSiO₄) and include Fe–Mg disordering on M2 and M1 sites. Ilmenite in the system FeO–Fe₂O₃–TiO₂ is incorporated using the model in White et al. (2000).

Mixing parameters for pyroxenes and garnet were refined in the Monte Carlo fit, at the same time as melt parameters, using experimental equilibria as described below.

A simple thermodynamic model for melts

The melt model is based on the basaltic model of Jennings & Holland (2015). However, the new model is somewhat more complex, having been modified to allow for the enormous compositional range spanned by granitic, tonalitic and basaltic melts. A Temkin-type model is adopted, in which the melt structure can be considered to consist of mixing sites, analogous to a solid solution, but with the number and occupancy of sites varying as a function of melt composition. The sites, and the mixing units contributed to each site by the melt end-members, are shown in Table 1. On the site designated F, random mixing of the nine species AlSi₂, AlSi, Si^p, ${\text{Si}}_2^{\mathrm{o}}$⁠, ${\text{Si}}_4^{\mathrm{q}}$⁠, Cr, Ti, Fe³⁺ and CaAl₂Si takes place, representing the networking behaviour of silicate melts. The three different types of Si units reflect the presence of chains (Si^p), isolated tetrahedra (⁠${\text{Si}}_2^{\mathrm{o}}$⁠) and 3-dimensional network (⁠${\text{Si}}_4^{\mathrm{q}}$⁠) elements. Random mixing of Mg, Fe, Al and Ca are assumed on separate, octahedral-like, M-sites. The A-sites are occupied by the alkali elements Na and K, which form large cations. A hypothetical site V, of multiplicity 2, is composed of just vacancies □ and H. This makes it possible to represent H₂O in a way that maintains compatibility with earlier work on granitic melts (e.g. Nicholls, 1980; Holland & Powell, 2001).

Table 1:

Melt end-members used in this study

	M	F	A	V
jdL	–	AlSi2	Na	□2
kjL	–	AlSi2	K	□2
slL	Al	AlSi	–	□2
woL	Ca	Sip	–	□2
foL	Mg4	${\text{Si}}_2^{\mathrm{o}}$	–	□2
faL	Fe4	${\text{Si}}_2^{\mathrm{o}}$	–	□2
qL	–	${\text{Si}}_4^{\mathrm{q}}$	–	□2
ekL	–	Cr	–	□2
tiL	–	Ti	–	□2
hmL	–	Fe3+	–	□2
h2oL	–	–	–	H2
ctL	–	CaAl2Si	–	□2

Cations are allocated to M, F and A sites. See text for details.

The melt model uses 12 end-members (Table 1; jdL = NaAlSi₂O₆, kjL = KAlSi₂O₆, woL = CaSiO₃, slL = Al₂SiO₅, ctL = CaAl₂SiO₆, foL = Mg₄SiO₈, faL = Fe₄SiO₈, qL = Si₄O₈, ekL = CrO_1.5, tiL = TiO₂, hmL = FeO_1.5, h2oL = H₂O). The end-member ctL is a fictive species (the Gibbs energy of ctL is constructed from ctL = woL + slL-1/4qL) whose amount in the melt is found by solving the equilibrium relation for the internal melt reaction ctL = woL + slL-1/4qL. The model is treated as a slightly asymmetric solution (van Laar), with ideal entropic mixing on sites, and a generally non-ideal enthalpy of mixing between pairs of end-members.

Compared with the Jennings & Holland (2015) model, the changes are: 1) diopside liquid has been replaced by wollastonite liquid; 2) sillimanite liquid is introduced; 3) quartz, forsterite and fayalite liquid end-member formula sizes have doubled; 4) K-jadeite and H₂O liquid are incorporated; and 5) an additional ctL end-member is introduced. The end-members and their model structural elements are as shown in Table 1. The formulae for qL, foL and faL were doubled so that relations in more polymerised granitic melts could be modelled well. The fictive species ctL was introduced to lower the Gibbs energy in some melt composition regions, and was essential in modelling plagioclase phase relations in granitic melts.

The model was developed in two stages. Firstly, the H₂O-absent system was calibrated using the experimental constraints discussed in Jennings & Holland (2015), as well as many additional ones given in the Appendix. Secondly, H₂O-bearing equilibria were added using the constraints from hydrous experiments listed in the Appendix. The fitting made use of the melting temperatures and melt compositions of KLB-1 peridotite (Takahashi, 1986; Takahashi et al., 1993; Davis et al., 2009, 2011), G2 quartz–rutile eclogite (Pertermann & Hirschmann, 2003a, 2003b), Leg 45–395A-8–1-9 basalt (Fujii & Kushiro, 1977), RE46 Icelandic basalt (Yang et al., 1996), Mix1G pyroxenite (Hirschmann et al., 2003), GS104–2-1 basalt (Tormey et al., 1987) and ARP74 basalt (Fujii & Bougault, 1983) to improve the Gibbs energy contribution of each end-member in the H₂O-absent system. To extend the model to more Si-rich compositions, the melting temperatures of albite–quartz and sanidine–quartz as well as the dry R1 granite composition used by Whitney (1975) are incorporated in the fitting together with the many other equilibrium constraints listed in the Appendix.

For bulk compositions containing iron, a bulk composition ferric: ferrous ratio needs setting via application of an appropriate Fe³⁺/∑Fe ratio. For 1 atm experiments (e.g. RE46 basalt) Fe³⁺/∑Fe = 0·09 was used, based on QFM-buffered experiments. For high-pressure piston-cylinder experiments values of Fe³⁺/∑Fe = 0·033, or about one third of that for QFM were used in the absence of good experimental evidence. A value of 0·033 was also used in Jennings & Holland (2015) for mantle peridotite. The fitted melt parameters are not affected significantly by varying these ratios, as the hmL end-member is determined from the experimentally determined Fe³⁺/∑Fe ratio in melts by Kress & Carmichael (1991), Jayasuriya et al. (2004) and O’Neill et al. (2006).

In the Monte Carlo fitting, the parameters that were varied were the H^mod terms of all end-members and the regular solution interaction energies (Ws) among them. After the initial stages where pyroxene terms were tuned, only melt parameters were varied in the final fits. Some interaction energies were initially taken as known, using values from binary and ternary subsystems, and then allowed to vary to see if any adjustments were needed. Some Ws benefitted from a modest pressure-dependence.

Silicate-bearing aqueous fluid

Below 10 kbar, it is sufficient to treat aqueous fluid as pure H₂O. However, as discussed above, at higher pressures aqueous fluid is able to dissolve significant quantities of silicate material. A new aqueous fluid model was, therefore, needed in order to incorporate constraints on wet melting from above 10 kbar (up to 25 kbar in some instances). The new aqueous fluid model allows fluid composition to depart from pure H₂O by incorporation of small amounts of the end-members from the melt model, excluding the fictive end-member ctL. Unlike the melt model, in which mixing takes place on multiple sites, the aqueous fluid model involves molecular mixing. The H₂O end-member and each melt end-member is treated as a single molecule that mixes with the other molecular end-members on a single site. The fluid is considered to be a dilute mixture in which all melt components mix ideally amongst themselves (e.g. W_jdL–kjL = 0), but each end-member mixes non-ideally with H₂O (W_jdL–h2oL ≠ 0, W_kjL–h2oL ≠ 0, …,). The non-ideal interaction energies must be calibrated.

Values for W_qL–h2oL, W_jdL–h2oL, W_kjL–h2oL, W_woL–h2oL, and W_foL–h2oL were found by initially fitting to wet melting experiments for quartz, albite, sanidine, wollastonite and forsterite, respectively. These values were then used as starting values in the Monte Carlo fitting of melting equilibria (see Appendix) and adjusted. W_slL–h2oL was determined from the fitting process. The remaining Ws were assigned arbitrary values similar to those of W_foL–h2oL. The H^mod terms were adjusted in the Monte Carlo fitting of the H₂O-bearing system. Although this simple model is crude, it does at least allow for some solubility of silicates in the fluid, increasing with pressure, and appears to yield solubilities that are consistent with measurements to modest (10–12 kbar) pressures: Kennedy et al. (1962) for quartz–H₂O; Goranson (1931, 1932); Morey & Hesselgesser (1951); Burnham & Jahns (1962); and the data of Stewart in Clark (1966), Yoder et al. (1957), Lambert et al. (1969) for albite–H₂O and san–H₂O; and additionally Burnham (1960) and Anovitz & Blencoe (1997, 1999) for albite–sanidine–quartz–H₂O. Examples for quartz–H₂O and albite–H₂O are shown in Fig. 1. This aqueous model should be used only as an approximate vehicle to lower the H₂O activities in calculations and not to estimate mineral solubilities. Although this approximation is an improvement over assuming pure H₂O, a more realistic model awaits development.

Fig. 1

Calculated and experimental compositions of aqueous fluids and melts coexisting with (a) quartz and (b) albite. Experimental data from Kennedy et al. (1962) for quartz, and from Goranson (1931, 1932), Morey & Hesselgesser (1951), Burnham & Jahns (1962) and Stewart (in Clark, 1966). The simple model fails to generate critical points where fluid and melt compositions approach one another.

Although the water contents of melt are not perfectly modeled for pure albite and pure quartz wet melting, the water contents of haplogranitic melts (Qz₂₈Ab₃₈Or₃₄) as given by the experiments of Holtz et al. (1992, 1995) are in good agreement with the aqueous silicate melt and solute-bearing fluid model used here (Fig. 2).

Fig. 2

Calculated solubility of H₂O in haplogranite granite melt (Qz₂₈Ab₃₈Or₃₄ wt %) compared with the experimental data of Holtz et al. (1995).

EXAMPLE CALCULATIONS

Mantle melting using KLB-1 as a bulk composition

The supra-solidus phase relations are an improvement over the earlier models of Jennings & Holland (2015), and show a much simpler sequence that matches experiments well. The new models include K₂O and TiO₂ in the melts, and they now agree with the liquidus temperatures of Takahashi (1986) and Takahashi et al. (1993). The change in end-members used, together with the new spinel model, has eliminated the appearance of small amounts of spinel which persisted up to liquidus temperatures in the earlier work. Figure 3 shows the revised P–T pseudosection for KLB-1 composition. The widening of the garnet + spinel peridotite field at low temperature is more marked, but this is highly sensitive to both ferric iron and chromium content. Orthopyroxene disappears along the solidus at pressures near 30 kbar rather than 40 kbar as in the earlier model, a change stemming from the improved fitting of garnet and pyroxene compositions as a function of pressure and temperature. As noted earlier in Jennings & Holland (2015), the exact pressure of opx disappearance remains highly sensitive to bulk composition variations.

Fig. 3

P–T pseudosection calculated for KLB-1 in the KNCFMASTOCr system. Phases: g, garnet; sp, spinel; pl, plagioclase; ol, olivine; cpx, clinopyroxene; opx, orthopyroxene; liq, liquid. In comparison with the earlier study of Jennings & Holland (2015) the spinel–garnet transition widens more dramatically at low temperatures and spurious spinel does not appear at high temperatures. Note that the model, as that of Jennings & Holland (2015), predicts a negative slope for the solidus just below the plagioclase to spinel transition pressure.

The models predict the sensitivity of the solidus to Na₂O and K₂O (Fig. 4) rather well in comparison to the 15 kbar experiments (cf. Hirschmann, 2000). The diagram shows the good agreement between the calculated and experimental solidus of KLB-1 and reproduces the trends in lowering of melting temperatures with addition of alkalis in experimental studies. The dashed lines were constructed by addition of K₂O and Na₂O to the bulk composition.

Fig. 4

Lowering of solidus temperature of KLB-1 at 15 kbar with addition of Na₂O and K₂O. The circle is the calculated KLB-1 solidus and squares are experimental compositions: MPY, T’Q (Robinson & Wood, 1998); KLB1, HK-66, PHN-1611 (Hirose & Kushiro, 1993; Kushiro, 1996).

Melting of Icelandic basalt RE46

Yang et al. (1996) performed experiments on melting and crystallization of an Icelandic basalt at ambient pressure with oxygen fugacity controlled by the QFM buffer. They documented the crystallization sequence liquid–olivine–olivine + plagioclase–olivine + plagioclase + clinopyroxene over the interval 1245°C down to 1198°C. A calculated pseudosection for RE46 (with Fe³⁺/∑Fe = 0·09) is shown in Fig. 5, along with the experimental temperatures for incoming of olivine, plagioclase and pyroxene at 1 bar pressure. The model shows the change of liquidus phase from olivine to pyroxene with increasing pressure, a feature commonly reported in experimental studies of such basalts (e.g. Fujii & Kushiro, 1977; Fujii & Bougault, 1983).

Fig. 5

P–T pseudosection calculated for RE46 Icelandic basalt in the KNCFMASTOCr system. Phases: pl, plagioclase; ol, olivine; cpx, clinopyroxene; ilm, ilmenite; liq, liquid. The black bars indicate incoming of olivine, plagioclase and clinopyroxene with decreasing T in the experiments of Yang et al. (1996).

Melting of G2 eclogite at 30 kbar

Pertermann & Hirschmann (2003a)[COMP: pls tag Pertermann & Hirschmann (2003a) as a reference] have determined melting relations for an eclogite at 30 kbar. Melting begins at 1300°C and rutile and then quartz melt out by 1360°C and both garnet and clinopyroxene remain to the liquidus around 1500°C. Calculations with the new model reproduce this behaviour quite well (Fig. 6a). On an up-temperature path at 30 kbar, the model predicts that quartz melts out at 1330°C and garnet outlasts cpx, whereas in the experiments garnet disappears just before cpx, but the difference is probably within combined experimental and model calculation errors. The compositions of the liquid are shown in Fig. 6b and are in good agreement with experiments, except for FeO and CaO that fall slightly low and high respectively in the calculations.

Fig. 6

Calculations for eclogite bulk composition G2 (Pertermann & Hirschmann, 2003a). (a) Modebox (calculated mineral modal proportions) at 30 kbar. The experimental decreasing-temperature sequence (Pertermann & Hirschmann, 2003a) shows liq–cpx–gt–q–ru–(liq) differing from calculations only in cpx being the liquidus phase. (b) Compositions of melt in eclogite bulk compositions. MgO ranges from 2–8 wt % as temperature ranges from 1300–1500°C. Filled circles: G2 bulk composition (Pertermann & Hirschmann, 2003a).

Melting of high-Al basalt composition TA-1

Baker & Eggler (1983) performed experiments on an Aleutian high-Al basalt, bulk composition AT-1 in the range 0–30 kbar. The phase relations were somewhat complex and difficult to decipher in the experiments. At 1 atm the crystallization sequence was determined as plagioclase (1250°C), olivine (1150°C) and cpx (1110°C), with olivine only stable below around 10 kbar (Fig. 7a). The solidus was determined at 2 kbar and 910°C. Just below 20 kbar pressure there is a change from plagioclase to garnet on the liquidus. Figure 7a delineates the incoming and outgoing of garnet, pyroxene, plagioclase and olivine, as proposed by Baker & Eggler (1983) on the basis of their experiments.

Fig. 7

Experiments on high-Al basalt AT-1. (a) Experiments: open symbols: noble-metal-encased samples; filled circles: graphite-encased samples. Figure adapted from Baker & Eggler (1983). (b) Calculated boundaries (from Fig. 8) for AT-1 to compare with experiments. The dashed line connects the pigeonite-cpx solvus crest to the cpx-out boundary where garnet joins. Although their diagram, reproduced here, shows the liquidus at 1250°C, their data at 1 atm show it to be at 1270°C, almost identical to the calculations. The more reduced graphite-encased samples and more oxidized noble-metal-encased samples may differ slightly in their phase boundaries.

Predictive calculations with the new model (neither the experiments of Baker & Eggler (1983), nor any other high-Al basalts were used in the fitting) reproduce these features rather well and can help make sense of the complexities. In Fig. 7b we show the equivalent phase limits for comparison with the interpretation of Baker & Eggler (1983). Near 20 kbar, the plagioclase-out and garnet-out curves intersect and run sub-parallel to the clinopyroxene-out curve, rather than all three curves intersecting at a point as in Fig. 7a. A complex region involving pigeonite is located between 4 and 9 kbar. The experimental diagram of Baker & Eggler (1983) has a region marked as possible opx, which, according to our calculations, may reflect inverted pigeonite on cooling. The calculations match the experimental solidus at 2 kbar, the loss of olivine at around 10 kbar, and the appearance of garnet on the liquidus at around 18 kbar. The complete pseudosection is shown in Fig. 8.

Fig. 8

Calculated P–T pseudosection for AT-1 high-Al basalt. Phases: g, garnet; cpx, clinopyroxene; pl, plagioclase; ru, rutile; q, quartz; ilm, ilmenite; pig, pigeonite; opx, orthopyroxene, liq, liquid. The main features were compared with the experiments of Baker & Eggler (1983) in Fig. 7. Interesting features are: the solvus top for cpx and pig at around 1175°C and the field of pigeonitic pyroxenes in the region 1100–1230°C and 4–13 kbar.

Wet melting of ksp-q and ab-q to high P

As an example of the new melt model applied to more silica (and H₂O)-rich compositions, the melting of albite + quartz and sanidine + quartz in both wet and dry situations is shown in Fig. 9. As a consequence of allowing dissolved silicate in the aqueous fluid, the wet melting curves are modeled well up to 25 kbar. In contrast, by assuming pure H₂O fluid, the curves deviate to lower temperatures progressively above 10 kbar. Application to higher pressures awaits a more sophisticated model for silicate solubility.

Fig. 9

Calculated wet and dry melting of feldspar + quartz to high pressures. (a) Albite + quartz melting: experimental data from Schairer & Bowen (1947, 1956), Tuttle & Bowen (1958), Shaw (1963), Luth et al. (1964), Luth (1969) and Boettcher & Wyllie (1969). (b) Sanidine + quartz melting: experimental data from Schairer & Bowen (1947, 1956), Tuttle & Bowen (1958), Shaw (1963), Luth et al. (1964), Luth (1969), Boettcher & Wyllie (1969), Lambert et al. (1969), Bohlen et al. (1983) and Huang & Wyllie (1974, 1975).

Granite wet/dry melting. Modelling of Whitney (1975) experiments

Whitney (1975) made pioneering measurements on melting of a haplogranitic melt composition at 2 and 8 kbar as a function of water content. Figure 10a shows Whitney’s data and his interpretation of phase equilibria at 2 kbar for water contents up to 12 wt % H₂O. Figure 10b shows the calculations for Whitney’s granite composition at 2 kbar for comparison. The melting temperatures are a few degrees low, but in acceptable agreement, and the water contents at saturation agree well. The experimental location of the pl + ksp + q + liq field at low temperatures is only poorly constrained. The new model appears to perform remarkably well in covering the range in compositions from dry peridotite melting at temperatures of up to 1900°C to wet granitic compositions below 600°C (see Figs 3, 8–10).

Fig. 10

Comparison of (a) experiments of Whitney (1975) and (b) calculations for haplogranite at various water contents at 2 kbar.

Wet melting of tonalite 101

Application of the model to wet melting of natural compositions is illustrated with respect to the melting of tonalite 101 from the Sierra Nevada batholith in experiments by Piwinskii (1968) and later at higher pressures by Stern et al. (1975). A calculated pseudosection for tonalite 101 with 15 wt % H₂O (as in the experiments) is presented in Fig. 11a. We have only used the mineral phases observed in the experimental studies with addition of a small amount of ilmenite to satisfy the bulk composition. Solidus and liquidus are shown with bolder lines. With 15 wt % H₂O, nominally water-saturated, the solidus is in general agreement with the schematic interpretation of Stern et al. (1975), Fig. 11b, which was based on Piwinskii’s (1968) data at 1–3 kbar. Because very small amounts of magnetite are involved at the calculated liquidus, an ‘effective liquidus’ would correspond to the first appearance of plagioclase. The presence of magnetite and/or ilmenite is highly sensitive to both the bulk Fe³⁺/Fe²⁺ ratio and to the value for the melt W_hmL–h2oL interaction energy. The amounts of magnetite and ilmenite are small (less than 1%) and their stability limits quite uncertain. Because the bulk composition contains 15 wt % H₂O, the magnitude of W_hmL–h2oL has a large effect on magnetite stability, and the value used was chosen to make magnetite stable at the liquidus so as to agree with the experiments of Piwinskii (1968) on several similar granites and granodiorites. Reducing W_hmL–h2oL from 60 to 50 kJ/mol would have the effect of just taking magnetite stability from the liquidus. Because of the uncertainties attending the magnetite and ilmenite stabilities, Fig. 11c has been drawn to show the non-oxide phase relations.

Fig. 11

(a) Calculated P–T pseudosection for tonalite 101 (Piwinskii, 1968) with 15 wt % H₂O. (b) Interpretation of Stern et al. (1975) for tonalite 101 phase relations based loosely on Piwinskii (1968) data at 1–3 kbar and high-pressure data of Stern et al. (1975). Note that Piwinskii (1968) does not observe pyroxenes in his study. (c) Phase limits for principal phases taken from (a) to facilitate comparison with (b). See text for discussion. (d) Phase relations in the range 4–6 kbar for 101 bulk composition with 0·8 wt % H₂O to show dramatic differences due to water undersaturation, particularly development of pyroxenes.

Figure 11a also shows that fluid is exhausted at higher pressures and temperatures (above the dashed curves in Fig. 11a, c). This, and the sensitivity of accessory and minor phases to small variations in bulk composition, especially in such H₂O-rich conditions, suggests difficulties might be anticipated in experimental equilibration of such assemblages and in the interpretation of the run products. Several features are different when compared with the experimental phase diagram: biotite and hornblende are less stable and the presence of pyroxenes is limited to rather low pressures, although in the experiments of Piwinskii (1968) pyroxenes were not observed, although they have been inferred in the interpretation of Stern et al. (1975) in Fig. 11b.

The lower temperature breakdown of biotite and hornblende in the calculations may be related to inadequacies in the activity models for these phases, or to the possible presence of fluorine in the natural materials used in the experiments. It may also be possible that capsule leakage in nominally water-saturated experiments is partly responsible for differences between calculated and experimental phase relations. In calculations with only 0·8 wt % H₂O in the range 4–6 kbar (Fig. 11d) the assemblages are fluid-absent just above the solidus and biotite and quartz stabilities are considerably extended relative to fluid saturation, as is the liquidus temperature. In addition, two pyroxenes are present over a large range of the melting interval. Excess H₂O is unlikely to be a common scenario for the generation of tonalitic rocks and the common occurrence of two pyroxenes and hornblende is typical in those generated under granulite facies conditions (Palin et al., 2016). The diagrams using an earlier melt model in Palin et al. (2016) are not recalculated here as the current model reproduces them quite closely, albeit with differences at high temperatures because the new melt model includes ferric iron and titanium and, therefore, allows calculations up to the liquidus.

Melting of average pelite in KFMASH

As a final example, we compare the melt model with the earlier calculations of White et al. (2014) for an average KFMASH pelite composition. The relationships among the phases garnet–cordierite–biotite–sillimanite–K-feldspar–orthopyroxene–quartz and melt are superimposed as dashed lines on Fig. 12, which is taken from Fig. 4 of White et al. (2014). Phase relations were calculated with the new melt and fluid solubility model, but keeping all other phase properties as in White et al. (2014), so that the effect of the new models may be compared. The position of the solidus below 6 kbar is virtually unaffected, and the principal differences to be seen are the enhanced upper stability of K-feldspar by around 20°C, and the small lowering of the solidus temperature at high pressures.

Fig. 12

Calculated P–T pseudosection for KFMASH pelite from figure 4 of White et al. (2014), with fields from the new melt and fluid model superimposed as dashed lines for comparison. Note that the ksp-out boundary has moved up-temperature by around 20°C.

The mixing properties of hydrous cordierite required adjustment in this study and were recalibrated using the garnet–cordierite Fe–Mg exchange equilibria of Aranovich & Podleskii (1983) and Perchuk & Lavrent’eva (1983). The large W (9 kJ) for mixing of hydrous Mg-cordierite with anhydrous Fe-cordierite, as in White et al. (2014), has now been set to 2 kJ and the value for Fe–Mg mixing has been lowered to 2·5 kJ (from 8 kJ). These changes are necessary to preserve agreement with the compositions of coexisting garnet and cordierite in the experiments. Use of the new cordierite model in Fig. 12 would cause an upward shift in pressure for fields containing cordierite and garnet by around 0.5 kbar.

OXYGEN FUGACITY OF MELTS

As outlined in Jennings & Holland (2015) the oxygen fugacity of dry melts may be estimated quite well if the Fe³⁺ to Fe(total) ratio is known. The model allows calculation of the Fe³⁺/ΣFe ratio of melts through the internal equilibrium among melt end-members 4 hmL + qL = 2 faL + O₂ and it was found that, given the interaction energies in the melt which had been determined from the experimental phase equilibrium data, only one parameter needed refining. This parameter was the enthalpy increment (H^mod) to the hmL end-member, which was fitted to the measurements reported in Kress & Carmichael (1991), Jayasuriya et al. (2004) and O’Neill et al. (2006). The quality of the fit may be seen in Fig. 13. The experimental melts measured in Kress & Carmichael (1991) range in logfO₂ from -8.5 to -0.6 and in temperature from 1200 to 1630°C. The calculations require solving for the proportion of the internal end-member ctL and must be done iteratively. We offer a program (fO2melt) that accepts melt compositions (including an estimate of Fe₂O₃) and pressure and temperature and which outputs logfO₂ as well as density, entropy and Gibbs energy of the melt.

Fig. 13

Plot of calculated against observed logfO₂ for melts experimentally determined by Kress & Carmichael (1991), Jayasuriya et al. (2004) and O’Neill et al. (2006). Open circles: anhydrous melt experimental data. Filled circles: melt experimental compositions with 10 wt % H₂O added. Note the offset increases slightly with oxidation state.

Oxygen fugacity estimation may also be extended to hydrous melt compositions. Figure 13 shows that the effect of adding 10 wt % H₂O to all the melt compositions causes a shift in logfO₂ by 0·7 log units at oxidizing conditions and up to nearly 1·5 log units at more reducing conditions. The accuracy of estimating logfO₂ in wet melts is tempered largely by our poor knowledge of the mixing parameter for W_hmL–h2oL. In this study, a value of 60 kJ/mol was found to bring magnetite onto the liquidus for tonalite 101 (see above), but this is not well constrained. A smaller value would reduce the oxygen fugacity offset between dry melts and melts with 10 wt % H₂O. Regardless of this uncertainty, influx of H₂O into a magma would cause an increase in its intrinsic logfO₂ and removal of water the reverse. Other factors such as gain or loss of hydrogen might be significant, but just adding H₂O into a magma would require a reduction in the Fe³⁺ to Fe(total) ratio if the oxygen fugacity were to be buffered by, for example, magnetite + quartz + ferrous silicate at close to the QFM buffer.

SUMMARY

The revised melt model is remarkable in being able to simultaneously reproduce phase equilibrium features (pressures and temperatures of melting, assemblage field boundaries, etc.) in compositions as diverse as lherzolite, basalt, eclogite, tonalite, granite and pelite under wet and dry conditions. It is based on a simplified hypothetical pseudo-structure (framework elements and interstitial cations in the melt) coupled with a nearly symmetrical solution model among a chosen end-member set. It may be unreasonable to expect perfect agreement with experimental phase equilibrium data given the over-simplified assumptions used in the model. Nevertheless, the model appears capable of describing a very wide range of features in diverse bulk compositions and should provide a useful tool for investigating a variety of petrological problems. It is hoped that this model will be modified in future as more is understood about the thermodynamics of mixing in molten silicate solutions.

These melt models have their limitations. The range in compositions currently encompasses the more common natural rock suites (peridotite through basalt and andesite to granite/rhyolite and pelitic compositions) but the models have not been tested on highly alkaline rocks. Additionally, the uncertainty in dealing with solubility of silicate in aqueous fluid limits use of the models to less than 20 kbar, with caution extendable to 30 kbar. Oxidation states of H₂O-bearing melts may be assessed approximately, but refinement awaits further calibration.

ACKNOWLEDGEMENTS

We thank Paul Asimow for his perceptive and helpful comments, and Eleanor Jennings and John Maclennan for valuable discussions. Any remaining errors of omission or commission are entirely ours.

FUNDING

This work was supported by a Natural Environment Research Council grant (NE/ J021539/1) to John Maclennan and T.J.B.H.

APPENDIX

The thermodynamic mixing properties of solid and melt solutions are provided here and in Jennings & Holland (2015). The updated thermodynamic dataset, tc-ds633.txt, a recent version of THERMOCALC (for OS X, Linux and Windows) and example files for running phase equilibrium calculations may be downloaded from http://www.esc.cam.ac.uk/research/research-groups/research-projects/tim-hollands-software-pages/thermocalc. Also available is a small program fO2melt taking rock analyses to yield logfO2, density, and entropy of melts.

Melt model

The melt model is discussed in the text (see Table 1) and is described in terms of 10 compositional variables (wo, jd, kj, sl, fo, fa, hm, ek, ti, h) representing the 11 external end-members (qL, woL, jdL, kjL, slL, foL, faL, hmL, ekL, tiL, h2oL), where the amount of qL is found by difference. The thermodynamic description is complemented by addition of an additional internal (associate) species CaAlSi₂O₆ (ctL). The 12 end-member proportions (with ∑ p_i = 1) are

$$\begin{array}{l}{p}_{\text{qL}}=1-s(1+3/4y)+y\\ p_{\text{woL}}=wo(1+3/4y)-y\\ p_{\text{jdL}}=jd(1+3/4y)\\ p_{\text{kjL}}=kj(1+3/4y)\\ p_{\text{slL}}=sl(1+3/4y)-y\\ p_{\text{foL}}=fo(1+3/4y)\\ p_{\text{faL}}=fa(1+3/4y)\\ p_{\text{hmL}}=hm(1+3/4y)\\ p_{\text{ekL}}=ek(1+3/4y)\\ p_{\text{tiL}}=ti(1+3/4y)\\ p_{\mathrm{h}2\text{oL}}=h(1+3/4y)\\ p_{\text{ctL}}=y\end{array}$$

where s = ek + fa + fo + hm + jd + kj + sl + wo + ti + h and y = p_ctL. As discussed in the text, the proportion of the ctL end-member (y) is determined from the internal equilibrium relation ctL = woL + slL - 1/4qL along with all the other mineral–melt equilibria in the system.

In terms of the A, M, F and V sites of Table 1 the numbers and sum of species elements on each site are:

$$\begin{array}{cc}\mathrm{A} \text{sites}: \hfill & \text{Na}=jd; \mathrm{K}=kj;\hfill \\ \hfill & \sum \mathrm{A}=jd+kj\hfill \end{array}$$

$$\begin{array}{cc}\mathrm{M} \text{sites}: \hfill & \text{Mg}=4fo; \text{Fe}=4fa; \text{Ca}=wo;\text{Al}=sl;\hfill \\ \hfill & \sum \mathrm{M}=4fo+4fa+wo+sl\hfill \end{array}$$

$$\begin{array}{cc}\mathrm{F} \text{sites}:\hfill & {\text{AlSi}}_2=p_{\text{jdL}}+p_{\text{kjL}}; \text{AlSi}=p_{\text{slL}};{ \text{Si}}^{\mathrm{p}}=p_{\text{woL}};{ \text{Si}}_2\mathrm{o}=p_{\text{foL}}+p_{\text{faL}};{ \text{Si}}_4\mathrm{q}=p_{\text{qL}}; \text{Cr}=p_{\text{ekL}}; \text{Ti}=p_{\text{tiL}},\hfill \\ \hfill & {\text{CaAl}}_2\text{Si}=y;{ \text{Fe}}^{3+}=p_{\text{hmL}};\hfill \\ \hfill & \sum \mathrm{F}=1-p_{\mathrm{h}2\text{oL}}=1-h(1+3/4y)\hfill \end{array}$$

$$\begin{array}{cc}\mathrm{V} \text{sites}:& \mathrm{H}=2h;\hspace{0.17em}\square =2(1-h); \sum \mathrm{V}=2\end{array}$$

The ideal activities are then given by mixing on the various sites as:

$$\begin{array}{ccccc}{a}_{\text{kjL}}& =& \frac{K}{\sum A}\hspace{0.33em}\frac{{\mathrm{A}\mathrm{l}\mathrm{S}\mathrm{i}}_2}{\sum F}{\left(\frac{\mathrm{\square }}{\sum V}\right)}^2& =& \frac{p_{\text{kjL}}}{\sum F}{(1-h)}^2\\ a_{\text{jd}}& =& \frac{Na}{\sum A}\hspace{0.33em}\frac{{\text{AlSi}}_2}{\sum \mathrm{F}}{\left(\frac{\mathrm{\square }}{\sum V}\right)}^2& =& \frac{p_{\text{jdL}}}{\sum F}{(1-h)}^2\\ a_{\text{woL}}& =& \frac{\mathrm{C}\mathrm{a}}{\sum M}\hspace{0.33em}\frac{{\mathrm{S}\mathrm{i}}^{\mathrm{p}}}{\sum F}{\left(\frac{\mathrm{\square }}{\sum V}\right)}^2& =& \frac{wo}{\sum M}\hspace{0.33em}\frac{p_{\text{woL}}}{\sum F}{(1-h)}^2\\ a_{\text{slL}}& =& \frac{\mathrm{A}\mathrm{l}}{\sum M}\hspace{0.33em}\frac{\mathrm{A}\mathrm{l}\mathrm{S}\mathrm{i}}{\sum F}{\left(\frac{\mathrm{\square }}{\sum V}\right)}^2& =& \frac{sl}{\sum M}\hspace{0.33em}\frac{p_{\text{slL}}}{\sum F}{(1-h)}^2\\ a_{\text{foL}}& =& \frac{{\mathrm{M}\mathrm{g}}^4}{{\left(\sum M\right)}^4}\hspace{0.33em}\frac{{\mathrm{S}\mathrm{i}}_2^{\mathrm{o}}}{\sum F}{\left(\frac{\mathrm{\square }}{\sum V}\right)}^2& =& \frac{{\left(4fo\right)}^4}{{\left(\sum M\right)}^4}\hspace{0.33em}\frac{p_{\text{foL}}+p_{\text{faL}}}{\sum F}{(1-h)}^2\\ a_{\text{faL}}& =& \frac{{\mathrm{F}\mathrm{e}}^4}{{(\sum M)}^4}\hspace{0.33em}\frac{{\mathrm{S}\mathrm{i}}_2^{\mathrm{o}}}{\sum F}{\left(\frac{\mathrm{\square }}{\sum V}\right)}^2& =& \frac{{\left(4fa\right)}^4}{{\left(\sum M\right)}^4}\hspace{0.33em}\frac{p_{\text{foL}}+p_{\text{faL}}}{\sum F}{(1-h)}^2\\ a_{\text{qL}}& =& \frac{{\mathrm{S}\mathrm{i}}_4^{\mathrm{q}}}{\sum F}{\left(\frac{\mathrm{\square }}{\sum V}\right)}^2& =& \frac{p_{\text{qL}}}{\sum F}{(1-h)}^2\\ a_{\text{ekL}}& =& \frac{\mathrm{C}\mathrm{r}}{\sum F}{\left(\frac{\mathrm{\square }}{\sum V}\right)}^2& =& \frac{p_{\text{ekL}}}{\sum F}{(1-h)}^2\\ a_{\text{hmL}}& =& \frac{{\text{Fe}}^{3+}}{\sum F}{\left(\frac{\mathrm{\square }}{\sum V}\right)}^2& =& \frac{p_{\text{hmL}}}{\sum F}{(1-h)}^2\\ a_{\text{tiL}}& =& \frac{\mathrm{T}\mathrm{i}}{\sum F}{\left(\frac{\mathrm{\square }}{\sum V}\right)}^2& =& \frac{p_{\text{tiL}}}{\sum F}{(1-h)}^2\\ a_{\text{ctL}}& =& \frac{\mathrm{C}\mathrm{a}{\mathrm{A}\mathrm{l}}_2\mathrm{S}\mathrm{i}}{\sum F}{\left(\frac{\mathrm{\square }}{\sum V}\right)}^2& =& \frac{y}{\sum F}{(1-h)}^2\\ a_{\mathrm{h}2\text{oL}}& =& {\left(\frac{\mathrm{H}}{\sum V}\right)}^2& =& h^2\end{array}$$

Non-ideal interactions are expressed in terms of regular solution energies W_i–j = a + bT + cP in kJ for end-members i and j, as in Table A1. Approximate uncertainties are given in parentheses. Only the non-zero interaction energies are listed.

The enthalpies of end-members are those in the Holland & Powell (2011) dataset, with increments expressed as a difference (H_i^mod = H_i − H_{i, ref}). Where an end-member is made from a linear combination of others H_i^mod is relative to that combination (Table A2). As before, the increment is expressed as H_i^mod = a_i + b_iT + c_iP in kJ, and the van Laar parameter α as a ratio relative to 1·0. Approximate uncertainties from Monte Carlo runs are a_i (0·2), b_i (0·0002), c_i (0·002) α_i (0·05).

The hmL (Fe₂O₃) component was calibrated as discussed in Jennings & Holland (2015). Enthalpy was fitted to measured Fe³⁺/ΣFe of basalts with low K₂O and low P₂O₅ from Kress & Carmichael (1991), where logfO₂ was calculated using the all-melt equilibrium ½ qL + 4 hmL = faL + O₂, using the thermodynamic data and activity model presented here.

Plagioclase

Plagioclase feldspars are treated here using the ternary feldspar model from Holland & Powell (2003).

Olivine

Olivines are treated as simple two-site symmetric solid solutions (M1M2SiO₄) between forsterite (MgMgSiO₄), fayalite (FeFeSiO₄), ordered cfm (MgFeSiO₄) and monticellite (MgCaSiO₄), with non-ideal mixing terms as follows:

$$\begin{array}{ll}{W}_{\text{mo}–\text{fa}}=24 \text{kJ}\hfill & W_{\text{fa}–\text{fo}}=9 \text{kJ}\hfill \\ W_{\text{mo}–\text{fo}}=38 \text{kJ}\hfill & W_{\text{fa}–\text{cfm}}=4.5 \text{kJ}\hfill \\ W_{\text{mo}–\text{cfm}}=24 \text{kJ}\hfill & W_{\text{fo}–\text{cfm}}=4.5 \text{kJ}\hfill \end{array}$$

The W terms are derived from W_Fe–Mg = 4·5 kJ per site following Wiser & Wood (1991). The other parameters were adjusted to fit the olivine–monticellite solvus of Warner & Luth (1973). The Fe–Mg ordering is assumed to be unimportant, so the Gibbs energy of the cfm end-member is given as just G_cfm = (G_fo+G_fa)/2 with no increment.

Clinopyroxene

Clinopyroxene is modelled after Green et al. (2012b) and Jennings & Holland (2015) with additions for Ti (Ca-buffonite) and K (K-jadeite). The end-members used are given in Table A3, with ideal activities being given by a mixing-on-sites model. The thermodynamic data were adjusted to give the small Ti and K contents seen in KLB-1 and in granulite facies metabasites.

The non-ideality is expressed as an asymmetric van Laar model, using the formalism of Holland & Powell (2003). Ws, H_i^mod and asymmetry parameters are given in Tables A4 and A5. Ws involving Fe, Fe³⁺, Na and Cr end-member non-idealities are determined using the approach in Powell et al. (2014). Symbols and values are explained under the melt section above.

Orthopyroxene

Orthopyroxene is also modelled after Jennings & Holland (2015), with addition of Mg-buffonite and ortho-jadeite. The thermodynamic data were adjusted to give the small Ti and Na contents seen in KLB-1 and in granulite facies metabasites. The end-members used are given in Table A6. The non-ideality is expressed as an asymmetric van Laar model, which is given in Tables A7 and A8, with symbols and values as for clinopyroxene.

Garnet

Garnets are modelled as in Jennings & Holland (2015) and Miller et al. (2016), with addition of a Ti end-member (tig) and by including asymmetric mixing. The model interaction energies for garnet are as follows (in kJ):

W_py–alm = 4·0 + 0·10P

W_py–gr = 45·4–0·010T + 0·04P

W_py–andr = 107·0 + 0·01T + 0·036P

W_py–knor = 2·0

W_py–tig = 0

W_alm–gr = 17·0–0·010T + 0·10P

W_alm–andr = 65·0–0·010T + 0·039P

W_alm–knor = 6·0 + 0·01P

W_alm–tig = 0

W_gr–andr = 2·0

W_gr–knor = 1·0 − 0·010T + 0·18P

W_gr–tig = 0

W_andr–knor = 63·0 − 0·010T + 0·18P

W_andr–tig = 0

W_knor–tig = 0

The asymmetry parameters (van Laar) are grossular (2·5) and andradite (2·5), with all others set at unity.

The interaction energies have been modified from those in Jennings & Holland (2015) and adjusted in the Monte Carlo fitting of garnet compositions in the assemblages used.

Spinel

Spinel is modelled as a two-site solid solution in which order–disorder may occur between the octahedral M and tetrahedral T sites. The end-members in Mg–Fe–Al–Fe³⁺–Cr–Ti–O used are as follows

	T	M
sp	Mg	Al2
isp	Al	MgAl
hc	Fe	Al2
ihc	Al	FeAl
nmt	Fe	${\text{Fe}}_2^{3+}$
imt	Fe3+	FeFe3+
pcr	Mg	Cr2
qnd	Mg	MgTi

There are four compositional parameters (x = Fe/(Fe + Mg), y = Fe³⁺/(Al+ Fe³⁺), c = X_{Cr, M}, t = 2X_{Ti, M}) and three order parameters (Q1 = X_{Mg, T} − X_{Mg, M}, Q2 = X_{Fe, T} − X_{Fe, M}, Q3 = X_{Fe3+, T} − X_{Fe3+, M}). The activity model is similar to that of Sack & Ghiorso (1991a, 1991b), but we reduce the entropy of mixing on the M site by a factor of 2 as reasoned in Holland & Powell (1996), following Carpenter & Salje (1994).

Ideal activities are calculated from

$$\begin{array}{l}{a}_{\text{sp}}=X_{\text{Mg},\mathrm{T}}{X}_{\text{Al},\mathrm{M}}{a}_{\text{hc}}=X_{\text{Fe},\mathrm{T}}{X}_{\text{Al},\mathrm{M}}\hspace{1em}\\ a_{\text{nmt}}=X_{\text{Fe},\mathrm{T}}{X}_{\text{Fe}3+,\mathrm{M}}\hspace{1em}{a}_{\text{pcr}}=X_{\text{Mg},\mathrm{T}}{X}_{\text{Cr},\mathrm{M}}\\ a_{\text{isp}}=2X_{\text{Al},\mathrm{T}}{(X_{\text{Mg},\mathrm{M}})}^{1/2}{(X_{\text{Al},\mathrm{M}})}^{1/2}\\ a_{\text{ihc}}=2X_{\text{Al},\mathrm{T}}{(X_{\text{Fe},\mathrm{M}})}^{1/2}{(X_{\text{Fe},\mathrm{M}})}^{1/2}\\ a_{\text{imt}}=2X_{\text{Fe}3+,\mathrm{T}}{(X_{\text{Fe},\mathrm{M}})}^{1/2}{(X_{\text{Fe}3+,\mathrm{M}})}^{1/2}\\ a_{\text{qnd}}=2X_{\text{Mg},\mathrm{T}}{(X_{\text{Mg},\mathrm{M}})}^{1/2}{(X_{\text{Ti},\mathrm{M}})}^{1/2}\end{array}$$

Monte Carlo methods were used to determine the spinel mixing parameters so that the following experimental constraints were met:

1. Ordering of simple end-members (spinel, hercynite, magnetite and magnesioferrite) using data from Harrison & Putnis (1996) and Andreozzi et al. (2000).

2. In MAS, the equilibrium py + fo = opx + sp was made to pass through the brackets of Danckwerth & Newton (1978) and Gasparik & Newton (1984).

3. The solvus between magnetite and hercynite (Turnock & Eugster, 1962) was fitted.

4. The complex ‘quad’ solvus between Al-poor spinels and Al-rich spinels in the system Mg–Al–Fe–Fe³⁺ was fitted to experimental tie lines in Lehmann & Roux (1986).

5. The Fe–Mg distribution between olivine and spinel for Al–Cr spinels of Engi (1983) were fitted.

6. Fitting of corundum–spinel tie lines from Oka et al. (1984).

7. Fitting to activity of FeCr₂O₄ to the data of Petric & Jacob (1982) for the FeCr₂O₄– Fe₃O₄ binary.

8. The solvi between Fe₂TiO₄– FeAl₂O₄ and Mg₂TiO₄– MgAl₂O₄ binary spinels (Muan et al., 1972; Boden & Glasser, 1973).

9. The top of the ulvospinel–magnetite solvus (Price, 1981).

The final set of interaction energies (in kJ/mol) used is

W_sp–isp  −8·2

W_sp–hc  3·5

W_sp–ihc  −13·0

W_sp–nmt  43·2

W_sp–imt  49·1

W_sp–pcr  −5·0

W_sp–qnd  22·5

W_isp–hc  4·4

W_isp–ihc  −6·0

W_isp–nmt  36·8

W_isp–imt  20·0

W_isp–pcr  14·0

W_isp–qnd  21·5

W_hc–ihc  −8·2

W_hc–nmt  18·1

W_hc–imt  49·0

W_hc–pcr  −19·0

W_hc–qnd  35·1

W_ihc–nmt  −4·0

W_ihc–imt  7·6

W_ihc–pcr  −11·0

W_ihc–qnd  9·0

W_nmt–imt   18·1

W_nmt–pcr  11·9

W_nmt–qnd  62·2

W_imt–pcr  −6·4

W_imt–qnd  24·3

W_pcr–qnd  60·0

The increments (H_i^mod values, in kJ/mol) to the dataset Gibbs energies for end-members are given as:

$$\begin{array}{c}\text{isp}:H_i^{\text{mod}}=23.6-0.00576303T\\ \text{ihc}:H_i^{\text{mod}}=23.6-0.00576303T\\ \text{nmt}:H_i^{\text{mod}}=0.00576303T\\ \text{imt}:H_i^{\text{mod}}=0.30\\ \text{qnd}:H_i^{\text{mod}}=-30.0\end{array}$$

In the above sp, isp, hc, and ihc are made ordered before adding the H_i^mod values, which are for disordering. For nmt and imt the H_i^mod values are added to the equilibrium Gibbs energies from the dataset.

Experimental constraints used for the melt model

1. Beginning of melting of plagioclase peridotite at 1 bar, with a temperature bracket of 1120–1180°C (Takahashi et al., 1993).

2. Beginning of melting at 30 kbar of garnet + olivine  + clinopyroxene + orthopyroxene, with a bracket of 1480–1510°C (Takahashi et al., 1993; Davis et al., 2011).

3. Beginning of melting of olivine + garnet + clinopyroxene at 43 kbar (1600–1710°C), 60 kbar (1720–1800°C) and 65 kbar (1740–1800°C) from Takahashi et al. (1993). Pressures for the above constraints were corrected by -5% relative to the nominal experimental values.

4. The compositions of the melt for (2) above was assumed to be that determined by Davis et al. (2011) for initial melt from KLB-1 at 30 kbar. In terms of end-member proportions, the melt composition brackets used were p_woL (0·34–0·37), p_jdL (0·13–0·16), p_slL (0·12–0·15), p_faL (0·05–0·07), p_foL (0·17–0·19), p_hmL (0·0–0·05), p_tiL (0·015–0·06), p_ekL (0·002–0·0033), p_kjL (0·034–0·040), with p_qL making up the remainder, and not imposed.

5. Beginning of melting for G2 (Pertermann & Hirschmann, 2003a, 2003b) at 30 kbar (1280–1320°C), with composition brackets p_woL (0·16–0·20), p_jdL (0·24–0·28), p_slL (0·10–0·14), p_faL (0·019–0·035), p_tiL (0·13–0·15, p_kjL (0·03–0·06). Also, rutile-out (1305–1350°C) and quartz-out (1315–1370°C).

6. MORB Leg 45-395A-8-1-9 basalt (Fujii & Kushiro, 1977), olivine-in (1200–1225°C) and plagioclase-in (1160–1190°C) at 1 bar. At 9 kbar, cpx-in (1250–1280°C) and plagioclase-in (1205–1235°C).

7. RE46 Icelandic basalt (Yang et al., 1996) at 1 bar: olivine-in (1220–1250°C), plagioclase-in (1205–1230°C) and clinopyroxene-in (1180–1210°C). Composition brackets for melt at 1200°C were p_woL (0·42–0·45), p_jdL (0·08–0·11), p_slL (0·16–0·22), p_foL (0·086–0·10), p_faL (0·045–0·06).

8. Liquidus for KLB-1 (Takahashi et al., 1993). Olivine-in at 1 bar (1780–1810°C) and 28 kbar (1900–1935°C). Pressure corrected by -5%.

9. MIX-1G pyroxenite (Hirschmann et al., 2003). Melt composition with garnet + clinopyroxene at 25 kbar and 1500°C: p_woL (0·35–0·41), p_jdL (0·09–0·13), p_slL (0·10–0·24), p_faL (0·04–0·06) ), p_tiL (0·02–0·04, p_kjL (0·001–0·003).

10. MORB (Tormey et al., 1987) at 1 bar. Olivine-in (1165–1220°C).

11. MORB (Fujii & Bougault, 1983). Olivine-in at 1 bar (1240–1285°C), 6·8 kbar (1250–1300°C), plagioclase-in at 6·8 kbar (1210–1290°C), and clinopyroxene-in at 1 bar (1150–1180°C) and 13·5 kbar (1325–1355°C). Pressures corrected by −5%.

12. NAS equilibrium albite + quartz/tridymite = liquid (Schairer & Bowen, 1956; Luth, 1969): 1 bar (1020–1070°C), 8 kbar (1140–1180°C), 12·3 kbar (1200–1230°C), 20 kbar (1270–1320°C) and Bell & Roseboom (1969): 33 kbar (1360–1380°C).

13. KAS equilibrium sanidine + quartz/tridymite = liquid (Schairer & Bowen, 1956; Luth, 1969): 1 bar (980–1000°C), 10 kbar (1090–1230°C), and 20 kbar (1290–1320°C).

14. KNCAS granite R1 dry melting (Whitney, 1975). 2 kbar (1130–1200°C) and 8 kbar (1160–1220°C).

15. KNAS dry granite. Huang & Wyllie (1975): 6·5 kbar (1030–1060°C), 16 kbar (1125–1150°C) and 26·5 kbar (1175–1225°C).

16. NASH system. Albite + quartz + fluid = liquid (Luth et al., 1964). 1 kbar (760–812°C), 2 kbar (730–770°C), 4 kbar (685–695°C), 8 kbar (655–695°C), 10 kbar (637–677°C), 15 kbar (610–640°C), 25 kbar, with jadeite (660–700°C).

17. KASH system. Sanidine + quartz + fluid = liquid (Huang & Wyllie, 1974; Goldsmith & Peterson, 1990). 1 kbar (780–805°C), 2 kbar (750–770°C), 4 kbar (725–745°C), 7 kbar (707–727°C), 10 kbar (700–730°C), 15 kbar (695–730°C), 20 kbar (670–710°C), 25 kbar (670–720°C).

18. SH system. Quartz/cristobalite + fluid = liquid (Kennedy et al., 1962; Stewart in Luth, 1976). 1300°C (0·6–0·9 kbar), 2 kbar (1095–1120°C), 5 kbar (1020–1050°C), 10 kbar (1015–1050°C).

19. KASH system (Huang & Wyllie, 1974). Sanidine–sillimanite–muscovite–quartz–fluid–liquid 5·9 kbar (700–725°C). Sanidine–sillimanite–quartz–muscovite–liquid (fluid-absent): 10 kbar (772–785°C). Sanidine–quartz–muscovite–fluid–liquid: 10 kbar (705–720°C).

20. NKASH system. Sanidine–albite–quartz–sillimanite–fluid–liquid: 4 kbar (620–660°C).

21. CMAS system. Invariant points from Green et al. (2012a). Anorthite + forsterite + clinopyroxene + orthopyroxene + spinel + liquid, 8·5 kbar (1310–1350°C). Garnet + forsterite + clinopyroxene + orthopyroxene + spinel + liquid, 25 kbar (1530–1570°C).

22. CMASH system (Thompson & Ellis, 1994). Zoisite + kyanite + quartz + amphibole + fluid = liquid, 12 kbar (750–820°C), 20 kbar (750–800°C). Zoisite + quartz + amphibole + fluid = clinopyroxene + liquid, 9 kbar (750–870°C), 20 kbar (800–900°C). Amphibole + fluid = garnet + clinopyroxene + orthopyroxene + liquid, 20 kbar (850–900°C). Amphibole + fluid = anorthite + clinopyroxene + orthopyroxene + liquid, 10 kbar (896–950°C). Anorthite + quartz + wollastonite + fluid = liquid, 7 kbar (750–810°C).

23. KNASH haplogranite. Huang & Wyllie (1975). Albite  + sanidine + quartz + fluid = liquid. 10 kbar (615–650°C) and 15 kbar (580–650°C).

24. CASH system (Boettcher, 1970). Anorthite + q + wollastonite + fluid = liquid. 7 kbar (790–820°C).

25. Granite 104 (Stern et al., 1975). Plagioclase-out, 3 kbar (690–720°C).

26. KCNASH system. R1 granite (Whitney, 1975). Plagioclase-out. 2 kbar (880–920°C) for 4 wt % H₂O.

27. Sierra Nevada Granite 104. Piwinskii (1968). Solidus 2 kbar (650–700°C).

28. KFMASH melting, Carrington & Harley (1995). Invariant bi–cd–g–san–q–opx–sill–liq at 8·4–9·5 kbar and 865–900°C.

29. Water solubility in haplogranite melt, 2 kbar (800°C, 6 wt % H₂O), 2 kbar (1150°C, 5·5 wt % H₂O), 5 kbar (800°C, 10·3 wt % H₂O), 5 kbar (1150°C, 10·7 wt % H₂O), 8 kbar (800°C, 13·7 wt % H₂O)

In addition to these melting equilibria, several subsolidus equilibria were included in the Monte Carlo runs to ensure consistency with experimental data on pyroxene and garnet compositions in CMAS. These were

1. CMS system clinopyroxene + orthopyroxene (Mori & Green, 1975; Lindsley & Dixon, 1976; Perkins & Newton, 1980; Nickel & Brey, 1984; Carlson & Lindsley, 1988). Brackets used were 1 bar (925°C), 15 kbar (850, 1100°C), 20 kbar (1100, 1400°C), 30 kbar (1200°C), and 50 kbar (1400°C).

2. CMS system clinopyroxene + orthopyroxene + pigeonite univariant at 15 kbar and 1470°C (Lindsley & Dixon, 1976), di_opx (0·05–0·15), en_cpx (0·3–0·55), en_pig (0·55–0·75).

3. Coexisting clinopyroxene + orthopyroxene + pigeonite in CFMS ‘quad’ system at 15 kbar and 1000°C. (Turnock & Lindsley, 1981; Lindsley, 1983). x_opx (0·57–0·67), Ca_opx (0·0–0·2), x_cpx (0·5–0·65), en_cpx (0·2–0·45), x_pig (0·52–0·70), en_pig (0·60–0·85).

4. CMAS Garnet + clinopyroxene + orthopyroxene (Nickel et al., 1985): 30 kbar, 1400°C, ts_opx (0·06–0·12), ts_cpx (0·07–0·13), di_opx (0·04–0·1), en_cpx (0·29–0·41), Ca_gt (0·11–0·15). Also Perkins & Newton (1980): 24 kbar, 1500°C, ts_opx (0·12–0·2), ts_cpx (0·19–0·22), Ca_gt (0·12–0·15); 20 kbar, 900°C, ts_opx (0·02–0·07), ts_cpx (0·02–0·05), di_opx (0·01–0·025), en_cpx (0·03–0·075), Ca_gt (0·12–0·15); 16 kb, 1000°C, ts_opx (0·08–0·12), ts_cpx (0·07–0·10), di_opx (0·02–0·35), en_cpx (0·07–0·105), Ca_gt (0·12–0·15); 40 kbar, 1000°C, ts_opx (0·012–0·0265), ts_cpx (0·011–0·020), di_opx (0·01–0·03), en_cpx (0·05–0·10), Ca_gt (0·12–0·15).

5. CMAS garnet + forsterite + clinopyroxene + orthopyroxene + spinel (Jenkins & Newton, 1979): 16 kbar, 1000°C, ts_opx (0·08–0·13), ts_cpx (0·07–0·10), Ca_gt (0·12–0·16).

6. CMASCr garnet + forsterite + orthopyroxene + spinel (Klemme & O’Neill, 2000): 33 kbar, 1400°C, ts_opx (0·10–0·16), Cr_opx (0·03–0·05), Cr_gt (0·08–0·20), Cr_sp (0·2–0·5).

Thermodynamic data

Updates to thermodynamic data relevant to this study are provided in Table A9.

Table A1:

Interaction energies for melt end-members (kJ)

WqL–slL = 9·5–0·10P WqL–woL = −10·3 WqL–foL = −26·5– 3·12P WqL–faL = −12·0– 0·55P WqL–jdL = −15·1– 0·13P WqL–hmL = 20·0 WqL–ekL = 0 WqL–tiL = 24·6 WqL–kjL = −17·8– 0·05P WqL–ctL = −14·6 WqL–h2oL = 17·8– 0·61P WslL–woL = −26·5 + 0·85P WslL–foL = 2·2 WslL–faL = 2·5 WslL–jdL = 16·8 WslL–hmL = −5·0 WslL–ekL = 0 WslL–tiL = 15·2– 0·04P WslL–kjL = 7·0 WslL–ctL = 4·0 WslL–h2oL = 23·7–0·94P WwoL–foL = 25·5 + 0·11P WwoL–faL = 14·0 WwoL–jdL = −1·2 WwoL–hmL = 0 WwoL–ekL = 0 WwoL–tiL = 18·0 WwoL–kjL = −1·1 WwoL–ctL = 9·5 WwoL–h2oL = 40·3– 0·86P WfoL–faL = 18·0 WfoL–jdL = 1·5 WfoL–hmL = 0

WfoL–ekL = 0 WfoL–tiL = 7·5 WfoL–kjL = 3·0 WfoL–ctL = −5·6 WfoL–h2oL = 9·4– 1·58P WfaL–jdL = 7·5–0·05 WfaL–hmL = −30·0 WfaL–ekL = 0 WfaL–tiL = 6·7 WfaL–kjL = 10·0 WfaL–ctL = −6·5 WfaL–h2oL = 9·2– 1·58P WjdL–hmL = 10 WjdL–ekL = 0 WjdL–tiL = 16·5 + 0·14P WjdL–kjL = −5·9 WjdL–ctL = 7·6 WjdL–h2oL = −8·3– 0·06P WhmL–ekL = 0 WhmL–tiL = 0 WhmL–kjL = 10·0 WhmL–ctL = 0 WhmL–h2oL = 60·0– 0·66P WekL–tiL = 0 WekL–kjL = 0 WekL–ctL = 0 WekL–h2oL = 30·0– 0·66P WtiL–kjL = 9·0 WtiL–ctL = 0 WtiL–h2oL = 30·0– 0·60P WkjL–ctL = −5·6 WkjL–h2oL = −0·1 + 0·22P WctL–h2oL = 17·3 + 0·05P

Table A2:

${\mathit{H}}_{\mathrm{i}}^{\text{mod}}$ increments (in kJ) to Gibbs energies of end-members in the dataset (tc-ds633) revised from Holland & Powell (2011) and van Laar parameters (α)

Himod	a	b	c	α	ref
qL	0·22	0	−0·059	1·0	4 qL
slL	6·20	0	−0·318	1·2	sillL
woL	−0·45	0	−0·114	1·4	woL
foL	8·67	0	−0·131	2·4	2 foL
faL	13·70	0	−0·055	1·0	2 faL
jdL	12·19	0	−0·089	1·2	abL-qL
hmL	3·30	0	−0·032	1·0	0·5 hemL
ekL	24·85	0	0·245	1·0	0·5 eskL
tiL	5·58	0	−0·489	1·0	ruL
kjL	11·98	0	−0·210	1·0	kspL-qL
ctL	−108·30	0·055	0·053	1·0	woL+silL-qL
h2oL	3·20	−0·0039	0·00087	1·0	h2oL

Table A3:

Clinopyroxene end-members in the model

	M2	M1	T
di	Ca	Mg	Si2
cfs	Fe	Fe	Si2
cats	Ca	Al	AlSi
crdi	Ca	Cr	AlSi
cbuf	Ca	½Mg½Ti	AlSi
cess	Ca	Fe3+	AlSi
jd	Na	Al	Si2
cen	Mg	Mg	Si2
cfm	Fe	Mg	Si2
kjd	K	Al	Si2

Table A4:

Interaction energies for mixing of clinopyroxene end-members (kJ)

Wi–j	a	b	c	Wi–j	a	b	c
di–cfs	25·8	0	−0·03	cats–cfm	27·0	0	−0·1
di–cats	13·0	0	−0·06	cats–kjd	6·0	0	0
di–crdi	8·0	0	0	crdi–cess	2·0	0	0
di–cess	8·0	0	0	crdi–cbuf	0	0	0
di–cbuf	0	0	0	crdi–jd	3·0	0	0
di–jd	26·0	0	0	crdi–cen	52·3	0	0
di–cen	29·8	0	−0·03	crdi–cfm	40·3	0	0
di–cfm	20·6	0	−0·03	crdi–kjd	3·0	0	0
di–kjd	26·0	0	0	cess–cbuf	0	0	0
cfs–cats	25·0	0	−0·1	cess–jd	3·0	0	0
cfs–crdi	38·3	0	0	cess–cen	57·3	0	0
cfs–cess	43·3	0	0	cess–cfm	45·3	0	0
cfs–cbuf	0	0	0	cess–kjd	3·0	0	0
cfs–jd	24·0	0	0	cbuf–jd	0	0	0
cfs–cen	2·3	0	0	cbuf–cen	0	0	0
cfs–cfm	3·5	0	0	cbufvcfm	0	0	0
cfs–kjd	24·0	0	0	cbuf–kjd	0	0	0
cats–crdi	2·0	0	0	jd–cen	40·0	0	0
cats–cess	2·0	0	0	jd–cfm	40·0	0	0
cats–cbuf	0	0	0	jd–kjd	10·0	0	0
cats–jd	6·0	0	0	cen–cfm	4·0	0	0
cats–cen	45·2	0	−0·35	cen–kjd	40·0	0	0
				cfm–kjd	40·0	0	0

Table A5:

${\mathit{H}}_{\mathrm{i}}^{\text{mod}}$ increments to Gibbs energies, and asymmetry parameters α for mixing, of clinopyroxene end-members

Himod	a	b	c	ref	α
di	–	–	–	–	1·2
cfs	2·1	−0·002	0·045	fs	1·0
cats	–	–	–	–	1·9
crdi	−4·90	0	0	cats + kos - jd	1·9
cess	−3·45	0	0	cats + acm - jd	1·9
cbuf	1·75	−0·0012	−0·005	cats + ½per +½ ru - ½cor	1·9
jd	–	–	–	–	1·2
cen	3·5	−0·002	0·048	en	1·0
cfm	−1·6	−0·002	0·0465	½(en + fs)	1·0
kjd	−3·75	0	1·189	jd + san - abh	1·2

H_i^mod = a + b T + c P, in kJ.

Table A6:

Orthopyroxene end-members in the model

	M2	M1	T
en	Mg	Mg	Si2
fs	Fe	Fe	Si2
fm	Fe	Mg	Si2
odi	Ca	Mg	Si2
mgts	Mg	Al	AlSi
cren	Cr	Mg	AlSi
obuf	Mg	½Mg½Ti	AlSi
mess	Mg	Fe3+	AlSi
ojd	Na	Al	Si2

Table A7:

Interaction energies for mixing of orthopyroxene end-members (kJ)

Wi–j	a	b	c	Wi–j	a	b	c
en–fs	7·0	0	0	fm–obuf	0	0	0
en–fm	4·0	0	0	fm–mess	12·0	0	0
en–odi	32·2	0	0·12	fm–ojd	35·0	0	0
en–mgts	12·5	0	−0·04	odi–mgts	75·5	0	−0·84
en–cren	8·0	0	0	odi–cren	20·0	0	0
en–obuf	0	0	0	odi–obuf	0	0	0
en–mess	8·0	0	0	odi–mess	20·0	0	0
en–ojd	35·0	0	0	odi–ojd	35·0	0	0
fs–fm	4·0	0	0	mgts–cren	2·0	0	0
fs–odi	25·54	0	0·084	mgts–obuf	0	0	0
fs–mgts	11·0	0	−0·15	mgts–mess	2·0	0	0
fs–cren	10·0	0	0	mgts–ojd	7·0	0	0
fs–obuf	0	0	0	cren–obuf	0	0	0
fs–mess	10·0	0	0	cren–mess	2·0	0	0
fs–ojd	35·0	0	0	cren–ojd	−11·0	0	0
fm–odi	22·54	0	0·084	obuf–mess	0	0	0
fm–mgts	15·0	0	−0·15	obuf–ojd	4·0	0	0
fm–cren	12·0	0	0	mess–ojd	−11·0	0	0

Table A8:

${\mathit{H}}_{\mathrm{i}}^{\text{mod}}$ increments to Gibbs energies, and asymmetry parameters α for mixing, of orthopyroxene end-members

Himod	a	b	c	ref	α
en	–	–	–	–	1·0
fs	–	–	–	–	1·0
fm	−6·6	0	0	½(en + fs)	1·0
odi	−0·1	0·000211	0·005	di	1·2
mgts	–	–	–	–	1·0
cren	−25·9	0·0155	0·05	mgts + kos - jd	1·0
obuf	3·35	−0·0051	−0·0061	cats + ½per +½ ru - ½cor	1·0
mess	4·80	0	−0·089	mgts + acm - jd	1·0
ojd	18·8	0	0	jd	1·2

H_i^mod = a + b T + c P, in kJ.

Table A9:

Updates to thermodynamic data of Holland & Powell (2011) relevant to this study

• rutile liquid and hematite liquid added (see Jennings & Holland, 2015) • eskolaite liquid added (see Jennings & Holland, 2015) • K-jadeite added (see Jennings & Holland, 2015) • new entropy for knorringite (see Jennings & Holland, 2015) • changed volume for woL to fit 1 bar melting better • almandine entropy changed to 337 J/K, and ilmenite entropy to 107·5 J/K (Dachs et al., 2012) • new volume and compressibility for eskolaite (Kantor et al., 2012); also enthalpy from Klemme et al. (2000) • added Cr-en and Cr-di tschermak endmembers (see Jennings & Holland, 2015) • fitted picr + cor = sp + esk equilibrium from Chatterjee et al., l984 • refitted geik heat capacity • added qandilite (qnd) with data of O'Neill & Scott (2005) • added 2esk + en = 2picr + 2q equilibrium (Klemme & O'Neill, 1997) • hercynite and spinel updated for the new disordering model • changed heat capacity and other values for h2oL • andradite entropy changed to 324 J/K (Geiger et al., 2017) • pyrope entropy changed to 267·0 J/K (was 269·5)