Evolving magma temperature and volatile contents over the 2008-2018 summit eruption of Kīlauea Volcano.

Magma rheology and volatile contents exert primary and highly nonlinear controls on volcanic activity. Subtle changes in these magma properties can modulate eruption style and hazards, making in situ inference of their temporal evolution vital for volcano monitoring. Here, we study thousands of impulsive magma oscillations within the shallow conduit and lava lake of Kīlauea Volcano, Hawai'i, USA, over the 2008-2018 summit eruptive sequence, encoded by "very-long-period" seismic events and ground deformation. Inversion of these data with a petrologically informed model of magma dynamics reveals significant variation in temperature and highly disequilibrium volatile contents over days to years, within a transport network that evolved over the eruption. Our results suggest a framework for inferring subsurface magma dynamics associated with prolonged eruptions in near real time that synthesizes petrologic and geophysical volcano monitoring approaches.

INTRODUCTION

Kīlauea volcano, Hawai‘i, USA, is one of the most active, best-monitored, and best-studied volcanoes on Earth (), serving as a focal point for volcanologic research (). However, resolving in situ variation in subsurface magma dynamics remains a challenge at Kīlauea and volcanoes globally (). The 2008–2018 Kīlauea summit eruption represents an opportunity to address this knowledge gap. The eruption involved a persistent lava lake in the Halema‘uma‘u summit vent and multiple subsurface magma intrusions and East Rift Zone eruptions, ending with a spectacular caldera collapse sequence representing the highest historical sustained eruption rate at Kīlauea (–). Previous studies suggested the main Kīlauea shallow summit magma plumbing system during this time consisted of the 1- to 2-km-deep Halema‘uma‘u reservoir and the 3- to 5-km-deep South Caldera reservoir (Fig. 1) (, ). The Halema‘uma‘u reservoir and overlying lava lake were continuously connected () by a ∼10-m-wide conduit (). Magma passed through the summit en route to the East Rift Zone, although the nature of hydraulic connections between the summit reservoirs, rift zone, and deeper magma sources is not well known (, ).

Fig. 1.

Kīlauea map and magma dynamics model.

(A) Map including the Halema‘uma‘u vent, inferred shallow magma storage zones, GNSS stations, and seismometers used in the VLP catalog (). (B) Typical lava lake activity on 13 February 2017 U.S. Geological Survey (USGS). (C) Seismic waveform from a VLP conduit-reservoir resonance event along with a model solution for reference fixed parameter inversion results forced with a Gaussian pressure perturbation (fig. S1). UTC, universal time coordinated. (D) Conduit-reservoir resonance model with approximate 2018 magma system geometry; black arrows illustrate vertical sloshing of the stratified magma column. ASL, above sea level. (E) Magmastatic depth profiles from piecewise linear total (dissolved plus exsolved) volatile mass fractions at a uniform temperature of 1200°C.

A wide range of data, interpreted using physical and chemical models, inform this picture of magma dynamics. Transport geometry is constrained primarily through inversion of seismic and geodetic data (, , ). Continuous gravity data are only available over limited time segments but constrain the density of magma in the lava lake and suggest temporal variation of up to 1500 kg/m3 (). Analysis of erupted products provides limited temporal and spatial resolution but suggests that Halema‘uma‘u magma consists of near-liquidus (1150° to 1300°C) crystal-poor basalt outgassed in CO2 with respect to the primary mantle magma (, ). Subsurface magma volatile contents are also indirectly informed by continuous gas emissions (, , ). These analyses suggest substantial disequilibrium outgassing or mechanical decoupling of gas bubbles from melt because of continuous convecting and outgassing (). However, geochemical and geophysical data are rarely combined in a quantitative manner.

Very-long-period (VLP) seismicity, with energy concentrated at periods above 2 s, has the potential to help unify these diverse constraints. VLP seismicity is prevalent at many volcanoes and often inferred to represent transient magma flow (), thus directly probing magma properties and transport geometry in ways not readily obtainable by other geophysical analyses. VLP signals are part of a spectrum of oscillatory motions that can result from impulsive or continuous forcing of magma transport structures (, ), but the VLP band is advantageous because it is less sensitive to path distortions from heterogeneous earth structure than shorter period signals.

Multiple resonant modes have been identified at Kīlauea, but the dominant VLP signal is from “conduit-reservoir” resonance, in which stratified magma in the conduit and lake sloshes in and out of the underlying reservoir (Fig. 1) (, , ). This resonance occurs sometimes as continuous tremor but most often as discrete minutes-long events triggered both from the lake surface (such as via rockfalls from the crater walls) and from depth (, ). Oscillation-restoring forces are from gravity and magma reservoir elasticity, while damping is from viscous drag on the conduit walls. Resonant period is primarily sensitive to conduit length and bulk magma density/density stratification (). Decay rate is quantified by quality factor (the ratio of energy stored per cycle over energy lost per cycle) and is primarily sensitive to conduit radius and apparent magma viscosity. In the shallow Halema‘uma‘u magma system, where melt composition does not vary much in time or space and where crystal contents are low (, , ), magma density is primarily controlled by porosity, and magma viscosity is primarily controlled by porosity and temperature. In chemical equilibrium, gas mass fraction (hence porosity) depends on total volatile mass fraction and pressure-dependent solubility of dominant volatile species (H2O, CO2, and sulfur) ().

VLP seismicity at Kīlauea thus reflects evolving magma thermal and chemical state as well as transport structures. Over the 2008–2018 Kīlauea eruption, thousands of conduit-reservoir resonance events provide an unprecedented record of time-evolving subsurface magma transport.

Approach: Inferring magma properties from geophysical data

Figure 2 outlines our workflow. We first conduct kinematic elastic inversions between 2008 and 2018 of continuous Global Navigation Satellite System (GNSS) ground deformation data (figs. S3 and S4) () for shallow magma reservoir pressure histories. In particular, Halema‘uma‘u reservoir pressure constrains magma column density in the overlying summit lava lake and conduit. Summit deformation at Kīlauea is complex: To resolve Halema‘uma‘u reservoir pressure, we build on constraints from previous geodetic studies (, , ) and include three known deformation sources ().

Fig. 2.

Inversion approach.

(A) Simplified flowchart of methods and data input/output. Additional constraints on GNSS inversions are from previous geodetic studies (, , , ). Additional constraints on VLP magma resonance inversions are from previous modeling (), gravity data (), and geochemical (gas and ejecta) data (, , , ). (B to F) Conduit-reservoir resonance period and quality factor, plus conduit bottom pressure, as a function of the parameters varied to fit Kīlauea VLP seismic and geodetic data. Variations in lava lake elevation and (assumed uniform) radius are prescribed from measurements (, ). Dashed black lines indicate default values used in the other plots.

We next use a perturbation approach to model transient flow associated with conduit-reservoir magma resonance (Fig. 1) (), extending previous analyses (, ). We treat fluid properties of the multiphase magma as functions of magmastatic pressure (an approximation given slow exchange flow within the conduit/lava lake ()), temperature, and vertically stratified total volatile mass fractions (CO2 + H2O; Fig. 1 and fig. S2), neglecting crystals and assuming an average melt composition based on 2008–2010 Halema‘uma‘u samples (, , –). We use this model to invert for magma properties from Halema‘uma‘u reservoir pressure, lava lake elevation and areal extent (, ), and the resonant period and quality factor of VLP seismic events cataloged over 2008–2018 by () (Fig. 2) ().

Resolving time evolution of shallow magma properties at Kīlauea is a long-standing challenge (, , ). We focus on shorter-term changes in multiphase magma properties by assuming a fixed magma system geometry based on previous inversions (, , ). Four additional assumptions are made to facilitate unique inversions for magma properties (Supplementary Text) (): (i) Temperature is spatially uniform in the conduit and lake. This is justified because the conduit undergoes quasi-steady exchange flow/mixing (), and the lake contributes negligibly to viscous damping. (ii) Magma in the conduit/lake has a fixed total (dissolved + exsolved) H2O/CO2 mass ratio. Volatile composition could vary over time but is unconstrained in our model without additional data, so we fix volatile ratios based on erupted products and gas emissions (, , ). (iii) Total volatile mass fraction varies linearly with depth (Fig. 1) subject to stable stratification, which should be approximately valid for the largely quiescent magma column. (iv) Total volatile mass fraction at the lake surface is constant. While there is known to be some variation in porosity near the lake surface from continuous gravity data (), these data are not available over most of the time span. In addition, our model exhibits minimal sensitivity to density stratification within the lake; it is primarily sensitive to average density (which controls the magmastatic pressure load of the lake on the conduit).

We test different fixed parameter combinations and conduct an a posteriori assessment of these assumptions. The magma properties we invert for are (i) magma temperature, (ii) conduit average total volatile mass fraction Xavg, and (iii) total volatile mass fraction stratification (difference between conduit top and conduit bottom) ΔX. We note that while the magma temperature parameter is applied to the whole magma column, the model is primarily sensitive to conduit temperature. We also note that because of the trade-offs between volatile contents at the bottom and top of the lava lake, ΔX should be considered to represent a general volatile stratification over the whole magma column (conduit and lava lake).

RESULTS

For our reference fixed parameters, Fig. 3 shows the timeline of GNSS inversion results and VLP magma resonance inversion results, along with other data. Shaded regions in Fig. 3 show the envelope of inversion results obtained by varying individual fixed parameters over feasible ranges, as detailed in the Supplementary Text (fig. S5). Evolution of magma system geometry, which is not considered in our inversions, is more likely to affect trends in inversion results over long (year or more) time scales. In particular, inversion results with the reference fixed parameters are likely not reliable in 2009 to early 2010 and mid-2011 (Discussion). On short time scales, noise in input data likely contributes to scatter and outliers in the inversion results. We thus focus most analysis on temporally averaged values and, in particular, on the relative variability in these values over time scales of a year or less rather than their absolute value at a given time. Figure 4 shows amplitude spectra, coherence, and phase lags between data sets with 95% significance thresholds (Supplementary Text). Additional analyses are shown in figs. S6 to S8.

Fig. 3.

Time-series data and inversion results.

Inverted relative changes in magma properties are from our reference fixed parameters (Fig. 1 and table S1). Dots represent individual VLP seismic events, bold lines are 30-day moving averages, while vertical green lines are East Rift Zone eruptions (solid), summit intrusions (dashed), and slow-slip events (dotted) (). (A) VLP seismic event resonance period and quality factor (). (B) Lava lake elevation and mean radius (, ) (C) GNSS inverted reservoir pressure changes, set to zero at the 7 March 2011 lava lake draining. Shaded areas indicate possible variation with different South Caldera reservoir geometries tested (Supplementary Text). (D) Inverted conduit magma temperature, with MgO thermometry for comparison (, ). The shaded area indicates possible variation with all fixed model parameter values tested (Supplementary Text). (E and F) Inverted conduit total volatile contents, with 30-day moving average SO2 emissions for comparison (, ) and possible variation shown in shaded areas. Values from 2009 to early 2010 are unreliable because of exact solutions not being obtainable with the fixed parameters chosen.

Fig. 4.

Wavelet amplitude spectra and coherence.

(A) Amplitude spectra of resonance properties (), lava lake elevation (, ), SO2 emissions (, ), GNSS inverted Halema‘uma‘u (HMMR) and South Caldera (SCR) reservoir pressures, and VLP magma resonance inverted magma properties. (B) Magnitude squared coherence colored by phase lag. The gray area is beneath the 95% significance threshold. Positive phase lags indicate that the second variable trails the first. Data before December 2011 were excluded from this analysis.

As expected for an open-vent magma system, Halema‘uma‘u reservoir pressure is well correlated with lava lake elevation over time scales from days to about a year (Figs. 3 and 4) (, ). Strong coherence between Halema‘uma‘u and South Caldera reservoir pressures over time scales of days to months (Fig. 4 and fig. S6) suggests that magma is often transferred between the reservoirs, although the anticorrelation implies hydraulic disequilibrium. This could indicate an intermittent connection, consistent with the unsteady connectivity inferred during hours- to days-long “deflation-inflation” events (, , ). We are not aware of any other settings where a consistent anticorrelation is observed between different magma reservoirs at the same volcano, although intermittent hydraulic connections have been inferred between Kīlauea and Mauna Loa (), as well as at other volcanoes such as Soufriére Hills () and Etna ().

Different fixed parameters affect the absolute value of inverted magma temperature, but the pattern of relative temporal variation is robust, and the magnitude of such changes varies by less than ∼20°C (Fig. 3 and fig. S5). Inverted temperature is primarily sensitive to conduit radius; decreasing radius by 10 m (to 5 m) uniformly increases temperatures by ∼60°C, while increasing radius by 10 m (to 25 m) uniformly decreases temperatures by ∼40°C. Conduit magma temperatures span the full 1150° to 1300°C range of Halema‘uma‘u magma storage temperatures previously estimated from ejecta geothermometry (, ), although it is difficult to make a direct comparison given uncertainty in the depths and/or time scales recorded by geothermometers.

On time scales from days to months, temperature exhibits up to 100°C variation (Fig. 3), corresponding to up to an order of magnitude variation in magma viscosity (figs. S2 and S8). Temperature and resonant quality factor are strongly correlated (fig. S6), which suggests that temperature is a primary driver of variations in magma viscosity. The dominance of temperature is unexpected because porosity has previously been proposed as a likely source of variation in VLP quality factor () and is known to vary significantly as bubbles rise and accumulate (, ).

Different fixed parameters affect the inverted absolute value of Xavg by up to ∼1 weight % (wt %), but the pattern of relative temporal variation is robust and the magnitude of such changes varies by less than ∼0.4 wt % (Fig. 3 and fig. S5). Similarly, different fixed parameters affect the inverted absolute value of ΔX by up to ∼1 wt %, but the pattern of relative temporal variation is robust, and the magnitude of such changes varies by less than ∼0.2 wt % (Fig. 3 and fig. S5). Over most of the timeline Xavg is greater than the inferred primary magma volatile mass fraction of 1 to 2 wt % , a notable accumulation particularly because some of the primary CO2 may have already been lost at depth (, , ). In addition, ΔX is mostly similar to or larger than inferred primary magma volatile mass fraction. Together, these indicate substantial departures from equilibrium outgassing, with an accumulation of volatiles in the upper conduit and lava lake.

On time scales of days to months, Xavg varies by up to ∼0.6 wt %, and ΔX varies by up to ∼1 wt % (Fig. 3). That this temporal variation is similar to the inferred primary magma’s total volatile mass fraction of 1 to 2 wt % (, ) suggests strong variations in the outgassing regime (). The only volatile species with continuous emission measurements that can be compared with ΔX and Xavg is SO2. SO2 has roughly similar solubility to H2O in mafic melts () and so will approximately trade-off with H2O in our model. SO2 emissions exhibit strong variation (an order of magnitude or more) on time scales from days to years (, ). We do not observe consistently strong coherence between ΔX or Xavg and SO2 emissions (fig. S6), although several pronounced increases in either ΔX or Xavg do correspond to increases in SO2 (e.g., April 2015, January 2016, October 2016, and August 2017). Inconsistent coherence could partly reflect the high uncertainty in SO2 emission data, although we note that gas emissions from the lava lake surface will not necessarily directly correlate with the amount of volatiles accumulated in the magma column. The strong in-phase coherence between Halema‘uma‘u reservoir pressure (or lava lake elevation) and ΔX on time scales of less than 90 days (Fig. 4) suggests that volatiles build up in the upper conduit/lake as magma accumulates in the Halema‘uma‘u system rather than maintaining a steady volatile mass balance through the shallow magma column. This could reflect an increase in volatile flux (e.g., from magma recharge), but could also be caused by less efficient outgassing through the lava lake as it fills.

DISCUSSION

Halema‘uma‘u magma mass balance

Maintaining a persistent lava lake for a decade requires a remarkable thermal and mechanical balance. Relatively constant magma supply from depth is needed to drive continuous convection, but supply must be countered by sufficient outflux to prevent conditions leading to violent eruption. Ground deformation and VLP seismicity provide a quasi-continuous probe of magma properties that facilitates interrogation of the multiscale processes maintaining (and modulating) this balance within the Halema‘uma‘u reservoir during an extended eruption.

In general, magma reservoir pressure can change even without any magma input due to gas exsolution and (to a lesser extent) crystallization. However, because the low-viscosity mafic melt and open-vent structure of Halema‘uma‘u facilitates gas escape, reservoir pressurization has been inferred to reflect accumulation of melt due to changes in either influx (e.g., recharge from the South Caldera reservoir or deeper storage regions) or outflux (e.g., to the East Rift Zone) (, ). For example, the inferred causes of the May 2015 summit intrusion, the 2018 eruption, and the prevalent hours- to days-long deflation-inflation summit deformation events are months of increased magma influx (, , ), months of reduced magma outflux (), and transient restrictions of magma influx or outflux (, , ). However, the general controls on magma mass balance over days to years are unknown. The 60- and 130-day period spectral peaks in Halema‘uma‘u reservoir pressure (also apparent in temperature, ΔX, and Xavg) (Fig. 4) may indicate dominant time scales for such changes in influx-outflux (). Quasi-periodic deformation and/or eruptive activity on similar time scales has also been observed at other volcanoes (, ).

We might expect magma recharge to increase conduit temperature, although this would depend on the temperature and influx of recharging magma and also its path through the ∼4 km3 of near-liquidus magma in the Halema‘uma‘u reservoir (, ). The inferred 2011–2012 average magma supply rate of ∼109 kg/day () would permit complete exchange with the ∼1010 kg of magma in the conduit and lava lake over a week. However, if this injected magma were uniformly mixed with the magma in the reservoir (∼1013 kg assuming a density of 2500 kg/m3) at a 100°C temperature difference, the mixture temperature would only increase by ∼0.01°C/day (neglecting latent heat and outflow). Given the poor coherence between Halema‘uma‘u reservoir pressure (or lava lake elevation) and temperature (Fig. 4), we expect that melt injected into the Halema‘uma‘u reservoir generally either was not appreciably hotter than existing magma and/or was not directly routed to the conduit.

One prominent exception that could exemplify an influx of hotter melt from depth is the persistent ∼100°C increase in temperature 6 months before the March 2011 Kamoamoa fissure eruption. There was no corresponding increase in volatile mass fractions, potentially due to deeper separation and upward flux of volatiles over the preceding months of elevated volatile mass fractions. Temperature then dropped by ∼100°C in the months leading up to the eruption, which we expect relates to lava lake downwelling rather than magma influx/outflux, as discussed in the next section. Another potential example of hot melt influx is the ∼90°C increase in temperature between the May 2012 slow-slip event on Kīlauea’s south flank décollement and the October 2012 intrusion, although there was also no corresponding increase in volatile mass fractions. The temperature increase supports previous suggestions that slow-slip events are linked to magmatism (), although we do not see similar temperature increases immediately following the 2010 or 2015 slow-slip events.

It is less obvious what changes in magma properties might be expected from decreased magma outflux, so we use the 2018 eruption as a case study. The months of pressurization preceding the eruption are accompanied by a decrease in magma temperature and increase in Xavg, but these do not clearly stand out from the background variation over the preceding year (Fig. 3). The lack of clear changes in magma properties is consistent with the idea that the 2018 eruption was triggered by decreasing outflux rather than by recharge () and, by extension, suggests that outflux does not necessarily drive notable changes in shallow magma properties. The May 2014 and May 2015 intrusions were also preceded by a month of Halema‘uma‘u reservoir pressurization without other clearly associated changes in the summit magma system. The lack of clear changes in magma properties would seem to suggest they were induced by decreased magma outflux, although at least in 2015, changes in East Rift Zone lava effusion were not apparent (, , ). The June 2014 and May 2016 Pu‘u‘Ō‘ō vent openings were not preceded by notable pressurization of the shallow summit magma system, suggesting they were not primarily caused by increased melt flux from the summit but rather by processes along the rift zone.

Shallow magma dynamics

Our results illuminate shallow fluid dynamic processes underlying a persistent lava lake. Observed covariation of parameters in our inversions suggests that volatile mass fraction and temperature in the conduit and lava lake vary in ways not always directly related to Halema‘uma‘u reservoir magma influx/outflux. We infer that such variation occurs because of unsteady exchange flow between the conduit and Halema‘uma‘u reservoir (), as well as because of changing convective efficiency in the lava lake and/or surface crust dynamics (which influence the outgassing rate and efficiency of heat loss to the atmosphere and host rock) (, ).

The negative correlation on time scales of months or less between Xavg and temperature (Fig. 4 and fig. S6) likely reflects such dynamics, because relatively poor coherence with Halema‘uma‘u reservoir pressure (or lava lake elevation) indicates neither Xavg nor temperature is primarily driven by magma mass balance. Simple thermal arguments suggest likely causes of temperature variation. Atmospheric heat exchange at the lake surface will be dominated by radiative heat flux , where ϕr is ∼1 gigawatt (GW) for lake surface area A≈ 104 m2, thermal emissivity ϵ ≈ 0.8, Stefan-Boltzmann constant σ = 5.7×10−8 W m−2 K−4, and average surface temperature Ts≈ 700°C (). Heat flux to the host rock depends on hydrothermal circulation, but can be approximated with an effective thermal conductivity ϕc = keΔT/ΔL, where ϕc is 10 to 1000 W/m2 for ke of 2 to 20 W m−1 C−1 () and temperature gradient ΔT/ΔL of 10 to 100°C/m (). Total heat transfer rate Φ from the conduit and lake (surface area ∼105 m2) and from the Halema‘uma‘u reservoir (surface area ∼107 m2) is 1 to 100 megawatt (MW) and 0.1 to 10 GW, respectively. Neglecting latent heat, average temperature of a magma mass M will decrease as dT/dt = Φ/(cpM). For specific heat cp ≈ 1000 J kg−1 K−1, average temperature of the ∼1010 kg of magma in the conduit and lake could decrease by ∼10°C/day, whereas average temperature of the ∼1013 kg of magma in the Halema‘uma‘u reservoir would only decrease by ∼0.01 to 1°C/month. We thus expect the prevalent temperature drops of 100°C or more that occur over days to weeks represent downwelling of magma that cooled in the upper lava lake. Episodic downwelling suggests episodically decoupled convection cells in the lava lake rather than a convective regime that settles persistently into one of the configurations previously proposed (, ). This mechanism likely explains the ∼100°C temperature drop preceding the March 2011 Kamoamoa fissure eruption, where a changing convective regime is perhaps related to the rapidly filling lava lake and/or high short-term (hours to days) variability in lava lake elevation during this time. In some other cases, rapid lava lake draining might also induce downwelling of cool magma. This downwelling could explain the days-long temperature decreases accompanying the October 2012 and May 2014 intrusions, although if so, it is interesting that the 2015 intrusion did not cause a temperature drop.

An evolving magma plumbing system geometry

Given a consistent open hydraulic connection between the Halema‘uma‘u reservoir and lava lake, the weakening coherence between them over years or longer (Fig. 4) could represent changes either in the magma column density or in the relation between reservoir pressure and ground deformation (a function of geometry and poroviscoelastic rock properties). Our fixed geometry inversions test the former and show that for a range of feasible fixed parameter values (fig. S5), very high values of Xavg and/or ΔX are required over some portion of the timeline (e.g., 2009 through mid-2010 for reference parameters). These volatile contents would correspond to a foam in the upper conduit and lava lake with an average porosity in excess of 90%. Available constraints from gravity data () suggest average porosity in the lava lake of only up to 70%, so the higher values inferred at early times are likely unrealistic. We thus expect subsurface magma plumbing system geometry evolved over time, which could also contribute to the weak coherence between inverted South Caldera and Halema‘uma‘u reservoir pressures over long time scales (fig. 4).

Changes in conduit length (reservoir-roof depth) of ∼10 m or changes in conduit radius of ∼1 m could measurably affect VLP resonance period and quality factor at Kīlauea (fig. S5). Such changes might occur gradually because of processes such as viscous deformation of the host rock, thermal/mechanical erosion, or crystallization. Geometry could also change abruptly because of host rock failure or opening/closing of hydraulically connected dikes/sills. To fit the low VLP periods in 2009–2010 with realistic volatile contents, a ∼100-m-higher reservoir roof elevation (510 instead of 410 m above sea level, which is within estimated uncertainty ()) and/or strongly tapered conduit (e.g., top radius <5 m and bottom radius >15 m) is required (fig. S9). It is unlikely that the roof of an ellipsoidal reservoir would have grown downward this much over year time scales because of crystallization, so it may have been shallower throughout the eruption. In this case, the drastic change in VLP periods over the early part of the eruption likely represents an evolving conduit geometry due to some combination of a widening upper conduit and a change in conduit length due to a changing dip angle and reservoir attachment depth. A shallow dike/sill above the main Halema‘uma‘u reservoir could have also impacted the resonance (); this would potentially be consistent with some seismic inversions (, ), but such additional source complexity is not needed according to other seismic and geodetic inversions (, , ).

Toward a new generation of volcano monitoring

Resolving the dynamics of subsurface magma transport is a grand challenge that dictates hazard forecasting efficacy as well as connections between active volcanic processes and the geologic record. Inferring relative changes in magma properties over days to months by identifying the fluid origin of VLP seismic events represents a concrete step toward unifying the inversion of geophysical and geochemical data. In particular, we have resolved temperature changes of over 100°C that likely reflect both convective overturns and magma recharge. We have also resolved stratified volatile profiles that represent a highly disequilibrium outgassing regime. Volatile contents vary by over 1 wt % on time scales from days to months, revealing an unsteady shallow volatile mass balance. We have also inferred an evolving magma system geometry, highlighting the need to develop models and data sets that can deconvolve changing fluid properties from changing transport pathways.

Incorporating additional data would yield even more precise constraints on multiphase magma properties and their depth variation. For example, continuous gravity data would provide independent constraints on magma density in the lake. Video of lake surface oscillations could independently constrain vertical motions of the lake and triggering mechanisms of VLP events. In addition, surface gas emission data could constrain volatile stratification and outgassing/convective regimes if combined with models for gas flux through the magma column.

Similar VLP events have been detected at Vanuatu and Erebus volcanoes (, ) and are expected at open-vent volcanoes generally (), suggesting that this type of analysis could be adapted to improve near real-time monitoring at other eruptions. These data will inform basic volcano science and lead to better understanding of physical controls on volcanic eruptions.

MATERIALS AND METHODS

GNSS inversions

To obtain time series of pressure change in the Halema‘uma‘u reservoir, we must consider other known sources of ground deformation at the Kīlauea summit: the South Caldera reservoir (, ), 2015 intrusion (), and steady slip along the south flank décollement (fig. S3) (). We assume a temporally fixed geometry for the three magma reservoirs (Fig. 1 and Supplementary Text) but constrain the 2015 intrusion to be an active deformation source only over May 13 to 17 (). We adopt the 2-km-deep 4-km3 ellipsoidal Halema‘uma‘u reservoir geometry and 3-GPa rock shear modulus from (), consistent with other studies (, , , , ). We assume a horizontal centroid location of the South Caldera reservoir based on inversions of (); depth and geometry are less well constrained, so we choose a reference 20-km3 sphere centered 4 km deep and test different values based on published ranges (, , ). We fixed the 2015 intrusion geometry following ().

Reservoir pressures are found using linear least square inversions (Supplementary Text) of daily average surface position solutions from the Nevada Geodetic Laboratory () for GNSS stations within a few kilometers of the reservoirs (Fig. 1), corrected for steady background south flank slip with the multicomponent dislocation model of () (figs. S3 and S4). We use an approximate solution for deformation associated with a pressurized ellipsoid in an elastic half space () for each of the three magma bodies.

Conduit-reservoir magma oscillation model

We model VLP seismic events as small amplitude, isothermal, and incompressible oscillatory magma flow within a lava lake–conduit–reservoir system. The model is extended from () to include inertial effects in the lava lake and experimentally constrained models for multiphase magma properties (Supplementary Text). We consider an inclined radially symmetric magma column encompassing the lava lake and conduit, underlain by a reservoir within elastic rocks (Fig. 1).

The magma column prior to VLP events is assumed magmastatic, justified because fluid particle velocities associated with resonance are larger than background exchange flow (). During VLPs, viscous drag is determined from shear stress at the magma column wall where a no-slip velocity condition is enforced. With z and r distance parallel and perpendicular to the magma column axis (a function of conduit dip from horizontal θ), linearized conservation of momentum (primed variables) around a background state (bars) is

Here, 〈u′〉 is the cross-sectionally averaged conduit-parallel fluid particle displacement (so the orientation of 〈u′〉 is a function of θ), v′ is the conduit-parallel fluid particle velocity, 〈v′〉 is the cross-sectionally averaged v′ (the time derivative of 〈u′〉), ρ is the magma density, p′ is the pressure perturbation, μ is the dynamic viscosity, and R is the conduit radius. Conservation of mass is , where subscript 0 indicates evaluation at the bottom of the magma column (Fig. 1).

We assume equilibrium joint solubility of CO2 and H2O in Halema‘uma‘u composition melts () as a function of pressure and gas composition () (Supplementary Text and fig. S2). We neglect other volatile species as they have generally lower concentrations and/or poorly constrained solubility at Kīlauea (, ). We assume ideal gas behavior and consider melt density a function of pressure, temperature, and composition (). Melt viscosity μ(z) is assumed to be a function of temperature and dissolved H2O (). The impact of bubbles on apparent magma viscosity depends on the magnitude of capillary forces (). For expected strain rates of ∼10−1 s−1 associated with slow exchange flow in the conduit, bubbles less than ∼10 cm across will increase apparent viscosity approximately according to (fig. S2), where is the background magma porosity ().

For conduit-reservoir resonance, pressure at the base of the magma column is (), where Cr is the total storativity of the reservoir (reservoir volume change per unit pressure increase). The Halema‘uma‘u reservoir assumed here corresponds to a “buoyancy-dominated” limit where reservoir pressure changes have a negligible effect on the magma column during VLPs (Supplementary Text) (). Pressure at the top of the magma column is , where subscript H indicates evaluation at the top of the magma column, and Pex(t) is the external forcing (Fig. 1). This system is equivalent to a driven harmonic oscillator with frequency-dependent damping and exhibits exponentially decaying oscillations in response to an impulsive forcing (fig. S1). We find the resonant period and quality factor by solving numerically for the free response of the system (Supplementary Text).

VLP seismic event inversions

We assume a temporally fixed magma plumbing system geometry, except for lava lake radius and surface elevation, which are interpolated from measurements (Supplementary Text) (, ). We choose reference fixed parameters based on previous constraints where available. Where minimal constraints are available, we test a range of values and select combinations that produce feasible inversion results over most of the timeline, as detailed in the Supplementary Text. We approximate the lava lake and conduit as cylinders, with a reference conduit radius of 15 m and conduit dip of 90° from horizontal (Fig. 1 and table S1).

We conduct inversions using the conduit-reservoir resonance model for the three free parameters (temperature, Xavg, and ΔX) from the three target values for each VLP seismic event: conduit bottom (Halema‘uma‘u reservoir top) pressure, resonance period, and resonance quality factor (Fig. 2). We use an iterative nonlinear trust-region-reflective solver to find the combination of free parameter values that minimizes misfit Ewhere vertical bars indicate absolute value, asterisks indicate observed/target values, Q is the resonance quality factor, ω is the resonance angular frequency, and is the magmastatic pressure at the bottom of the conduit (top of the reservoir). To prevent unfeasible solutions, we impose bounds on the search space such that volatile mass fraction at all depths is between 0 and 7 wt % and temperature is between 900° and 1600°C. In most cases, there is an exact solution (E = 0), although for some VLP events (e.g., in 2009 and early 2010), exact solutions do not exist for the reference parameters, and the solver will find a local minimum instead. Grid searches indicate that the misfit spaces are convex, so the solver is finding unique global minima and/or unique exact solutions. Time-series analysis methods used to interpret inversions are detailed in the Supplementary Text.