Diverse explanations exist for the large-volume catastrophic eruptions that formed the Bishop Tuff of Long Valley in eastern California, the Bandelier Tuff in New Mexico, and the tuffs of Yellowstone, Montana, USA. These eruptions are among the largest on Earth within the last 2 Myr. A common factor in recently proposed petrogenetic scenarios for each system is multistage processing, in which a crystal mush forms by crystal fractionation and is then remobilized to liberate high-silica liquids. Magma evolves in the lower crust in earlier phases. We have tested these scenarios quantitatively by performing phase equilibria calculations (MELTS) and comparing the results with observed liquid (glass) and phenocryst compositions. Although comparison of tuff samples from each ignimbrite reveals distinct phenocryst compositions and proportions, the computed results exhibit a remarkable degree of congruity among the systems, pointing to some underlying uniform behavior relevant to large-volume silicic ignimbrites. Computed liquid compositions derived from more than ∼25% fractional crystallization of the parental melt in the deep crust are marked by SiO2 concentrations several weight per cent too low compared with the observed compositions, suggesting a limit on the extent of magma evolution by crystal fractionation in the deep crust. In all cases, the phase equilibria results and related considerations point to evolution dominated by crystal fractionation of a water-saturated mafic parental melt at shallow depths (∼5 km). Parental melt compositions are consistent with those of observed regional primitive basalts erupted prior to ignimbrite eruption for each system in each region. Fractional crystallization of water-rich mafic melt at shallow levels leads inherently to destabilization near thermodynamic pseudoinvariant points at around 800°C within the melting interval close to, but above, the solidus. For each system, the magmas evolve to states of high exsolved H2O volume fraction even at 5 km depth, eventually exceeding the criterion for magma fragmentation of ∼60 vol. % near the pseudoinvariant point temperature. Copious exsolution and possible expulsion of fluid occurs at this temperature, where the solid fraction in the magma changes almost discontinuously (isothermally) to significantly higher values. This instability mechanism acts as an eruption trigger by generating a gravitationally unstable arrangement of low-density, water-saturated magma beneath a thin (several kilometres) crustal lid. The trigger mechanism is common to fractional crystallization scenarios based on a variety of conditions, when crystallized solids and/or exsolved fluids are fractionated from residual melt isobarically (constant pressure) or isochorically (constant volume). In a single system, differences in liquid compositions resulting from constant volume versus constant pressure crystallization and expulsion versus retention of exsolved H2O are small compared with those arising from variations in initial water concentration, lithostatic pressure, and oxygen fugacity. It is these latter quantities that lie at the crux of the commonality in large-volume ignimbrite-forming eruptions, with a reasonable range of metamodel parameters. Scale analysis provides thermal timescales for fractional crystallization, including age ranges for discrete crystal populations. For the Bishop Tuff, the overall timescale for the Bishop magma body is >1 Myr. For the Yellowstone Tuffs, calculated thermal timescales are consistent with recurrence intervals of ∼600 kyr between successive caldera collapses. Although it is recognized that petrogenetic processes other than perfect fractional crystallization play a role in ignimbrite petrogenesis, by emphasizing common features the uniqueness of each system can be brought into better focus by sound and quantitative analysis.
The mean age of the Earth’s continental crust is ∼2 Ga. Crustal growth has occurred throughout geological time, perhaps episodically (Taylor & McLellan, 1985; Condie, 2000; Hong et al., 2004), by lateral accretion along destructive plate margins, and by mafic magma underplating of continental crust (e.g. formation of large igneous provinces; Coffin & Eldholm, 1994). Many geohazards, mineral, and geothermal resources are associated with crustal magmatism. Understanding crustal magmatism on Earth—its underlying, unifying processes, consequences, tempo and evolution—is an important and practical problem central to various aspects of Earth evolution.
In this study, we sketch a metamodel for crustal magmatism. This represents a convenient framework for classifying the plethora of possible evolutionary paths pertinent to crustal magmatism in batholithic, island arc and continental environments. We then focus on one of the most spectacular ways continental crust can grow. Specifically, we quantitatively test proposed petrogenetic scenarios using constraints from phase equilibria for three Quaternary large-volume ignimbrite deposits in the western USA. In particular, we study (1) the Bishop Tuff of Long Valley in eastern California, USA (Hildreth & Wilson, 2007), (2) the Bandelier Tuff in New Mexico, USA (Rowe et al., 2007), and (3) the tuffs of Yellowstone, Montana, USA (Bindeman & Valley, 2001; Bindeman et al., 2008). Upon identifying disparate aspects of the existing scenarios, we adopt the parsimonious approach of finding the very simplest model consistent with the observed available phase relations, volcanological, and geochemical data for each system. We find that the data do not demand models that incorporate complicated temporal sequences of events. For each system, models dominated by fractional crystallization can explain most of the observations. We advocate that a simple scenario that explains most observations is advantageous compared with one that is more complicated and often impossible or very difficult to test quantitatively because of too many unresolved degrees of freedom. Although it is very likely that processes other than fractional crystallization contribute to giant ignimbrite petrogenesis, we argue that fractional crystallization is the underlying most quantitatively significant process. Hence, by quantifying fractional crystallization, we establish a ‘reference model’ for each system. Deviations between the reference model and observations then allow one to evaluate quantitatively more exotic petrogenetic processes. This approach is consistent with the metamodel briefly described below.
The genetic link between large-volume ignimbrites, Cordilleran batholithic rocks (the diorite–tonalite–granodiorite–granite association), and intermediate to silicic volcanic rocks (the basaltic andesite–andesite–dacite–rhyolite association) of subduction zones along active continental margins has long been explored by petrologists and geochemists (e.g. Smith, 1960, 1979; Shaw, 1965; Chappell & White, 1974, 2001; Gill, 1981; Lipman, 1984; Vernon, 1984; Hildreth & Moorbath, 1988; Barbarin, 1990; Atherton, 1993; Atherton & Petford, 1993; Pitcher, 1993; Brown, 1994; Chappell, 1996; Cobbing, 2000; Annen & Sparks, 2002; Bachmann et al., 2002; Jellinek & DePaolo, 2003; Annen et al., 2006). In detail, however, each petrotectonic environment—indeed each magmatic system—possesses a unique blend of textural, petrological, rheological, and geochemical attributes, the interpretation of which has sparked significant controversy during the past century. Attempts to develop a unified conception of crustal magmatism date to the very beginnings of igneous petrology (e.g. Iddings, 1909; Grout, 1918; Bowen, 1928; Kennedy & Anderson, 1938). These controversies have continued, in modified form, in more recent times (Chappell & White, 1974; Brown, 1979; Clarke, 1992; Atherton, 1993), as petrologists and geochemists have attempted to understand the relative importance via deconvolution of fractional crystallization, partial melting, magma mixing, liquid immiscibility, and other multistage processes (e.g. melt rejuvenation or ‘defrosting’) pertinent to petrogenesis. Regarding giant ignimbrite-forming systems, very active debate surrounds the relative roles of partial melting versus fractional crystallization versus melt rejuvenation and magma mixing in generating large volumes of rhyolitic melt within the upper crust (depth <10 km).
Most investigators now recognize the significant role of ascent and intrusion of mantle-derived mafic magma in the transport of heat into the crust. Advected heat provides the energy to power the observed large volumes (10–1000 km3) of low-crystal-content magmas that reside for long periods (∼104–107 years) within the upper to middle continental crust. An important question is whether support for long-term lifetimes of intermediate to silicic systems in continental and island arc settings requires heat delivered by continuing advective transport of basaltic magma generated within the mantle. If so, then crustal magmatic systems are best viewed as diathermal systems. The specific (per unit mass of magma) heat requirement () for replacing ambient upper crustal rock at temperature Tc with evolved magma at Tm is where Cp is the average specific isobaric heat capacity of crust (c) or magma (m), T is temperature, and Δhm is the magma crystallization enthalpy. For typical values ( = 1000 J/kg K, = 800 J/kg K, Δhf = 400 kJ/kg, Tm = 1300 K, Tc = 600 K, ρc = 2700 kg/m3, ρm = 2500 kg/m3; Spera, 2000), the energy per unit mass requirement is ∼1·2 MJ/kg. This estimate provides a basic constraint on any mechanism for the generation of large volumes of evolved magma: heat of ∼1 MJ for each kilogram of evolved magma produced needs to be supplied over characteristic timescales of ∼105–107 years to focused regions of relatively small volume (∼1000 km3). This requirement is consistent with observed rates of crustal magma generation (10−4–10−2 km3/a) in terrestrial magmatic systems (Smith, 1979; Spera & Crisp, 1981; Raia & Spera, 1997; Petford & Gallagher, 2001; White et al., 2006) and observed volumetric rates of magma generation within the mantle of ∼10−2–10−1 km3/a. Crustal radioactivity and frictional melting are too small by orders of magnitude to meet the energy requirement, as is the transport of heat by geothermal fluid migration (e.g. decompression and ascent of geothermal fluids). On the other hand, thermodynamic calculations show that the heat liberated by cooling and crystallization of basaltic melt from its liquidus to its solidus is 0·8–1·3 MJ/kg of magma. This is sufficient to power crustal igneous systems. The range of heat available encompasses different mantle melt compositions (including volatile contents), pressure (which affects the solidus to liquidus crystallization interval), and oxygen fugacity. Simply put, the melt isentrope of ∼0·4 K/km is significantly smaller than typical crustal geothermal gradients of 20–40 K/km. This implies significant vertical transport of heat during magma ascent. The point is that the energy for crustal magmatic systems ultimately depends on a mantle source.
Although certainly important, energetics are not the whole story. The material source of erupted magma evidently is important. There are two end-member possibilities, with petrological reality residing somewhere between. The first is anatexis of existing subsolidus country rock. In this scenario, mafic magma simply acts as a heat-transfer fluid. Intruded magma cools and crystallizes and delivers its sensible and latent heat to the crust, thereby triggering partial melting of pre-existing crust. In this extreme case, the generated crustal magma derives no material from the intruded mafic melt. At the other extreme is the scenario of fractional crystallization, where perfect closed-system fractionation of a mafic parental melt generates the evolved magma. The compositional lever is approximately 10:1. That is, evolved melt is generated by crystallization of ∼90% or more of parental mafic melt. On the spectrum between these extremes are scenarios that can be relatively simple or fiendishly complex. A simple example is fractionally crystallizing parental magma contaminated by anatectic liquids generated by partial melting of surrounding country rock. A more complex example is a multistage (polybaric) one in which parental mafic melt undergoes significant lower crustal crystallization involving generation of a mush region and associated evolved silicic melt that is later remobilized by fresh intrusion of mafic magma. There are innumerable possible scenarios in which the importance of melt recharge, assimilation, and crystallization (isobaric or polybaric) play roles of varying quantitative importance spatially and temporally. A problem with complex models is that they are difficult to falsify. For example, a polybaric model based on assimilation of partial melts derived from various crustal source compositions, variable amounts of recharge of fresh melt of varying composition, and crystallization that is neither fractional nor equilibrium and that occurs at a range of depths (polybaric) is virtually impossible to model in detail; there simply are too many variables and poorly known initial conditions to account for. Such scenarios are commonly based on a limited subset of data, such as the isotope signatures of a limited subset of components that stabilized at an unknown stage in the evolution of a particular magma body, rather than on phase equilibria considerations of the components that make up the bulk of a system. In complex systems there is always the danger of presenting a ‘just-so’ scenario based on too few data and too many degrees of freedom. The primary purpose of this study is to quantitatively investigate the consequences of the very simplest scenario for the origin of silicic magma in the upper crust: perfect fractional crystallization of a single parental melt at a single depth (under isobaric or isochoric conditions) along an oxygen buffer, where the components that make up the bulk of the system are accounted for. We recognize that just because a scenario is complex does not mean it is incorrect; similarly, because a model is simple and testable, with limited degrees of freedom, does not mean it is correct. However, simple models have the advantage of being testable and of providing a reference state. That is, in comparing observations with predictions, the shortcomings of the model can be distinguished. More complex scenarios can be developed starting from a simple model by adding on complications sequentially and testing the consequences of each additional process. For example, if a simple isobaric fractionation model does not reproduce the observations well, one can modify the reference model by allowing for polybaric crystallization. If isotopic data indicate some degree of crustal contamination, then the reference model can be modified to include assimilation. If there is strong evidence of magma recharge, then the reference model can be modified to allow input of fresh parental-type magma to the evolved liquid. Our basic assumption is that even when complicating processes are taking place, the underlying process of fractional crystallization remains important. Fractional crystallization is the background process-prism in a sense, with which other processes can be profitably viewed.
We apply this approach to three well-documented natural systems. The underlying hypothesis is that these systems possess some degree of commonality and we test this hypothesis via phase equilibria. The basic question is, if we assume the very simplest model, how large is the deviation between observed melt and phenocryst compositions and the model predictions?
A METAMODEL FOR CRUSTAL MAGMATISM
Most petrologists agree that three essential elements of crustal magmatism are: (1) generation of a primary melt by partial melting of the mantle (peridotite or peridotite–eclogite mixture); (2) processing, by crystal fractionation and other processes, of the primary melt during transport through the mantle, thereby generating parental melts; (3) further processing of the parental melts during magma ascent through, or storage within, the crust, ending in the generation of intermediate to silicic melts. The term ‘processing’ is used here sensu lato. It may include classical (Bowen-style) fractional crystallization (FC) under polybaric or isobaric conditions or along other thermodynamic paths (e.g. isentropic decompression, isochoric crystallization), development of an exsolved fluid phase that may be expelled from or retained within the magma body, assimilation (A) of melt generated by partial melting of country rock or of stoped blocks, recharge (R) of ‘fresh’ parental melt into a pre-existing magma body, and mixing (M) of recharge magma with pre-existing magma that has undergone some in situ evolution (Fig. 1). These and other processes can take place at virtually any depth within crust of widely varying composition, temperature, thermophysical properties, and local stress state to generate an almost infinite array of possible temporal sequences, leading to a wide variety of end states illustrated schematically in Fig. 1. Understanding the relative importance of the aforementioned mechanisms, the depths at which they occur, and their temporal order, forms the essence of the petrogenetic problem.
The ability to gather information on small scales by, for example, crystal isotope stratigraphy (e.g. Wolff & Ramos, 2003), crystal size distribution analysis (e.g. Cashman & Marsh, 1988; Cashman, 1993), high-resolution zircon geochronology (e.g. Reid, 2008), and study of glass (melt) and fluid inclusions in phenocrysts (e.g. Anderson et al., 2000; Webster et al., 2003; Cannatelli et al., 2007) has provided significant constraints for testing various petrogenetic hypotheses. Fortunately, there are now many hundreds of detailed case studies of plutonic, hypabyssal, and volcanic systems. Phase equilibria models (e.g. MELTS), experimental studies of specific systems (e.g. Rutherford & Hill, 1993; Rutherford et al., 1998; Nekvasil et al., 2004), and trace element models (DePaolo, 1981; Bohrson & Spera, 2007) can be applied routinely. Below we highlight a few of the many paths implicit in the metamodel of Fig. 1.
In some cases, the parental melt may ascend directly to the surface (left-hand side of Fig. 1) and erupt with relatively little modification. Quaternary xenolith-bearing, small-volume alkali basalt lava and scoria deposits from the Great Basin, USA of North America result from this branch of the metamodel. Such magmas ascend at rates of ∼0·01–1 m/s, traversing the crust relatively rapidly (1 day to 1 month), insufficient for marked changes in melt composition by assimilation or magma mixing. The compositional evolution of the melt that does occur in this case is driven by fractional crystallization, perhaps by the mechanism of flowage crystallization and reaction along conduit–mantle contacts (Irving, 1980; Spera, 1980; Kelemen, 1990). Along another possible branch of the metamodel, the intruded melt may solidify completely in the subsurface, liberating the heat needed to trigger voluminous production of in situ melts derived by partial melting of local crust. The shadowy area of transition from migmatite to granite encapsulated within the right-hand side of Fig. 1 recognizes the blurred boundary between magma and country rock. Along this branch, evolved melts form by partial melting of pre-existing crust with no additions from mantle-derived melts. In another alternative, when only a portion of the intruded melt solidifies, residual melt may ascend to shallower levels and undergo further evolution (bottom and middle portions, Fig. 1). Recently, the idea of remobilization (‘defrosting’) of a mostly crystallized mush by heat brought into the system via addition of ‘hot’ recharge magma has been suggested (e.g. Bachmann & Bergantz, 2003). This idea links recent material addition with the energy needed to remobilize melt from a previously existing porous network. In general, recharge can occur at any time after an initial intrusion of magma, and mixing can either go to completion, forming hybridized magmas, or remain incomplete (mixed magmas) until the moment of eruption and hence be observed in the surface deposit.
Some key elements governing magma evolution during storage periods include the amount of heat exchange between magma and host environment, as this governs the extent of fractional crystallization and assimilation, the depth (pressure) of magma processing, which exerts a significant effect on phase equilibria, and the volumetric rate of magma (and heat) recharge. The orientation and magnitude of the principal stress tensor is important in controlling the form magma bodies acquire and hence the probability of eruption. Also important is the influence of a supercritical fluid phase if one develops (‘second boiling’) or the incorporation of volatile-bearing wall-rock during magma evolution. The thermophysical state of the country rock surrounding the magma body, especially the spatial variation and anisotropy of permeability, modulates heat transfer across the boundary between magma and the surrounding hydrothermal system.
Given the large number of variable parameters defining a given system, it is not surprising that so many seemingly distinct scenarios have been proposed. Although there has been great interest and development of magma dynamics since the 1970s, the obvious should be explicitly stated: magma dynamics is ultimately built upon a foundation of phase equilibria and chemical thermodynamics.
The metamodel provides a framework for systematic quantitative evaluation of the most important processes affecting the evolution of particular magmatic systems. As a first step, in this study we focus on closed-system evolution. Our results show that this may not be an unreasonable approach. For the Bishop and Yellowstone Tuffs, we have found that a major contribution from open-system processes is not required to account for either crystal age discrepancies within a single deposit or phase relations. During fractional crystallization, solids can stabilize and persist over varying timescales. Although isotopic evidence (discussed below) points to some role for magma contamination by assimilation or magma mixing in each system, we find that data available to support the role of assimilation are generally linked with phases (sanidine, quartz, accessory phases) that, according to phase equilibria, stabilized near the end of closed-system crystallization. The timing of possible magma contamination via stoping or by addition of wall-rock anatectic melts versus the chronology of fractional crystallization remains an important, but unresolved issue. As noted below, in isochoric crystallization the mean magma pressure first falls below the local lithostatic value and then increases beyond the lithostatic value near the solidus. One might expect assimilation to be most important during time intervals where the magma pressure is less than the lithostat. These intervals are as yet sparsely represented in the isotopic database.
GEOLOGICAL SETTING AND PETROGENETIC SCENARIOS
In Electronic Appendix 1 (available for downloading at http://www.petrology.oxfordjournals.org/) we present a synopsis of the petrological settings and petrogenetic scenarios previously offered for three large-volume Pleistocene ignimbrite deposits in western North America. These scenarios are evaluated in a subsequent section, followed by our own analysis of the petrogenesis.
The petrogenetic hypotheses evaluated here are encapsulated schematically within the context of the metamodel in Fig. 1. These hypotheses share many fundamental similarities. Most significantly, crystal fractionation is of primary importance; it is the mechanism that generates erupted high-silica rhyolite compositions in each case. Remobilization of the fractionation products leads to liberation of the rhyolite melts. For the Bishop Tuff, Hildreth & Wilson (2007) considered crystal fractionation accompanied by melt extraction and transport to a shallower reservoir. For the Bandelier and Yellowstone systems, respectively, Rowe et al. (2007) and Bindeman et al. (2008) considered crystal fractionation followed by rejuvenation or remelting. The Bandelier rhyolites undergo minimal compositional modification during remobilization. The Yellowstone rhyolites are hydrothermally altered prior to remelting. In general, assimilation of crustal partial melts is not considered important in rhyolite production, although precursor magmas, such as the pre-caldera mafic lavas associated with all three systems, are hypothesized to have undergone assimilation, among other processes. It is interesting to note that although trace element and isotopic evidence testifies to the occurrence of contamination in each system, the implications and importance of these patterns are weighted very differently in each proposed scenario. In the following sections we explore the extent to which each scenario can explain the observed phenocryst and major element liquid (glass) compositions in each system.
APPLICATION OF MELTS TO NATURAL SYSTEMS
Phase equilibria constraints on major element changes in magmas during crystallization form a fundamental starting point in the discussion of magma compositional evolution. MELTS is a thermodynamic model of crystal–liquid equilibria that uses experimentally derived data on the compositions of coexisting solid(s) and liquids at specified temperatures, pressures, and oxygen fugacities to calibrate models for the compositional dependences of thermodynamic potentials for phases in equilibrium. As heat is extracted and the temperature drops in a system, phase identities, compositions, and proportions are calculated. The MELTS algorithm is based on classical equilibrium thermodynamics and has been extensively reviewed (Ghiorso & Sack, 1995; Asimow & Ghiorso, 1998). Below we provide a brief summary of points that are relevant to the present study. From a practical perspective, using MELTS to model the compositional evolution of cooling magma requires specification of the initial state of the system and the constraints under which evolution proceeds. Initial conditions define the starting temperature, temperature step, pressure, and parental magma composition, including an initial water concentration. In this study we have chosen the liquidus temperature as the starting temperature; MELTS computes a liquidus temperature based on the specified initial liquid composition and pressure. An end temperature is selected by comparison of MELTS results and observed data. A system ferric/ferrous ratio must be defined, either from FeO/Fe2O3 analyses or from total Fe, in which case selection of an oxygen fugacity distributes iron in the liquid and solids between FeO and Fe2O3 according to the specified temperature and pressure. The constraints specify a reaction path in which the system is closed or open to mass transfer within the magma body (fractional or equilibrium crystallization) or from wall-rock to the magma body (assimilation). Pressure may be held constant or may vary along some P–T path. An oxygen fugacity constraint path may be defined. These and a number of additional issues must be considered when using MELTS to explore the consequences of magmatic evolution in natural systems.
Thermodynamic properties database
Compared with data from our investigation of some seven natural explosive volcanic systems, including those studied here, MELTS effectively predicts oxide concentrations for liquids during fractional crystallization, from the liquidus down to low melt fraction (∼0·04).
However, the predicted CaO and K2O concentrations are systematically displaced from the observed data by up to ∼3 wt %. The discrepancy for CaO can probably be attributed to under-stabilization of clinopyroxene (M. Ghiorso, personal communication). The problem with K2O may be due to a lack of experimental data to calibrate the thermodynamic model that describes alkali feldspar–liquid equilibria (M. Ghiorso, personal communication). The calibration of the activity of the K-component in the liquid may give a value that is too high, resulting in overprediction of the stability of alkali feldspar that is manifest in a reduction of liquid K2O concentrations.
Lack of accounting for CO2
MELTS currently does not have a thermodynamic model for CO2 solubility. However, the solubility of CO2 in magmas is minimal (Holloway, 1976).
H2O solubility model
The water solubility model incorporated in MELTS is adequate; it captures the experimental observation that water solubility increases as silica increases from basalt to rhyolite. There are issues that Gordon Moore and Mark Ghiorso have identified in calc-alkaline systems, but we are not dealing with calc-alkaline systems here. MELTS can deal with the effect of water dilution on the stability of nominally anhydrous minerals.
Amphibole and biotite models
MELTS does not perform well when amphibole and biotite are primary phases, as a result of problems with the related solid solution models. There are simply not enough good phase equilibria data to constrain liquid compositions when amphibole and biotite are modally important phases, as in calc-alkaline magmatic systems. However, amphibole and biotite are present in minor amounts (e.g. 1 vol. %, maximum) in the systems we are studying.
THERMODYNAMIC AND MASS-BALANCE CONSTRAINTS ON MAGMA EVOLUTION: METHOD TO ESTABLISH THE BEST CRYSTALLIZATION SOLUTION
Figure 2 shows a classification of closed-system thermodynamic and system mass-transfer paths associated with magma solidification, modified from Fowler & Spera (2008). These solidification paths represent end-member scenarios for crystallization of crustal magma bodies and their physical evolution when crystallization without significant recharge or assimilation is the dominant mechanism of evolution. Phases (precipitating crystals and/or exsolving fluid) can either remain in chemical potential equilibrium with residual melt (equilibrium crystallization) or be isolated from it, precluding further reaction with melt and/or other crystals (fractional crystallization). In terms of thermodynamic limits, solidification can take place either at constant pressure (isobaric) or constant volume (isochoric). For isobaric solidification, the independent variables include temperature, system bulk composition, oxygen chemical potential, and magma internal pressure, Pmagma. Isobaric crystallization corresponds to a magma chamber within weak country rock that expands or contracts perfectly to accommodate changes in the system volume accompanying phase transitions and maintaining fixed pressure. During isochoric solidification (Ghiorso & Carmichael, 1987), volume rather than magma pressure is the independent variable. Because phase change is generally accompanied by a change in system density, the isochoric constraint leads to variations in magma pressure along the crystallization path because the country rock provides a perfectly rigid container. This case is approximated when wall-rock deformation is too slow to accommodate volume change as a result of crystal and fluid bubble formation by ductile flow and the criterion for brittle failure in the country rocks is not exceeded.
Over 150 MELTS calculations have been performed and compared with available data for natural systems to develop bounds on the parameters that governed production of the Bandelier, Bishop, and Yellowstone Tuffs. Information about MELTS can be found in the Electronic Appendix. In the model scenario, mafic parental melt containing a specified amount of dissolved water crystallizes at fixed depth within weak (isobaric crystallization) or rigid (isochoric crystallization) crustal host-rock. Parental melts represent mafic precursor melts taken from published data that are specific to each of the systems studied. Each calculation begins above the liquidus, based on the parental melt composition and pressure (i.e. fixed pressure for isobaric calculations; initial pressure for isochoric calculations). Magma is assumed to be in lithostatic equilibrium with host-rock at the initiation of crystallization. We explored a variety of potential parental melt compositions over a grid of lithostatic pressures (0·05, 0·1, 0·2, 0·3, and 0·5 GPa), oxygen fugacity (QFM – 2 to QFM + 3, where QFM is the quartz–fayalite–magnetite buffer), and parental melt dissolved water concentration (dry up to H2O saturation at the liquidus). We imposed oxygen buffer constraints on the liquid oxidation state. Only total iron has been used as an input; the FeO–Fe2O3 redox state of evolving liquids reflects ferric/ferrous ratios along the buffer. Using the estimated initial parameters and calculation constraints, we computed phase equilibria for six of the possible eight solidification paths shown in Fig. 2: (1) isobaric solidification (crystals and H2O fluid bubbles retained in magma); (2) isochoric solidification (crystals and bubbles retained); (3) isobaric fractionation (crystals removed but bubbles retained); (4) isochoric fractionation (crystals removed but bubbles retained); (5) isobaric fractionation (crystals and bubbles removed); (6) isochoric fractionation (crystals and bubbles removed). In all cases, calculations end at temperatures approaching the solidus, generally when the fraction of melt (fm) remaining is ∼0·05.
From assessment of the results, we then chose a ‘best’ value for the initial pressure, oxygen buffer, and initial water content via comparison of predictions with observed pumice, glass, and bulk-rock major element compositions, and the identity and composition of all phenocryst phases. We found, in general, that varying the assumed parental melt composition within the bounds of possible parental melts based on regional forerunners leads to a restricted set of possible outcomes (i.e. liquid line of descent and phenocryst compositions). In this study, all calculations are based on a single representative parental melt composition for each system.
Below, we evaluate existing petrogenetic scenarios for the Otowi (∼400 km3) and Tshirege (∼200 km3) members of the Bandelier Tuff, the ∼600 km3 Bishop Tuff, and the 2500, 300, and 1000 km3 Yellowstone high-silica rhyolite tuffs. Site-specific data and interpretations published over the last ∼25 years provide starting conditions and geological constraints. We define reference states for the phase equilibria, major element, and physical property evolution of each system that can be tested through further collection of geochemical and volcanological data. For each tuff, we describe the results of phase equilibria calculations to illustrate the initial pressure (isochoric crystallization) or fixed pressure (isobaric crystallization), initial H2O concentration, and redox conditions that governed the liquid lines of descent and phase relations under isobaric and isochoric conditions. We find that fractional crystallization (FC) involving removal of exsolved water bubbles, starting from a relatively ‘wet’ parental basalt (3–3·5 wt % dissolved water), near water saturation at the liquidus temperature, and low starting pressure provides a good model for derivation of the explosive rhyolitic eruptive rocks. We have shown previously (Fowler & Spera, 2008) that a similar scenario works well to explain the generation of the 40 ka Campanian Ignimbrite in central Italy. For one system (the Bandelier Tuff), we present a sensitivity analysis that demonstrates how changes in intensive variables influence phase equilibria. All plotted major element data are presented on an anhydrous basis.
Evaluation of existing scenario
Here we evaluate the scenario of Rowe et al. (2007) for the petrogenesis of the Bandelier Tuff. Those workers identified melting of hybridized crust in the form of buried, warm intrusions, generated during pre-6 Ma Jemez Mountains Volcanic Field (JMVF) activity as an important mechanism in Bandelier Tuff petrogenesis. As evidence they cited Nd and Pb isotope similarity between the Bandelier Tuff and pre-caldera mafic Paliza Canyon volcanic rocks interpreted to have crystallized from crustally contaminated, mantle-derived melts. Rowe et al. (2007) rejected crystal fractionation based on a scarcity of evidence for a large mass of cumulate mafic residue beneath the JMVF and an absence of intermediate magmas that are isotopically similar to the Bandelier Tuff and Paliza Canyon volcanic rocks. However, investigating the composition of the crust beneath the JMVF is not straightforward; thermal effects related to continuing magmatic activity and probable extensive modification of the crust over the past ∼20 Myr may complicate recognition of cumulate deposits. Additionally, the results of fractional crystallization calculations presented below suggest that intermediate-composition lavas related to the Bandelier Tuff may be absent at the surface because they were un-eruptible.
We have compared Bandelier Tuff trace element concentrations with predictions from trace element modelling. The trace element models incorporate phase proportions and compositions determined from phase equilibria. Rowe et al. (2007) discussed the special challenges associated with tracking trace elements in high-silica rhyolites, where the presence of small amounts of accessory phases such as allanite, chevkinite, monazite, and zircon can complicate understanding of large ion lithophile element, high field strength element, and rare earth element partitioning. In addition, published crystal–liquid partition coefficients vary widely for many trace elements, sometimes over several orders of magnitude. Thus we consider only Sr and Ba, which extend to distinctly low concentrations in the Bandelier Tuff (∼2 ppm; Kuentz, 1986; Balsley, 1988) and have relatively restricted published partition coefficient ranges, leading to distinct liquid concentrations for some competing end-member petrogenetic scenarios.
First we performed MELTS isobaric batch melting calculations on Paliza Canyon samples (Rowe et al., 2007), with oxygen fugacity constrained at the QFM buffer and initial water set just below the concentration at which water exists as a free phase near the solidus. A pressure of 0·1 GPa was chosen for consistency with the shallow depth of the Paliza Canyon Formation (Rowe et al., 2007); small positive or negative pressure variations (i.e. + ∼0·1 GPa, – ∼0·05 GPa) do not lead to significant changes in the melting solid phase assemblage. After fixing the Sr and Ba concentrations in the bulk source and calculating bulk mineral–melt partition coefficients based on the phase equilibria results and mineral–melt partition coefficient values from the Geochemical Earth Reference Model (GERM) database (http://earthref.org/GERM/), we forward-modelled the consequences of partial melting by numerical solution to the batch partial melting equation. Source trace element abundances used in the calculations are from Rowe et al. (2007). To explore a wide range of potential source compositional variability, we chose as starting compositions samples with higher-silica, higher-Sr concentration (sample JM9307, 778 ppm Sr; Rowe et al., 2007) and lower-silica, lower-Sr (sample JM93124, 264 ppm Sr; Rowe et al., 2007) concentration. These samples have 1185 and 1345 ppm Ba, respectively. Accessory phases such as allanite, chevkinite, monazite, and zircon are not included within the MELTS database, but we considered their influence on Sr and Ba evolution by performing calculations in which we manually included 1 and 3% by mass in the computed solid phase assemblages of a hypothetical phase with a mineral–melt Kd for Sr and Ba of 100.
The computed solid phase assemblage for the lower-Sr, lower-SiO2 Paliza Canyon source rock (Rowe et al., 2007) consists of apatite, amphibole, plagioclase, alkali feldspar, quartz, rhombohedral oxide, and spinel. The higher-Sr, higher-SiO2 Paliza Canyon source rock composition differs in the presence of diopsidic clinopyroxene and orthopyroxene. Biotite and fayalite are absent.
The minimum observed Bandelier Tuff bulk Sr concentration is ∼2 ppm and the majority (111/134) of published Bandelier Tuff whole-rock Sr concentration values (Kuentz, 1986; Balsley, 1988; Spell et al., 1990; Wolff et al., 2005) are lower than 20 ppm. Under conditions chosen to maximize the coincidence between the calculations and the observed Bandelier Tuff data (i.e. using the Paliza Canyon intrusion parent with the lower initial Sr concentration along with the maximum mineral–melt Kd values included in the GERM database of ∼20 for plagioclase and alkali feldspar), the lowest computed melt Sr concentration is ∼80 ppm (no hypothetical accessory phase included). Addition of 3% by mass of a hypothetical accessory phase with a mineral–melt Sr partition coefficient of 100 leads to a minimum melt Sr concentration of ∼40 ppm. Minimum Ba concentrations for the two cases approach 100 ppm compared with minimum observed Bandelier Tuff Ba contents of ∼2 ppm. For lower Sr and Ba mineral–melt Kd values (4), minimum calculated Ba values are ∼400 ppm at the lowest melt fraction. For the higher initial Sr concentration source composition, the minimum Sr concentration is ∼160 ppm and the minimum Ba concentration is ∼1500 ppm. These results suggest that partial melting of Paliza Canyon intrusions may not have been the dominant petrogenetic mechanism for the Bandelier Tuff.
Kuentz (1986) proposed that fractional crystallization near a multicomponent invariant point generated the Otowi Member mineral assemblage and element variations, and we investigated this possibility through calculation of diverse crystallization paths presented below. Kuentz (1986) speculated that fractionation generated vertical gradients in composition, viscosity, and density, and led to accumulation of a volatile-rich cap at the top of the magma reservoir that subsequently generated the Guaje pumice deposit. We do not attempt trace element calculation of fractional crystallization as a result of the relatively more complex role of accessory minerals that potentially stabilize at a range of temperatures and melt fractions. However, we note that even in the absence of accessory minerals, fractional crystallization starting at JMVF mafic compositions, using reasonable partition coefficient values from GERM, can generate Sr and Ba concentrations similar to the lowest observed concentrations in the Bandelier Tuff. Because of the physical and chemical similarity between the Otowi and Tshirege members (Smith & Bailey, 1966), we use a single MELTS model to describe both Bandelier Tuff members.
Bandelier Tuff isobaric and isochoric MELTS crystallization calculations
In our examination of crystallization scenarios, we first performed numerical experiments for equilibrium crystallization (crystals, fluids, and melt remain in chemical equilibrium throughout the crystallization history). Low observed crystal contents and predicted trends on MgO–oxide variation diagrams that are in poor agreement with Bandelier Tuff data suggest that equilibrium crystallization was not an important petrogenetic process for the Bandelier Tuff, so the related discussion is brief. We follow with presentation of crystal fractionation calculation results, beginning with selected results of sensitivity tests to illustrate potential compositional variability introduced by independently varying initial conditions or constraints (e.g. pressure, redox condition, initial H2O content of parental melt). Somewhat surprisingly, upon conducting an extensive set of numerical experiments, we found that varying the thermodynamic constraint (isobaric or isochoric) and the mass-balance constraint (fractionation of solids or fractionation of solids + fluid) has only a small influence on liquid lines of descent. The sensitivity analysis is based on the case of isobaric crystallization with fractionation of both exsolved fluid and precipitated crystals; the variables are the initial water concentration, the chosen oxygen buffer, and the initial pressure.
We have explored parental melt compositions corresponding to a number of JMVF basalts; however, the results reported here are based on sample EA3 from Wolff et al. (2005), a tholeiite erupted within the El Alto mafic volcanic field in the northern part of the JMVF. El Alto lavas are petrographically and chemically similar to mafic lavas in the nearby Cerros del Rio volcanic field dated at 2·48 and 2·33 Ma (Woldegabriel et al., 1996), but may be slightly older; one El Alto flow has been dated at 3·2 Ma (Baldridge et al., 1980). The use of a slightly different parental anhydrous bulk composition from the JMVF does not significantly alter the conclusions based on the El Alto parental magma composition.
Bandelier Tuff isobaric and isochoric equilibrium crystallization MELTS calculations
Bandelier Tuff isobaric and isochoric equilibrium crystallization calculations are based on the best-fit model constraints (parental melt composition and its initial dissolved water concentration, oxygen buffer, and lithostatic pressure; described below). Despite isochoric-case variations in pressure (Fig. 1, Electronic Appendix), the resultant compositional paths are broadly similar (Figs 2–4, Electronic Appendix). Pressure first decreases, then increases, as temperature falls (Fig. 1, Electronic Appendix). Temperature is shown as a nondimensional variable () based on the liquidus and solidus temperatures specific to each composition (). The most important observation is that predicted trends agree poorly with Bandelier Tuff data. A marked discrepancy is the deficiency in calculated liquid SiO2 (Fig. 2a, Electronic Appendix); SiO2 attains a maximum of ∼70 wt %, ∼6–9 wt % lower than the bulk of Bandelier Tuff data. Also, calculated K2O and Al2O3 liquid lines of descent do not intersect the Bandelier Tuff data field (Fig. 2b and d, Electronic Appendix). Because MELTS generates equilibrium mineral compositional data for each step as temperature falls, a further tool for evaluating the MELTS calculation results is comparison of the observed and predicted phenocryst compositions. Significantly, alkali feldspar, anorthoclase, and quartz, all recorded phases in the Bandelier Tuff, are absent in the equilibrium crystallization models. The predicted compositional range of feldspar (An76–40) overlaps only with An-rich (An54) and intermediate Bandelier Tuff plagioclase (Fig. 3, Electronic Appendix). Clinopyroxene compositions agree fairly well with those of the more magnesian phenocrysts (<Fs33; Fig. 2, Electronic Appendix); however, isobaric-case clinopyroxene compositions are restricted to <Fs20 (Fig. 4, Electronic Appendix). Calculated orthopyroxene, although not well constrained by observations because of a scarcity of published data, does overlap partly with the observed range (Fig. 4, Electronic Appendix). Fayalitic olivine is not predicted; olivine is restricted to <Fa68. Based on these results, we reject equilibrium crystallization as an important mechanism for the evolution of the Bandelier Tuff. In the next step, we systematically evaluate model parameter sensitivity for fractional crystallization.
Effects of varying the initial water concentration: fractional crystallization
We examined the influence of parental melt dissolved water concentration by performing isobaric (∼0·13 GPa) fractional crystallization (crystals and exsolved fluid removed) calculations with oxygen fugacity constrained to follow the QFM buffer and 0·5, 1, 2, 3, and 3·5 (saturation) wt % H2O. For 0·5, 1, and 2 wt % initial H2O, MgO vs K2O and Al2O3 trends are respectively ∼1·5 and >6 wt % lower than the Bandelier Tuff data field at low MgO (Fig. 3b and d). In contrast, MgO vs Na2O trends (Fig. 3c) exceed observed concentrations (0·5 and 1 wt % H2O) or do not coincide with the bulk of the data (2 wt % H2O). Maximum SiO2 concentrations are higher for the 0·5, 1, and 2 wt % initial H2O cases (∼80 wt %) than for the 3 and ∼3·5 wt % cases (∼78 wt %), more closely approximating the range of SiO2 observed in the Bandelier Tuff (Fig. 3a).
These variations are due mainly to the behaviour of feldspar (Fig. 4). Increasing initial water concentration results in progressively lower feldspar saturation temperatures, higher total feldspar abundances (0·52 at 0·5 wt % initial H2O, 0·58 at ∼3·5 wt % H2O), and plagioclase destabilization in favour of relatively Al2O3-depleted, SiO2-rich alkali feldspar and anorthoclase (Fig. 4). The appearance of quartz at initial H2O ≥2 wt %, consistent with observations of quartz phenocrysts in the Bandelier Tuff, is an additional factor leading to lower liquid SiO2 contents at higher initial H2O concentrations (Fig. 4). Although the calculated ratio of plagioclase to alkali feldspar plus anorthoclase falls as initial H2O rises, liquid K2O concentrations are lower at lower initial H2O because alkali feldspar becomes more potassic (Fig. 5a; 0·5 wt % H2O: Or46–59; 1 wt % H2O: Or55–58; 2 wt % H2O: Or49–59; 3·0 wt % H2O: Or40–54; 3·5 wt % H2O: Or44–55). The higher initial H2O cases provide the best agreement between observed and predicted feldspar, but none of the models fully reproduces the range of observed anorthoclase compositions. Anorthoclase phenocrysts that do not coincide with model predictions (Fig. 5a; ∼Or32–40) are mostly interpreted as cognate (Caress, 1994).
Pyroxene predictions for the 0·5 wt % H2O case best coincide with observed Fe-rich clinopyroxene crystals (Fig. 5b), but only predictions based on the 3·0 and ∼3·5 wt % initial H2O cases overlap Fe-rich orthopyroxenes (Fig. 5b). However, we show below that Fe-rich pyroxene is stable at slightly lower pressure (0·10 GPa), at initial water concentrations ≥3 wt %. All input water concentration values yield fayalitic olivine (>Fa80) at low melt fraction, plus spinel (Fig. 5c) and rhombohedral oxide (Fig. 5d). The compositions of these phases are minimally dependent on initial H2O concentration.
Bandelier Tuff oxide data correlate poorly with MELTS model trends for parental melt dissolved water concentrations of 0·5, 1, and 2 wt % H2O, suggesting that the concentration of dissolved water was relatively high, from ∼3·2 to 3·5 wt %, approaching water saturation at the liquidus at 0·1–0·15 GPa.
Effects of varying pressure: fractional crystallization
To explore the influence of pressure, we compared phase equilibria calculation results for isobaric fractional (crystals plus fluid) crystallization at systematically varying pressure (0·05, 0·10, ∼0·13, 0·20, 0·30, and 0·50 GPa), with oxygen fugacity and parental melt dissolved water concentration fixed at QFM and 3·5 wt % H2O, respectively. This pressure range corresponds to equilibration depths from 1 to 15 km, roughly. Figures 6–8 show the calculated liquid major element trends, phase proportions, and solid phase compositions.
At 0·2, 0·3, and 0·5 GPa, MELTS liquid and solid phase oxide trends depart significantly from observations (Fig. 6). The 0·3 and 0·5 GPa cases have leucite, garnet, and muscovite; trace quantities (<1% by mass) of leucite and garnet are also present at 0·2 GPa (Fig. 7). These phases are not observed in the Bandelier Tuff. Below 0·2 GPa, the agreement between observed and model solid phases is broadly consistent; however, at 0·05 GPa, quartz, a ubiquitous phase in the Bandelier Tuff, is absent (Fig. 7).
As lithostatic pressure rises, the stabilization temperature and mass fraction of plagioclase decreases progressively as more silica-rich phases stabilize (Fig. 7), so that increasing pressure leads to lower liquid SiO2 and higher Al2O3 concentrations compared with the Bandelier Tuff data field (Fig. 6). At 0·05 GPa, the total proportion of plagioclase is twice as high as in the 0·5 GPa case (Fig. 7). At 0·2, 0·3, and 0·5 GPa, liquid SiO2 concentrations achieve maxima of at least 5 wt % below the Bandelier Tuff data field, then plunge with falling temperature as quartz appears (Fig. 6). An increasing mass proportion of clinopyroxene at the expense of olivine is an additional factor leading to liquid SiO2 depletion with increasing pressure (Figs 6 and 7). Na2O and K2O trends are inconsistent with Bandelier Tuff data at higher pressure; at 0·2 and 0·3 GPa, a relative abundance of leucite, anorthoclase, and potassic alkali feldspar (Fig. 7) leads to Na2O and K2O concentrations that are significantly lower than the Bandelier Tuff data field (Fig. 6). Na2O concentrations are higher than the bulk of Bandelier Tuff data at 0·5 GPa as a result of relatively low plagioclase and alkali feldspar Na2O contents (Fig. 8a) and a jump in the mass fraction of garnet to ∼0·1 at the expense of an equivalent mass of plagioclase (Fig. 7). Olivine >Fo75 is a stable phase at 0·2, 0·3, and 0·5 GPa—these higher pressure models do not reproduce the range of fayalitic olivine compositions documented in the Bandelier Tuff (Fig. 8b). The agreement between model and observed clinopyroxene is fairly good at 0·5 GPa, but orthopyroxene is absent (Fig. 8b). Orthopyroxene is present at 0·2 and 0·3 GPa, but model clinopyroxene is restricted to diopside (Fig. 8b). At all pressures, model and observed compositional data are in good agreement for spinel and rhombohedral oxide (Fig. 8c and d).
Liquid SiO2 concentrations are highest at 0·05 and 1·0 GPa, coinciding with the full range of Bandelier Tuff SiO2 (Fig. 6). However, the relatively high abundance of plagioclase at the lowest pressures (0·05 and 1·0 GPa; Fig. 7) leads to liquid Al2O3 concentrations at low melt fraction that are up to ∼7 wt % lower than the observed data (Fig. 6). At 0·20 GPa, where plagioclase is relatively less stable, liquid Al2O3 concentrations exceed those of the Bandelier Tuff data field by >2 wt %. Melt Na2O concentrations coincide only with the most sodic Bandelier Tuff data (Fig. 6) as a result of predicted alkali feldspar compositions that are restricted to >Or48 and an abundance of sodium-poor plagioclase (Fig. 8a). The 0·10 and 0·20 GPa cases provide the closest agreement with measured feldspar compositions (Fig. 8a). At 1· 0 GPa, the predicted clinopyroxene composition extends towards hedenbergitic Bandelier Tuff pyroxene data (Fig. 8b); orthopyroxene and fayalitic olivine are present at the lowest pressures (Fig. 8b).
The misfit of oxide data, coupled with the mismatch between observed and predicted phases, eliminates the possibility that the Bandelier Tuff magma evolved at pressures of 0·2 GPa or higher or at pressures lower than 0·1 GPa. We performed an additional calculation at 0·13 GPa (see below). Based on the results we conclude that pressure between ∼0·1 and 0·13 GPa yields reasonably good agreement between the model and Bandelier Tuff compositions. This pressure range corresponds to a depth of 4–5 km in the middle part of the upper crust.
Effects of varying oxygen fugacity: fractional crystallization
We performed isobaric fractional crystallization calculations based on removal of crystals and fluid, with pressure and initial water concentration set at ∼0·13 GPa and ∼3·5 wt %, respectively, to examine the effects of varying oxygen fugacity (QFM – 2, QFM – 1, QFM, QFM + 1, QFM + 2, and QFM + 3). We found that Bandelier Tuff liquid major element trends (Fig. 9) are not well constrained by oxygen fugacity, but oxygen fugacity does have a strong influence on the relative composition and stability of Fe2+- and Fe3+-bearing solid phases (Fig. 10).
At QFM, QFM – 1, and QFM – 2, fayalitic olivine is stable at low melt fraction, consistent with observations on the Bandelier Tuff. Olivine is stable at QFM + 1 and QFM + 2, but only at high melt fraction and >Fo80 (Fig. 11a). Olivine is absent at QFM + 3.
Figure 11a summarizes model pyroxene compositions. At QFM + 1 to QFM + 3, orthopyroxene En contents exceed observed values by up to 25–35%. Measured data, although sparse, are reasonably approximated at QFM. Orthopyroxene is absent at QFM – 1 and QFM – 2. The QFM to QFM – 2 buffers provide the closest agreement between model and Bandelier Tuff clinopyroxene compositions. As mentioned above (Fig. 8a) but not shown here, clinopyroxene is more Fe rich at slightly lower pressure (0·1 GPa), overlapping with the more hedenbergitic Bandelier Tuff clinopyroxene phenocrysts.
Plotted in Fig. 11b and 11c are spinel and rhombohedral oxide model results. Variations between the models are not significant, but the observed spinel compositions are in closest agreement with the QFM-based model results. Calculated feldspar compositions are relatively insensitive to model oxygen fugacity (Fig. 11d).
We conclude that the fugacity of oxygen during crystallization was below QFM + 1. QFM provides the closest agreement between the calculated and actual Bandelier Tuff compositions. For simplicity, all results presented from this point forward are based on QFM.
Bandelier Tuff isobaric and isochoric fractional crystallization
Based on the results above, we conclude that the Bandelier Tuff evolved via fractional crystallization of a water-rich (∼3 wt % up to a liquidus saturation concentration of 3·5 wt % H2O) parental melt along the QFM buffer at c. 4–5 km depth. The QFM buffer constraint is consistent with the findings of Warshaw & Smith (1988), who concluded that the Bandelier Tuff crystallized at or below the QFM buffer from analysis of compositional variations of ferromagnesian silicates within the Tshirege Member. We use these results as a foundation for further calculations where we assume (a) isobaric solidification, crystals fractionated and H2O bubbles retained; (b) isochoric solidification, crystals removed but H2O fluid bubbles retained; (c) isobaric solidification, crystals and H2O bubbles fractionated; (d) isochoric solidification, crystals and H2O bubbles fractionated. These models are all based on values of ∼3·5 wt % and 0·13 GPa for parental melt dissolved water concentration and fixed (isobaric) or initial (isochoric) pressure.
The four calculations are all broadly consistent with the observed phenocrysts, and the temperature versus phenocryst proportion diagrams (Fig. 12) share many significant features. Olivine, the universally predicted liquidus phase, is not observed; its absence can be explained by loss via crystal settling in relatively low-viscosity magma (hot, nearly aphyric and water-rich). Following olivine precipitation in succession are clinopyroxene, spinel, plagioclase, and apatite. As temperature continues to decrease in cases (a), (b), and (c), orthopyroxene, anorthoclase, alkali feldspar, rhombohedral oxide, and quartz appear in succession with falling temperature. In case (d), rhombohedral oxide is not stable. Also, orthopyroxene appears after alkali feldspar and anorthoclase, and is followed by leucite, quartz, and amphibole. Leucite phenocrysts are not observed in the Bandelier Tuff. Biotite and hornblende are observed, but biotite is calculated only under isochoric conditions, whereas hornblende is stable only in case (d). However, minor positive or negative variations in input lithostatic pressure yield both phases (Fig. 7). The minor phases zircon, allanite, and chevkinite, also observed, are not included within the MELTS database, and are therefore not predicted.
Major element trends are similar in all four calculations (Fig. 13); however, in case (d), leucite crystallization leads to lower K2O at low melt fraction. In addition, the concentration of SiO2 is ∼2·6 wt % lower than in the other three cases in case (d), where the pressure at low melt fraction is relatively high, as a result of a higher total proportion of quartz (∼6·5%) relative to cases (a), (b), and (c) (∼3·5% quartz).
Feldspar, pyroxene, and olivine compositions are shown in Fig. 14. The agreement between predicted and observed phases is good, considering the computed trend assumes perfect fractional crystallization of solid phases for an inferred parental melt composition and water content at either a single pressure (isobaric) or an initial pressure (isochoric). Predicted phenocryst compositions generally overlap, but also reflect minor differences in the modelled reaction paths for each system. At low melt fraction in case (d), plagioclase is up to ∼5–10% less potassic and alkali feldspar is up to ∼30% more potassic than the observed data and cases (a), (b), and (c) (Fig. 14a). Also in case (d), pyroxene compositions are generally more iron-rich, overlapping with the most Fe-rich Bandelier Tuff pyroxenes (Fig. 14b). Fayalitic olivine is stable in all cases except (d) (Fig. 14b).
Figure 15a–e shows the variation of temperature versus physical properties for the four fractional crystallization scenarios for each of the large-volume eruptions investigated. Although we focus on the Bandelier Tuff at present, it should be noted that the other systems share many features. As noted by Fowler & Spera (2008), marked physical changes occur along the liquid line of descent in each case; changes are particularly striking near the lowest temperatures [710°C, = 0·02 in cases (a) and (b); 700°C, = 0·03 in case (c), 670°C, = 0·04 in case (d)]. These points correspond to multicomponent pseudoinvariant points where the temperature is constant (within a few degrees) during heat extraction. We designate these points as temperature Tinv. The isochoric constraint [cases (c) and (d)] leads to negative and positive variations in pressure during cooling (Fig. 15a), with kinks corresponding to the appearance of new phases. In spite of general similarities in the liquid lines of descent for the isobaric and isochoric cases, isochoric versus isobaric crystallization leads to significant differences in the evolution of some physical properties, primarily driven by differences in pressure shown in Fig. 15a. In particular, pressure varies during isochoric crystallization and that variation does affect the phase relations.
In the isochoric cases, the first interval of crystallization begins with the appearance of olivine and ends with water saturation. The higher density of olivine compared with melt results in an initial internal pressure decrease in the melt as temperature decreases (Fig. 15a). The pressure decrease is relatively steep, even at minor degrees of crystallization, because the density of olivine is significantly higher than that of the melt. The lower system pressure causes water to saturate at higher temperature (1130°C; = 0·96) than in the isobaric cases (1010°C; = 0·67). During subsequent crystallization, olivine is replaced by clinopyroxene, then spinel, and the system pressure decrease becomes more gradual, reaching a minimum P/Po of ∼0·8. Although spinel and clinopyroxene are denser than the coexisting melt, continuing exsolution of water, a low-density phase, ‘buffers’ or moderates the pressure variation. The onset of plagioclase stability at 1010°C ( ≈ 0·7) coincides with an abrupt pressure rise. Plagioclase is the dominantly stable phase during the succeeding crystallization interval, saturating at ∼1000°C after ∼20% crystallization along all four reaction paths. This interval represents a significant extent of crystallization—plagioclase makes up ∼37–40% of the calculated final phase assemblages—and internal magma pressure increases for the Bandelier, Bishop, and Yellowstone systems as a result of the high proportion and low density of Ab-rich plagioclase relative to coexisting melt. All three systems attain their initial pressures after rather minor plagioclase precipitation, at 34–38% of magmatic crystallization. The final crystallization interval begins with steepening pressure increase driven by alkali feldspar stabilization and increased H2O exsolution [800°C, = 0·23 in case (c); 780°C, = 0·24 case (d)]. At the appearance of alkali feldspar, phase types and phase proportions for the two Bandelier Tuff isochoric cases are broadly similar to each other, although pressure is higher in case (d) [P/Po = 1·23 in case (c); P/Po = 1·79 in case (d)]. Continuing crystallization leads to stabilization of leucite in case (d), as the system achieves P/Po = 2·13. Leucite is not an observed phase in the Bandelier Tuff. Its presence affects calculated major element trends insignificantly because of its small proportion and late appearance. Following the appearance of leucite, varying stability of leucite, quartz, exsolved H2O, and alkali feldspar with decreasing temperature drives abrupt positive and negative pressure variations of up to P/Po ∼ 0·03 towards the end of the crystallization interval. The maximum calculated internal magmatic overpressure towards the conclusion of isochoric crystallization is P/Po = 0·31. This high pressure near the lowest melt fractions is the source of differences in the low-melt fraction phase assemblages between case (d) and cases (a), (b), and (c). The retention of H2O bubbles and crystallization of phases less dense than melt leads to a flattening off, then decrease in pressure in case (c). In this case, the maximum pressure during the final phase of crystallization is P/Po = 1·25, with a final pressure of P/Po = 1·16.
In the isobaric cases [(a) and (b)], melt saturates with respect to H2O (Fig. 15b) at 1010°C ( = 0·69), at ∼4·3 wt % dissolved H2O (Fig. 15c). The concentration of dissolved H2O continues to increase to ∼4·7 wt % H2O as heat is extracted and crystals precipitate. At ∼900°C ( ≈ 0·44), dissolved H2O levels off as the mass of anhydrous solid phases produced at each temperature step decreases, until at Tinv there is a small rise in H2O content, from ∼4·13 wt % to ∼4·25 wt %. In the isochoric cases, melt saturates with respect to H2O at 1130°C ( = 0·96; 3·52 wt % dissolved H2O). The dissolved water content increases to a maximum of 5·14 wt % in case (c) and 11·3 wt % in case (d) (Fig. 15c).
In Fig. 15b, variation in the volume fraction (θ) of supercritical fluid along the crystallization paths is shown. It is important to recognize that at temperatures at and just below Tinv, ∼90 wt % of the original melt has crystallized. The viscosity of a hypothetical homogeneous multiphase mixture with ∼90 wt % crystals is enormous (∼1015 Pa s) and is clearly not relevant to the products of the Bandelier eruptions, which have crystal contents of <5 vol. % to a maximum of 32 vol. % in rare samples. On the other hand, the explosive Bandelier eruptions clearly did involve very large volume fractions of fluid. The fluid volume fraction is accordingly computed based on the two-phase (melt + fluid) sub-system. The high final volume fraction of fluid in the magma, ∼87 vol. % just below Tinv, has important dynamical implications (e.g. see Blake, 1984). Fluid volume fractions exceeding ∼60% correspond to the volume fraction at which magma fragmentation occurs and are consistent with pyroclastic plinian eruptions. The upward movement of magma initiates further exsolution of fluid from fluid-saturated melt, and upon ascent, a positive feedback will come into play such that the magma viscosity continues to decrease. This, in turn, increases the rate of ascent and hence the rate of decompression, and generates even further fluid exsolution, thereby leading to positive feedback.
Figure 15d shows the variation of melt density as a function of temperature. Melt densities are similar in all four cases because the mass fraction of fluid remains rather small. During temperature decrease from the liquidus, the density of melt is nearly constant, then begins to decrease as clinopyroxene and spinel join the fractionating assemblage at 1030°C and 1020°C, respectively ( ≈ 0·75). Melt density levels off at ∼900°C ( ≈ 0·42), then increases slightly as temperature decreases further, primarily reflecting the effect of the melt becoming increasingly Na2O- and water-poor. The decrease in melt density, similar to the effects of decreasing viscosity, enhances the proclivity for crystal–melt separation.
Finally, Fig. 15e shows the relationship between the volume fraction of bubbles in the magma (θ), as a function of the initial wt % H2O, melt fraction, and lithostatic pressure. Along with Fig. 15f, this figure shows dramatically why silicic large-volume explosive eruptions originate at shallow depth. At pressures equivalent to mid-crust depths (∼15 km), the volume fraction of exsolved fluid never exceeds about 15 vol. % regardless of the dissolved water content of the parental melt or the extent of crystallization. Even with extensive crystallization the volume fraction of exsolved water never rises above 15 vol. %. In contrast, at low pressure there are many initial water content–melt fraction (fm) combinations that lead to volume fractions of fluid that exceed the critical limit for fragmentation of ∼60 vol. %.
Hildreth & Wilson (2007) advocated a multistage model in which mantle-derived hydrous basalts partially crystallize in the deep crust, incorporating wall-rock anatectic melts of diverse trace element and isotope chemistry to form H2O-rich, and therefore, low-viscosity, low-density, rhyolite to basalt residual products (Annen & Sparks, 2002; Annen et al., 2006). During subsequent ascent, differences in the relative slopes of liquidii and the adiabat lead to attainment of superliquidus states and resorption of any entrained xenoliths and crystals from the deep crust. Extensive crystallization begins in the shallow crust as decompressing melt batches intersect their H2O-saturated liquidii. Associated large viscosity increases lead to stalling and generation of a quartz- and feldspar-rich dacite to rhyodacite mush many kilometres thick (Hildreth, 1981; Hildreth & Moorbath, 1988; Annen & Sparks, 2002; Annen et al., 2006; Hildreth & Wilson, 2007). Upward extraction (e.g. Bachmann & Bergantz, 2004) of fractionation products—H2O-rich, rhyolitic crystal-poor interstitial melt—via upward intrusion as dykes and pods, forms a crystal-poor, density-stratified, H2O-rich melt reservoir overlying the mush zone. Gas-saturation pressures in melt inclusion-bearing quartz crystals indicate that the melt reservoir roof was 5–6 km beneath the paleosurface and the deepest magma withdrawn trapped melt inclusions at 10–11 km (Wallace et al., 1999).
We first performed MELTS phase equilibria calculations based on numerous pre- and post-caldera basalt to dacite parental melt bulk compositions erupted between 3·8 and 2·5 Ma, from Bailey (2004). Basalt Mg-number [≡ atomic Mg/(Mg + Fe2+)] values extend from 58 to 80. Results presented here are based on a single representative parent magma composition, a pre-caldera basaltic trachyandesite from the Basin and Range Province (sample b11, Mg-number = 63; Bailey, 2004). The basalts are hypothesized to have evolved isobarically in sills injected at 20–30 km depth (Conrad discontinuity to Moho). Accordingly, we set initial H2O at its liquidus saturation concentration for pressure fixed at 0·5 and 0·75 GPa. Oxygen fugacity follows the QFM + 2 buffer (Hildreth, 1979: Bishop Tuff between QFM + 1 and QFM + 2). The solidification mode is fractional crystallization, which provides the optimal condition for generating silicic melt. We neglect assimilation of wall-rock anatectic melts. Varying fixed pressure from 0·05 to 0·75 GPa with all other conditions the same yields calculated melt compositions that are nearly identical up to residual melt silica concentrations (anhydrous) of ∼60 wt % (fm ≈ 0·75). This is because magnesian olivine and clinopyroxene plus spinel are the only stable phases. Continued fractionation at fm < 0·75 at 0·5–0·75 GPa leads to divergence in residual melt compositions and residual melts that are dissimilar to the field of Bishop Tuff major element data. At 0·5 GPa, crystallization to near-solidus conditions leads to maximum SiO2 concentrations of 67 wt %, with a predicted phase assemblage that does not correspond to the observed phase assemblage (i.e. garnet, leucite, and muscovite are predicted). Above 0·5 GPa, the discrepancies increase. We suggest that a melt silica concentration of ∼60 wt % SiO2 (∼25% fractional crystallization) represents an approximate upper limit on deep crustal fractional crystallization. That is, ≥75% or more of the crystallization must take place at shallower levels. Because fractionation to 60 wt % SiO2 generates identical liquids at all pressures (all other conditions identical) from ∼0·05 to 0·75 GPa, we employ the simplifying assumption of a single pressure along the entire liquid line of descent in all subsequent calculations of shallow-level fractionation.
According to Hildreth & Wilson (2007), superliquidus residual melts from lower crustal fractional crystallization rose to varying levels, then underwent isobaric fractional crystallization, forming an intermediate mush zone several kilometres thick. The high-silica rhyolite residual melts, transported subsequently to an overlying melt reservoir, are compositionally equivalent to the Bishop Tuff. At the initiation of mush zone fractional crystallization, we assume that there is no free water phase; we fix the initial dissolved H2O at its saturation concentration. Fluid-saturation pressures in quartz-hosted melt inclusions (Wallace et al., 1999) suggest that the deepest Bishop Tuff magma withdrawn trapped melt inclusions at 10–11 km, so we infer that the underlying mush zone in which fractionation occurred extended to 10–11 km at its shallowest levels (∼0·25 GPa). The mush zone is hypothesized to be several kilometres thick, but we present results for mush zone fractionation based on the lowest possible fixed pressure value of 0·25 GPa in all calculations; this generates the best possible correspondence between Bishop Tuff data and MELTS predictions for the Hildreth & Wilson (2007) scenario.
Significantly, maximum calculated melt silica concentrations are lower than the Bishop Tuff data field and Al2O3 concentrations are generally higher than the data field for parental compositions with less than ∼57 wt % SiO2 (Fig. 5, Electronic Appendix). There are pre-caldera dacites from Bailey (2004) that have >57 wt % SiO2 and these do produce the Bishop Tuff conmposition when fractionated at the top of the mush (0·25 GPa). We have shown, however, that if these liquids started as basalts from the mantle and fractionated to dacite in the deep crust, as Annen et al. (2006) and Hildreth & Wilson (2007) suggested, then they would not attain Bishop Tuff liquid (glass) compositions. Fractionation below the top part of the hypothesized mush zone—increasing the pressure above 0·25 GPa—only increases the discrepancy between predictions and observations. The poor overlap between predicted and observed melt compositions for shallowest-level mush zone crystallization at P ∼0·25 GPa motivates us to reject the mush model of rhyolite extraction for the Bishop Tuff. We find that fractionation at P ≤ 0·25 GPa is required, broadly consistent with the model of Anderson et al. (2000).
Bishop Tuff isobaric and isochoric MELTS crystallization calculations
The best-case MELTS models for Bishop Tuff petrogenesis are based on either isobaric or isochoric fractionation of precipitated solids and exsolved H2O, or fractional crystallization of solids only, starting with a dissolved water concentration of 3·3 wt %, oxygen fugacity fixed at the QFM + 2 buffer, and a lithostatic pressure of ∼0·12 GPa (∼5 km depth). The representative parent melt composition is a pre-caldera basaltic trachyandesite erupted between 3·8 and 2·5 Ma within the Basin and Range Province from Bailey (2004).
Figure 16 shows solid phase identities and abundances and the temperature at which water saturates as a function of temperature for cases (a) isobaric solidification, crystals removed and H2O fluid bubbles retained; (b) isochoric solidification, crystals removed but bubbles retained; (c) isobaric solidification, crystals and bubbles removed; (d) isochoric solidification, crystals and bubbles removed. In all four cases, olivine is the liquidus phase, followed by clinopyroxene, spinel, plagioclase, orthopyroxene, biotite, apatite, rhombohedral oxide, and alkali feldspar and quartz. In case (d), leucite appears prior to quartz. Similar to the Bandelier Tuff, predictions generally match observations for major phases, except for olivine, which is predicted but not observed. Allanite, zircon, and pyrrhotite are observed minor phases that are not included in the MELTS database and are therefore not predicted.
Calculated liquid major element trends reflect removal of these phases and are illustrated in Fig. 17. The four trends are similar to each other; however, in case (d), destabilization of sanidine in favour of leucite at low melt fraction leads to maximum SiO2 concentrations near the conclusion of the crystallization interval that are up to ∼2·75 wt % lower than in cases (a), (b), and (c).
Plotted in Fig. 18a are observed and predicted clinopyroxene and orthopyroxene compositions in the pyroxene quadrilateral. Bishop Tuff data, all from Hildreth & Wilson (2007), are shown as grey circles. For all cases, calculated clinopyroxene compositions are similar to each other and agree well with observed data, considering the computed trend assumes perfect, immediate removal of fractionated phases. However, the computed orthopyroxene compositional trend does not correspond well to observed phenocryst compositions, being ∼35 mol % richer in the enstatite component. The results of parameter sensitivity tests show that varying oxygen fugacity, lithostatic pressure, and parental melt dissolved water concentration does not improve the coincidence between calculated and observed orthopyroxene.
Figure 18b shows Bishop Tuff plagioclase and alkali feldspar compositions plotted along with the computed trends. Observed plagioclase crystals correspond reasonably well to the low-temperature part of the MELTS plagioclase trend, but are ∼5–7 mol % too high in dissolved KAlSi3O8. For alkali feldspar, the observed Or content varies from Or70 to Or63. Case (d) best reproduces the observed data range (Fig. 18b), with models (a), (b), and (c) extending towards the lowest-Or phenocryst compositions (∼Or60). The trend towards higher Or in case (d) is due to a dramatic pressure increase at low melt fraction (up to ∼0·34 GPa, P/Po = 2·83; Fig. 15a). Similar to the Bandelier Tuff models described above, increasing pressure yields more potassic alkali feldspar compositions. Figure 19 shows the results of isobaric crystal fractionation calculations with respect to solids and fluid, illustrating the influence of pressure on feldspar compositions. Bishop Tuff feldspars are bracketed by pressures up to 0·3 GPa. At 0·5 GPa, calculated alkali feldspar compositions approach ∼Or80.
Figure 18 also illustrates the range of observed and predicted spinel (Fig. 18c) and rhombohedral oxide (Fig. 18d) ternary compositions in FeO–MgO–Fe2O3 and FeO–MgO–TiO2 normalized coordinates, respectively. Bishop Tuff spinel and rhombohedral oxide (ilmenite) phenocrysts correspond fairly well to the low-temperature portions of the model trends.
MELTS predicts biotite saturation, consistent with observations. However, because so few experimental data are available to calibrate activity–compositions in relevant di- and tri-octahedral micas, a quantitative comparison between observed and predicted biotite compositions is not meaningful.
In summary, for most phases there is good correspondence between the identity and composition of phases observed and predicted, assuming isobaric or isochoric fractional crystallization with removal of precipitated crystals and exsolved H2O, oxygen fugacity along the QFM + 2 buffer, and a parental melt initial (dissolved) water concentration of ∼3·3 wt % H2O. At a lithostatic pressure of 0·12 GPa (depth ∼5 km), similar to 0·14 GPa fluid-saturation pressures of melt inclusions in some early erupted quartz phenocrysts (Wallace et al., 1999), the major discrepancies are the compositions of predicted versus observed orthopyroxene and feldspar. Feldspar compositions are bracketed at lithostatic pressure values up to ∼0·3 GPa (∼10 km depth) and are best described by case (d). These conditions are not inconsistent with the 5–6 km magma body thickness that is implied by maximum gas-saturation pressures of 0·26 GPa for some late-erupted melt inclusions (Wallace et al., 1999), as mentioned by Anderson et al. (2000).
Evolving melt and magma properties as a function of temperature are indicated in Fig. 15. Figure 15a shows the variation of pressure for the isochoric cases [(c) and (d)]. At the initiation of the crystallization interval, high-temperature crystallization of olivine, clinopyroxene, and spinel, all denser than melt, is associated with initial decrease of magma internal pressure to a minimum of 0·1 GPa (P/Po = 0·83). Melt saturates with respect to H2O at 1150°C ( = 0·94) in the isochoric cases [(c) and (d)] and at 1030°C ( = 0·72) in the isobaric cases [(a) and (b)]. As in the Bandelier Tuff isochoric calculations discussed above, the subsequent appearance of plagioclase marks the beginning of a pressure increase that levels off, then reverses towards the end of the crystallization interval in case (c), where H2O bubbles are retained. In case (d), pressure continues to increase, rising to a maximum of ∼0·34 GPa (P/Po = 2·83) at the lowest melt fraction.
The appearance of alkali feldspar at ∼840°C, ≈ 0·3 [cases (a), (b), and (c)] and at 810°C, = 0·27 [case (d)], is universally associated with a steepening decrease in melt fraction and an increase in the mass (and volume) fraction of supercritical fluid (Figs 15b and 16), associated with the pseudoinvariant point, Tinv, 1. A second, more pronounced pseudoinvariant point develops when quartz joins the crystallizing feldspars (plagioclase and alkali feldspar) at Tinv, 2 ≈ 700°C in all cases (Figs 15b and 16). Compared with cases (a), (b), and (c), isochoric fractional crystallization with respect to solids + fluid [case (d)] yields a lower final exsolved supercritical fluid fraction [case (d): 0·74; case (a): 0·85; case (b): 0·79; case (c): 0·87; Fig. 15b] as a result of higher magma internal pressure at low melt fractions and related higher dissolved H2O concentrations. However, all of these volume fraction values for exsolved H2O exceed the critical value for magma fragmentation. Melt viscosity remains below 10 Pa s where H2O fluid is retained [i.e. cases (a) and (c); not shown]. Expulsion of H2O fluid [cases (b) and (d)] leads to a maximum in melt viscosity at Tinv, 2 in both the isobaric and isochoric cases, but the isobaric case value is significantly higher (∼1·1 × 106 Pa s), as a result of lower dissolved H2O and Na2O concentrations, than in the isochoric case (3·6 × 104 Pa s). It is also higher than in the Bandelier calculations for the same reasons. Melt density values along the liquid line of descent (Fig. 15d) are similar to those of the Bandelier Tuff. Trends for the four cases are similar to each other. Density initially decreases with falling temperature, reflecting rising dissolved water and decreasing Na2O concentrations in the melt, then increases at the lowest temperatures as the dissolved water content decreases.
Comparison with existing models
The feasibility of crystal fractionation can be assessed independently through calculation and comparison of the timescales necessary for solidification. We use a simple thermal model (Fowler et al., 2007) based on heat-transfer scaling to define the upper and lower bounds on the time frame of development of Bishop Tuff composition melt from a basaltic trachyandesite parent composition via fractional crystallization, then compare the results with independent age data.
Approximately 1 MJ/kg of heat must be removed during isobaric crystallization to pass from the melt at the liquidus temperature to the near-solidus condition of fm = 0·05 where MELTS calculations were stopped. A timescale can be derived by making simple assumptions regarding the rate of heat extraction per unit area of magma body–wall-rock contact (), the fraction of parental melt volume (fm) that differentiates to form Bishop Tuff magma of volume VBT, and the fraction (α) of that volume that erupts to form the Bishop Tuff (VEBT). Typically α lies between unity and 0·1. By definition, VEBT = α fmV. In the scale analysis we treat the volume of the magma body that undergoes crystallization as a connected region of volume V and surface area A. Area and volume are related according to A = KV2/3 where K is a dimensionless constant that depends on the shape of the magma reservoir. For cubical, spherical and disk-like magma volumes of differing aspect ratio, K varies between about 5 and 7. A scale time, τ, for fractionation is constructed by assuming that magma of volume V and density ρ loses heat at rate and that the total amount of heat removed is the difference in the enthalpy of the magma from the starting temperature, 1160°C, to the lowest melt fraction ( fm = 0·07) temperature of 700°C. This enthalpy, , is computed from MELTS and represents the total decrease in specific enthalpy (J/kg) of the magma during crystallization from the liquidus to the lowest fraction of remaining melt, fm, based on the starting volume of the parental magma. With these assumptions, the scale time is calculated by simple heat balance as
Equation (2) may be written explicitly in terms of the time since the start of fractionation (t) according tois a function of the melt fraction; and are respectively, the specific enthalpy (J/kg) at the liquidus and the specific enthalpy as a function of fm. is zero at the liquidus temperature and reaches a maximum at low temperature, at the lowest fm. Once the parameters are set, the enthalpy of the system, which comes from the phase equilibration calculation, can be related to the absolute time since the start of cooling and crystallization. The scale time, τ, is calculated for parameters chosen to provide upper and lower scale times as a function of the DRE (Dry Rock Equivalent) eruptive volume (VEBT) of the Bishop Tuff. The parameters r, , K, and fm are set equal to 2150 kg/m3, 0·87 MJ/kg, 5, and 0·07, respectively. Conductive heat flow measurements obtained within the upper 300 m of sedimentary fill at Long Valley caldera indicate values of 0·167 W/m2 near the western rim and 0·084 W/m2 near the eastern rim of the caldera (Lachenbruch et al., 1976). We adopt these values as minimum and maximum bounds on heat flow, respectively. For comparison, the regional background heat flow is 0·063–0·084 W/m2, whereas measured heat flow at Yellowstone, where a partially molten crustal magma chamber is hypothesized to reside beneath the surface (Husen et al., 2004), is 2 W/m2 (Morgan et al., 1977; Fournier, 1989). The timescale is inversely proportional to , so using higher values for the heat loss naturally leads to smaller scale times, τ. The fraction α of differentiated magma that erupted to form the Bishop Tuff from a possibly larger subjacent volume of undifferentiated magma is difficult to estimate but probably lies in the range 0·5–1 (White et al., 2006). Here we adopt maximal and minimal values of 0·8 and 0·4. In Fig. 20a, the scale time τ is shown as a function of the DRE eruptive volume of the Bishop Tuff (VEBT) for the range of estimated parameters. Scale times range from ∼1·58 Ma to ∼3·95 Ma for a Bishop Tuff volume of VEBT = 600 km3. Using a value for τ according to this simple thermal model allows approximation of an absolute timescale for the fractionation process, and hence, a scale to estimate the age of crystals. Figure 20b shows such a scale with τ equal to 1·58 Ma based on equations (1), (2), and (3). A τ of 1·58 Ma corresponds to the maximal present-day heat flow value. The implication of Fig. 20b is that each type of phenocryst has an age range.
Such timescale predictions must agree with isotopic ages. Consistent with our results, 87Rb–87Sr mineral and glass ages and Ar apparent ages well in excess of the ∼760 ka 40Ar/39Ar eruption age (Sarna-Wojcicki et al., 2000) for the Bishop Tuff have been identified as evidence that a long-lived (>1 Myr) silicic magma chamber was the source of the Bishop Tuff (e.g. Halliday et al., 1989; Christensen & DePaolo, 1993; van den Bogaard & Schirnick, 1995; Christensen & Halliday, 1996; Davies & Halliday, 1998). Bindeman & Valley (2002) interpreted large-scale δ18O homogeneity of the Bishop Tuff as evidence of its longevity (>105 years) and convective homogenization. On the other hand, they noted that isotopic zoning in quartz, trace element gradients in feldspars, and quartz and zircon crystal size distributions are more consistent with far shorter timescales (102–104 years). They reconciled this discrepancy in terms of recharge 103–104 years before eruption of the Bishop Tuff. In addition, recent U–Pb studies of accessory phases (Reid & Coath, 2000; Simon & Reid, 2005; Crowley et al., 2007; Simon et al., 2007) have revealed crystal populations with young, narrow age ranges that have been interpreted as evidence that Long Valley rhyolites, through eruption of the Bishop Tuff, evolved in <300 kyr. Bishop Tuff sanidines erupted late in the pyroclastic sequence range in Sr isotope composition from 0·70595 to 0·70603 and yield model ages >400 kyr older than the eruption. However, the corresponding scatter in Pb isotope values extends up to three times that of analytical uncertainty, leading Simon et al. (2007) to propose that the variation in Sr isotopes reflects an isotopically heterogeneous system instead of a long crystallization interval. Those workers suggested that glass might be differentially more affected by open-system processes, so that sanidine–glass model ages overestimate the range and absolute value of crystallization ages. They called on this phenomenon to explain a 150 kyr difference between a mean feldspar model age (1000 ka) and a zircon crystallization age (850 ka) from Simon & Reid (2005). They suggested that the extent of chemical relaxation of Sr and Ba heterogeneity in zoned Bishop Tuff sanidine supports a duration of no more than a few hundred thousand years for magma accumulation (Morgan & Blake, 2006). Relative to late Bishop Tuff sanidines, early erupted feldspars exhibit clustered Pb isotope compositions and greater Sr isotope heterogeneity, with 87Sr/86Sr in some crystals extending to >0·713 (Davies & Halliday, 1998). According to Simon et al. (2007), these feldspars could have crystallized after a long period (>1 Myr) of radiogenic Sr ingrowth in melts similar to the Bishop Tuff, or they could be antecrysts derived from earlier intrusions. Simon et al. (2007) concluded that these feldspars have ambiguous significance for magma residence times because some early erupted feldspars yield negative Sr model ages when coupled with their host melt. In addition, no similarly old (>1 Ma) crystallization ages have been identified for >50 early erupted non-xenocrystic zircons. We do not accept that the U–Pb accessory phase data (Simon et al., 2007) rule out a long residence time for the Bishop Tuff magma body. Although we agree that some magma–crust interaction is likely to have occurred (we do not ruminate here on the issue of crustal contamination), we suggest that differences in feldspar model ages and accessory phase crystallization ages may arise from distinctions in their stabilization behaviour. Zircon is not included in the MELTS database, but the model ages of >400 ka for Bishop Tuff sanidines erupted later in the pyroclastic sequence are approximately equivalent to sanidine stabilization times from our thermal scaling (Fig. 20b). However, we maintain that accumulation of the Bishop Tuff magma began earlier, well before sanidine stabilized, at >1 Ma. Different phases subsequently stabilized at different times prior to eruption. Isotopic age data on other phases are needed to test these models in detail.
Hindered settling calculations (see Fowler et al., 2007) that incorporate evolving physical properties from the phase equilibria models confirm that fractionation can occur within the calculated thermal timeframe, consistent with Anderson et al. (2000).
Bindeman et al. (2008), augmented by aspects of the study by Hildreth et al. (1991), rejected the scenario of a single, large-volume, supersolidus, mushy-state Yellowstone magma chamber that existed for over 0·2 Myr, periodically producing rhyolite (as in the study by Vazquez & Reid, 2002). They proposed that the source rocks for Yellowstone silicic volcanism cooled below the solidus prior to hydrothermal alteration and remelting to liberate high-silica rhyolite melts. Radiogenic isotope trends reflect isotopic differences between unrelated batches of remelted precursor volcanic and sub-volcanic rhyolite, variably altered at shallow levels (1–5 km) by heated meteoric waters to lower δ18O. Basaltic magma pulses provided the heat for melting of these precursor rocks, but not matter. The basalts are the products of deep crustal hybridization.
This is a particularly challenging scenario to test quantitatively. Bulk melting following initial crystallization does not lead to unique element concentration signatures. Additionally, using the expected isotopic effects of post-crystallization hydrothermal alteration as a progress variable in this regard is complicated by the unknown degree to which alteration (possibly involving leaching by hydrothermal fluids of highly radiogenic Sr and Pb from surrounding Archaean country rocks) progressed. A further complication involves the potential incorporation of unknown masses of hybridized basaltic recharge melt (Bindeman et al., 2008).
We concur that isotopic data, including zircon oxygen isotope values and U–Pb age data, respectively provide incontrovertible evidence for influx of meteoric water and/or other inherited material into the Yellowstone silicic magmas; however, we do not attempt to model these processes. The rhyolite intrusions were presumably subjected to minimal major element modification during hydrothermal alteration and bulk melting. We can therefore use major element data from the large-volume tuffs to constrain phase equilibria calculations of rhyolite generation during the initial cooling of hybridized basalt below the solidus (Bindeman et al., 2008). In the light of the Bishop and Bandelier results discussed above, we suspect that examining the consequences of crystallization may yield interesting insights. We define a best-case scenario for crystallization below.
Yellowstone Tuffs isobaric and isochoric MELTS crystallization calculations
Whole-rock major element and feldspar compositional variation for the Huckleberry Ridge, Mesa Falls, and Lava Creek tuffs is limited (Fig. 21), and can therefore be described by a single MELTS model. We have not evaluated the degree of compositional overlap for solid phases other than feldspar because of a scarcity of relevant compositional data. The parental magma composition is a pre-Huckleberry Ridge basalt (sample 8yc475; Christiansen, 2001).
Figure 22 shows mineral identities and abundances and the temperature at which H2O saturates for isobaric and isochoric fractional crystallization with respect to solids and solids + fluid. Olivine is the liquidus phase (Tliq ≈ 1140°C), followed by clinopyroxene, plagioclase, and spinel. Anorthoclase, alkali feldspar, rhombohedral oxide, apatite, fayalitic olivine, and quartz are also predicted, but their order of appearance varies for the four cases. Orthopyroxene stability is predicted for cases (a) and (d); biotite and leucite appear in case (d) only. These predictions broadly accord with observations for the Yellowstone Tuffs. Amphibole is observed, but not predicted. Biotite is predicted only in case (d). Leucite is predicted in case (d), but is not observed. Zircon, chevkinite, and allanite are present within the Yellowstone Tuffs in trace quantities, but are not included in the MELTS database, and therefore are not predicted.
Plots comparing major element trends versus MgO concentration for the best-case scenario are shown in Fig. 21. The calculated trends are generally similar to each other and overlap with the observed data. Na2O and Al2O3 reflect the greatest differences in the modelled reaction paths (up to 0·75 and 2·0 wt %, respectively) towards the conclusion of crystallization (Fig. 21c and d). The largest discrepancy between predicted and observed melt compositions is for CaO in all cases (Fig. 21e).
Plotted in Fig. 23a are predicted and observed ternary feldspar compositions. For alkali feldspar, observed compositions fall into one of two groups: a band from ∼Or62 to Or50 and a less frequent group representing samples from the Lava Creek Tuff at ∼Or40. The fractional crystallization models do not reproduce the observed alkali feldspar compositions at ∼Or40. The computed trend for case (d) (isochoric fractional crystallization with removal of both solids and fluid) describes only the highest-Or end of the high-Or group. However, agreement between observed data for the higher Or group and computed trends for cases (a), (b), and (c) is good (Fig. 23a). Other observed feldspar compositions include oligoclase, andesine, and anorthoclase, all of which correspond fairly well to the low-temperature portions of the computed compositional trends for each case. Few compositional data are available for other phases (e.g. pyroxene, spinel) in the Yellowstone Tuffs. We show predicted compositions in the pyroxene quadrilateral (Fig. 23b), but only two measured data points constrain the trends. At 0·14 GPa, the samples display significantly greater enrichment in Fe (∼Fs48) compared with the model trends for each case (maximum ∼Fs33). We have found that the model trends reproduce the data at slightly lower pressure (∼0·12 GPa). Figure 23b also shows that fayalite is predicted, consistent with observations. No measured data are available to constrain calculated spinel.
In summary, the best-case lithostatic pressure is ∼0·12–0·14 GPa (∼5 km depth). Oxygen fugacity is constrained to follow the QFM – 1 buffer, and the initial dissolved H2O concentration is ∼3·5 wt %.
Figure 15a shows a comparison of pressure versus temperature during isochoric crystallization based on fractionation of solids only and solids + fluid. Here, as in the Bandelier and Bishop systems, there is an initial decrease in pressure associated with the removal at high temperature of phases that are denser than the melt. Thereafter, the pressure flattens as a result of H2O fluid exsolution at 1050°C ( = 0·8, isobaric cases) and at 1130°C ( = 0·96, isochoric cases). With decreasing temperature, the two isochoric reaction paths diverge. Retention of H2O fluid bubbles combined with crystallization of phases denser than the melt leads to pressure decrease towards the conclusion of crystallization, whereas expulsion of H2O bubbles drives up internal magma pressure to a maximum of 0·33 GPa (P/Po = 2·35) at the lowest temperature.
In Fig. 15b–d, the variation of fluid volume fraction, dissolved water content, and melt density is shown versus temperature. Like the Bandelier Tuff and Bishop Tuff systems, the onset of quartz crystallization coincides with an abrupt increase in exsolved supercritical fluid and a decrease in melt fraction. H2O fluid retention [cases (a) and (c)] leads to melt viscosity values along the liquid line of descent of <10 Pa s; however, where H2O is expelled, melt viscosity achieves peak values of 1·6 × 104 Pa s [case (d)] and 1·1 × 106 Pa s [case (b)] before falling abruptly at Tinv (not shown graphically).
Comparison with existing models
Bindeman et al. (2008) proposed that at Yellowstone basaltic magma provided the heat for remelting shallow-level, earlier-erupted, hydrothermally altered volcanic rocks, liberating rhyolite melts. Our results show that the basalts evolved inevitably towards destabilization as they solidified. Although isotopic data show that inherited components entered the magma from the surrounding wall-rocks, we maintain that fractional crystallization was the dominant process governing the compositional and dynamic fate of the pre-eruptive Yellowstone large-volume rhyolite bodies.
Upper and lower scale times calculated above for Bishop Tuff solidification exceed 1 Myr. Similarly derived timescales for the Mesa Falls Tuff and Lava Creek Tuff must satisfy constraints provided by recurrence intervals of ∼600 kyr following the Huckleberry Ridge Tuff and Mesa Falls Tuff eruptions. Again we use equations (1)–(3) with parameters chosen to provide upper and lower scale times as a function of the DRE eruptive volume of the Mesa Falls and Lava Creek tuffs. Parameters α and K are respectively set at 0·4–0·8 and 5. The starting volumes of parental magma are VEMFT = 300 km3 (Mesa Falls Tuff) and VELCT = 1000 km3 (Lava Creek Tuff). Results from phase equilibria provide values of 2304 kg/m3, 0·86 MJ/kg, and 0·05 for average melt density, ρ, , and fm, respectively. A current power output of 4·5–6·0 GW corresponds to a heat flux of ∼1·5–2·1 W/m2 over the 2900 km2 Yellowstone caldera, ∼30–40 times the continental average (Fournier, 1989; Friedman & Norton, 2007; Lowenstern & Hurwitz, 2008). We adopt these values as minimum and maximum bounds on heat flow. In Fig. 24, scale times are shown as a function of DRE eruptive volume for the range of estimated parameters. Scale times range from ∼120 to ∼200 kyr for Mesa Falls Tuff and from ∼180 to ∼310 kyr for the Lava Creek Tuff. For the Huckleberry Ridge Tuff (not constrained by preceding large-volume eruption), scale times range from ∼250 to ∼420 kyr. These results, although based on a highly simplified model that neglects smaller-volume events following the three large-volume caldera collapses, show that fractional crystallization of each large-volume Yellowstone rhyolite plausibly occurred within <600 kyr. Fractionation timescales for the Mesa Falls and Lava Creek tuffs remain within 600 kyr for heat flux values as low as 0·5–0·6 W/m2.
Bindeman & Valley (2001) and Bindeman et al. (2008) suggested that the ubiquity of zircon in Yellowstone rhyolites—some zircons are included within sanidine, quartz, and magnetite crystals—testifies to the critical importance of zircon for interpreting Yellowstone rhyolite petrogenesis. Our results show that magnetite stability persists from high to low temperature, whereas sanidine and quartz are stable only at low melt fractions during fractional crystallization. Although zircon provides valuable age and petrographic information, its apparent absence over the bulk of the crystallization interval (there is no record of its inclusion in phases other than those mentioned here) leads us to wonder whether its role has been overestimated in this case, leading to overemphasis of the importance of bulk melting. Bindeman et al. (2008) assessed the plutonic versus volcanic (buried by caldera collapse) origin of zircons inherited from remelted source rocks through Rayleigh fractionation calculations of zircon U, Th, and Y. Concentrations of 3000 ppm (U), 4000 ppm (Th), and 8000 ppm (Y) are taken as corresponding to 40–60% crystallinity, the rheological limit for melt extractability from a crystal mush (Bachmann & Bergantz, 2004). According to Bindeman et al. (2008), these threshold values demarcate zircons inherited from volcanic precursor rocks (lower U, Th, and Y than the cut-off values) and plutonic or sub-volcanic rocks (higher U, Th, and Y than the cut-off values). We cannot test these calculations because zircon is not included within the MELTS database (i.e. we cannot model its evolving proportions during differentiation). It does not necessarily have a uniform distribution or stability along the liquid line of descent. Zircon U, Th, and Y concentrations may be strongly influenced by the stability of other phases with potentially nonuniform distribution (i.e. minor phases), variable stabilization behaviour, and varying affinity for zircon. Perusal of the GERM partition coefficient database shows that the Kd for U in high-silica rhyolite allanite can extend as high as 17. For Th, high-silica rhyolite allanite Kd values are as high as 648. The ranges for clinopyroxene, magnetite, and fayalite extend from incompatible to compatible, at 0·36–2·65, 0·08–13·1, and 0·28 to 5·0, respectively. Consideration of data for lower-silica rhyolites further broadens some of these ranges. For Y, rhyolite amphibole Kd values range from 4 to 45·2 and ilmenite values vary between 0·21 and 5·09. At present there are no age data with which to compare our results on phases other than zircon, but we are working on this issue.
COMPARISON OF LARGE-VOLUME SYSTEMS
For all systems, differences in the liquid lines of descent and in the compositions of precipitated crystals are relatively small, regardless of thermodynamic constraints (isobaric versus isochoric) and mass-transfer assumptions (fractionation of solids only or fractionation of crystals and H2O). The largest differences occur near the solidus for systems that evolve through isochoric fractional crystallization with respect to solids + fluid, because the pressure near the solidus climbs by a factor of two to three over its initial value (Fig. 15a). In natural systems, host-rock failure would most probably occur before such high anomalous pressures could be attained. The minimum internal magmatic underpressure (pressure less than the lithostatic stresses at the reservoir margins; P/Po) during olivine + clinopyroxene + spinel crystallization is ∼0·8 of the absolute lithostatic pressure. Interestingly, such a low-pressure condition could trigger periods of enhanced assimilation associated with addition of stoped blocks and anatectic country rock liquids into the magma body.
For the three systems, calculated phase proportions for the three distinct crystallization paths associated with each of the parental magma compositions share many significant features. In all three cases, crystallization begins with olivine followed by clinopyroxene, then spinel or plagioclase. The order of spinel and plagioclase depends on oxygen fugacity; higher oxygen fugacity stabilizes early spinel crystallization. At lower temperatures, anorthoclase (Yellowstone Tuffs and Bandelier Tuff) and/or alkali feldspar crystallize. Finally, quartz saturates very close to the solidus. Very small proportions of other phases also crystallize. For example, apatite, orthopyroxene, rhombohedral oxide (ilmenite), and biotite crystallize in all systems under all or many thermodynamic and mass constraints, and fayalite crystallizes very close to the solidus in the Bandelier and Yellowstone Tuffs. These predictions are consistent with observations for each system.
Important conclusions are that regardless of reaction path, eutectic-like crystallization occurs in all systems near the solidus. This eutectic-like behaviour is associated with quartz saturation and, in the Bishop Tuff, with the appearance of alkali feldspar. This eutectic-like behaviour leads to an abrupt decrease in melt fraction over a small temperature interval (a few degrees) and drives each system towards magma fragmentation and hence, dynamic instability under isobaric or isochoric conditions. Although the present study considers the simplest case of closed-system fractional crystallization and does not explicitly consider assimilation, the drop in magma pressure during the early stages of isochoric FC leads to a sizeable pressure gradient that can lead to mechanical conditions that favour assimilation.
VOLUMES OF CUMULATES
The phase equilibria models presented above reasonably describe the production of the Bandelier, Bishop and Yellowstone tuffs by fractional crystallization. However, if valid, these models imply that evolution of thousands of cubic kilometers of parental melt in the upper crust resulted in crystalline products of comparable volume beneath each volcanic field. Seismic velocity structures in the crust beneath all three systems are reasonably well resolved by active seismic experiments and tomographic studies (Hill et al., 1984; Ankeny et al., 1986; Kissling, 1988; Dawson et al., 1990; Steck & Prothero, 1994; Steck et al., 1998; Miller & Smith, 1999; Aprea et al., 2002; Husen et al., 2004), but the corresponding compositional structure of the crust and upper mantle remains poorly resolved. Although a number of high-velocity, possibly mafic, bodies have been identified in the lower crust–upper mantle beneath Long Valley and the Jemez Mountains volcanic field, geophysical studies have not identified the requisite cumulate deposits, either because large-scale fractional crystallization did not occur in the uppermost crust or because of extensive thermal and compositional modification of the crust and upper mantle during prolonged, continuing magmatism that complicates investigation into the compositional characteristics of the crust.
Even assuming that each system evacuated its entire liquid volume during caldera collapse, erupted volumes in each of our calculations represent ∼5% of a parental melt volume, so that the ∼600 km3 Bandelier and Bishop tuffs correspond to parental melt volumes of ∼1·2 × 104 km3. The three major Yellowstone tuffs of total volume 3·8 × 103 km3 suggest a total parental melt volume of 7·6 × 104 km3. The dimensions of the bodies implied by our models seem formidable, but their existence cannot be ruled out. Current regional topographic uplift and extension related to magmatic activity is associated with the Jemez Mountains volcanic field and Long Valley region, and is particularly well known at Yellowstone. The 0·6 km high, 600 km wide topographic swell centred on Yellowstone caldera is widely interpreted as being supported by a buoyant volume of material in the crust or mantle (Smith & Braile, 1994). Estimating the volume of material supported by the anomalous topography is complicated by the 16 Myr history of basaltic–silicic magmatism over 800 km along the Snake River Plain to Yellowstone. For example, the present bulge overlaps with the voluminous (4 × 103 km3), next oldest (∼6·5–4·3 Ma) Heise silicic volcanic field to the SW along the Snake River Plain (Morgan & McIntosh, 2005). A topographic swell is likely to have persisted throughout the history of Snake River Plain–Yellowstone volcanism (Parsons et al., 1994; Smith & Braile, 1994); its previous dimensions and the exact proportion of melt within the swell and the volume of erupted material during the history of Snake River Plain–Yellowstone volcanism are not well known. Nevertheless, its current proportions can provide a first-order constraint on the parental melt volume for the 2·05 Ma to present Yellowstone silicic volcanic field that can be compared with our phase equilibria results. A cone of radius 300 km and height 0·6 km has a volume of ∼5·7 × 104 km3. A cylinder with the same radius and height has a volume of ∼1·7 × 105 km3. Accommodation of a parental melt volume of 7·6 × 104 km3 implied by MELTS is not inconceivable then. At Long Valley, the picture is simpler. A topographic swell 100 km wide, currently centred on Mono Basin, has been attributed to magmatic activity (Hill, 2006). Assuming a feature of similar diameter has persisted for the duration of Long Valley magmatic activity, the cylindrical and conical volumes constrained vertically by the height of the Mono–Inyo Craters volcanic chain above the surrounding terrain (610 m; Bailey, 1989) are respectively ∼4·8 × 103 km3 and ∼1·6 × 103 km3. If we use the maximum topographic height of Long Valley caldera (3·6 × 103 m) above the surrounding terrain (∼2·2 × 103 m) as the vertical constraint (1·4 × 103 m), the corresponding conical volume is ∼3670 km3 and the cylindrical volume is ∼1·1 × 104 km3. Again, accommodation of a significant portion of parental magma is possible. There is mention of an epeirogenically uplifted topographic swell (Morgan & Swanberg, 1985; Eaton, 1986, 2008; Morgan, 2003) in the northern New Mexico region, but its dimensions are unknown.
We have sketched a metamodel that encompasses potential scenarios for crustal magmatic systems as a foundation for quantitative testing. The metamodel encapsulates the potential complexity of magmatic system evolution into the temporal sequence of events (e.g. fractional crystallization, assimilation, magma recharge, etc.) and concomitant prevailing environmental conditions (oxygen fugacity, pressure, fugacity of H2O, parental melt composition) during those events. This approach is distinct from complicated ‘just-so’ scenarios, each one specific to a particular system. The difference in these two approaches is that the former recognizes the unity of natural systems whereas the latter focuses upon the distinctiveness of each. Both approaches are useful in building towards a comprehensive understanding of the behaviour of natural systems. Whereas the former may miss important details the latter can miss the forest for the trees. The metamodel further serves as a framework for posing quantitative tests of petrogenetic hypotheses, such as those proposed for the ∼1·2 and ∼1·6 Ma Bandelier Tuff, the ∼0·7 Ma Bishop Tuff, and the ∼2·1, 1·3, and 0·64 Ma Yellowstone Tuffs. The most recent published scenarios for each system are associated with distinct aspects of the metamodel, forming a continuum that reflects the long-standing debate regarding the importance of crustal melting versus fractional crystallization as the essential mechanism of magma generation. At one extreme, the erupted material is derived from pre-existing crust. At the other extreme, the erupted material is derived mostly from crystallization of the parental magmas. A feature shared by all of the published scenarios is that solidification in the upper crust leads to an uneruptible mush. In the approach of Bindeman et al. (2008) (Yellowstone Tuffs), basaltic recharge leads to bulk remelting of the solidified mush and intra-caldera volcanic rocks (right side of Fig. 1). In that of Rowe et al. (2007) (Bandelier Tuff), recharge partially remelts young mafic intrusions that are solid, but still warm (right side of Fig. 1). In the approach of Hildreth & Wilson (2007) (Bishop Tuff), recharge occurs, but interstitial silicic melt separation takes place through extraction from the mush.
We evaluated each scenario through systematic comparison of predictions of multicomponent–multiphase equilibria models subject to varying thermodynamic constraints (isobaric or isochoric, fractionation or retention of fluid and/or crystals). For each system, parental melts are pre-caldera basalts that are spatially related to each large-volume tuff deposit. The results lead to the following conclusions.
Fractional crystallization can account for the bulk of variability in both liquid and solid phase compositions for each system. Although we do not claim that the match between predictions and observations is perfect, we do maintain that fractional crystallization is the dominant mechanism.
Predicted liquid lines of descent, phase proportions, and solid phase compositions for models based on either isobaric or isochoric crystallization, with either retention or fractionation of exsolved fluid, are similar. Although in general phase relations are not grossly dependent upon whether fractional crystallization occurs isobarically or isochorically, the magma pressure in strictly isochoric crystallization undergoes an interesting evolution. Initially, in isochoric crystallization the mean magma pressure falls below that of the lithostat. This condition could lend itself to magma contamination by the action of stoped block calving or by the percolation of wall-rock anatectic melts into the adjoining magma body. Either of these processes would probably be of limited extent in most cases, although one can envisage scenarios in which they could be quantitatively as important as fractional crystallization. Physical models can be constructed to quantitatively evaluate the importance of these assimilation mechanisms.
High volatile concentrations (Anderson, 1976) and shallow magma reservoirs—at ∼0·1–0·3 GPa (∼4–10 km depth; e.g. Bachmann & Bergantz, 2008)—are widely associated with large-volume explosive eruptions. Accordingly, we find for each system that parental melts saturated at the liquidus with ∼3·2–3·5 wt % H2O evolved by isobaric or isochoric fractional crystallization, mostly at lithostatic pressures of ∼0·1–0·14 GPa (∼4–6 km depth). The phase equilibria models allow up to 20% by mass of crystallization within the deep crust. This is because the nature and composition of the solid phases dominant at high melt fraction (olivine and clinopyroxene) vary little with varying lithostatic pressure. The nature and composition of other phases that stabilize at >20% crystallization do vary with varying lithostatic pressure. We expect that magma stagnating at depth or at low initial water concentration would freeze beneath the surface.
The consequence of evolution at shallow levels by fractional crystallization of a water-rich parental magma is dramatic variation along the liquid line of descent in magma physical properties, leading to inherent destabilization and explosive giant eruptions. A second stage of processing, as proposed in existing scenarios, is not required. In each system an invariant point temperature corresponds to a jump in supercritical fluid fraction accompanying the appearance of quartz and/or feldspar to values that exceed recognized fragmentation limits. In some systems (e.g. Campanian Ignimbrite, Fowler et al., 2007), the onset of quartz and/or feldspar stability corresponds to a dramatic increase in crystallinity. In the systems discussed here, the increase in crystallinity is less marked.
In isochoric models, the pressure of melt can increase considerably compared with the lithostatic pressure, and can exceed the tensile strength of the country-rock, leading to failure by crack propagation and concomitant magma decompression. The pressure increase is greatest (P/Po ≥ 2·35 for all systems) in the case that relatively compressible supercritical fluid is expelled from the melt as it forms. Retention of compressible supercritical fluid bubbles dampens the pressure increase. In any case, for initially H2O-rich magma hosted at shallow levels by weak or strong country-rocks, where melt and fluid are separated or remain in equilibrium, explosive eruption is an inevitable outcome of solidification. Several researchers have proposed that recharge events can rejuvenate crystal-rich mushes that are otherwise too viscous to flow or erupt (e.g. Murphy et al., 2000; Bachmann et al., 2002). Adding fresh (hot) recharge magma can remelt crystals and increase the magma volume, leading to over-pressurization, the mechanical failure of wall-rocks, and the triggering of an eruption. We point out that evidence of basaltic and andesitic recharge-driven eruption is difficult to test because mixing during magma withdrawal from a stratified body will juxtapose blobs of magma initially at very different depths and different initial temperatures (Spera, 1984; Trial et al., 1992). In any case, volumetrically significant masses of basaltic recharge are not observed in the systems discussed here.
We have calculated thermal timescales for magmatic evolution based on scale analysis, given rates of heat loss constrained by geothermal observations, the liquidus to solidus enthalpy difference, and the volume of the system. For the Bishop Tuff, reasonable values yield a timescale of fractional crystallization for the Bishop magma body of >1 Myr. This analysis can be applied to single crystal phases to provide estimates of durations of crystallization and therefore potential age ranges for discrete crystal populations. For example, a minimum estimate of the onset of sanidine stability is 400 kyr, a value that coincides with model ages from Simon et al. (2007). Those workers interpreted a difference in mean ages for sanidine and zircon crystals as evidence for open-system processes during Bishop Tuff petrogenesis. We recognize that radiogenic isotope data provide some evidence for assimilation of wall-rocks, but another source of variation is that crystals from a single magma body may record different stages of evolution. If our model for Bishop Tuff petrogenesis is valid, zircon, which is apparently included in solids stable at the lowest melt fractions (spinel, sanidine, and quartz), may provide information only about the very last stages of Bishop Tuff evolution. For the Yellowstone Tuffs, calculated thermal timescales are consistent with recurrence intervals of ∼600 kyr between successive caldera collapses.
Supplementary data for this paper are available at Journal of Petrology online.
We are very grateful to editor Marjorie Wilson and reviewers Catherine Annen and an anonymous reviewer for their constructive reviews. This work was supported by grants from the National Science Foundation (EAR-0838182 and ATM-0425059) and US Department of Energy (Geosciences) to F.J.S.