A Myriad of Melt Inclusions: A 3D Analysis of Melt Inclusions Reveals the Gas-Rich Magma Reservoir of Colli Albani Volcano (Italy)

Abstract

Melt inclusions are small parcels of magma trapped in crystals, which hold key information about pre-eruptive magmatic conditions, including volatile content and melt chemistry. We focus here on melt inclusions for a nuanced view of the magmatic pre-eruptive state of Colli Albani, a mafic-alkaline ignimbrite forming system in central Italy. Recent years have seen an increased amount of research surrounding the feasibility of using melt inclusions as tracers of pre-eruptive magma volatile content including, namely the concern about measuring trapped CO₂ present in vapour bubbles. Here we present synchrotron 3D tomographic scans of over 2000 melt inclusions from 35 pyroxene and leucite crystals from Colli Albani. We show that incorporating 3D information into melt inclusion geometry calculations allows for the development of a novel classification scheme, which we then use to draw inferences about the pre-eruptive evolution of the plumbing system of Colli Albani. We define six types of melt inclusions based on shape, crystallinity, and the characteristics of their vapour bubble. We also identify a strong variability in melt inclusion type proportions with stratigraphy, which ultimately is a reflection of variation in reservoir residence times, magma ascent rates, and tephra quench rates. Additionally, a large number of melt inclusions have large-volume bubbles, suggesting the magma reservoir was bubble bearing at the time of melt inclusion trapping. We suggest that this is essential to prepare the eruption of large volumes of the low-viscosity magma at Colli Albani.

INTRODUCTION

How magma accumulates, what pressure and temperatures it resides at, how and if it chemically evolves, and the role of volatiles, are all fundamental questions that serve to link pre-eruptive processes to eruptions of different magnitude and dynamics. These processes we cannot directly observe, so we focus on erupted products to draw conclusions about the state of the magma prior to eruption. One way to view these processes is by investigating melt inclusions (MI), which are small droplets of magma that get trapped in crystals sampling magma before its eruption (Kent, 2008; Lowenstern, 1995; Sobolev, 1996; De Vivo & Bodnar, 2003). While MI are extremely useful time capsules, they can be subject to a number of processes which can change their chemistry (Moore et al., 2015; Audétat & Lowenstern, 2013; Wallace, 2005). Additionally, due to the difficulty of sample preparation and advanced methods required for analysis, only a small number of MI are studied which can result in a biased view of the pre-eruptive magma storage conditions.

In an ideal world, all MI would quench to a single homogeneous phase (also called glassy) where one can simply measure the glass to gather all the information on the volatile content and inclusion chemistry. However, this is very rarely the case and instead MI often include vapour bubbles or are chemically altered during or after trapping via crystallization of the host phase in the MI walls in a process called post entrapment crystallization (PEC; Kent, 2008; Roedder, 1979; Danyushevsky et al., 2002; Moore et al., 2015; Audétat & Lowenstern, 2013; Wallace, 2005). PEC changes MI composition by depleting it in elements compatible with the host phase and inducing vapour bubble (VB) exsolution with changes to volatile solubility of the melt and by changing the MI pressure (Maclennan, 2017; Kent, 2008; Steele-MacInnis et al., 2011; Aster et al., 2016). Additionally, slow cooling can induce crystallization of daughter phases in the MI creating microcrystalline inclusions (Roedder, 1979; Danyushevsky et al., 2002). In the case of post-entrapment bubble formation, the VB sequesters volatiles, such as CO₂ and H₂O, which causes underestimated values of volatile species dissolved in the melt (Esposito et al., 2016). During cooling, differential thermal contraction of host mineral and melt will cause the pressure inside the MI to drop and thus volatiles that were soluble at high pressures will become insoluble and exsolve (Aster et al., 2016; Maclennan, 2017; Ferrero & Angel, 2018; Esposito et al., 2011). CO₂ is less soluble in melts than H₂O and its concentration in the melt decreases when vapour bubbles are formed (Ghiorso & Gualda, 2015). Several studies show that up to 90 wt % of the total CO₂ budget of MI can reside in the vapour bubble (Buso et al., 2022; Aster et al., 2016; Tucker et al., 2019). An accurate measurement of dissolved CO₂ requires measurements of the volatile content of the VB, which can be achieved by combining Raman spectroscopy and techniques to estimate the bubble volume within the MI (Aster et al., 2016; Moore et al., 2015; Wallace et al., 2015; Hartley et al., 2014; Hanyu et al., 2020). Typically, photomicrographs combined with the assumption of a spherical vapour bubble are used to estimate volume (Hanyu et al., 2020; Aster et al., 2016). However, if the MI or the host is opaque, it is difficult to measure VB dimensions. Furthermore, this study points to errors associated with volume reconstructions from photomicrographs as this is the largest error when determining CO₂ content from Raman spectroscopy, which can result in a 20% to 40% relative error of the CO₂ content of the MI (Hanyu et al., 2020; van Gerve et al., 2024). Developments in 3D X-ray computed tomography have shown to provide high-resolution volume measurements of melt inclusions and vapour bubbles, which have been used better constrain CO₂ budgets (Pamukcu et al., 2013; Pamukcu et al., 2015; Richard et al., 2019; van Gerve et al., 2024).

While there are a myriad of complications that exist for chemical analysis, textural evidence and volume fractions of VB can be utilized to reveal key information about the magma saturation state. Studies suggest that post entrapment shrinkage bubbles are typically 0.2 to 5 vol % of the inclusion (Hartley et al., 2014; Lowenstern, 1995; Lowenstern, 2003) and several studies (e.g. Moore et al., 2015; Steele-MacInnis et al., 2017; Hanyu et al., 2020; Lowenstern, 2003; Hartley et al., 2014) indicate that MI with large vapour bubbles (⁠|$>$| 10 vol %) are trapped from an already bubble bearing magma. Trapping magma that already has exsolved fluids is sometimes referred to as heterogeneous bubble trapping but in order to avoid confusion with heterogeneous bubble nucleation, we refer to this bubble trapping mechanism as exsolved volatile trapping and use the term ‘pheno-bubble’ to describe bubbles existing in magma chamber (Steele-MacInnis et al., 2017; Toramaru, 2014).

Here, we present a 3D X-ray computed tomography data set collected at Deutsches Elektronen-Synchrotron DESY of over 2000 MI. This large dataset allows us to link the macroscopic processes of a volcanic system with the microscopic textures of MI and put forward that looking at the textures of MI can provide additional precious information on pre-eruptive processes. This allows us to further our understanding of the relationship between melt inclusions, the magma reservoir processes responsible for forming them, and ultimately saturation state of deep seated magma and its role on explosivity.

Colli Albani

We present MI that are hosted in pyroxene (px) and leucite (lct) phenocrysts from the Colli Albani Caldera Complex. Colli Albani is located in central Italy 20 km SE of Rome and has a history of large-volume and explosive magmatism, making it unusual amongst other mafic-alkaline volcanoes (Giordano and the CARG Team, 2010). Colli Albani is currently in a period of quiescence, having not erupted in over 23 ka, though still exhibits sustained CO₂ degassing (⁠|$>$|4.2 × 109 mol year⁻¹) and uplift over the last 70 years (Trasatti et al., 2018; Chiodini et al., 2001; Todesco & Giordano, 2010; Chiodini et al., 2004; Giordano and the CARG Team, 2010; Chiarabba et al., 1997). The anomalous eruptive activity of Colli Albani provides a fascinating question for the volcanological community, and current literature suggests that elevated quantities of CO₂, combined with rapid magma ascent, contributes to the explosivity of Colli Albani magmas (Freda et al., 1997; Iacono Marziano et al., 2007; Bianchi et al., 2008; Mollo et al., 2010; Freda et al., 2011; Cross et al., 2014; Kleest & Webb, 2021; Jorgenson et al., 2024).

In this study, we investigate the Villa Senni eruptive unit (VSN) as it is the most recent of the large-volume ignimbrites at Colli Albani and is thus the best exposed. VSN erupted 18 km³ dense rock equivalent (DRE) at 355 ka. VSN is comprised of a basal fallout (VSN0) and two flow units: Tufo Lionato Ignimbrite (VSN1) and Pozzolanelle Ignimbrite (VSN2), which are separated by a co-ignimbrite breccia at some localities (VSN2b; Giordano and the CARG Team, 2010). We also analysed one sample which is from a pre-VSN fall deposit and one from Pozzolane Rosse, the largest ignimbrite of Colli Albani. A detailed study of clinopyroxenes of the same sample locations as this study has been done by Jorgenson et al. (2024) who suggests that the eruption of VSN is triggered by a rapid ascent of deep seated magma which blows through the upper crustal reservoir, as evidenced by high magnesium number (Mg/(Fe + Mg); |$>$| 0.8 and up to 0.94) and Cr₂O₃ (⁠|$>$|0.25 and up to 0.94 wt %) in VSN0. Further information on these units can be found in Giordano and the CARG Team (2010), Vinkler et al. (2012), and Jorgenson et al. (2024). By looking at MI textures from the VSN ignimbrite in 3D we are able to better understand the processes that lead to MI trapping, variability in quench and ascent rates, and the state of volatile exsolution, all of which help us to understand this unique system, as well as to draw general conclusions for the significance of MI in magmatic systems.

METHODS

Sample preparation and analysis

Samples were collected as bulk ignimbrite samples from the Villa Senni unit of the Colli Albani Caldera Complex. For more details on the sample locations, readers are referred to Jorgenson et al. (2024). Samples were then crushed and crystals of pxroxene and leucite were separated in an ethanol bath using a binocular microscope to hand pick mineral separates. Crystals were mounted in crystal bond, polished just below the crystal surface, and a transmitted light binocular and petrographic microscope were used to identify crystals with MI. Samples were scanned at the PETRA III beamline P05 (Wilde et al., 2016) operated by the Helmholtz Centre Hereon at Deutsches Elektronen-Synchrotron DESY. Crystals were mounted on a metal sample pin and adhered to the surface with UV hardening dental resin. In the experimental setup the full beam (ca. 6 × 3 mm²) illuminated the sample, 2500 projections were collected over an angular range of 360 degrees, around the longest axis. Different photon energies for the beam were tested, ranging from 13 to 40 keV. The upper bound was set to penetrate even the thickest sample (up to 4 mm), and the lower bound was set to provide the highest sensitivity in the full resolution depth of the camera for smaller samples.

Data processing and segmentation

Reconstructions, converting a set of 2D images into a 3D representation, were obtained immediately after data acquisition, with a binning factor of 2 (effective pixel size of 0.5–1 μm) exploiting a code by Moosmann et al. (2014) and van Aarle et al. (2016). Absorption contrast was sufficient to highlight MI in most samples with resolution on the order of 1 μm. Sample post processing samples was done with Image J (cropping) and Avizo software (segmentation and volume rendering; Schindelin et al., 2012; ThermoFisher Scientific, 2021). A median filter was applied to remove ‘salt and pepper’ noise on the image and reduce overlapping grey scale values. After denoising, a mask for the crystal volume was created. We then applied the ‘interactive thresholding’ tool to segment the inclusions. In many cases, interactive thresholding was not sufficient for accurate segmentation, and in this case, data were manually processed using the ‘brush’ tool. Then each inclusion was assigned a label using the ‘labelling’ tool. Vapour bubbles were added by using the ‘fill holes’ tool and greyscale threshold of the vapour phases (needed for VB on the edge of MI). Using the ‘arithmetic’ tool, we subtracted the filled MI and empty MI to get only the vapour phase. Finally, we used the ‘label analysis’ tool on both the filled MI and the vapour bubble. Parameters measured include volume, area, shape factor, length, width, breadth, thickness, and mean integral curvature. Length is the maximum distance across the object, whereas width is the minimum. Critically, because 3D geometries were often irregular, these two measurements are not necessarily orthogonal to each other. Breadth is defined as the largest distance between two parallel lines which touch the object but do not intersect it and is orthogonal to the length. Thickness is the largest distance that lies in an orthogonal plane to the length and breadth and touched the end points of the object. The shape factor is calculated as:

where 1 is a perfect sphere (ThermoFisher Scientific, 2021). From this data table, we visually inspected each label for artefacts, edges, and fractures. Afterwards, MIs were categorized, their number of bubbles and phases was determined, and we recorded any other pertinent information. We note that inclusions on the crystal edge are incredibly difficult to measure accurately with our technique, as the grey scale value for the vapour bubbles and air are the same so we do not include any MI along crystal edges in this analysis.

Segmentation of intracrystalline zoning with sufficient phase contrast was performed following the methods of (Lubbers et al., 2023) and a small python module built on top of the package scikit-image (Pedregosa et al., 2011): CTPy available at https://github.com/jlubbersgeo/ctpy. In brief, the workflow consists of the following steps: loading data in the form of a stack of 2D images, resampling the data to fit memory limits of the user’s computer, normalizing the data for processing, contrast stretching, denoising with a non-local means filter, and zone segmentation with a watershed algorithm. Readers are referred to Lubbers et al. (2023) and github documentation for further explanation.

Limitations

Although we benefit from a very high imaging resolution, we still need to consider partial volume effects (Ketcham & Mote, 2019; Kato et al., 2013). More specifically, since each 3D dataset is composed of voxels, it represents a discretized version of a real 3D volume. If the boundary between two materials is not at the exact edge of a voxel the resulting greyscale value will be a mixture of each material. While methods of sub-resolution feature quantification do exist (Ketcham, 2019), our greyscale thresholding used for segmentation and manually corrected afterwards makes it difficult to quantify the variability in volume from including and excluding the partial volume. We estimate the variability in the manual segmentation by repeating the segmentation of a scan three times. The results of this are outlined below. Following the works of Spavieri et al. (2018), we determine the minimum detectable size to be 1 μm (double the voxel size). The minimum quantifiable volume we determine is 15 μm³, based on a minimum quantifiable diameter which is five times the voxel size (Gualda & Rivers, 2006).

Segmentation reproducibility

Ideally, segmentation via thresholding alone would give reproducible segmented MI, but given the variability in the contrast of some inclusion, many needed manual corrections. Consequently, the segmentation is not automated and subject to user variability. To understand the reproducibility of our manual thresholding segmentation method, we repeated the segmentation and labelling process three times for a single scan—scan 0017. Scan 0017 is a px with 30 enclosed inclusions that cover a range of textures representative to our overall dataset. Due to poor phase contrast between host and inclusion, all 30 inclusions were not segmented each iteration. The average MI volume variability per inclusions is 428 and 30 μm³ for the VB. As the variability of MI volume ranges several orders of magnitude, we use the relative uncertainty to assess the measurement reproducibility. We calculate the percent uncertainty as:

where the range is the difference between the maximum volume and the minimum measured value.

The mean percent uncertainty is 10.5% for MI volume, 7.7% for VB volume, and 17.9% for the volume fraction (VB volume and volume fraction were only calculated on glassy inclusions). As shown in Fig. 1a, the percent uncertainty varies with MI volume. This result is expected as variation of a single voxel has a greater effect on an object with a total of 50 voxels versus 50 000 voxels. The mean percent uncertainty for MI volumes >1000 μm³ is 4% and for <1000 μm³ is 15%. For the vapour bubbles, the percent error is less systematic, perhaps due to the already small volumes. The percent error of G1 type inclusions (see section Glassy MI with a single VB; n = 17) is 9.2% for the MI, 7.5% for the VB, and 18.2% for the volume fraction. For GM type inclusions (See section Glassy MI with multiple VB; n = 6), the error is 6.9% for the MI, 6.9% for the VB, and 17.3% for the volume fraction. For G (See section Glassy MI; n = 6) and MC (see section microcrystalline inclusions; n = 1) type inclusions, the MI volume percent error is 20.6% and 3%, respectively.

RESULTS

We analysed and processed a total of 2045 MI from 35 crystals. Px and lct phenocrysts are from the VSN0, VSN1, and VSN2 subunits, as well as one px from Pozzolane Rosse and a fall deposit before after VSN. We scanned and segmented 10 samples from VSN0 (5 px and 4 lct, and 1 apatite), 14 samples from VSN1 (9 px and 5 lct) and 9 samples from VSN2 (6 px and 3 lct).

Type of MI

We observe a wide range of inclusion textures and also find this variability exists on both the individual crystal and sub-unit scale. Based on this observed variability, we separate inclusion into six types, broadly defined by three metrics - degree of crystallinity, shape, and distribution of vapour bubbles. We note that though these classifications are an oversimplification of MI complexity they aid in interpreting the link between MI textures, type of deposits, and the eruptive sequence.

Glassy MI (G)

The term glassy refers to a single homogeneous glass phase which we see represented in our scans as no variability in phase contrast within the inclusion. This inclusion type includes all MI, which are glassy and VB-free (G in Fig. 2). Spherical high density blebs are also common in glassy MI, which we assume are sulphide. Inclusions that are otherwise homogeneous other than the sulphide bleb (and VB for the other MI types) were consider to be glassy. Of all px hosted MI analysed, 33.0% are glassy and bubble-free; of the leucite-hosted MI measured, 15.6% are type G. We note that these inclusions may have a vapour bubble that is not resolvable as our maximum resolvable object is ~1 μm³ (Spavieri et al., 2018).

Glassy with single VB (G1)

These MI are glassy with a single spherical vapour bubble (Fig. 2, G1). Of the px measured, 30.6% are type G1; and of the leucites, 18.4% are type G1.

Glassy with multiple VB (GM)

These are glassy with more than one spherical VB. We found up to 159 bubbles in a single MI. These vapour bubbles are only found on the rims of the MI and never in the centre of the inclusion (Fig. 2, GM; note the bubbles are on the top rim of the MI). Of the px measured, 2.6% are type GM; and of the leucite host MI measured, 17.2% are type GM.

Glassy with irregular VB (Gi)

These are glassy MI with on or more vapour bubbles that are irregularly shaped (Fig. 2, Gi). We consider irregular shaped to be anything other than a spheroid or ellipsoidal shape. Of the px measured, 8.3% are type Gi; and of the leucites measured, 4.7% are type Gi.

Tube inclusions (T)

These are glassy inclusions with elongated tube or blade shapes (Fig. 2, T). We exclude these from the previous MI types as the variation in shape might be indicative of a different entrapment process. We define these MI to be tubes visually. The average length over breath is 5.5, and the average length over thickness is 9.7. Of the px measured, 13.6% inclusions are type T; and of the leucites measured, 7.0% are type T. We note that it may be possible that some T inclusions are not actually melt but crystals. Apatite and phlogopite inclusions have been chemically identified in several crystals. However, it is hard to definitely say via textural evidences (and especially without clear crystal habit), so all are considered melt inclusions for this manuscript.

Microcrystalline inclusions (MC)

These inclusions have one or more mineral phase within the glass portion (Fig. 2, MC). This melt type shows the most variability in vapour bubble distribution and shape, as well as the melt inclusion shape. Of the px measured, 11.9% are type MC; and of the lct, 37.1% are type MC. In many cases, but not all, these inclusions were connected to a crack in the crystal.

MI location within crystal

One way to look at the relationship between a MI and host is the respective distance from the core to rim and where it is with respect to zoning patterns. Several of our crystals were partially polished to ensure MI presence and to optimize the experimental setup, so we cannot investigate the MI location this way. However, in four crystals (px: 0018 and 0019, lct: 0077 and 0062), a clear core and rim is defined, and the MI is preferentially located near the core. There is no strong correlation between melt inclusion type and location of the inclusions in the crystal centre.

Additionally, there are many instances where MI are oriented along a plane in the host phase (scans 0023, 0021; 0018; 0019; 0048; 0010; 0007; 0064; 0056; 0017 in supplementary file S2). In some cases, this is along crystallographic orientation, and is often the case with T type inclusions. However, there are instances where there is a plane of MI in the crystal centre not associated with the crystal axis, and we suggest that these are from cracks which are later filled by melt. In such cases there, is no preferential melt inclusion type found along these planar alignments. Planar alignments of inclusions not oriented along crystallographic axes are similar to secondary assemblages in fluid inclusion studies (Roedder, 1979). Thus, we suggest that MI oriented along noncrystallographic planar features should be avoided as they might not be representative of the melt in equilibrium with the host crystal.

An alternative to looking at the relative location of the MI is to investigate MI location within crystal zones. Measurements were optimized to see contrast between melt and host phase and not for zoning of the host phase so not all crystals showed clear zoning patterns in these scans. However, we managed to successfully segment clear zones in seven px crystals and one apatite crystal (scans 0018; 0019_x2; 0019_x3; 0023; 0039; 0096; 0100; 0053, see supplementary file S2). Zones correspond to regions of different density, and correlating these zones with the geochemistry of the host phase is out of the scope of this work, but is a possible avenue for future work. Five px have very distinct core and rim and a concentration of MIs in the crystal core. The apatite and two px show zoning and a preference for melt inclusions to be in a single zone. In all scans, MI preferentially reside in one zone. We further discuss the melt inclusion relationship with the varying zone below.

There is no distinct spatial relationship between the MI type and location in the crystal (aside from T type inclusions as discussed above), so one cannot determine if an MI is primary or secondary from looking at a single inclusion. This stresses the need for 3D data to be obtained for MI studies.

MI type and volume

Figure 5 shows a relationship between MI type and size. G and G1 type inclusions have the smallest mean volumes, Gi and T type have intermediate volumes and GM and MC inclusions have the largest volumes (Table 1). Notably, leucites MI are generally smaller than px MI (mean volume of 2432 |$\pm$| 11% μm³ versus 6937 |$\pm$| 11% μm³ for the pyroxene), and have less variability in the distribution of the MI volume with a mean skewness and stdev of 8.2 and 8862 |$\pm$| 11% μm³ for the leucite and 9.9 and 32 948 |$\pm$| 11% μm³ for the px (Table 1). Within each individual crystal, there are many instances where the inclusions follow the same trends. The Gi and MC type inclusions are the largest, and the G and G1 types are the smallest (Fig. 5 and Supplementary file S2 scans 0007; 0018; 0056; 0064; 0096; 0100). GM types are often larger than G and G1 type inclusions (Fig. 5 and Supplementary file S2 0017; 0064; 0013; 0064). However, there are some crystals with notable exceptions (e.g. scans 0010; 0031; 0048 in supplementary file S2) where G1 inclusions are larger than GM inclusions.

Shapes

Shape is a difficult parameter to quantify outside of classical geometric shapes. As clearly seen in the 3D renderings, MIs are generally ellipsoidal in nature but can deviate from this shape (Figs 3 and 4 and Supplementary file S2). Fluid inclusion research separates inclusions on the basis of regularity (defined as a ratio between the perimeter squared divided by the area times 4|$\pi$| where above 1.75 a shape is ‘irregular’ and otherwise ‘regular’) and elongation (defined by a ratio between the major and minor axes where a shape is equant until this ratio is above 2; Bakker & Diamond, 2006). Our MI classification system, while more focused on VB distribution, utilizes elongation to classify T type inclusions. Instead of looking at regularity we can compare MI and VB shapes with shape factor and curvature metrics. We use shape factor to look at the shape of the MI (Equation 1). An MI with a shape factor of 1 is a perfect sphere and using the (Bakker & Diamond, 2006) shape classification system would be a ‘regular-equant’ inclusion. Figure 6 shows the shape factor of all the MI types. We also note that deviations from a perfect sphere may not only be due to the overall inclusion shape, but can also be due to surface roughness (e.g. a dimpled golf ball versus a smooth pingpong ball). Results of the shape factor are shown in Table 2.

As the T type inclusions were chosen on a shape basis, predictably none are close to the spherical value of 1 (Fig. 6 and Table 2). MC inclusions also deviate from the perfect shape factor of 1 (mean shape factor is 3.1 and 2.3 for px and lct; Table 2). GM and Gi type inclusions start to approach 1 with the mean px of 1.8 and 2.4 and leucite 1.3 and 6.8 respectively. G and G1 inclusions are the closest to a perfect sphere where G and G1 for px has a mean shape factor of 1.9 and 1.6 and 1.5 and 1.7 for leucite (Table 2). Between crystals, we clearly see px inclusions deviate further from a sphere than lct MI (Fig. 6). Notably, in some of Gi, GM, and MC inclusions, the inclusion is two or more interconnected spheroids, which maybe be evidence of a ‘necking-down’ process as described in Yao et al. (2020).

Curvature is the measure of how abruptly a curve deviates from a straight line, or can be considered as how much an object varies from being flat. Mean curvature is defined as the arithmetic mean of the two principle curves (Crane et al., 2013), H = (⁠|$\kappa$|1 + |$\kappa$|2)/2. The integral of the mean curvature for 3D objects is computed in Avizo software using this equation:

where |$\varSigma$| 3d_i is the sum of measures in a local 2 × 2 × 2 environment. Mean curvature does not take into account surface roughness, just the overall global geometry. However, it offers us some insights into the complexity of an object. Objects with an integral of mean curvature equal to 0 have a minimal surface area, such as a spherical droplet in a vacuum or the helical shape of our DNA (Gennes et al., 2004). In curvature space, we see that almost no inclusions have a minimal integral mean curvature, and thus none have a minimal surface energy. The px inclusions seem to show more variance than the leucite MI. G and G1 type inclusions seem to systematically have a smaller curvature than the other inclusion types. We suggest the more complex shapes as we empirically see are represented here in the curvature.

Vapour bubble volume fraction

To investigate the relationships between VB and MI volumes, we focus on G1 and GM type inclusions, which are less likely to be affected by post entrapment processes than Gi and MC type inclusions. Figure 7 shows that px have more G1 type inclusions than GM type inclusions where lct shows a more equal distribution. In terms of VB volume, G1 and GM VB hosted in px have larger VB than MI hosted in lct (px G1 mean volume 124 |$\pm$| 8% μm³ GM mean volume of 558 |$\pm$| 8% μm³ and lct G1 mean volume of 96 |$\pm$| 8% μm³, GM mean volume of 249 |$\pm$| 8% μm³; Table 1).

Both px and lct hosted MI have a wide-volume fraction range (px, 0.4–82 vol %; lct, 0.3–78 vol %). Lct G1 type inclusions have a slightly larger volume fraction than px (px mean vol frac = 0.14 |$\pm$| 15%, lct 0.16 |$\pm$| 15%; see Table 1). Px hosted GM inclusions have almost the same average volume fraction (Px mean vol frac = 0.062 |$\pm$| 15%, lct = 0.060 |$\pm$| 15%; see Table 1).

Vapour bubbles can represent both volume fraction trapping pheno-bubbles and/or post-trapping exsolution, thus we must determine what volume fraction is from an exsolved melt and what are the other factors. Bubble formation is directly linked to volatile solubility, which varies with pressure (Papale et al., 2006; Iacono-Marziano et al., 2012; Moussallam et al., 2015). Changes to MI internal pressure can occur from differential thermal contraction (DTC) between the host phase and MI and from postentrapment crystallization (PEC; Steele-MacInnis et al., 2011). PEC changes MI compositions by crystallizing host phase on the MI rim, and thus depletes the MI in host phase elements which changes the solubility of volatiles in the melt (Steele-MacInnis et al., 2011). DTC can be calculated using the Moore et al. (2015) calculator. They calculate the relative volume change of the glass (V/V0) as the reciprocal of the melt density normalized to the density at the trapping temperature. The volumetric change of the host phase is calculated using empirically derived thermal expansion of the host minerals, and the calculated volume proportion of the bubble is thus the total difference between the thermal contraction of the inclusion and the host. We use a composition based on GEOROC data of SiO₂ 45.17 wt %; TiO₂ 1.028 wt %; Al₂O₃ 14.183 wt %; FeO 8.676 wt %; MnO 0.143 wt %; MgO 5.65 wt %; CaO 10.689 wt %; Na₂O 1.323 wt %; K₂O 6.884 wt %; H₂O 3 wt %; CO₂ 0.8 wt %; a trapping temperature of 1200°C which is one of the maximal temperature estimates from Jorgenson et al. (2024); and a glass transition temperature of 748.4°C from Giordano et al. (2008). This results in a 4.7 vol % difference for clinopyroxene and 5.3 vol % for alkali feldspar, which we use as a proxy for lct. If we assume a 2% volume difference from PEC (Hanyu et al., 2020) then the total volume that can be from volatile exsolution postentrapment is 6.7 vol % for px and 7.3 vol % for lct. Using these limits, we find that in total, 62.6% of the G1 and GM type inclusions in px are above the limit, whereas for lct, 46.8% of the inclusions are above the limit. Data for all melt inclusion metrics measured can be accessed in Supplementary Table S1.

DISCUSSION

How to form each MI type?

While MI can provide a wealth of information regarding pre-eruptive processes and magma storage, they should be carefully selected to ensure inclusions are representative. Here, we consider each MI type and consider its provenance with respect to the aforementioned trapping mechanisms (Roedder, 1979). We stress that while textural variability of a single melt inclusion may reveal information about trapping, it is also important to place the melt inclusion in the spatial context of the host crystal. While MIs are often assumed to be primary, the identification of planar alignments of MI that are not aligned with crystallographic axes implies that it is possible to have secondary trapping of MI. Fluid inclusion studies rely heavily on the spatial context of the host, grouping inclusion in assemblages and separating them into primary or secondary trapping mechanisms (Bodnar et al., 2006), and we encourage melt inclusion studies to adopt these practices, as also suggested in Rose-Koga et al. (2021).

G type

G type inclusions are the smallest inclusions on average (Fig. 5). We suggest that these inclusions are less likely to be captured via a disequilibrium process and are more likely to be captured where crystal defects develop, owing to their small size (Roedder, 1979). These inclusions may form without a vapour bubble since their small volume favours rapid quenching before there is time for vapour bubble exsolution, and there is a smaller probability for vapour nucleation. Additionally, smaller inclusions have a larger internal pressure from surface tension, and thus the ability for volatiles to exsolve is smaller (Roedder, 1979; Tait, 1992). It is also possible that the lack of bubbles in these inclusions is due to a lack of volatiles in the melt or a VB that is too small to be resolved by our measurements. Spavieri et al. (2018) suggests that minimum detectable size of an object is double the voxel size, which in our scans is 0.49–0.98 μm, thus any bubble below 1 μm in diameter is not resolvable.

G1 type

Whether a MI is trapped from a melt with exsolved or dissolved volatiles can be gleaned from VB volume fraction (Moore et al., 2015). To simplify variable VB distributions within an MI, we consider several endmember scenarios, as outlined in Fig. 8, noting that this does not encompass all processes that can affect VB exsolution (Rasmussen et al., 2020; Aster et al., 2016; Danyushevsky et al., 2002; Maclennan, 2017; Bucholz et al., 2013).

Considering MI trapping in presence of excess volatiles, we show a situation in Fig. 8.1a where the MI traps several pheno-bubbles exsolved in the melt. The MI is then enclosed by the crystal and quenched with the bubbles in place. Here, we expect many bubbles randomly distributed throughout the MI and the time between exsolution and quenching to be rapid as otherwise bubble coalescence would occur, as is the case in Fig. 8.1b. In the case of an exsolved magma with bubbles which are large relative to the MI, (Fig. 8.1c), the volume fraction of the VB to MI is not representative of the rest of the magma. We note in the case of large pheno-bubbles in an exsolved melt, it is also possible that MI may not trap a bubble at all, and thus, the resulting MI would be volatile-free and not representative of the initial magma.

In the case of MI trapping from a melt with dissolved volatile, the generation the VB subsequently exsolves upon a pressure or temperature change of the MI. MI are shown to retain high pressures even when the external pressures are low (Steele-MacInnis et al., 2017), thus we suggest that it is most likely that exsolution happens due to a pressure drop during ascent in tandem with PEC + DTC. Heterogeneous bubble formation on the rims of the inclusion is likely the driving exsolution mechanism, as homogeneous nucleation requires a much larger over pressure (Gardner et al., 2023). Heterogeneous trapping would utilize the MI-crystal interface to nucleate either one (Fig. 8.2d) or several (Fig. 8.2a–c) bubbles on the inclusion rim. The nucleation of several small bubbles and subsequent coalescence (Fig. 8.2a) is seen in experiments from Mangan & Sisson (2000), Hanyu et al. (2020), and Drignon et al. (2021). Coalescence is a fast but not an instantaneous process, thus it is feasible to trap bubbles in the middle of coalescence as in Fig. 8.2b and further discussed below. However, it is possible that bubbles may not coalesce at all (Fig. 8.2c), which may happen due to either a too high activation energy to start to coalesce (bubbles to far apart or surface tension too high) or quenching that is too fast to initiate coalescence. Lastly, in a dissolved volatile-rich magma, a single VB may nucleate (Fig. 8.2d).

Clearly, G1 inclusions can be generated from trapping either a bubble bearing (Fig. 8.1b or c) or bubble-free melt (Fig. 8.2a or d) and to determine this we can evaluate the vapour volume fraction, further discussed below. We note that shapes of G1 inclusions are commonly close to spherical (px, mean of 1.6 and mode of 1.0; leucite mean of 1.7 and mode of 5.1), suggestive of a trapping mechanism which allows for a minimal energy and thus is more likely to be in equilibrium with the host. These inclusions are great candidates for measuring volatiles as they have only a single vapour bubble and are often larger than the G type inclusions (px mean G1 volume = 1609 |$\pm$| 11% μm³ and G = 1489 |$\pm$| 11% μm³, lct mean G1 = 4201 |$\pm$| 11% μm³ and G = 522 |$\pm$| 11% μm³).

GM type

Multiple VB are not commonly discussed in geochemical studies of MI, perhaps owing to the difficulties of measuring and reconstructing volatile contents from many bubbles, yet they are present in nature (Frezzotti, 2001; Steele-MacInnis et al., 2017; Cannatelli et al., 2016; Wallace et al., 2003; Rose-Koga et al., 2021) and multiple bubbles have been reported in MI re-homogenization experiments (Hanyu et al., 2020; Pintea, 2013; Drignon et al., 2021). Notably, Pintea (2013) calls these melt inclusions ‘foamy’. We find that GM inclusions are larger relative to G and G1 inclusions (Fig. 5) and are slightly more irregularly shaped (Fig. 6). This is more pronounced in px than lct hosted MI (GM px mean shape factor = 1.8, mode = 3.1; lct mean shape factor = 1.3, mode = 2.1). This finding is also reflected in the work of Yang & Scott (2002) who reports olivine and px hosted MI with multiple VB (which occupy 5–40 vol % of the MI) are often larger and more irregular than other MI. From this observation we suggest irregular shaped GM inclusions are more likely to be formed during a period of rapid crystal growth.

Frezzotti (2001) suggests multiple bubbles are trapped from a melt with exsolved fluids (as in Fig. 8.1a), yet experiments show multiple bubbles can form from heterogeneous nucleation and coalesce into a single bubble (Ohashi et al., 2022a; Drignon et al., 2021; Hanyu et al., 2020; Pintea, 2013). There are no GM inclusions with bubbles in the inclusion centre, pointing away from the hypothesis that all GM inclusions are trapped from a bubble bearing melt. Some of GM inclusions of this study show mid-coalescence bubbles connected by a thin neck (Fig. 9). We suggest that these MI bubbles have exsolved and began coalescence but have passed the glass transition temperature and quenched before coalescence was completed. In our samples we have many instances of bubbles trapped mid-coalescence (Fig. 9 and Supplementary File S2), however we note that films between two bubbles can be small and may not be resolvable with our technique (Castro et al., 2012; Nguyen et al., 2013; Ohashi et al., 2022a).

Coalescence is a three-step process of bubble approach, film drainage, and shape relaxation to a spherical shape (Toramaru, 2022). Bubble coalescence has been well studied in the context of bubbles in open system magmas and is shown to be a rapid process (Masotta et al., 2014; Ohashi et al., 2022a; Ohashi et al., 2022b; Castro et al., 2012). Masotta et al. (2014) shows experiments with bubble growth in basaltic melts on the timescales of seconds to minutes and Nguyen et al. (2013) finds that film drainage of a low-viscosity magma is on the order of magnitude of tens of seconds or even less. Toramaru (2022) shows theoretically that shape relaxation is on the order of seconds, which is also seen empirically by Masotta et al. (2014) and Ohashi et al. (2022a). Thus, all stages of the bubble coalescence process take seconds to minutes. This is confirmed by re-homogenization experiments of Hanyu et al. (2020), who found that 5-minute re-homogenization experiments of olivine hosted MI with multiple bubbles merged to become one single bubble (experiments ran at 1180–1280°C). We suggest that the presence of GM inclusions could be indicative of rapid quenching. However, GM inclusions may also be controlled by inclusion volume. GM inclusions are systematically larger than G1 inclusions (px mean GM volume = 16 866 |$\pm$| 11% μm³, px mean G1 volume = 1609 |$\pm$| 11% μm³, lct mean GM volume = 6230 |$\pm$| 11% μm³, lct mean G1 = 4201 |$\pm$| 11% μm³). An alternative to fast quenching is that larger volume inclusions are less likely to have bubbles close enough to approach each other and thus VB remain attached to MI walls without coalescing.

Proportionally, there are more GM inclusions below the DTC + PEC limit than G1 inclusions which we propose indicates nucleation of many small bubbles on the MI rim is favoured over nucleation of a single VB during rapid contraction. Proportionally, lct hosts have more GM inclusions than px. Notably, above 1000 μm inclusions with a large vapour fraction are absent, which we at suggest is because large MI have a higher change to decrepitate or crystallize.

Gi type

Many MI measured have irregular VB shapes, and were classified as Gi type inclusions. We suggest that irregular bubble shape is due to decrepitation, when the inclusion ruptures and loses CO₂ or H₂O to the external melt. This is a major process controlling the distribution of CO₂ and measuring decrepitated MI may result in CO₂ underestimations (Maclennan, 2017). Whether or not a MI decrepitates is ultimately due to pressure difference between the inclusion interior and exterior inclusion. Inclusion pressure is linked to size and shape, and smaller inclusions (diameters |$<$|10 μm) can reach higher pressures (⁠|$>$| 300 MPa) without decrepitation (Wanamaker et al., 1990; Campione et al., 2015). We do not see strong dependency on inclusion volume (Fig. 5). However, the stress state between a MI and host phase also plays an important role in if decrepitation occurs. Tait (1992) found that non-spherical (cylindrical) MI or MI with corners and irregular points will produce stress concentrations and make them more susceptible to decrepitation. Thus, upon ascent (in an isothermal regime), the change in pressure between MI and host phase causes induces cracks within the host phase leading to volatile loss and decrepitation. We suggest that Gi type MI are formed from MI decrepitation and, therefore, are not representative of the true melt or volatile content of initial melt that formed them. Naturally, these are best avoided for volatile studies.

MC type

MC inclusions are the largest of all inclusions types (px MC mean volume = 26 600 |$\pm$| 11% μm³, lct MC mean volume = 1348 |$\pm$| 11% μm³) and we suggest this follows the relationship between cooling rate and volume proposed by Roedder (1979) which calls these melt inclusions ‘mixed’. The time needed to quench a large inclusion is longer, thus allowing time for crystals to nucleate and grow within the inclusion. Experiments by Bodnar et al. (2006) found H₂O saturated conditions may also promote crystallization of a MI during cooling. Leucite crystals have the highest proportion of MC inclusions (px is 11.9% where lct is 34.7%). We suggest this could be a function of the shape or location of MC inclusions, which are frequently found in the centre of large crystals (scans 0077 and 0062). While MC are certainly harder to reconstruct geochemical information from, we note the works of Créon et al. (2018) who utilizes X-ray micro-tomography to reconstruct MC compositions.

T type

T type inclusion volume is variable, indicating that the main control of these inclusions is host crystal shape and habit and not cooling rate. The orientation of these inclusions is often along the crystallographic axis (see supplementary file S2). We suggest T type inclusions could be generated from cracks along cleavage planes in the crystal that have filled and healed at a later stage. In this instance, the melt phase is not representative of the melt the host crystal grew from. It is also possible that these MI formed via decrepitation, where a MI causes a crack as explained for the Gi type inclusions. This is evidenced by a mostly tube shape but with a thicker blob section of the inclusion. T type inclusions may also be closed embayments, of whose use and morphology is further discussed in Ruefer et al. (2021), Hosseini et al. (2023), and others. Additionally, these MI have the least amount of gas, which we suggest is due to melt infill from a primarily degassed magma.

Crystal zoning around MI

Crystal zoning is an important consideration in geochemical studies, so it stands that zoning should also be considered in MI studies (Ruth et al., 2018). While sample preparation for the synchrotron required partial polishing the crystal away to ensure MI presence, there are some samples which have resolvable preserved crystal rims and zoning from phase the contrast scans. Scans 0018, 0019_x2, 0019_x3, 0023, and 0096 are all well-preserved samples, with a clear crystal habit and a distinctly zoned core and rim (Fig. 10). It is important to note that while the observed zoning is due to density variation in the crystal, we have not quantitatively correlated with geochemistry at this moment in time. In these scans, it is clear that MI are predominately found in the core, irrespective of MI type. In some crystals where crystal habit is not preserved and there is no clear core and rim zones, we still find zoning and a preferred zone where the MI reside (scan 0039, 0100, and 0053; Fig. 10). Previous studies have denoted MI assemblages as zonal or azonal and suggest that zonal inclusions are always primary and azonal inclusions are only primary if there is no evidence of fracturing or mineral dissolution (Bodnar et al., 2006).

Notably, all zoned samples with an MI-rich core are from VSN2. The other samples are from VSN0 (0100 and 0053) and RED (0039). While we do not have enough crystals for a representative view, we can begin to speculate on what the zoning may indicate. The increased abundance of MI-zoned crystals in VSN2 versus VSN0 and VSN1 may be from differing processes in the magma reservoir that make the crystals of VSN2. One way to generate the MI-zoning patterns of VSN2 crystals would be to have magma rise quasi-adiabatically from a depth at which a magma is at a subliquidus temperature and exsolves no or very modest quantities of H₂O upon ascent, causing the magma to become overheated, and the crystals to partially resorb. The resulting irregularities on the mineral rims can be conducive to melt trapping, followed by a period of rapid crystal growth induced by water exsolution at shallow depths (i.e. equivalent to 200–100 MPa, Jorgenson et al., 2024). This also indicates a pause in the magma at shallow levels to allow time for resorption and subsequent crystal growth. We suggest the variation in zoned crystals in VSN1 and VSN0 versus VSN2 is indicative of varying residence time in the upper crustal reservoir.

MI and VB volumes from 2D data

Reconstructing CO₂ from vapour bubbles requires accurate volumetric measurements of both MI and VB. Propagated errors from volume measurements can result large errors recalculated CO₂ (Hartley et al., 2014; Bakker & Diamond, 2006). A common assumption used for volume estimates is that MI and VB are ellipsoidal; however, our results show that this is not always the case and average shape factor varies with MI type (Fig. 6). Bakker & Diamond (2006) show large variation in inclusion size and vapour fraction with varying inclusion shape, and stress the need for accounting for a third dimension when assessing vapour fraction of irregularly shaped fluid inclusions which they do utilizing a spindle stage. The shape of MI in this study are both regular and irregular where G type are the most spherical inclusions, where G1 and Gi types are also commonly spherical. GM and MC inclusions are more likely to be not spherical, and clearly T inclusions are not spheres. Here, we investigate variability of volume estimations assuming an ellipsoidal MI compared to the 3D volume estimations, allowing for a greater understanding of the limits of using 2D data for 3D volume reconstruction and to provide a guide for future studies.

To test these assumptions, we use G1 type inclusions as they are what are commonly used in MI studies. To calculate the volume we assume an ellipsoidal shape and use the length (longest measurement from Avizo) as the a axis and the breath (orthogonal to a) as the b axis. To determine the c axis we calculate it from the a and b axes using five different methods 1) c = b, 2) c = arithmetic mean, 3) c = geometric mean, 4) c = true thickness (the third orthogonal measurement from Avizo) and 5) c = b = a. As we can see in Fig. 11a and b, volume estimates using a calculated c axis and assuming a ellipsoidal shape are poor and overestimate volume, with errors of 100% or even 1000%. We suggest that these over estimations are due to deviations from perfectly ellipsoidal nature and from using the absolute maximal length and breadth estimates. The best fitting technique was to assume a perfect sphere with all axes equal to the breath. For the px this gives a mean, median, and modal % error of 67, 29, and − 100%, a std deviation of 430, and a skewness of 13. For the lct it gives a mean, median, and modal % error of 69, 55, and 20%, a std deviation of 67, and a skewness of 3 (also see supplementary table S3). The volume calculations using the measured thickness as the c axis predictably gives better results than recalculated values (Px: mean % error of 53, median % error of 38, modal % error − 99, std dev of 112, skewness of 11. Lct: mean % error of 72, median % error of 55, modal % error 79, std dev of 48, skewness 1).

Despite the wider distribution of error for the sphere estimation (⁠|$b=a=c$| axis), it gives best performance in volumetric estimates for both px and lct and, therefore, we recommend this for volume reconstruction from 2D images. However, we note this is only based on samples from Colli Albani and this method may be more robust in other systems. Vapour bubbles are more spherical than MI (Fig. 11), but their smaller size exacerbates volume reconstruction errors. As such % errors for VB are much larger than for MI (Fig. 11). Similar to MI best estimates are from assuming a sphere (⁠|$b=a=c$|⁠; Px: mean % error 166, median % error 48, modal % error − 100, std dev 589, skewness 8. Lct: mean % error 122, median % error 73, modal % error − 1, std dev 135, skewness 2). Again, we suggest MI volumes are recalculated using this method if only 2D data are available.

Overall we find that volume estimates from 2D sections are associated with large errors, with the lowest median percent error of 29%. We compare our findings with that of Tucker et al. (2019) who explores volume estimates using computer generated ellipsoids and measuring based off randomly intersected planes. They find that using the visible axes arithmetic mean to give the best estimate of true volume compared to the geometric mean and |$c=b$|⁠. However, their computed relative error for all methods can be very high; 36% on average for |$c$| equal to smallest axis, 5% for arithmetic mean, and −9% for geometric mean. Based on our findings and the findings of Tucker et al. (2019), we suggest that when possible 3D data are obtained for studies analysing CO₂. While 3D data are not always possible, there are other methods on may consider which account for the z-dimension including the spindle stage method of Bakker & Diamond (2006), which also accounts for inclusion shape, alternatively if recalculating from microscope images one could also consider a volume correction based on host crystal and phase measured (VB or MI).

Melt inclusions of Colli Albani reveal the evolution of bubble bearing magma

Proportions

Proportions of MI types between different subunits for both px and lct varies considerably (Fig. 12). The most striking variability is the proportions of px hosted T and MC inclusions in VSN2 compared to VSN0 and VSN1. Of the VSN2 px inclusions 22% are T type, where VSN0 and VSN1 have no T type inclusions. This contrasts with lct hosted MI where VSN2 has no T type inclusions and VSN0 and VSN1 are 10% and 6% respectively. The proportions of MC inclusions follow a similar trend, with the highest proportion in VSN2 px (18%, VSN0 4%, VSN1 3%) and higher in VSN0 and VSN1 in the leucite crystals (VSN0 22%, VSN1 54%, VSN2 18%). Generally lct crystals have a higher proportion of MC inclusions, inclusions which are associated with a slower cooling rate (Kent, 2008; Roedder, 1979). Within the px host it seems that a higher proportion of T type inclusions is coupled with a higher proportion of MC inclusions, and we suggest that the processes to trap these inclusions are linked. Additionally, some VSN2 crystals have a MI-rich core and a MI-free rim, where we do not see this MI zoning style in the VSN0 or VSN1 crystals (Fig. 10). The VSN2 px represent an endmember of MI assemblages, generated by slower cooling rates as evidenced by the high proportions of MC inclusions. Furthermore, these crystals have experienced resorption and rapid growth as evidenced by the MI free rims.

VSN0 px has a larger abundance of GM inclusions than VSN1 and VSN2. The proportion of GM inclusions in lct is larger than px, perhaps due to variability in trapping conditions and cooling rates. Shape factor of GM inclusions is more irregular than the other glassy inclusions, and leucite hosted GM inclusions are even less spherical than px hosted inclusions. Lct crystals of VSN commonly have a skeletal texture from rapid growth, which may encourage trapping of non-spherical MI and inhibit VB coalescence (Giordano and the CARG Team, 2010; Vinkler et al., 2012). Lct and px have similar thermal diffusivity, but px is slightly larger indicating a faster cooling rate (Kanamori et al., 1968; Hofmeister & Pertermann, 2008; Hofmeister & Ke, 2015), the variation is minimal and crystal size and shape also influences cooling rate (Vollmer, 2009). The presence of GM inclusions is indicative of a rapid MI quenching, which is also reflected in the proportion of VSN0 G type inclusions, the fastest cooling inclusions (Wallace et al., 2003; Roedder, 1979). VSN0 lct and px have the largest proportion of G inclusions (lct: VSN0 45%, VSN1 1%, VSN2 7% and px:VSN0 56%, VSN1 26%, and VSN2 33%). VSN0 crystals are from a fall deposit and VSN1 and VSN2 come from ignimbrite deposits. Fall deposits cool faster than ignimbrite deposits, which may drive the variability in MI proportions between VSN0 and VSN1/VSN2 (Thomas & Sparks, 1992; Wallace et al., 2003; Trolese et al., 2017). We can speculate that variation in ascent rate may also influence quench rate. Rapid decompression causes undercooling due to rapid water exsolution and subsequent increase in the liquidus temperature. Thus, a larger proportion of quickly quenched MI (G and GM) in VSN0 may be indicative of a faster ascent rate.

Overall, we suggest that MI proportions varies between VSN0 and VSN2 endmembers, and VSN1 lies somewhere between. VSN2 clearly has a larger abundance of T and MC type inclusions where VSN0 has a larger proportion of G and GM inclusions. We propose this variation of melt inclusion proportions reflects variability in time spent in the upper crust (varying proportion of MC inclusions) and possibly ascent rate (variation in GM inclusions). We note that as our samples are from bulk material, it is possible that the crystal population includes ante-crysts and xenocrysts; however, our large population of MI allows for us to overcome this issue. This hypothesis corroborates with findings from previous work on the geochemical populations of the VSN clinopyroxenes from Jorgenson et al. (2024). They find three crystal populations for VSN px: (1) px found predominately with a low crystallization temperature (as low as 839°C), patchy zoned without a rim, and relatively lower SiO₂, MgO, CaO, and higher in FeO, Na₂O, and MnO with respect to other VSN ignimbrite crystals; (2) high |$T$| (up to 1250 |${}^{\circ }$|C), not zoned, and with a high Mg# (⁠|$>$| 0.8 and up to 0.94) and Cr₂O₃ (⁠|$>$|0.25 and up to 0.94 wt %) which are predominately found in VSN0 and VSN1; and (3) variable zoning (patchy with and without zoned rims, sector zoned, and not zoned), with most temperature estimates 1000°C and lower Mg# and Cr₂O₃ than the second type, found in VSN2 and VSN1. They suggest that a deep seated mafic pulse of magma, characterized by high Mg#, Cr₂O₃, and temperature, blows through the upper crustal reservoir (and entrains some of the low |$T$| crystal) to erupt VSN0 and destabilizes the reservoir enough to erupt the main ignimbrite unit. Following MI textures and zoning patterns, we suggest that the VSN0 clinopyroxenes are mostly from the high |$T$| group of clinopyroxenes (1) as they are not zoned. The VSN2 clinopyroxenes are more likely to be from the variably zoned group (3) as the zoning patterns match this group.

Evidence of exsolution at depth

Volatiles play an important role in eruption processes especially in terms of magma bulk density and buoyancy, as well as influencing phase equilibria and crystallization (e.g. Edmonds & Woods, 2018; Anderson, 1995; Ghiorso & Gualda, 2015). Volatiles exsolve via first and second boiling, where broadly first boiling exsolution is due to a drop in the magma pressure and second boiling exsolution is due to crystallization (Edmonds & Woods, 2018; Townsend et al., 2019). Volatile exsolution lowers magma bulk density which in turn increases the volume and bulk compressibility and pressurizes the magma chamber (Townsend et al., 2019; Mastin et al., 2008). Degruyter et al. (2016) shows via thermo-mechanical modelling that with a higher bulk compressibility from exsolved volatiles the magma favours accumulation over eruption. In order to assess the state of magma exsolution of VSN magma we turn to variable proportions of VB volume fraction in G1 and GM inclusions.

As previously discussed, we utilized the (Moore et al., 2015) calculator for DTC combined with an approximate volume correction of 2 vol % to account for PEC. This resulted in an overall limit of DTC + PEC to be 6.7 and 7.3 vol % for px and lct (Fig. 7; Moore et al., 2015; Hanyu et al., 2020). Below this threshold we propose inclusions are trapped from a pheno-bubble free melt (Fig. 8.1a–c), indicating magma was either volatile free or had dissolved volatiles during trapping, as is the case for many inclusions. Of GM inclusions 64.3% px and 86.4% lct are below the PEC + DTC limit, where 35.1% and 46.8% of the px and lct G1 inclusions are below the limit. The inclusions above the limit (62.6% GM + G1 px, 35.1% GM + G1 lct) suggest trapping MI from a magma with exsolved volatiles (Fig. 8.2a–d). G1 px above the limit show a wide variation of volume fraction, with some MI above 70 vol % (Fig. 13). We suggest that these large-volume fraction MI are indicative of an excess fluid phase in the pre-VSN magma reservoir.

While all subunits have some MI with a volume fraction above the estimated limit, the proportion varies with stratigraphy, where 44.9% of the VSN0 MI are above the limit (27/50 px and 4/19 lct). VSN1 has 56.8% inclusions above the limit (176/281 px and 8/43 lct), and VSN2 has 63.7% inclusions above the limit (125/197 px and 19/29 lct; Fig. 13). Variability in MI volume fraction reflects variability in magma exsolution state at the time of crystallization. In this respect, VSN0 and VSN2 are endmembers, trapping from bubble-poor and bubble-rich melts, similar to the endmember variability seen in MI type proportions. For VSN0, we suggest that the reason for a lack of high volume fraction MI is that many MI are trapped before fluid exsolution can take place (Fig. 8.1a–c), indicating rapid magma ascent from a deep source. VSN2, on the other hand, has many high volume fraction MI suggestive of crystallization from a bubbly magma reservoir. This is further corroborated by the presence of fluid inclusions in several crystals which were observed in the 3D scans as voids. While the progressive increase in volatile fraction from bottom to top of the unit may seem counter intuitive, it corroborates with previous works of Jorgenson et al. (2024) who suggest that cpx erupted from VSN0, and in part VSN1, are likely from a deeply sourced mafic rapid pulse where time for magma exsolution, crystal growth and MI trapping would be minimal. Their PT estimates and textural data indicate that the VSN2 crystals were formed in a shallower magma reservoir, which we suggest is the main reservoir of bubble bearing magma. The increasing volume fraction with stratigraphy also agrees with findings of Vinkler et al. (2012) who finds an increase in the vesicularity and median bubble size in the juvenile material of the VSN eruption. Furthermore, the variation in volume fraction may be indicative of a slower ascent rate, which has been suggested by Vinkler et al. (2012), evidenced by syn-eruptive changes in vesicularity, bubble number density, and an increase in leucite microlites (Shea et al., 2009).

Evidence of a magma with excess volatiles in the shallow crust reservoir prior to the VSN eruption can also lead us to speculate on the unique set of conditions that lead to the Colli Albani ignimbrite eruptions. The mafic-alkaline nature of the magma (i.e. low viscosity) does not favour the accumulation of large volumes of magma, which are required for the Villa Senni caldera forming eruption (Giordano and the CARG Team, 2010). Thermo-mechanical modelling from Degruyter et al. (2016) and Townsend et al. (2019) suggest that with progressive magma recharge whether eruption or accumulation occurs is influenced by several factors, including magma injection rate, crustal viscous relaxation, initial reservoir volume, and presence of volatiles. Given a constant recharge rate and identical initial chamber volume, magma with exsolved volatiles is more likely to accumulate magma than a reservoir without exsolved volatiles due to magma compressibility (Degruyter et al., 2016). From our results, we can infer that the magma of CA is gas rich, and exsolved volatile in the reservoir may be what allows for accumulation of such a large quantity of low-viscosity magma without erupting. Furthermore, exsolved volatiles lend the magma to become more buoyant which can contribute to the eruption of large volumes of magma (Caricchi et al., 2014; Sigmundsson et al., 2020). This combined with the proposed fast ascent rate of Colli Albani magma (Jorgenson et al., 2024; Vinkler et al., 2012; Campagnola et al., 2016) allows a fresh perspective on the eruptibility of large volumes of low-viscosity magma.

CONCLUSION

Take-home messages for the MI community

Melt inclusions are tool to understand pre-eruptive magma. However, MI research is arguably one of the more difficult petrological endeavours as preparation is time consuming, difficult, and involves a many step process. In addition, once the data have been collected, there are a myriad of corrections and processes that must be accounted for (Aster et al., 2016; Maclennan, 2017; Kent, 2008; Gaetani et al., 2012; Moore et al., 2015; Audétat & Lowenstern, 2013; Wallace, 2005). While our findings are based on a very large number of MI, a similar study could be conducted with a smaller number of MI, as long as variability for host crystals is accounted for. In this work, we use the findings of our textural study to provide new insights for MI research:

1. We have developed a classification scheme based on shape, crystallinity, and vapour phase, which has allowed us to compare different inclusions from this study. These six categories are glassy bubble free (G), glassy with a single bubble (G1), glassy with multiple bubbles (GM), glassy with irregular bubbles (Gi), microcrystalline (MC), and tube shaped (T). While these classifications may not be relevant for all MI studies we hope the addition of nomenclature is helpful to the MI community.

2. MI should be considered with respect to their location in the crystal, as MI can be strongly associated with mineral zones (Fig. 10) or found in planar alignment indicating a secondary process capturing MI. Following works of Bodnar et al. (2006), Roedder (1979), Yang & Scott (2002), Rose-Koga et al. (2021), Esposito et al. (2018), and literature from the fluid inclusion community, MI should be considered as zonal or azonal and grouped into assemblages rather than a single parcel representative of the entire crystal. Additionally, clear planar alignments of MI that do not follow zoning patterns are best avoided as these are unlikely to be in equilibrium with the host (secondary inclusions).

3. Traditionally, MI studies avoid inclusions with multiple bubbles though they are clearly apparent in natural rocks (Frezzotti, 2001; Steele-MacInnis et al., 2017; Cannatelli et al., 2016; Wallace et al., 2003; Hanyu et al., 2020; Pintea, 2013; Drignon et al., 2021; Rose-Koga et al., 2021). While properly constraining the vapour phase in GM inclusions is more difficult, they are viable options for measuring volatiles and geochemistry. Additionally, GM inclusions are likely indicative of a faster cooling rate which may limit some other MI considerations (PEC, diffusive loss of H₂O, etc.)

4. Following our analysis of recalculated MI and VB volumes (Fig. 11) we strongly suggest that studies reconstructing total MI CO₂ using the VB obtain 3D data. Otherwise, we suggest users take a conservative approach and assume a spherical shape with the smaller of the two axes measured.

Colli Albani MI reveal a bubble bearing magma reservoir

Our comprehensive study of 1996 MI of the VSN eruption provides us with a novel view into the pre-eruptive state of the magmatic reservoir. By separating MI into distinct types we have been able to look at varying proportions, revealing stratigraphic variability. VSN0 shows a larger proportion rapidly quenched MI (G and GM type) where VSN2 shows a larger proportion of slowly cooled MI (MC inclusions). VSN1 has variable proportions of these MI types, which indicates a progressive slowing of quench rate from VSN0 to VSN2, which may be attributed to variation of deposit type or possibly magma deceleration. Additionally, VSN2 crystals show a strong zonation of MI-rich cores and MI-poor rims indicating a period of resorption and growth, which is markedly different than the VSN0 and VSN1 crystals.

VB volume fractions also reveal key information about the state of exsolution of the magma reservoir prior to MI trapping. Volatile-rich and glassy-type inclusions (G1 and GM) have a large-volume fraction, which indicates pheno-bubble trapping. The proportion of G1 and GM inclusions above the volume limit for DTC + PEC varies with stratigraphy, where VSN0 has the lowest and VSN2 has the highest (VSN0 44.9%, VSN1 56.8%, VSN2 63.7%). This is suggestive of varying exsolution states of the magma during MI trapping and crystallization. We suggest lower proportion of MI with a high volume fraction in VSN0 is indicative of magma from deep and may point to a rapid ascent. VSN2, on the other hand, is extremely volatile rich with volume fractions up to 78%, suggestive of a trapping pheno-bubble. We suggest this points to the magmatic reservoir of CA to be rich in exsolved melts prior to the eruption. As magma reservoirs with exsolved fluids are more likely to accumulate than erupt magma given a fixed magma input (Degruyter et al., 2016; Townsend et al., 2019), this may point to how such a large volume of low-viscosity magma can accumulate prior eruption. Given this information, in conjunction with works of Jorgenson et al. (2024), Vinkler et al. (2012), and Campagnola et al. (2016) who all suggest rapid ascent, we gain a better understanding that to create a large-volume mafic-alkaline eruption, it is necessary to not only have rapidly ascending magma, but to have a magma with enough exsolved fluid to be able to accumulate a buoyant magma that has the power to create a VEI 6 eruption.

Competing interests

No competing interest is declared.

Author contributions statement

L.C. and M.S. conceived the experiments. C.J., L.C., and M.S. wrote the beam time proposal. G.G. helped with sample collection. C.J. did the sample preparation. C.J., M.S., Gi.F., Ge.F., T.K., F.W., and G.W. all assisted with measurements and image reconstruction. C.J. did the sample preparation, image segmentation, and wrote the manuscript. L.C., G.W., G.G., and M.S. reviewed and edited the manuscript.

Supplementary Data

Supplementary data are available at Journal of Petrology online.

Acknowledgements

We acknowledge the work of Julian P. Moosmann for the experimental support (Hereon). We would like to thank Peter Westenberger from ThermoFisher who greatly aided in the image segmentation workflow in Avizo. We thank Jordan Lubbers who greatly aided in the image segmentation and edits of the manuscript. We thank Rosario Esposito and Victoria Smith for their feedback on this research. We also thank Terry Plank, Glen Geatani, and Ayla Pamucku for their insightful thoughts and discussions which greatly benefited the manuscript. We acknowledge DESY, a member of the Helmholtz Association for granting beamtime at beamlines P05 (PETRA III) operated by Helmholtz-Zentrum Hereon. This research was partly supported by the Maxwell computational resources operated at DESY. We would like to thank two anonymous reviewers and the associate editor Adam Kent for their feedback.

Funding

CJ and LC received funding from the Swiss National Science Foundation (Grant No. 200021_184632). GW acknowledges funding through an early postdoc mobility fellowship from the Swiss National Science Foundation.

Data Availability

Data underlying this article are available in the online supplementary materials.

References