Disequilibrium response to tapping crustal magma reveals storage conditions

Despite their importance, estimates of the pressure, temperature, saturation state, geometry and location of magmatic storage regions vary widely for even the most-studied individual volcanic systems. Geothermobarometers using mineral–mineral chemistry or phase equilibria have been applied to volcanic materials to constrain magmatic origins, but these methods have large uncertainties (about 50–200 MPa or about 2.5–10.0 km) arising from sparse laboratory constraints, limits on analytical precision, assumptions of local equilibrium and interdependence between temperature and pressure^2,3. Furthermore, magmatic reservoirs are challenging to image with geophysical methods owing to limitations in resolution and poorly constrained relationships between lithology and geophysical signals, resulting in typical uncertainties on magma depth of approximately 0.5–10.0 km (ref. ⁴). As a result, the scientific community lacks consensus on even the fundamentals of the spatial distribution of melt in the crust⁵.

Deep drilling in hydrothermal fields offers the unique potential for placing tight constraints on the location, temperature, pressure and chemistry of melt stored in natural volcanic systems. Hydrothermal drilling has occasionally intersected magma: dacite in the Puna geothermal field (Hawaii), trachyte at Menengai (Kenya) and rhyolite at Krafla (Iceland). At Krafla, the KJ-39 and Iceland Deep Drilling Project-1 (IDDP-1) boreholes directly intersected magma at roughly 2,500 and at 2,104 m depth, respectively^6,7, which was not anticipated on the basis of coarse geophysical imaging before drilling. Retrospectively, a magma body was recognized during reanalysis of magnetotelluric data⁸ and from seismic imaging^9,10, which have still been unable to resolve questions about its lateral and vertical extent.

Silicic glass fragments were recovered from the IDDP-1 borehole, which were quenched through interaction with drilling fluids^7,11. The precisely known recovery depth, and temporally and spatially constrained ascent and quench, makes them ideally suited to resolve the unknowns of magmatic storage, as the melt has not been subject to the complex ascent processes that afflict the products of volcanic eruptions. Indeed, the glass chemistry has sparked discussion about the origins of the magma from (1) partially molten, hydrothermally altered basaltic crust¹² or (2) mantle-derived basalts evolved by means of fractional crystallization¹³ and about the degree of crustal assimilation^1,12,14,15. Despite the direct sampling of the glass, established equilibrium-based methods to determine pressure and temperature have yielded wide constraints for these magmas (Fig. 1a): various geothermometers yield different estimates and uncertainties from two-pyroxene equilibration between 920 °C and 940 °C (refs. ^1,16) and 890 °C and 910 °C (refs. ^12,17), whereas the phase assemblage places the weak constraint between 800 °C and 950 °C (ref. ¹). Pressure estimates from projection on the haplogranitic ternary suggest <50–100 MPa, whereas Rhyolite-MELTS yields 44 ± 11 to 47 ± 32 MPa depending on oxidation state, all of which require an assumption of volatile content and fail to model this system in that they produce quartz + plagioclase ± orthopyroxene, when the observed phases are plagioclase + two clinopyroxenes¹⁸. The pressure of magma is often, instead, determined by the volatile (H₂O and CO₂) concentrations of crystal-hosted melt inclusions, which are sensitive to pressure and temperature and less susceptible to decompression-induced changes than the residual melt (quenched to glass). Melt inclusions are arguably inappropriate to investigate the pressure and saturation state of the IDDP-1 magma for three reasons: (1) the few crystals present in the IDDP-1 chips frequently show dissolution (melt embayments and rounding)¹, so a meaningful, reliable population has not been isolated; (2) the pressure/depth of melt inclusions are usually determined on the basis of the assumption of volatile saturation and are, by their very nature, contained within individual crystals and thus separated from mineral pairs, which could provide an independent pressure determination; and (3) are produced preferentially during disequilibrium crystal growth and so systematically oversample non-equilibrated conditions during transport or perturbations from background³ and therefore are poorly suited to investigate equilibrium storage. Instead, we are left with the measured volatile contents of the glass, which correspond to saturation pressures between about 35 MPa and 45 MPa (Fig. 1a and Extended Data Fig. 1), below the lithostatic pressure (about 50–57 MPa) and above the hydrostatic pressure of the well (about 16 MPa)¹. These measurements have been interpreted to suggest that the magma is either (1) stored and degassed to equilibrium at a pressure less than lithostatic owing to interaction with the hydrothermal system^1,12 or (2) originally undersaturated^14,19.

a, Temperature (red) with depth in the borehole⁶ versus geothermometric constraints from augite–pigeonite and clinopyroxene–orthopyroxene–plagioclase–magnetite–ilmenite and Rhyolite-MELTS¹ plotted at the recovery depth, with a secondary axis showing lithostatic pressure at the corresponding depth and saturation pressure (blue, from VolatileCalc^1,19; error bars indicate 1σ uncertainty). b,c, Although a trans-crustal mush arrangement of magma storage (b) is typically favoured⁵, classic views of a basalt-underplated continuous, large rhyolite source (c) have also been proposed at Krafla^12,13,41 and are not distinguishable based on present geophysical observations^8,9,10. d,e, Backscattered electron images of IDDP-1 glass chips that show the variability in vesicularity. Scale bars, 1 km (b); 500 μm (d,e).

Although the drilling fluids rapidly cooled the magma, it was still subject to decompression following intersection by the drill string¹⁴, resulting in remobilization; the magma flowed 8 m up the well in 9 min (ref. ²⁰). So we ask, can the IDDP-1 glass be used as a direct record of storage conditions? Here we expand a bubble growth model for H₂O (ref. ²¹) to include CO₂, coupled to a model for water species interconversion²². As well as the dynamics of bubble growth, a direct output of this model is the residual volatile content in the melt/glass between bubbles from which we constrain the magma storage conditions of pressure, temperature and volatile composition against the measured glass chips.

The glass fragments contain total water of 1.3–2.0 wt%, with an average of 1.8 wt%, consistent across several studies^{1,11,13,19,23,24}, and 50–200 ppm CO₂ (refs. ^1,19). Vesicularity in the quenched glasses is low, with most chips having <6 vol%, although occasionally up to about 15 vol% (ref. ¹⁴). Bubble sizes are 1.5–75.0 μm, with an increase in bubble size with drilling time¹⁴. Bubble number densities range between 10^11.7 and 10¹⁵ m⁻³, which are inferred to have nucleated during drilling-induced decompression at rates of 10⁵–10⁷ Pa s⁻¹ (refs. ^14,25). OH/H₂O_m ratios are between 1.68 ± 0.45 and 2.19 ± 0.37, increasing over time¹⁹, and are slightly lower in more vesicular, 2.07 ± 0.20, than in less vesicular, 2.13 ± 0.38, fragments¹.

We explore which pressure–temperature (P–T) paths best reproduce the volatile chemistry and vesicularity of the IDDP-1 glass. Although the cooling from the drilling fluid should be rapid, we use a thermal model to seek paths consistent with measurements of the natural glass; geospeedometry through differential scanning calorimetry indicates that the IDDP-1 magma cooled through the glass transition regime (T_g) at 7–80 °C min⁻¹, with a T_g of about 480 °C (ref. ²⁶). We begin with one-dimensional cooling in a planar geometry from an initial temperature of 900 °C by means of conduction and forced convection at the surface with steam (400 °C and 16 MPa), in which drilling fluid is in direct contact with the melt. At a planar interface, cooling rates far exceed those measured in the glass. However, at larger length scales, thermal diffusion is inefficient, such that we can only reproduce the cooling rates in a narrow region between 0.3 mm and 1.3 mm from the interface (Extended Data Fig. 2). To produce a nearly homogeneous glass, similar cooling rates must be sustained across a large portion of the sampled material. This is consistent with fracturing or fragmentation of the material during cooling (Fig. 1d,e) to enhance the surface area and enable cooling from several directions.

If we assume that magma fragments (by any mechanism) during the onset or early progression of cooling, we should model the thermal evolution instead using a spherical geometry. Under these conditions, we reproduce the measured cooling rates at T_g in the interior of melt fragments of radius 9–20 mm (Extended Data Fig. 3), a spacing common in perlitic glass^27,28. This results in a timescale for cooling from storage to T_g (at 480 °C) in about 4 min, with a decelerating cooling rate.

Now we turn to the pressure path experienced by the magma. During drilling, the magma should experience rapid decompression, reflected by high bubble number densities^14,25. However, the volatile composition and vesicularity of the glass changed only modestly over 9 h of drilling¹⁴, suggesting that decompression of the magma was also progressive. We suggest that the decompression could be offset at increasing depth by the pressure drop generated by viscous resistance to flow up the borehole. We estimate the relevant length scale using the Hagen–Poiseuille equation:

using a borehole radius, R, of 0.15 m (12.5″ diameter), a melt viscosity, μ, of 3.2 × 10⁵ Pa s for magma at 900 °C and 1.8 wt% H₂O, a vesicularity, ϕ, of 1–5%, bubble number density of 10¹¹–10¹³ m⁻³ (ref. ²⁵) and a rate of bubble growth of 5 μm s⁻¹, resulting in a ∂ϕ/∂t ≈ 1–10% s⁻¹. We find a length scale, H, of about 1.5–4.0 m over which the pressure increases back to the investigated storage conditions (35–55 MPa). This length scale would become shorter directly proportional to an increase in viscosity. This has two implications for magma response: (1) the growth of bubbles in a partially confined space acts to buffer the pressure of magma reservoirs, such that perturbations from drilling probably only have a local (centimetres to metres) region of influence, and consequently (2) most magma may continue to experience near-storage pressures until the onset of cooling.

The model outcomes above suggest that magma conditions are modulated by cooling and decompression fronts, whose rates of propagation will be closely linked to, and driven by, fragmentation. The shear rate for magma flowing in the conduit is about 1–3 s⁻¹, insufficient to produce large-scale shear-driven fragmentation at this viscosity²⁹. Similarly, rapid decompression can result in explosive fragmentation but requires large pressure drops for low-porosity magmas³⁰. Instead, we suggest that non-explosive fragmentation from thermal shocking may be the prevalent fragmentation mechanism, similar to fuel–coolant reactions and the production of hyaloclastite and other subaqueous products³¹. We use the planar thermal simulations above and a thermal expansivity of about 6 × 10⁻⁶ °C⁻¹ to find a region of high strain rate in a propagating front that moves at about 30 μm s⁻¹. Any crack opening would depend on the tensile strength of the material and be assisted by fluid injection and volume expansion within the cracks. This crack opening results in an ‘unzipping’ of the magma, leading to progressive decompression.

Because we do not simulate the exact fragmentation dynamics, we consider a range of possible decompression rates from the initial volatile saturation pressure (35–55 MPa) evolving towards hydrostatic (16 MPa) over 10 s, 30 s, 1 min, 5 min, 30 min and 1 h, corresponding to decompression rates 5.3–2,900 kPa s⁻¹ (Extended Data Table 1). The initial conditions are specified by choosing a starting volatile saturation pressure and molar ratio of water and CO₂ in the coexisting vapour phase, which together set the total water and CO₂ volatile contents. With a vesiculation model, we cannot distinguish between the two competing hypotheses of an undersaturated magma (which produces a first volatile phase when the system pressure drops to the saturation pressure) and magma stored saturated, but at less-than-lithostatic pressure with a negligibly small proportion of coexisting vapour. Figure 2 shows results with initially 1.8 wt% H₂O_t and 111 ppm CO₂ at 45 MPa and 900 °C and a OH/H₂O_m ratio of 4.92 (at equilibrium).

a–c, Vesicularity along different decompression timescales (a) and cooling timescales along rapid (30 s) (b) and slow (5 min) (c) decompression paths. The shaded regions represent the typical ranges of vesicularity measured in the IDDP-1 glass chips, with darker shading representing greater abundance among measured chips. Line colour represents the timescales for decompression (blue) and quenching offset (red).

Source data

Complete decompression from 45 MPa and 900 °C and 35 MPa and 920 °C to 16 MPa produces final vesicularities of 26.2% and 25.6% and water concentrations of 1.34 wt% and 1.32 wt%, respectively (Fig. 2), both more vesicular and drier than the measured glass. At magmatic temperatures, bubble growth is rapid and only at the highest decompression rates (<5 min) is there a substantial kinetic delay in which bubble growth continues for 20–30 s after the end of decompression. This suggests that quenching of the glass must begin before the magma is completely decompressed.

We simulate decompression interrupted by cooling at some time offset, ∆t_quench. When the timescales for decompression are short compared with cooling (Fig. 2b) or when cooling occurs late, the magma nearly completely degases, but cooling produces resorption of bubbles; although decompression reduces solubility, cooling does the reverse and the decompressed magma can resorb up to approximately 10 vol% vesicles during cooling. However, to produce final vesicularities within the measured range, cooling must be early or rapid compared with decompression (Fig. 2c), in which case the final vesicularity is sensitive to the relative timing of cooling onset and the glass locks in a low vesicularity and high water content.

To find the likely storage conditions and decompression paths, we restrict the timing of cooling onset to produce the desired final vesicularity of about 5%. Because water is the volumetrically dominant species, the final water concentration of the glass and the vesicularity are closely related (Fig. 3a,b,e,f,i,j). Together they help constrain the initial water content (assuming minimal loss to the surroundings), but there remain a large number of possible P–T paths that are suitable.

a–l, Vesicularity (a,e,i), total water (b,f,j), OH/H₂O_m ratio (c,g,k) and CO₂ concentrations (d,h,l) in the melt/glass along different decompression paths from 35 MPa (a–d), 45 MPa (e–h) and 55 MPa (i–l) initial saturation pressure with syn-decompression quenching whose onset is chosen to reproduce the total vesicularity.

Source data

Similarly, the OH/H₂O_m ratio is primarily dependent on the cooling path. Accordingly, all simulations produce a narrow range of values between 1.75 and 1.76 (Fig. 3c,g,k), which corresponds to an equilibrium temperature of roughly 500 °C. Because the cooling rate close to the glass transition was constrained by the measured glass transition temperatures and the thermal modelling, the OH/H₂O_m ratio has little further discriminatory power but is consistent with the geospeedometry²⁶.

Conversely, owing to the large differences in water and CO₂ diffusivity (one to two orders of magnitude), the final H₂O/CO₂ ratio of the glass is highly path-dependent. On initial decompression, both H₂O_m and CO₂ diffuse into the bubbles, but the slower transport of CO₂ increases the H₂O/CO₂ ratio in the bubble, which in turn changes the solubility of both species in the melt. This has the overall effect of delaying rapid bubble growth in the first approximately 10 MPa of decompression compared with the pure water case. On cooling, we see rapid water resorption but no marked change in the CO₂ dissolved in the melt, both because resorption is kinetically limited and owing to the changing H₂O/CO₂ ratio in the gas, which may allow CO₂ to continue moving into the bubble. As a result, the CO₂ in the final glass is more sensitive to the initial vesiculation than to resorption.

Because the CO₂ has preferentially been lost to the bubbles, previous estimates using a mixed solubility model and only the glass chemistry¹ would have resulted in an underestimate of storage saturation pressure. Although we can produce reasonable glass chemistry along rapid decompression and early cooling from 45 MPa, decompression from lithostatic pressure (55 MPa) using a modestly higher initial CO₂ concentration of 170 ppm and 1.8 wt% water better encompasses the range of measured volatile concentrations in the glass chips (Fig. 4). We cannot produce any glasses with more than about 50 ppm CO₂ from 35 MPa. Starting from a cooler (about 850 °C) initial temperature would have only a modest effect on starting CO₂ (121 versus 108 ppm at 45 MPa and 1.77 wt% H₂O) and so this method cannot place tighter constraints on the starting temperature compared with existing geothermometers and experimental studies³².

To preserve sufficient CO₂ in the melt, we require that the characteristic timescale for decompression must be ≲5–10 min, similar to the characteristic timescale for cooling, and cooling onset must be approximately synchronous with or slightly trailing (≲3–4 min) decompression (Figs. 4 and 5). If decompression occurs over the 1.5–4.0 m characteristic length scale found above, we find propagation rates of 15–80 cm s⁻¹, resulting in a fragmentation level that is nearly stationary or rises up the borehole with the flow of magma. The magma cannot have undergone slow decompression for tens of minutes to hours during approach of the drill string, suggesting that the roof rock is capable of isolating the magma from the surrounding hydrothermal system, at least on short timescales.

Calculated paths (blue lines) of average H₂O and CO₂ in the melt along decompression from magma storage (sketched in red) at lithostatic pressure to the final glass (blue markers), which span the field of measured glasses^1,19 (grey markers; error bars indicate 1σ uncertainty). Black lines show the saturation pressure⁴² at 900 °C. The inset highlights the trajectory of the volatile saturation path during decompression at different rates and synchronous resorption.

Source data

Sketch of processes during magma perturbation by drilling.

In our simulations, diffusion is rapid enough to almost completely erase the gradient in both total water content and speciation except in the very last stages of cooling while CO₂ either minimally resorbs or even continues to diffuse into the bubble (Extended Data Fig. 4). This suggests that resorption may be both pervasive^33,34,35 and difficult to detect in the final glass, in which a notable gradient towards the bubble is conventionally used as the primary indicator of resorption^35,36,37. This also has implications for our interpretation of the origin of the magmatic volatiles; the δD, δ¹⁸O and total water contents have been interpreted to evidence mixing between mid-ocean ridge basalts and meteoric water, along with crustal assimilation by assuming that the isotopic signatures are not substantially altered by exsolution owing to the low final vesicularity¹⁴. However, vesiculation followed by rapid rate-limited resorption contributes to fractionation, in which the heavier isotopes are enriched in the bubble during exsolution, and cannot diffuse as rapidly during resorption^38,39,40.

Our study highlights the rapid, disequilibrium changes that overprint equilibrium conditions and the importance of carefully accounting for the magma physico-chemical response during quenching; we demonstrate that glass is not a direct representation of storage, even under rapid cooling. During drilling, magma at Krafla responded to pressure and temperature perturbations over extremely short timescales (seconds to minutes) and small length scales (centimetres to metres), closely tied to fragmentation driven by thermal shock. By correcting for this rapid response using our multiparametric model, we produce the first, to our knowledge, robust constraint of magmatic conditions, confirming storage of volatile-saturated magma at lithostatic pressure. The ability to invert for magmatic conditions using drilling samples during rapid cooling and quenching at depth, however, raises the challenge inherent in reconstructing storage conditions from volcanic products that have traversed kilometres of crust over hours to weeks. Beyond the validation of magmatic conditions, the thermo-physico-chemico-rheological-kinetic model developed here provides a first rigorous tool to engineer magma response and enable safe access through drilling by forward modelling to find the optimal drilling conditions (for example, borehole geometry, penetration rate, fluid composition) to inhibit magma ascent up the borehole.