Silicic magma reservoirs in the Earth’s crust

Determining the shape, size, depth, and state of magma bodies in the Earth’s crust as well as how they evolve over time remain key issues for several Earth science communities. With a better determination of these variables, petrologists could construct more accurate chemical models of differentiation processes; magma dynamicists could postulate causes for magma migration, storage, and interaction at different levels within the crust; volcanologists could better predict the style and volume of upcoming volcanic eruptions; and geochronologists could construct an evolutionary trend with greater precision. Although our knowledge on magmatic systems has expanded greatly over the past decades, much remains to be done to constrain the multiphase, multiscale processes that are at play in those reservoirs and the rates at which they occur.

Previous attempts to summarize the state of knowledge in this field (e.g., Smith 1979; Lipman 1984; Sparks et al. 1984; Marsh 1989a; de Silva 1991) have greatly helped crystallizing ideas on the contemporaneous state of knowledge and possible alleys for future directions to explore. Since then, many more publications have attempted to draw magma reservoirs geometries and internal complexities [see Fig. 1 for a potpourri of such schematic diagrams for several examples of large caldera systems in the U.S.A., as well as the more recent reviews by Petford et al. (2000), Bachmann et al. (2007b), Lipman (2007), Cashman and Sparks (2013), Cashman and Giordano (2014), and de Silva and Gregg (2014)]. At this stage, it is clear that magmas are complex mixtures of multiple phases (silicate melt, H₂O-CO₂-dominated fluid, and up to 10 different mineral phases in some systems) with very different physical properties. For example, viscosity contrast between water vapor and crystals can span more than 20 orders of magnitude, and even vary significantly within one given phase (e.g., several orders of magnitude for silicate melt), as composition- $P-T-fH2O-fO2$ vary (Mader et al. 2013). The way these magmas move and stall in the crust will therefore reflect the complex interaction between all these phases and the typically colder wall rocks.

Over the last decades, an impressive body of work has accumulated on the characterization of magma mainly through: (1) high-precision, high-resolution geochemical work [with ever more powerful instruments, including the atom-probe, e.g., Valley et al. (2014) and nanoSIMS in the recent past; e.g., Eiler et al. (2011); Till et al. (2015)], and (2) experimental petrology (e.g., phase diagrams, trace element partitioning between phases, isotopic fingerprinting), culminating in powerful codes for thermodynamic modeling (Ghiorso and Sack 1995a; Boudreau 1999; Connolly 2009; Gualda et al. 2012). This work will continue to greatly improve the constraints that we place on these magmatic systems for years to come, but we argue that the field will have to take into account kinetic, disequilibrium processes, which are commonplace, and can be severe in the world of magmas.

Thermodynamics and kinetics are both valuable to understand magmatic processes, and kinetics is bound to become even more important in the following decades as we strive to better constrain the timescale of dynamical processes associated with magma transport and evolution in the crust. This Centennial Volume of American Mineralogist appears as an appropriate landmark to summarize the recent advances on the subject and bring the multiple facets of this problem together to better frame the future directions of research.

The most exciting scientific challenges are, we believe, those that are driven by enigmas arising by observations of Nature. Below are a few examples of key questions or long-standing observations that relate to magma reservoirs and magmatic differentiation:

1. The fundamental dichotomy in supervolcanic deposits (also called ignimbrites or ash-flow tuffs), which show, in most cases: (a) compositionally and thermally zoned dominantly crystal-poor units [sometimes grading to crystal-rich facies at the end of the eruption (Lipman 1967; Hildreth 1979; Wolff and Storey 1984; Worner and Schmincke 1984; de Silva and Wolff 1995)], and (b) strikingly homogeneous crystal-rich dominantly dacitic units (the “Monotonous Intermediates” of Hildreth 1981), known to be some of the most viscous fluids on Earth (Whitney and Stormer 1985; Francis et al. 1989; de Silva et al. 1994; Scaillet et al. 1998b; Lindsay et al. 2001; Bachmann et al. 2002, 2007b; Huber et al. 2009; Folkes et al. 2011b; see also Bachmann and Bergantz 2008b; Huber et al. 2012a for reviews).

2. The striking resemblance in chemical composition (for major and trace elements) between the interstitial melt in crystal-rich silicic arc magmas and the melt-rich rhyolitic magmas (Bacon and Druitt 1988b; Hildreth and Fierstein 2000; Bachmann et al. 2005).

3. The ubiquitous observation of magma recharges from depth prior to eruptions, and their putative role on the remobilization of resident (often crystal-rich) magmas (“rejuvenation”), for large and small eruptions (e.g., Sparks et al. 1977; Pallister et al. 1992; Murphy et al. 2000; Bachmann et al. 2002; Bachmann and Bergantz 2003; Wark et al. 2007; Molloy et al. 2008; Bachmann 2010a; Cooper and Kent 2014; Klemetti and Clynne 2014; Wolff and Ramos 2014; Wolff et al. 2015).

4. The long-noticed compositional gap in volcanic series, coined the Bunsen-Daly Gap after R. Bunsen and R. Daly’s early observations (Bunsen 1851; Daly 1925, 1933), and confirmed in many areas worldwide since (Chayes 1963; Thompson 1972; Brophy 1991; Thompson et al. 2001; Deering et al. 2011a; Szymanowski et al. 2015) even where crystal fractionation, which should lead to a continuous melt composition at the surface, dominates (Bonnefoi et al. 1995; Geist et al. 1995; Peccerillo et al. 2003; Macdonald et al. 2008; Dufek and Bachmann 2010; Sliwinski et al. 2015).

5. The complex relationship between silicic plutonic and volcanic units with possible hypotheses (not mutually exclusive in plutonic complexes) ranging from plutons being “failed eruptions” (crystallized melts; Tappa et al. 2011; Barboni et al. 2015; Glazner et al. 2015; Keller et al. 2015) to plutons being “crystal graveyards” (crystal cumulates; Bachl et al. 2001; Deering and Bachmann 2010; Gelman et al. 2014; Putirka et al. 2014; Lee and Morton 2015). Related to this volcanic-plutonic connection is the corollary that granites sensu stricto (evolved, high-SiO₂, low-Sr composition) appear to be less common than their volcanic counterparts (e.g., Navdat database, Halliday et al. 1991; Hildreth 2007; Gelman et al. 2014).

6. The observation that silicic plutons, when fully crystallized, contain little volatile left (typically <1 wt% H₂O that is bound up in the minerals; e.g., Caricchi and Blundy 2015) despite the fact they started with volatile concentrations above 6 wt%, implying significant loss of volatile at different stages of the magma bodies’ evolution [e.g., by first and second boiling, pre-, syn-, and post-eruption degassing, see for example Anderson (1974); Wallace et al. (1995); Candela (1997); Lowenstern (2003); Webster (2004); Wallace (2005); Bachmann et al. (2010); Blundy et al. (2010); Sillitoe (2010); Baker and Alletti (2012); Heinrich and Candela (2012)].

7. The well-known “Excess S,” highlighted by the much higher S content released during eruptions (measured by remote sensing or estimated from ice-core records) that can be accounted for by the melt inclusion data [i.e., petrologic estimate, see, for example, reviews by Wallace (2001); Costa et al. (2003); Shinohara (2008)].

8. The remarkable cyclic activity observed large-scale volcanic systems, which culminates in caldera-forming eruptions (“caldera cycles”) that can repeat itself several times, e.g., the multi-cyclic caldera systems such as Yellowstone (Hildreth et al. 1991; Bindeman and Valley 2001; Christiansen 2001; Lowenstern and Hurvitz 2008), Southern Rocky Mountain Volcanic field (e.g., Lipman 1984, 2007; Lipman and Bachmann 2015), Taupo Volcanic Zone (e.g., Wilson 1993; Wilson et al. 1995; Sutton et al. 2000; Deering et al. 2011a; Barker et al. 2014), High Andes (e.g., Pitcher et al. 1985; Petford et al. 1996; Lindsay et al. 2001; Schmitt et al. 2003; de Silva et al. 2006; Klemetti and Grunder 2008), Campi Flegreii (e.g., Orsi et al. 1996; Civetta et al. 1997; Gebauer et al. 2014), and Aegean arc (e.g., Bachmann et al. 2012; Degruyter et al. 2015).

All those questions/observations require an overarching concept of magma chamber growth and evolution that ties together these seemingly disparate concepts. Before attempting such a task, first this briefly outlines the methods that are typically used to infer the state of these reservoirs, highlighting some of their strengths and limitations.

The foundation of volcanology is motivated by observations of ongoing and past volcanic activity, with the purpose of understanding the underlying processes that govern the evolution of magmatic systems on Earth and other bodies of the Solar System (Wilson and Head 1994). However, large-scale eruptions are rare (Simkin 1993; Mason et al. 2004) and as volcanologists/petrologists we rely much on a forensic approach: study magmatic rocks (both plutonic and volcanic) that formed/erupted in the past, and try to reconstruct the conditions (P, T, $fO2$⁠, $fH2O$⁠, $fSO2$⁠, …) that prevailed in the reservoirs at the time of formation or eruption [see, for example, a review of techniques by Putirka (2008) or Blundy and Cashman (2008)]. Of course, volcanic and plutonic lithologies do not provide the same type of information, and could/should be complementary. Volcanic rocks carry limited spatial context for the reservoir, but provide an instantaneous snapshot into the state of the magma body just before the eruption. In contrast, plutonic rocks present a time-integrated image of the magma accumulation zone, often with a history spanning multimillion years (e.g., Greene et al. 2006; Walker et al. 2007; Schoene et al. 2012; Coint et al. 2013), for which we can tease out some information about sizes and shapes of magmas bodies. Much of the geochemical data acquired over the last century is now tabulated in large databases such as Georoc (http://georoc.mpch-mainz.gwdg.de/georoc/), EarthChem (http://www.earthchem.org), and Navdat (http://www.navdat.org), providing a remarkable resource to analyze global geochemical problems (see recent papers of Keller and Schoene 2012; Chiaradia 2014; Gelman et al. 2014; Glazner et al. 2015; Keller et al. 2015).

Measuring gases coming out of active volcanoes provides us with some key information about the state of magmatic systems (e.g., Giggenbach 1996; Goff et al. 1998; Edmonds et al. 2003; Burton et al. 2007; Humphreys et al. 2009). Due to their low density and viscosity, volcanic gases can escape their magmatic traps and provide “a telegram from the Earth’s interior,” as laid out by one of the pioneers of gas measurements, Sadao Matsuo (e.g., Matsuo 1962). For example, the abundance of chemical species like He, CO₂, or Cl can help determine the type of magma that is degassing (Shimizu et al. 2005), and the potential volume that is trapped in the crust (Lowenstern and Hurvitz 2008). As mentioned above, the amount of S released during eruptions is also key in estimating the presence of an exsolved gas phase in the magma reservoir prior to an eruption (Scaillet et al. 1998a; Wallace 2001; Shinohara 2008). As discussed by Giggenbach (1996), the composition of volcanic gases measured at the surface integrates both deep, i.e., magma reservoir processes, and shallow processes within the structure of the edifice, and possibly the interaction with a hydrothermal system. As such it is often challenging to directly relate gas composition to the conditions of shallow magma storage (e.g., Burgisser and Scaillet 2007) and a better monitoring strategy is generally to look for relative changes in gas compositions and temperature rather than focusing on absolute values.

With the exception of gas measurements, imaging active magma reservoirs mainly involves geophysical methods. Several techniques are now available, and provide different pictures of active magmas bodies. Such techniques include: (1) seismic tomography [both active and passive sources (see for example Dawson et al. 1990; Lees 2007; Waite and Moran 2009; Zandomeneghi et al. 2009; Paulatto et al. 2010), and ambient-noise tomography (e.g., Brenguier et al. 2008; Fournier et al. 2010; Jay et al. 2012)]; (2) magneto-telluric surveys (e.g., Hill et al. 2009; Heise et al. 2010); (3) deformation studies using strain-meter, GPS, and satellite data (e.g., Massonnet et al. 1995; Hooper et al. 2004; Lagios et al. 2005; Newman et al. 2006; Hautmann et al. 2014); and (4) muon tomography on steep volcanic edifices (Tanaka et al. 2007, 2014; Gibert et al. 2010; Marteau et al. 2012).

Geophysical methods are based on the inversion of signals transmitted through the crust and measured from or close to the surface. The choice of methods and, in several cases, the frequency band of the signals studied depends on the targeted resolution and depth of interest. Geophysical methods probe variations in physical properties caused by the presence of partially molten reservoirs in the crust and provide two- or three-dimensional images of these systems. These inversions however do not provide an unequivocal picture of the thermodynamical state of the magma storage region, because several variables (melt fraction, temperature, exsolved gas content, composition, connectivity of the various phases) affect the elastic, magnetic, and electric properties of magmas. Similarly, surface deformation signals are not only function of the state, shape, orientation, and size of a magma body, but also depend on the crustal response to stresses around the body (through anelastic deformation and movement along weakness planes such as fractures and faults; Newman et al. 2006; Gregg et al. 2013).

In summary, geophysical surveys mainly provide constraints about the depth, shape, and extent of a partially molten (or hot) region under an active volcanic edifice. The many surveys using inSAR, GPS field campaignsm or ambient-noise tomography where, as for gas monitoring, the focus lies on relative changes rather than the effective state of the reservoir have significantly increased our ability to monitor volcanic unrest (see recent review by Acocella et al. 2015). They efficiently highlight rapid (decadal or less) changes in the mechanical state of magmatic systems, which is generally difficult to achieve using petrological data sets. However, the spatial scales probed by geophysical imaging techniques are limited by the resolution of the method considered. Under most circumstances, the spatial resolution remains greater or equivalent to the kilometer scale (Waite and Moran 2009; Farrell et al. 2014; Huang et al. 2015), which is in stark contrast with the extremely high resolution of geochemical and petrological data. As such, geophysical inversions can be thought of as upscaling filters of the state of heterogeneous multiphase reservoirs. A fundamental assumption underlying these inversions is that spatial resolution is commensurable with a Representative Elementary Volume (REV) where the physical properties of the heterogeneous medium can be (linearly) averaged out into a single effective value (a more detailed discussion of the limitations imposed by this assumption is provided below).

Starting with N.L. Bowen at the turn of the 20th century, many geoscientists have applied experimental studies to study the thermodynamics that control the assemblage of phases in magmas with different bulk compositions. For nearly a century, a large amount of data has been collected on crystal-melt stability and composition as a function of several variables (P, T, fugacities of volatile elements and oxygen, compositions of magmas, cooling rate…). Much of these data are now available in databases, some in web-based portals, as is the case for the geochemical data (see for example Berman 1991; Holland and Powell 2011 or LEPR, http://lepr.ofm-research.org/YUI/access_user/login.php). The data sets range from crystal-melt equilibria in the mantle melting zones (e.g., see Poli and Schmidt 1995 to lower crustal (e.g., Muntener and Ulmer 2006; Alonso-Perez et al. 2009) and upper crustal conditions in many different settings (e.g., Johnson and Rutherford 1989; Johannes and Holtz 1996; Scaillet and Evans 1999; Costa et al. 2004; Almeev et al. 2012; Gardner et al. 2014; Caricchi and Blundy 2015).

For shallow magma reservoirs in arcs, a third phase, volatile bubbles, is also commonly present. This exsolved volatile phase does not only affect the thermo-physical properties of magmas (see below) but can also impact greatly the partitioning of some volatile trace elements and therefore influence the chemical fractionation that takes place in magma reservoirs. Laboratory experiments have played a major role in determining the solubility of volatile phases (H₂O, CO₂, S, CL, F…) as a function of pressure, temperature, and melt composition [see the review by Baker and Alletti (2012) and references therein]. Several species, such as Cl and S compounds, have a strong influence on the fractionation and transport of precious metals out of magma reservoirs to the site of ore deposits and have therefore received a thorough attention (Scaillet et al. 1998a; Williams-Jones and Heinrich 2005; Zajacz et al. 2008; Sillitoe 2010; Fiege et al. 2014).

Chemical equilibrium between the melt and newly formed crystals or crystal rims is often a decent assumption, but this assumption breaks loose when changes taking place in the reservoir are more rapid than the equilibration timescale (see Pichavant et al. 2007). Such rapid changes in reservoirs’ conditions can include: (1) magma recharges and partial (incomplete) mixing; (2) efficient physical separation of mineral and/or bubble-melt; and (3) during rapid pressure drops associated with mass withdrawal events (eruptions, dike propagation out of the chamber). In these particular cases, kinetics controls the mass balance between the different phases; one therefore needs to consider diffusion of elements in silicate melts and minerals, as well as mineral and bubble growth or dissolution. The diffusive transport of chemical species in silicate melts is a strong function of: (1) viscosity (which is a proxy for the degree of polymerization of the melt); (2) dissolved water content; and (3) to a lesser extent, oxygen fugacity (redox conditions in the magma affect speciation). We bring the attention of the interested readers to the review of Zhang and Cherniak (2010) for additional information on this topic.

The presence of kinetics is readily observed in mineral zoning patterns because solid-state diffusion is generally orders of magnitude slower than diffusion in silicate melts. Over the last decade, several models have been developed to use the diffusion of cations in silicate minerals and retrieve timescales relevant to the dynamic evolution of magmas before and during eruptions. These models rest heavily on experiments where the kinetics of diffusion as a function of temperature and composition is measured [see review of Zhang and Cherniak (2010) and reference therein]. Diffusion modeling in minerals informs us on the interval of time between a significant change in thermodynamic conditions of the magma (mixing from recharges for example) and the time at which the minerals transition across the closure temperature (e.g., eruption) where subsequent diffusion can be ignored. To date, studies have explored the diffusion of a growing quantity of cations in hosts such as quartz, plagioclase, pyroxene, biotite, and olivines (for more mafic systems than the ones considered here; see Fig. 2 for published examples).

Stable isotopes can also provide information regarding the thermodynamic conditions (equilibrium fractionation between different phases) and kinetics (kinetic fractionation) in the magma. Because of improved accuracy and spatial resolution of analytical measurements, there is a growing interest in using different stable isotope systems, such as O, S, Ca, as well as other major elements, to fingerprint the oxidation state of magmas (Metrich and Mandeville 2010; Fiege et al. 2015), crustal assimilation (e.g., Friedman et al. 1974; Taylor 1980; Halliday et al. 1984; Eiler et al. 2000; Bindeman and Valley 2001; Boroughs et al. 2012), or mixing processes between magmas with different compositions (Richter et al. 2003; Watkins et al. 2009a).

The migration, convection, and deformation of magmatic mixtures in the crust have direct implications on the rate and efficiency of fractionation processes, on the cooling rate of magma reservoirs, and eruption dynamics (e.g., Gonnermann and Manga 2007; Lavallée et al. 2007; Cordonnier et al. 2012; Gonnermann and Manga 2013; Parmigiani et al. 2014; Pistone et al. 2015). The response of magmas to differential stresses is complex, owing to the multiphase nature of these systems. Several research groups have performed deformation experiments to better constrain rheological laws pertinent to magmas (Spera et al. 1988; Lejeune and Richet 1995; Caricchi et al. 2007; Ishibashi and Sato 2007; Champallier et al. 2008; Cimarelli et al. 2011; Mueller et al. 2011; Vona et al. 2011; Pistone et al. 2012; Del Gaudio et al. 2013; Laumonier et al. 2014; Moitra and Gonnermann 2015). These studies have focused on diverse aspects of the rich non-linear dynamics of the rheology of suspensions, such as: (1) the effect of the crystal volume fraction; (2) the effective shear rate; and (3) the effect of polydisperse crystal size distributions and various crystal shapes. These results clearly show that suspensions become very stiff as crystal content approach or exceed 40 to 50% of the volume of the magma, although they provide little information about the onset and magnitude of a yield stress for crystal-rich suspensions. Even more challenging is the development of physically consistent rheology laws for three-phase magmas because of the additional complexity of the capillary stresses (coupling bubbles, melt, and crystals), bubble deformation, three phases lubrication effects and possible coalescence (Pistone et al. 2012; Truby et al. 2014).

Magma rheology experiments have been synthesized to yield different empirical (Lavallée et al. 2007; Costa et al. 2009) and semi-empirical (Petford 2009; Mader et al. 2013; Faroughi and Huber 2015) models for the deformation of magmas under shear flow conditions. The mechanical behavior of magmas remains a rich field for future investigations, and several important challenges remain to constrain the rate at which magma flows in reservoirs and up to the surface. They include:

• How to treat multiphase rheology when phase separation takes place concurrently?

• How do bubbles interact with crystals suspended in viscous melts, and how does this behavior depend on the strain rate and volume fractions of each phases?

• Can we constrain the yield strength of crystal-rich magmas under magma chamber conditions (imposed stresses rather than imposed strain-rate)?

• By rheology, we often restrict ourselves to measure the shear viscosity using experiments or models that relate stress and deformation under very simple geometries. As viscosity is a tensor, is it possible to measure its anisotropy for non-linear multiphase systems with complex crystal shapes under various strain rates and conditions pertinent to magma chambers?

The rheology of a magma body should therefore be considered as a dynamical quantity, which can vary significantly (orders of magnitude) because of changes in stress distributions, variations in crystal-melt-bubble phase fractions and potentially many other factors.

Several other physical properties have a resounding impact on the evolution of magmas in the crust. Magmas are far from thermal equilibrium in the mid to upper crust, and heat transfer in and out of these bodies controls their chemical and dynamic evolution to a great extent (Bowen 1928; McBirney et al. 1985; Nilson et al. 1985; McBirney 1995; Reiners et al. 1995; Thompson et al. 2002; Spera and Bohrson 2004). In that context, the heat capacity, latent heat, and thermal conductivity of silicate melts/minerals and CO₂-H₂O fluids are important to establish the timescale over which these reservoirs crystallize in the crust and by extension the vigor of convective motion in these systems. It is important to emphasize that thermal convection in the classical sense is not directly applicable to magma chambers under most conditions, because density contrasts between the different phases far exceeds those by thermal expansion (e.g., Marsh and Maxey 1985; Bergantz and Ni 1999; Dufek and Bachmann 2010). Additionally, the rate of crystallization directly impacts the cooling rate of magmas (due to non-linear release of latent heat during crystallization, e.g., Huber et al. 2009; Morse 2011), and this latent heat buffering becomes dominant as silicic magmas approach eutectic behavior near their solidus (Gelman et al. 2013b; Caricchi and Blundy 2015). Some recent studies have tested how the temperature dependence of thermal properties influence the cooling and crystallization timescales of magmas in the crust and showed that these effects are sometimes non-trivial (Whittington et al. 2009; Gelman et al. 2013b; de Silva and Gregg 2014).

An elegant approach to study the complex multiphase fluid dynamics that takes place in magma reservoir is to resort to laboratory tank experiments (see for example, McBirney et al. 1985; Nilson et al. 1985; Jaupart and Brandeis 1986; Martin and Nokes 1988; Shibano et al. 2012). Experiments provide the means to investigate some aspects of the complex non-linear dynamics that prevail in these systems, however, in most cases, it is difficult or impossible to satisfy dynamical similitude with real magmatic systems. Nevertheless, experiments reveal interesting feedbacks that may have been overlooked and can serve as a qualitative assessment of how processes couple to each other. One of the main directions for laboratory fluid dynamics experiments is to study how magmatic suspensions behave at low Reynolds number during convection, phase segregation, and magma recharges with mixing (e.g., McBirney 1980; Sparks et al. 1984; McBirney et al. 1985; Davaille and Jaupart 1993; Koyaguchi et al. 1993). Here, the term suspension is used loosely to describe both suspension with solid particles (crystals) and bubbly emulsions.

As discussed throughout this review [and clearly alluded to by early petrologists such as R. Daly as far back as the 1910s, e.g., Daly (1914)], magma bodies are open reservoirs that exchange mass and energy with their surroundings. Chemical heterogeneities and zonations then reflect either phase separation by gravity (e.g., McBirney 1980, 1993; Hildreth 1981; Sparks et al. 1984) or incomplete mixing between magmas with different compositions (Eichelberger 1975; Dungan et al. 1978; Blake and Ivey 1986; Eichelberger et al. 2000). Because buoyancy forces related to compositional differences (including presence of various phases of variable densities) exceed thermal buoyancy effect by several orders of magnitude, mechanical convection and stirring are dominated by chemical variations, bubble rise, and/or crystal settling. Strong spatial contrasts in bulk viscosity because of mixing between magmas with different compositions, contrasts in crystal content or thermal gradients also tend to stabilize against hybridization or magma mixing (Sparks and Marshall 1986; Turner and Campbell 1986; Koyaguchi and Blake 1989; Davaille and Jaupart 1993; Jellinek and Kerr 1999). Mixing in magma reservoirs can also be accelerated during eruptions, where substantial stress differentials can be generated by decompression or roof collapse (Blake and Campbell 1986; Blake and Ivey 1986; Kennedy et al. 2008).

Convection strongly influences phase separation, which in turn controls chemical differentiation. Several groups have investigated the effect of crystals (Jaupart et al. 1984; Brandeis and Jaupart 1986; Brandeis and Marsh 1990; Koyaguchi et al. 1993; Shibano et al. 2012), bubbles (Huppert et al. 1982; Thomas et al. 1993; Cardoso and Woods 1999; Phillips and Woods 2002; Namiki et al. 2003), and shapes of reservoirs (e.g., de Silva and Wolff 1995) on convective motion. The presence of crystals or bubbles controls the gravitational energy potential that is converted to kinetic energy during convection, as experiments confirm that they impact convection even at relatively small volume fractions. Obviously, suspended phases affect convection, but convection itself also influences crystal settling or bubble migration. For example, Martin and Nokes (1988, 1989) used a series of experiments to quantify the timescale over which a convecting suspension of crystals clears out because of Stokes settling. Hydrodynamic interactions between the different phases can significantly reduce the settling velocity and affect phase separation even at low-particle (bubbles or crystals) volume fractions (Faroughi and Huber 2015, see Fig. 3).

Numerical studies enable us to isolate specific dynamical feedbacks and to provide a framework to understand the complex temporal evolution of these systems, which can often be challenging to infer from natural samples or laboratory experiments where dynamic similarity is difficult to satisfy. Numerical approaches based on different sets of starting assumptions and governing equations have been used to address questions such as:

• The chemical evolution of open-system magma reservoirs undergoing assimilation, crystal fractionation and magma evacuation.

• The thermal evolution and longevity of magma reservoirs in the crust.

• The relative importance of crustal melting vs. crystal fractionation of more mafic parents in driving the evolution of magmas toward silicic compositions.

• The pressure variations in and around magma chambers and stability of reservoirs following recharges.

• Gas exsolution and how it impacts the dynamics in shallow magma reservoirs.

The first category of studies summarized below is based on box models where the reservoirs are considered homogeneous and internal structures are not explicitly accounted for. These models involve ordinary differential equations that depend solely on time (no spatial dependency) and by extension assume that homogenization takes place over a shorter timescale than the processes considered. In other words, in these models, the different phases (melt, crystals, sometimes bubbles) are assumed to remain in thermal and chemical equilibrium with each other over the scale of the whole reservoir.

Since the seminal work of H. Taylor and D. DePaolo (Taylor 1980; DePaolo 1981), where a box model was constructed to infer the relative importance of assimilation and crystal fractionation processes on the trace elements and isotopic record, many studies have expanded on the method to relate chemical tracers to the actual dynamical processes that control the evolution of magmatic systems. The scope of these models ranges from non-linear reactor models (Bonnefoi et al. 1995) investigating the development of bimodal compositions from a single magmatic system to constraining the effect of melt extraction on the chemical signature of the residual cumulate and high-silica melts (Gelman et al. 2014; Lee and Morton 2015) as well as focusing on the development of crystal zonations in a well-mixed reservoir (Wallace and Bergantz 2002; Nishimura 2009). These models provide insights into the chemical response of reservoirs to various fractionation processes, but do not explicitly treat the source of these processes (e.g., cooling and crystallization leading to gravitational settling or the reheating of the surrounding crust and its partial melting/assimilation).

As a first step toward answering these limitations, W. Bohrson and F. Spera developed, in a series of papers (Spera and Bohrson 2001, 2004; Bohrson and Spera 2007; Bohrson et al. 2014), a model that couples the mass balance of trace elements and isotopes in a reservoir where cooling and assimilation processes are parameterized. They also added a thermal energy equation for the reservoir. The energy balance in open magma chambers is quite complex and involves more than just sensible and latent heat contributions in response to the heat loss to the surrounding crust. Even within the assumption that the reservoir remains homogeneous and the phases remain in equilibrium at all times, enthalpy is constantly transferred from mechanical work to heat and vice versa. As an example, the injection of fresh magma into a reservoir affects the energy budget not only by providing a possible heat source but also exerts mechanical work (pressure changes) that affects the stability of the phases present in the magma. These feedbacks become even more important when shallow volatile-saturated systems are considered (Huppert and Woods 2002; Degruyter and Huber 2014; Degruyter et al. 2015). It is therefore important to extend the energy conservation statement in these models to include mechanical work and phase changes when possible [hence introducing a consistent two-way coupling between the hot magma reservoir and the surrounding crust (de Silva and Gregg 2014; Degruyter and Huber 2014)].

The assumption of thermal and chemical equilibrium at the scale of the reservoir is useful and allows us to draw interesting and informative conclusions on the general trends that govern magma chamber evolution. However, this assumption is untenable under most conditions. Introducing spatial heterogeneities in the magma body and therefore allowing for some degree of disequilibrium between the phases (at least down to the scale of the model’s resolution, where local equilibrium is generally still assumed between phases) is necessary to get realistic estimates of the temporal scales over which magma bodies evolve. More specifically, a spatiotemporal description of magma reservoirs provides the opportunity to test for the conditions under which efficient mixing is no longer possible (Huber et al. 2009) and chemical differentiation by mechanical separation (flow) occurs. We will divide the spectrum of numerical models geared to address these questions into: (1) thermal models on one hand and (2) coupled thermal and mechanical models on the other hand.

Thermal models are concerned only with one aspect of the energy balance and consider only sensible and latent heat. These models involve various levels of complexity, from equilibrium crystallization, i.e., the amount of crystallization over each time step is that predicted from equilibrium thermodynamics (e.g., Annen et al. 2008; Leeman et al. 2008; Annen 2009; Gelman et al. 2013b; Karakas and Dufek 2015) to the model of C. Michaut and coworkers where crystallization is considered kinetically limited (Michaut and Jaupart 2006). Equilibrium crystallization models rely on a constitutive equation that relates temperature and crystallinity for a given choice of magma composition and pressure, generally constrained experimentally or using thermodynamic models such as MELTS (Ghiorso and Sack 1995b). Because no momentum conservation is considered, crystals and melt are static and bound to each other, i.e., no phase separation and by extension no chemical differentiation by crystal fractionation is possible. Heat is therefore transported only by diffusion within and out of the magma reservoir. These models have been used extensively to study the longevity of active magma bodies in the crust as a function of the rate of magma injections. While they suggest high average rates of magma injections to counteract the heat lost to the crust, they overlook the role of mechanical processes and volatile exsolution on the energy balance in response to repeated magma recharges and sporadic evacuation events (Degruyter et al. 2015).

Another type of model accounts for the coupled momentum and enthalpy balance in magma reservoirs and has been used to model the thermal evolution as well as chemical differentiation (by crystal fractionation and assimilation) in magmatic systems. In these models, the momentum balance is either parameterized (Huber et al. 2009) or requires a set of coupled continuum equations that allows the different phases to separate from one another (Dufek and Bachmann 2010; Gutierrez and Parada 2010; Lohmar et al. 2012; Simakin and Bindeman 2012; Solano et al. 2012). In those latter models, the different phases are coupled in the momentum equation through drag terms. The first numerical models applied to mixing of chemical heterogeneities (single-phase fluids) in magmas by convective stirring where developed by C. Oldenburg and colleagues (Oldenburg et al. 1989; Spera et al. 1995), highlighted the different rates of mixing by normal and shear strains. Later, Huber et al. (2009) discussed a method to characterize the efficiency of mixing by convective stirring in time-dependent systems where cooling and crystallization dampen the efficiency of convection over time. Multiphase mixing by sinking crystal plumes initiated at thermal boundary layers has been studied extensively as the source of convection in magmas (Bergantz and Ni 1999; Dufek and Bachmann 2010; Simakin and Bindeman 2012). The focus of these models is to look at crystal fractionation (Bergantz and Ni 1999; Dufek and Bachmann 2010; Solano et al. 2012) and crustal assimilation (Simakin and Bindeman 2012) as sources of chemical differentiation for magmas emplaced in the crust.

At this stage, a model that involves a thermo-mechanical solution to the crust-magma body systems and that also solves for the internal fluid mechanics and chemical evolution of the reservoir is yet to be completed [see Gregg et al. (2013) and Degruyter and Huber (2014) for preliminary attempts to do this]. Some of the processes associated with the mechanical response to these reservoirs to mass addition or withdrawal can operate over much faster timescales than the thermal and chemical evolution of the reservoir and coupling these various processes proves to be challenging.

The mechanical and chemical interaction between the many phases that coexist in magma reservoirs is complex. It is generally parameterized with simple constitutive empirical relationships (e.g., drag terms, interphase heat transfer coefficients). Because models are only as good as the assumptions they rely on, it is important to revisit the validity of these constitutive models that introduce the dynamical coupling between the different phases in continuum models. In truth, these interaction terms involve length scales that are far beyond the resolution of continuum models, and it is the interplay between crystals, melt, and, possibly, gas bubbles that drive most of the dynamics in magma bodies. One example for such an empirical law is the mechanical separation between crystals and the melt by gravity, which controls chemical differentiation. Below a threshold crystallinity, gravitational settling laws based on Stokes settling are generally assumed valid (Martin and Nokes 1989). However, it assumes that crystals are far apart and never interact hydrodynamically and that there is no return flow of melt to conserve the mass flux over any horizontal surface in the magma body, i.e., the magma body has an infinite volume. The mechanics of multi-particles gravitational settling is quite complex though and requires an account of both short- and long-range hydrodynamic corrections even at crystal volume fractions significantly below 10% (Koyaguchi et al. 1990, 1993; Faroughi and Huber 2015). Another example lies in a proper way to account for the thermal energy balance in multiphase media (crystals+melt for example). At the continuum scale, the two general approaches commonly used assume either that. over the resolution of the numerical model, the two phases are in thermal equilibrium (local thermal equilibrium conditions or LTE) as in Simakin and Bindeman (2012) or that the contrast in thermal properties requires a thermal lag or inertia between the phases parameterized with inter-phase heat transfer coefficients (local thermal non-equilibrium conditions or LTNE) as in Dufek and Bachmann (2010). Although the latter thermal model (LTNE) is more realistic and provides better results (the current state-of-the-art), it has several issues that sometimes precludes its use for multiphase heat transfer (e.g., Karani and Huber 2015). There is therefore a need to study multiphase magmas at the discrete scale and develop new continuum-scale conservation laws that do a more careful job when filtering out heterogeneities during volume averaging (see Frontiers section).

A possible solution to resolve the actual thermal and mechanical interactions in multiphase magmas is to resort to granular (discrete) scale models as a complement to these continuum scale studies. The fundamental difference between continuum and granular scale models is that the latter explicitly solve for the mass (including chemical exchanges), momentum (e.g., drag), and energy (e.g., heat transfer) exchange between the different components through their interfaces, while the former involves a parameterization of these interactions or, in some cases (e.g., effective medium theory), relies on the definition of an effective medium with hybrid (parameterized) properties. The granular-scale approach offers a significant advantage by solving a more realistic model and allows us to study complex feedbacks arising from the interaction of multiple phases, but at the cost of heavy computational requirements.

Several “granular” multiphase models have been developed over the last decade to study magma dynamics. For example, multiphase fluid dynamics modeling where crystals are introduced using the discrete element method (DEM) coupled to Stokes flow solvers has been applied to magma chamber dynamics recently by Furuichi and Nishiura (2014) and Bergantz et al. (2015), considering mostly how crystal melt mechanical interactions affect settling in magma chambers. Granular-scale models, using the lattice Boltzmann technique, have also permitted to study the migration of exsolved volatile bubbles in magma chambers at high crystal content (see Fig. 4; Parmigiani et al. 2011, 2014; Huber et al. 2012b). As mentioned above, discrete models (for the dispersed phase) are valuable in that they provide an accurate description of the actual physics that governs the mass, momentum, and energy exchange between phases, but they are generally limited by computer power to small computational domains and require upscaling strategies to extrapolate the dynamics to the scale of the reservoir. Upscaling multiphase dynamics is and will remain a great challenge in the years to come to better constrain the chemical and dynamical evolution of magma bodies in the crust (see Frontiers topic at the end of this review).

After reviewing most of the techniques and tools used to study magma reservoirs in the Earth’s crust, the next section will focus on a specific model of magma reservoir evolution, the so-called “Mush Model” (Bachmann and Bergantz 2004; Hildreth 2004; Huber et al. 2009), strongly influenced by the pioneering efforts of many previous researchers (e.g., Smith 1960, 1979; Hildreth 1981; Lipman 1984; Bacon and Druitt 1988b; Brophy 1991; Sinton and Detrick 1992 to cite a few). Although we fully admit that this choice is strongly dictated by our previous experience, we believe it helps providing a framework to many of the questions we listed in the introduction. Obviously, the model needs many refinements, and we will try to bring out the caveats as we see them. We are also aware of the controversial aspects that remain to be resolved (e.g., Glazner et al. 2008, 2015; Tappa et al. 2011; Simakin and Bindeman 2012; Rivera et al. 2014; Streck 2014; Wotzlaw et al. 2014; Keller et al. 2015), but hope that such a review can provide a stepping stone to bring the discussion forward.

Since volcanic rocks are, at times, relatively crystal poor [particularly the compositional extremes, basalts and rhyolites, whereas intermediate magmas are often more crystal rich (Ewart 1982)], researchers have tended to draw magma reservoirs as dominantly liquids [big pools of liquid magmas for both mafic and silicic units, see Yellowstone reservoir depicted on Figure 1 (Fig. 9–14 of McBirney 1993; Irvine et al. 1998; Maughan et al. 2002)], which easily caught on in the general public’s views of magma reservoirs. However, as magmas are emplaced in contact with a colder crust, heat loss to the surrounding crust limits the time magmas remain in a mostly liquid state in these reservoirs. Hence, many researchers argued that solidification zones [or fronts, (Marsh 1981, 1989b, 2002; Marsh and Maxey 1985; Gutierrez et al. 2013; Fig. 5)] are likely to develop quickly, and form a crystal-rich buffer zone between the hottest part, most liquid part of the reservoir, and the subsolidus wall rocks. Since, many lines of evidence have suggested that magmas are dominantly kept as high-crystallinity bodies in the Earth’s last few tens of kilometers. These high-crystallinity bodies are often referred to as “crystal mushes.” Miller and Wark (2008) defined crystal mush as a “mixture of crystals and silicate liquid whose mobility, and hence eruptibility, is inhibited by a high fraction of solid particles.” (Miller and Wark 2008). This is the case for subvolcanic reservoirs (Hildreth 2004; Bachmann et al. 2007b) in large silicic systems, but also at Mid Ocean Ridges (Sinton and Detrick 1992). The supporting evidences for the importance of crystal mushes stem from petrological, geophysical, geochemical, and geological observations, as well as thermal modeling. These are detailed below.

As the temperature contrast with the surrounding crust diminishes, the cooling rate also decrease and magmas spend more time at temperature close to their solidus (Marsh 1981; Koyaguchi and Kaneko 1999; Huber et al. 2009). The waning of convection also plays a role decreasing the cooling rate with decreasing T. As stirring occurs, it steepens the temperature gradient at the edge of the body, and tends to hasten the cooling. Since convection only occurs at low crystallinity, it enhances the cooling near the liquidus (see Fig. 6 of Huber et al. 2009). Finally, latent heat released during crystallization provides energy that needs to be removed before subsequent cooling is possible. If the temperature vs. crystallinity diagram is relatively linear, latent heat is released equally over the cooling interval between the liquidus and solidus, and does not favor a prolonged state at either the low- or high-crystallinity ranges. However, the temperature vs. crystallinity relationship is quite complex and non-linear for many types of magmas, particularly those with near-eutectic points. For example, evolved high-SiO₂ magmas rapidly reach the haplogranite eutectic and are thus characterized by a strongly non-linear phase diagram [e.g., dacitic magmas reach near-eutectic conditions at 40–50 vol% crystals (Bachmann et al. 2002)]. The latent heat is mostly released when the magma approaches this eutectic condition and provides an additional thermal buffering effect, i.e., slows down the cooling rate significantly from that point onward [a process called “mushification” by Huber et al. (2009)]. Magma recharges can provide a sufficient amount of enthalpy to slow or even reverse the cooling trend generally observed in these reservoirs, which can prolong their existence at low-melt fraction (Annen and Sparks 2002; Annen et al. 2006; Leeman et al. 2008; Annen 2009; Gutierrez et al. 2013; Gelman et al. 2014), see Figure 6 for an attempt to illustrate the temperature-time evolution in the warmest parts of the reservoirs. The maturation of a magmatic field through repeated intrusions of magmas appears to also play an important role in priming the thermal state of the crust to host long-lived active magma reservoirs (e.g., Lipman et al. 1978; de Silva and Gosnold 2007a; Lipman 2007; Grunder et al. 2008).

The eruption of large, homogeneous crystal-rich units, such as the Fish Canyon Tuff (Whitney and Stormer 1985; Lipman et al. 1997; Bachmann et al. 2002), the Lund Tuff (Maughan et al. 2002), or La Pacana Tuff (Lindsay et al. 2001) have also been pivotal in the development of the Mush Model (see below). Such “monotonous intermediates” (Hildreth 1981) are common in large magmatic provinces, and unambiguously document the presence of large, homogeneous zones of crystal-rich magmas in the upper crust, potentially remobilizable by eruption. Their homogeneity is striking, and requires a relatively thorough stirring prior to eruption, likely triggered by a recharge and remobilization prior to eruption [see below, Huber et al. (2012a) and Parmigiani et al. (2014) for reviews].

Large, low-seismic velocity zones and conductive MT areas beneath active volcanoes or volcanic provinces also suggest the presence of near-solidus magmas bodies. Although the resolution remains poor (hundreds of meters to kilometers), the anomalies (both in terms of electric conductivity and seismic velocities) and the size of those bodies are not supportive of dominantly molten bodies (e.g., Dawson et al. 1990; Lees 1992; Weiland et al. 1995; Steck et al. 1998; Zollo et al. 1998; Masturyono et al. 2001; Zandt et al. 2003; De Natale et al. 2004; Husen et al. 2004; Hill et al. 2009; Waite and Moran 2009; Heise et al. 2010; Farrell et al. 2014; Ward et al. 2014). Estimates of melt range from a few percent to nearly 50% melt in some areas of the anomalies (Heise et al. 2010; Farrell et al. 2014). In addition to those seismic and MT anomalies, Bouguer gravity anomalies under older volcanic fields also suggest the accumulation of low-density, silicic rocks in the upper crust (e.g., SRMVF, Fig. 7, Plouff and Pakiser 1972; Lipman 1984; Drenth and Keller 2004; Drenth et al. 2012; Lipman and Bachmann 2015).

The presence of crystal-rich, cumulate blocks erupted in volcanic units (Wager 1962; Arculus and Wills 1980; Heliker 1995; Ducea and Saleeby 1998), as well as large plutonic masses in crustal sections with clear crystal accumulation zones (Voshage et al. 1990; Greene et al. 2006; Jagoutz and Schmidt 2012; Otamendi et al. 2012; Coint et al. 2013) suggest the presence of mush zones within the crust in active magma provinces. Cumulate signatures are not only seen in deep crustal zones, but also in shallow plutonic bodies, both using bulk chemistry (McCarthy and Groves 1979; Bachl et al. 2001; Miller and Miller 2002; Gelman et al. 2014; Lee and Morton 2015) and/or textural features (Beane and Wiebe 2012; Graeter et al. 2015).

Over the next sections, we will try to revisit the questions asked in the introduction, and see how mush zones (or the “mush model”) can help accounting for the main observations that we have laid out.

First of all, we observe many silicic crystal-poor deposits at the surface of the Earth (see Hildreth 1981; Lindsay et al. 2001; Mason et al. 2004; Huber et al. 2012a). Hence, the extraction of high-SiO₂, cold (<800 °C), viscous melt from crystalline residue in large quantities (several hundreds of km³) clearly happens in many cases. In some locations, these crystal-poor pockets of erupted magmas are homogenous (Dunbar et al. 1989; Ellis and Wolff 2012; Ellis et al. 2014), in others they grade into crystal-rich, typically hotter intermediate magmas at the end of the eruption (upper parts of the deposits; Lipman 1967; Hildreth 1981; Bacon and Druitt 1988a; Wolff et al. 1990; Deering et al. 2011b; Huber et al. 2012a; Lipman and Bachmann 2015; Fig. 8).

Generating these eruptible silicic pockets is challenging, as extracting such viscous melt [up to 10⁵–10⁶ Pa s, even with several wt% dissolved volatiles (Scaillet et al. 1998b)] from a network of millimeter-sized crystals is, at best, sluggish (McKenzie 1985; Wickham 1987) and the process can be easily disrupted. For example, convection currents should not occur, as this will significantly stir the whole magma reservoir, and lead to re-homogenizing of the crystal-melt mixture (i.e., impeding crystal-melt separation). Hence, a possibility, particularly in long-lived, incrementally built, sill-like magma bodies, is that melt is slowly extracted by gravitational processes (hindered settling, microsettling, compaction) as the magma reaches the rheological lock-up [~50 vol% (Brophy 1991; Thompson et al. 2001; Bachmann and Bergantz 2004; Hildreth 2004; Solano et al. 2012)], an extraction process potentially enhanced by gas filter-pressing (Sisson and Bacon 1999; Bachmann and Bergantz 2006; Ellis et al. 2014; Pistone et al. 2015). This melt extraction following lock-up in large mush zones has the advantage of: (1) preventing stirring and mixing of the crystal-melt mixture by convective currents, and (2) keeping the crystal-poor melt pocket in a warm environment, leading to significantly reduced crystallization rate and more time for crystal-liquid separation. Moreover, a sill-like geometry (see Fig. 9) optimizes extraction as it reduces the average vertical distance traveled by the viscous melt. The interstitial melt extraction from mushes seems not to be restricted to evolved magmas, but may also occurred at more mafic, less viscous compositions, such as basalts and andesites (Dufek and Bachmann 2010).

Petrological observations in many large, zoned ignimbrites support such a magma reservoir model, with a cap of crystal-poor material immediately underlain by a more crystal-rich zone. Some well-exposed deposits such as the 7700 BP Crater Lake, or 1912 Katmai eruptions show abrupt changes in crystallinity, from nearly crystal free early in the eruption, to 40–50 vol% crystals late in the sequence, with a very similar melt composition throughout the stratigraphy (Bacon and Druitt 1988a; Hildreth and Fierstein 2000; Bacon and Lanphere 2006). Other large ignimbrites present more gradual variations, as exemplified by famous Bishop and Bandelier Tuffs (Hildreth 1979; Hildreth and Wilson 2007; Wolff and Ramos 2014) among many others (see Hildreth 1981 and Bachmann and Bergantz 2008a for lists of these deposits). Such zoned ignimbrites are best explained by erupting a melt-rich cap followed by partial remobilization and entrainment of a silicic cumulate from the mushy roots of the system, chocking the eruption (see Smith 1979; Worner and Wright 1984; Deering et al. 2011b; Bachmann et al. 2014; Sliwinski et al. 2015; Wolff et al. 2015; Evans et al. 2016). This partial cumulate remobilization is likely a consequence of hotter recharge at the base of the silicic cap, as suggested by textural, mineralogical, and geochemical features indicative of mixing and reheating (Anderson et al. 2000; Wark et al. 2007; Bachmann et al. 2014; Wolff and Ramos 2014; Forni et al. 2015).

Other processes, such as incomplete/partial mixing between different upper crustal magma pockets and recharge (e.g., Dorais et al. 1991; Mills et al. 1997; Eichelberger et al. 2000; Bindeman and Valley 2003) or co-erupting physically separated and chemically different melt-rich lenses (Gualda and Ghiorso 2013a), have also been suggested to explain the ubiquitous chemical zoning in ignimbrites. Although we do not question that such processes are likely playing a role in generating chemical complexities in the erupted deposit, these hypotheses fail to account for the presence of co-magmatic crystal-rich cumulate blocks that erupt concurrently with the crystal-poor parts of the reservoirs [see also Ellis et al. (2014) for a case in a hot-dry rhyolitic system, and “The unzoned, crystal-poor to intermediate deposits” section below]. Hence, we argue that magma recharge reaching the base of one or several crystal-poor caps above a cumulate zone (as depicted in Fig. 9) is the most likely reservoir geometry to explain such deposits.

A second conspicuous group of ignimbrites are the “monotonous intermediates,” consisting of homogeneous crystal-rich, typically dacitic units, without any evidence for significant crystal accumulation (Hildreth 1981; Lindsay et al. 2001; Maughan et al. 2002; Folkes et al. 2011c; Huber et al. 2012a). Such units typically display corrosion textures on the low-temperature minerals (quartz and feldspars; Fig. 10) and pervasive reheating prior to eruption (Bachmann and Dungan 2002; Bachmann et al. 2002). Smaller eruptions of crystal-rich magmas [e.g., Montserrat (Murphy et al. 2000), Kos Plateau Tuff (Keller 1969; Bachmann 2010a), Taupo Volcanic Zone (Molloy et al. 2008), Cascade arc (Cooper and Kent 2014; Klemetti and Clynne 2014)] also show a very similar reheating and partial melting signal following recharge from deeper, hotter magmas (although melting evidence can be subtle or even nonexistent in some cases). The physical processes acting upon this remobilization process have been explored by several papers (Bacon and Druitt 1988b; Druitt and Bacon 1989; Couch et al. 2001; Huber et al. 2010a, 2011; Burgisser and Bergantz 2011; Karlstrom et al. 2012). The reheating signal and the partial melting signature in minerals favor a model whereby the reactivation results from the combination of melting and pore pressure increase in the mush in response to melting and the slow and progressive disaggregation of the weakened mush, a process coined as thermo-mechanical remobilization by Huber et al. (2010a). This model is also consistent with several striking characteristics (whole-rock homogeneity, high crystallinity, partial corrosion of minerals) that such deposits display (Huber et al. 2012a; Cooper and Kent 2014; Parmigiani et al. 2014).

Thermo-mechanical remobilization of locked areas following recharge is likely to be a very common process in incrementally built magma reservoirs. In most cases, such remobilization events would not lead to an eruption, but could potentially be inferred from surface deformation, change in gas outputs, and/or seismic signals related to the emplacement of the recharge. When partial melting of crystal-rich, mostly anhydrous material, occurs, the new melt released is dry, and will impact the rheological response and phase diagram of the newly produced mixture (Evans and Bachmann 2013; Caricchi and Blundy 2015; Wolff et al. 2015). As drier melt increases the temperature of mineral precipitation (e.g., Johannes and Holtz 1996), the amount of melting will tend to be limited in most cases (Huber et al. 2010a; Wotzlaw et al. 2013). Hence, the condition for eruption during remobilization might only be reached when significant heat and volatiles (mainly H₂O) are injected by the incoming recharge. Exactly how much heat and volatiles need to be injected strongly depends on several factors (including relative size of recharge and mush, geometry, compositions of magmatic end-members, average crystallinity of the mush, …), which will vary from case to case and are difficult to predict.

There is a third type of ignimbrite that commonly appears in volcanic sequences, particularly in areas where magmatism is relatively dry (hot spot or rift environments such as Yellowstone—Snake River Plain or Taupo Volcanic Zone). In such cases, homogeneous ignimbrites with crystal contents on the order of 10–20% are commonplace (Dunbar et al. 1989; Nash et al. 2006; Ellis and Wolff 2012; Ellis et al. 2013). Ubiquitous in such units are millimeter- to centimeter-sized glomerocrysts of pyroxene-plagioclase-oxides (Fig. 11). These glomerocrysts are much more mafic (when analyzed in bulk) than the bulk rocks they are found in, and entirely lack sanidine and quartz, while they are found as phenocrysts in the host units. However, the chemical compositions of their pyroxenes and plagioclase crystals are identical to the individual phenocrysts around them. These glomerocrystic pyroxene and plagioclase also show very evolved REE pattern, suggesting growth from rhyolitic melts (Ellis et al. 2014).

Within the framework of the mush model, such deposits are best understood as pockets of melt-dominated regions in relatively hot/refractory mushes built in such dry magmatic environments. Upon recharge from below, these dry mushes are not likely to melt as easily as the ones that contain sanidine and quartz (as found in wetter, colder subduction zone environments). Hence, the zoning produced by partial melting of cumulates does not happen, and the deposits remain homogenous. As gas sparging can also play an important role in mush defrosting (Bachmann and Bergantz 2006; Huber et al. 2010b), mushes in magmatically dry environments are more difficult to remobilize because exsolved volatiles fail to reach the critical volume fraction required for efficient heat transfer by advective transport (Huber et al. 2010b), and produce mostly this remarkable type of unzoned, crystal-poor to intermediate deposit (Ellis et al. 2014; Wolff et al. 2015).

Compositional gaps, or the paucity of erupted products within a range of compositions, are commonplace in volcanic sequences and called the Bunsen-Daly gaps, in homage to early work from Bunsen (1851) and Daly (1925) on Iceland and Ascension Island. Such gaps can be seen in different elements (major and trace), in all series, including tholeiitic trends (e.g., Thompson 1972) and subduction zones (Brophy 1991; see Fig. 12). Possible hypotheses to explain these gaps dominantly revolve around three main lines of reasoning: (1) two main magma types (one mafic, one felsic/silicic) are generated in the upper parts of our planet, mostly by melting mantle and crustal components (e.g., Bunsen 1851; Chayes 1963; Reubi and Blundy 2009); (2) two magma types are generated by liquid immiscibility (e.g., Charlier et al. 2011); and (3) magma evolution is dominated by crystal fractionation from mafic parents (generating a continuum of melt composition) but due to some mechanical trapping (potentially enhanced by nonlinear temperature-crystallinity relationships), magma reaching the surface are dominantly clustered at some compositions (Marsh 1981; Grove and Donnelly-Nolan 1986; Brophy 1991; Bonnefoi et al. 1995; Francalanci et al. 1995; Grove et al. 1997; Freundt-Malecha et al. 2001; Thompson et al. 2001; Dufek and Bachmann 2010; Melekhova et al. 2013).

Although these hypotheses can all play a role in different magmatic systems around the world, the mechanical trapping of melt, as seen within the context of the mush model, should be dominant in large magmatic provinces where such gaps are obvious (ocean islands, subduction zones in thin oceanic crust, or thin continental crust such as the Aegean and Taupo volcanic zones) and crustal melting/liquid immiscibility are likely not efficient. The key mechanical reasoning for trapping melt is the following: at low melt fractions (0 to ~50 vol% crystals), melt and crystals are not easily separable, due to: (1) the stirring effect of magma convection, and (2) the short time that these magmas spent at low crystallinities (Marsh 1981; Koyaguchi and Kaneko 1999; Huber et al. 2009; Dufek and Bachmann 2010). Hence, melts are dominantly extracted from the crystal cargo only after significant crystallization, showing then a composition that is very different from the mixture composition (Deering et al. 2011a). As an example of a geological setting where the amount of crustal contamination should be very limited (Sliwinski et al. 2015), the composition of erupted magmas in Tenerife displays three principal modes in SiO₂ (mafic, intermediate, and felsic magmas) likely indicating melt compositions released from: (1) the mantle (the most primitive), (2) the lower crustal differentiation zone (intermediate magmas), and (3) the upper crustal differentiation zone (the most felsic magmas). Although it was claimed by Melekhova et al. (2013) that the magnitude of the gaps should diminish with time (on the basis of purely thermal calculations), this does not appear to be true in magmatic provinces around the world (e.g., Tenerife; Sliwinski et al. 2015).

How compositionally and temporally related (i.e., same age range, see section on timescales below) volcanic and plutonic units fit onto a common framework has been actively debated for nearly two centuries [discussion started by James Hutton at the 18th century, and pushed further in classical books by Lyell (1838), Daly (1914), and Bowen (1928); and papers or memoirs such as Read (1948, 1957), Buddington (1959), Hamilton and Myers (1967), and Pitcher (1987)]. Several recent review papers summarize the main issues at play (Bachmann et al. 2007b; de Silva and Gosnold 2007a; Lipman 2007; Glazner et al. 2015; Keller et al. 2015; Lipman and Bachmann 2015). Focusing on felsic/silicic magmas (mafic plutons are clearly mostly considered cumulates for decades; e.g., Wager et al. 1960), the main hypotheses revolve around whether intermediate to silicic plutons are mostly the expression of crystallized melts (“failed eruptions,” i.e., have not undergone significant crystal-liquid separation), or cumulate leftovers (“crystal graveyards”) from volcanic eruptions (those two end-member hypotheses not being mutually exclusive in plutonic complexes containing multiple facies). Textural analysis of plutonic lithologies and trace element geochemistry (in both bulk rock and minerals) should be the two main lines of arguments to differentiate between these competing hypotheses.

1. If crystal accumulations occurs in areas of the crust: compatible trace elements should show significant increase in their concentrations, while incompatible elements should remain low (e.g., Mohamed 1998; Bachl et al. 2001; Wiebe et al. 2002; Kamiyama et al. 2007; Deering and Bachmann 2010; Turnbull et al. 2010).

2. Plutonic rocks should show some mineral orientations or preferential alignment/foliation (e.g., Wager et al. 1960; Arculus and Wills 1980; Shirley 1986; Seaman 2000) related to magmatic processes (e.g., hindered settling and/or compaction), which will tend to orient anisotropically shaped crystals as they accumulate.

Despite several attempts to look in detail at this issue, this topic remains controversial (see de Silva and Gregg 2014; Glazner et al. 2015; Lipman and Bachmann 2015). The root of the controversy may largely stem from the fact that the amount of crystal/melt segregation is less important in evolved (silicic) magmas compared to more mafic compositions [i.e., more trapped melt component in the intermediate to silicic crystal cumulates, due largely to viscosity difference (Bachmann et al. 2007b)]. Hence, silicic units have more subtle geochemical or textural signatures of crystal accumulation (or melt loss) than their mafic counterparts, and can easily be overlooked (see Gelman et al. 2014, and Fig. 13). However, although recent publications involving compilations of geochemical data from large databases (e.g., Glazner et al. 2015; Keller et al. 2015) cannot clearly differentiate volcanic from plutonic compositions at the felsic/silicic end of the spectrum, other compilations do see significant variabilities with volcanic units being on average richer in SiO₂ (Lipman 1984, 2007; Halliday et al. 1991; Gelman et al. 2014; Deering et al. 2016). Recent studies focusing on specific field examples show that several intermediate to silicic plutonic lithologies have high-compatible/low-incompatible element concentrations in bulk rock (although such units clearly still have significant trapped melt components; e.g., Bachl et al. 2001; Barnes et al. 2001; Deering and Bachmann 2010; Lee and Morton 2015; Fig. 13). Additionally, detailed textural data, using electron backscattered electron diffraction (EBSD) techniques, on non-metamorphosed granitoids (i.e., observed mineral orientation must be magmatic) indicate crystal alignment compatible with compaction/hindered settling, even for very viscous magma compositions (Beane and Wiebe 2012; Graeter et al. 2015).

While high-SiO₂ rhyolites are not rare in the volcanic record, including several very large units (>100 km³; Lipman 2000; Christiansen 2001; Mason et al. 2004), granites sensu stricto appear much less commonly. The upper crustal plutonic record is dominantly granodioritic/tonalitic in composition (Taylor and McClennan 1981; Rudnick 1995). When compiling data from large geochemical databases, such as the NAVDAT (http://www.navdat.org), one notices that low-Sr rocks are, relatively speaking, much more abundant in the volcanic realm than in the plutonic record (Parmigiani et al. 2016), an observation already discussed a quarter of a century ago on a much smaller database (Halliday et al. 1991; Cashman and Giordano 2014).

A possible explanation for the unusually high volcanic/plutonic ratio of evolved magmas is the fact that rhyolites are generated and stored in the upper crust (Tuttle and Bowen 1958; Gualda and Ghiorso 2013b; Ward et al. 2014; Lipman and Bachmann 2015), where volatiles can exsolve in quantity. If exsolved magmatic volatiles can accumulate in such melt pools, it will render them highly eruptible (increasing the gravitational potential energy of the magma) and will tend to promote more voluminous eruptions (Huppert and Woods 2002). As a consequence it will also lead to a deficit in the plutonic record. Using multiphase numerical simulations of exsolved volatile transport in magma reservoir, Parmigiani et al. (2016) have shown that bubbles tend to accumulate in crystal-poor environments. In crystal-rich environments, the volume available for the exsolved volatile phase is significantly reduced (porosity) by crystal confinement and promotes bubble coalescence and the formation of buoyant fingering pathways. Once the pathways are established, the very viscous melt does not limit the migration efficiency of the vapor phase and transport is efficient. In contrast, in crystal-poor environment, fingering channels are not stable and break into bubbles under the action of capillary forces. Discrete bubbles can only rise as fast as the viscous melt can be moved out of their way, which significantly reduces the migration efficiency of exsolved volatiles in crystal-poor magmas. As a result, bubbles tend to accumulate in the convecting rhyolitic cap, providing additional potential energy to drive or sustain future eruptions.

The accumulation of a buoyant volatile phase in these high-SiO₂ melt pools favors large eruption, but it is unlikely to prevent entirely the formation plutonic counterparts as it only brings these magmas closer to a critical state. In fact, zones of evolved high-SiO₂ granites do seldom occur in plutonic series. They appear to be mostly present in the upper parts of plutons, near the roof. Typical examples of such evolved pockets have been described in plutons from Nevada (Bachl et al. 2001; Barnes et al. 2001; Miller and Miller 2002), Klamath Mountain, California (Coint et al. 2013), Sierra Nevada Batholith, California (Putirka et al. 2014), and Peninsular Range Batholith, Mexico (Lee and Morton 2015).

The chemical composition of the atmosphere is tied to some extent to magmatic degassing (e.g., CO₂ and H₂O cycles). As magmas are dominantly forming plutons [(plutonic/volcanic ratios are typically at least 5:1 to 10:1, and possibly even greater in some areas (Crisp 1984; White et al. 2006; Ward et al. 2014), the cycle of degassing and outgassing in plutons exerts an important control on the chemistry of the surface reservoirs on Earth. The protracted periods of storage and slow crystallization at high crystallinity predicted by the mush model are prone to lead to particularly efficient magma degassing at shallow depths. In such slowly cooled crystal-rich environments, exsolved magmatic volatile build up in the mush by second boiling (exsolution produced by the crystallization of dominantly anhydrous minerals), and progressively fill more of the constricted pore space leading to a more efficient outgassing (Parmigiani et al. 2011; Huber et al. 2012b, 2013). Once the pore space, and more specifically the pore throats, are significantly reduced in size (down to the scale of a few micrometers), capillary stresses can even become large enough for gas pathways to deform locally the granular structure in the mush and volatile-filled “dikes” may finalize the outgassing at large crystal fractions (>80% vol; Holtzman et al. 2012; Oppenheimer et al. 2015), possibly leading to the formation of the ubiquitous aplitic dykes seen in plutons (e.g., Jahns and Tuttle 1963; Candela 1997). In either case, the combination of high-crystal content and increasing exsolved volatiles volume fraction can lead to the formation of efficient degassing pathways even during quiescent periods (sometimes potentially preserved as residual porosity in granites, e.g., Dunbar et al. 1996; Edmonds et al. 2003; Lowenstern and Hurvitz 2008).

The mass balance of volatiles released during volcanic eruptions can provide valuable insights into the state of shallow magma reservoirs prior to the eruption. As such, the mismatch in S mass balance between petrological inferences (i.e., comparing the S content of pre-degassed melt inclusions and fully degassed matrix glass) and remote sensing spectroscopy or sulfur output estimated from ice core data can provide additional clues about the dynamics that prevail in magma chambers and more specifically about the factors that control volatile exsolution. The notion of “Excess S” refers to the underestimation of the sulfur output from petrological constraints (melt inclusion record) during an eruption. This Excess S paradox has been evident since the eruption of El Chichon in Mexico (Luhr et al. 1984) and described in detail for the eruption of Mt. Pinatubo in 1991. At Pinatubo, remote sensing methods estimated the S output to nearly 20 Mt (e.g., Soden et al. 2002), while the petrological method provided an estimate smaller by a factor of 10 (Gerlach et al. 1996). Since these two eruptions, the release of S in excess to petrologic estimates seems to be more the rule than the exception, especially for silicic magmas in arcs (see Fig. 14 and Wallace 2001; Shinohara 2008).

Although multiple hypotheses have been suggested to account for this Excess S (see Gerlach et al. 1996 for a review), the presence of a S-rich exsolved volatile phase in the magma reservoir seems to be the most commonly accepted postulate (Gerlach et al. 1996; Scaillet et al. 1998a; Wallace 2001; Shinohara 2008). Crystallization-driven (second boiling) exsolution in a shallow magma reservoir can generate significant excess S where the S trapped by vapor bubbles during crystallization is not accessible to melt trapped subsequently as inclusion in crystal hosts. Although second boiling is bound to play an important role on generating Excess S in crystal-rich magmas, in crystal-poor silicic eruptions the deficit in S in melt inclusions requires an alternative process. In light of the previous section, accumulation of exsolved (S-rich) volatile bubbles in the erupted, low-crystallinity portions of magma reservoirs would provide a closure to the problematic missing S in crystal-poor silicic magmas, These bubbles (poorly constrained, but likely up to 20–30 of percent by volume) can be extremely rich in S, due to the very high affinity of S for the exsolved volatile phase (Zajacz et al. 2008). Hence, the mush model, predicting efficiently degassed crystal-rich zones and gas-charged, highly eruptible pockets of magma, provides a consistent framework to explain the missing source of S necessary to close the volatile budget for crystal-poor units.

Caldera-forming events typically happen after long periods of maturation in long-lived magmatic provinces during which andesitic volcanism dominates for up to several millions of years (Lipman 2007; de Silva and Gosnold 2007b; Grunder et al. 2008; Deering et al. 2011a) (Fig. 15). During this andesite-dominated period, the crust initially warms up to become more amenable to host silicic magma reservoirs in the mid to upper crust (Jellinek and DePaolo 2003; de Silva and Gregg 2014). Once this mature stage is reached, it is not rare to see several caldera-forming events occurring in relatively short time spans [~2–3 Myr; multicyclic caldera systems (Hildreth et al. 1991; Graham et al. 1995; Christiansen 2001; Lipman and McIntosh 2008; Lipman and Bachmann 2015)]. During such caldera cycles, magmas erupting shortly (~5–30 kyr) prior to the caldera-forming event typically appear to be compositionally similar to the climactic event (e.g., Bacon and Druitt 1988b; Bacon and Lanphere 2006; Bachmann et al. 2012), whereas the post-caldera magmas appear more mafic and drier (Shane et al. 2005; Bachmann et al. 2012; Gelman et al. 2013a; Barker et al. 2014; Fig. 15). These petrological shifts can be explained by the fact that the leftover, crystal-rich (mush) portion of the magma system following the eruption will be strongly affected by the caldera-forming event (Hildreth 2004). As those eruptions typically lead to significant decompression of the underlying crustal column (potentially up to 50–100 MPa, accounting for the overpressure necessary for the eruption to start and the redistribution of mass following the removal of much of the eruptible parts of the reservoir, see Bachmann et al. 2012), some rapid and significant degassing and crystallization is expected in the mush, particularly if it is near-eutectic and volatile-saturated (Bachmann et al. 2012). For mushes that were volatile-saturated prior to the decompression, the caldera eruption will likely promote the rapid formation of gas channels resulting in deep-gas release. Hence, the leftover mush is expected to become more crystalline (decompression-driven crystallization) and more degassed than it was before the caldera eruption (Degruyter et al. 2015). It will take another long maturation stage for the system to be primed for a new caldera event with a shallow, gas-charged, evolved magma body.

How long does a magma body stay above the solidus in the Earth’s crust (i.e., potentially active) remains one of the key questions in magmatic petrology and volcanology. The longevity, of course, impacts the extent of magmatic differentiation (including ore formation) and the eruptive volume that is available at any given time. As discussed above, magmatic reservoirs are mostly kept as high-crystallinity regions above the solidus when active, but some controversies remain as to how long the system remains in such a mush state (with estimates from a few thousand to more than 1 000 000 years). The crux of the debate lies in the fact that different analytical techniques provide different answers, and likely date different processes. Zircon dating of silicic plutons and volcanic rocks record extended crystallization histories, ranging from tens of thousands to several millions of years (Reid et al. 1997; Brown and Fletcher 1999; Lowenstern et al. 2000; Reid and Coath 2000; Charlier et al. 2003, 2005, 2007; Schmitt et al. 2003, 2010b; Coleman et al. 2004; Vazquez and Reid 2004; Bacon and Lowenstern 2005; Simon and Reid 2005; Bindeman et al. 2006; Matzel et al. 2006; Bachmann et al. 2007a; Miller et al. 2007; Walker et al. 2007; Reid 2008; Claiborne et al. 2010; Klemetti et al. 2011; Storm et al. 2012, 2014; Walker et al. 2013; Wotzlaw et al. 2013, 2014; Chamberlain et al. 2014b). In contrast, crystal size distributions (CSD; Pappalardo and Mastrolorenzo 2012), and diffusion timescales (see Druitt et al. 2012 and the compilation of data in Cooper and Kent 2014) are typically much shorter (typically months to a few hundreds of years).

As most elements used in diffusion timescale modeling move relatively fast (e.g., Fe-Mg in mafic minerals, Ti in quartz, Sr and Mg in plagioclase, Li in zircon …), it is not surprising that preserved zoning patterns in such elements indicate relatively short timescales. Using simple scaling, one can provide a rough estimate of maximum timescale that diffusion re-equilibration would provide if one knows the size of the crystals and the diffusion coefficients (time ~R²/D, where R is the radius and D the diffusivity). Elements that are diffusing slower (REE, Ba) can be used to obtain longer timescales (e.g., Morgan and Blake 2006; Turner and Costa 2007), but, in such cases, periodic resorption, concurrent crystal growth, and 3D effects associated with the finite size of the host lead to significant complications in estimating time. For example, growth of feldspars in the Yellowstone magma reservoirs (Till et al. 2015) demonstrably led to “diffusion-like” profiles that look similar for elements that are diffusing at very different rates. Similarly, CSD studies typically assume linear growth rates (no pauses in crystallization, no resorption event), with values obtained from experiments (hence typically faster than what actually happens in silicic magma reservoirs). Hence, crystal size distributions provide a lower bound, and sometimes strongly underestimated, timescale.

The short timescales obtained by diffusion modeling are best interpreted as reactivation/remobilization/mixing timescales (e.g., Martin et al. 2008; Matthews et al. 2012; Till et al. 2015). They, however, do not provide a timescale of magma reservoir growth and storage in the crust. For such information, the zircon age distributions or bulk U/Th/Ra/Pb ages likely provide more reliable estimates (Reid et al. 1997; Cooper and Reid 2003; Bachmann et al. 2007a; Turner et al. 2010; Cooper and Kent 2014; Guillong et al. 2014). The difference between eruption age (estimated using the young zircon population or the Ar/Ar age when available) and the oldest co-genetic zircons (i.e., antecrysts, see Bacon and Lowenstern 2005; Miller et al. 2007 for a definition) is typically at least 10⁴–10⁵ years (see Costa 2008; Reid 2008; Simon et al. 2008; Bachmann 2010b for reviews on the topic). Xenocrysts (foreign crystals coming from unrelated wall rock lithologies) can also occur, but, in some cases, they can be isolated and removed from the data set (as they can be significantly older, plot off the Concordia and/or have different trace element chemistries; Miller et al. 2007; Lukács et al. 2015). Of course, the presence of older but co-genetic zircons (antecrysts) does not mean that the system remained above the solidus for the whole time; some parts of the system could have reached the solidus, and later recycled by reactivation events (Schmitt et al. 2010a; Folkes et al. 2011a). Hence, physical models are needed to address this issue.

Numerical models of the thermal state of such systems usually agree with relatively long timescales of tens to hundreds of thousands of years, particularly when they include incremental recharge at relatively high fluxes (e.g., Annen et al. 2008; Gelman et al. 2013b). Moreover, it is very costly, both in terms of volatiles and energy, to reactivate sub-solidus plutons; nearly all volatile elements (including CO₂ and H₂O; Caricchi and Blundy 2015) and all latent heat has been lost to the surrounding. If the subsolidus material has to be melted again, the thermal budget required involves not only the sensible and latent heat loss, but also a large fraction of energy wasted on reheating solid material that will not undergo melting (Dufek and Bergantz 2005). The assimilation of these devolatilized sub-solidus portions of the body will also induce a depletion of volatile content in the reservoir. From a mechanical standpoint, entrainment/recycling of material that is not completely solid is much more likely, as it will disaggregate more easily, particularly if some local melting at grain boundaries is involved (Beard et al. 2005; Huber et al. 2011).

In summary, the periodic recharge of magma at reasonably high fluxes (>10⁻⁴–10⁻³ km³/yr; Annen et al. 2008; Gelman et al. 2013a, 2013b; Karakas and Dufek 2015) is likely to maintain at least parts of silicic magma reservoirs in a mush state for periods extending to at least several hundreds of thousands of years, even in the upper part of the crust (8–10 km depth). Shorter timescales, obtained with other techniques (particularly diffusion modeling), mainly refer to periods of reactivation following significant heat and mass addition into the reservoirs (Cooper and Kent 2014; Till et al. 2015). Hence, crystal/liquid separation is likely to have enough time to occur, despite being sluggish, to some extent in those reservoirs, driving differentiation toward more evolved and more eruptible crystal-poor magma pods (e.g., Bachmann and Bergantz 2004, 2008c; Hildreth 2004; Lee et al. 2015). We note that those thermal models are lacking the mechanical part of the energy budget (overpressurization during recharge, depressurization during eruption), which can have a significant impact on the stability and duration of reservoirs (Degruyter and Huber 2014; Degruyter et al. 2015) and the feedbacks with crustal rheology (Gregg et al. 2013; de Silva and Gregg 2014). For example the exsolution of volatiles affect the effective thermal properties of the magma by reducing its thermal capacity and conductivity (Huber et al. 2010b) and consuming latent heat. We note, however, that adding this mechanical part to the models is unlikely to significantly shorten the lifetimes of the reservoirs.

Over the last paragraphs, we have tried to assemble many important questions in volcanology/petrology into a coherent framework. The Polybaric Mush Model (Fig. 9) provides a context to explain several seemingly unconnected observations, such as the geophysical imaging of crystal-rich zones in magmatically active zones at different levels in the crust, the chemical differentiation of magmas dominated by fractional crystallization (with some assimilation, explaining some of the variability in isotopic data) of mafic parents, the production of compositional gaps in volcanic series, the presence of cumulates in the plutons records, and the volatile mass balance in shallow reservoirs. However, there are many questions that remain to be answered, and the next section tries to outline some of the challenges that our community faces in the years to come.

At the most fundamental level, the goal is to continue developing a more integrated picture between physical models, geophysical/geodetical inversions and petrology/geochemistry. In particular, we should tend to obtain a multiphase and multiscale (time and space) picture of magmatic systems, from the formation of the primary melts to their final resting place in the crust or at the surface (including interactions with the atmosphere). To reach such an ambitious goal, we would need, among many things:

1. To revisit the information that can be extracted from the inversion of geophysical data sets. The best resolution that one can obtain from a magma body at depth is much lower than most structures of interest in these dynamical environments. As such, it is important to understand how the implicit filtering operates and how it impacts the set of questions that can be logically addressed with these studies. For example, if the target of the analysis is to estimate melt fraction from v_p/v_s tomography, even with good experimental calibration, one is left to wonder what the estimates of a melt fraction at the scale of hundreds of meters to even kilometers really mean. Is the outcome of the inversion to be understood as a volumetric average (as is commonly assumed)? The effective mechanical properties of heterogeneous media, which is what one measures through inversions, are generally quite different from a volumetric average. For complex, but structured, heterogeneous media, the effective properties inferred from elastic properties for example, result from a complex non-linear optimization process. For instance, it is quite possible that the effective medium elastic properties that one retrieves from tomography are dominantly controlled by heterogeneities that are volumetrically secondary. We certainly need more effort to relate the best-fit results of these inversions to the physical reality of complex multiphase magma bodies. This starts with a better understanding of the role of subgrid scale structures on the effective properties at the resolution of the inversion.

2. To continue developing and improving our interpretation of tracers for rates of processes in magmatic systems. For example, it is necessary to improve precision and spatial resolution of mineral geochronology; and do more diffusion modeling of measured gradients to determine timescales, and focus on what they mean. It also becomes possible to use stable and radiogenic isotopes to constrain thermodynamics and kinetics, although the most substantial effort in this direction so far relates to mafic magmas (Richter et al. 2003; Watkins et al. 2009b; Huang et al. 2010; Savage et al. 2011). We also need better constraints on possible growth rates of minerals and bubbles. Dissolution events are transparent to the present chronometers and leaves hiatus in the temporal evolution of magma bodies that are impossible to quantify at the present time. Can we find kinetic tracers that would provide information about the duration and frequency of resorption/dissolution events?

3. To better constrain disequilibrium/kinetic effects in magma reservoirs. Most of the petrology and trace element partitioning among phases in magma reservoirs assumes at least a local equilibrium. The mere fact that crystallization (and sometimes dissolution) takes place requires a finite degree of disequilibrium, even over small (crystal sizes) scales. To relate dynamical processes to the rock record, a kinetic framework is necessary. There is ample information to mine from heterogeneities and disequilibrium; these clues will become crucial in testing dynamical models and establishing what set of processes controls the temporal evolution of these reservoirs over different spatial and temporal scales. Shall the high-T geochemistry and petrology community follow the low-T geochemistry community and invest time in developing reactive transport modeling for trace elements and isotopes in magma bodies? The chemistry can become very complicated, but the idea even in the context of oversimplified models is worth pursuing.

4. To focus on eruption triggering mechanisms. One of the fundamental societal goals in volcanology is to understand the causes (triggers) that drive volcanic eruptions. There are several possible factors that can influence the state of a shallow reservoir and influence the timing of an eruption. Triggers can be external [e.g., caused by a regional earthquakes, roof collapse, or magma recharges from deeper (Sparks et al. 1977; Pallister et al. 1992; Watts et al. 1999; Eichelberger and Izbekov 2000; Gottsmann et al. 2009; Gregg et al. 2012)] or internal [e.g., buoyancy, pressure buildup during recharge or gas exsolution (e.g., Caricchi et al. 2014; Degruyter and Huber 2014; Malfait et al. 2014]. Fundamentally, the notion of trigger, which is related to an event or a process that has a causal link to a subsequent eruption, remains quite loose. Establishing a causality between a trigger and an eruption is much more challenging than establishing a correlation between the two.

5. To continue providing better models of long-standing issues such as the “room problem.” Although the “room problem” is alleviated by having long-term incremental addition in a mainly ductile crust, pushing material (a) to the side, (b) down, as dense cumulates from the lower crust founder back into the mantle (Kay and Mahlburg Kay 1993; Jull and Kelemen 2001; Dufek and Bergantz 2005; Jagoutz and Schmidt 2013) and (c) upward, as upper crustal reservoir dome up during resurgence (e.g., Smith and Bailey 1968; Lipman et al. 1978; Marsh 1984; Hildreth 2004), the question continues to drive controversies (e.g., Glazner and Bartley 2006, 2008; Clarke and Erdmann 2008; Paterson et al. 2008; Yoshinobu and Barnes 2008).

At the present time, many questions remain to be answered about the migration, emplacement, evolution, and ultimately eruption of silicic magmas. The main challenges stem from the fact that: (1) we have access to only a limited set of indirect observations, in most cases difficult to relate to the actual processes that control the chemical and dynamical evolution of these magmas, and (2) most magmas never reach the surface (trapped as plutons within the crust).

The development of a self-consistent model for silicic magmas residing in the upper crust is complex and requires a joint effort across multiple disciplines. The mush model presented in this review provides, we believe, a testable hypothesis from which better models can be developed. In spite of what remains to be done to understand the fate of magmas in the crust, the mush model offers a self-consistent paradigm that is supported by several independent observations from geophysics, geochemistry, and petrology. For instance the development of mush zones is consistent with:

• Thermal models of incrementally growing magma bodies and geophysical data on active systems require that such reservoirs are dominantly crystal-rich on long timescales (>100 kyr; Marsh 1981; Koyaguchi and Kaneko 1999; Bachmann and Bergantz 2004; Cooper and Kent 2014).

• Chemical differentiation can occur internally by melt-crystal separation and produce shorter-lived eruptible pockets (with melt content higher than 50 vol%) from time to time. Melt extraction can largely occur when convection currents wane down, i.e., when mixing by stirring becomes less efficient than gravitational settling. Phase separation should be most efficient at intermediate (40–70%) crystallinities.

• Melt extraction in slowly cooled mush zones (receiving periodic enthalpy additions by recharges and thermally buffered by latent heat at high crystallinity) helps explaining the longevity of magma centers/lifetime of volcanic fields (as inferred from zircon age populations) and the comparatively short duration over which magmas are mobile and prone to erupt (Cooper and Kent 2014).

• The presence of cumulate plutonic lithologies in well-exposed upper crustal batholiths (e.g., Bachl et al. 2001; Barnes et al. 2001; Putirka et al. 2014) and as cognate clasts in silicic volcanic units (Gelman et al. 2014; Lee and Morton 2015; Wolff et al. 2015) confirms the close spatio-temporal link between fossil mush zones (“crystal graveyards”) and their extracted melt-rich complements.

Despite our attempt to cite the largest amount of previous work (we believe the reference list is extensive with nearly 500 citations), we are certainly not giving enough due credit to the amazing work that people have done over the years in our topic of interest. We apologize for this in advance. Most of these ideas put forward here have been nucleating in the womb of our own biased view. We have them here for discussion, and trust our community that some or most of them will be challenged. We look forward to it. We thank Mikito Furuichi, Naomi Matthews, Marc-Antoine Longpré, and Katy Chamberlain for providing figures that we use to illustrate some recent advances in our field, and Keith Putirka for motivating us to write this paper. Comments from Ben Ellis, Wim Degruyter, and Peter Lipman on a previous version of the manuscript helped better crystallize some ideas that were too tangled up in a mushy discussion. We are indebted to Charlie Bacon and Shan de Silva as formal reviewers, for the insightful and constructive comments, as well as to the very efficient editorial handling by Fidel Costa.