The longevity of giant magma bodies in the Earth’s crust prior to eruption is poorly constrained, but recognition of short time scales by multiple methods suggests that the accumulation and eruption of these giant bodies may occur rapidly. We describe a new method that uses textures of quartz-hosted melt inclusions, determined using quantitative three-dimensional propagation phase-contrast X-ray tomography, to estimate quartz crystallization times and growth rates, and we compare the results to those from Ti diffusion profiles. We investigate three large-volume, high-silica rhyolite eruptions: the 240 ka Ohakuri-Mamaku and 26.5 ka Oruanui (Taupo Volcanic Zone, New Zealand), and the 760 ka Bishop Tuff (California, USA). Our results show that (1) longevity estimates from melt inclusion textures and Ti diffusion profiles are comparable, (2) quartz growth rates average ∼10−12 m/s, and (3) quartz melt inclusions give decadal to centennial time scales, revealing that giant magma bodies can develop over notably short historical time scales.


Supereruption deposits suggest that giant pools of magma (>450 km3; Sparks et al., 2005) have accumulated in the crust and erupted repeatedly throughout Earth’s history. Supereruptions are potentially devastating to life and infrastructure, and understanding these magmatic systems is critical to our ability to prepare for such events. The time scales over which giant magma bodies accumulate is crucial to constrain and has been a topic of considerable study; however, no consensus has been reached, and results vary widely with the method and/or mineral phase used (e.g., Bishop Tuff, California, USA, zircon geochronology: 103–106 a; Crowley et al., 2007; Reid and Coath, 2000; oxygen isotopes: 104 a; Bindeman and Valley, 2002; thermal modeling: 103–105 a; Gualda et al., 2012a; Gelman et al., 2013; textures and diffusion chronometry: 103 a; Gualda et al., 2012a). This discrepancy exists in part because different timekeepers in volcanic deposits record different time scales, and the resolution of these time scales varies by method. We introduce a new method to estimate time that is based on textures (sizes, shapes, and positions) of melt inclusions hosted in quartz crystals. The power of this approach is that it elucidates the duration that a magma existed as an eruptible melt-rich body, rather than a crystal-rich uneruptible one. Ultimately, to establish successful volcanic hazards assessment, monitoring, and mitigation programs, it is imperative that the former time scales be constrained.

We (1) assess quartz crystallization times and growth rates from melt inclusion textures, derived from quantitative analysis of X-ray tomography of quartz crystals, for the paired Ohakuri-Mamaku (Taupo Volcanic Zone, New Zealand, ∼100 km3; Gravley et al., 2007; Bégué et al., 2014), Oruanui (Taupo Volcanic Zone, New Zealand, ∼530 km3; Wilson, 2001), and Bishop Tuff (∼1000 km3; Wilson and Hildreth, 1997) supereruptions; (2) compare times and growth rates from melt inclusions with those from Ti diffusion chronometry in the same crystals; and (3) assess the significance of the results in regard to the final accumulation of giant magma bodies in the crust, using rhyolite-MELTS modeling (Gualda et al., 2012b) to evaluate the completeness of the time record preserved in quartz.


As magmas evolve, crystallization occurs and round blebs of melt (liquid) can be trapped in crystals (melt inclusions). Round melt inclusions are disequilibrium features, and equilibrium is approached as diffusion acts to fully facet them (i.e., develop the negative crystal shape of the host mineral; Beddoe-Stephens et al., 1983; Chaigneau et al., 1980; Manley, 1996; Gualda et al. 2012a); this is a spontaneous process, as the free energy reduction in making curved surfaces flat more than offsets the rise from increasing the total surface area of the inclusion (Manley, 1996; Wortis, 1988). The ability for such a disequilibrium feature to be preserved in erupted material is a function of kinetics: at magmatic temperatures, diffusion of many elements is fast relative to geologic time, but it effectively stops upon eruption and cooling. Thus, melt inclusions record passing time through shape maturation, changing from round to faceted, and the extent to which they are faceted can be used to quantitatively estimate crystallization time scales and, more broadly, the longevity of a melt-rich magma body.

Quartz is a particularly useful mineral for such analysis because it is a major phase in many silicic rocks, and the faceting rate is limited only by the diffusion rate of Si in melt, for which diffusion coefficients at relevant temperatures have been determined experimentally (Baker, 1991). The underlying theory and equations for melt inclusion faceting in quartz were developed and implemented by Gualda et al. (2012a) to assess the longevity of the Bishop Tuff; however, a limited number of inclusions were studied, and shapes were characterized only qualitatively, using two-dimensional (2-D) sections or projections from optical microscopy. Here we employ 3-D propagation phase-contrast X-ray tomography to image quartz crystals and their inclusions, which provides quantitative metrics for their textures (Pamukcu et al., 2013). We also use these images to document inclusion positions within crystals, which allows us to calculate crystal growth rates, a poorly constrained parameter for mineral crystallization in magmas.

To validate the faceting method, we compare results from Ohakuri-Mamaku melt inclusions to those determined using Ti diffusion chronometry at zone boundaries near or containing the selected melt inclusions. Diffusion chronometry also relies on kinetic markers of time: a zone boundary represents a concentration gradient, and the severity of this feature is lessened by diffusion. Thus, the slope of a compositional profile across a zone boundary gradually becomes shallower. This process stops at eruption, so the steepness of this slope today can be used to estimate the boundary residence time (i.e., the time a zone boundary spent at magmatic conditions). Quartz is also an excellent phase for diffusion chronometry (e.g., Wark et al., 2007; Matthews et al., 2012): zoning is primarily due to Ti (Spear and Wark, 2009; Leeman et al., 2012), and the diffusion rate of Ti in quartz is isotropic, independent of composition or the presence of other elements, and has been measured experimentally (Cherniak et al., 2007).


Melt Inclusion Faceting and Diffusional Relaxation Times

We calculate melt inclusion faceting times assuming that (1) faceting is a constant-volume process (i.e., total inclusion volume does not change); (2) transport is limited by diffusion through the melt; and (3) inclusions are initially round and the equilibrium shape is a hexagonal bipyramid. These are conservative assumptions that lead to maximum estimates for the inferred time scales; if grain boundary diffusion prevailed (potentially much faster than volume diffusion; Dohmen and Milke, 2010) and/or the initial and equilibrium shapes were less extreme (i.e., partially faceted initially, fully faceted with rounded corners), our calculated times would be shorter.

While Gualda et al. (2012a) calculated the complete (full) faceting time for a given inclusion, we calculate the time represented by the current shape of the inclusion (determined using tomography; see the GSA Data Repository1), which may not be fully faceted. Ultimately, we calculate the time that was required to diffuse the volume of the inclusion that currently protrudes from a best-fit ellipsoid that represents the initial shape of the inclusion (Fig. 1). Faceting time is primarily controlled by the volume of diffused material (ΔV) and Si diffusivity in the melt (Fig. DR1 in the Data Repository), the main control of which is temperature. Fully faceted inclusions record minimum times because an inclusion cannot change shape after it is fully faceted; thus, partially faceted inclusions provide more meaningful estimates, and we use somewhat conservative parameters (see Gualda et al., 2012a) to yield maximum faceting time estimates. We assess how faceted an inclusion is by comparing the time required to achieve the current shape with that required for full faceting (Fig. 2; Fig. DR1).

Cathodoluminescence (CL) intensity has been suggested as a good proxy for Ti contents in quartz (Cherniak et al., 2007; Spear and Wark, 2009; Gualda et al., 2012a; Leeman et al., 2012), and we verified this in a subset of our crystals by comparing CL images with Ti maps from X-ray microprobe (see the Data Repository). We find that this correlation holds for our samples (Fig. DR2), at least for the more pronounced variations in CL intensity relevant for diffusion chronometry (i.e., some fine-scale CL zoning is unrelated to Ti); thus, we calculate times on profiles derived from CL images. We apply a 1-D diffusion model to these profiles and find the best-fit error function, from which the diffusional length scale (L) is estimated (following Gualda et al., 2012a). Time (t) is then calculated using the relationship t = DTiL2/4, where DTi is the diffusion coefficient of Ti in quartz. We assume that profiles are initially step functions, so resultant times are maximum estimates.

Errors on our calculated times are 125% for faceting and 200% for diffusion chronometry (Gualda et al., 2012a). Importantly, although we use a single temperature for these calculations (see the Data Repository), we assume uncertainties of ±30 °C (2σ) in our error estimation (Gualda et al., 2012a), which encompasses the range of reasonable temperatures over which crystallization may take place in rhyolitic magmas. Furthermore, these systems are nearly invariant and characterized by feldspar saturation prior to quartz (Gualda and Ghiorso, 2013), so we expect modest temperature changes of much less than the 60 °C we assume here (Gualda et al., 2012b).

Quartz Growth Rates

We use the simple relationship of Growth Rate = Length/Time to estimate average quartz growth rates. Time is that estimated from melt inclusion faceting or diffusion chronometry. Length for melt inclusion–based estimates is determined during image processing: host quartz crystals can be processed in the same way as inclusions, and melt inclusions can be repositioned inside the host. The distance from the centroid of the inclusion polyhedron to the nearest face in the crystal polyhedron can then be determined (Fig. 1). For zoning-based estimates, Length is obtained by measuring the distance from the relevant zone boundary to the nearest edge of the crystal.


To compare results from faceting and diffusion chronometry, we studied a subset of Ohakuri-Mamaku crystals using both methods (Fig. 2; Fig. DR3; see the Data Repository). Time scale and growth rate estimates from faceting and diffusion chronometry are comparable within and between individual crystals: both approaches suggest that quartz in these magma bodies grew over short time scales (∼101–102 a; Figs. 2 and 3; Fig. DR4) and that average growth rates are ∼10−12 m/s (Fig. 1; Fig. DR3). The consistency of the results from the two approaches suggests that the faceting method is robust.

Given this, we consider faceting times for all three eruptions in more detail (Fig. 3). Most inclusions in Ohakuri-Mamaku crystals are partially faceted, giving time scales on the order of ∼101–102 a. In both Oruanui and Bishop crystals, a number of small inclusions are fully faceted, but many larger ones are partially faceted and give time scales of ∼102 a. Bishop inclusions give the longest times (<600 a), followed by Oruanui (all but a few give <300 a). Ohakuri-Mamaku inclusions give considerably shorter times, with most clustering at <50 a; however, most of the partially faceted inclusions from Bishop and Oruanui give times of <100 a (Fig. 3). Growth rates are also impressively similar in all three systems, with a near-normal distribution peaking at ∼10–11.5 to 10–12.5 m/s (Fig. 1). The consistency of the results, both between samples and methods, suggests that an average growth rate of ∼10−12 m/s is appropriate for quartz in various large to giant systems. These are one to two orders of magnitude faster than those determined previously for Bishop Tuff quartz using diffusion chronometry (Gualda et al., 2012a), but those estimates are likely limited by the lower spatial resolution of the CL images that were used.

In the three systems we studied, the time scales from melt inclusions are strikingly short; most suggest times of 101–102 a, up to 103 a. (Fig. 3). Considering the errors on these times, which are conservative, maximum estimates from the inclusions do not exceed 2 ka. This is consistent with crystal textures: using our most conservative growth rate (10–12.5 m/s; Fig. 1), a quartz crystal 2 mm in diameter would grow in ∼100 a. Even allowing for 100% error on the growth rates (errors on measured growth distances are likely much smaller than those on crystallization times, so the latter dominate the errors on calculated growth rates), growth times for a 2 mm crystal would be <200 a. Such large crystals are not common in these rocks, so this result is commensurate with the suggestion of decadal to centennial time scales for the crystallization of these bodies.

In interpreting these results, three important points should be noted:

(1) The system may have waxed and waned over time, such that the calculated times accrued over a longer period of time. Using a combination of geochronology and diffusion chronometry, Cooper and Kent (2014) demonstrated that crystals may only spend a fraction of their total growth history in melt-rich magmas (i.e., significant time may be spent at or near subsolidus conditions); thus, faceting or diffusional relaxation may take place over multiple time intervals. We emphasize that, for eutectoid systems like these, quartz would likely be fully resorbed during remobilization events, so quartz should lack long memory (in contrast to other minerals, like zircon and even feldspars). Furthermore, we focus on the longevity of the magma bodies that feed giant eruptions (i.e., the residence time of large pools of crystal-poor magma in the upper crust since the final remobilization events that resulted in their accumulation). If prior faceting or diffusional relaxation has taken place, our calculated time scales are overestimates of this residence time.

(2) Quartz may not be in the liquidus assemblage for these magmas. Quartz is unlikely to be the first phase to crystallize in silicic magmas (Gualda and Ghiorso, 2013), unless it is part of the liquidus assemblage with other minerals. Assuming that magma solidification is limited by heat loss to the crust, we can use heat budgets derived from rhyolite-MELTS to assess the fraction of the total crystallization time that is recorded by quartz (Gualda and Ghiorso, 2011; see the Data Repository). For Oruanui, simulations on average whole-rock compositions from phase 10 (final eruptive phase; Wilson, 2001) pumice clasts show that ∼7 wt% crystals (primarily plagioclase) would have crystallized before the onset of quartz crystallization (Fig. 4). Assuming that heat loss to the crust is relatively constant, the change in enthalpy of the system (ΔH) between the onset and end of crystallization can be used as a relative measure of time. In this case, a total ΔH of –42 J/g is released between plagioclase and quartz saturation, and another –14 J/g is released until 13 wt% crystals (maximum crystallinity in Oruanui pumice; Wilson, 2001) is achieved (Fig. 4). This suggests that quartz times may underestimate the longevity of the erupted Oruanui magma by no more than a factor of three; late-erupted Bishop Tuff simulations are similar (factor of two; Gualda et al., 2012b). Early-erupted Bishop Tuff and Ohakuri-Mamaku simulations suggest these magmas were saturated in quartz and feldspars from the outset (Wark and Watson, 2006; Gualda et al., 2012b); thus, quartz times represent the entire duration of magmatic crystallization in these cases. Consequently, we conclude that all systems studied here exhibited a final remobilization stage that led to the accumulation of large volumes of magma that resided and crystallized in the shallow upper crust for a few decades to (at most) centuries prior to eruption.

(3) Our current data set supports a number of melt extraction and assembly models, from multiple discrete magma batches (Cashman and Giordano, 2014) to single giant pools (e.g., Wilson et al., 2006), as well as a patchwork of discrete very large bodies (Gualda and Ghiorso, 2013; Bégué et al., 2014).

These results are in stark contrast to those from some studies that suggest that giant crystal-poor silicic magma bodies accumulate over 104–106 a (Gelman et al., 2013; Reid and Coath, 2000; Bindeman and Valley, 2002). Others have proposed much shorter time scales (102–103 a; Crowley et al., 2007; Gualda et al., 2012a; Wilson et al., 2006; Allan et al., 2013), and our results are consistent with these; however, our results also suggest that even these shorter time scales may be maximum bounds, with final accumulation and crystallization of the melt-rich magma body more likely taking place on decadal to centennial time scales.

We emphasize that these short time scales can be reconciled with many of the longer estimates, and the discrepancy may be merely a matter of perspective: longer time scales, particularly those recorded in long-lived, resilient zircon crystals, likely represent the protracted history of mushy root zone development, solidification and partial melting events, and priming of the crust that allows for the eventual accumulation of giant melt-rich magma bodies (Cooper and Kent, 2014; Gualda et al., 2012a). Meanwhile, the short, decadal to centennial time scales derived here reflect the final accrual of buoyant, unstable, crystal-poor magma parcels in the shallow crust.

Portions of this work were performed at GeoSoilEnviroCARS (Sector 13), Advanced Photon Source (APS), Argonne National Laboratory (Illinois, USA). GeoSoilEnviroCARS is supported by the National Science Foundation (NSF)–Earth Sciences (EAR-1128799) and U.S. Department of Energy–GeoSciences (DE-FG02-94ER14466). This research used resources of the Advanced Photon Source, a U.S. Department of Energy (DOE) Office of Science User Facility operated for the DOE Office of Science by Argonne National Laboratory under Contract No. DE-AC02-06CH11357. Funding was provided by NSF grant EAPSI-1209584 to Pamukcu, and NSF grants EAR-1321806, EAR-1151337, EAR-0948528, and two Vanderbilt Discovery Grants to Gualda.

1GSA Data Repository item 2015316, expanded methods, Figure DR1 (schematic of faceting time plots, times for individual eruptive units and different temperatures), Figure DR2 (Ti map versus CL image, plots of beam intensity versus CL profiles), Figure DR3 (CL images, measurements, and calculations for Ohakuri-Mamaku samples), Figure DR4 (effects of melt inclusions on CL zoning), and Table DR1 (calculations for Oruanui and Bishop inclusions), is available online at www.geosociety.org/pubs/ft2015.htm, or on request from editing@geosociety.org or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.