Silicic magma systems are of great scientific interest and societal importance owing to their role in the evolution of the crust and the hazards posed by volcanic eruptions. MELTS is a powerful and widely used tool to study the evolution of magmatic systems over a wide spectrum of compositions and conditions. However, the current calibration of MELTS fails to correctly predict the position of the quartz + feldspar saturation surface in temperature, pressure and composition space, making it unsuitable to study silicic systems. We create a modified calibration of MELTS optimized for silicic systems, dubbed rhyolite-MELTS, using early erupted Bishop pumice as a reference. Small adjustments to the calorimetrically determined enthalpy of formation of quartz and of the potassium end-member of alkali feldspar in the MELTS calibration lead to much improved predictions of the quartz + feldspar saturation surface as a function of pressure. Application of rhyolite-MELTS to the Highland Range Volcanic Sequence (Nevada), the Peach Spring Tuff (Arizona–Nevada–California), and the late-erupted Bishop Tuff (California), using compositions that vary from trachydacite to high-silica rhyolite, shows that the calibration is appropriate for a variety of fluid-bearing silicic systems. Some key observations include the following. (1) The simulated evolutionary paths are consistent with petrographic observations and glass compositions; further work is needed to compare predicted and observed mineral compositions. (2) The nearly invariant nature of silicic magmas is well captured by rhyolite-MELTS; unusual behavior is observed after extensive pseudo-invariant crystallization, suggesting that the new calibration works best for relatively small (i.e. <50 wt %) crystallization intervals, comparable with what is observed in volcanic rocks. (3) Our success with rhyolite-MELTS shows that water-bearing systems in which hydrous phases do not play a critical role can be appropriately handled; simulations are sensitive to initial water concentration, and although only a pure-H2O fluid is modeled, suitable amounts of water can be added or subtracted to mimic the effect of CO2 in fluid solubility. Our continuing work on natural systems shows that rhyolite-MELTS is very useful in constraining crystallization conditions, and is particularly well suited to explore the eruptive potential of silicic magmas. We show that constraints placed by rhyolite-MELTS simulations using late-erupted Bishop Tuff whole-rock and melt inclusion compositions are inconsistent with a vertically stratified magma body.
Silicic magma systems are the source of fascinating magmatic phenomena. They also play an important role in the origin and evolution of the crust and their eruptions pose significant hazards to society. One of the challenges in understanding these systems is the fact that many of the relevant processes occur at depth, and on timescales that far exceed human lifetimes. Unable to directly observe these phenomena, petrologists resort to a combination of observations on active systems, study of the record preserved in rocks, and modeling of relevant processes.
One powerful modeling tool, often used to understand the chemical and phase assemblage evolution of magmatic systems, is the thermodynamic modeling software MELTS (Ghiorso & Sack, 1995; Asimow & Ghiorso, 1998). MELTS has been extensively applied to understand the evolution of relatively dry mafic–intermediate magmas (e.g. Ghiorso, 1997; Gaetani et al., 1998; among many others). Application to wet intermediate magma compositions is limited by the lack of appropriate thermodynamic models for hydrous mafic silicates, particularly amphibole and biotite. Inasmuch as the evolution of silicic magmas is dominated by felsic minerals and a fluid phase, there is no major limitation for the application of MELTS to silicic systems. However, a very limited number of relevant experiments on silicic systems were available at the time of calibration of MELTS, and the current calibration fails to correctly predict important features of these systems. Here, we present a modified calibration of MELTS, dubbed rhyolite-MELTS, optimized for fluid-bearing silicic systems. The fundamental goal of the modification is for rhyolite-MELTS to correctly predict the quartz + feldspar saturation surface in temperature, pressure, and composition space. In addition to a modified calibration, rhyolite-MELTS employs a more efficient algorithm for detecting feldspar and pyroxene phase saturation, and implements the rhombohedral oxide solution model of Ghiorso & Evans (2008).
Our calibration is based on compositions determined and crystallization conditions inferred for early erupted Bishop Tuff [California (CA); Bailey et al., 1976] magma. Importantly, rhyolite-MELTS differs from MELTS only in the stability of quartz and sanidine, with the result that rhyolite-MELTS inherits the virtues and limitations of MELTS in other respects. We test the quality of the new calibration by comparing predictions by MELTS and rhyolite-MELTS for the quartz + feldspar saturation surface using experimental data available in the Library of Experimental Phase Relations (LEPR; Hirschmann et al., 2008). We also test the effectiveness of the new calibration in retrieving the position of the quartz + feldspar saturation surface using examples from the Highland Range Volcanic Sequence [Nevada (NV); Faulds et al., 2002; Colombini et al., 2011] and the Peach Spring Tuff [NV–Arizona (AZ)–CA); Young & Brennan, 1974], with compositions ranging from trachydacite to high-silica rhyolites with distinct compositional traits. We show that rhyolite-MELTS yields much improved results for the position of the quartz + feldspar saturation surface as a function of pressure. Our applications highlight how rhyolite-MELTS can be used to constrain crystallization conditions, and how it can be used to explore differences in behavior owing to varying initial fluid contents and timing of fluid saturation. We conclude by showing how rhyolite-MELTS can be used to retrieve significant information on the heat of cooling and crystallization and the crystallization conditions of the Bishop Tuff magma body.
STATEMENT OF THE PROBLEM
This contribution emerged from our desire to apply MELTS to evolved silicic systems, in particular rhyolitic magmas. Simulations for silicic systems using the current calibration of MELTS are affected by significant limitations.
We take the early erupted Bishop Tuff as an example (Fig. 1). We use as the initial composition the average early erupted bulk-pumice composition of Hildreth (1979). MELTS calculations assume water saturation at 200 MPa [which yields 5–6 wt % dissolved water, in agreement with the estimates of Wallace et al. (1999)], with oxygen fugacity buffered at Ni–NiO, and are performed under equilibrium conditions. Sanidine is the liquidus phase at >800°C, with magnetite-, plagioclase- and quartz-in within the temperature interval 725–715°C (Fig. 1a).
The problems with these calculations are well portrayed in projections onto the quartz–albite–orthoclase (Qz–Ab–Or) ternary. We follow the projection scheme of Blundy & Cashman (2001) to project anorthite-bearing compositions onto the ternary, a scheme that is suitable for metaluminous compositions (i.e. normative corundum < 1, normative acmite = 0). The advantage of using this projection is that natural and calculated compositions can be readily compared with the experimentally determined quartz + sanidine cotectic curves, known for a variety of relevant pressure conditions (for details, see Blundy & Cashman, 2001). In our early erupted Bishop Tuff example (Fig. 1), given that the initial composition lies very close to the 200 MPa cotectic (the pressure used in the simulation), we would expect quartz and sanidine to saturate at very similar temperatures, which is not observed. Instead, simulated compositions (filled circles in Fig. 1b) move away from the Ab–Or join, leading to unrealistic melt compositions with >82 wt % SiO2 (anhydrous basis).
We argue that the stability field of quartz is underestimated, and this is a systematic problem with the current calibration of MELTS. The root of the problem is that the liquid regular solution model used by MELTS is unable to appropriately describe the activity of SiO2 in the melt phase for silicic melt compositions. Correction of the problem would require adopting non-regular solution models for the melt. The large temperature interval between sanidine and quartz saturation may also be due to an excessively large stability field for alkali feldspar. It is noteworthy that, in their study of the Campanian Ignimbrite, Fowler et al. (2007) found that K2O values in the calculated melts are systematically underpredicted, suggesting that the stability field of sanidine may be overpredicted by MELTS. At the time the MELTS calibration was developed there were no available experiments with which to calibrate the potassium end-member of alkali feldspar, and calibration was performed in an internally consistent fashion using experimental results for leucite. However, the data were assumed to refer to pure leucite (KAlSi2O6), when actual leucite has significant amounts of the sodium end-member (up to 10 mol % NaAlSi2O6). As a result, not only is the stability field of sanidine potentially overestimated, but also, the calculated compositions are likely to be excessively potassic.
The purpose of this contribution is to remediate these two problems and make the use of MELTS for silicic systems practical. Since the original calibration of MELTS (Ghiorso & Sack, 1995), significant advances have been made towards developing more sophisticated solution models for the melt phase (e.g. Ghiorso, 2004), and more experiments now exist involving the potassium feldspar end-member of alkali feldspar (Holtz et al., 2001; Koester et al., 2002; Scaillet & Macdonald, 2003; Patiño-Douce, 2005). However, any modifications involving a new model for the melt phase, as well as the incorporation of new experiments into the reference experimental database, would require recalibration of MELTS over the entire spectrum of natural compositions. Such an effort would be hampered by the relative scarcity of relevant experiments for intermediate compositions, and it is not clear that it would be possible to obtain a calibration that is superior to the existing one.
Given these difficulties, we take a different approach. Rather than recalibrating MELTS over the entire spectrum of compositions, we search for small corrections to the stability of quartz and the potassium end-member of alkali feldspar that lead to more appropriate calculation of the quartz + feldspar saturation surface for silicic systems. Although the modified calibration we propose has some limitations, we argue here that it leads to much improved results and can be used to gain significant insight into the evolution of silicic magmas. Importantly, no modification is made to any of the other phases in the MELTS database, so that predicted phase relations other than those involving quartz and alkali-feldspar are unaffected.
THE MODIFIED CALIBRATION
We use early erupted Bishop Tuff pumice as a reference in searching for a modified calibration. The similarity in composition of bulk-pumice and melt inclusions (Fig. 1b) highlights the nearly invariant (or eutectoid) nature of early erupted Bishop magma (Anderson et al., 2000), and suggests crystallization along the quartz + sanidine + plagioclase saturation surface, consistent with the observed phenocryst assemblage. These findings are in complete agreement with previous studies. The homogeneity of phenocrysts in major element composition (Hildreth, 1979) implies nearly isothermal, closed-system crystallization. Melt inclusion trace and volatile compositions are consistent with closed-system crystallization at pressures of ∼175 MPa under fluid-saturated conditions (∼5·5 wt % H2O, ∼100 ppm CO2; Wallace et al., 1995, 1999; Anderson et al., 2000), consistent with the location of bulk-pumice and melt inclusion compositions projected onto the Q–Ab–Or ternary (Fig. 1b).
Accordingly, we use an average melt inclusion composition (Table 1) calculated from the data of Anderson et al. (2000), which MELTS calculates to be water-saturated at 175 MPa, under Ni–NiO oxygen fugacity conditions. We find that simultaneous saturation in quartz, sanidine and plagioclase can be achieved by modest changes in the enthalpies of formation of quartz and the potassium end-member of alkali feldspar; the required changes are within or of the same order as the uncertainties associated with fitting the solution model used by MELTS to experimental data. Specifically, the average error in fitting the MELTS solution model is ∼2 kJ mol−1, and the required changes of −1·3 kJ mol−1 SiO2 for quartz and 3·4 kJ mol−1 KAlSi3O8 for alkali feldspar represent small amounts of residual energy. Importantly, the changes made are in the direction we would expect from the considerations above; that is, a negative change in enthalpy leads to larger stability fields for quartz, whereas a positive change for alkali feldspar leads to smaller fields of stability for sanidine.
|Early erupted Bishop Tuff||Late-erupted Bishop Tuff||Highland Range (HRL 21)||Highland Range (HRL 27)||Peach Spring Tuff (KPST 01D)|
|Early erupted Bishop Tuff||Late-erupted Bishop Tuff||Highland Range (HRL 21)||Highland Range (HRL 27)||Peach Spring Tuff (KPST 01D)|
Using this modified calibration we find that at 175 MPa our bulk composition is saturated in quartz, sanidine, plagioclase and fluid at 759·9°C, with magnetite joining the crystallizing assemblage at ∼758·6°C (Fig. 2a). Crystallization is essentially invariant, with almost 30 wt % crystallization in the first 1·5°C of cooling, followed by an additional 60 wt % crystallization in the following 1·7°C, with a total of 90 wt % crystallization in slightly more than 3°C cooling. As expected, the calculated composition of the melt phase is effectively constant (Fig. 2d–g). Feldspar compositions vary over a narrow range of compositions (sanidine: Ab49–45Or49–54An2–1; plagioclase: Ab76–78Or17–13An7–9; see Fig. 2b and c), in agreement with the uniformity of phenocryst compositions in Bishop pumice (Hildreth, 1979). Feldspar compositions compare reasonably well with observed compositions (sanidine: Ab36Or63An1; plagioclase: Ab80Or14An6; Hildreth, 1979), particularly for plagioclase, although observed sanidine compositions are somewhat more potassic than calculated ones. The absolute crystallization temperatures are 30–40°C higher than the temperatures of ∼720–730°C estimated by Hildreth (1979) based on Fe–Ti oxides [for a discussion of the potential role of post-eruptive re-equilibration of Fe–Ti oxides, see Ghiorso & Evans (2008)] and experimental evidence (Scaillet & Hildreth, 2001), which could explain much of the discrepancy between observed and calculated sanidine compositions.
Inspection of our results (Figs 1b and 2a) shows that after substantial invariant crystallization (e.g. >50 wt %), the rate of crystallization with temperature starts to decline; this is an issue that we observe in all of our simulations, both with MELTS and rhyolite-MELTS. This behavior is unexpected, and we believe it to be an artifact of the calculations. We attribute this problem to inaccuracies in the calculated phase proportions or phase compositions, which eventually cause the calculated melt compositions to drift off the invariant; it follows that calculated sub-solidus assemblages are also inaccurate. We thus expect rhyolite-MELTS to be most successful when applied to volcanic systems, in which crystal contents are well below 50 wt %.
Although sanidine compositions, absolute temperatures, and sub-solidus phase assemblages calculated by our modified calibration of MELTS are potentially somewhat inaccurate, we are encouraged by its ability to retrieve the phase assemblage and phase properties of an invariant magma such as the early erupted Bishop Tuff. The success of our modified calibration will be, however, ultimately measured by how successful rhyolite-MELTS is in providing useful information for the evolution of silicic magmas in general. In the following sections, we compare the two calibrations with available experiments. We then discuss initial applications of rhyolite-MELTS aimed at (1) testing the quality of the proposed modified calibration, and (2) highlighting the potential brought by the application of rhyolite-MELTS to silicic volcanic systems.
COMPARISON WITH EXPERIMENTAL DATA
A relatively small number of experiments are now available for silicic systems that can be used to test our modified calibration. We extracted from LEPR (Hirschmann et al., 2008) all experiments that met the following conditions: (1) pressure <1 GPa; (2) SiO2 >55 wt % (anhydrous basis); (3) quartz-saturated; (4) anhydrous or with H2O as the only volatile species (no CO2-bearing experiments available; Cl- and/or F-bearing experiments excluded); (5) bulk composition not within the Na2O–K2O–SiO2 ternary. The last condition is necessary given that MELTS was designed to take advantage of the simplifications that arise in describing the energetics of mixing in complex chemical systems such as natural silicate melts (Ghiorso et al., 1983); as such, MELTS is ill-suited for extrapolation to simpler systems (e.g. pure silica liquids, binary and ternary systems), so we limit testing and applications to compositions with at least five chemical components (e.g. CaO–Na2O–K2O–Al2O3–SiO2).
A total of 89 experiments in LEPR meet all five conditions above, and they can be used to test rhyolite-MELTS; 84 of them are saturated in at least one feldspar. Of these, 45 are water-saturated, all of which are also saturated in at least one feldspar. Only two experiments are simultaneously saturated in quartz, plagioclase, alkali feldspar, and water, which underscores the difficulty in calibrating MELTS for silicic systems, and justifies our choice of using a well-studied natural example for calibration. Most of the silica-rich, water-undersaturated experiments correspond to pelite melting experiments that yield very low melt fractions (e.g. Patiño-Douce, 2005); water contents in the melt are typically calculated by difference from electron microprobe glass analyses; as such these experiments are much less reliable than the water-saturated ones, and we focus on the water-saturated experiments in our discussion below.
For the 89 experiments, we calculate the chemical potential of silica in the melt phase and in quartz under the experimental conditions; we refer to the difference in these two chemical potentials as the residual. Because all selected experiments are quartz-saturated, all residuals should be zero. Comparison of the residuals calculated using MELTS and rhyolite-MELTS (Fig. 3) shows that MELTS yields positive residuals for most experiments, particularly for water-saturated experiments (Fig. 3a and b); positive residuals correspond to quartz-undersaturation, consistent with the lack of quartz in our MELTS simulations for early erupted Bishop Tuff compositions. Residuals for rhyolite-MELTS (Fig. 3c and d), on the other hand, scatter around zero, as expected for quartz saturation. Importantly, the experiments span the compositional spectrum from dacite to high-silica rhyolite; no obvious effect of silica content on residuals is seen for the water-saturated experiments. This is a remarkable result, given that the experiments were not used in the calibration; the much better residual distribution achieved using rhyolite-MELTS strongly suggests that the new calibration substantially improves calculation of the quartz + feldspar saturation surface.
APPLICATIONS TO SILICIC SYSTEMS
Highland Range Volcanic Sequence
The Colorado River Extensional Corridor (southern Nevada–northwestern Arizona) is characterized by large crustal blocks tilted during the opening of the Rio Grande Rift, which expose extensive vertical sections of the crust in a region where magmatism played a major role at 15–19 Ma (Faulds et al., 2001). In these blocks, thick sections of granitic intrusions are exquisitely exposed at the surface, along with potentially correlative volcanic sequences (Miller & Miller, 2002). Age correlations, geochemical similarities, and geographical proximity between the upper, more silicic portion of the Highland Range sequence (Colombini, 2009; Colombini et al., 2011) and the Searchlight Pluton (Bachl et al., 2001; Miller & Miller, 2002) suggest that the Highland Range is the volcanic counterpart of the Searchlight plumbing system (Faulds et al., 2002). The rock record in the Highland Range sequence reveals that the compositional evolution of the system from trachydacite to high-silica rhyolite was accompanied by a shift in eruptive style, from effusive to explosive (Colombini, 2009; Colombini et al., 2011).
We model the evolution of two magma compositions from the Highland Range using rhyolite-MELTS (following Vaum et al., 2009), a high-silica rhyolite (sample HRL 21) and a trachydacite (HRL 27), for which we have whole-rock and glass compositions (Table 1). Further details of these samples have been given by Colombini et al. (2011).
The case for the high-silica rhyolite (HRL 21) has much in common with the early erupted Bishop Tuff example presented above. Although the crystallization pressure is only relatively constrained (less than ∼200 MPa, based on correlations with the nearby Searchlight Pluton; see Bachl et al., 2001; Miller & Miller, 2002), the relative position of whole-rock and glass compositions in the Qz–Ab–Or ternary (Fig. 4a) suggests cotectic crystallization at ∼150 MPa. Consistent with the considerations above, MELTS calculations show that quartz saturation is underestimated and calculated melt compositions move away from the cotectic towards the quartz apex; rhyolite-MELTS calculations, on the other hand, correctly predict quartz + sanidine saturation, with observed glass compositions corresponding to pseudo-invariant melt compositions attained upon plagioclase saturation, in accordance with the coexistence of quartz, sanidine and plagioclase in the actual rock.
The bulk composition of the trachydacite (HRL 27) falls well within the quartz-undersaturated field (Fig. 4b). We expect crystallization conditions to be similar to those of HRL 21, so we simulate crystallization at 150 MPa and variable amounts of water. Plagioclase is the first felsic phase in both MELTS and rhyolite-MELTS simulations, which drives melt compositions away from the Ab apex in the ternary and toward the glass compositions (Fig. 4b). The next felsic phase to crystallize is sanidine. In MELTS, sanidine becomes stable in less evolved melt compositions, which causes the melt phase composition to evolve toward the Qz apex, bypassing the observed glass composition. In the rhyolite-MELTS simulations, sanidine becomes stable only in more evolved melt compositions, causing these to approach the observed glass composition much more closely. It is also apparent that once again MELTS underpredicts quartz saturation, causing the melt phase to evolve toward unreasonably high silica concentrations (Fig. 4b); the same does not happen with rhyolite-MELTS, which calculates nearly invariant crystallization upon reaching quartz saturation for melt compositions falling between the 100 and 200 MPa cotectic curves, as expected for crystallization at 150 MPa. The observed phenocryst assemblage, consisting of plagioclase + sanidine with no quartz, is consistent with the evolution predicted by rhyolite-MELTS.
In combination, these simulations reinforce the notion that quartz is underpredicted and sanidine is overpredicted by MELTS, lending confidence to the corrections suggested based on early erupted Bishop Tuff compositions. Rhyolite-MELTS is not only more successful in predicting melt evolution paths consistent with observed phenocryst assemblages and glass compositions, but it also seems to much better predict the quartz + feldspar saturation surface as a function of pressure.
Peach Spring Tuff
The Peach Spring Tuff (southwestern USA; Young & Brennan, 1974) is a giant ignimbrite (>500 km3), exposed over an area of >35 000 km2 in SE California, NW Arizona, and southern Nevada (Buesch & Valentine, 1986; Glazner et al., 1986), formed by a super-eruption at 18·8 Ma (Nielson et al., 1990; Ferguson et al., in preparation). The Peach Spring Tuff has a distinctive phenocryst assemblage, dominated by sanidine with little to no quartz and relatively abundant titanite (Glazner et al., 1986; Pamukcu et al., in preparation). Recent work suggests that the Peach Spring Tuff erupted from the Silver Creek Caldera, Arizona (Ferguson, 2008; Ferguson et al., in preparation); intracaldera deposits reveal a history of heating and mush remobilization that contrasts with a record of monotonic cooling and crystallization preserved in the high-silica rhyolites typical of the outflow deposits (Pamukcu et al., in preparation).
Our simulations use a high-silica rhyolite composition representative of the outflow sheet exposed in Kingman, Arizona (Table 1; for sample description and location, see Gualda et al., 2010), which corresponds to one of the most evolved compositions in the deposit (Carley, 2010; Pamukcu et al., in preparation). As in the case of the Highland Range, the crystallization pressure is poorly constrained; however, inspection of the Qz–Ab–Or ternary (Fig. 5) suggests crystallization at <250 MPa. As part of a study of the eruptive potential of giant magma bodies (Carley, 2010; Carley et al., in preparation), we are conducting simulations over a broad range of pressure (150–350 MPa) and water contents (2–7 wt % H2O); only representative results aimed at testing the new calibration are shown in Fig. 5
It is known from experiments that the quartz + feldspar saturation surface migrates towards more quartz-rich compositions with decreasing pressure (for reviews see Johannes & Holtz, 1996; Blundy & Cashman, 2001). Our simulations using the Peach Spring high-silica rhyolite composition show that rhyolite-MELTS captures well the expected change in quartz saturation with pressure (Fig. 5b). There is good agreement between the positions of the quartz + sanidine saturation surface calculated by rhyolite-MELTS and the experimental cotectic curves (as sketched by Blundy & Cashman, 2001) for lower pressures (150–200 MPa). The agreement is not nearly as good for higher pressures (300–350 MPa), with rhyolite-MELTS suggesting smaller changes in the position of the saturation surface with increasing pressure. Interestingly, at higher pressure (e.g. 350 MPa), quartz saturates in less evolved melts, with the effect that the melt phase has a more extended evolution along the quartz + feldspar saturation surface (Fig. 5b).
Importantly, the quartz + feldspar saturation surface is consistently calculated irrespective of initial dissolved water content (Fig. 5c). However, the temperature of feldspar saturation is affected by the initial water content, and, accordingly, calculated feldspar compositions vary significantly as a function of initial water contents. For instance, for 2 wt % H2O, sanidine saturates with composition Or58Ab40 at 892°C, and plagioclase saturates with composition Ab70An20 at 805°C; in contrast, for 7 wt % H2O, sanidine appears at 754°C and is more potassic (Or79Ab21), whereas plagioclase starts crystallizing at 735°C but is similar in composition (Ab70An23).
It is noteworthy that the saturation sequence suggested by our rhyolite-MELTS simulations (sanidine → quartz → plagioclase) is not what would be expected from the phenocryst assemblage (sanidine > plagioclase > quartz). Although the changes in the stability of quartz and sanidine implemented in rhyolite-MELTS are justified based on the limitations of MELTS and its calibration, and are well constrained by the early Bishop Tuff example, we do not find similar grounds for a modification of plagioclase stability. Nonetheless, what our data show is that even with this potential limitation, rhyolite-MELTS successfully calculates the quartz + feldspar saturation surface as a function of pressure for a variety of water contents, which makes it significantly superior to MELTS, and opens new possibilities for the study of silicic systems.
Late-erupted Bishop Tuff
Since the seminal work by Hildreth (1979), the Bishop Tuff (Bailey et al., 1976; Hildreth, 2004) has been considered an archetypical example of a magma body vertically stratified in both temperature and composition (see also Hildreth & Wilson, 2007). Evidence for thermal stratification stems from Fe–Ti oxide (Hildreth, 1979; Hildreth & Wilson, 2007), oxygen isotope (Bindeman & Valley, 2002) and Ti-in-quartz thermometry (Wark et al., 2007); accompanying compositional zonation is inferred from phenocryst (Hildreth, 1979) and melt inclusion (Wallace et al., 1999) compositions; vertical stratification has been implied from the stratigraphy of the deposit (Wilson & Hildreth, 1997; Hildreth & Wilson, 2007), with some support from melt inclusion H2O–CO2 barometry (Wallace et al., 1999). In particular, Wallace et al. (1999) favored crystallization at ∼250 MPa for late-erupted melt inclusions, in contrast to crystallization at ∼175 MPa for early erupted melt inclusions.
To test rhyolite-MELTS, we performed simulations under a variety of pressures (100–350 MPa) with both bulk-pumice and melt inclusion compositions representative of the late-erupted Bishop Tuff (Table 1). Application of the H2O–CO2 fluid solubility model of Papale et al. (2006) shows that the combination of relatively low H2O contents (3·7 wt %) and high CO2 (450 ppm) results in fluid saturation for the conditions of interest, reinforcing the original conclusions of Wallace et al. (1999). Accordingly, all simulations shown here are conducted under water-saturated conditions, which are achieved by using H2O concentrations in excess of those measured. In effect, our approach causes water activity in the fluid phase to be kept constant—as expected for fluid-saturated conditions—but overestimated. Experimental evidence (see Johannes & Holtz, 1996) indicates that saturation temperatures are affected by reduced water activity in the fluid; for XH2O = 0·6 calculated using the Papale et al. (2006) model, temperatures would be underestimated by ∼40°C. The effect on the melt (and crystallizing assemblage) is expected to be small, however, given the extremely low solubility of CO2 in the melt (e.g. Tamic et al., 2001). This is confirmed by experiments (e.g. Johannes & Holtz, 1996), which show that phase compositions and proportions—including the position of the cotectic surfaces and the ternary minimum—are not affected, lending support to our approach.
The projection of the late-erupted bulk-pumice composition of Hildreth (1979; Table 1) into the Qz–Ab–Or ternary (Fig. 6) plots close to the 300 MPa cotectic; this is consistent with initial crystallization in the feldspar field, under pressures of 250 MPa (or lower), with quartz saturation only after some amount of crystallization. Interestingly, melt inclusions (Anderson et al., 2000; Table 1) are somewhat more silicic than the bulk-pumice, and they plot along cotectic curves of ≤200 MPa pressure (Fig. 6).
We use rhyolite-MELTS simulations for a late-erupted average melt inclusion composition to explore the variation in crystallization trajectories as a function of pressure, and we seek conditions that lead the initial bulk-pumice composition toward the melt inclusion compositions. Our simulations show once again that rhyolite-MELTS captures well the effect of pressure on the quartz + feldspar saturation surface (Fig. 6b), in a manner consistent with the projection of Blundy & Cashman (2001). For 250 MPa and above, the liquidus phase is quartz, which drives the melt composition toward the Ab–Or join, until sanidine saturation is reached, from which point the melt composition evolves along a path subparallel to the cotectic curves; nearly invariant crystallization is observed upon plagioclase saturation. For 100 MPa, crystallization starts in the feldspar field, with sanidine as the liquidus phase, until quartz saturation—and nearly invariant crystallization—is reached very close to the trace of the 100 MPa cotectic. For 175 MPa, quartz and sanidine appear at ∼762°C, and plagioclase appears at only 5°C lower, at 757°C; simultaneous crystallization of sanidine, quartz and plagioclase in the presence of a fluid phase is consistent with the observed phenocryst assemblage and measured melt inclusion volatile contents. Just as in the Peach Spring Tuff case, the agreement between the position of the quartz + feldspar saturation surface calculated by rhyolite-MELTS and the traces of the cotectic curves on the Qz–Ab–Or ternary is excellent for low pressure (i.e. <200 MPa), but not for the 300 MPa cotectic, suggesting that the 300 MPa cotectic curve lies too far into the feldspar field. Surprisingly, our simulations suggest fairly strongly that melt inclusion compositions are not consistent with crystallization under relatively high pressures (i.e. close to 250 MPa), but favor instead crystallization under ∼175 MPa.
Simulations using the bulk late-erupted pumice composition reveal a strong dependence of the crystallization path on pressure (Fig. 6c). For 250 and 350 MPa, sanidine is the first felsic phase to appear, followed first by quartz and then plagioclase; for 175 MPa, sanidine is followed by plagioclase and then quartz; and for 100 MPa, plagioclase is the first felsic phase, followed by sanidine then quartz. This results in contrasted evolution paths: for the higher pressures (250 and 350 MPa), melt compositions evolve toward the Qz–Ab join, away from the observed melt inclusions; in contrast, for lower pressures (100 and 175 MPa), melt compositions evolve toward the Qz–Or join, toward the melt inclusion compositions (Fig. 6c). In this sense, our simulations using bulk-pumice compositions also suggest crystallization under pressures lower than ∼200 MPa, consistent with our simulations with melt inclusion compositions. The same conclusion is drawn based on the comparison of calculated melt compositions with measured melt inclusion compositions (Fig. 7), as particularly well shown in the Or–Ab–An ternary (Fig. 8a) and in select binary diagrams of major oxides (Fig. 8b). As expected, silica concentrations in the melt are constrained by the pressure of crystallization, with increasing SiO2 values with decreasing pressure (Fig. 8c). Although the calculated compositions are not as silicic as the observed melt inclusion compositions, calculated and measured compositions differ by less than 3% (<1% for 100 and 175 MPa, ∼2% for 250 MPa, ∼3% for 350 MPa).
Interestingly, the temperature of nearly invariant crystallization in simulations using bulk-pumice compositions varies significantly with pressure (781°C at 100 MPa, 756°C at 175 MPa, 743°C at 250 MPa, 729°C at 350 MPa; Fig. 8d); this is in fact a well-known property of hydrous magmas that derives from the negative Clapeyron slope of the fluid-saturated solidus. This would imply that vertical stratification of water-saturated, nearly invariant magmas would lead to a modest (∼15°C) inverse temperature gradient; the effect of reduced water saturation would be sufficient to offset this negative gradient if a gradient in fluid composition existed, resulting in a possible modest positive thermal gradient (∼25°C using average melt inclusion compositions). The onset of co-crystallization of magnetite and ilmenite takes place at temperatures fairly close to or at the invariant (invariant at both 100 and 175 MPa, 5°C above invariant at 250 MPa, 30°C above invariant at 350 MPa). In this sense, our simulations are largely inconsistent with the high temperatures retrieved from Fe–Ti oxides for the late-erupted Bishop pumice (Hildreth, 1979; Hildreth & Wilson, 2007), unless very low pressures (<150 MPa) prevailed during crystallization. The much smaller temperature gradients inferred from Ti concentrations in quartz crystal interiors (Wark et al., 2007) can be reconciled with our simulations.
Surprisingly, our simulations are in many respects inconsistent with the prevailing view of a vertically stratified magma body for the Bishop Tuff. Crystallization of late-erupted Bishop pumice at 250 MPa is inconsistent with the melt inclusion compositions, which themselves suggest crystallization at ∼175 MPa, a pressure that is identical to that inferred for the early erupted Bishop magma. We note that whereas Wallace et al. (1999) interpreted melt inclusion compositions to reflect crystallization under relatively high pressures (i.e. close to 250 MPa), single estimates span a wide pressure interval between 110 and 280 MPa (Fig. 9a), displaying a unimodal distribution with average of 178 MPa (see also Anderson et al., 2000). The temperature record also seems to be problematic, in that vertical stratification of fluid-saturated, nearly invariant magmas would lead to modest (<30°C) temperature gradients, much smaller than what is retrieved from Fe–Ti oxide thermometry, but consistent with the record of Ti-in-quartz (Wark et al., 2007). We note that the Fe–Ti oxide thermometry record has been called into question before, in particular because of the lack of agreement between temperatures retrieved from Fe–Ti and Fe–Mg equilibria (Frost & Lindsley, 1992; Ghiorso & Evans, 2008). Finally, the advocated continuum in composition between early and late-erupted compositions is not entirely borne out by melt inclusion major, trace, and volatile elements (Fig. 9b–d) or by phenocryst major and trace element compositions (Fig. 10), which reveal bimodal distributions suggestive of two discrete magma compositions. In this context, it seems relevant to note that the early erupted and late-erupted pumice are inferred to have erupted from separate vents (Hildreth & Mahood, 1986; Wilson & Hildreth, 1997), potentially tapping lateral, rather than vertical, variations. All in all, derivation of Bishop pumice from two discrete, laterally juxtaposed magma bodies of rhyolitic magma is favored by phase equilibria considerations using rhyolite-MELTS, and can be reconciled with much (if not all) of the existing evidence. We recognize that our analysis is somewhat limited in scope, and may be affected by sampling bias; but it does cast significant doubt on the widely held view of an archetypical vertically stratified magma body for the Bishop Tuff.
Application of the phase equilibria calculation software MELTS to silicic magma compositions shows that its current calibration fails to capture important first-order features of these systems; specifically, the stability field of quartz is underestimated, whereas the stability field of sanidine is overestimated.
We have developed a new calibration, dubbed rhyolite-MELTS, based on the well-constrained crystallization conditions and the ‘eutectoid’ nature of the early erupted Bishop Tuff magma. Small changes in the enthalpies of formation of quartz and the potassium end-member of alkali feldspar lead to nearly invariant crystallization of an average, early erupted melt inclusion composition under 175 MPa pressure and water-saturated conditions.
Comparison of data for the early erupted Bishop pumice reveals potential issues with the calibration, as follows.
Absolute temperatures of ∼760°C calculated by rhyolite-MELTS for nearly invariant crystallization are 30–40°C higher than estimated from thermometry and experiments; we expect this effect to be systematic.
The agreement between calculated and observed mineral compositions is adequate, but calculated sanidine compositions are more sodic than observed ones; this can, at least in part, be explained by the bias in absolute temperature.
After extensive nearly invariant crystallization, minor inaccuracies in the calculated phase proportions or compositions cause the calculated melt composition to drift away from the invariant; it follows that rhyolite-MELTS is best suited for volcanic systems with <50 wt % crystals.
Simulations using magma compositions representative of the Highland Range Volcanic Sequence (NV), the Peach Spring Tuff (AZ–NV–CA), and the late-erupted Bishop Tuff (CA) show the following.
Rhyolite-MELTS correctly calculates the position of the quartz + feldspar saturation surface and reveals nearly invariant behavior for systems with initial compositions ranging from dacite to high-silica rhyolite.
The new calibration captures well the dependence of the position of the quartz + feldspar saturation surface on pressure; specifically, rhyolite-MELTS correctly calculates the shift toward more silicic compositions with decreasing pressure.
Simulations using bulk-rock and glass (matrix or melt inclusion) compositions can be very useful in constraining the crystallization conditions of silicic systems, in particular crystallization pressure. This is particularly well shown by our simulations using late-erupted Bishop Tuff compositions, which are inconsistent with crystallization at pressures >200 MPa, casting doubt on the vertical stratification of the Bishop Tuff magma body.
Rhyolite-MELTS is available for public use and can be downloaded from: http://melts.ofm-research.org/.
Material support for this investigation was provided by the National Science Foundation grants to G.A.R.G. (EAR-0948528) and M.S.G. (EAR-0948734). Support to R.V.L. was provided by the Vanderbilt Undergraduate Summer Research Program. This work was partly supported by the National Science Foundation through a Graduate Research Fellowship to T.L.C.
Discussions with the MESSY group at Vanderbilt helped improve the paper. Reviews by James Gardner, Claudia Cannatelli and an anonymous reviewer are greatly appreciated.