Magmatic Storage and Volatile Fluxes of the 2021 La Palma Eruption

The 2021 La Palma eruption (Tajogaite) was unprecedented in magnitude, duration, and degree of monitoring compared to historical volcanism on La Palma. Here, we provide data on melt inclusions in samples from the beginning and end of the eruption to compare the utility of both melt and fluid inclusions as recorders of magma storage. We also investigated compositional heterogeneities within the magmatic plumbing system. We found two populations of olivine crystals: a low Mg# (78–82) population present at the beginning and end of eruption, recording the maximum volatile contents (2.5 wt % H2O, 1,800 ppm F, 700 ppm Cl, 3,800 ppm S) and a higher Mg# (83–86) population sampled toward the end of the eruption, with lower volatile contents. Despite their host composition, melt inclusions share the same maximum range of CO2 concentrations (1.2–1.4 wt %), indicating olivine growth and inclusion capture at similar depths. Overall, both melt and fluid inclusions record similar pressures (450–850 MPa, ∼15–30 km), and when hosted in the same olivine crystal pressures are indistinguishable within error. At these mantle pressures, CO2 is expected to be an exsolved phase explaining the similar range of CO2 between the two samples, but other volatile species (F, Cl, S) behave incompatibly, and thus, the increase between the two olivine populations can be explained by fractional crystallization prior to eruption. Finally, based on our new data, we provide estimates on the total volatile emission of the eruption.


Introduction
The Canary Islands are part of the Canary hotspot track of the Central Eastern Atlantic and consist of seven main volcanic islands stretched over 500-km, sitting on Jurassic-aged oceanic crust (Carracdeo & Troll, 2016).These volcanic islands record a general westward age progression with La Palma and El Hierro subaerial volcanism being the youngest at ∼1.7 and 1.1 Ma, respectively (Carracedo et al., 2001).Multiple models of formation have been argued as an explanation for the existence of the Canary Islands, including a mantle plume, edge-driven convection, or a tectonically controlled fracture (Anguita & Hernán, 2000;Carracedo et al., 1998;Duggen et al., 2009;Geldmacher et al., 2005;King & Ritsema, 2000;Zaczek et al., 2015).In recent years, a combination of a mantle plume edge-driven convection toward the "cold" African craton explains the apparent weak age progression of the Canary Islands as well as historic volcanism on both sides of the archipelago (Manjón-Cabeza Córdoba & Ballmer, 2022;Negredo et al., 2022;Troll & Carracedo, 2016).Although isotope systematics suggest the presence of recycled materials in the Canary Islands source (Day & Hilton, 2020;Day et al., 2009Day et al., , 2010Day et al., , 2022;;Thirlwall et al., 1997;Walowski et al., 2019), the role of recycled volatiles, specifically CO 2 , has been underinvestigated as a potential driver of eruptive activity (Aiuppa et al., 2021;Allison et al., 2021;Taracsák et al., 2019).

The 2021 La Palma (Tajogaite) Eruption
The 2021 Tajogaite-eruption on La Palma displays many characteristics of a typical Canary historical eruption mentioned above.Prior to the start of volcanic activity on 19 September 2021, unrest under the CVR had been detected following seismic swarms in 2017 and 2018 (Torres-González et al., 2020).An increase in helium isotopic ratios and CO 2 emissions indicated a mantle-sourced magmatic intrusion, placed at 25 km 4 years prior to the start of eruptive activity (Torres-González et al., 2020).Following the start of seismic unrest in 2017, La Palma experienced nine seismic swarms prior to the final seismic event which began 8 days before the eruption start (Fernandez et al., 2022;Longpré, 2021).Pre-eruptive activity starting on 11 September has been interpreted as magmatic injection or underplating in a sill-like body at 7-10 km, which then transitioned to vertical dike motion as the path was cleared for vent opening on 19 September (D'Auria et al., 2022;Del Fresno et al., 2023;De Luca et al., 2022).Propagation of both sills and dikes led to the shallowing of earthquake hypocenters immediately prior to eruption (as also observed in the 2011-2012 Tagoro eruption off El Hierro) (López et al., 2012;Martí et al., 2013).Dike complexes, like those which fed the 2021 eruption, are exposed in eroded sections of the rift zones and are the main mechanisms of oceanic island growth (Carracedo, 1994;Carracedo et al., 2011;Delcamp et al., 2012).
When compared to the last two previous on-island eruptions of Cumbre Vieja in 1949 (378 days, volume of 44 × 106 m 3 ) and 1971 (245 days, volume of 46 × 106 m 3 ), the 2021 eruption was twice as long (85 days) and emitted an order of magnitude more material with a total volume of 253.9 ± 6.9 × 106 m 3 (Civico et al., 2022) (Figure 1c).Ash emission from the eruption was significant and long-lasting, reaching all Canary Islands, mainland Europe, and even Barbados (Gaston et al., 2022).Plume heights reached up to 8.5 km with an average of 3.5 km above sea level visible from Tenerife (Carracedo et al., 2022).Explosivity of the 2021 eruption was also significant and variable, with multiple vents active simultaneously each displaying various behavior ranging from pulsating strombolian blasts to effusive lava flows (Bonadonna et al., 2022;Carracedo et al., 2022).
The 2021 lavas exhibit a compositional hybrid nature.Initially, they displayed a tephritic composition characterized by low MgO (∼6 wt%) and TiO 2 (∼4 wt%).As the eruption progressed, the lavas became increasingly mafic, reaching a consistently basanitic composition with higher MgO (∼8 wt%) and lower TiO 2 (3.7 wt%) around day 20 (Day et al., 2022).Recent investigations into Sr and Os isotopes have identified the tapping of a more primitive reservoir during progression of the eruption (Day et al., 2022;Gonzalez-Garcia et al., 2023;Ubide et al., 2023).The process of the initial eruption involving a more fractionated magma is not unique to the 2021 event.In both the 1949 and 1971 eruptions, lower-MgO amphibole-bearing products erupted before basanites, suggesting fractionation as a widespread process beneath the Cumbre Vieja Ridge (Barker et al., 2015;Klügel et al., 2000).Although exhibiting temporal and compositional variability, a homogeneous asthenospheric magma source with a constant degree of 2.5%-3% partial melting has been proposed as the primary feeding source for the eruption, with modification by fractional crystallization (Day et al., 2022;Gonzalez-Garcia et al., 2023).
Because of the explosive activity and deposition of substantial quantities of tephra containing large, euhedral olivine phenocrysts, the 2021 eruption provides ideal conditions for studying melt inclusions due to the small size of tephra clasts, and, therefore, fast quench time minimizing post-eruptive modification (Lloyd et al., 2012).Melt inclusions, tiny droplets of melt trapped during crystal growth, are a powerful tool to understand the drivers and evolution of magmatic eruptions (Allison et al., 2021;Cashman et al., 2017;DeVitre, Gazel, et al., 2023;Hauri, 2002;Kent, 2008;Lowenstern, 1995;Metrich & Wallace, 2008;Moore et al., 2015Moore et al., , 2021;;Rasmussen et al., 2018;Roedder, 1979;Wieser et al., 2021) (Figure 2).Though ideally able to preserve valuable information about magmatic processes prior to eruption, melt inclusions are subject to several, simultaneously occurring, postentrapment processes.These processes include olivine crystallization on melt inclusion walls and subsequent Mg-Fe diffusive re-equilibration, also known as post-entrapment crystallization (PEC) (Danyushevsky et al., 2000;Esposito et al., 2012).H 2 O loss from trapped melt upon ascent and clast-cooling (Lloyd et al., 2012), melt inclusion vapor bubble formation (Moore et al., 2015;Rasmussen et al., 2020;Wallace et al., 2015), and carbonate crystallization on vapor bubble walls (Moore et al., 2015(Moore et al., , 2021;;Venugopal et al., 2020) are also common post-entrapment processes.Additionally, analytical challenges can also obscure the true volatile contents of melt inclusions.Uncalibrated Raman analysis of melt inclusion vapor bubbles within the CO 2 miscibility gap (ρCO 2 = 0.2-0.7 g/cm 3 ) below the critical point of CO 2 (∼31.1°C)can lead to underestimates of CO 2 density in melt inclusion vapor bubbles leading to errors in saturation pressure calculations of up to 6 km (DeVitre, Dayton, et al., 2023).Here, we use the most recent best practices of melt-inclusion analysis to ensure the most accurate estimates of saturation pressure by hand-picking single crystals from fine-grained tephra and avoiding melt inclusions with anomalously high bubble-volume fractions that indicate obvious heterogenous entrapment (Steele-MacInnis et al., 2017).We also avoided melt inclusions with carbonates on the bubble walls and measured melt inclusions in the third dimension using Raman spectroscopy.Lastly, we analyzed melt-inclusion vapor bubbles above the critical point of CO 2 using calibrated Raman spectroscopy (DeVitre, Dayton, et al., 2023;DeVitre, Gazel, et al., 2023;DeVitre et al., 2021;Hanyu et al., 2020;Mironov et al., 2020;Rose-Koga et al., 2021).
From our previous work using fluid-inclusion barometry (Dayton et al., 2023), we constrained the magma storage at mantle depths prior to the 2021 eruption (15-27 km).In an attempt to comprehensively examine the geochemical processes and investigate the drivers of volcanism in the Canary Islands, we have now analyzed melt inclusions (glass plus vapor bubbles) and matrix glasses from the start (September, LM6) and end (December, LM0) of the 2021 eruption.We obtain volatile budgets of the mantle source that feeds the 2021 eruption and further constrain the magmatic plumbing system of Cumbre Vieja by identifying heterogeneities in volatile concentrations during the eruptive event.Specifically, we report volatile contents (H 2 O, CO 2 , fluorine, chlorine, sulfur) and major and trace element data on melt inclusions as well as whole rock data on the stratigraphic section used in Dayton et al. (2023) (Figure 1d).We analyze multiple melt inclusions (up to six) within the same olivine crystal (e.g., Figure 2b) to elucidate heterogenous entrapment (cf.Steele-MacInnis et al., 2017) of exsolved CO 2 fluid during melt inclusion capture, which can result in erroneously high CO 2 contents and saturation pressures if not properly identified.We compare our results to previously published fluid inclusion data from the same samples (Dayton et al., 2023), and analyzed both melt and fluid inclusions within the same crystals (Figures 2c and 2d) to constrain magmatic processes and directly compare melt and fluid inclusion to investigate their utility and reliability as recorders of magma storage location.

Materials and Methods
Samples were collected from a stratigraphic section containing the entirety of 2021 eruptive deposit (Figure 1d).Sample LM6 represents the first 3 days of eruptive activity (September 19 to ∼21) (Romero et al., 2022), while LM0 represents the final eruptive unit (∼December).Olivine crystals were hand-picked, and individually polished to identify melt and fluid inclusions.CO 2 densities of melt inclusion vapor bubbles and fluid inclusions were collected via Raman spectroscopy in the Department of Earth and Atmospheric Sciences at Cornell University.Volatiles in exposed glass of melt inclusions were analyzed by secondary ion mass spectrometry (SIMS) in the Northeast National Ion Microprobe Facility at Woods Hole Oceanographic Institution, followed by major elements in the glass and host olivine by electron microprobe analyzer (EMPA) at the University of Massachusetts Amherst in the Department of Geosciences.Finally, both the glass and trace elements in olivine were collected by laser ablation induced coupled plasma mass spectrometry (LA-ICP-MS) in the Department of Earth and Atmospheric Sciences at Cornell University.Whole rock tephra data was collected via X-ray fluorescence (XRF) and induced coupled plasma mass spectrometry (ICP-MS) in the Peter Hooper GeoAnalytical Lab at Washington State University.Details on converting melt and fluid inclusion pressures are discussed at length in Supporting Information S1 and Data Set S1, along with more details on all materials and methods.To calculate volatile fluxes, we use an estimates of the total emitted volume of volcanic material, including subaerial deposits and proximal fallout plus the volcanic edifice (0.2539 ± 0.0069 km 3 ), constrained through high-resolution surface modeling of the eruption (Civico et al., 2022) and the volatile contents of melt inclusions and matrix glass, we calculated volatile budgets of SO 2 and CO 2 associated with the 2021 eruption using the petrologic approach (Sharma et al., 2004).As melt inclusion volatile compositions are variable through the eruption, and there are uncertainties in other parameters such as erupted volume, melt crystal fraction and lava density, we utilize normal distributions to assess the spread and error associated with all inputs (Figures S3-S6 in Supporting Information S1, see Data Set S1).We use a density of eruptive material value of 2,403 kg/m 3 calculated from the DensityX (Iacovino & Till, 2019) results of this work on matrix glass compositions (2,700 kg/m 3 ) and the average vesicle modal percent of 11% (Gisbert et al., 2023).We include a density error of 170 kg/m 3 to include the spread of vesicularity values (Gisbert et al., 2023) and DensityX error (Iacovino & Till, 2019).We use a crystal density of 3,010 kg/m 3 calculated from the point counting averages in Table S1 of Ubide et al. (2023).We use a normal distribution, truncated at zero, to obtain a crystal mass fraction of 0.25 ± 0.13, to account for the wide spread of crystallinity estimates within the eruption spanning various grain size fractions (Day et al., 2022;Hornby et al., 2023;Pankhurst et al., 2022;Ubide et al., 2023) (Figure S3 in Supporting Information S1).We utilize a Monte Carlo simulation with 10,000 duplicates to account for uncertainty in vesicularity across all eruptive products.We use an eruption-wide average melt inclusion (3,062 ± 500 ppm) and matrix glass (345 ± 53 ppm) concentrations for SO 2 estimates.A CO 2 budget was calculated using the same methods as for SO 2 with inputs of 12,040 ± 4,553 ppm and 11 ± 7 ppm for melt inclusion and matrix glass concentrations.A HF and HCl budget cannot be calculated as fluorine and chlorine do not degas from the matrix glass in sufficient quantities to be resolvable from analytical noise (Figures 5a  and 5c).We use a normal distribution for all measured volatile contents (melt inclusions and matrix glasses) and truncate any negative CO 2 matrix glass values to zero (Figures S4 and S5 in Supporting Information S1).

Results
Melt inclusions are hosted within olivine crystals with varying core compositions.To simplify, we discuss these two compositional populations based on the M Mg# with one group characterized by lower Mg# (∼78-82) present at both the start and end of the eruption, and a higher Mg# population (∼83-86) present only in the sample at the end of the eruption (LM0) (see Figure 2 of Dayton et al. (2023)).This spread of olivine Mg# is not unique to the 2021 eruption and has been identified in other eruptive suites over the course of La Palma island growth (e.g., Taburiente volcano, Nikogosian et al., 2002).Melt inclusions trapped within the higher Mg# populations yield olivine-liquid thermometry temperatures of ∼1190°C compared to the lower Mg# population with temperatures of ∼1160°C (Putirka, 2008, Equation 22).Regardless of host composition, after post-entrapment crystallization (PEC) corrections (details in Supporting Information S1 and Data Set S1), melt inclusions plot along with lava and tephra whole rock as basanites with ∼5 wt% total alkalis (Figure 3a).The matrix glass, especially at the start of the eruption (LM6) is significantly more fractionated (Figures 3a and 3b) as expected due to continued magmatic evolution and reflecting the abundance of observed microlites in analyzed tephra fragments and ash particles (cf.Hornby et al., 2023).PEC-corrected melt inclusions and whole rock lava span MgO contents of 5.5-8.5 wt% with the most mafic inclusions recorded at the end of the eruption (Figure S2 in Supporting Information S1).The FeOt versus MgO variations suggest a liquid line of descent controlled mostly by olivine crystallization up to around ∼5.5 wt%  Bas et al., 1986) including melt inclusion and matrix glass from samples LM6 (September) and LM0 (December).Tephra (this work) and lava whole rock data (Day et al., 2022) span the entire Tajogaite eruption.(b) CaO/Al 2 O 3 versus SiO 2 diagram identifying the "geochemical kink" in November 2021 and transition from fractionated, more evolved magma representing the start of the eruption to the more primitive composition present through the end of the eruption (Day et al., 2022;Gonzalez-Garcia et al., 2023;Ubide et al., 2023).
MgO, although clinopyroxene and plagioclase have also been observed.Clinopyroxene is also observed as inclusions toward the crystal cores of the lower Mg# olivine population, recording cotectic crystallization during the growth of the most evolved olivine (Figure 2e).Clinopyroxene inclusions are present more infrequently within the lower Mg# rims of the high Mg# olivine crystal population (Figure 2a).The transition between relative fractionated melts recorded at the early stage of the eruption into a more primitive basanite magma (Day et al., 2022;Gonzalez-Garcia et al., 2023;Ubide et al., 2023) is also tracked through melt inclusion CaO/Al 2 O 3 (Figure 3b).The most mafic melt inclusions record higher CaO contents, likely reflecting trapping prior to significant clinopyroxene fractional crystallization.Inclusions of magnesiochromite are much more common than pyroxene inclusions and are present within both olivine crystals and melt inclusions (Figure 2).
When our new whole rock tephra data and lava data from previous studies (Day et al., 2022) are compared, the overall trace-element patterns are remarkably similar.The results show depletions in fluid-mobile elements (e.g., K, Pb) and enrichment in high field strength elements (e.g., Nb, Ta, Figure 4a) and light rare earth elements (La/Yb > 10, Figure 4b), which are characteristic of oceanic island basalts (Hofmann, 2007).That said, the range of trace element variations observed within the melt inclusions is larger (e.g., La/Yb 20-40, Nb/Yb 20-40) than what is recorded in the lavas and tephras (La/Yb 30-35, Nb/Yb 35-45) (Figures 4b and 4c) (Maclennan et al., 2003).
Our vanadium-in-olivine estimates of melt inclusion oxygen fugacity of the 2021 eruption span +0.75 to +1.5 FMQ, trending slightly lower than those estimates found by Day et al. (2022) Canil, 2002).This is likely because the fO2 values of Day et al., 2022 are derived from bulk rock measurements of V, while our V values of in situ melt inclusion measurements from LA-ICP-MS are higher (318 ppm for lava whole rock data vs. 353 ppm V for melt inclusions).Nevertheless, estimates of the V values of olivine from both this work and Day et al. (2022) are very similar, being 5.7 versus 6 ppm, but our fO2 values for the 2021 eruption agree when the uncertainty of the V-in-olivine method is considered (∼±0.5 ΔFMQ) (Nicklas et al., 2022).Thus, compared to previous work on the Canary Islands, the 2021 La Palma eruption is slightly less oxidized, but still records an overall oxidized mantle source compared to mid-ocean ridge basalts (Moussallam et al., 2019;Nicklas et al., 2022;Taracsák et al., 2022).
Volatile concentrations of melt inclusions also display a systematic variation between olivine populations with the more evolved melt inclusions population (Mg# 78-82) containing higher concentrations of volatile elements (fluorine, chlorine, sulfur, H 2 O) (Figure 5).H 2 O concentrations noticeably separate the samples, with inclusions from the start of the eruption (LM6) recording ∼1.5-2.5 wt% H 2 O and those from the end of the eruption (LM0) showing ∼1.0-1.5 wt% H 2 O.Despite inclusions from LM0 being hosted in both low Mg# (78-82) and high Mg# (high 83-86) olivine, H 2 O contents of melt inclusions from the end of the eruption contain similar H 2 O contents regardless of Mg# of the host olivine core (Figures 5a-5d).Sulfur concentrations also divide each sample with early inclusions having between ∼3,000 and 4,000 ppm and late inclusions retaining between ∼2,000 and 3,500 ppm of sulfur (Figure 5b).However, the volatile divide between low and high Mg# olivine hosted melt inclusions is less noticeable for sulfur and chlorine than it is for fluorine (Figure 5a).Aside from the variations in H 2 O, two populations are apparent when fluorine contents are investigated regardless of the sample.Lower Mg# hosted melt inclusions typically contain ∼1,400-1,900 ppm  (Pearce & Norry, 1979) diagram identifying all eruptive materials as retaining Oceanic Island Basalt signatures.fluorine while higher Mg# crystals, all from LM0, contain ∼1,000-1,400 ppm fluorine.Interestingly, there are two populations of chlorine contents within the low Mg# hosted melt inclusions (550-750 ppm and <550 ppm), but there is only one population of chlorine contents within the high Mg# population of LM0 (Figure 5c).CO 2 contents of melt inclusion glasses are slightly higher for sample LM6, resulting in an apparent positive correlation with H 2 O (Figure 5d).That said, once the CO 2 within the melt inclusion vapor bubbles are accounted by mass balance, the CO 2 contents in melt inclusions from both samples range from 0.5 to 1.4 wt% CO 2 regardless of Mg# of the host population when inclusions that are potentially heterogeneously co-trapped are excluded (see Section 4.2 and Figures 6 and 7).

2021 La Palma Eruption Volatile Record and Canary Island Volatile Variability
As previously shown, H 2 O concentrations in melt inclusions from the same locality, even within single samples, can display considerable variations in H 2 O content (Hauri, 2002;Lloyd et al., 2012;Portnyagin et al., 2008).This has been found to depend on the grain size selected for inclusion analysis due to rapid H 2 O diffusivity upon clast cooling, and several studies have employed experimental techniques to reverse and constrain these changes (Bucholz et al., 2013;Gaetani et al., 2012;Portnyagin et al., 2008;Rasmussen et al., 2020;Wallace et al., 2015).To avoid these complications, we selected olivine crystals coated only with a thin veneer of glass sieved from a  6).Isobars are calculated using Magmasat in VESIcal for a representative composition (PEC-corrected LM6 G3 MI4) with the Moho (black bold line) at 14 km (Ghiorso & Gualda, 2015;Iacovino et al., 2021;Ranero et al., 1995;Tenzer et al., 2013).fine grain size fraction of 0.5-2 mm to ensure rapid quenching, following the best practice recommendations of Lloyd et al. (2012).

Geochemistry, Geophysics, Geosystems
We find that melt inclusions in lower Mg# (78-82) crystals erupted during the opening (sample LM6) of the eruption record the highest H 2 O content (up to 2.5 wt%).To properly evaluate the diffusive loss or gain of H 2 O during ascent, we plotted our H 2 O contents against other volatiles (e.g., fluorine, chlorine, sulfur) and measured the distance from the melt inclusion to the crystal edge along the a-axis, in which H 2 O diffuses the fastest (Barth et al., 2019;Ferriss et al., 2018;Le Voyer et al., 2014).We did not find major differences in melt inclusion volume or distance to crystal edge in either sample.As fluorine did not degas from the co-erupted matrix glasses (Figure 5a), we assume that fluorine did not diffuse from melt inclusions as concentrations remained in equilibrium with their carrier melt (Bucholz et al., 2013;Le Voyer et al., 2014;Portnyagin et al., 2008), but note the lack of consensus on the behavior of fluorine in olivine (Koleszar et al., 2009).When H 2 O is plotted against fluorine, we find minimal H 2 O loss in inclusions at the start of the eruption (LM6), with a few exceptions.We notice that inclusions very close to the edge of crystals along the a-axis (<30 μm) likely gained some H 2 O (Figure 5a) based on multiple inclusion measurements in single crystals as those closest to the grain edge have slightly elevated H 2 O contents (∼0.2 wt%) (e.g., LM6 Grain 2).We interpret this slight H 2 O gain as melt inclusion rehydration and incomplete re-equilibration during magma ascent as olivine crystals from LM6 encountered the evolved, hydrous, amphibole-bearing tephrite melt (see LM6 matrix glass composition, Figure 3) that was emitted during the first days of the eruption (Day et al., 2022;Pankhurst et al., 2022;Ubide et al., 2023).These low Mg# crystals were the primary crystal cargo during the start of the eruption (sample LM6) but were also sampled together with a lower H 2 O and higher Mg# (83-86) population toward the end of the eruption (both within sample LM0) (Figure 5).As mentioned, the melt inclusions from the end of the eruption are hosted in two populations of olivine, but regardless of Mg#, all the late erupted crystals contain melt inclusions which record concentrations of ∼1.0-1.5 wt% H 2 O.Given the lower fluorine, sulfur, and chlorine concentrations of the melt inclusions within the high Mg# olivine population (red circles, Figure 5), we interpret these inclusions as initially containing lower water concentrations, rather than being fully controlled by water loss.We identify the lower Mg# population that erupted at the end of the eruption (blue circles, Figure 5) as genetically related to the initial lower Mg# crystal cargo sampled at the beginning of the eruption in sample LM6 (blue squares).The lower H 2 O in the lower Mg# population from the end of the eruption is likely controlled by diffusive hydrogen loss, which occurs on the scale of hours in olivine within the Mg# range (78-82) (Barth et al., 2019).We explain this relationship among the lower Mg# olivine population seen at the end of the eruption by having resided in the volcanic plumbing system for a longer time and thus interacting with new, more H 2 O poor magma replenishments within the plumbing system.The influx of the drier, more primitive magma (recorded by melt inclusions with high CaO/Al 2 O 3 in the higher Mg# olivine crystals) (Figure 3) drove diffusive re-equilibration of the H 2 O contents of the lower Mg# hosted inclusions that erupted at the end of the eruption.
Overall, we see positive correlations between H 2 O and fluorine, chlorine, and sulfur (Figures 5a-5c) in these melt inclusions, with lower Mg# crystals containing inclusions higher in fluorine, chlorine, and sulfur regardless of sample.Generally, incompatible volatile elements are higher in the lower Mg# olivine populations, suggesting volatile increase connected with fractional crystallization.Our analyzed melt inclusions do not record evidence of sulfur saturation, identified by the presence of sulfide blebs (Hartley et al., 2017), except in two heterogeneously entrapped inclusions with anomalously low H 2 O (∼0.5 wt%) and high chlorine contents (∼1,000 ppm) (Figure 5; LM0 Grain 1 Inclusion 1 and LM0 Grain 23 Inclusion 1; see Supporting Information S1 and Data Set S1).Observed sulfide and ubiquitous magnesiochromite inclusions are identified through Raman spectroscopy (Lafuente et al., 2015).
Our fluorine contents of melt inclusions from the 2021 eruption align with those from the 1949 eruption (Duraznero vent), Barranco Fagundo deposits (of the main Taburiente shield volcano), and other Holocene cones on La Palma (Kirstein et al., 2023;Walowski et al., 2019), but we have found higher H 2 O (up to 2.5 wt%) agreeing with phase equilibrium experiments (Fabbrizio et al., 2023).Compared to the most volatile rich inclusions of the El Hierro 2011-2012 eruption, which contain over 3 wt% H 2 O, 5,000 ppm sulfur, and 2,000 ppm fluorine (Longpré et al., 2017) and other El Hierro inclusions (Taracsák et al., 2019), the inclusions of the 2021 La Palma eruption are slightly less volatile rich.Tenerife melt inclusions from 1705 Volcán de Fasnia (TF-16-08) and Montaña Bilma (TF-16-09A) eruptions (DeVitre, Gazel, et al., 2023) also contain very similar H 2 O, CO 2 (in melt inclusion glass), fluorine, chlorine, and sulfur contents as the 2021 La Palma eruption in similar olivine host populations, indicating that similar volatile systematics may apply on an inter-island scale.
We found half the amount of chlorine in our samples for the 2021 eruption (∼300-800 ppm, with later trapped inclusion containing ∼1,100 ppm) versus 600-1,900 ppm found in previous studies of melt inclusions from earlier eruptions on La Palma (Kirstein et al., 2023;Walowski et al., 2019).An overestimate in chlorine concentrations by the SIMS at Edinburgh Ion Microprobe Facility (EIMF) was noted by Taracsák et al. (2019) when analyzing melt inclusions from El Hierro, which may explain the chlorine discrepancy between melt inclusion populations from older La Palma eruptions as SIMS data was collected at EIMF (Kirstein et al., 2023;Walowski et al., 2019).To evaluate this discrepancy, the chlorine contents of our melt inclusions were also measured using EPMA, with results in agreement with SIMS data from WHOI (Supporting Information S1 and Data Set S1).Inclusions from the 2011-2012 eruption of El Hierro (Longpre et al., 2017), also measured at WHOI, record chlorine contents (∼600-1,300 ppm) higher than those from La Palma but lower than those measured at EIMF (150-1,700 ppm; Taracsák et al., 2019).Chlorine contents are especially important when comparing F/Cl ratios, with recent work (Hartley et al., 2022) finding average F/Cl ratios of melt inclusions from the 2021 eruption of ∼1.9, while we find average F/Cl of ∼3 in our melt inclusions.
Despite lower fluorine, chlorine, and sulfur when compared to previous La Palma and El Hierro eruptions, the glass phase of the 2021 melt inclusions contains some of the highest dissolved CO 2 of the Canary Islands.Maximum CO 2 contents reach up to 6,000 ppm versus 4,000 ppm in the 1949 eruption products and Barranco Fagundo deposits (Walowski et al., 2019) and ∼3,500 ppm in El Hierro eruption products (Longpré et al., 2017;Taracsák et al., 2019).Maximum CO 2 contents of melt inclusion glasses may also be due to compositional variations among eruptions, resulting in differences in CO 2 solubility as well as differences in CO 2 partition within vapor bubbles (Moore et al., 2015;Shishkina et al., 2014).

2021 La Palma Eruption Volatile Record and Canary Island Volatile Variability
While a lack of vapor bubble data from previous La Palma eruptions hinders comparisons between eruptions, variations in total CO 2 contents in our combined vapor bubble and glass data can be compared to identify changes within the 2021 eruption.Almost all selected melt inclusions from both sample sets had CO 2 detectable by Raman spectroscopy in vapor bubbles with CO 2 densities ranging from 0.19 to 0.51 g/cm 3 .On the other hand, nondecrepitated fluid inclusions retained CO 2 densities between ∼0.74-0.98 g/cm 3 (Figure 6a) (Dayton et al., 2023).Specifically, LM0 Grain 2 Inclusion 2 (pictured in Figures 2a and 2b) had no detectable CO 2 in the vapor bubble.Other inclusions with a similar morphology (rounded, smooth, and unfaceted) have also been found to not contain CO 2 in the vapor bubble, as these vapor bubbles formed during quenching and CO 2 did not have time to diffuse into the vapor bubble (cf.Maclennan, 2017;Wieser et al., 2021).This shallower, later-trapped inclusion (LM0 Grain 2 Inclusion 2, 15.8 ± 0.4 km) was large enough to allow for SIMS data collection and comparison with a more primitive, deeper-trapped inclusion (LM0 Grain 2 Inclusion 1, 19 ± 1.5 km) trapped within the same crystal (Figures 2a and 2b).
As evident by inclusions with large bubble-volume percent and the presence of fluid inclusions in phenocrysts (Dayton et al., 2023), heterogenous entrapment of exsolved CO 2 fluid occurs in some melt inclusions (e.g., Steele-MacInnis et al., 2017).From our selected melt inclusions, vapor bubbles occupied between ∼3%-20% of the melt inclusion volume (minus co-trapped oxide volumes; Figures 6b and 6c).We find that inclusions with a bubblevolume percent above 7.5 ± 1.5 for LM6 and 10 ± 1.5 for LM0 (plotted with opaque symbols on Figures 6b and  6c) and total reconstructed (vapor bubble plus glass) PEC-corrected CO 2 above ∼1.4wt% indicate heterogenous entrapment.We note that vapor-bubble percent alone cannot identify heterogenous entrapment and inclusions must be measured in all three dimensions to obtain the best of vapor-bubble-volume percent (DeVitre, Gazel, et al., 2023;Mironov et al., 2020).If the third dimensions of these melt inclusions were assumed from measurements of the exposed two dimensions, vapor bubble percents of these inclusions would be significantly lower with an average of 2.3% for this sample set.Whereas when the shortest, dimension of these melt inclusions are measured, in this case with Raman spectroscopy (see Data Set S1), this sample set has an average vapor bubble percent of 6.5%.
In addition to anomalously high reconstructed PEC-corrected CO 2 contents (bubble plus glass), saturation pressures, and bubble volume percent, we find that inclusions with more than ∼75% CO 2 in the vapor bubble for LM0 and 70% in the vapor bubble for LM6 also heterogeneously entrapped melt inclusions in these samples.We find that vapor bubbles of inclusions within the dehydrated, low Mg# hosted melt inclusions of LM0 erupted at the end of the eruption (Figure 5) contain the highest percentage of CO 2 .We note that heterogeneous entrapment is more common in crystals erupted during the end of the eruption within LM0 (Figure 6b), potentially indicating more exsolved CO 2 fluid as the eruption progressed.
Heterogeneous entrapment of melt inclusions can further identified by analyzing multiple inclusions within the same inclusion assemblage in the same crystal.LM6 Grain 10 contains six melt inclusions analyzed with all analytical techniques (Raman, SIMS, EPMA, LA-ICP-MS) (Figure 7a).Inclusion 4 is heterogeneously entrapped and records a spurious saturation pressure of 900 MPa, while the five other melt record pressures of ∼700 MPa when total CO 2 contents are reconstructed (bubble plus glass).Upon further investigation, LM6 Grain 10 Inclusion 4 has a bubble volume percent of 7% and a CO 2 density of 0.47 g/cm 3 , while other five inclusions have bubble volume percent between 3.4% and 5.7% and CO 2 densities between 0.33 and 0.35 g/cm 3 .
Unfortunately, identifying heterogenous melt inclusions from previous studies is difficult without accurate measurements of melt inclusions in all three dimensions.When inclusions from the La Palma 2021 eruption are compared to those from other studies of Island melt inclusions that reported 2 in melt inclusion vapor bubbles, the La 2021 eruption melt inclusions contains more CO 2 (up to 1.4 wt%) than those of El Hierro (up to ∼1.3 wt%, Taracsák et al., 2019) and Tenerife (up to 1.3 wt%, DeVitre, Gazel, et al., 2023), though these highest estimates of other studies may be inaccurate as heterogenous entrapped inclusions cannot be identified.Nevertheless, this indicates that a deeper magma storage reservoir fed the 2021 La Palma eruption as recorded through higher CO 2 contents.

Comparing Melt and Fluid Inclusions in Single Crystals
To investigate if melt and fluid inclusions record equivalent pressures within the 2021 magma reservoir, we identified grains which contained both melt and fluid inclusions within a single crystal to directly compare preserved pressures from both inclusion types in grains with similar histories (Figure 7b).Overall, grains containing both melt and fluid inclusions are rare.Even if both types of inclusions are within the same crystal, inclusions must be suitable for comparisons.Fluid inclusions must not be decrepitated, so they retain their original trapping pressures (unlike the decrepitated fluid inclusion in Figure 7a) (Campione, 2018;Campione et al., 2015;Wanamaker et al., 1990) and melt inclusions must not be entrapped (Figures 2c and 7a) to preserve their true trapping pressure (Steele-MacInnis et al., 2017).To complicate this comparison further, we find that melt inclusions within grains that also have fluid inclusions are commonly heterogeneously entrapped as crystal growth may occur near areas with exsolved CO 2 fluid (ex.LM0 Grain 25 Inclusion 2) (Figure 2c).
Despite fluid and melt inclusions being trapped within the same grain, careful consideration must also be taken to ensure the inclusions are petrogenetically related and they belong to the same fluid-inclusion assemblage (FIA) (Goldstein & Reynolds, 1994) or in this case, inclusion assemblage.This is displayed in LM6 Grain 28 (Figure 2d, note these images are focus-stacked to allow photographing of inclusions at different depth within the same crystal in a single image), which contains multiple inclusions assemblages.Melt inclusions 3 and 4 belong to the same assemblage recording pressures of 735-756 ± 60 MPa (25.5-26 ± 2 km).While a ∼3 µm fluid inclusion trapped within a different inclusion assemblage (inclusion 5; Figure 2d), cross-cutting below the inclusion assemblage that contains inclusions 3 and 4 within this same crystal records a pressure of 590 ± 22 MPa (17.6 ± 0.7 km).
On the rare occasion that both suitable melt and fluid inclusions are trapped within the same inclusion assemblage (Goldstein & Reynolds, 1994), a comparison between inclusion types can be evaluated.In LM0 Grain 30 (Figure 7b), both fluid inclusions (inclusions 2, 5, 6, and 7) and a melt inclusion (inclusion 8) are trapped within a single inclusion assemblage and record pressure between 560 and 600 MPa (∼20 km).This indicates that inclusions trapped within the same assemblage may record similar pressures and single crystals can record multiple events at different depths by analyzing multiple assemblages and different inclusion types.By analyzing both melt and fluid inclusions within crystals, we are also allowed a glimpse into a single crystal's journey through a magma chamber, identifying magmatic processes at different depths.

Reconstructing Magmatic Plumbing System
Overall, melt inclusions retain slightly higher pressures than fluid inclusions, but average depths and pressures are within error of each other (839 ± 75 MPa or 29 ± 2.5 km for the deepest melt inclusion vs. 780 ± 27 MPa or 27 ± 0.9 km for the deepest fluid inclusion; Dayton et al., 2023).This may be due to the higher uncertainty of saturation pressure calculations using a solubility model for melt inclusions or as the fluid inclusion densities of Dayton et al. (2023) assume a pure CO 2 fluid using the equation of state of Span and Wagner (1996) (Figure 8).As the quantity of H 2 O in CO 2 fluid inclusions could not be directly determined, we re-calculate fluid-inclusion pressures and depths using the molar fraction of H 2 O (XH2O) values of melt inclusions with MagmaSat (Ghiorso & Gualda, 2015;Iacovino et al., 2021).We use average values of XH2O of 0.15 for LM6 and 0.05 for LM0 and utilize the mixed CO 2 -H 2 O equation of state of Duan and Zhang (2006) republished by Yoshimura ( 2023) and implemented into DiadFit (Wieser & DeVitre, 2023).If XH2O of melt inclusions is applied to fluid inclusions within the same samples, we find that fluid inclusions for LM6 record entrapment depths 6 km (186 MPa) deeper and those from LM0 would be placed 2 km (54 Mpa) deeper (see Supporting Information S1 and Data Set S1).Despite these changes in entrapment pressure, melt and fluid inclusions still align within error.We note that as more work is needed on the diffusion of H 2 O out of CO 2 fluid inclusions at magmatic temperatures, we discuss the fluid-inclusion data as published in Dayton et al. (2023).There is a lack of fluid inclusions retaining pressures from the 2021 eruption above 780 ± 27 MPa (27 ± 0.89 km), which may indicate that CO 2 did not significantly exsolve deeper than ∼27 km, corresponding with the deepest depths of earthquake epicenters recorded during the 2021 eruption (D'Auria et al., 2022;Dayton et al., 2023).
In addition to intact fluid inclusions, decrepitated fluid inclusions were analyzed.Inclusions with any optically identifiable indications of decrepitation (Campione, 2018;Campione et al., 2015;Wanamaker et al., 1990) were excluded from Dayton et al. (2023) but are now evaluated further (Figure 8).Some fluid inclusions with optically identifiable features of decrepitation still retain densities of ∼0.7 g/cm 3 , indicating decrepitation and equilibration within the deep mantle chamber prior to further ascent.Generally, decrepitated fluid inclusions retain two groups of separate pressure ranges in addition to the deep mantle magma reservoir.
One group of decrepitated fluid inclusions retained CO 2 densities between 0.2 and 0.6 g/cm 3 , corresponding to pressures between 60 and 340 MPa (2.5-12.5 km) and aligned with the upper earthquake swarm within the crust (Figure 8) (D'Auria et al., 2022;Del Fresno et al., 2023).The second group of decrepitated fluid inclusions displayed ultra-low densities between 0.001 and 0.05 g/cm 3 , with the lowest densities right at the detection limit of our Raman spectroscopy method (DeVitre et al., 2021).These ultra-low-density inclusions indicate fluid inclusion decrepitation upon ascent and preserve pressures of under 10 MPa (or less than 0.5 km).The preservation of CO 2 within decrepitated fluid inclusions indicates rapid olivine fracture healing following loss of initial trapping pressure (Wanamaker et al., 1990).This allows decrepitated inclusions to record additional, shallower, magma storage depths as they rapidly re-equilibrate to their evolving storage conditions.

Geochemical Evolution of the 2021 La Palma Eruption
Despite relatively uniform compositions, subtle geochemical changes in eruptive products can serve as indicators of magmatic processes at depth.While in a TAS diagram (Figure 3a) (Le Bas et al., 1986) erupted products (excluding matrix glasses) generally plot together as tephrites/basanites, changes in the magma supply are apparent in lavas and melt inclusions when investigating CaO/Al 2 O 3 versus SiO 2 for the duration of the eruption (Figure 3b).Though whole rock composition of the 2021 lavas identifies a generally homogenous mantle source with input of a young HIMU recycled oceanic lithosphere and relatively uniform degree of partial melting (2.5%-3%) (Day et al., 2022), variations in Os isotopes identify fractional crystallization or even a degree of crustal contamination, specifically in eruptive products from the first 20 days of eruption.This differentiation of early erupted magma results in the evolved compositions of the matrix glasses from the first days of the eruption (LM6) and the early eruption of amphibole.In addition, this is in line with the slight re-hydration of melt inclusions within olivine erupted during the first days of the eruption (cf.Day et al., 2022;Ubide et al., 2023) (Figures 5a and  et al., 2022) indicate two main reservoirs.Fluid inclusions (Dayton et al., 2023;this work), and melt inclusions plot together below the Moho at 14 km (Ranero et al., 1995;Tenzer et al., 2013).Heterogeneously entrapped melt inclusions are plotted in gray.Decrepitated fluid inclusions in the olivine plot in two groups, one between 2 and 14 km aligned with the shallower earthquakes and upper magma chamber within the crust, and another population recorded depths of ∼0.3 km, which we interpret as representing decrepitation upon final magma ascent.9).Notably, melt inclusions within the low Mg# olivine population also record this pre-eruptive fractionation through lower CaO/Al 2 O 3 ratios (Figure 3b) and inclusions of clinopyroxene (Figures 2a and 2e).
In situ analysis uncovers a more detailed picture of distinct melt batches exhumed during the eruption.In situ Sr isotope data traces two distinct magma batches, and identifies a colder, more radiogenic, evolved crystal mush at the beginning of the eruption (high 87 Sr/ 86 Sr in andesine) and a less radiogenic, mafic melt (lower 87 Sr/ 86 Sr in labradorite microcrysts) (Ubide et al., 2023).Melt inclusion compositions along with isotope studies (Day et al., 2022;Sandoval-Velasquez et al., 2023;Ubide et al., 2023) identify the exhumation of this new batch of less fractionated magma in November 2021 (Day et al., 2022;Ubide et al., 2023) recorded by the change in CaO/ Al 2 O 3 versus SiO 2 of lava around day 50 of the eruption and is sampled by melt inclusions within the high Mg# olivine population in LM0 (Figure 3b).
Our melt inclusion data further constrains this process to specific depths, and although trace element variation between samples is small, volatile contents within melt inclusions trace fractional crystallization and pulses of hotter primitive magma replenishment (Figures 5 and 9).Melt inclusions that record this fractionation processes with elevated H 2 O, sulfur, chlorine, and fluorine (hosted within low-Mg# olivine) record pressures spanning the entire deep magma reservoir (17.6 ± 1.3-29 ± 2.5 km), indicating that fractionation occurs throughout the entire deep chamber system.As the eruption began to be overtaken by the more primitive magma supply (Figure 9), we find that these lower-Mg# crystals containing inclusions of the early appearing derivative magmas become less frequent as this presumably older crystal mush is evacuated from the deep chamber and becomes entrained with the higher-Mg# crystals recording the emplacement of more primitive melt.Our analyzed melt inclusions recording this primitive melt record pressures between 14.9 ± 1.4 and 28.1 ± 2.7 km and align with fluidinclusion pressures from the same sample (15-27 km) (Dayton et al., 2023), tracing this primitive magma through the entire deep magma chamber.Melt inclusions hosted within more evolved crystal cargoes that crystallized from a more H 2 O-rich evolved magma likely underwent diffusive re-equilibration with the more primitive, H 2 O-poor replenishing magma, resulting in the significantly lower H 2 O contents of the melt inclusions in the lower Mg# olivine of LM0 (blue circles, Figure 5) compared with the same population at the beginning of the eruption within LM6 (blue squares, Figure 5).This natural dehydration has implications for melt-inclusion studies which characterize the volatile contents of an eruption using a single sample, as our results show that melt inclusion H 2 O contents can vary significantly between different samples (Bucholz et al., 2013;Gaetani et al., 2012;Portnyagin et al., 2008).Prior to eruption, magma fractionates within both the shallow and deep reservoirs, creating the evolved tephrites and basanites.Following magma replenishment, this evolved composition, along with more evolved low Mg# olivine phenocrysts retaining the more volatile enriched melt inclusions, was cleared out of the magma chamber.Following replenishment, the primitive lower volatile-rich composition begins to dominate the eruption, entraining the more primitive higher Mg# olivine phenocrysts and melt inclusions.

Volatile Budgets of the 2121 La Palma Eruption
Our Monte Carlo simulations suggest that over the span of the 2021 eruption, 2.4 ± 0.4 Mt of SO 2 would have been emitted from the degassing (determined from the difference between melt inclusion and matrix glass composition) of 0.25 km 3 of eruptive materials.This is higher, but within error of the 1.84 ± 0.92 Mt estimates from satellite-constrained emissions (Milford et al., 2023;Valade et al., 2019) (Figure 10a).This discrepancy can be reconciled in a number of ways.First, it is possible that satellites underestimate the emission rate if the cloud fraction is high and plume altitude is low (Esse et al., 2023).Second, it may be that some proportion of the erupted lava did not experience as much degassing as the matrix glasses measured here (e.g., lava transported in closed tubes).We address this scenario by performing the same calculation assuming that the average amount of degassing of 0.25 km 3 of erupted material was 85% and 70%.If only 70% of the accounted for erupted volume is Figure 10.Atmospheric emissions of (a) SO 2 calculated for the Tajogaite eruption using the petrographic method of Sharma et al. (2004).Errors in input parameters were propagated using Monte Carlo simulations and budgets are reported for three scenarios: full degassing (melt inclusion-matrix glass), 85% degassing, and 70% degassing of erupted material using the 0.25 km 3 volume estimates of Civico et al. (2022).Reported satellite data and error on SO 2 emissions are from Milford et al. (2023) and Valade et al. (2019).(b) CO 2 budget calculated using the same methods as for SO 2, but for scenarios where for 0.25 km 3 (actual volume; Civico et al., 2022), 0.75 km 3 , and 1.25 km 3 of magma degassed.Black lines indicate the CO 2 estimates and error of Burton et al. (2023).
Geochemistry, Geophysics, Geosystems 10.1029/2024GC011491 degassed, SO 2 emissions from the petrographic method (2.0 ± 0.4 Mt) and satellite estimates of 1.84 ± 0.92 Mt (Milford et al., 2023) more closely align (Figure 10a).Third, there are also other sources of uncertainty that may resolve offsets between satellite and petrological methods.The presence of fumarolic deposits containing native sulfur and sulfates (as well as carbonates, fluorides, and chlorides, Campeny et al., 2023;Martínez-Martínez et al., 2023), and the formation of sulfate crystals on fine volcanic ash (Hornby et al., 2023) during the 2021 indicate that some amount of emitted sulfur was sequestered prior to detection by satellites.Using the 100% degassing of 0.25 km 3 of the eruptive volume scenario, it is noteworthy that the La Palma 2021 eruption's emission of 2.4 ± 0.4 Mt of SO 2 pales in comparison to the 2018 lower East Rift Zone eruption of Kīlauea which emitted ∼10 Mt of SO 2 (Kern et al., 2020).Despite the higher sulfur concentrations of the La Palma melt inclusions (3,062 ± 500 ppm) and matrix glasses (345 ± 53 ppm) compared to Kīlauea 2018 melt inclusions (400-900 ppm) and matrix glasses (100-900 ppm) (Lerner et al., 2021), the Kīlauea 2018 eruptive volume (0.9-1.4 km 3 DRE) was significantly larger than the 2021 La Palma eruption (0.21-0.24 km 3 DRE) (Dietterich et al., 2021).For the La Palma 2021 eruption to release as much SO 2 as the Kīlauea 2018 eruption, only 0.84-0.96km 3 DRE would be needed to emit 9.7 ± 1.8 Mt (Figure S6 in Supporting Information S1).
Using the same approach as for SO 2 , we find 5.4 ± 1.0 Mt of CO 2 were released in La Palma during the 2021 eruption when assuming degassing of the 0.25 km 3 of erupted volume (Figure 10b).These are minimum values for CO 2 , as any volatile emissions from magma which ascended to shallower reservoirs but did not erupt would not be accounted for using this method.We note that our calculated values are five times lower than other previously published values, which find 28 ± 14 Mt of CO 2 emitted for the 2021 eruption based on plume CO 2 / SO 2 ratios and solubility modeling (Burton et al., 2023) (Figure 10b).Assuming the same CO 2 concentrations as found in melt inclusions and matrix glasses, we find that degassing ∼1.25 km 3 of magma is needed to release 27.0 ± 5.1 Mt of CO 2 to align with the estimates of Burton et al. (2023) (Figure 10b), indicating that 80% of magma needed to satisfy this budget did not erupt.
Discrepancies clearly exist among SO 2 and CO 2 emissions during the 2021 eruption, as measured SO 2 values via satellite (Milford et al., 2023) and those of CO 2 from solubility models (Burton et al., 2023) would require degassing of significantly different quantities of magma (<0.25 km 3 for SO 2 vs. ∼1.25 km 3 for CO 2 ), if they are to be directly compared.If five times the measured eruptive volume, or 1.25 km 3 , of La Palma 2021 magma degassed, 12.2 ± 2.2 Mt of SO 2 , or about five times what was detected via the satellite method (1.84 ± 0.92 Mt SO 2 , Milford et al., 2023), would have been emitted (Figure S6 in Supporting Information S1).This indicates that the satellite method (Milford et al., 2023) and petrologic method (Sharma et al., 2004) would severely underestimate SO 2 emissions, or CO 2 /SO 2 ratios and solubility modeling (Burton et al., 2023) may overestimate CO 2 emissions if the same amount of magma was unlikely, completely degassed of all volatiles.Each method has weaknesses, with the satellite method not as effective for low-plume altitude eruptions, and the petrographic method unable to account for previously degassed CO 2 .Despite these discrepancies, if 28 Mt of CO 2 was released by the 2021 La Palma eruption, it would have accounted for an unexceptional 0.0008% of the 34.9 Gt total global CO 2 emitted in 2021 (Liu et al., 2022), whereas 0.00015% of global CO 2 if the 5.4 ± 1.0 Mt of the petrographic method is correct.

Conclusions
• The melt inclusions and olivine host compositions in Tajogaite tephras provide chemical insight into a deep magma chamber filled with an evolving fractionating magma prior to eruption.This pre-eruptive fractionation enriched incompatible volatiles recorded by melt inclusions in low Mg# olivine and allowed cotectic olivine and pyroxene crystallization prior to the eruption.A pulse of primitive magma evidently perturbed this crystallizing deep magma reservoir, recorded by melt inclusions in higher Mg# olivine, leading to the eruption of a mixed crystal record.• Despite their differences, we find that both melt and fluid inclusions can reliably identify magma storage depths.Though melt inclusions can further elucidate distinct chemical changes in magma source, if rapid magma depth information is needed, fluid inclusions are a more rapid, precise, accurate, and more costeffective tool to locate previous magma storage levels if they are present within samples.• We find utility in decrepitated fluid inclusions, and stress that careful petrographic observations can allow their use to identify multiple (shallower) levels of magma storage.
• Complicating melt inclusion work, we find that H 2 O contents of melt inclusions within the same olivine population can vary by 1.5 wt% within a single eruption, indicating that previous studies which may have only selected one sample per-eruption may incompletely characterize the diversity of magmatic water contents.• We note that volatile budgets from petrological, satellite, and ground-based methods are inherently different and must be comprehensively evaluated to understand the true volatile budgets of volcanic eruptions.This work was funded by NSF Grants EAR 2318614 and 2216738 to EG and 2217371 to PW. VRT acknowledges support from the Swedish Research Council.FJPT and MA thank to the LAJIAL (PGC2018-101027-B-I00, MCIU/AEI/FEDER, EU) and MESVOL (SD RD 1078/2021 LA PALMA) research projects for financial support of their fieldwork.We thank Brian Monteleone and Mike Jercinovic for assistance with SIMS and EPMA analysis.Discussions with Terry Plank, Glenn Gaetani, Sandro Aiuppa, and Paul Wallace improved this manuscript.We thank Jake Lowenstern and Jacqueline Dixon for their thoughtful reviews which improved this manuscript.We are grateful to PEVOLCA for enabling access to the Exclusion Zone of La Palma during the 2021 eruption for sample collection.

Figure 1 .
Figure 1.Geographic overview and sample collection.(a) The Canary Islands, with La Palma circled.(b) Simple geologic stages of La Palma, with the younger Cumbre Vieja ridge extending south from the older and extinct Taburiente shield volcano.(c) Historical volcanism of La Palma modified from Carracedo et al. (2022).Sample locations marked with a white star.(d) Stratigraphic section collected from the middle of road LP-212 on 17th January 2022 in Las Manchas (28.606552°N, 17.878136°W) following the conclusion of the 2021 eruption.The blue square and red circle indicate samples with melt inclusion data included in this work.Green stars indicate samples with fluid inclusion data published by Dayton et al. (2023).Dates of samples are derived from the stratigraphic column of Romero et al. (2022).

Figure 2 .
Figure 2. Melt Inclusion Photomicrographs (note: these are focus stacked images to allow inclusions at different depths to be visible within a single image).(a) LM0 Grain 2 with a clinopyroxene (CPX) inclusion and two analyzed melt inclusions in the grain center.(b) Magnified view of inclusions one and two in panel (a).(c) LM0 Grain 25 containing both fluid and melt inclusions.This grain is polished down the a-axis to give a "side" view of melt inclusions.Inclusion 2 is heterogeneously entrapped.(d) LM6 Grain 28 contains multiple inclusion assemblages recording different pressures.(e) LM0 G35 grain contains both a melt inclusion and clinopyroxene (CPX) inclusion trapped prior to oxide crystallization.(f) LM6 Grain 10 Inclusion five.Typical morphology of melt inclusions selected for analysis.

Figure 3 .
Figure 3. Major element (wt%) composition of melt inclusions and eruptive materials of the 2021 eruption.(a) Total Alkalis (Na 2 O + K 2 O) versus SiO 2 diagram (LeBas et al., 1986) including melt inclusion and matrix glass from samples LM6 (September) and LM0 (December).Tephra (this work) and lava whole rock data(Day et al., 2022) span the entire Tajogaite eruption.(b) CaO/Al 2 O 3 versus SiO 2 diagram identifying the "geochemical kink" in November 2021 and transition from fractionated, more evolved magma representing the start of the eruption to the more primitive composition present through the end of the eruption(Day et al., 2022;Gonzalez-Garcia et al., 2023;Ubide et al., 2023).

Figure 4 .
Figure 4. Trace elements from melt inclusions and eruptive materials (tephra and lava) of the Tajogaite eruption.Lava data from Day et al. (2022).(a) Spider diagram normalized to pyrolite mantle (McDonough & Sun, 1995) displaying the compositional homogeneity of eruptive materials.(b) Dy/Yb versus.La/Yb diagram indicating low degree of melt and relative homogeneity among magmas in the garnet peridotite field.(c) Th/Yb versus Nb/Yb discrimination(Pearce & Norry, 1979) diagram identifying all eruptive materials as retaining Oceanic Island Basalt signatures.

Figure 5 .
Figure 5. Volatiles (Fluorine, Sulfur, Chlorine, and CO 2 vs. H 2 O) of Melt Inclusions and Matrix Glass.(a) Fluorine contents decrease with decreasing H 2 O over the course of the eruption.Inclusions in low Mg# olivine from the start of the eruption (LM6) closest to the crystal rims (<30 μm) display some H 2 O gain and inclusions trapped in low Mg# crystals erupted toward the end of the eruption (LM0) display the most H 2 O loss.Matrix glasses retain a majority of initial fluorine contents.(b) Sulfur decreases with decreasing H 2 O in both inclusions and matrix glasses.(c) There is also a decrease of chlorine with H 2 O, but the two inclusions with 0.5wt% H 2 O have chlorine contents around ∼1,100 ppm, indicating later trapping of degassing melt.(d) Melt inclusion total CO 2 (vapor bubble + glass) and melt inclusion glass-only CO 2 (transparent symbols) versus H 2 O. Gray symbols are melt inclusions (bubble + glass), which are interpreted as heterogeneously entrapped (Figure6).Isobars are calculated using Magmasat in VESIcal for a representative composition (PEC-corrected LM6 G3 MI4) with the Moho (black bold line) at 14 km(Ghiorso & Gualda, 2015;Iacovino et al., 2021;Ranero et al., 1995;Tenzer et al., 2013).

Figure 6 .
Figure 6.Evaluation of CO 2 density, melt inclusion volumes, and heterogeneous entrapment.(a) CO 2 density (g/cm 3 ) of melt inclusion vapor bubbles and non-decrepitated fluid inclusions from this work and Dayton et al. (2023) versus inclusion diameter (μm).(b) Melt inclusion volume versus bubble volume for sample LM6 (start of the eruption) sized and colored by total PEC corrected CO 2 .The inclusions circled in gray are interpreted to be heterogeneously entrapped.(c) Melt inclusion volume versus bubble volume for sample LM0 (end of the eruption) sized and colored by total PEC corrected CO 2 .The inclusions circled in gray are interpreted to be heterogeneously entrapped.

Figure 7 .
Figure 7. Multiple Inclusions in Single Crystals (note: these are focus stacked images to allow inclusions at different depths to be visible within a single image) (a) LM6 Grain 10 contains six analyzed melt inclusions with all inclusions preserving approximately the same H 2 O contents (∼2.2-2.3 wt%).Five of the six inclusions retain saturation pressures of ∼700 MPa, or ∼24.5 km and a CO 2 vapor bubble densities of 0.33-0.35g/cm 3 .Inclusion four is heterogeneously entrapped recording a saturation pressure of ∼900 MPa, ∼31.6 km and a CO 2 vapor bubble density 0.47 g/cm 3 .(b) LM0 Grain 30 with four fluid inclusions and one melt inclusion.All inclusions retain entrapment pressures within error.

Figure 8 .
Figure 8. Cross-section of the magmatic plumbing system under Cumbre Vieja.Deep and shallow earthquakes (D'Auriaet al., 2022)  indicate two main reservoirs.Fluid inclusions(Dayton et al., 2023; this work), and melt inclusions plot together below the Moho at 14 km(Ranero et al., 1995;Tenzer et al., 2013).Heterogeneously entrapped melt inclusions are plotted in gray.Decrepitated fluid inclusions in the olivine plot in two groups, one between 2 and 14 km aligned with the shallower earthquakes and upper magma chamber within the crust, and another population recorded depths of ∼0.3 km, which we interpret as representing decrepitation upon final magma ascent.

Figure 9 .
Figure 9. Schematic diagram of the Tajogaite eruption based on the results of this study.Prior to eruption, magma fractionates within both the shallow and deep reservoirs, creating the evolved tephrites and basanites.Following magma replenishment, this evolved composition, along with more evolved low Mg# olivine phenocrysts retaining the more volatile enriched melt inclusions, was cleared out of the magma chamber.Following replenishment, the primitive lower volatile-rich composition begins to dominate the eruption, entraining the more primitive higher Mg# olivine phenocrysts and melt inclusions.