Explosive Activity on Kīlauea's Lower East Rift Zone Fueled by a Volatile‐Rich, Dacitic Melt

Magmas with matrix glass compositions ranging from basalt to dacite erupted from a series of 24 fissures in the first 2 weeks of the 2018 Lower East Rift Zone (LERZ) eruption of Kīlauea Volcano. Eruption styles ranged from low spattering and fountaining to strombolian activity. Major element trajectories in matrix glasses and melt inclusions hosted by olivine, pyroxene and plagioclase are consistent with variable amounts of fractional crystallization, with incompatible elements (e.g., Cl, F, and H2O) becoming enriched by 4–5 times as melt MgO contents evolve from 6 to 0.5 wt%. The high viscosity and high H2O contents (∼2 wt%) of the dacitic melts erupting at Fissure 17 account for the explosive Strombolian behavior exhibited by this fissure, in contrast to the low fountaining and spattering observed at fissures erupting basaltic to basaltic‐andesite melts. Saturation pressures calculated from melt inclusion CO2‐H2O contents indicate that the magma reservoir(s) supplying these fissures was located at ∼2–3 km depth, which is in agreement with the depth of a dacitic magma body intercepted during drilling in 2005 (∼2.5 km) and a seismically imaged low Vp/Vs anomaly (∼2 km depth). Nb/Y ratios in erupted products are similar to lavas erupted between 1955 and 1960, indicating that melts were stored and underwent variable amounts of crystallization in the LERZ for >60 years before being remobilized by a dike intrusion in 2018. We demonstrate that extensive fractional crystallization generates viscous and volatile‐rich magma with potential for hazardous explosive eruptions, which may be lurking undetected at many ocean island volcanoes.

1. What was the cause of the high explosivity observed in the eruptions at the western end of F17? 2. What was the most probable parent magma(s) of the variable melt compositions erupted in the first few weeks of the 2018 eruption? 3. At what depth were these magmas stored in the LERZ?
We show that fissures erupting in the first few weeks of the 2018 eruption tapped a stored magma body which experienced variable amounts of fractional crystallization at 2-3 km depth. Based on a comparison of incompatible trace element ratios from 2018 samples and historic activity in this section of the LERZ, we show that the 2018 eruption tapped a magma body with very similar chemistry to that supplying the 1955 and 1960 eruption. The large amounts of fractional crystallization experienced by dacitic melts erupted at F17 drove up melt H 2 O contents to unprecedented levels at Kīlauea (∼2 wt%) which, combined with an increase in melt viscosity, accounts for the explosive activity observed at this fissure.

Chronology of the 2018 LERZ Eruption
The 2018 LERZ eruption marked the end of a 35-year-long period at Kīlauea Volcano during which activity was focused predominantly at the Puʻuʻōʻō vent, on the Middle East Rift Zone (the name change from PuʻuʻŌʻō signifies that the eruption has ended; Hawaiʻi Board on Geographic Names). Following inflation of the summit and Puʻuʻōʻō during March and April 2018 (Neal et al., 2018;Patrick et al., 2020), Puʻuʻōʻō vent collapsed on 30 April, accompanied by the propagation of a dike downrift into the Leilani Estates subdivision on the LERZ.
Following the formation of ground cracks, the first eruptive fissures opened on 3 May. Fourteen additional fissures opened along a linear trend in the first week (F1 to F15; Figure 1b). These first fissures, classified as Early Phase 1 by Gansecki et al. (2019), were characterized by whole-rock compositions with ∼4-5 wt% MgO, erupted as relatively small-volume spatter mounds and lava flows. These lavas are thought to represent material which was intruded during previous eruptive episodes in the LERZ (e.g., 1790, 1840, 1955, or 1960Gansecki et al., 2019;Figure 1c). The injection of the dike beneath the Leilani Estates subdivision at the start of May 2018, is thought to have disturbed these resident melts and forced them to the surface (e.g., by increasing overpressure, changing the crustal stress state, or thermal rejuvenation).
Following an eruptive pause, activity recommenced on 12-18 May (termed Late Phase 1), with more fissures opening up along the same linear trend (Fissures 16,18,20,and 22). These later fissures erupted magma compositions with higher bulk MgO contents (5-6 wt% vs. 4-5 wt% in Early Phase 1), which have been interpreted to represent mixing between the early rift-stored component, and the more MgO-rich magma supplied by the dike (Gansecki et al., 2019). Between 17 and 27 May (Phase 2) effusion rates increased and whole-rock MgO contents increased again to become very similar to those erupted at Puʻuʻōʻō vent (6-7 wt%), indicating that the vast majority of melts stored in the LERZ had been flushed from the system (Gansecki et  showing the difference in eruption style between the basaltic-andesite eastern end (low fountaining and spattering, LHS), and the more silicic andesitic western end (strombolian explosions and gas jetting, RHS). (e) Low spatter mounds produced during the first day of the eruption at F3. Photos in panels (d and e) from USGS. Panel (b) was adapted from Patrick et al. (2019), and panel (c) from Pietruszka et al. (2021), using geospatial data from Trusdell et al. (2006) and Zoeller et al. (2020 Gansecki et al., 2019) and previous Kīlauea eruptions (black triangles, Wieser et al., 2019) are shown, along with the whole-rock composition of the dacitic magma drilled in 2005 (Teplow et al., 2009). Two fractional crystallization models run in alphaMELTS for MATLAB (Rhyolite-MELTS v1.2.0, Antoshechkina & Ghiorso, 2018;Ghiorso & Gualda, 2015;Gualda et al., 2012) are shown, with Fe 3+ /Fe T, initial = 0.15 (unbuffered), P = 650 bars. Model A uses an initial H 2 O content of 0.1 wt% (red line, best recreates the onset of Fe-Ti oxide saturation), while Model B uses a more realistic initial H 2 O content of 0.5 wt% (black dashed line). As the rhyolite-MELTS model fails to recreate the olivine-only fractionation assemblage at >6.8 wt% (Wright & Fiske, 1971), the red dashed line shows a reverse fractionation crystallization model run from the starting composition in Petrolog3.1.3 (Danyushevsky & Plechov, 2011). The fractionation trajectory for Zr is calculated assuming incompatible behavior.  It has been suggested, based on the presence of high forsterite olivines incorporated within Phase 3 lavas, and the bimodal distribution of melt inclusion saturation pressures at 1-2 and 3-5 km, that this magma originated from the two magma reservoirs identified by geophysical imaging beneath Kīlauea's summit (Lerner et al., 2021;Wieser et al., 2021), with a hydraulic connection linking these reservoirs to Ahuʻailāʻau (Gansecki et al., 2019;Neal et al., 2018;Patrick et al., 2019). Conversely, Pietruszka et al. (2021) suggest that major and trace element differences between Phase 3 lavas and those erupted at Kīlauea's summit and Puʻuʻōʻō prior to the onset of the 2018 eruption rules out a summit origin, and instead suggest that this MgO-rich magma had accumulated in the Middle East Rift Zone (MERZ) downrift of Puʻuʻōʻō over ∼10 years prior to 2018. The exact source of Phase 3 magma is beyond the scope of this study, which focuses on the early phase lavas erupted between 3 May to 16 May, so we refer to Phase 3 magma as "dike-supplied" to avoid ambiguity.

Analytical Methods
Spatter, reticulite, and lava samples were collected from nine separate fissures (Fissures 4,5,8,10,11,12,13,17,and 20) on Kīlauea's LERZ during the 2018 eruption and a follow-up field campaign in 2019 (Table S1 in Supporting Information S1). Care was taken to select the most rapidly quenched water and air-quenched material to minimize H + diffusion out of melt inclusions (Gaetani et al., 2012). Specifically, we selected the smallest and most visibly glassy clasts available in each location. Clast size was 0.5-8 cm and water-quenched material was favored where available. Samples were jaw-crushed, and olivine, orthopyroxene, clinopyroxene, and plagioclase crystals were picked under a binocular microscope. Crystals were individually mounted on glass slides using Crystalbond and ground down to the level of target melt inclusions. Care was taken to prioritize melt inclusions not visibly connected to the outside of the crystal, but particularly in plagioclase-hosted melt inclusions, it was often difficult to determine the degree of connectivity due to complex melt networks. An indication of whether each melt inclusion was fully isolated or connected to the carrier melt is provided in Data Set S1 in Supporting Information S2.
All vapor bubbles within melt inclusions were analyzed for CO 2 using Raman spectroscopy following the methods described by Wieser et al. (2021). Unlike the olivine-hosted melt inclusions erupted in the Phase 3 between 28 May and 1 August, none of the vapor bubbles in melt inclusions from the early fissures (3-16 May) investigated in this study produced measurable CO 2 peaks, indicating that the concentration of CO 2 in their vapor bubbles were below the detection limit of Raman spectroscopy (∼0.02 g/cm 3 ; Wieser et al., 2021). Following exposure of a target melt inclusion at the surface, crystals were mounted in epoxy stubs in groups of 30-40 inclusions, and polished using progressively finer diamond pastes. Samples were analyzed for H 2 O, CO 2 , and F (as well as MgO and SiO 2 for normalization) using the Cameca IMS-7f GEO at the NERC Ion Microprobe Facility at the University of Edinburgh in two analytical sessions (July 2019 and January 2020). A variety of glass standards were analyzed with a range of Si contents to convert counts to element concentrations, and to account for matrix effects on ion yields as a function of melt SiO 2 content. Detailed analytical information is available in the Figures S1, S2 and Table S2 in Supporting Information S1.
Following SIMS analysis, the Au coat was removed by polishing on a 1 μm pad (which also helped to reduce the size of SIMS pits and remove SIMS sputter) and a C coat was applied. Matrix glasses, melt inclusions, and the composition of the mineral host ∼30-100 μm away from each melt inclusion were analyzed using a Cameca SX100 EPMA at the Department of Earth Sciences, University of Cambridge following the analytical techniques described in Wieser et al. (2021) (also see Text S1 and Table S3-S9 in Supporting Information S1). Care was taken to analyze away from the SIMS pit where possible.
Trace elements were measured in melt inclusions and matrix glasses using a Photon Machines G2 193 nm excimer laser system equipped with a HelEx II 2-volume cell coupled to an Agilent 8800 ICP-MS/MS) at the School of Environment, Earth and Ecosystem Sciences, The Open University. Depending on the size of the melt inclusion and the number of microlites in matrix glasses, ablation spots were between 20 and 50 μm. Analyses were conducted following the techniques described by Jenner et al. (2015) and Wieser et al. (2020). Comparison of measurements of BCR-2G run as secondary standard to long-term laboratory preferred values is shown in Data Set S2 in Supporting Information S2.

Correction for Post-Entrapment Crystallization
Olivine-hosted melt inclusions were corrected for the effects of post-entrapment crystallization (PEC) using the "Olivine MI" tool in Petrolog3.1.1.3 (Danyushevsky & Plechov, 2011), which requires users to enter the measured major element composition of the melt inclusion, the host Fo content (both taken from EPMA analyses), and an estimate of the initial FeO t content of the melt inclusion prior to PEC. The initial FeO t content was estimated by determining a relationship between olivine Fo content and melt FeO t contents by calculating the equilibrium olivine Fo content for co-erupted matrix glasses (see Text S2 and Figure S3 in Supporting Information S1). While olivine-hosted melt inclusions erupted in the later stages of F8 (28 May-1 August) have experienced up to 35% PEC , the maximum amount of PEC experienced by olivines erupted at the earlier fissures examined in this study is 10% (median = 4%).
Unlike for olivine, there is no clear consensus as to the best way to reconstruct orthopyroxene, clinopyroxene, and plagioclase-hosted melt inclusion compositions for the effects of PEC. We discuss the approach taken in this study in detail in Text S2 in Supporting Information S1. Briefly, for clinopyroxene we examine the degree of disequilibrium between each melt inclusion and its host clinopyroxene crystal using the four equilibrium tests described by Neave et al. (2019): Fe -Mg , Di-Hd, En-Fs, and CaTs. Melt inclusions lie within the equilibrium field for all parameters except Fe -Mg , ( Figure S4 in Supporting Information S1, likely resulting from the fact that the concentration of FeO and MgO are far more sensitive than other major elements to the crystallization of small amounts of clinopyroxene on the wall of the melt inclusion). If a melt inclusion has crystallized clinopyroxene on its walls, the composition of that clinopyroxene must be added back into the measured melt inclusion composition to restore the composition at the time of entrapment. Similarly, if the inclusion was heated up and clinopyroxene dissolved following entrapment, that clinopyroxene must be removed from the measured inclusion composition. We add/subtract the composition of the host clinopyroxene until the melt inclusion and host crystal are in Fe -Mg equilibrium using Equation 35 of Wood and Blundy (1997). The median amount of clinopyroxene addition needed to obtain Fe -Mg , equilibrium was 4 wt% (1σ = 5%; Figure S5 in Supporting Information S1). Corrected melt inclusion compositions meet all four equilibrium tests, and lie closer to major element trajectories defined by co-erupted matrix glasses than uncorrected compositions ( Figures S4 and S5 in Supporting Information S1).
Many orthopyroxene-hosted melt inclusions also lie out of Fe -Mg , equilibrium with the composition of the host crystal ( Figure S6 in Supporting Information S1). However, unlike for clinopyroxene, adding/subtracting the composition of the host crystal to reach Fe -Mg equilibrium results in a worse fit to the major element trajectories defined by matrix glasses in CaO, and TiO 2 versus MgO space ( Figure S7 in Supporting Information S1). Given that most orthopyroxene-hosted melt inclusions lie close to the major element trajectory of matrix glasses, we leave their compositions uncorrected.
The vast majority of plagioclase-hosted melt inclusions lie very close to the compositional trend defined by matrix glasses, with the exception of six melt inclusions erupted at F17, which are offset to significantly higher FeO ( Figure 2b) and lower CaO contents at ∼6 wt% MgO ( Figure S8 in Supporting Information S1). These offsets are indicative of crystallization of ∼15 wt% plagioclase on the walls of the melt inclusion following entrapment, likely due to extensive cooling between the formation of these relatively high anorthite crystals (∼An 60 ) and their incorporation in the significantly cooler F17 melt composition (in equilibrium with ∼An 45−50 ). These six melt inclusion compositions were corrected by adding the composition of the host crystal back into the measured melt inclusion composition to obtain the best fit to matrix glass data in MgO versus Al 2 O 3 space ( Figure S9 in Supporting Information S1). All other plagioclase-hosted melt inclusions were left uncorrected.

Major and Volatile Elements
The composition of melt inclusions and matrix glasses erupted between 3 and 16 May define coherent trends in major element space, spanning basaltic to dacitic compositions (∼7-0.5 wt% MgO, 48-68 wt% SiO 2 ; Figure 2a). The compositions of matrix glasses are distinctly clustered, with the largest group at ∼4 wt% MgO (erupted at F4, F5, F8, F10-12 during Early Phase 1 on 3-9 May), a second cluster at ∼1 wt% MgO (erupted at F17 on 14 May), and a third cluster at ∼6-7 wt% MgO erupted in Phase 3 at Ahuʻailāʻau (F8) after 28 May. Lavas erupted from F13 following its re-activation (F13-react) plot between the Early Phase 1 and F17 cluster, while lavas erupted at F20 on 16 May plot between the Early Phase 1 cluster and the higher MgO Phase 3 cluster.
Melt inclusion and matrix glass Cl and F concentrations increase with decreasing melt MgO contents, with the most volatile-rich dacitic melt inclusions erupted at F17 having Cl and F concentrations of ∼1,000 ppm (3-5 times higher than observed in melt inclusions erupted in Phase 3; Figures 3a and 3b). F concentrations show significantly more scatter than Cl concentrations. The upper limit of melt inclusion H 2 O contents also increases with fractional crystallization, reaching 1 wt% in the early Phase 1 lavas with ∼4 wt% MgO, and up to 2 wt% in melt inclusions erupted at F17 with 0.5-1 wt% MgO (Figure 3d). These  . For Cl and F, models A and B produce very similar trajectories, as these models predict similar amounts of solid fractionated. The solubility of pure CO 2 and H 2 O calculated using the MagmaSat solubility model (Ghiorso & Gualda, 2015) implemented in VESIcal  for the best-fit to the observed major element path (see Figure S11 in Supporting Information S1) is shown with a blue dashed line.

Trace Elements
As discussed further in Section 4.4, Nb/Y is not extensively fractionated during crystallization at Kīlauea, so can be used to identify distinct magma batches . Early Phase 1 matrix glasses from F4, F5, F8, and F10-12 (green symbols) have Nb/Y ratios of ∼0.6-0.8, which overlap with the composition of matrix glasses erupted at F17 (Figure 4a). In contrast, matrix glasses erupted at F13 following its reactivation and F20 have lower Nb/Y ratios, plotting closer to the composition of Phase 3 glasses. Overall, the trends defined by matrix glasses measured by LA-ICP-MS in this study are very similar to those determined by whole-rock XRF measurements , magenta dots, Figure 4a).
Unlike the temporal variations seen in matrix glass and whole-rock trace element compositions, there are no apparent inter-sample differences in the Nb/Y ratios of melt inclusions erupted from early fissures, or relationships between Nb/Y ratios and host crystal An# or Mg# in these samples (Figures 4c and 4d). As for matrix glasses, the vast majority of melt inclusions from Phase 3 (orange colors; Figures 4c and 4d) lie to distinctly lower Nb/Y  Pietruszka & Garcia, 1999a, 1999b, LERZ lavas from 1790, 1823, 1840, 1955, 1960(red dots, Marske, 2010, early phases of the 1955 LERZ eruption (yellow colors), and later phases (cyan colors, Helz & Wright, 1992;Norman, 2005). Analyses of matrix glass from the 1960 Kapoho eruption from  and whole-rock measurements from Puʻuʻōʻō are also shown (Gansecki et al., 2019;Garcia et al., 2021).

Matrix Glasses (this study)
ratios, with only a small number of melt inclusions erupted on 28 May overlapping with the Nb/Y ratios of melt inclusions from early fissures (particularly those hosted in lower forsterite olivines).

Modeling Fractional Crystallization
The coherent trends in major and volatile element trajectories defined by melt inclusions and matrix glasses indicate that the dominant cause of chemical variability in erupted melts during the first 2 weeks of the 2018 eruption was fractional crystallization. Previous studies have shown that Kīlauea melts fractionate only olivine (+minor Cr-spinel) at >6.8 wt% MgO, defining a relatively flat trajectory for all elements versus MgO (Wright & Fiske, 1971). Next, clinopyroxene joins the liquidus, followed by plagioclase after only 5-15°C cooling, as well as minor orthopyroxene (Helz & Wright, 1992). The crystallization of plagioclase accounts for the rapid increase in melt FeO T as MgO contents drop from 6 to 4 wt% ( Figure 2b). The appearance of plagioclase and clinopyroxene on the liquidus also causes the slope of relatively incompatible elements like TiO 2 and Zr versus MgO to greatly increase, because a larger mass of solid is fractionated for a given decrease in melt MgO content. At ∼4 wt% MgO, Fe-Ti oxides saturate (Helz & Wright, 1992), driving the concentration of FeO T and TiO 2 down, and SiO 2 up. Concurrently, P/Nd begins to decline (along with P 2 O 5 ; Figure S10 in Supporting Information S1), indicating the onset of apatite crystallization (visible as microphenocrysts in matrix glasses from F17).
Although it is likely that the melt fractionating within the LERZ was a previously stored component, the starting composition of models was set as the average major-element composition of glasses erupted in Phase 3. This is because the vast majority of literature data for previous eruptions comprises whole-rock compositions, which may be affected by crystal addition. The starting composition used for the modeling lies within the range of whole-rock compositions from various early fissures from Gansecki et al. (2019), justifying this choice. The initial Fe 3+ /Fe T ratio was set at 0.15 (Lerner et al., 2021;Moussallam et al., 2016), and models were run without a buffer, allowing for changes in oxidation state following Fe-Ti oxide saturation. The pressure was set at 650 bars, based on the range of saturation pressures measured in Early Phase 1 melt inclusions (see Section 4.5). The trajectory of trace and volatile elements during fractional crystallization was modeled using the mass of solid phases produced, and the fractional crystallization equation (C final = C initial × F D−1 , where F is the amount of melt remaining, and D is the bulk partition coefficient).
In detail, we consider two separate model trajectories in MELTS. Model A (red line) was run with an initial H 2 O concentration of 0.1 wt%. Model A is a good match to the trajectory of major elements from ∼6.8 to 4 wt% MgO, capturing the sharp increase in FeO t and TiO 2 contents with decreasing MgO. At lower MgO contents, Model A is a good fit to FeO t , TiO 2 , and SiO 2 (Figures 2a-2c) but overestimates CaO, and underestimates Al 2 O 3 ( Figure  S11 in Supporting Information S1) indicating that the ratio of clinopyroxene versus plagioclase fractionating in this model is incorrect. Similarly, Model A does not recreate the prominent downtick in P 2 O 5 (and P/Nd) caused by apatite crystallization, because MELTS currently only accounts for a hydroxy-apatite phase (rather than the more stable fluorapatite observed in natural magmas; Rooney et al., 2012). The failure to correctly model apatite saturation may also account for some of the overestimation of melt CaO contents at <4 wt% MgO ( Figure S11 in Supporting Information S1). However, despite these major-element discrepancies, the concentration of Zr and Cl (assuming complete incompatibility) in Model A provides a good fit to the composition of melt inclusions and matrix glasses. This indicates that MELTS model results are accurately predicting the amount of solid being fractionated, even if predictions of relative phase proportions are slightly inaccurate.
If MELTS models are initialized with water contents similar to those in melt inclusions with ∼6 wt% MgO (H 2 O i = 0.5 wt%, Model B, Figure 2, dashed black lines), MELTS saturates Fe-Ti oxides significantly earlier and plagioclase later than Model A, such that Model B fails to achieve the prominent TiO 2 and FeO enrichment seen in erupted glass compositions. This has been noted previously for Kīlauea by Garcia et al. (2003). Fortunately, the total mass of solids predicted by model A and model B are reasonably similar, so these models predict similar trajectories for incompatible elements such as Zr (Figure 2d). We also generate a third major element path for fractional crystallization, using MELTS model A at >4 wt% MgO, and a best fit through measured glass and melt inclusion data for <4 wt% ( Figure S11 in Supporting Information S1). The major element differences between these three models have a relatively small effect on volatile solubility ( Figure S12 in Supporting Information S1) and calculations of melt viscosity (Figure 5b). This justifies the use of Model B to investigate changes in volatile systematics and viscosity during fractional crystallization, even though this model doesn't fully recreate the observed major element systematics ( Figure 2).
Rhyolite-MELTS V1.2.0 fails to reproduce olivine-only fractionation at >6.8 wt% MgO, so this section of the fractionation path was modeled using the "reverse crystallization" tool in Petrolog3 (Danyushevsky & Plechov, 2011) from the specified starting composition.

Volatile Element Evolution During Crystallization
MELTS modeling demonstrates that the prominent enrichment in melt inclusion Cl concentrations with decreasing MgO results from extensive fractional crystallization of Cl-poor crystals ( Figure 3a). F data shows significantly more scatter than Cl (Figure 3b), but approximately follows the fractional crystallization trend. Notably, six plagioclase-hosted F17 melt inclusions with higher MgO contents (∼4-6 wt%) lie significantly above the fractional crystallization line for F, indicating that these inclusions may have diffusively re-equilibrated with the significantly more F-rich carrier melts in which they were erupted ( Figure 3b). Apparent open-system behavior of F in olivine-hosted melt inclusions has been noted previously by Koleszar et al. (2009) and Portnyagin et al. (2008), while Neave et al. (2019) attribute F enrichment in plagioclase-hosted melt inclusions to diffusive pile up during rapid growth. We favor diffusive re-equilibration here, because many F17 melt inclusions have lower F contents than co-erupted matrix glasses, and it is only the most MgO-rich melt inclusions which are expected to have equilibrated with more F-rich, MgO-poor melts that show notable enrichment above the fractional crystallization trend. Two of these melt inclusions also show elevated H 2 O contents (Figure 3d), which is consistent with these more primitive melt inclusions equilibrating with the more MgO-poor, F-and H 2 O-rich F17 carrier melt prior to eruption (with F re-equilibration appearing to occur faster than H 2 O re-equilibration in plagioclase).
Accounting for these complexities with F, the relative similarity of the trajectories of matrix glasses and melt inclusions indicates that the amount of F and Cl degassed upon eruption is sufficiently small that it cannot be resolved petrologically by comparing melt inclusions and matrix glasses. In contrast, H 2 O and CO 2 degas strongly upon eruption, shown by the significantly lower concentrations in matrix glasses relative to melt inclusions (Figures 3c and 3d).   (1987) MgO thermometer. The red line shows the viscosity calculated using Giordano et al. (2008) for the best-fit model to measured major elements ( Figure S11 in Supporting Information S1), with temperatures calculated using Helz and Thornber (1987), and H 2 O contents estimated for a parameterization of the relationship between H 2 O and MgO shown as a red line in Figure 3d.
in a mixed H 2 O-CO 2 fluid at crustal pressures (Gerlach et al., 2002). This means there is a competing influence between their incompatible behavior in solid phases and their loss to an exsolving fluid phase. To demonstrate this trade-off, we use the solubility model MagmaSat (Ghiorso & Gualda, 2015) implemented in the Python3 tool VESIcal  to calculate the solubility of pure-CO 2 fluids and pure-H 2 O fluids at 650 bars for the synthetic major element path which provides the best fit to measured glass and melt inclusion data ( Figure  S11 in Supporting Information S1), using temperatures from the MgO geothermometer of Helz and Thornber (1987). At low pressures where non-ideality between H 2 O and CO 2 is negligible, the concentration of either volatile species will not exceed the pure solubility limit. The solubility of pure CO 2 is relatively low (∼400 ppm), and plots close to the measured CO 2 contents of melt inclusions (dashed blue line, Figure 3c). In contrast, the solubility of pure H 2 O is much higher, plotting significantly above the concentration of melt inclusions (dashed blue line, Figure 3d).
In reality, H 2 O-CO 2 exsolves as a mixed fluid in magmatic systems. At the low pressures considered here, where mixing between H 2 O and CO 2 is close to ideal, the behavior of mixed fluids is best understood considering Henry's Law (Lowenstern, 2001). Namely, the addition of H 2 O to a pure CO 2 fluid lowers the partial pressure of CO 2, and therefore lowers its solubility in the co-existing melt phase.  Figure 5a). This decrease in the proportion of CO 2 in the fluid causes the CO 2 solubility predicted by MELTS model B (black dotted line, Figure 3c) to decrease more rapidly than the model for pure CO 2 solubility. MELTS model B tracking mixed volatile solubility during fractional crystallization effectively brackets the upper limit of the majority of melt inclusion data, providing strong evidence that melts were vapor-saturated at depth, and that the solubility model MagmaSat used in Rhyolite-MELTS is effectively capturing changes in volatile solubility during fractional crystallization.
A significant proportion of the total amount of CO 2 in the system is lost to the vapor phase during fractional crystallization, partially because XCO 2 >65 wt% across most of the fractionation interval, and because there was so little CO 2 in the system to start with (∼385 ppm at 6.7 wt% MgO). For example, fractionation of just 0.02 wt% fluid with XCO 2 =80 wt% removes 160 ppm of CO 2 from the system. In contrast, the higher initial mass of H 2 O in the system combined with high XCO 2 ratios mean that very little H 2 O is lost to the fluid; in the example above only 40 ppm of H 2 O is lost to the vapor compared to 160 ppm of CO 2 . Notably, this is only ∼1% of the initial amount of H 2 O in the system, which is well within analytical error of SIMS measurements ( Figure 5b) Figure 6). Some samples (e.g., F10; Figure 6b) show tightly clustered volatile contents, which plot very close to fractional crystallization trends (Figures 3c and 3d). Others (e.g., F13-react, Figure 6a) show considerable spread to lower CO 2 contents, approximately following the trajectory of modeled degassing paths. This indicates that many melt inclusions were formed or sealed during ascent, trapping a progressively degassing magma. Horizontal spread to variable H 2 O contents may reflect diffusive re-equilibration of melt inclusions with degassing carrier melts, or crystals being entrained into melts with higher H 2 O than those they grew from (e.g., the 6 plag-hosted MI from F17). F17 melt inclusions show the most scatter, which can be recreated by a combination of the trajectory shown by MELTS model B, degassing, and H + re-equilibrium ( Figure 6d). However, without measurements of D/H ratios in melt inclusions, the exact contribution of degassing, fractional crystallization, and re-equilibration of H 2 O is difficult to determine. In general, we note that H + loss through plagioclase is likely slower than in olivine (Johnson & Rossman, 2013). Additionally, we observe that multiple melt inclusions hosted within single crystals erupted at F17 have a wide range of H 2 O-CO 2 contents, and there is no relationship between H 2 O contents and the position of the melt inclusion within the crystal, or the size of the inclusion (inconsistent with H-loss by diffusion, Figures S13-S15 in Supporting Information S1). The observed variability within individual crystals is more consistent with melt inclusions being trapped at various points along a degassing path. The influence of degassing is considered further in Section 4.5.

Implications for Eruption Style Variations
The higher water contents of F17 melts relative to basaltic and basaltic-andesite melts (resulting from extensive fractional crystallization) means that a significantly larger volume of volatiles is exsolved as these melts ascend to the surface (Figure 7c). We demonstrate this by tracking the ascent of magma along a closed-system degassing path in alphaMELTS for MATLAB for a composition representative of the basaltic material erupted at the early fissures (MgO = 4.3 wt%, SiO 2 = 50.8 wt%, and H 2 O = 0.95 wt%) and a composition representative of the dacitic melt erupted at F17 (MgO = 1 wt%, H 2 O = 1.78 wt%; Figure 7). For both melt compositions, these models predict that the volume of exsolved volatiles relative to the total volume of the system (fluid + melt + crystals) is very low at 650-400 bars (Figure 7a), and degassing only produces measurable changes in the dissolved concentration of CO 2 in the melt (Figure 7b). This is also apparent from the degassing paths in Figures 6a-6d. At Figure 6. (a-c) Example degassing paths for three fissures, with open and closed system degassing paths calculated using MagmaSat (Ghiorso & Gualda, 2015) in VESIcal    To compare these ascent paths, we divide the volume of volatiles exsolved during ascent of the dacite by the volume of volatiles exsolved during ascent of the basaltic (Figure 7a). The higher initial H 2 O content of the dacite, combined with the fact that H 2 O begins to exsolve at higher pressures means that at ∼200 bars, the volume of volatiles in the dacitic melt is 8-9 times higher than in the basaltic melt (Figure 7c). The larger volatile volume at depth, combined with higher melt viscosity, will promote bubble nucleation and coalescence, which is a key driver of Strombolian activity (Jaupart & Vergniolle, 1988). By the time these melts reach the surface, where almost all of the initial H 2 O in both melts has exsolved, the volume of volatiles exsolved in the dacitic melt is approximately twice that for the rift-stored basalts.

Pass filter
Eruption style is also influenced by melt viscosity. First, we model changes in melt viscosity during fractional crystallization (Figure 5c), and then assess changes in viscosity during ascent using the viscosity model of  Giordano et al. (2008) model. The dashed line shows viscosity using the temperature from MELTS, and the solid line shows temperature calculated from the MELTS MgO content using Helz and Thornber (1987). The yellow star shows the viscosity ratio between the dacitic F17 lavas and the basaltic/basaltic-andesite Early Phase 1 lavas reported by Soldati et al. (2021), using temperatures they calculate from whole-rock MgO contents using Helz and Thornber (1987). The cyan star uses the viscosity parameters (A, B, C) of Soldati et al. to calculate viscosity at the temperatures for the MgO contents considered here for the dacitic melt (1 wt% MgO) and basaltic melt (4.3 wt% MgO) using Helz and Thornber (1987). The lower MgO content of our F17 glass relative to their whole-rock measurement result in a lower temperature, so a higher viscosity.   Giordano et al. (2008) (Figure 7d). As viscosity is highly sensitive to melt H 2 O content, we calculate viscosity for the major and volatile-element trajectories from MELTS Model B, using MELTS temperatures and temperatures from Helz and Thornber (1987). Because Model B is a relatively poor fit to major element data at <4 wt% MgO, we also calculate the viscosity for the best fit to measured major and volatile element data (see Figure S11 in Supporting Information S1) using the MgO geothermometer of Helz and Thornber (1987) for melt temperature, and H 2 O contents from a best fit to MELTS model B. The two viscosity models using Helz and Thornber (1987) are very similar, predicting a rapid increase in melt viscosity below 4 wt% MgO (Figure 5b) as a result of SiO 2 enrichment in the melt following the onset of Fe-Ti oxide fractionation ( Figure 2a) and a steady decline in the temperature of the melt. The model using MELTS temperatures is similar until ∼2 wt% MgO, after which it rises to higher values because MELTS predicts significantly lower temperatures than Helz and Thornber (1987).
Prior to the onset of degassing upon ascent, the dacitic melt is ∼10-20 times more viscous than the basaltic melt ( Figure 7d). During ascent along a closed-system degassing path, the discrepancies in viscosity between the basaltic and dacitic melt are enhanced, so at the surface, the dacitic melt is 70-130 times more viscous (Figure 7d; Soldati et al., 2021). The relative viscosities calculated in our models are very similar to viscosities obtained from direct measurements of 2018 lavas at atmospheric pressure (yellow and blue stars on Figure 7d; Soldati et al., 2021) particularly when the Helz and Thornber (1987) MgO thermometer is used for both comparisons. The increasing contrast in viscosity between the basalt and dacite during ascent to the surface results from the fact that the dacitic melt has more H 2 O to lose (which drives a rapid increase in viscosity). It is also worth noting that these calculations only track the viscosity of the melt phase. While the proportion of crystals relative to melt is similar for F17 versus early Phase 1 lavas (∼39% vs. 42%; Gansecki et al., 2019), F17 samples have completely different crystal shapes, with a larger number of smaller, elongated plagioclase crystals which may further amplify these differences in viscosity (see BSE images in Data Set S3, p. 78-80 vs. p. 82-84 in Supporting Information S2). The larger volume of bubbles in F17 melts will also affect viscosity, although this is difficult to quantify because it depends greatly on the amount of shear rate and strain of bubbles (Llewellin et al., 2002;Manga & Loewenberg, 2001).
It is well recognised that differences between Hawaiian style lava fountaining (exhibited weakly by the early fissures) and Strombolian activity (exhibited strongly by the F17 andesite) are controlled by magma rise speed, viscosity, and volatile content (Houghton et al., 2016;Houghton & Gonnermann, 2008;Wilson & Head, 1981). Strombolian eruptions occur when bubbles can coalesce and rise faster than the surrounding melt, resulting in the bursting of large bubbles at the surface along with relatively minimal volumes of melt (see Movie S1, Wilson & Head, 1981). Separated two phase flow is promoted at higher gas/magma ratios and for higher melt viscosities, which tends to reduce magma rise rates (Gonnermann & Manga, 2012). The 10-100 higher viscosities of F17 melts compared to more primitive basaltic melts at <300 bars (Figure 7d), combined with the higher exsolved volatile fraction (Figure 7c), mean that melt rise speeds would have been significantly slower for a given conduit width and pressure gradient, favoring bubble rise and coalescence (e.g., Parfitt & Wilson, 1995). Thus, the high explosivity observed at F17 is likely a direct consequence of extensive fractional crystallization driving up melt viscosity and H 2 O contents, causing a larger volume of bubbles to exsolve (and coalescence) at greater depths in the conduit at a deeper depth, along with likely feedbacks between a higher viscosity and lower ascent rates in the dacitic melt (further favoring bubble coalescence).

Identifying Parental Magmas
Major and volatile element systematics indicate that melts erupted in the first 2 weeks of the 2018 eruption underwent extensive (and variable amounts) of fractional crystallization during rift-zone storage. However, using major elements systematics alone it is difficult to determine whether these melts were derived from fractionation of a single parent magma, or a number of different magma bodies, because the major element contents of different Kīlauea eruptions are relatively similar at a given MgO content (Helz & Wright, 1992). In contrast, through time, Kīlauea eruptions show clear variations in incompatible trace element ratios such as Nb/Y, La/Yb, and Zr/Y (Figure 4b), interpreted to represent heterogeneity in the Hawaiian mantle plume, as well as variations in the degree of mantle melting (Pietruszka & Garcia, 1999a, 1999b. It is probable that the parent magma(s) for the variably evolved melts erupted in 2018 were intruded into the rift zone and erupted at the surface during a previous eruptive episode (the most recent eruptive events occurred in 1790, 1823, 1840, 1955, and 1960. The 1790 and early 1955 fissures are of particular interest, because these fissures were located in the same area of the LERZ as the 2018 eruption (Figure 1c; Helz, 2008;Moore, 1992;Wright & Fiske, 1971). Similarly, the 1960 eruption occurred a short distance downrift, so presumably its magma supply passed through the section of the rift zone where the 2018 eruption occurred. Here, we use trace element ratios to determine the composition of historic eruptions on the LERZ and compare these to our measurements from 2018 eruptive products to deduce the most probable parent magma(s).
Unlike previous work using trace elements to track magma batches at Kīlauea in lavas which have undergone mostly olivine crystallization (Pietruszka & Garcia, 1999b;Wieser et al., 2019), the melts analyzed here have experienced fractionation of a number of phases (olivine, clinopyroxene, orthopyroxene, plagioclase, ± ilmenite, magnetite, and apatite). The presence of these other phases in addition to olivine provides a potential mechanism for these trace element ratios to be fractionated during crystallization, which might obscure attempts to compare rift-stored melt compositions to magma batches defined by analyses in the literature of mostly olivine-saturated samples from previous LERZ eruptions and Kīlauea's summit. To assess which element ratios are least fractionated by extensive crystallization at Kīlauea, we use the trace element systematics of glasses from drill cores taken through the Kīlauea Iki lava lake (Greaney et al., 2017). This lava lake experienced extensive fractionation of a single magma batch across a range of MgO contents very similar to the lavas examined here, with a very similar mineral assemblage (Helz, 1980;Helz & Thornber, 1987). Thus, evaluation of changes in trace element contents with decreasing MgO in these lava lake samples provide a unique opportunity to establish whether trace element ratios are fractionated across this differentiation interval. Between 6 and ∼0 wt% MgO, Kīlauea Iki matrix glass Nb, Yb, and Y all increase by a factor of 2-3, while La and Zr behave more incompatibly; increasing by a factor of 3-4 times ( Figure S16 in Supporting Information S1). This causes the Nb/Y ratio to remain almost constant during extensive fractionation, while La/Yb increases slightly (∼20%), and Zr/Y more than doubles ( Figure S17 in Supporting Information S1). Pietruszka et al. (2021) also demonstrate that Nb/Y is resistant to fractionation using available partitioning data for relevant phases. In addition to being the least fractionated ratio, it is also advantageous to use Nb/Y for comparisons because it is the only trace element ratio which has been reported by a number of different studies investigating historic lavas on the LERZ (because both Nb and Y can be measured accurately and precisely by XRF, unlike La and Yb).
Previous work has shown that Nb/Y ratios in historic summit and rift eruptions follow a prominent trough-peaktrough shape between 1800 to present Marske, 2010;Pietruszka & Garcia, 1999a, 1999b. The highest Nb/Y ratios are observed in lavas from the mid 1900s, declining toward the values measured during the final stages of the Puʻuʻōʻō eruption (Figure 4b; Gansecki et al., 2019). Nb/Y ratios from whole rock analyses of previous LERZ eruptions follow a similar trend to summit eruptions (Helz & Wright, 1992;Marske, 2010;Norman, 2005). Matrix glasses erupted in Early Phase 1 in 2018 have much more similar Nb/Y ratios to the 1955-1960 eruptions than melts erupted in the period between 1790 and 1840 CE (characterized by lower Nb/Y ratios; ∼0.5-0.6, see also Pietruszka et al., 2021). This indicates that melts erupted in 1955 or 1960 are a more probable parent than those erupted in 1790 and 1840.
Two main phases of the 1955 eruption have been identified, with an early phase (1955E) thought to tap more MgO-poor magmas already present within the rift zone (5-5.7 wt % MgO Whole-rock ), while later phase (1995L) tapped more MgO-rich magmas (6.2-6.8 wt% MgO Whole-rock ) interpreted to reflect mixing of the 1955E lava with a more MgO-rich, summit-derived component (Helz & Wright, 1992). This is analogous to the temporal evolution from MgO-poor to MgO-rich observed in the 2018 eruption. Using the analyses of Marske (2010), along with additional analyses of 1955 lavas from Helz and Wright (1992) and Norman (2005), we show that the mean and range of Nb/Y ratios in lavas erupted in the early and late phase of 1955 are a good match to the matrix glass compositions erupted in Early Phase 1 in 2018 (Figure 4a). Using Nb/Y ratios alone, it is not possible to definitively distinguish between the 1955E, 1955L, and 1960 lavas, but it is notable that the 1955E vents are geographically extremely close to the 2018 vents (Figure 1c).
Whole-rock Nb/Y ratios (pink dots on Figure 4a; Pietruszka et al., 2021) of lavas erupted in the first 20 days of the 2018 eruption form a mixing trend between the "1955-1960 like" Nb/Y ratios exhibited by Early Phase 1 lavas, and the lower Nb/Y ratios measured in Phase 3. This trend is also seen in our matrix glass analyses of F13-react and F20. These mixing trends indicate that LERZ-stored melts were progressively flushed out by the more MgOrich, lower Nb/Y component supplying Phase 3 (Gansecki et al., 2019). Interestingly, while some of the more MgO-rich samples erupted at F17 (WR MgO = 3-4 wt%, Gansecki et al., 2019;Pietruszka et al., 2021) show evidence for mixing with the lower Nb/Y component from Phase 3 ( Figure S18 in Supporting Information S1), our sample (glass MgO = ∼0.5 wt%, WR MgO = 2.3-2.6 wt%) has glass Nb/Y ratios which entirely overlap with these Early Phase 1 glasses (Figure 4a). This overlap indicates that melts erupted in Early Phase 1 and at F17 were likely derived from a common parent magma, and that our F17 sample experienced minimal mixing with magma erupted in Phase 3. Thus, it seems likely that the initial activity at F17 was triggered by an increase in overpressure, change in crustal stress, or thermal rejuvenation resulting from the intrusion of a dike containing Phase 3 magma into this section of the rift zone, rather than as a direct consequence of mixing between different magma compositions. Subsequent mixing of magmas indicated by whole-rock Nb/Y data from Gansecki et al. (2019) and Pietruszka et al. (2021; Figure S18 in Supporting Information S1), likely helped to liberate more of this stored material, fueling the large lava flows with more MgO-rich compositions that erupted from F17 until 25 May (not sampled in this study).
Plagioclase-hosted melt inclusion Nb/Y ratios show no correlation with the anorthite (An) content of the host (Figure 4c), which is consistent with a common parental magma for all rift-stored melts. Similarly, melt inclusion Nb/Y ratios in olivine, clinopyroxene and orthopyroxene from early fissures show no correlation with host Mg# (Mg/[Mg + Fe] atomic; Figure 4d), and are clearly distinct from the significantly lower Nb/Y ratios in matrix glasses and the majority of melt inclusions in Phase 3. This indicates that very few inclusions trapped melts that had undergone significant mixing with Phase 3 melts. Interestingly, three plagioclase-hosted melt inclusions from early F8, F5, and F20 plot with F17 melt inclusions in Nb/Y-An space, and 6 inclusions from F17 have higher MgO contents and major element contents indicating post-entrapment cooling. This supports the idea that more MgO-poor pockets of melt exist in close proximity with more MgO-rich bodies in the rift zone magma body, in order for more MgO-rich fissures to scavenge more An-poor crystals, and for more MgO-poor fissures to scavenge more An-rich crystals. The five olivine-hosted melt inclusions erupted on 28 May in Phase 3 (F8) with higher Nb/Y ratios (∼0.7; Figure 4d) are hosted in olivines with some of the lowest forsterite contents of Phase 3, and have melt inclusion TiO 2 contents overlapping with Early Phase 1 melt inclusions and glasses ( Figure S19 in Supporting Information S1). This indicates that these 5 crystals may have also grown from rift-stored melts with higher Nb/Y ratios, and were then recycled into the lower Nb/Y Phase 3 melts.
Overall, melt inclusion Nb/Y ratios indicate that melt inclusions erupted in the first 2 weeks of 2018 likely crystallized from a single parent magma body (or a series of magma bodies with similar parental magmas), and that very few melt inclusions crystallized following mixing between rift-stored and Phase 3 melts (although crystal recycling between different stored magma batches did occur).

Pre-Eruptive Magma Storage Depths
Melt inclusion H 2 O and CO 2 contents provide further constraints on the nature of magma storage within the LERZ. The increase in H 2 O and decrease in CO 2 with decreasing MgO content (Figures 3c and 3d) provide strong evidence that melts were volatile-saturated throughout the fractionation interval considered here. Thus, estimates of the melt inclusion temperature, major elements, H 2 O and CO 2 contents at the time of entrapment may be used to calculate the pressure at which those dissolved volatile contents would be saturated, and, by extension, the pressure at which the melt inclusion was trapped within its host crystal. This is termed the "saturation pressure," or "entrapment pressure." By estimating the crustal density profile, these saturation pressures can be converted into depths within the volcanic edifice.
As discussed in Section 4.2, H 2 O-CO 2 systematics indicate that a number of melt inclusions were likely trapped during ascent, after degassing had begun. While these still yield useful saturation pressures (indicating the depths at which inclusions are forming on the way to the surface), we also filter the data set to only consider melt inclusions trapped from magmas prior to the onset of degassing on ascent to identify magma storage reservoirs.
To remove melt inclusions trapped during degassing, we exclude inclusions with CO 2 contents lying more than 30% below the CO 2 content predicted by MELTS model B at the MgO content of each melt inclusion. To remove inclusions which have undergone H 2 O degassing or diffusive re-equilibration, we filter inclusions more than ±30% from MELTS model B. An example of this filtering process is shown in Figures 6d-6f for F17 melt inclusions. Melt inclusions meeting both criteria are marked with red crosses (see also Figure S20 in Supporting Information S1).
We use the solubility model MagmaSat (Ghiorso & Gualda, 2015) implemented in the Python3 tool VESIcal (v.0.01; Iacovino et al., 2021) to calculate entrapment pressures from PEC-corrected major element, H 2 O, and CO 2 contents of each melt inclusion, and a temperature estimated using the Helz and Thornber (1987) MgO thermometer. We choose MagmaSat because this model provides the best fit of available models to experimental data with basaltic to dacitic compositions, as a result of its thermodynamic nature and extensive calibration data set (Wieser et al., 2022). Additionally, use of this model provides consistency with our models of volatile solubility during fractionation and ascent, as MagmaSat is the solubility model used in Rhyolite-MELTS v.1.2.0.
Filtered early Phase 1 melt inclusions show a distinct clustering of saturation pressures at ∼0.48-0.8 kbar, with four inclusions with pressures up to 1.15 kbars (Figure 8a). Filtered F17 melt inclusions show a remarkable overlap with Early Phase 1, mostly spanning 0.55-0.8 kbar, with two higher pressure inclusions (Figure 8b). Only two melt inclusions from F13-react and five from F20 pass our filters (Figures 6a, 6c and 8c-8d, Figure S20 in Supporting Information S1), but these also yield very similar pressure distributions. Unfiltered data for all samples stretches to significantly lower pressures, indicating that a number of melt inclusions were likely trapped during ascent (Figures 8e-8h). It is worth noting that melt inclusions from F13-react and F20 samples show clusters at 0.01-0.03 wt% CO 2 , which could indicate that a large number of inclusions were trapped during degassing, or that these fissures were tapping a slightly shallower magma reservoir.  (Ghiorso & Gualda, 2015) implemented through VESIcal . (a-d) shows filtered melt inclusions (within ±30% of MELTS model B for H 2 O, and −30% for CO 2 , see Figures 6e and 6f). Panels (e and h) shows all melt inclusions. Each gray line on the violin plot represents a melt inclusion, and the width of the "violin" represents the clustering of data at each pressure. On the RHS, the individual saturation depths for each inclusion are shown with symbol shapes and colors corresponding to the host phase. The depth at which the dacitic magma was drilled at PGV is shown as a yellow star, and the depth of the low V p /V s region (V p /V s = 1.66 vs. 1.78 background) imaged by P. Cooper and Dustman (1995) is shown in yellow.
Depth ( Using a crustal density of ρ = 2,600 kg/m 3 based on the relationship between drill depth and pressure presented from borehole drilling results at Puna Geothermal Venture from Teplow et al. (2009), melt inclusion saturation pressures can be converted into magma storage depths. Filtered melt inclusions cluster between ∼2 and 3 km (Figure 8), which aligns remarkably well with the depth at which an even more silicic dacitic magma was drilled at the PGV injection well KS-13 in 2005 (depth of 2,488 m, ∼650 bars pressure, yellow star on Figure 8). There is also an interesting correspondence between our saturation depths and the depth of a low V p /V s anomaly imaged in this area between Puʻuhonuaʻula and Puʻulena Crater (Figure 1b, ∼2km depth, vertical extent of 1.5 km, P. Cooper & Dustman, 1995). This anomaly has been interpreted as a geothermal reservoir with steam-filled cracks. The lateral dimensions of this body (2.5 km parallel to the rift zone, 2 km perpendicular to the rift zone) means that it directly underlies the fissures erupted in the first 2 weeks of the 2018 eruption. It seems highly probable that the magma body being tapped by the 2018 eruption lies directly beneath this geothermal reservoir, providing a heat source.

Temporal Evolution of Kīlauea's LERZ
We present a schematic model summarizing the evolution of Kīlauea's LERZ, based on our findings and previous work on historic lava samples (Figure 9). The early fissures of the 1955 eruption, which are located very close to the 2018 eruption site (Figure 1c), tapped melts with whole-rock contents of ∼5-5.7 wt% MgO. These melts were likely resident in the rift zone since at least since 1924 (when a dike swarm propagated past the 2018 eruption site into Kapoho) or 1790 (the last eruption where magma breached the surface in this section of the LERZ; Figure 1c). Alternatively, 230 Th-226 Ra dating of 1955E lavas suggests that clinopyroxenes and plagioclases were resident for at least 550 years prior to eruption. The later phases of the 1955 eruption focused at fissures uprift of the 2018 eruption site tapped more MgO-rich magmas (6.2-6.8 wt% MgO Whole-rock ) which have been interpreted to reflect renewed supply from the summit reservoir (Helz & Wright, 1992; Figure 9a). This resupply may have propagated downrift, helping to replenish reservoirs depleted during the early phases of the 1955 eruption. In 1960, activity commenced downrift of the 1955E fissures, first trapping stored magma with chemical similarities to the 1955 lavas (6.1-6.6 wt% MgO Whole-rock ), followed by the appearance of progressively higher MgO contents associated with supply of material from Kīlauea's summit reservoir (7-13 wt% MgO Whole-rock , Russell & Stanley, 1990). As in the later phases of the 1955 eruption, the 1960 summit resupply may also have replenished the larger reservoir complex residing beneath the 2018 eruption site (as the 1960 magmas are very similar in terms of Nb/Y and major elements; Russell & Stanley, 1990, Figure 9b). This pulse of activity in 1955-1960 left substantial amounts of magma at depth, which continued to cool and fractionate over the next 60 years as activity shifted to Kīlauea's summit, South West Rift Zone and Middle East Rift Zone. In 2005, a dacitic pocket of magma was encountered during geothermal drilling (Figure 9c).
The propagation of a dike downrift from Puʻuʻōʻō on 30 April 2018 disturbed these rift-stored melts, with magma erupting at the surface beginning on 3 May (Figure 9d). After a short pause after 9 May, erupted melt compositions between 12 and 18 May begin to show mixing between LERZ-stored and dike-supplied melts (e.g., in matrix glass Nb/Y ratios, Figure 4a, Gansecki et al., 2019, Figure S18 in Supporting Information S1), although erupted crystals are still dominated by geochemical signatures indicative of the LERZ-stored component (Figures 4c and 4d). Dacitic melts also began erupting at F17 and Fissure 13 shows mixing between these dacitic melts and the more MgO-rich LERZ-stored component. By the time activity focused at Ahuʻailāʻau (F8) on 28 May, almost all LERZ-stored crystals (and melts) had been flushed out by dike stored material (Figures 4d and 9e, Figure S19 in Supporting Information S1).
It is noteworthy that Early Phase 1 lavas have similar MgO contents of early 1955 lavas, so if the 1955E lavas are the parent, the portion of the reservoir tapped by Early Phase 1 fissures must have been sufficiently thermally buffered that minimal fractionation took place over ∼60 years. Relatively slow cooling rates within a thermally buffered magma body are consistent with ∼500 years residence times inferred for the early 1955 lavas by K. M. Cooper et al. (2001). Alternatively, the 1955E reservoir may have been topped up by more MgO-rich magma in the later phases of the 1955 eruption or the 1960 eruption (Figure 9b).
The larger amount of fractional crystallization experienced by F17 melts erupting on 14 May and the dacitic magma encountered during drilling in 2005 compared to melts erupted between 3 and 9 May is intriguing. The similarity of melt inclusion saturation pressures and drilling depths between all three samples, and similarity of Nb/Y ratios between Early Phase 1 and F17 from these samples indicate that the melts supplying all three bodies had a similar parent magma, and were stored at a similar depth. We suggest that the F17 dacitic melt may formed on the periphery of the larger magma body located in this region, allowing more cooling to occur in a specific amount of time (Figure 9d). This interpretation is supported by the fact that F17 is offset ∼220 m NE from the trend defined by other fissures. Pockets of dacitic melt may also have formed within regions of enhanced hydrothermal cooling. Deconvolving the nature of the evolution of this complex LERZ magma body in space and time will likely require detailed thermal modeling, and isotopic data to provide a more precise way to identify different magma batches with more resolution than can be achieved using incompatible trace element ratios.

Global Occurrence of More Silicic Melts in Basalt-Dominated Volcanic Region
Analyses of erupted matrix glasses and melt inclusions from the 2018 eruption demonstrate that extensive crystallization of an ocean island basalt causes substantial enrichment in incompatible volatile and trace element species (Cl, F, and H 2 O). In particular, the increase in melt H 2 O contents combined with a rapid increase in melt viscosity following Fe-Ti oxide fractionation means that these silicic magmas have the potential to display more explosive eruption styles relative to the more typical Hawaiian style exhibited by basaltic to basaltic-andesite magmas.
The eruption of dacitic melt in 2018 at Kīlauea provides support to the growing volume of literature indicating that andesitic-rhyolitic magmas are more common than we thought in predominantly basalt-dominated volcanic regions. For example, Stock et al. (2020) report low An plagioclase crystals (down to An 48 ) in basaltic lava and reticulite dominated by a higher An plagioclase population (An 78−82 ) from the 2015 eruption of Volcan Wolf in the Galápagos Archipelego. They also identify low An plagioclase crystals in gabbroic nodules from the 1968 Fernandina eruption. Fractionation models indicate that these plagioclase crystals must have grown from basaltic trachy-andesite and trachy-andesitic melts. Additionally, resorbed quartz has been identified in tephra from the 2015 eruption of Wolf. This indicates that the subvolcanic plumbing system at these two volcanoes which erupt monotonous basaltic compositions must contain highly fractionated, silicic melts at depth not yet observed at the surface.
Rhyolites comprise approximately 10% of erupted lavas in Iceland (Jónasson, 2007), and have been noted in a variety of settings among mid-oceanic ridges (Wanless et al., 2010). While the continuum of melt inclusion major element contents from basalt to dacite, and good fit of major, volatile and trace element contents to MELTS fractional crystallization models indicate that Kīlauea dacites result from extensive fractional crystallization, silicic magma occurrences in Iceland and MOR have been variably attributed to fractional crystallization, crystal-melt segregation, and partial melting of hydrothermally altered crust based on isotopic, volatile, major and trace element data (Elders et al., 2011;Geist et al., 2021;Masotta et al., 2018;Rooyakkers, Stix, Berlo, Petrelli, Hampton, et al., 2021;Wanless et al., 2010).
Interestingly, Rooyakkers, Stix, Berlo, Petrelli, and Sigmundsson (2021) highlight the fact that silicic magmas have been encountered four times in the last two decades years during borehole drilling: dacite at Kīlauea (Teplow et al., 2009), trachyte in Menegai, Kenya (Mbia et al., 2015), and rhyolite on two occasions at Krafla (Elders et al., 2011;Mortensen et al., 2010). While the Menegai and Krafla drilling encounters were a little less surprising in terms of the composition encountered compared to Kīlauea (because rhyolitic and syenitic magmas are better represented in surface deposits at these locations), It is noteworthy that all of these encounters intercepted silicic magmas at very similar depths (2,488 m at Kīlauea, ∼2,000 m in Kenya, 2,571 m and 2,104 m at Krafla), all in areas of active hydrothermal exploration. Although these silicic magmas may be generated by different mechanisms, all require high heat fluxes at shallow levels, either to melt/rejuvenate existing material, or to allow fractional crystallization to proceed without the magma simply freezing in the crust. At Kīlauea, this high heat flow likely results from the focus of repeated intrusion over hundreds of years in this particular section of the rift zone, perhaps resulting from a "jog" with two offset branches creating a local extension stress field (Kenedi et al., 2010). This is somewhat analogous to the higher focusing of melt, and therefore higher occurrence of silicic magma in overlapping spreading centers of mid-oceanic ridges (Kent et al., 2000).
While 2018 marked the first documented occurrence of dacitic melt erupting at the surface, given that 70% of Kīlauea's surface is <500 years old, and 90% is <1,100 years old (Holcomb, 1987), these melts could erupt on a centurial-millennial basis, and still be absent in surface exposures. This is particularly true given that dacitic melts are most likely to form in parts of the volcano characterized by high intrusion rates, to maintain the high heat fluxes at shallow levels required for melts to fractionate to dacitic compositions without freezing in the crust. In turn, high intrusion rates mean these sites are most likely to be rapidly resurfaced by volcanic activity. For example, 75% of the surface of the Lower-East Rift Zone is <500 years old (Moore, 1992). If magmas with silicic compositions were able to evade sampling at Kīlauea, one of the world's best-studied ocean island volcanoes, prior to drilling in 2005, it seems highly likely that these compositions are present at other ocean island volcanoes. Specifically, based on the need for high heat flow to generate silicic magma bodies, and the increasing number of drilling encounters, the anomalously hot areas of ocean islands are not only the most likely to contain more silicic melts at shallow depths, but will also be the preferred targets for hydrothermal exploration.
When the dacitic body was drilled at Kīlauea in 2005 at 2,488 m depth, it was noted that the glass did not show vesiculation (Teplow et al., 2009), leading to inferences that the magma was volatile-poor. Given that our melt inclusion measurements demonstrate that F17 melts were H 2 O-rich, we instead hypothesize that the low vesicularity resulted from the fact that a relatively small mass of volatiles is exsolved during fractionation from a basalt to dacite at this depth, meaning there is limited potential for a large, gas-rich cap to develop on top of these bodies (Figure 5a). Upon ascent, significant quantities of H 2 O only begin to exsolve from a dacitic melt composition at ∼300 bars. Thus, as long as the drilling fluid quenches these silicic melts at higher confining pressures, there would be limited potential for large amounts of volatile exsolution and catastrophic well blowout (Teplow et al., 2009). However, if allowed to rise to shallow levels prior to quenching, silicic melts could present a hazard to geothermal drilling in basalt-dominated settings.
In addition to being encountered during drilling, it appears that the eruption of silicic melts at the surface may be closely coupled to the injection of basaltic dikes into the region. At Kīlauea, based on the limited evidence for mixing between the dacite and the dike-supplied magma in the initial products of F17, we suggest that dike injection increased overpressure or changed the stress state in the magma body, causing the eruption. Similarly, Rooyakkers, Stix, Berlo, Petrelli, and Sigmundsson (2021) examine the rhyolitic magma producing the Víti maar at Krafla, Iceland, and conclude that these melts show very little textural or chemical evidence for mingling/ mixing with the basaltic "trigger." This contrasts with explosive eruptions like Askja, 1875 which show clear evidence textural evidence for pre-eruptive magma mixing (Sparks et al., 1977). Clearly, the potential of silicic magma bodies to generate explosive eruptions should be considered during hazard assessments in a wide variety of ocean islands and other basalt-dominated settings (Rooyakkers, Stix, Berlo, Petrelli, & Sigmundsson, 2021), both in terms of natural eruptions initiated by a dike, and drilling encounters.

Conclusions
Basaltic to andesitic magmas erupted during the first 2 weeks of the 2018 eruption (3-16 May) were likely formed from variable amounts of crystallization of a magma body at ∼2-3 km depth within the LERZ. The similarity in Nb/Y ratios in melt inclusions from these early eruptive fissures and a range of host crystal chemistry implies that crystallization proceeded from a single magma parent body (or a series of bodies with near-identical trace element chemistry). Comparison of Nb/Y ratios to historic LERZ lavas indicates that melts erupted in 1955-1960 are the most likely parent to the early 2018 lavas (see also Helz, 2008;Pietruszka et al., 2021;Teplow et al., 2009). Extensive crystallization of a section of this larger magma body (perhaps on the periphery or in a region with enhanced hydrothermal cooling) produced a dacitic melt composition highly enriched in incompatible elements such as Cl, F, Zr, and H 2 O. Combined with an increase in magma viscosity with increasing SiO 2 content and dropping temperatures, this H 2 O-enrichment accounts for the explosive strombolian behavior exhibited by the eruptive fissure tapping this melt (F17) without requiring external sources of volatiles such as groundwater. Although the high viscosities of these dacitic melts mean they are unlikely to erupt spontaneously, the 2018 eruption shows that they may be triggered by the injection of basaltic intrusion. The volatile-rich nature of these bodies should also be considered when prospecting for geothermal wells; if the melt erupted at F17 had been drilled and allowed to depressurize to ∼300 bars prior to quenching, large amounts of vesiculation may occur.