Modelling of crystal fractionation in magmas predicts that the content of compatible trace elements such as strontium (highly enriched in plagioclase, the dominant mineral phase in the mid to upper continental crust) drops in the residual melt and should correlate with the SiO content4. In geochemical datasets, such as the North American Volcanic and Intrusive Rock Database (NAVDAT; http://www.navdat.org), the distribution of strontium for a given SiO range between volcanic and plutonic suites presents a striking conundrum. For dacitic/granodioritic compositions (intermediate SiO content), the frequency of rock samples with similar strontium contents is nearly identical between volcanic and plutonic rocks; for high-SiO magmas, there is a clear trend towards many more low-strontium compositions in volcanic units than in plutonic units. Clearly, geochemical datasets such as NAVDAT may carry sampling biases, but we argue that the number of samples considered (1,989 volcanic rock samples with 65–69% SiO ; 2,355 volcanic rock samples with 73–77% SiO ; 2,018 plutonic rock samples with 65–69% SiO ; and 1,647 plutonic rock samples with 73–77% SiO ; see Fig. 4a) and the prominence of the low-strontium mode in high-SiO rhyolites imply that magmas depleted in strontium are more commonly represented in the eruptive than in the intrusive record. We assume that these low-strontium magmas are formed by interstitial/residual melt extraction from dacitic crystal-rich mushes, and segregate into liquid-dominated caps15, 16, 33, 34 that are consequently more prone to erupt than to stall in the crust.
Estimates of volatile mass balance in magmas stored in the crust are generally based on the volatile content dissolved in melt inclusions35. However, melt inclusion data are limited by several factors, including: (1) leakage of volatiles (loss through volume diffusion, cracks or cleavage surfaces); and (2) an inability to record volatiles trapped in hidden reservoirs (exsolved bubbles or sulfide phases). The excess sulfur paradox is a clear consequence of such limitations. Sulfur degassing from the melt during magma ascent in the conduit does not contribute to the excess sulfur because it is typically accounted for in the petrologic estimate from melt inclusions2. The processes that govern the unbalanced sulfur budget are the build-up of an abundant sulfur-rich MVP produced by crystallization-driven exsolution, the transport and accumulation of sulfur-rich MVP from deeper untapped portions of the magmatic system, and, in some cases, the breakdown of sulfides or anhydrite before an eruption1, 2, 36, 37. In the context of the large excesses of sulfur released during the explosive eruption of crystal-poor rhyolitic caps, the contributions to the total sulfur mass budget from crystallization-driven exsolution and sulfide/anhydrite breakdown are bound to be tenuous, and call for an efficient migration and accumulation of MVP exsolved deeper down (by second boiling in the crystal-rich mushy roots of the magmatic system; see Fig. 4b).
At high crystallinity, buoyant bubbles are likely to deform along the direction perpendicular to gravity and, therefore, experience a significant hydrostatic pressure drop. Once this pressure drop is high enough to invade a pore throat, drainage is initiated and the blobs of MVP migrate vertically. The formation of anisotropic MVP clusters along the direction of gravity requires a confinement from the crystal phases to work against interfacial tension. Interfacial tension will tend to make bubbles spherical, whereas gravity will provoke the horizontal expansion of bubbles when their ascent is obstructed. Crystal confinement is therefore key to the development and stability of mobile MVP fingers in a crystal mush. In Extended Data Fig. 2, we show snapshots (Extended Data Fig. 2a–c show initial conditions and Extended Data Fig. 2d–f show the steady-state MVP distribution) for three numerical calculations conducted using a multiphase lattice Boltzmann model (that is, the interparticle potential Shan–Chen method21, 38, 39, 40). We implemented this model using the open-source Palabos library (http://www.palabos.org) and ran the simulations on the supercomputer clusters at Georgia Tech, the Swiss Federal Institute of Technology Zurich (ETHZ) and CSCS-Switzerland (the lattice Boltzmann method and code validations are described in detail below). In these calculations, the pore volume fraction of MVP and the size of pore throats are identical but with different crystallinities (crystallinity (1–ϕ) increases to the right). We see that increasing the spatial confinement (crystallinity) leads to enhanced coalescence and the formation of stable fingering. The run at highest crystallinity (Extended Data Fig. 2c, f) maximizes the vertical bubble pressure drop and is the only one that leads to the formation of a continuous fingering feature across the domain. In the other calculations, bubbles are mechanically trapped by capillary and viscous forces. A high MVP volume fraction and high crystallinity favour the formation and stability of viscous fingers, because they prevent the growth of Rayleigh–Plateau instabilities41, 42, 43 that are responsible for the break-up of fingering.
In essence, fingering pathways, once established, require little displacement of the viscous melt in the porous medium and reduce therefore the rate of energy dissipation in the melt. This results in an increase of MVP discharge in the mush. In Extended Data Fig. 2g, we report the results of a set of calculations (78 in total) conducted with the same porous medium geometry (porosity 0.4) but a varying initial spatial distribution and volume fraction of MVP. The porous medium is made up of spherical pores, each connected to six cylindrical throats along each dimension (in three dimensions). Throat radii are randomly generated to introduce a random distribution of capillary entry pressure for the MVP invasion in neighbour pores. The calculations at a given MVP pore volume fraction are repeated with a different initial distribution of MVP in the pore space.
These calculations clearly show that the MVP discharge through the porous medium increases with the MVP pore volume fraction. At low MVP volume fraction (red region), MVP remains distributed as discrete bubbles trapped in the medium because of capillary and viscous forces. The discharge is negligible. At intermediate MVP pore volume fraction (green region), coalescence becomes important and makes the formation of percolating fingering pathways of MVP possible, which leads to a sharp increase in MVP discharge. However, in this region, the connectivity of fingering pathways depends on the initial distribution of MVP. In this regime, the runs that do not yield an efficient MVP discharge often display intermittent formation and destruction of fingering pathways, leading to successive periods of short-lived efficient transport and periods of capillary and viscous trapping of MVP bubbles. Above a certain pore volume fraction (blue region), the MVP always forms and sustains percolating fingering pathways and the MVP migration rate is fast.
The rising velocity of an isolated bubble through an infinite stagnant fluid can be described by the law derived by Hadamard and Rybczynski (reviewed in refs 44 and 45). However, when it comes to finding the velocity of a single bubble rising inside a cloud of bubbles, the dynamics becomes more complex because bubbles interact hydrodynamically with each other and with the ambient melt. For example, at low Reynolds number, the rising velocity of a trailing bubble aligned with another (lead) bubble along their direction of motion is greater than that of an individual bubble (the Smoluchowski effect; see ref. 46), while misaligned bubbles experience a greater viscous drag because of the melt return flow. Recently, Faroughi and Huber18 characterized both local and non-local bubble interactions theoretically, and proposed a new hindrance function, F(Ψ,λ), which represents the ratio of the migration velocity of a bubble in a suspension to that of the same bubble in a bubble-free melt. The relative velocity of bubbles in a suspension at low Reynolds number is controlled by the balance between buoyancy and viscous stresses. The presence of bubbles decreases the hydrostatic pressure by a factor (1−Ψ), whereas the presence of a cloud of MVP bubbles dispersed in the melt affects the effective shear viscosity of the magma17. The general expression for the hindrance function is18:
where represents the drag caused by the return flow of melt, parameterized as:
Ψ and Ψ are, respectively, the bubble volume fraction and the random close packing limit for spherical bubbles; β = 0.45, being a geometrical proportionality constant determined experimentally17; captures Smoluchowski’s effect:
and expresses the change in momentum diffusivity via:
Finally, under the assumption that bubbles are inviscid relative to the melt, we obtain the relative bubble velocity:
where U and U are, respectively, the bubble velocity in the suspension and its Stokes ascent velocity. Equation (1) is plotted against experimental data over a wide range of particle volume fractions in Extended Data Fig. 1. We carried out experimental studies of bubble migration by using water injected at the top of a tank filled with silicon oil (Extended Data Fig. 3). The localized and fixed injection points at the top of the tank (water is denser than silicon oils, so buoyancy is reversed compared with the typical situation in a magma reservoir) mimic the localized point sources that will transfer the MVP from the mush to the cap. The experimental set-up allows local bubble plumes to form, where hydrodynamic interactions introduce a smaller penalty to bubble buoyant migration. It is expected that bubble plumes (‘vents’) will form out of the mush in heterogeneous magma bodies; it was therefore necessary to validate equation (1) against this set of experimental data for our MVP cap suspension model (see Extended Data Fig. 1 inset). Note the significant decrease in MVP flux as the bubble fraction increases in a suspension; this contrasts strongly with the results shown in the section ‘MVP migration in the crystal-rich mush’, where the MVP flux increases significantly with increasing volume fraction in a porous mush.
Crystal-poor caps are prone to convect, especially when buoyant bubbles are fluxed in from below. Convective motion will affect the migration of buoyant bubbles47. At low Reynolds numbers, the overall motion of bubbles can be decomposed as a vectorial sum between the imposed convective motion and the buoyant phase separation calculated above. The behaviour of bubbles is determined by the ratio of these two velocity components (sometimes parameterized as a Stokes number47). For small (millimetre-size) bubbles in a silicic magma, one can assume that bubbles remain highly coupled to the convective flow motion, except when the flow decelerates in boundary layers next to the edges of the reservoir. Thus, we adapt the model derived by Martin and Nokes48 and also used by Dufek and Bachmann34 for crystal suspensions to calculate the residence time of bubbles in the convecting cap.
The main differences between our calculations and those presented in refs 34 and 48 are: (1) in our calculations, the segregating phase comprises buoyant bubbles with free-slip conditions at the interface between bubbles and melt; and (2) we use the hindrance function derived above to correct for the presence of other bubbles, which can significantly affect the buoyant bubbles’ ability to migrate in magmas. The model we obtain for the mass (here volume) conservation of bubbles in the cap therefore reads:
where F(Ψ,λ) is the hindrance function calculated from equation (1), H is the thickness of the crystal-poor layer, and q is the volumetric flux of MVP coming from the mush.
We first solve equation (2) with q = 0, and retrieve a characteristic residence time for bubbles in a convecting magma. In Extended Data Fig. 4, we show the solution to this differential equation under magmatic conditions. We find that increasing the initial volume fraction of bubbles in a convecting magma has a positive impact on accumulation—that is, at a higher volume fraction, bubbles remain trapped in the convective motion longer because of the hindrance to phase separation. Moreover, the decay rate of the bubble fraction that remains suspended in the convecting magma no longer follows an exponential law18, 48, because of the nonlinear dependence of the MVP ascent velocity on the MVP volume fraction. We also calculate the residence time of bubbles with two arbitrary sizes over a wide range of dynamic shear viscosities of the melt, for dilute (Ψ = 0.01) and high (Ψ = 0.3) volume fractions (Extended Data Fig. 5a). We determine the residence time as the half-life of bubbles in the cap, Ψ(t ) = 0.5Ψ (see ref. 19).
where Q = q/U is a dimensionless sourcing term. We solve this equation to find the volume fraction of bubbles that can accumulate in the convecting layer, Ψ . The equation is nonlinear because of the hindrance function, and can admit more than one root. The physically meaningful solution is plotted in Extended Data Fig. 5b, and shows that the accumulated MVP volume fraction increases monotonously with the influx of MVP from the mush. Interestingly, equation (3) does not admit a real solution for injection rates that are greater than 15% of the Stokes final velocity of a 2-mm-diameter bubble in an infinite pool of melt with a dynamic viscosity of 106 Pa s (Extended Data Fig. 5b). Because these injection rates are quite modest, we expect that accumulation of bubbles up to a few tens of per cent in crystal-poor layers in magma chambers is possible. At higher injection rates, accumulation is still possible and likely to occur. However, the lack of a steady solution to our simple convecting suspension model implies that the multiphase dynamics will probably depart from that of a convecting suspension. We hypothesize that, as the volume fraction of bubbles in the crystal-poor cap increases, more complex processes may arise and lead, for example, to massive Rayleigh–Taylor overturns47.
We explain the accumulation of MVP in crystal-poor horizons of magma reservoirs by the formation of continuous MVP fingers in crystal-rich environments21, 40 and their break-up at the crystallinity transition between crystal-rich and crystal-poor magmas. This break-up of MVP fingers results in a significant change in the viscous dissipation regime.
We investigate this scenario numerically using a rather simplified geometry. We model the complex geometry of the crystal mush at the pore scale as a capillary tube that opens in a crystal-free/solid-free environment (Extended Data Fig. 6a, b). This is a simple proxy for the more realistic mush–cap transition, but it captures its essential ingredients: the dynamics of two immiscible fluids through a change in spatial confinement, where the low viscosity fluid is non-wetting and buoyant. We justify this approximation with the finding21 that the transport of immiscible fluids in a porous medium becomes mostly similar to an annular flow once the percolating pathway for the non-wetting fluid is reached. In our numerical calculations, a constant influx of MVP and a fixed pressure for the melt are set at the bottom boundary (inlet), while the top boundary (outlet) absorbs the outfluxing MVP and maintains a fixed pressure for the melt. The sides are periodic boundaries.
The competition between viscous, buoyancy, capillary and inertial forces controls both MVP transport and the breaking of continuous MVP fingering at the crystalline transition between crystal-rich and crystal-poor environments (bubble pinch-off frequency and volume49, 50, 51). Because this balance operates at the pore scale, we resort to pore-scale multiphase flow calculations to study the formation and destruction of fingering pathways in a heterogeneous medium. The force balance can be described with three dimensionless numbers, the Archimedes (Ar), Bond (Bo) and Reynolds (Re) numbers. Ar, Bo and Re represent, respectively, the ratio between buoyancy and viscous forces (equation (4)), the ratio between buoyancy and capillary forces (equation (5)) and the ratio between inertia and viscous forces (equation (6)):
where Δρ = ρ − ρ (ρ and ρ are the densities of MVP and melt), g is the acceleration due to gravity, μ the dynamic viscosity of the melt, D the bubble diameter and u the MVP average pore velocity. A rough estimate of these dimensionless numbers in shallow and highly evolved magmatic systems leads to Ar ≪ 1, Re ≪ 1 and Bo ≈ 0.1–1, we obtain a Bo of the order of approximately 0.1 and we force Re and Ar to be lower than unity. Therefore, our results can serve as good first-order estimates for MVP accumulation in crystal-poor environments. The numerical method described in the ‘Lattice Boltzmann for two-phase fluid flows’ section limits us to relatively small viscosity contrasts compared with those expected in magmatic systems. Once pathways of MVP are established in the mush, the melt plays a passive role and does not affect the ascent of the MVP. The same is not true for the suspension, where the viscosity of the melt controls the rate of energy dissipation; as such, we expect accumulation to become more efficient as the viscosity contrast between the wetting and the less viscous non-wetting fluid increases. We decided to use our numerical model to test whether bubbles are likely to accumulate under less optimal conditions, that is, when the viscosity contrast is 1/20 ≤ λ ≤ 1. We found that bubbles accumulate in the crystal-poor region even when the two fluids share the same viscosity (λ = 1), and that the accumulation potential increases as the viscosity contrast becomes more pronounced (Fig. 3).
The lattice Boltzmann method (LBM) solves a discretized version of the continuum Boltzmann equation52, 53. Based on statistical mechanics, the LBM focuses on the mechanical interaction of an ensemble average distribution of particles f (x,t), and retrieves mass and momentum conservation (Navier–Stokes) equations from the statistical moment of the Boltzmann equation.
The LBM has been extended to multicomponent (MC) immiscible fluid flows. Among others, the MC Shan–Chen (SC) model38, 54 is often applied because of: (1) its straightforward implementation; and (2) the numerical stability of the algorithm in complex geometries such as porous media. In this work, we use the SC model extended by ref. 39, which allows us to model immiscible fluids characterized by notable viscosity contrast. These improvements result from an explicit formulation of the forcing term acting on the particle distribution functions and the use of a multi-relaxation-time (MRT) collision procedure. Below, we describe the improved algorithm briefly; for more details, see refs 39 and 55.
The explicit evolution rule for the particle distribution function f α(x,t) with a single-relaxation-time (SRT) collision operator, Ω α, can be written as:
where τ is the relaxation time for fluid A and B (α = A,B) and relates to the fluids viscosity; f eq,α is the equilibrium distribution function; and f F,α is the explicit forcing term56.
The left-hand side of equation (7) is generally referred to as the streaming of f α values from the lattice node x to one of its neighbours x+e ; the right-hand side (the collision operator Ω α) describes the exchange of momentum between the colliding f values. In equation (7), e are a set of velocity vectors connecting nearest neighbour nodes (the spatial discretization of the lattice). Here we use the D3Q19 lattice—a three-dimensional lattice in which each node is connected to 19 neighbours. Lattice velocities e and weights w for a D3Q19 lattice can be found in ref. 57. The equilibrium distribution function and the explicit forcing term in equation (7) read respectively:
Here, c is the lattice speed of sound and ueq is the fluid mixture velocity defined as:
where ω = 1/τ and the statistical moments ρ (density) and ρ u (momentum) are calculated respectively as:
The forcing vectors F contain several contributions, notably cohesion (particle–particle), adhesion (particle–wall) and bulk (for example, gravity and buoyancy) forces. The cohesion forces, responsible for the physical separation between immiscible components, are calculated as:
where α and β are the two complementary phases and Gc is a free parameter that is used to tune the interfacial tension between the two fluids. The magnitude of the repulsive force applied by fluid B on fluid A at the node x (and vice versa) depends on the density gradient of fluid B (for example, ). The evaluation of is critical for the stability of the calculations. High-density gradients (thin fluid–fluid interfaces) require an extended neighbourhood to reach the required accuracy58. However, a better evaluation of density gradients comes at the price of an increase in computational time (especially in three dimensions). See refs 55, 59 for a detailed description of how to include adhesive and bulk forces. Here, in order to keep the numerical performance acceptable, we calculate the density gradients using the nearest neighbours only.
In order to improve the stability and accuracy of the SRT SC algorithm described above, we use an MRT collision operator, Ω MRT. Then, the linear collision operator is re-cast into the space of velocity moments m = M × f = (m , m , m … m , m ) (where M is the transformation matrix57); next, the relaxation parameter of each moment is adjusted individually to improve numerical stability. Ω MRT can be written as:
In this equation, Sα are diagonal matrices where the 19 diagonal components represent the relaxation parameter for each moments of f α. As suggested in ref. 57, we use:
The 19 components of the vectors mα, meq,α and meF,α can be calculated respectively as:
The stability of the algorithm depends mainly on the choice of repulsion constant (Gc) and its correspondent value at solid wall nodes (Gwall, used to introduce wetting forces). Here we want to deal with a highly non-wetting MVP phase. The non-wetting behaviour of MVP affects its dynamics both in the porous medium (higher capillary entry pressures) and at the transition between crystal-rich and crystal-poor environments (pinch-off dynamics).
In order to validate the MRT SC multicomponent algorithm that we use to model the capillary finger formation and the pinch-off dynamics (Figs. 2 and 3), we test our model with two benchmarks. The first test is an annular Poiseuille flow, where the non-wetting fluid A is located in the centre of the pipe such that r < R , and the wetting fluid B is placed in the outer ring such that R ≤ r ≤ R (where R is the radius of the pipe flow). Both fluids are accelerated by the same bulk force Fb. For the case of the two-phase Poiseuille profile problem, an analytical solution exists:
where ν is the kinematic viscosity of either fluid. In Extended Data Fig. 7a–c, we compare the analytical and numerical solutions for three different viscosity ratios (λ = 1/5, 1/10, or 1/20).
The second validation test is a three-dimensional calculation of the equilibrium shape of a drop of fluid A embedded in fluid B and in contact with a flat solid surface. The goal of this validation is to reproduce the correct equilibrium (static) contact angle between the fluids and solid phases for different wetting properties (Extended Data Fig. 7d–f).