Recent discoveries of isotopically diverse minerals, i.e., zircons, quartz, and feldspars, in large-volume ignimbrites and smaller lavas from the Snake River Plain (SRP; Idaho, USA), Iceland, Kamchatka Peninsula, and other environments suggest that this phenomenon characterizes many silicic units studied by in situ methods. This observation leads to the need for new models of silicic magma petrogenesis that involve double or triple recycling of zircon-saturated rocks. Initial partial melts are produced in small quantities in which zircons and other minerals undergo solution reprecipitation and inherit isotopic signatures of the immediate environment of the host magma batch. Next, isotopically diverse polythermal magma batches with inherited crystals merge together into larger volume magma bodies, where they mix and then erupt. Concave-up and polymodal crystal size distributions of zircons and quartz observed in large-volume ignimbrites may be explained by two or three episodes of solution and reprecipitation. Hafnium isotope diversity in zircons demonstrates variable mixing of crustal melts and mantle-derived silicic differentiates. The low δ18O values of magmas with δ18O-diverse zircons indicate that magma generation happens by remelting of variably hydrothermally altered, and thus diverse in δ18O, protoliths from which the host magma batch, minute or voluminous, inherited low-δ18O values. This also indicates that the processes that generate zircon diversity happen at shallow depths of a few kilometers, where meteoric water can circulate at large water/rock ratios to imprint low δ18O values on the protolith. We further review newly emerging isotopic evidence of diverse zircons and their appearance at the end of the magmatic evolution of many long-lived large-volume silicic centers in the SRP and elsewhere, evidence indicating that the genesis of rhyolites by recycling their sometimes hydrothermally altered subsolidus predecessors may be a common evolutionary trend for many rhyolites worldwide, especially in hotspot and rift environments with high magma and heat fluxes. Next, we use thermomechanical finite element modeling of rhyolite genesis and to explain (1) the formation of magma batches in stress fields by dike capture or deflection as a function of underpressurization and overpressurization, respectively; (2) the merging of neighboring magma batches together via four related mechanisms: melting through the screen rock and melt zone expansion, brittle failure of a separating screen of rocks, buoyant merging of magmas, and explosive merging by an overpressurized interstitial fluid phase (heated meteoric water); and (3) mixing time scales and their efficacies on extended horizontal scales, as expressed by marker method particle tracking. The envisioned advective thermomechanical mechanisms of magma segregation in the upper crust may characterize periods of increased basaltic output from the mantle, leading to increased silicic melt production, but may also serve as analogues for magma chambers made of dispersed magma batches. Although not the focus of this work, dispersed magma batches may be stable in the long term, but their coalescence creates ephemeral, short-lived eruptable magma bodies that erupt nearly completely.


The origin of silicic magmas is a traditional topic of igneous petrology that has fundamental significance for the origin of the first continental crust (Hutton, 1794). It was recognized early that fractional crystallization of basaltic magmas could only generate small, difficult to extract, interstitial silicic liquids (Bowen, 1928; Tuttle and Bowen, 1958), recently advocated for the Snake River Plain and Yellowstone (western USA) plume by McCurry et al. (2008) and Whitaker et al. (2008). An alternative mechanism, partial melting of the lower basaltic crust (Leeman et al., 1992; Dufek and Bergantz, 2005) has problems related to heat availability, conduction, and the large amounts of basaltic magma required, often producing excessive amounts of associated cumulates (e.g., Beard and Lofgren, 1991; Petford and Gallagher, 2001; Ducea, 2001). Models involving filter pressing, or reactivations of melts from cumulate piles and static or dynamic mushes, have advantages (e.g., conservation of latent heat and space), but also significant limitations on the physical mechanisms of melt extraction and the slow rates of these processes (Wickham, 1987; Philpotts et al., 1996; Bachmann and Bergantz, 2004; Huber et al., 2009; Streck and Grunder, 2008; Bindeman et al., 2008b).

Many factors that militate against the production of silicic magmas, including slow extraction, the room problem, and heat, are eliminated by a hypothesis of simple recycling of originally silicic or intermediate protoliths, which may have been formed by any or all of these slow and incremental processes. If rhyolites are thought to be produced by such recycling, estimates of the total amount of silicic magma generated in a given magmatically active area, the amount of basalt required, and the amount of heat required become significantly smaller.

A large-volume silicic eruption provides a snapshot of an instant in the temporal evolution of a magma body; repeated eruptions over millions of years give insights into petrogenesis at depth (Lipman, 2007). However, the debate over whether plutons and ignimbrites provide comparable records of magma genesis is ongoing, especially with reference to arc and collisional environments, where batholiths may be assembled over millions of years (Coleman et al., 2004; Glazner et al., 2004; de Silva and Gosnold, 2007) and magma bodies are kept near solidus.

The study of the plume-related Snake River Plain (SRP; Idaho, USA) silicic volcanism (Fig. 1) has an advantage over other areas of voluminous silicic volcanism because the North American plate is moving at a speed of 2–3 cm/yr over the Yellowstone hot magmatic zone (e.g., the hotspot and/or propagating rift; Pierce and Morgan, 2009), creating a conveyor-belt record of silicic magma genesis as opposed to more stationary examples in convergent margins and subduction arcs, where volcanism occurs in a localized area and younger magmatic episodes obscure their predecessors. At least six major nested caldera clusters have developed over the past ∼15 m.y. (Fig. 1); current activity is centered at Yellowstone. Each of the volcanic centers is the site of several major caldera-forming eruptions and overlapping calderas, although individual caldera boundaries are often obscured in the older 15–11 Ma centers.

Recent work (Figs. 1 and 2) has demonstrated that Picabo, Heise, and Yellowstone, the three most recent caldera clusters, share common features of magmatic evolution for the past 10.5 m.y. Eruption of the first major tuff occurred after an ∼2 m.y. hiatus from the last major eruption in the previous caldera center. The initial Huckleberry Ridge Tuff of Yellowstone erupted at 2.1 Ma, after the 4.4 Ma Kilgore tuff, which finished activity at the Heise center to the southwest, and the oldest Heise tuff, the 6.6 Ma Blacktail Creek tuff, erupted ∼2 m.y. after the eruption of the 8.3 Ma Pocatello tuff, the youngest tuff of the Picabo center (Fig. 1) (e.g., Christiansen et al., 2007; Morgan and McIntosh, 2005; Drew et al., 2013). These relationships become more obscure in the older and less-well-exposed central SRP centers (Fig. 1), including the Twin Falls–Bruneau-Jarbidge centers, which also exhibited lingering silicic magmatism as young as 8–6 Ma (see Ellis et al., 2013, and references therein), which were also likely influenced by coeval rifting (Konstantinou et al., 2013), and are now buried by Quaternary SRP basalts. Ongoing investigation of the earliest ca. 15 Ma volcanism before the Bruneau-Jarbidge center (Colon et al., 2013), however, demonstrates a gap of ∼1–2 m.y. between scattered post–Columbia River Basalt volcanism and more focused central SRP volcanism. The ∼2 m.y. gap is sufficiently long to generate the first voluminous silicic magma body at the beginning of each caldera cluster cycle, comparable to the Huckleberry Ridge Tuff batholith, ∼5000–10,000 km3 in volume, assuming an erupted versus residual magma ratio of 1:3. To produce these rhyolites, we therefore favor slow incremental assembly of magma by fractional crystallization of SRP basalt and melting of the lower crust in a hot zone (via sill emplacement; Annen and Sparks, 2002), with insignificant contribution from the upper, highly radiogenic crust, as is indicated by the modestly radiogenic values of rhyolites (Leeman et al., 1992; McCurry and Rodgers, 2009), or more rarely with an exclusively upper crustal source (e.g., the Arbon Valley Tuff Member; Drew et al., 2013). Usual or slightly elevated hotspot basaltic intrusion rates of 0.001–0.01 km3/yr (Pierce and Morgan, 2009) are sufficient to explain the required volumes of silicic magma to build the first magma bodies that are normal δ18O and have 60%–88% mantle contribution to their radiogenic isotopes (Leeman et al., 1992; McCurry and Rodgers, 2009; Nash et al., 2006). The rest of this paper deals with the recycling of these first magma bodies as a chief mechanism of rhyolite petrogenesis.


The advent of microanalytical isotopic techniques, i.e., single crystal laser fluorination, isotopic dilution-thermal ionization mass spectrometry, single zircon analyses with high precision, and ion microprobe analysis with worse precision but superb spatial resolution, allows a novel look at igneous petrogenesis. The rapid and accelerating rate of discoveries of isotopically diverse zircons, quartz, feldspars, and other minerals in large-volume ignimbrites and smaller volume lavas worldwide suggests that this phenomenon characterizes many silicic units studied by in situ methods. We here summarize recent results of these investigations with emphasis on the SRP.

Zircon and Quartz δ18O Diversity

Early evidence of zircon and quartz δ18O zoning and diversity in a few Yellowstone lavas (Bindeman and Valley, 2001) has been observed in more than half of units studied in the Columbia River Basalt–SRP large igneous province (Figs. 1 and 3), and similar evidence is being obtained for other areas (Fig. 4). This diversity variably characterizes zircon δ18O, Hf isotopes, and sometimes older U-Pb ages from nearly half of studied SRP tuffs and lavas erupted since ca. 15 Ma: Yellowstone (Bindeman et al., 2008b; Watts et al., 2012; Vazquez et al., 2009), Heise (Watts et al., 2011), Picabo (Fig. 3; Drew et al., 2013), Bruneau-Jarbidge and Twin Falls (Cathey et al., 2011), and the early post–Columbia River Basalt volcanism (Colon et al., 2013). The diversity of the refractory mineral zircon, the mineral with the slowest oxygen diffusion rates, is not surprising, but diversity in δ18O quartz crystals is also present (Fig. 3).

In the SRP, SW Nevada, and in Iceland (Figs. 3 and 4) there is common association of diverse zircon and quartz crystals with low δ18O magmas, indicating their inheritance during remelting of buried and hydrothermally altered volcanic rocks and their subvolcanic units (e.g., Carley et al., 2011). However, there are also diverse zircons in several normal δ18O to moderately high δ18O rhyolites of the SRP, Iceland, Kamchatka Peninsula, and Toba (Fig. 4). This suggests that nearly all studied low δ18O magmas worldwide, and some nominally normal δ18O magmas, display δ18O diverse zircons and sometimes other minerals. Zircons are in equilibrium only in the Fish Canyon tuff, likely a reflection of long residence of this magma body (Bachmann and Bergantz, 2004). The isotopic fingerprinting of the upper crust by low δ18O meteoric waters values thus serves as a tracer of the processes of remelting and recycling of hydrothermally altered rhyolites. If hydrothermal alteration happens at the seafloor using 0‰ seawater, a spread from low to high δ18O values will result (Grimes et al., 2013; Wanless et al., 2010; Bindeman et al., 2012), potentially leading to high and low δ18O magma batches. Finally, advocated processes may occur in nature even if the precursor rocks have not been hydrothermally altered with respect to δ18O, particularly in the lower and middle crust.

High δ18O, Normal δ18O, and Abundant Low δ18O Magmas in the SRP

Throughout the past decade there has been rapid discovery of low δ18O magmas in the SRP and worldwide (e.g., Bindeman et al., 2007, 2012; Boroughs et al., 2005; Ellis et al., 2010; Drew et al., 2013), and these rocks have gained increasing attention as probes into upper crustal petrogenesis. Nearly all nested caldera clusters in the SRP have produced thousands of cubic kilometers of rhyolites depleted in δ18O (the volume of these magmas in the central SRP is 10,000–30,000 km3; Fig. 1; Bonnichsen et al., 2008; Ellis et al., 2010, 2013), and individual units are as voluminous as 1800 km3 (Morgan and McIntosh, 2005). We advocate that the origin of these magmas requires a two- to three-step process: (1) hydrothermal alteration of the protolith by low δ18O meteoric waters at large water/rock ratios; (2) remelting and convective mixing during which low and variable in δ18O oxygen is homogenized in a newly formed magma batch; and (3) formation of new crystals representing batch-averaged δ18O composition, and repetition of 1 and 2. Washington State University researchers favor the idea that step 1 significantly predates step 2 (e.g., alteration preceding central SRP volcanism occurred during formation of the Eocene Idaho batholith; Boroughs et al., 2005; Ellis et al., 2013), while the University of Oregon group favors auto development of magmas from normal to low δ18O values immediately before low δ18O eruptions via successive caldera collapse and rift burial (Bindeman et al., 2007; Watts et al., 2011; Drew et al., 2013). Shallow crustal fingerprinting by low δ18O meteoric water may also be accomplished via tectonically induced water-rock interaction during syngenetic rifting, core complex formation, and alteration along normal faults (Konstantinou et al., 2013; Drew et al., 2013; Colon et al., 2013). Conditions for remelting are optimized when the rocks are brought down closer to the heat source during crustal burial processes of caldera collapses, rifting, or both. Rifting and burial due to extension may have created additional hydrogeological conditions for hydrothermal fluid circulation particularly in Iceland, where burial and alteration of rocks were followed by remelting (Martin and Sigmarsson, 2010; Bindeman et al., 2012; Figs. 3 and 4). It seems highly unlikely that low δ18O magmas are produced in any other way than by remelting of hydrothermally altered low δ18O protoliths (for a review of old models, see Taylor, 1974, 1986).

The remelting process happens in shallow depths available to meteoric water percolation, thus requiring a shallow heat source to generate thousands of cubic kilometers of magma in the upper crust. It further reaffirms, however, that even if remelting occurs effectively in the shallow crust, magmatic recycling in the lower and middle crust should be even more prevalent, given higher temperatures and melt-segregation efficiency (Annen et al., 2006; Dufek and Bergantz, 2005). Furthermore, this process must be relatively rapid, faster than mineral-diffusive and recrystallization time scales or currently the most precise geochronologic method, ID-TIMS U-Pb dating (102-104 yr; Wotzlaw et al., 2014; Rivera et al., 2014). Both of these statements contradict the currently prevailing view that silicic magma are generated slowly (because it takes effort; see above), and in the lower crust (where it is hotter), or via filter-pressing extraction from the mushy batholiths.

Insights from Zircon Hf Isotopes and Whole-Rock Sr-Nd-Hf-Pb Isotopes

In situ radiogenic isotope analyses provide further insights. The Pb and Sr isotope diversity of individual feldspar crystals within single units (Watts et al., 2012) fingerprint the isotopic sources of Yellowstone lavas. Diverse Sr and Pb isotopes in individual feldspars are a well-known phenomenon worldwide (Davidson et al., 2001), suggesting that many crystals are xenocrysts and antecrysts.

The Hf isotopes within individual zircons record crust/mantle melt proportions of each magma batch from which the zircons crystallized (Figs. 3 and 5). In two extreme cases when a silicic differentiate of the SRP crystallizes high (mantle like) εHf zircon mixes with a pure partial melt of silicic crust (ultralow εHf) zircons, the resulting final magma may have extremely diverse εHf zircons (see Fig. 3B). While the evidence for diverse Hf isotopes in zircons is growing, preliminary data suggest that both homogeneity and heterogeneity of Hf isotopes in zircons is common in SRP magmas (Picabo: Drew et al., 2013; Yellowstone: Stelten et al., 2013). Nonetheless, sporadic ultralow-Hf zircons (Fig. 3; down to –47 εHf units; Drew et al., 2013) suggest addition of melts of Archean crust during the late stages of magma segregation; in most cases ultralow εHf zircons have reset U-Pb ages, suggesting that remelts of the Archean crust appeared and then crystallized new zircons before mixing with other magmas.

Isotopically Diverse Crystalline Cargo: Inheriting Isotopic Variations from Protolith

Gaining insights into the range of isotopic values and their spatial distribution within a protolith undergoing melting are important for our future modeling of the efficiency of melting and mixing and the mechanisms of incorporation of these parameters into zircon and other minerals (Fig. 4).

Ranges of δ18O and εHf Variations in Protolith

Ranges. Given the –14‰ to –19‰ δ18O value of meteoric water, it is theoretically possible to generate a volume of rock (or a vein) with comparable δ18O values if the temperature of exchange were >500 °C, meteoric water was isotopically unshifted at these depths and temperature, the water/rock ratio was high, and/or the isotope exchange was 100%. Investigation of hydrothermally altered intracaldera tuffs and ore deposits around the world suggests that such conditions are rarely realized in nature. Isotopic ranges in and around shallow plutons, subvolcanic intrusions (e.g., Taylor, 1974, 1986), and intracaldera tuffs (McConnell et al., 1997; John et al., 2011) demonstrate higher and narrower δ18O ranges (–2‰ to +7‰). Thus, the lack of observed –18‰ to –14‰ values in zircons and quartz (Fig. 3) is partly a reflection of incomplete water-rock isotope exchange of their protolith, shifted waters, and a lower temperature of hydrothermal alteration. Another argument is purely hydrologic; in order to expect a high water/rock ratio, high porosity and permeability (open cracks) are required, and this is best achieved above the brittle-ductile transition (e.g., Norton et al., 1984). For silicic volcanic rocks these temperatures are likely to be below ∼400–450 °C (Simpson, 1985). Given Δ18Orock-water fractionation (e.g., Cole and Charraborthy, 2001), rocks are several permil higher in δ18O than equilibrium water at ∼400 °C, further reducing the lowest possible δ18O rock value in hydrothermally altered rock.

Spatial δ18O distribution. Does the protolith typically contain a frequently spaced, centimeter- to meter-scale network of isotopic variations such as ultralow δ18O veins in a normal δ18O matrix, or ultralow εHf Archean crustal slivers in high εHf matrix? The isotopic variations that we observe in zircons and other minerals will then have a highly localized (centimeter to meter scale) origin. Partial melting of such areas will generate lenses of low-δ18O melt as is observed inside granodiorite blocks of the climactic Mt. Mazama eruption (e.g. Bacon et al., 1989) or in remelted Icelandic crustal xenoliths (Gurenko et al., 2013). If such high-amplitude variations span across wide regions undergoing melting (e.g., 8‰ δ18O variation over >40 km for the Kilgore tuff; Watts et al., 2011), this will simplify not only the efficiency of subsequent mixing, but may also help with mass and heat balance problems of extraction of small-volume, zircon-saturated melt batches with diverse zircons. This extraction may be facilitated by dikes and sills generated by external tectonic forces (e.g., Norton et al., 1984). Natural evidence, however, suggests that high-amplitude variations with limited isotopic magnitudes of several permil are also accompanied by larger scale (10–103 m) variations across greater isotopic magnitudes, but the detailed structure of this variation still needs to be studied by isotope mapping. Mapping of bullseye alteration patterns over plutons (e.g., Taylor, 1974, 1986), remote sensing of alteration products, and isotopic investigation of alteration in intracaldera tuffs (McConnell et al., 1997; Livo et al., 2007; John et al., 2011) provide concentric patterns of δ18O variation (10–103 m scale). Moreover, numerical modeling of meteoric fluid flow through porous media in one dimension (Bowman et al., 1994) and two dimensions (Norton and Taylor, 1979) suggests predominance of continuous isotopic fronts with diffusive-hydraulic dispersion around fast pathways, isotopically shifted waters, and overlapping, time-integrated fronts leading to laterally smoothed δ18O distribution on a scale of meters to kilometers (Hayba and Ingebritsen, 1997).

Isotopic Reequilibration During Remelting: Lessons from Skaergaard

Figures 3 and 4 demonstrate the large isotopic diversity of zircons; however, among thousands of zircons studied by in situ methods for U-Pb age, only a few older xenocrysts survive. This suggests that any partial melts of older crust were heated above zircon saturation (e.g., Fig. 5). In all of the subsequent melting history, as shown in Figure 6, all processes lead to averaging out of isotopic differences, causing solution reprecipitation of zircons and other minerals. Prolonged residence leads to oxygen isotope equilibrium between coexisting minerals, and further homogenization of hydrothermally induced heterogeneity.

An important insight into the partial remelting of hydrothermally altered rocks with and without zircons is provided by the Skaergaard intrusion in East Greenland (Taylor and Forester, 1979; Norton and Taylor, 1979; Bindeman et al., 2008a). Following crystallization and hydrothermal alteration of the Skaergaard intrusion, the mafic Basistoppen sill had intruded 100 m above the hydrothermally altered siliceous Sandwich Horizon. The melting of Sandwich Horizon at near eutectic compositions generated partial melt inside of the matrix of rock of moderately diverse δ18O (Bindeman et al., 2008a; Wotzlaw et al., 2012). These low-degree (<20%) silicic partial melts are difficult to extract; however, minerals in them become isotopically equilibrated and zircons partially reset their U-Pb ages and acquire oxygen and hafnium isotopes reflecting each minute magma batch. In the remelted Sandwich horizon of Skaergaard, isotope diversity of zircons occurs on the 10–102 m scale (Bindeman et al., 2008a), a scale we use for modeling. The Skaergaard example is evidence that there are no extreme δ18O values in the partial melts of variably hydrothermally altered rocks, as isotopic differences induced by heterogeneous isotopic alteration (Taylor and Forester, 1979) are further averaged out by melting process in line with the reasoning above (Figs. 5 and 6B; stage III). The Sandwich horizon recrystallized 20 k.y. later (Wotzlaw et al., 2012), without extraction of this partial melt, which only underwent limited intercumulus fluid flow (Tegner et al., 2009). If the Skaergaard intrusion had another sill intrusion soon after Basistoppen, and/or if the Basistoppen sill intruded closer to the Sandwich horizon, it could have resulted in greater degree of partial melting approaching the extraction window of ∼50% melt (Marsh, 1981; Wickham, 1987; Dufek and Bachmann, 2010; stage V in Fig. 6B) and would have segregated silicic melt. Furthermore, if the sill intruded directly under the Sandwich Horizon, partial melting of this laterally confined eutectic area would have been more efficient, and large fraction partial melt extraction would have been thermomechanically driven by the viscous convection in the lower layer. This process is not static and involves melting through the separating rocks by hot magma pulses (e.g., Simakin and Bindeman, 2012), mixing and diking between magma batches, and segregation by buoyancy; these processes are thermomechanically modeled in the following. This set of processes is called bulk cannibalization (Bindeman, 2008; Watts et al., 2011) and is similar to the reactive bulk assimilation model of Beard et al. (2005) and Bacon and Lowenstern (2005).

Batch Assembly, Preeruptive Mixing of Melts, and Diverse Crystalline Cargo Leading to Diverse Crystals

The advanced stages of segregation of 10–100 m magma batches lead to further averaging out of isotope diversity of magma, but bring together isotopically-diverse crystal cargo with diverse zircons. Because diffusion of oxygen and Hf is the slowest in zircons, the zircons provide the best memory effect of the described multistage segregation, while other minerals undergo greater reequilibration. Most important, the majority of zircons are out of isotopic equilibrium with the final melt they are erupted in (Figs. 3 and 4), which suggests that a preeruptive mixing episode was efficient and isotopically homogenized the melt.

Textural Evidence of Recycling in Crystals

Abundant bright cathodoluminescence rims in quartz and zircon provide evidence of recycling of as much as half of the crystal mass due to remelting processes, and the time between melting and mixing on the batholithic scale is measured in centuries (Gualda et al., 2010; Pamukcu and Gualda, 2010). The observation by Poldervaart (1956) that zircons in granites are have log-normal size distributions is supported by investigations of zircons and other phenocryst phases in single pumice clasts from major tuff units (Fig. 7; Bindeman, 2003). In Simakin and Bindeman (2008) we analytically computed that the abundant concave-down crystal size distribution of zircons and quartz may be explained by only 1–3 remelting episodes that recycle the outermost 10%–30% of quartz or zircon radii (Fig. 7). Delayed nucleation leads to concave-up distribution and often a lack of the smallest zircon and quartz crystal sizes.

Time Scales as Recorded by Crystals

Isotopic diversity in zircons occurs within error of U-Pb ion microprobe age dating (Watts et al., 2011; Drew et al., 2013) or even ID-TIMS age dating with the precision of several thousand years (Wotzlaw et al., 2014; Rivera et al., 2014), but mineral-diffusive time scales suggest even shorter time scales of mixing. Because oxygen diffusion in quartz is >104 times more rapid than in zircons (Watson and Cherniak, 1997; Cole and Chakraborty, 2001), the preservation of isotopically zoned crystals in many units suggests short times (e.g., hundreds of years) between segregation, mixing, and eruption. In many other units quartz is isotopically annealed due to residence times >∼100 yr (e.g., Bindeman and Valley, 2001; Bindeman et al., 2006). Furthermore, Cathey and Nash (2004) and Ellis and Wolff (2012) presented evidence for compositionally (Fe-Mg) polymodal pyroxene crystals in the Cougar Point Tuffs (Table 1), which are largely unzoned and record very high (900–950 °C) temperatures. Fe-Mg interdiffusion in pyroxenes is even more rapid than oxygen diffusion in quartz, and pyroxene crystals are typically small; the occurrence of such crystals can help further constrain time scales of extraction from the hot, once recycled, low δ18O crystal batches or mushes that underwent very rapid preeruptive mixing.

Tasks for Numerical Modeling

The microanalytical isotopic evidence mentioned here directs the numerical efforts. We aim to quantify our above observations on (1) efficacy of melting and formation of magma batches, (2) efficacy of merging together of initially separated batches, (3) efficacy of mixing by convection and other processes involved in merging the magma bodies with diverse crystal cargos, and (4) the time scales of these processes, as well as the conditions when magma mixing is interrupted and magma bodies erupt.


For this work we used proprietary written and relatively simple two-dimensional finite element code to model of convection of incompressible fluid at low Reynolds numbers (plug flow approximation), valid for small magma sills with vertical dimensions of as much as ∼30 m. Additional calculations were performed on larger grids covering 100–1000 m to analyze the mechanics of large-scale batch assembly, from small magma batches through slow viscous deformation of the frame rocks; numerical details are described in the following. We also performed several more roof melting + mixing experiments (started in Simakin and Bindeman, 2012; technical details provided therein). Pressure was eliminated from the Navier-Stokes equations with the penalty formulation (King et al., 1990). Material transport is modeled by Lagrange markers with 10× finer resolution in each spatial dimension. The typical grid was 180 ×180. The code was tested with the analytical solution for the convection in the liquid with sharp viscosity contrasts (Trybitsyn at al., 2006). Further technical details and material properties of magmas and country rocks are listed in Table 2.


Formation of Larger Magma Bodies

Melting and generation of small magma batches for hotspot and rift environments were numerically modeled in Simakin and Bindeman (2012) and continued in this work (Fig. 8; see Supplemental Files 11, 22, 33). We model multistage generation of eruptable high-silica rhyolitic magma bodies (Fig. 5). The first is generation of silicic melts in the hot zone near the Moho or lower crust (e.g., Annen et al., 2006) by differentiation of basalt and melting the lower mafic crust; these dry and near-liquidus silicic melts rise up rapidly, carrying heat and undergoing further adiabatic decompression and overheating in pressure-temperature-X compositional space. They are then intruded into the upper crustal silicic lithologies defining rheological or density traps, i.e., hydrothermally altered portions of the silicic batholith or under low δ18O silicic caldera fill (Figs. 5 and 6). Formation of these small, ephemeral magma bodies proceeds rapidly, in which roof-rock melting of isotopically distinct protolith is accompanied by mixing via convective melting; the roof collapses, leading to inheritance of partial melts and isotopically diverse crystalline cargo, in which δ18O and εHf in zircons survive the longest (e.g., Fig. 2). The melting processes of silicic, near eutectic, often glassy rocks are incredibly rapid and energy efficient (Fig. 8; Simakin and Bindeman, 2012), and serve as an energy sink, with minimal heat dissipation to the environment.

We thus consider several first-order thermomechanical mechanisms of growth and differentiation of larger volume magma bodies in the upper crust, assembled from these melts, and those that were produced in situ. These shallow magma bodies should be at least as large and as hot as the volume and temperature of the observed erupted ignimbrites, i.e., 102–103 km3, but likely more than twice that given that some residual magma will remain in the crust.

Expansion of the Melting Zone

With addition of heat and mass from below, the upper crustal magma body may simply grow by expansion of the melting zone if magma addition is faster than cooling. Continuous supply of new hot silicic differentiates of basalt to the base leads to rapid upward expansion of the melt zone and sinking of separating screen rocks. High versus low basaltic power outputs (Lipman et al., 1972; Hildreth, 1981; Lipman, 2007; de Silva, 1989) determine particular temperature-time trends and the subsequent lifetimes of each of these magma bodies, pertinent to the isotopic reequilibration of inherited crystalline cargo (Figs. 5 and 6).

Merging by Melting Through Separating Screen Rocks

We envision that areas of melting and shapes of magma bodies are initially complicated (e.g., Figs. 5 and 6) and may be originally determined by country-rock composition, areas of higher glass content (so less latent heat is needed), and hydrated areas of hydrothermal preconditioning. In addition, hydrothermal redistribution of silica and alkalies will create local silica oversaturation conditions, favoring lower temperature eutectic melting in selected areas. Furthermore, screens of higher melting temperature rocks may separate low-temperature eutectic rhyolite lithologies. Figure 8 and Supplemental File 2 (see footnote 2) demonstrate such a scenario of amalgamation by two melt zones via the melting through of separate screen rocks with different silica contents. Stoping of these screen areas will lead to the release of gravitational energy, enhancing convective mixing. The thermomechanics of convection play a more important role than silica oversaturation (as also observed in Simakin and Bindeman, 2012).

Physics of Batch Assembly in Stress and Density Fields

We consider simple numerical models with two vertically separated and offset magma sills (Figs. 9 and 10 and Supplemental Files 44, 55, and 66) and explore conditions of when and how these sills can merge, and exchange and mix their crystalline cargoes. In the following we employ scaling analyses, numerical experiments, and consider relevant physical aspects of large magma body assembly from batches in a stress field.

Formation of Magma Sills and Batches

Each magma batch forming in situ by melting, or those arriving from below into the zone of accumulation, can mechanically interact with neighboring batches if they are close enough to each other, or with a large dominant magma reservoir. This interaction occurs through the thermomechanical influence of the stress field around each magma sill that controls the generation and direction of dike propagation (Fig. 9A). While there are multiple studies of the interacting vertical dikes (e.g., Takada, 1994a, 1994b), there are limited analytical or numerical simulations of the mechanical interaction of a magma reservoir and a rising dike. Field observations (e.g., Horsman et al., 2009) from some shallow silicic intrusions of ∼1 km in diameter document that when small magma batches are emplaced successively, magma bodies internally inflate then spread laterally, forming a series of silicic sills emplaced on top of each other, creating voluminous silicic stocks.

Subsequent analysis of the interaction between magma batches and a rising dike often employs consideration of the principle stress (σ1) directions around a magma reservoir, defining approximate trajectories for diking via hydraulic fracturing of rocks (Jellinek and DePaolo, 2003; Simakin and Ghassemi, 2010; Karlstrom et al., 2012). This analysis leads to the observation that for the overpressurized deep spherical magma reservoir (equally inelastic, viscoelastic, or viscous media), σ1 directions are such that the greatest compressive stress is radial to the magma body and σ3 is tangential (Fig. 9). For the underpressurized magma reservoir, principal stress directions are transposed and σ1 becomes tangential. Our numerical model shows (Fig. 9B) that an overpressurized magma body attracts a dike, causing it to propagate radially and coalesce with the magma body, thus causing the magma body to grow.

When an underpressurized magma body deflects a dike, it may produce a separate sill (Fig. 9B) capable of growing into a separate, independent magma body with a continuing supply of magma. Pressure lower than lithostatic may develop in a magma body with a low supply flux during the solidification of magma due to the volume effect of crystallization and escape of exsolved fluid phases, and due to extensional tectonics. Numerical simulations of viscoelastic elliptic chambers (see Simakin and Ghassemi, 2010) demonstrated that at magma pressures lower than lithostatic pressures, expressed through the bottom-parallel σ1 stress directions, the rock envelope presses into the chamber and undergoes radial extension. These arguments support accumulation of newly arriving batches of magma into a single growing magma body in sites of active magma generation and high intrusion rate. At reduced magma supply rates, crystallization and therefore underpressurization of magma will promote formation of separate magma batches, but underpressurization conditions may even exist at some high supply rates (Gregg et al., 2013). An additional factor promoting separate magma batch formation may be preexisting oblique faults that serve as structural traps for rising magmas and relieve stress from overpressured magma bodies (Fig. 9B; Simakin and Ghassemi, 2010). This is particularly relevant for Basin and Range faults and extension along the SRP (Fig. 1), or any other silicic magma generation in rift zones.

Merging of Vertically Stacked Magma Batches

Sills that were closely spatially emplaced are separated by bridges (or screens) of country rocks, or crystal mush if emplacement happens into a cooling and crystallizing batholith (Figs. 10–13). If certain mechanical conditions of interaction are met, these bridges can fail, causing sills to merge to form bigger magma bodies (e.g., Schofield et al., 2010, 2012a, 2012b).

We consider two basic controls on the mechanical interaction between sills: overpressure (Fig. 10), and the density contrast between magma and country rocks (Figs. 11–13). In the modeling approach used here, when slow motion of country rocks is simulated, dimensions of computational domain increase to tens of kilometers at the same grid size. In our simulations the viscosity contrast between deforming rocks and magma was set to 104.5, similar to the maximum viscosity contrast 105 used by Bittner and Schmeling (1995) for the modeling of the Rayleigh-Taylor instability of the melting front in the underplated silicic crust; this value is less than the real ratio between hot country rocks (1017–1018 Pa·s) and silicic magma (105–106 Pa·s). However, further increase in the viscosity ratio (above 104–105) has no effect on the numerical solution for slow deformation of rocks. It allows us to portray large-scale slow features (this does not require accounting for the inertia term in the Navier-Stokes equation) instead of modeling the fast processes with the 10 cm scale resolution used for melting processes. Consideration of a far greater viscosity ratio or equally high Rayleigh number (Ra) requires a supercomputer approach (Longo et al., 2012) for implementation of the inertia terms and integration of the small and large time and spatial scales, but does not change our outcome for slowly deforming rocks.

We also simplify our analysis also by neglecting complex and fast thermally activated processes inside sills, such as convection, melting, and crystallization, thus concentrating only on the viscous (slow) deformation of the rocks.

Viscoelastic materials (Maxwell viscoelasticity) behave elastically at fast deformation and flow as viscous material at slow deformation (Gerya, 2009; Kaus and Podladchikov, 2006). Scaling analysis defines the Debora number, De = η/Gt0, as a measure of viscoelastic behavior, in which η/G is the characteristic time of viscous relaxation of the deviatoric stresses, G is the shear modulus, and t0 is the time scale of the process. When De < 1, mechanical evolution of the viscoelastic solid can be approximated by viscous flow. The viscous flow approximation of the deforming crust is likely the correct case, as is suggested by geodetic data on surface deformation above an overpressurized magma chamber, data that include the volume of the rocks in the thermal aureole (Newman et al., 2006). Newman et al. (2006) suggested that the threshold viscosity for liquid behavior of the thermal aureole is 3 × 1015 Pa·s (G = 5 GPa), and De ≤ 0.26 for several months (global positioning system and interferometric synthetic aperture radar observations). The viscous approximation used in the following is therefore most likely valid for slow deformation of the rocks around magma sills below the brittle-ductile transition depth, typically at 400–450 °C.

Effects of Overpressure

We investigated the effects of overpressure on a pair of vertically stacked sills (Fig. 10), with the deeper sill undergoing gradual infill at a constant volumetric rate. For this purpose we use a vertical conduit of rectangular shape (Fig. 10) and apply a parabolic vertical magma ascent velocity (with zero values at the walls) at the inlet. This sill inflates with time, inducing viscous stress. At a given filling rate, a quasi-steady-state stress field in the country rocks around the magma sill arises at the very beginning of deformation, and changes slowly with the expansion of the sill, often leading to dike initiation. When the level of stress, a product of the second invariant of strain rate and viscosity, exceeds the yield stress of country rocks (Fig. 10), the screen rock between magma batches fails. Shear strain localization is facilitated by means of dynamic power law (DPL) rheology (proposed by Sleep, 2002, and developed by Simakin and Ghassemi, 2005, and Simakin, 2014, and references therein) to produce sharp, fracture-like shear zones (Fig. 10). The DPL model treats viscosity as a dynamic parameter accounting for the effect of the texture evolution and described by a system of ordinary differential equations on the material Lagrange particles. As an asymptotic limit at a fixed strain rate (ė),the DPL model yields an ordinary power-law viscosity law [η = a(T)/ėm] where T is temperature, a is a some temperature-dependent parameter, and m is an exponent coefficient used in description of power laws (e.g., Becker, 2006).

Calculations were performed in nondimensional form. Since density is uniform at the overpressure approximation (i.e., we consider magma overpressure over the lithostatic), the gravity term can be dropped by setting the Rayleigh Number, Ra, to 0. Stresses and viscosity are scaled with some pressure (P0) and time scales of deformation (t0): 
In these calculations nondimensional viscosities of the frame rocks are kept at forumla and magma forumla with yield stresses at forumla. In the following text we will use symbol forumla for the scaled yielding stress. Then, for example, at P0 = 1 MPa, the physical value of forumla will be 50–200 MPa. Time scales can be found by setting rock viscosity forumla. With the relatively low viscosity of the country rocks ηr = 3.0 1018 Pa·s, close to the value estimated for the thermal aureole of the magma chamber underlying the Long Valley Caldera (Newman et al., 2006), and pressure fixed at P0 = 1 MPa we get a time scale for deformation t0 = 1018.5–6-5 = 107.5 s, or ∼1 yr.

In our experiments the normalized sill expansion rate dV/V0/dτ was ∼0.01 or 0.01/t0 = 0.01 yr–1 for the above parameters. By increasing P0 we proportionally increase the rock yield strength forumla and the physical expansion rate of the magma. In other words, these nondimensional results can be interpreted equally with the choice of a high filling rate and high forumla, or a low filling rate and small Sy, as these are linearly interrelated. The linear scale, l0, remains arbitrary. Possible limits of l0 are constrained by the absolute intrusion rate (km3/yr). Results of the numerical experiments with varying Sy and a constant filling rate (0.01) show that bridge failure occurs at Sy < 150 (150 MPa at P0 = 1 MPa). This confirms that the order of magnitude scaling estimate of the viscous stress in the bridge is correct, since 150 MPa is close to the typical shear strength of the rocks (Barton, 2013).

Overpressure defines the amplitude of the induced stress, and for the expanded two-dimensional (2-D) circular shell around a magma body it was shown to be (e.g., Jellinek and DePaolo, 2003; Simakin and Ghassemi, 2010): 
where R1 and R2 are internal and external (R1 < R2) radii of the 2-D circular cavity filled with liquid. The external radius in these approximations is defined as an elliptic envelope of deformation, which coincides with the brittle-ductile transition (Newman et al., 2006). In the case of two cavities next to each other, we can expect the level of stress in the separating screen rock to be: 
where σb is the stress component normal to the boundary (radial for circle inclusion) of the expanded sill, and α is some geometric factor <1, which depends on the particular geometry of the system; α can be determined analytically or numerically (Simakin and Ghassemi, 2010). For a circular cavity in the shell α = 1 – (R1/R2)2 or (R1/R2)2 = 0.5, which corresponds to equal volumes (areas in 2-D) of magma and rocks in the thermal aureole, α = 0.5. For an elliptic chamber and viscoelastic rocks the numerically estimated value of this coefficient is 0.32 (recalculated from Simakin and Ghassemi, 2010). By setting α = 1/3, P0 = 1MPa, and dV/dτ/V0 = 0.01, we get an estimate of 333 MPa for the viscous stress in the bridge.

These numerical experiments (Fig. 10) demonstrate how the process of thin bridge failure between closely spaced magma lenses proceeds in 2-D. In particular, they show that the bridge is initially faulted at one end at the very beginning of deformation and then rotates upward onto and into the upper sill as a handle or trapdoor, eventually giving way for lower magma to flow into the upper sill. The thicker bridge between magma lenses yields two fractures, leading to the extrusion of a part of the bridge in a piston-like fashion into the upper magma body (see Fig. 10); we call this mechanism an internal caldera collapse. Our analysis also demonstrates that widely spaced magma batches can coexist for a long time without interaction, and homogenize their crystalline cargo and whole-rock chemical and isotopic composition on time scales of cooling. For example, Wilson and Charlier (2009) demonstrated coeval, probably lateral, reservoirs separated by ∼20 km at the Taupo volcanic zone (New Zealand). The numerical results here were obtained exclusively due to an overpressure response in the viscous media without considering density terms, which we consider next.

Influence of Magma Buoyancy

When sills are intruded into sufficiently soft country rocks or into a crystal mush, the density difference between the magma and surrounding media leads to slow deformations due to buoyancy. In the set of stacked magma lenses (sills), there are two possible mechanisms for buoyant rise: (1) the large-scale diapiric rise of one or two magma bodies, thinning and disrupting the screen rocks between, causing the merging of these into a single magma body; and (2) development of Rayleigh-Taylor instabilities of the plastic screen rocks separating the sills. Our numerical experiments (Figs. 11–13) show that diapir-like motion always prevails in the set of the stacked sills, as opposed to Rayleigh-Taylor instability of the separating screen rock. Diapir-like motions are driven by horizontal density gradients. For diapiric batch rise, the maximum flow rate during approximately one nondimensional unit of time varies slowly, and the vertical position of the screen midpoint rises linearly above the ascending flow of the lower sill (Fig. 12B). For sufficiently weak screen rocks, viscous stresses overstep the shear strength (yield stress in the frame of the viscoplastic rheology), resulting in sharp, fracture-like shear zones splitting the screen rock and surrounding rocks, thus facilitating further deformation.

Numerical models of the buoyant rise can be further rationalized with simple scaling analysis. With the nondimensional Navier Stokes equation, we use time scale t0, length scale l0 (size of the element grid), and viscosity scale η0. Pressure (stress) scale is derived as P0 = η0/t0. With choice of the time scale equal to t0 = η0/Δρgl0 Rayleigh number (multiplier in the scaled gravity term) becomes unity and P0 = Δρgl0. If Ra is set larger than unity (for practical considerations), then t0 = t0Ra. Thus gravity flow develops identically with the fixed initial geometry for the various density contrasts and contrasts in viscosities and linear sizes, but on different time scales. The characteristic diapir velocity for a sill semi-width W can be compared with corresponding Stokes rise velocity Us of the sphere (3-D) or circle (2-D): 
where for a solid sphere kS = 2/9, for a solid circle kS = 1/4, and for a gas 3-D bubble kS = 1/3; all values are for slow flow at small Re numbers. Our numerical experiment is not identical to any of these geometrically idealized models. While tentatively setting kS = 1/3, for the half-width of the sill, W1/2, we rearrange the equation for Us in nondimensional form as: 
Numerical nondimensional maximum flow velocity at the beginning of the experiment (ηr0 = 105, Ra = 120, forumla) can be approximated with W1/2/l0 = 21.2 while the circle inscribing both sills has the nondimensional radius 32. Because low-density material fills only part of this large circle, correspondence seems to be satisfactory.

Dimensional velocity of the diapiric flow (Ud) strongly depends on the physical parameters. In a crystal mush environment with at ηr = 1013 Pa·s, l0 = 40 m, Δρ = 300 kg/m3, and geometry as in Figure 11, we obtain Ud = 6 m/day with maximal rise rate of as much as 4 m/min during plastic failure. Sills intruded into stiffer rocks with ηr = 1018 Pa·s, a value typical for a thermal aureole of an intrusion, will undergo diapiric rise at a velocity of 2 cm/yr; i.e., a diapiric component can only contribute little to the commonly observed deformation rates of as much as several meters per year associated with sill intrusion. The main factor of uplift in this case is magma overpressure.

Cooling and Deformation Rates

Another important time scale involved is that of the cooling and solidification of magma bodies. Unlike salt diapirs, driven by a density instability without temperature contrast with country rocks, magma sills cool and solidify; therefore, in order to be effective, the time of the screen-rock instability and the ability to fail between magma sills by diapiric flow (tb = hb/Us, where hb is screen or bridge thickness) should be less than the characteristic conductive heat dissipation time t1: t1 = hsill2/kT, where hsill is half-thickness of the sill, and kT is thermal conductivity. A ratio of these two time scales gives: 
Note that the last multiplier in the above expression is a quasi-Rayleigh number, which defines the effectiveness of slow motion versus heat loss. Equation 6 demonstrates that in order for deformation to remain self-similar, a factor of 10 increase in the sill half-width W1/2 requires a viscosity increase by a factor of 1000. The condition of the buoyant rise time scale over the conductive cooling time scale can be thus expressed as this nondimensional quasi-Rayleigh number multiplied by the geometric factor. The geometric factor allows comparison of nondimensionalized parameters for different geometric configurations. 
The geometric factor ratio β is proportional to the product of the inverse aspect ratio (W1/2/hsill = 5) and the ratio of the sill half-width to the bridge width (hsill/hb = 1). The resulting value of β here is 0.067.

We can compute the required viscosities that are required for the gravity instability to be important on the cooling time scales (3–300 yr) inferred for the relatively small, 10–100-m-thick magma sills considered in this work. For Ra = 2 at β = 0.067 it yields small viscosities in the range of 1012–1015 Pa·s that are probably typical for crystal mushes (Lavalee et al., 2007) or for some nonwelded pyroclastic rocks and hyaloclastites. For rock viscosity 1018 Pa·s, a value typical for rocks near the brittle-ductile transition and in thermal aureoles of intrusions (e.g., Newman et al., 2006), and to satisfy Equations 6 and 7 for Ra = 2 and Δρ = 300 kg/m3, the sill sizes are W1/2 = 2165 m, and W1/2 = 3702 m, if the Ra = 10. Therefore, a significant increase in size is required to compensate for an increase in viscosity, and so only less viscous bridges are capable of deformation due to buoyancy for geologically realistic thicknesses and cooling time scales.

Initial smooth deformations are accelerated and modified after plastic failure of the surrounding rocks. For the geometry presented in Figure 10 we find that the threshold nondimensional yield strength, forumla providing for plastic failure of the bridge, rises linearly with Ra and is 18 at Ra = 80. The dimensional value of Sy can be recalculated via natural scale for pressure as: 
With Δρ in the range 300–500 kg/m3 and l0 in the range 40–100 m, we estimate Sy as 10–120 kPa. Such low values of Sy are compatible only with shear strength of the nonwelded pyroclastics at the effective mean stress approaching zero (σn, eff = 10 kPa, Sy = 30–350 kPa, Cecconi et al., 2010) or crystal mushes (Sy = 0.14 kPa for the solid fraction of isometric crystals, 50% and less; Hoover et al., 2001; or 0.3–4.3 kPa; Philpotts and Carroll, 1996).

Our analysis demonstrates that buoyancy does not play a significant role in the merging of small (10–102 m) magma batches, unless they are already located in a weak mush. The brute force buoyancy approach in this case requires that magma batches are already of magma reservoir size, >∼15km, in order to cause diapiric rise, but such sizes may cause caldera collapse (e.g., Ramberg, 1970; Geyer et al., 2006; de Silva et al., 2006; Gregg et al., 2012).

Overpressure Versus Buoyancy—What Controls the Merging Processes and Stability of Multibatch Magmatic Systems?

Figures 10 through 12 demonstrate that overpressure due to magma production and supply (which causes sill formation on the first place) is primarily responsible for the merging of magma sills. This is achieved through mechanical failure and dike focusing across separating screen rocks connecting dikes, followed by the initiation of magma coalescence. Density difference between rock and magma may cause rapid coalescence of magma lenses only if they are very close, longer lasting, and larger than the country rock screens, which have low viscosities and attributes of partially molten rocks, crystal mushes, or nonwelded pyroclastic rocks. Buoyant merging may be important at later stages of magma coalescence on a caldera scale, and may be one of the causes of some caldera-forming eruptions (Malfait et al., 2014; Caricchi et al., 2014), although the mechanisms for achieving the excess buoyancy postulated remain unclear. Overpressure or external tectonic controls (Allan et al., 2012) were proposed as an explanation for the Taupo volcanic zone (New Zealand) eruptions.

The analysis here demonstrates that widely separated magma batches may exist without interaction. In 2-D, stresses around circular cavities decay as 1/(R/Rc)2, where Rc is the cavity radius, in agreement with theoretical and numerical results (e.g., Jellinek and DePaolo, 2003; Simakin and Ghassemi, 2007). The stresses are thus reduced ∼9 times when the distance between batches is increased from 1 to 4 times their radius. This has implications for the long-term stability of widely dispersed magma batches with distinct crystal cargoes. The magma batches may exist without merging, crystallizing new zircons with rejuvenated U-Pb ages and batch-averaged εHf and δ18O (Fig. 4).

Overpressure in the Interstitial Fluid Phase

Interstitial hydrothermal fluid in country rocks may play an important role in low-pressure upper crustal conditions, in porous rocks, and especially in buried tuffs. Fluid pressure may rise rapidly in the thermal aureoles around magma lenses upon their emplacement or redistribution, leading to mechanical failure of the screen rocks in the fast elastic mode. Heating of the pore fluid increases its pressure almost proportionally to the temperature (e.g., according to gas law) at a constant volume of pores, or faster if there is a phase change. Real rock expands and fluid migrates in response to the internal pressurization, and the net effect is determined by the rock poroelastic properties and permeability. The direct proportionality of fluid pressure with respect to temperature and its dependence on other parameters at the impermeable boundary of the heated porous half-space was provided by McTigue (1986). After simplifications, 
where G is the shear modulus, ν and νu are the drained and undrained Poisson’s ratio, respectively, jo is the porosity, c and kT are the pressure and temperature conductivities, aT,fl is the fluid expansion coefficient, and ΔT is the temperature increment. Inserting typical parameter values from Detournay and Cheng (1993) and calculating aT,fl for an ideal gas at T = 800 °C, we get a fluid pressure in the range of 16–470 MPa (e.g., typical for sandstone and marble). The fluid pressure increase is the greatest for rocks with low permeability and high porosity, such as impermeable pumice, poorly consolidated, clay-rich sediments, or hydrothermally altered rocks. Intrusion of sills in such sediments at depths of <2 km results in fluidization of the rocks (Reid, 2004; Buttinelli et al., 2011; Schofield et al., 2010), sometimes giving rise to clastic dikes (Svensen et al., 2006). Even at greater depths of 6–8 km, common for magma reservoirs, pressurization effects might reach tens of megapascals; the details depend on the poroelastic properties and permeability of the country rocks. We stress that the rocks with large isolated pores, characteristic of vesiculated pyroclastic materials, are candidates for failure due to the strong pressurization and mechanical weakening at the thermal aureoles of sills, which we model. This makes realistic scenarios of sills merging due to buoyancy (see Figs. 11 and 12), as rocks with low viscosity and Sy can be interpreted as requiring fluid overpressure-degraded cataclastic materials.

Mixing Inside and Between Reservoirs

We next consider mixing effectiveness within merged batches. Many papers in the past 50 yr have documented abundant natural examples of efficient magma mixing and mingling on short time scales of magma chamber refill, or on mineral dissolution time scales for a variety of compositions. Examples including mixed and unmixed ignimbrites are also abundant (for recent numerical treatment of this problem, see Huber et al., 2009, 2012; for a review of diffusion time scales, see Cooper and Kent, 2014). For example, buoyant overturn caused by bubbles demonstrates the efficiency of this process on the order of hours or days (Ruprecht et al., 2008), and work in the Taupo volcanic zone demonstrates effective rhyolite mixing on eruptive time scales (Wilson and Charlier, 2009; Allan et al., 2013).

Mixing Within a Single Magma Batch

In these sets of numerical experiments we studied the dynamics of the magma reservoir with fine spatial resolution, down to several centimeters, sufficient to resolve fine flow structures in the compositional field and thus evaluate goodness of mixing (Figs. 13 and 14; see Supplemental File 77). This effort continues from Simakin and Bindeman (2012), where coupled Navier-Stokes and heat transfer equations were solved with realistic phase diagrams and proper viscosity models. We observed in these and previous melting experiments that convective melting induces vigorous rates of compositional convection, significantly exceeding the intensity of thermal convection in the solidifying magma chamber. This is essentially due to the high efficiency of in situ melting in generating density instabilities, as compared with the far-field control of chamber cooling. Although somewhat counterintuitive, it is thus the melting process that leads to convection and homogenization inside of magma bodies (Fig. 14; see Supplemental Files 1, 2, and 3 [see footnotes 1, 2, and 3]). We observed that any roof heterogeneity in density or compositional character (e.g., higher melting temperature blocks) will release gravitational energy and induce mixing upon downward movement through the magma body.

Here we pay greater attention to the motion and mixing in the horizontal direction responsible for the homogenization in sills with large aspect ratios (D/H), relevant for magma bodies parental to large-scale ignimbrite eruptions. Thermocompositional plumes chaotically descending from the roof cause rotations of the multiphase fluid, and are characterized by a spotty vorticity distribution. We observe that due to relatively small scale stochastic flow, liquid particles caught by successive vortexes rotate in different directions. Their horizontal motion and horizontal displacement can be approximately represented by random jumps with steps proportional to a vortex diameter in size. An example of the distribution of the absolute difference between initial and final horizontal coordinates, i.e., 
of the magma markers inside a model sill of 7 × 24 m is presented in Figure 13, recorded in 5 and 17 days from the beginning of the numerical experiment. Absolute displacements are averaged over a subgrid with size equal to one-half of the velocity elements (one-quarter by area).

Trajectories of 6 markers initially placed in the filled square are mapped after 8 days in Figure 15A. Markers first move coherently and then spread across the studied magma domain. In Figure 15B, absolute horizontal displacement δXm is plotted as a function of time for the averaged trajectories of 25 markers, originally located in the same square. During the first 2 days, this parameter oscillates around 1 m, which corresponds to the lifetime of the individual plume (vortex), where particles move vertically with small periodic horizontal displacements. An abrupt rise of δXm corresponds to the moment of pulling apart the magma parcel. Later, the δXm increases, yielding oscillations with periods of 2.5–2 days and average rates of ∼0.5 m/day.

Dependence δXm on time (t) after 10 days is nearly linear with a rate of 0.5 m/day (Fig. 15B; R2 = 0.969). However, square root dependence on time is expected on the basis of physical 1D Markov chains model with diffusive stochastic jumps (Levin et al., 2008). If so, the diffusional relationship allows the definition of the marker diffusion coefficient: 
where δx is the size of the jump each δt time units, and D is the effective diffusion coefficient. By using δXm (t) data for time >4 days we can derive the square root of time approximation with diffusion coefficient, D = 4 m2/day. These are obviously very fast rates as compared to the chemical diffusion of SiO2 or ZrO2 in rhyolitic magma, 10–9 cm2/day (Baker, 1991; Harrison and Watson, 1983).
Another way to treat mixing is to apply average separation between the markers approach (e.g., Ruprecht et al., 2008), as implemented for magma chambers and mantle mixing problems. Average separation is defined as the average distance between initially neighboring pairs of material points in some domain Ω with NΩ points: 
This parameter δR as a function of time was calculated for 16,000 markers initially located inside the 1 × 1m square shown in Figure 15A, and results of the numerical run are plotted in Figure 15C. Until ∼4 days, separation almost does not change; in the time interval 4–11 days it undergoes exponential stretch, and demonstrates square root dependence on time afterward. Ruprecht et al. (2008) treated separation as an exponential function of time, using a Lyapunov parameter formulation with exponent defined as: 

This temporal variation of the exponent parameter in our experiments is displayed in the Figure 15C inset. The Lyapunov exponent attains a maximum plateau with a level of 0.6 day–1 that characterizes the intensity of convective mixing in our case. Obviously λ decreases with time, since convective velocities are finite while constant λ implies growth of average separation with an exponentially increasing rate. We observe in our runs that exponential mixing with the Lyapunov exponent is only a transient phenomenon.

The random walk model with either square root time or linear dependence allows for extension of the results obtained in the restricted space domain toward the geologically more relevant larger dimensions and larger magma bodies. In particular, retaining the same dynamic parameters with a simple increase in the size of magma chamber results in a linear regression prediction of δX = 1 km in 2 yr. The diffusional model predicts 1 km dispersion in 684 yr. However, as dispersion (mixing) occurs in parallel with melting (Figs. 8, 13, and 14), horizontal parameters are also affected by the overall rates of displacement of the convective body of the magma. During 2 yr of roof melting, the melting front will displace the solid-melt boundary by 6–20 m, given modeled melting rates, um, of 4–12 m/yr; during 684 yr of melting, the magma body may migrate 300–1000 m, comparable with overall vertical chamber size, further contributing to homogenization.

It is likely that this model of displacement due to repeated rotations in random directions in descending plumes provides a minimum proxy for chemical mixing and horizontal homogenization in magma bodies. This is because it neglects the contribution of (1) the curved trajectories of the descending interacting plumes, and (2) occasional larger scale flows (magmatic avalanches), as seen in some runs (See Supplemental Files 1, 2, and 3 [see footnotes 1, 2, and 3]) and as widely speculated to occur in nature (e.g., layered intrusions; Wager and Brown, 1967). In the latter case, horizontal flows are occurring with significant horizontal velocities of meters per hour, much faster than either the linear or ∼t1/2 horizontal spreading rates mentioned here. However, compositional homogenization on a centimeter scale is ensured by the slower and smaller scale motions. Efficiency of mixing is also dependent on the amplitude and spatial distribution of tracer heterogeneity (e.g., Figs. 14, 16, and 17), whether large scale or small scale. The hand-specimen isotopic homogeneity is achieved when isotopically distinct individual minerals such as zircons with diverse εHf and δ18O (Figs. 4 and 14) are brought together from areas that were originally tens of meters apart.

In summary, our numerical results demonstrate that thorough compositional mixing is easy, and proceeds on time scales faster than cooling, but slower than mineral diffusive time scales. The process generates homogeneous whole-rock composition with isotopically diverse crystalline cargo, including zircons.

Comparison with Natural Data on Mixing

There are attempts to study convective homogenization in a magma reservoir in experiments with real silicate melts. De Campos et al. (2011) experimentally performed mixing of siliceous and basaltic magmas in a high-T cell. The experimental setup with independently rotated noncoaxial cylinders generated a chaotic (nonperiodic) flow field (Fig. 17). Mingling in the experiments at the large viscosity contrast of a factor of ∼1000 proceeds through formation of the continuously twisted magma layers rather than separated drops. The textures produced are similar to the commonly observed interlayering of dark and light bands of glasses in obsidians and pumices. Geologic observations (e.g., Slabunov, 1995) show us that partially melted xenoliths of amphibole-bearing schist in an andesitic melt are characterized by elliptic shapes in section and a unimodal distribution of the their short axes with median values of 1.5–2 cm, a mean of ∼5 cm, and a maximum of 18 cm. In our numerical runs we observe comparable size distributions of the partially molten roof material scattered in the lower convecting rhyolite (see Fig. 16, inset). This means that our experiments correctly reproduce the initial stage of homogenization, continuing with the fast diffusive dissolution of the smallest drops. Our results contradict earlier published conclusions of the small efficiency of the convection in a magma reservoir and its inability to achieve real mixing (Gutierrez and Parada, 2010). Both natural data (Vazques et al., 2009; Girard and Stix, 2010) and numerical models (Huber et al., 2009; Burgisser and Bergantz, 2011) demonstrate high efficiency of mixing, extraction, and differentiation in crystallizing rhyolites, especially when relative crystal-melt separation promotes local scale mixing.

We thus consider that the large-scale merging of magma batches that we modeled here creates many possibilities for compositional convection and mixing. A three-layer melting experiment in Figure 8 further illustrates how compositional heterogeneities are efficiently mixed when magma batches are merged together, in this case by melting through the screen rock between.


Case for Short-Lived Batch Magma Bodies in the Stress Field of the Upper Crust

Local high-resolution seismic tomography (e.g., Fig. 18) provides information about the depths, sizes, and approximate shapes of possible magma bodies under Yellowstone and elsewhere (Beachly et al., 2012). Miller and Smith (1999) and Husen et al. (2004) imaged two separate regions of low Vp under the Sour Creek (Yellowstone) resurgent domes at the same depth of 7 km (Fig. 18). The isotopically diverse crystal cargoes in eruptive products best argues that that the supposed magma chambers in the SRP and elsewhere are short lived and ephemeral, like those shown in Figure 18, rather than a single convectively homogenized batholith. Provided enough basalt and heat flux is maintained, dispersed magma bodies can be stable in this form, crystallize, and be hydrothermally altered, or rapidly reactivated. Thermomechanical stability analysis of widely dispersed melt pockets is outside the scope of this paper, but instability can be triggered by a greater supply of basalt (thermal reactivation), or by changing stress conditions (tectonic reactivation). Our modeling demonstrates that formation and mixing of large, eruptable magma volumes of 102–103 km3 may proceed rapidly by either amalgamation of these numerous bodies or continuous growth of one or several large ones, consuming the others.

Accounting for the Amount of Silicic Magma Produced

The estimates of the volcanic/plutonic ratio for the erupted rhyolites vary widely in any volcanic area, from ∼10% erupted, 90% remaining (Smith, 1979), to 1:2 (Bonnichsen et al., 2008), or even 1:1, as we favor here due to efficient near-eutectic melting and extraction. These estimates affect the corresponding amount of heat and basalt required to produce the observed amounts of silicic melt. Basalt differentiation or partial melting of the mafic crust (slow process; Fig. 6), or a combination of the two of these processes, is a strong function of the melting efficiency, i.e., how much heat is wasted for heating of the country rock prior to and during partial melting (e.g., Annen et al., 2006; Dufek and Bergantz, 2005). Also, many such models do not take into account temperature-dependent thermal diffusivities and/or conductivities, which, if incorporated, would reduce the heat contribution needed from basalts by a factor of ∼2 (Whittington et al., 2009). Hot zone models require significant amounts of mafic cumulate to be produced at the end (Fig. 19), often in excess of the observed thickness of the lower crust, requiring that these cumulates delaminated and sank into the mantle. Thus the room problem exists not only for silicic plutons in the crust, but also for parental cumulates and restites in the lower crust and upper mantle. In addition, the material versus heat contribution of basalt is one of the longest-discussed problems of igneous petrology (e.g., Hildreth and Moorbath, 1988; Grunder, 1995; and many others).

However, all of these mass, space, and heat constraints are significantly relaxed if silicic rocks are produced by remelting and recycling of already existing silicic magma bodies, which were formed by any of the above mechanisms. As we advocate here, such remelting is easy and fast (fast process; Fig. 5), and leads to the sequestration of silicic melts into the upper crust. For example, in the Andes, the U-Pb zircon evidence suggests that nonerupted parts of older eruptions in the same magmatic system or sequence are commonly recycled (e.g., Folkes et al., 2011), and a similar process must also characterize the multicyclic San Juan volcanic field of Colorado, USA (Lipman, 2007), although we note that these systems appear not to be fingerprinted by low δ18O values, indicating recycling of previously erupted and altered material (e.g., Folkes et al., 2013). Much recent work has been done on Taupo volcanic zone in New Zealand, also pointing to zircon recycling and short residence times of magma segregation (Charlier et al., 2008; Wilson and Charlier, 2009; Allan et al., 2013): oxygen isotopic evidence on recycling there is still to be investigated.

What Remains in the Crust After Recycling?

The volcanic-tectonic processes that we envision to characterize the evolution of an SRP-type environment leave behind basaltified, topographically reduced (Fig. 1), and bimodal terrain in the wake of the plume. Internally, it should lead to the complex intracrustal structures of dikes and sills, not common in typical clean granite plutonic examples such as those of the Sierra Nevada (California). However, plutons represent significantly prolonged histories after cessation of volcanism (Paterson and Vernon, 1995; Zak et al., 2007). In particular, large screens of intracaldera roof rocks may sink to the bottom of a lower magma body, as well as result in mafic magma recharge. Furthermore, plutons may undergo doming and upward motion on the order of 105–106 m.y. time scales, creating internal density redistribution (Ramberg, 1970; Petford, 2003; Schubert et al., 2013). Here we emphasize an important conclusion: observations from plutonic examples significantly postdate the stages of magmatic evolution documented by the volcanic evidence (or processes inferred for the volcanic magma chamber), and thus the observational records may not be compatible (e.g., Glazner et al., 2004; Coleman et al., 2004) (e.g., the geophysically imaged gabbro or gabbro-dioritic sill that remains in the crust after passage of the Yellowstone plume; Husen et al., 2004; DeNosaquo et al., 2009). Discordance between plutonic and volcanic records should be expected because plutonic rocks represent significant posteruptive thermal evolution. So, links between plutonic and volcanic records should not be discarded just because there is chronological discordance.

Predicting Yellowstone’s Future from Studying Its Past

The causes of the production and eruption of magmas in the SRP (Fig. 1), have direct implications for the 620 ka Yellowstone caldera, and in particular on whether it represents the last of a series of eruptions in a caldera cluster (e.g., Watts et al., 2012) and a crustal block that has been remelted and recycled from top to bottom, and thus exhausted its eruptive potential, or if generation and eruption of large volumes of silicic magma is still possible. In other words, it is important to estimate if the current (past 600 k.y.) voluminous rhyolitic eruptions at Yellowstone represent the path of a dying cycle, and if a new Huckleberry Ridge Tuff–size magma reservoir is currently being slowly established in southwest Montana with anticipated completion in 1–2 m.y. If so, the next eruption is to be normal δ18O, assembled slowly and with difficulty from the lower crust and differentiation of SRP basalt, with limited zircon diversity.

The low δ18O magmas with diverse zircons (Figs. 2–4) that appeared at the end of each caldera cluster evolution (and throughout the central SRP) for at least the past 12 m.y., and likely as early as 15 Ma (Colon et al., 2013) are identified as derived by recycling of buried silicic (zircon saturated) hydrothermally altered volcanic and plutonic predecessors from the same crustal block. Relative similarities in the radiogenic isotope ratios and trace elemental abundances in early normal δ18O and subsequent low δ18O rhyolites (Hildreth et al., 1991; Hanan et al., 2008; Christiansen and McCurry, 2008; Drew et al., 2013) suggest that the low δ18O rhyolites are derived by recycling and remelting of the normal δ18O rhyolites after they became depleted in 18O by high-temperature alteration, which did not disturb, or minimally affected elemental and other isotopic abundances. Given that Yellowstone has already erupted >600 km3 of recycled low δ18O rhyolites with diverse crystal cargo for the past 620 k.y., we consider that it is on a dying cycle, and that the crustal block under Yellowstone has been depleted in silicic melt, which was sequestered upward and recycled.

Big Picture View of Silicic Magma Petrogenetic Models and Crustal Evolution by Sequestration

The large-scale recycling model of rhyolite petrogenesis in the upper crust that we propose here (Fig. 19) serves as a previously underappreciated, important end-member to existing models of formation of silicic magmas in the hot zone in the lower crust during successive mafic sill intrusions (Hildreth and Moorbath, 1988; Huppert and Sparks, 1988; Annen and Sparks, 2002; Dufek and Bergantz, 2005; Annen et al., 2006), and the resulting assembly of large silicic magma bodies or batholiths in the middle and upper crust (Schaltegger et al., 2009). Extraction of the partial silicic melt from lower crustal mafic protoliths (slow and difficult process above and in Fig. 7) is facilitated by the slow coeval deformation (e.g., Spiegelman, 2003; Rabinowicz and Vigneresse, 2004; Simakin and Talbot, 2001). The batholiths may be initially mushy or fully crystalline (e.g., Coleman et al., 2004; Bachmann and Bergantz, 2004; Allan et al., 2013), or may not be batholiths at all, but an agglomeration of segregated silicic sills and stocks (Fig. 19; Horseman et al., 2009). Our model of upward sequestration (fast and easy process) then requires that these already segregated areas are remelted more than once, each time accompanied by further fractionation that increases crustal differentiation indices (e.g., increasing Rb/Sr, Eu anomalies, and SiO2, decreasing MgO). Each time, melting of siliceous lithologies proceeds rapidly; the 50% melt extraction window (e.g., Fig. 5) is achieved quickly as the heat of intruding hot silicic differentiates is very efficiently consumed by the eutectic melting, without much heat loss on rapid melting time scales. Basaltic magma may also intrude the upper crust along the magma pathways (Fig. 19; e.g., Fowler and Spera, 2010), delivering heat but offsetting the progress of sequestration. Rhyolitic bodies are then segregated in large volumes by upward sequestration, and our modeling demonstrates easy assembly, amalgamation, and effective mixing between disparate magma bodies.

We point to the observation that bimodal basaltic-rhyolitic provinces predominantly occur within extensional and hotspot environments with high basaltic power outputs, but also in continental arc areas with delamination features. Rising basaltic magmas are trapped by density filters of silicic magma batches at different levels in the crust, leading to a Daly Gap (i.e., a lack of intermediate compositions; Seligman et al., 2014). The described process (Fig. 19) furthermore leads to crustal sequestration of silicic material in the upper crust and mafic cumulates in the lower crust. We further speculate that because high basaltic output conditions were predominant on the early Earth, the recycling and sequestration processes were responsible for the segregation of the upper silicic and lower mafic crust, with significant components of thermomechanical lower crustal convection (e.g., Gerya, 2009). A dynamic view of silicic magma petrogenesis by sequestration is a must.

Supported by National Science Foundation EAR-CAREER-0844772 and EAR-1049351, and Russian Foundation for Basic Research grant 07-05-00629. We thank Matthew Loewen, Dana Drew, Dylan Colon, and Taras Gerya for comments, Chris Folkes, Graham Andrews, Ben Ellis, and two anonymous reviewers for their extensive comments, and Shan de Silva for helpful editorial handling.

1Supplemental File 1 is a generic melting by hot rhyolite intruded into cold, solid rhyolite, both with realistic phase diagrams (see Simakin and Bindeman, EPSL, 2012). If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES00969.S1 or the full-text article on www.gsapubs.org to view Supplemental File 1.
2Supplemental File 2 is a simulation of the bridge failure by “melting it through” without considering gravity or stress fields leading to subsequent mixing. The lower magma is a hot rhyolite with T=870°C, the upper is cold, near solidus rhyolite with T=730°C, the green bridge is solid, higher melting T. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES00969.S2 or the full-text article on www.gsapubs.org to view Supplemental File 2.
3Supplemental File 3 is a melting run emphasizing two refractory blocks in the roof rock, which have higher melting temperature (e.g. slivers of pre-caldera ancient rocks). They release potential energy upon fall, promote mixing, and carry down frozen-in rhyolitic material. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES00969.S3 or the full-text article on www.gsapubs.org to view Supplemental File 3.
4Supplemental File 4 shows batch merging in gravitational and stress fields in viscoplastic crust, as described in Figs. 9 and 10 in the text of the paper. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES00969.S4 or the full-text article on www.gsapubs.org to view Supplemental File 4.
5Supplemental File 5 shows batch merging in gravitational and stress fields in viscoplastic crust, as described in Figs. 9 and 10 in the text of the paper. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES00969.S5 or the full-text article on www.gsapubs.org to view Supplemental File 5.
6Supplemental File 6 shows batch merging in gravitational and stress fields in viscoplastic crust, as described in Figs. 9 and 10 in the text of the paper. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES00969.S6 or the full-text article on www.gsapubs.org to view Supplemental File 6.
7Supplemental File 7 is a temporal evolution of vorticity distribution, which is responsible for mixing. If you are viewing the PDF of this paper or reading it offline, please visit http://dx.doi.org/10.1130/GES00969.S7 or the full-text article on www.gsapubs.org to view Supplemental File 7.