INTRODUCTION

Magma bodies spend the majority of their supra-solidus existence in the crust as high-crystallinity mushes; in such a high-crystallinity state, the lower degree of thermal disequilibrium with their surroundings, the inability to convect and latent heat buffering efficiently reduce their cooling rate (Marsh, 1981 ; Koyaguchi & Kaneko, 1999 , 2000 ; Huber et al. , 2009 ). It is therefore expected that new magmatic inputs rising from the mantle and deep crust will often interact with crystal mushes. This interaction between a mature, high-crystallinity magma body and fresh intrusions plays a central role in the construction and protracted existence of large magmatic systems (Mahood, 1990 ; Reid et al. , 1997 ; Lipman, 2007 ).

As magma bodies grow in the crust, they perturb the local stress field and act as a focusing lens for new injections of magma (Karlstrom et al. , 2009 ). When two magmas enter into contact, the dynamics and the amount of mixing between them depends on their respective viscosity, density and the injection rate (Huppert & Sparks, 1980 ; Jellinek et al. , 1999 ; Snyder, 2000 ). In the ubiquitous case where the mush is more evolved than the intrusion, the more mafic and dense magma input spreads underneath the mush and acts as a heating plate with little chemical mixing as a result of density stratification and rheological contrast (Couch et al. , 2001 ). As the lower magma cools and partially crystallizes, it exsolves and releases buoyant volatiles that will migrate upwards to the contact between the two magma bodies. The rate of volatile injection into the mush is controlled by the cooling rate and convective motion of the underplating magma (Eichelberger, 1980 ; Huppert et al. , 1982 ; Cardoso & Woods, 1999 ). The collection of volatiles at the magma interface (rheological barrier) leads to the formation of a buoyant instability that can eventually overcome the rheological barrier and enter the mush once the volume fraction of gas and therefore the buoyancy of the foamy layer reaches a critical value (Eichelberger, 1980 ).

Bachmann & Bergantz ( 2006 ) hypothesized that the volatiles exsolved from an underplating magma would enhance the remobilization of crystal mushes because of their ability to increase the vertical heat transfer. Their numerical model does not solve for the enthalpy coupling between the two magmas (the intrusion was treated as a fixed temperature boundary condition) and therefore does not solve for the injection rate of volatiles from the intrusion to the mush. The crystal mush was also not allowed to melt during the course of their calculations.

Building on this earlier model of gas sparging, this study addresses the following questions:

1. What features control the mobility of volatiles through a mush?

2. Can the vertical transport of volatiles enhance or prevent the defrosting (partial melting of the crystalline framework) of a mush (see Mahood, 1990 )?

3. What is the effect of the mush composition on the thermal evolution of both magma bodies?

4. What explains the paucity of crystal-rich ignimbrites in ‘dry’ (extensional, hotspot) tectonic environments?

To address these questions, we have developed a one-dimensional numerical model that solves for the combined thermal evolution of the two magma bodies. We also solve for the exsolution rate and transport of volatiles from the intrusion into the mush. Enthalpy conservation in the mush allows us to calculate the evolution of the porosity during partial melting of the mush locally and accounts for changes in permeability by using a parameterization of the temperature–crystallinity relationship. The calculations span a range of pressures (from 1 to 4 kbar at the level of the intrusion), a range of compositions for the mush (through temperature–crystallinity parameterization) and a range of initial water contents for the intrusion (from 1· 5 to 6 wt % water).

PHYSICAL MODEL

We solve for the heat and volatile mass transfer between a hot intrusion and a crystal mush (see Table 1 for mathematical symbols and Fig. 1 for a schematic illustration of the physical model). For simplicity, the intrusion is emplaced instantaneously and the volume ratio between the mush and the intrusion is kept constant at a ratio of 4:1. The colder crystalline mush is initially set at a uniform crystallinity (55 vol. %) and temperature (which depends on the crystallinity–temperature relation used). As the underplating magma cools at the contact with the mush (cooling from below the intrusion is neglected here, for simplicity), it partially crystallizes, can reach volatile saturation and starts exsolving a volatile phase. The volatile phase consists mostly of water, as we assume that a large fraction of any CO ₂ in the magma has already degassed before reaching the depth of the intrusion (≤4 kbar). The whole-rock composition of the intrusion is assumed to be andesitic and we use a different crystallinity–temperature relationship for the mush to understand the importance of its composition on the thermal evolution of both magmas. The mush is assumed to be initially saturated with volatiles and we neglect any volatile dissolution or exsolution that is expected to arise as a consequence of both reheating and partial melting of the mush. Furthermore, we assume that volatile transfer from the andesitic intrusion to the mush carries heat but does not change the stability of the different minerals (i.e. we do not track the activity of water in the volatile phase). Neglecting changes in the stability field of different mineral phases in the mush upon the injection of volatiles is a reasonable assumption as the kinetics associated with the growth or dissolution of minerals caused by changing the exsolved gas composition in the mush are much slower than melting by heat transfer.

Fig. 1.

Open in new tab Download slide

Schematic representation of the coupled mush–intrusion system. The two magmas exchange heat and the mush is subjected to injections of volatiles exsolved from the underlying intrusion.

THE UNDERPLATING MAGMA

We assume that the crystal mush provides the main resistance to heat flow in the underplating magma–crystal mush system. We consider only heat transfer between the intrusion and the mush and do not account for cooling through the surrounding wall-rocks. Therefore, all the calculations that follow provide underestimates of the cooling rate, and as a consequence, the exsolution rate for the intrusion. We simplify the heat balance for the underplating magma:

(1)

where the first term on the right-hand side is sensible heat and the second term is latent heat; the subscript i refers to the intrusion, H _i , L _i , c _i , ρ _i and χ _i are respectively the thickness, the latent heat of crystallization, the specific heat, the density and the crystallinity of the underplating magma body and q _out is the heat transfer from the intrusion to the mush. This heat balance equation assumes a well-mixed intruding magma with negligible thermal gradients; this is a relatively good approximation for a convecting magma body with strong temperature-dependent viscosity. The dynamic viscosity of the magma in the intrusion depends both on temperature and crystallinity (Dingwell et al. , 1993 ):

(2)

where Q is the activation energy (fixed here at 500 kJ/mol), R is the ideal gas constant, T _ini is the reference temperature (equal to the temperature during emplacement) and χ _cr = 0·5 is the critical crystallinity at which the intrusion behaves like a solid. To calculate q _out we need an expression for the thermal gradient at the interface between the two magmas and thus the boundary layer thickness. The convection of a fluid with strongly temperature-dependent viscosity is parameterized in terms of the Rayleigh number:

(3)

where ρ is the average density of the convecting fluid, α is the thermal expansion coefficient, κ is its thermal diffusivity, H is the thickness of the fluid layer, and g is the acceleration due to gravity. This definition of Ra implies that the boundary conditions for the intrusion are prescribed temperatures, which is a good approximation as long as the time step in the calculations is much shorter than the cooling timescale. In all the following calculations, Ra never decreases below 10 ⁵ in the intrusion. The assumption of a well-mixed intrusion magma is therefore valid. The temperature difference responsible for convection (Δ T _e has given by Davaille & Jaupart ( 1993 )):

(4)

The boundary layer thickness δ between the two magmas depends on the Rayleigh number of the andesitic magma body:

(5)

where δ ₀ and Ra ₀ are the initial boundary layer thickness and Rayleigh number. The heat exchange between the two magmas becomes

(6)

where k _cond is the thermal conductivity at the interface and χ _mush is the crystallinity of the mush. q _out is the thermal coupling between the two magmas and requires knowledge of the thermal gradient on the other side of the interface (in the mush). Figure 2 shows schematically the heat balance at the interface between the two magmas described by equation ( 6 ). Equation ( 6 ) solves for the conjugate heat transfer between the intrusion and the mush, which is limited by the medium with the slower heat transfer, in our case the mush because the intrusion is able to convect (Carrigan, 1988 ).

Fig. 2.

Open in new tab Download slide

Schematic description of heat transfer at the interface between the intrusion and the mush. The heat flux out of the mush is balanced by changes in sensible heat (temperature) and melt fraction in the mush at the interface.

It is, however, important to note that equations ( 1 ) and (6 ) relate the phase diagrams of the two magmas through the variation of crystallinity inside the andesitic intrusion and at the interface with the mush. For the andesitic magma, we use MELTS with the major element composition listed by Parat et al. ( 2008 ) for the Huerto andesite (San Juan Volcanic Field, Colorado) for a range of initial water contents (from 1·5 to 6 wt % H ₂ O) and for pressures of 1, 2 and 4 kbar). We also use MELTS to parameterize the exsolution of volatiles (water only in this case; see below) as a function of temperature (see Fig. 3 for an example of crystallinity–temperature and exsolution parameterization for an underplating intrusion at 2 kbar containing initially 3·5 wt % H ₂ O).

Fig. 3.

Open in new tab Download slide

(a) Calculation of the crystallinity–temperature relationship for an andesite [major element composition of the Huerto andesite from Parat et al. ( 2008 )] with 3·5 wt % H ₂ O at 2 kbar. (b) Calculation of the volume fraction of exsolved water for the same magma composition and same pressure conditions. Both calculations were run with MELTS (Ghiorso & Sack, 1995 ). We use linear fits over a range of temperature from 750 to 850°C to parameterize the crystallization and exsolution of water in the intrusion.

The absence of CO ₂ in the MELTS calculation influences both the saturation temperature and the stability field of the various crystalline phases, and is expected to lead to more crystallization of anhydrous phases at higher temperature as CO ₂ decreases water activity. Because we are focusing on comparing the dynamics of the coupled system for different mush–intrusion conditions, we neglect these effects in the following calculations. We nonetheless obtain very reasonable saturation temperatures; for the calculation at 4 kbar and 4·5 wt % H ₂ O, the melt reaches H ₂ O saturation at about 830°C, which is consistent with the experiments of Parat et al. ( 2008 ).

The density of the exsolved water is calculated with the Modified Redlich–Kwong (MRK) equation of state of Halbach & Chatterjee ( 1982 ). We use a fit to the results obtained with the MRK equation to avoid inverting the equation of state for the density at each node and for each iteration (see Fig. 4 ). We obtain

(7)

for 600 ≤ T ≤ 900°C and 300 ≤ P ≤ 4000 bar. The density of the volatile phase (water) in the intrusion varies by a factor of eight between the calculations at 1 and 4 kbar (below and above the supercritical point for water).

Fig. 4.

Open in new tab Download slide

Parameterization of the equation of state for water. We use the modified Redlich–Kwong model of Halbach & Chatterjee ( 1982 ).

We assume that the transport of exsolved volatiles after the saturation point can be described by a multiphase Darcy equation:

(8)

where k and k _r are the permeability of the magma and the relative permeability for the volatile phase, respectively, and Δρ is the density contrast between the melt and the volatile phase. Gas transport in a crystal slurry, as opposed to a rigid mush, is not expected to exhibit a Darcian behavior. Instead, at low Reynolds number (ratio of inertia to viscous forces), equation ( 8 ) has the general form of the Stokes’ ascent velocity for a dilute suspension of bubbles, where k and k _r would then represent an average bubble radius. Moreover, we expect that the presence of at least 25–30 vol. % of crystals and the large viscosity contrast between the bubbles and the melt (∼10 ^–8 ) justifies this simplification. The closure of equation ( 8 ) requires expressions for both the permeability and the relative permeability. The permeability is related to the crystallinity by the Carman–Kozeny relation (Bear, 1972 ):

(9)

where the constant A is set to 2 × 10 ^–12 m ² , leading to a permeability k = 10 ^–12 m ² at χ _i = 0·5. This choice of permeability–porosity relationship is consistent with the range of porosities observed in our calculations (0·45–0·75). Many empirical and theoretical expressions have been derived for the dependence of the relative permeability on the volatile volume fraction; these commonly reduce to power laws with exponent n ∼ 4 (e.g. Brooks & Corey, 1964 ). We use the following expression:

(10)

where S _g is the volume fraction of exsolved volatiles. Equation ( 10 ) neglects the effects of irreducible and residual saturation. As a correction to the non-purely Darcian behavior of the gas in the crystal-slurry, we allowed gas transport below the residual saturation limit. This simplification is not valid for the mush, where the viscosity of the magma is much larger as a consequence of its dacitic composition and high crystallinity.

Crystal mush

The mush above the intrusion initially has a crystallinity of 55% and experiences both reheating and injection of volatiles from the intrusion. We solve for volatile mass conservation:

(11)

where ρ _g = ρ _H2O , χ _mush is the crystallinity of the mush, S _g is the pore volume fraction occupied by the volatiles and u _{d, z} is the z -component of the Darcy velocity for the volatile phase and z is parallel to gravity. The boundary condition at the interface between the mush and the intrusion requires matching the volatile flux with equation ( 8 ). We assume here that the melt is volatile-saturated and that no further degassing or dissolution of volatiles occurs within the mush. The flux of volatiles is obtained from

(12)

where g is the acceleration due to gravity and is counted positive upwards, k and k _mush and k _{r mush} are respectively the permeability and relative permeability of the mush for the volatile phase and μ _H2O is the dynamic viscosity of the volatile phase, which, for simplicity, was fixed at 10 ^–5 Pa s (Assael et al. , 2000 ).

The permeability of the mush evolves with time because of the changes in porosity associated with partial melting of the mush. We use the same permeability–porosity relationship as equation ( 9 ). The relative permeability accounts for the residual saturation; that is, the volume fraction of the pore space under which the mobility of the non-wetting phase is drastically reduced because of capillary effects:

(13)

The residual saturation can be understood as a threshold value below which the mobility of the non-wetting phase is negligible (see Fig. 5 ). It therefore plays a central role in volatile transport in the mush. Owing to the lack of experimental constraints on the residual saturation for exsolved water hosted in a feldspar-dominated mush saturated with a silicate melt, we use a value of S _residual = 0·2, consistent with many natural systems (Bear, 1972 ).

Fig. 5.

Open in new tab Download slide

Schematic representation of capillary effects responsible for the residual saturation of a non-wetting phase. (a)–(c) show three snapshots of the motion of a bubble (non-wetting phase) through a constriction calculated with the lattice Boltzmann method. The buoyancy force works against bubble deformation and the displacement of the wetting phase next to the solid. (d) A train of buoyant bubbles ascending through a saturated porous medium. At low saturation, the non-wetting phase is dispersed as bubbles, each of which is working against the forces described in (a)–(c). (e) At higher saturation (above the residual saturation) the mobility of the non-wetting phase is increased by several orders of magnitude because the connected phase mobility now is controlled by the viscosity of the nonwetting phase (about 8–10 orders of magnitude smaller for silicic magmas and exsolved water).

The crystallinity of the mush is related to temperature using a parameterization with a single free parameter b :

(14)

T _sol and T _liq are the solidus and liquidus temperatures, respectively set to 700 and 950°C, and b is a free parameter used to set different crystallization behaviors for the magma. For example, as b decreases to zero, the magma tends to crystallize mostly over a narrow range of temperature just above the solidus (see Fig. 6 ).

Fig. 6.

Open in new tab Download slide

Crystallinity–temperature relationships for the mush from equation ( 14 ) for different values of b . As b decreases towards zero, the system tends to crystallize more just above the solidus temperature. The crystallinity–temperature relation for the Fish Canyon magma body (dacite) calculated with MELTS is shown as a reference. The Fish Canyon magma is consistent with an exponent b < 0·4.

To follow the evolution of the crystallinity of the mush, we solve the energy balance that includes heat diffusion, absorption of latent heat during the partial melting of the mush and advection of heat by the volatile phase:

(15)

In equation ( 15 ), the specific heat of the volatiles c _g is assumed constant (3880 J/kg K; Lemmon et al. , 2003 ), ρ _H2O is calculated with the MRK equation of state [equation ( 7 )], and ρ _s is the average density of the crystalline phase subjected to melting. Equation ( 15 ) assumes that the volatile phase is locally in thermal equilibrium with the melt and crystals, which assumes a low Peclet number (i.e. Pe = u _{d, z} R /κ < 1, R being the average pore radius and κ the average thermal diffusivity). The maximum Peclet number achieved in the simulations is about 10 ⁻² for a pore size of about 1 mm. The thermal equilibrium assumption is, thus, a good approximation. The density, specific heat and thermal conductivity of the solid–melt–volatiles mixture are calculated with

(16)

(17)

(18)

where the subscripts s, m, g refer to solid, melt and volatiles, respectively. Table 1 lists the parameters used in the physical model.

Numerical model

We solve for the one-dimensional gas transport and enthalpy balance between the two magmas using an iterative explicit finite-difference method. We use a nested iterative scheme that can be divided in two groups, ( 1 ) the enthalpy balance group and ( 2 ) the volatile transport group. We use an iterative solver to calculate the heat transport and thermal coupling between the two magmas, the crystallinity profile in the mush, the average crystallinity of the intrusion and the permeability profile in the mush (enthalpy group). Once the variables converge to within a tolerance (relative difference between two consecutive iteration steps is lower than 10 ⁻¹² typically), we fix them and solve iteratively for the volatile transport, relative permeability of the mush, and thermal and mechanical properties of the volatile–melt–solid mixture, until convergence is reached. We then repeat the iterative process for the two parts of the scheme consecutively before proceeding to the next time step. Table 2 lists the major input variables and the output that is calculated.

We use similar thermal and mechanical properties, as well as similar initial conditions for each calculation. We conducted four calculations with different power-law exponents for the crystallinity–temperature relationship [equation ( 14 ), b = 0·4, 0·5, 0·7 and 1 for each intrusion initial water content (1·5, 2·5, 3·5, 4·5 and 6 wt % H ₂ O) and each pressure of emplacement (1, 2 and 4 kbar)].

All the calculations have a 2 km thick mush discretized with a resolution of 20 m. We set the initial crystallinity of the mush to 0·55. Because we compare calculations with the same initial crystallinity in the mush but with different temperature–crystallinity relationships, the initial temperature of the mush is not the same from one run to another. The initial temperature of the 500 m thick intrusion is set to 850°C, about 20°C above the water saturation conditions of the Huerto andesite composition at 4 kbar (Parat et al. , 2008 ). The calculations are run for 1000 years with a time step of a few hundred seconds to ensure stability. The mush is assumed to contain no exsolved volatiles initially. Each run is duplicated with an identical test run except that in the test run the volatiles exsolved from the intrusion are not allowed to penetrate into the mush. These control runs allow us to isolate and quantify the impact of volatiles on the overall thermal evolution of the system.

RESULTS

For each run, we keep track of the temporal evolution of the intrusion in terms of average temperature, crystallinity and volume fraction of exsolved volatiles. Figure 7 shows the temporal evolution of the average temperature of the intrusion and the temperature just above the interface between the two magmas (in the mush) for an intrusion emplaced at 1 kbar, containing 4·5 wt % H ₂ O and a mush with b = 0·4. The initial temperature difference is about 116°C and is quickly reduced to about 10°C over c . 10 years because of the high heat flux from the intrusion to the mush. After about 10 years, the temperature in the intrusion and in the lowest part of the mush both decrease, because of the cooling of the intrusion. The temperature difference at the interface between the two magmas is reduced to about 1°C after about 1000 years. Figure 8 shows the evolution of the average crystallinity and exsolved volatile fraction in the intrusion for the same calculation (1 kbar, 4·5 wt % H ₂ O). At 1 kbar, the intrusion contains about 23 vol. % volatiles at the emplacement temperature (850°C). The volatile volume fraction is quickly reduced because of the rapid transfer to the mush [volatile mobility depends strongly on the volatile volume fraction from equations ( 8 ) and (10 )].

Fig. 7.

Open in new tab Download slide

Average temperature of the intrusion and the intrusion–mush boundary with time. The andesitic intrusion is emplaced at a pressure of 1 kbar and contains 4·5 wt % water. The mush crystallinity–temperature is set with a power-law exponent b = 0·4.

Fig. 8.

Open in new tab Download slide

Evolution of the average crystallinity and volume fraction of exsolved volatiles in the intrusion. The characteristics of the mush and the intrusion are similar to those in Fig. 7 .

We also monitor the evolution of porosity, temperature and pore volume fraction occupied by the volatiles with depth in the mush. Figure 9 shows results after 1000 years for a mush with a temperature–crystallinity relationship with b = 0·4 and for an intrusion emplaced at 1 kbar with 4·5 wt % H ₂ O. Figure 9a shows the porosity (or extent of melting) and Fig. 9b the temperature in the mush. The vertical distance is normalized by the thickness of the mush (2 km). The volume fraction of volatiles in the pore space ( Fig. 9d ) exhibits a more interesting behavior. First, volatiles rise through more than half of the mush (more than 1 km) and much further than the melting front (see Fig. 9a ). Second, the maximum saturation of 0·2 is consistent with our choice of residual saturation, as the volatiles become much more mobile above this threshold. The profile, however, shows a decrease in pore saturation between distances of 0·1 and 0·4. Figure 9c shows the profile of the product of the porosity and the pore volume fraction of volatiles, and shows that the volume of volatiles monotonically decreases until reaching a plateau consistent with 20% of the pore space at 45% porosity all the way to the volatile propagation front. This illustrates that the reduction of S _g between distances 0·1 and 0·4 from the intrusion is the consequence of melting and increase of pore space at constant volatile volume ( Fig. 10 ). The mobility of the volatile phase is controlled by the local volume of volatiles in the pore space [equation ( 10 )].

Fig. 9.

Open in new tab Download slide

Results obtained after 1000 years; the characteristics of the mush and the intrusion are similar to those in Fig. 7 . (a) Profile of porosity ϕ with depth in the mush (zero is the mush–intrusion boundary and one is the top of the 2 km thick mush). (b) Temperature profile. (c) Local volume fraction of volatiles (=ϕ S _g ). (d) Volume fraction of the pore space occupied by volatiles, S _g .

Fig. 10.

Open in new tab Download slide

Schematic description of the dilution effect experienced by the volatile fraction during the increase of pore space associated with melting. Case (a) shows an example where the melting rate is faster than the injection rate, leading to a dilution effect of the volatile phase in the pore space ( S _g decreases with time). Case (b) shows the opposite situation, where the volatile injection rate is high enough to compensate for the increase in pore space associated with melting. In this case, the local pore fraction occupied by the volatile phase increases with time and can potentially reach the residual saturation threshold over which it becomes mobile and buoyantly rises in the mush.

At 1 kbar, the low density of the volatile phase (0·1–0·2 kg/m ³ ) allows it to reach residual saturation efficiently and propagate through the mush ahead of the melting front. As the melting front progresses behind the volatile front, the pore space increases and more volatiles are required to keep the melting region at residual saturation. If the increase in pore volume by melting is locally faster than the flux of volatiles, S _g decreases below the residual saturation value (see Fig. 10 ), leading to the formation of two large-scale slug regions of volatiles separated by a low hydraulic conductivity region (below residual saturation). The volatile slug that detaches ahead of the melting front (distance from the intrusion greater than 0·2 in Fig. 9b ) carries little excess enthalpy beyond the main melting front. As the crystallinity–temperature power-law exponent b increases the segmentation of the volatile-rich part of the mush is reduced ( Fig. 11 ), as the increase in pore volume owing to melting is more equally distributed over a wide range of temperature.

Fig. 11.

Open in new tab Download slide

Comparison between runs with different crystallinity–temperature power-law exponent b . The intrusion is emplaced at 1 kbar with 4·5 wt % H ₂ O.

Figure 12 shows the profiles for porosity, temperature, S _g and ϕ S _g for a similar calculation in which the intrusion is emplaced at 4 kbar. In this case, the melting front propagates ahead of the volatile front and the features observed in Fig. 9 are not present. The melting front has not propagated as much as for the 1 kbar case—this is mostly the consequence of low volatile mobility. The thermal conductivity of water-rich volatiles ( k _H2O = 0·3; Lemmon et al. , 2003 ) is about an order of magnitude lower than for the melt or crystals. Volatiles thermally insulate the mush from the intrusion because of their low thermal conductivity (Bagdassarov & Dingwell, 1994 ). When the volatiles are mobile (when they reach the residual saturation), the heat they advect upwards can compensate for the negative impact they have on diffusive heat transfer. We will explore this effect in more detail later.

Fig. 12.

Open in new tab Download slide

Results obtained after 1000 years. The characteristics of the mush and the intrusion are similar to those of Fig. 7 , except that the pressure (depth) is 4 kbar at the mush–intrusion boundary. (a) Profile of porosity ϕ with depth in the mush (zero is the mush–intrusion boundary and one is the top of the 2 km thick mush). (b) Temperature profile. (c) Local volume fraction of volatiles (=ϕ S _g ). (d) Volume fraction of the pore space occupied by volatiles, S _g .

The rate of cooling and crystallization of the intrusion depends on three main parameters: ( 1 ) the pressure at which the intrusion is emplaced, as this affects the phase stability of the different crystallizing phases; ( 2 ) the water content of the intrusion, for the same reasons; ( 3 ) the power-law exponent b of the overlying mush. The last is important as it controls the initial temperature of the mush (we fix the initial crystallinity) and also controls the partitioning of heat delivered by the intrusion into sensible and latent heat.

Figure 13 compares the average temperature and crystallinity of the intrusion (at 1 or 2 kbar) and the porosity of the mush just above the contact with the intrusion after 1000 years. Each contour plot is shown as a function of the power-law exponent b for the mush and the initial water content of the intrusion. The cooling rate of the intrusion decreases with increasing b because of the lower initial temperature difference between the two magmas at high b . The slower cooling at high b is also a consequence of the lower latent heat buffering of the mush; that is, the heat from the intrusion is partitioned primarily into sensible heat in the mush, which reduces more efficiently the temperature difference between the magmas.

Fig. 13.

Open in new tab Download slide

The two columns compare the results for the average temperature and crystallinity of the intrusion after 1000 years for calculations in which the pressure at the intrusion is set to 1 and 2 kbar. The x -axis and y -axis of these contour plots are respectively the initial water content of the intrusion and the power law exponent b , defining the temperature–crystallinity relationship in the mush. The porosity of the mush just above the intrusion is also shown for the calculations at 1 and 2 kbar.

The dependence of the final temperature of the intrusion on its water content is more complicated. At 1 kbar, the final average temperature of the intrusion increases with water content. As volatiles exsolve and migrate toward the overlying mush they carry heat and reduce the temperature difference between the intrusion and the lowest part of the mush, and thus reduce the subsequent heat flow out of the intrusion. At 2 kbar, the situation is slightly different, with a peak of high average final temperature for the intrusion at about 2·5 wt % H ₂ O. For dry intrusions (below 2 wt % H ₂ O), the cooling and crystallization of the intrusion does not lead to exsolution of water. For intrusions containing between 2 and 3·5 wt % H ₂ O, the mass of volatiles injected into the mush is relatively small. Moreover, at 2 kbar the average density of the volatile phase is almost twice the density at 1 kbar. As a consequence, the small mass of volatiles injected into the mush does not reach residual saturation and the volatiles remain immobile just above the interface between the two magmas. As volatiles have a much lower thermal conductivity than both melt and crystals they thermally insulate the intrusion and prevent efficient cooling. This ‘thermos bottle effect’ (Carrigan, 1988 ) caused by a gas-rich layer between two magmas has been previously described by Bagdassarov & Dingwell ( 1994 ). At higher initial water contents in the intrusion, the greater exsolution and injection rate of volatiles into the mush allows residual saturation to be reached and the volatiles become mobile in the mush, advecting heat away from the interface between the magmas. The porosity of the mush at the contact with the intrusion is consistent with these observations.

We define three lengthscales to characterize the melting efficiency and volatile transport in the mush for each calculation. The first lengthscale d _m0 represents the melting front position in the mush measured from the contact with the intrusion in the absence of volatiles. It is obtained from the control runs where the volatile phase is not injected into the mush. The second lengthscale d _m1 is similar, but for runs where volatiles were allowed to penetrate into the mush. The last lengthscale d _g characterizes the distance traveled by volatiles in the mush from the contact with the intrusion. We define two dimensionless ratios R ₁ and R ₂ :

(19)

represents the efficiency of volatiles to be transported into the mush, where R ₁ > 1 indicates calculations where the volatiles are transported ahead of the melting front (e.g. Fig. 9 ). The second dimensionless ratio is given by

(20)

This is the ratio of the melting front distance when volatiles are accounted for to when they are not. R ₂ therefore isolates the effect of volatiles on the melting of the mush.

Figure 14 shows contour plots of R ₁ and R ₂ for three intrusion emplacement pressures as function of b and the initial water content of the intrusion. The results for R ₁ are similar for 1 and 2 kbar. Volatile transport is more efficient for larger values of b and higher intrusion water contents. The dependence of R ₁ on b is consistent with the observation of the effect of b on the volatile volume fraction shown in Fig. 11 . At larger values of b the decrease in volume fraction induced by the increase of pore space is reduced.

Fig. 14.

Open in new tab Download slide

The left column shows the efficiency of the volatiles to rise ahead of the melting front in the mush for different pressures. The right column compares the melting front position with and without volatiles for three pressures. (For further details, see the main text.)

At 2 kbar, we observe similar results, but shifted to higher intrusion water contents and lower R ₁ as expected because of the higher density of the volatiles at higher pressures. The results at 4 kbar are different. Only the most favorable cases, low b (faster cooling and exsolution rate from intrusion) and high water content in the intrusion allow the volatiles to propagate ahead of the melting front ( R ₁ > 1).

The contour plots of R ₂ ( Fig. 14 ) characterize the impact of volatiles on the partial melting of the mush. R ₂ > 1 represents calculations where volatiles enhance the melting of the mush, whereas R ₂ < 1 indicates that volatiles suppress melting. At 1 kbar, no calculations lead to R ₂ < 1. At low water contents (< 3 wt %), no volatiles are injected into the mush and R ₂ = 1. At higher pressure (2 and 4 kbar), the same result is obtained but shifted to higher water contents as expected. Nevertheless, a major difference appears. A region where R ₂ < 1 (i.e. where volatiles decrease the melting efficiency) is observed at lower water contents and low b .

DISCUSSION

Effect of melting

Partial melting increases the hydraulic conductivity of the mush as the permeability increases with increasing porosity (melt fraction). For multiphase flows in porous media, the effect of melting on the transport of each phase depends on their wetting characteristics. For a fixed volume of non-wetting phase, its saturation scales with

(21)

where 1 – χ is the local porosity of the mush. As a consequence, the dependence of the hydraulic conductivity of the non-wetting phase on the porosity becomes

(22)

The partial melting of the mush therefore has a negative impact on volatile transport. This effect becomes even more important when the local volatile volume fraction is just above the residual saturation as observed in Fig. 9 . As a consequence, for each calculation, as the lowermost part of the mush is melting and reaching higher porosities, more volatiles have to be injected to reach the residual saturation and become mobile. This effect is clearly observed in our results, as calculations with larger b (which reduces the amount of partial melting at high crystallinity) allow volatiles to propagate further.

The mobility of the volatiles thus depends on two competing effects: melting and injection rate from the intrusion. To better understand these effects we discuss a conceptual model with a single pore (see Fig. 10 ). The volume fraction occupied by volatiles in the pore increases because of volatile injection and decreases because of the increase in pore space by melting. Using mass balance for the volatiles, we obtain

(23)

where is the initial pore volume, and q _vol is the volumetric injection rate of volatiles. The relative importance of the source and sink terms for S _g can be cast in terms of a dimensionless number ζ:

(24)

where τ _m and τ _v are respectively the timescales for melting and volatile injection. For ζ > 1, volatiles can accumulate in the pore space, reach the residual saturation and become mobile in the mush. On the other hand, when ζ < 1, melting dominates and keeps the volume fraction of the volatiles below the residual saturation. Mobility of the volatile phase (i.e. ζ > 1) requires either a small melting rate or fast injection rate of volatiles into the mush, or a combination of both. Smaller melting rates occur for larger b values. However, the slower cooling rate at high b decreases the volatile injection rate as observed in Figs 13 and 14 . For R ₁ > 1, the negative effect of melting on the mobility of volatiles has to be compensated by the high injection rate (i.e. ζ > 1). At 4 kbar, the transport of volatiles is limited by the volume of volatiles injected into the mush and requires the more water-rich conditions and the largest initial temperature difference between the two magmas.

Effect of volatiles

Another important result from the calculations is the effect of volatiles on the thermal evolution of the two magmas. When the degassing of the intrusion is such that volatiles can travel upwards through the mush (ζ > 1), the volatiles advect heat and push the melting front forward. The excess enthalpy advected by the volatiles is limited, as substantial melting of the mush would locally drive the volatile phase under the residual saturation limit (ζ < 1). For cases where the volume of volatiles injected into the mush is too low to reach residual saturation in the lowermost part of the mush, volatiles affect only the overall heat transfer because of their low thermal conductivity. They thermally insulate the mush from the intrusion and, as a result, slow down the diffusive heat transfer, i.e. R ₂ < 1 (see Fig. 15 ). In our calculations, the thermal conductivity of the volatiles, the melt and solids are set respectively to 0·3, 1·4 and 2 W/m K. As a consequence, thermal insulation associated with the presence of volatiles includes two contributions: a direct contribution from the low thermal conductivity of the volatile phase and an indirect contribution if volatiles increase the degree of partial melting of the mush locally, replacing solids with melt with a lower conductivity (1·4 instead of 2 W/m K; Murase & McBirney, 1973 ; Bagdassarov & Dingwell, 1994 ).

Fig. 15.

Open in new tab Download slide

Melting scenarios in dry and wet magmatic provinces. In the dry case, volatiles stall at the interface between the intrusion and the overlying mush and thermally insulate both magmas. In the wet case, the mobility of volatiles enhances the upward propagation of the melting front and allows the defrosting of the mush.

The magnitude of the thermal insulation associated with a volatile-rich immobile layer depends on the choice for the averaging of the mixture properties [equation ( 16 )]. We use mass averages because the product of the specific heat and the local temperature has to be integrated over the mass of the body to provide the sensible enthalpy fraction. This choice of averaging, however, provides a crude approximation for the thermal conductivity of a mixture. A more accurate estimation of the mixture thermal conductivity would be obtained using the bounds derived by Hashin & Shtrikman ( 1963 ). According to this model, the lower and upper bound for the thermal conductivity are both overestimated by respectively up to 10–20% with our mass averaging and therefore predict a stronger thermal insulation by volatiles. The calculation of these upper and lower bounds is straightforward but generates discontinuities in the mixture thermal conductivity owing to its stronger dependence on the volatile volume fraction. The discontinuity arises from the step function distribution of volatiles controlled by the residual saturation. For the sake of numerical simplicity, we thus use mass averages for the thermal conductivity. The extent of melting generated by the advection of volatiles will remain mostly unaffected by this choice, as the thermal conductivity of the mixture affects only the heat diffusion term. We emphasize that the overall trend of the results will remain similar, but the magnitude of R ₁ and R ₂ would be different.

Effect of convection

In all calculations, we did not account for the possible onset of convection within the mush once the crystallinity decreases below the rheological lock-up (χ < 0·50) in a layer thick enough to become unstable. The presence of volatiles injected from the intrusion can also lower the average density at the base of the defrosted mush and contribute to gravitational instabilities. Convection will increase the overall heat transfer from the intrusion to the mush and accelerate the cooling of the intrusion. As a result, the injection rate of volatiles to the mush is expected to increase, but the delivery rate of volatiles to the melting front in the mush can be substantially slowed down as exsolved bubbles have to separate efficiently from the convecting magma. The separation of the bubbles from the convecting magma depends on the viscosity of the magma and the average bubble size (Cardoso & Woods, 1999 ; Namiki et al. , 2003 ). Convective currents will be limited for the parameter range we considered, as efficient melting scenarios require small b values and layers of high melt fraction (above 60–65%) barely extend over a few tens of meters. Moreover, with time, the lowermost part of the mush cools with the intrusion after the efficient initial reheating, leading to decreasing melt fraction in the lowest part of the mush after a few tens of years.

Questions revisited

In the Introduction we listed several questions relative to the thermal evolution of the two magma bodies. In the light of the results of this study, we revisit these questions.

What features control the mobility of volatiles through a mush?

There are three major controls on the mobility of volatiles through a mush. If the advection of volatiles is associated with the advection of enthalpy, then their mobility is controlled by ( 1 ) the flux of volatiles, ( 2 ) the ambient pressure and ( 3 ) the crystallinity–temperature relationship in the mush. As the mobility of volatiles in the mush is controlled by the volume fraction they occupy [see equation ( 13 )], the flux of volatiles and the pressure constrain the volumetric rate of accumulation of volatiles locally, whereas the crystallinity–temperature relationship of the mush describes the increase in pore space resulting from the flux of enthalpy associated with the volatile flux. To summarize, the mobility of volatiles is controlled by the volume fraction they occupy in the mush. Partial melting of the mush decreases the mobility of volatiles (assuming volatiles are the non-wetting phase in the melt–volatile–solid system), whereas lower pressure increases their mobility.

Can the vertical transport of volatiles enhance or prevent the defrosting (partial melting of the crystalline framework) of a mush?

Our results show that once the volatiles are mobile, the excess enthalpy they advect away from the mush–intrusion interface enhances the defrosting of the mush. This, however, requires relatively low pressures, high permeability and/or water-rich (saturated) intrusions. If the intrusion is unable to deliver enough volatiles to the mush, the low thermal conductivity of the volatiles thermally insulates the two magmas and therefore prevents melting in the mush.

What explains the paucity of crystal-rich ignimbrites in ‘dry’ (extensional or hotspot) tectonic environments?

Large crystal-rich dacitic eruptions (‘Monotonous Intermediates’) are common in continental arc settings, with examples found in the high Andes (Lindsay et al. , 2001 ), and in the Oligocene magmatic ‘flare-up’ in the Western USA (Lipman, 2007 ; Christiansen & McCurry, 2008 ). On the other hand, such crystal-rich intermediate ignimbrites are yet to be found in hotspot or continental rifts (e.g. Yellowstone–Snake River Plain area; see Fig. 16 ). This observation has led some researchers to speculate that mushes do not exist in the source regions of crystal-poor rhyolites in those tectonic settings (Christiansen, 2005 ; Streck & Grunder, 2008 ).

Fig. 16.

Open in new tab Download slide

Total alkalis–silica diagram showing the wide compositional gap present in the dry, hotspot-related Yellowstone Snake River Plain (YSRP) units, whereas the wet, calc-alkaline series of the Oligocene Great Basin display a continuous range of compositions from basaltic andesite to rhyolite. We argue that this compositional gap is a consequence of the ineruptability of mushes with intermediate compositions in the YSRP because of their lower volatile content, leading to ( 1 ) a delay in volatile exsolution and ( 2 ) inefficiency of gas-driven defrosting (gas sparging) of rheologically locked mushes (modified from Christiansen & McCurry, 2008 ).

On the basis of this study, we propose an alternative scenario: dry mushes are nearly impossible to erupt. To appear in the volcanic record, a locked crystal mush needs to be defrosted. In hotspot or extensional magmatic environments, the low initial volatile content of the underplating intrusions results in slow gas transfer with a negative (insulating) impact on the remobilization of the mush ( R ₂ < 1, see Fig. 14 ). In continental arc settings, however, an intrusion initially containing more than 4 wt % of volatiles will generate a high enough injection rate of volatiles for them to become mobile in the mush (ζ > 1), enhancing defrosting ( Fig. 17 ).

Fig. 17.

Open in new tab Download slide

The different regimes of mush defrosting induced by volatile exsolution from an underplating intrusion at 2 kbar. b is the exponent of the crystallinity–temperature relationship [equation ( 14 )].

What is the effect of the mush composition on the thermal evolution of the two magma bodies?

The mush composition controls its crystallinity–temperature relationship. At a crystallinity of the mush above 0·5, the partitioning of enthalpy between sensible and latent heat during the heating of the mush by the intrusion depends on its composition. For more evolved magmas, the changes of crystallinity resulting from small additions of enthalpy can be significant ( b ≤ 0·5), whereas for magmas further away from the haplogranitic eutectic composition the changes in crystallinity can be substantially smaller and enthalpy is mostly absorbed by sensible heat (temperature change). The latter case is characterized by ζ > 1 and volatiles can rise and move beyond the melting front provided that the intrusion contains enough volatiles. For mushes with a composition closer to the eutectic, the increase in pore space associated with the heat diffused from the intrusion can, in some cases, prevent the volatiles from ascending or, when enough volatiles are available, lead to the formation of separate gas slugs as indicated in Fig. 11 .

CONCLUSIONS

In this study, we calculate the enthalpy coupling between a cooling and degassing intrusion and an overlying crystal mush. We quantify the partial melting of the mush for different crystallinity–temperature relationships, at different pressures and with different initial water contents for the intrusion. We show that partial melting of the mush reduces the mobility of the volatiles because the ability of volatiles to ascend through the mush depends on the volume fraction they occupy in the pore space. This effect couples the fate of volatiles exsolved from the underplating intrusion to the composition of the mush.

Our calculations demonstrate that volatiles can play two opposite roles in the thermal evolution of the mush depending on their mobility in the mush. At high volatile fluxes, the enthalpy they advect away from the interface between the two magmas favors remobilization. In contrast, low volatile fluxes partially prevent the reactivation of the mush by insulating it from the intrusion. These opposite effects at low and high volume fraction of volatiles can explain the paucity of erupted crystal-rich ignimbrites in dry tectonic environments where mushes of intermediate composition cannot be rejuvenated and erupted once they reached their rheological lock-up point.

FUNDING

C.H. and M.M. were supported by NSF EAR 0608885 and the Larsen Fund, and O.B. was supported by NSF EAR grant 0809828.

ACKNOWLEDGEMENTS

We thank our colleagues Joe Dufek and George Bergantz for unfailing support and numerous fruitful discussions on the topic of this paper. Insightful comments by an anonymous reviewer, F. Spera, J. Wolff, and editor Wendy Bohrson are gratefully acknowledged.

REFERENCES