A Data Driven Approach to Investigate the Chemical Variability of Clinopyroxenes From the 2014–2015 Holuhraun–Bárdarbunga Eruption (Iceland)

The Holuhraun–Bárdarbunga (Iceland) eruption lasted approximately 6 months. Magma propagated laterally through a 40 km long dyke, while Bárdarbunga caldera was collapsing. This event was intensely monitored, providing an opportunity to investigate the relationships between eruption dynamics and erupted products. Whole rock and melt inclusion data do not show chemical variations of magma during the eruption. Nevertheless, zoning patterns in clinopyroxene suggest temporal variations of intensive parameters during crystallization. We investigated the chemical zoning of clinopyroxene using a data driven approach, on major and trace elements analyses from lava flow lobes emplaced during the eruption. We applied hierarchical clustering (HC) to identify compositional groups based on major and trace element chemistry. This analysis identifies five compositional groups, which can be associated with specific petrographic features. One cluster represents the chemistry of hourglass sectors, two constitute the oscillatory zoned mantle of the crystals, one cluster corresponds to a seldom present bright rim (in back scattered electron images) in the outer portions of the crystals, and a last one represents most of the outer rims. HC applied to trace elements also identifies five compositional clusters, which highlight progressively more evolved clinopyroxene compositions from the core to the rim of the crystals. Two of the clusters identified with trace elements corresponds to major element clusters. All together the data suggest that the chemical zoning in the inner portions of the clinopyroxene crystals was generated by crystallization in the magma reservoir and interaction between hot magma propagating through the dyke and unerupted magma cooling within the dyke. The fraction of zones produced by interaction with colder portion of the magma residing within the dyke dropped during the eruption, potentially signaling the thermal maturation of the dyke. Some of the analyses reveal that relatively close to the eruption time (i.e., the outer portions of the crystals) the dyke intercepted a lens of low temperature magma with a chemical composition that is distinguishable from the 2014 to 2015 Bárdarbunga eruption. Our approach can provide insights on the evolution of deep processes occurring during long-lasting eruptions by combining the analysis mineral chemistry of erupted products with multiparametric monitoring signals.

The Holuhraun-Bárdarbunga (Iceland) eruption lasted approximately 6 months. Magma propagated laterally through a 40 km long dyke, while Bárdarbunga caldera was collapsing. This event was intensely monitored, providing an opportunity to investigate the relationships between eruption dynamics and erupted products. Whole rock and melt inclusion data do not show chemical variations of magma during the eruption. Nevertheless, zoning patterns in clinopyroxene suggest temporal variations of intensive parameters during crystallization. We investigated the chemical zoning of clinopyroxene using a data driven approach, on major and trace elements analyses from lava flow lobes emplaced during the eruption. We applied hierarchical clustering (HC) to identify compositional groups based on major and trace element chemistry. This analysis identifies five compositional groups, which can be associated with specific petrographic features. One cluster represents the chemistry of hourglass sectors, two constitute the oscillatory zoned mantle of the crystals, one cluster corresponds to a seldom present bright rim (in back scattered electron images) in the outer portions of the crystals, and a last one represents most of the outer rims. HC applied to trace elements also identifies five compositional clusters, which highlight progressively more evolved clinopyroxene compositions from the core to the rim of the crystals. Two of the clusters identified with trace elements corresponds to major element clusters. All together the data suggest that the chemical zoning in the inner portions of the clinopyroxene crystals was generated by crystallization in the magma reservoir and interaction between hot magma propagating through the dyke and unerupted magma cooling within the dyke. The fraction of zones produced by interaction with colder portion of the magma residing within the dyke dropped during the eruption, potentially signaling the thermal maturation of the dyke. Some of the analyses reveal that relatively close to the eruption time (i.e., the outer portions of the crystals) the dyke intercepted a lens of low temperature

INTRODUCTION
The 2014-2015 Bárdarbunga-Holuhraun (Iceland) eruption lasted for approximately 6 months and the emission of lava occurred in the Holuhraun region, approximately 40 km NE of the Bárdarbunga volcano (Sigmundsson et al., 2014). Lava outflow was accompanied by progressive collapse of the Bárdarbunga caldera, with a striking correspondence between the rate of caldera volume increase and the volume of erupted magma. This suggests that a reservoir associated with the Bárdarbunga plumbing system was progressively drained during the eruption . Earthquake swarms tracked the propagation of a horizontally extended sub-vertical dyke through which the magma moved discontinuously first toward the SE and then to the NE from the Bárdarbunga region to the eruptive fracture (Sigmundsson et al., 2014;Gudmundsson et al., 2016;Woods et al., 2019). This eruption was intensely monitored, which makes it an excellent case study to attempt to relate monitoring signals and the chemistry of the erupted products. Petrologic investigations find that the source of magmas might be heterogeneous, but homogenization processes within the magma reservoir masked both source heterogeneity and any potential evolution of magma chemistry during the eruption (Halldórsson et al., 2018;Hartley et al., 2018).
Several previous studies suggest that magma and melt inclusion composition is broadly constant over the course of the eruption (Browning and Gudmundsson, 2015;Gauthier et al., 2016;Geiger et al., 2016;Ruch et al., 2016;Ilyinskaya et al., 2017;Bali et al., 2018;Halldórsson et al., 2018;Hartley et al., 2018;Woods et al., 2019). We focused our analysis on clinopyroxene (cpx), which are characterized by different patterns of sector and concentric oscillatory zoning. We collected samples of lava flow lobes emplaced during the eruption, except for November 2014, for which we did not manage to reach the lobe formed during this period (Figure 1). We present an approach that combines classic petrography and geochemistry with unsupervised learning [Hierarchical clustering (HC)] to identify how chemical zoning in minerals capture the complexity and evolution of the Bárdarbunga-Holuhraun volcanic eruption.

Electron Probe Micro Analyzer (EPMA)
Major element analyses were collected using a JEOL JXA 8200 superprobe at the University of Geneva using 15 keV acceleration voltage, 20 nA beam current and a beam diameter of about 1 micron. Each element was measured acquiring for 30 s on the corresponding wavelength position and 30 s on each of the two background positions located on both sides of the peak. We used forsterite as a standard for Mg and fayalite for Fe, pyrophanite for Mn and Ti, wollastonite for Si and Ca, andradite for Al, Chromium oxide for Cr, and albite for Na. The uncertainty for each element was estimated by measuring the standards as unknowns and accepting only calibrations for each elements for which the ratio between the net counts per second during calibration and measurement was within of ±1%. We additionally performed the peak search for each element every 20 analyses to minimize the impact of any drift during prolonged analytical sessions.

Laser Ablation-Inductively Coupled-Mass Spectrometer (LA-ICP-MS)
Trace element analyses and high-resolution trace element mapping were performed by laser ablation-inductively coupled plasma-mass spectrometry (LA-ICP-MS) at the Department of Physics and Geology (University of Perugia), utilizing a Teledyne Photon Machine G2 laser ablation system coupled to a Thermo iCAP-Q ICP-MS (Petrelli et al., 2016a,b).
Trace elements analyses were performed on single spots of 25 µm and rasters characterized by widths of 10 µm. Data were collected for Li, Be, P, S, Cl, Sc, Ti, V, Cr, Mn, Co, Ni, Cu, Gd, Rb, Sr, Y, Zr, Nb, Cs, Ba, Hf, Ta, Pb, Th, U, plus Rare Earth Elements (REEs), with a repetition rate and fluence of 8 Hz and 4 J/cm2, respectively. The NIST 610 reference material, Si concentrations determined by EPMA, and USGS BCR2G glass were used as the calibration standard, internal standard, and quality control, respectively. Data Reduction was carried out using the Iolite v.3 software package (Paton et al., 2011). Under the reported analytical conditions, precision and accuracy are typically better than 10% (Petrelli et al., 2016a,b).
Trace elements were imaged following the line rastering technique proposed by Ubide et al. (2015). A spot size of 10 µm, a scan speed of 2 µm/s, a fluence of ∼4 J/cm 2 and a repetition rate of 10 Hz were utilized, respectively. When creating crystal maps, the number of analytes was reduced to Sc, Ti, V, Cr, Mn, and Zr to increase the efficiency of the acquisition. Crystal maps and transects were produced using the CellSpace tool (Paul et al., 2012). Finally, compositional zonation relationships were further interpreted using "Images from selections" with Monocle tool (Petrus et al., 2017).

Hierarchical Cluster Analysis
The procedure was performed using the library "cluster" within the freeware software R (R Core Team, 2017). Major element   analyses of cpx with totals higher than 98 wt.% (n. 809; but using all 876 analyses did not change the result of the analyses) were first normalized for each element to values between 0 and 1 to give them the same weight and avoid false results. The normalized values were used to calculate the Euclidean distance (d ij ) between each analytical spot as: where x to z are all measured elements and i and j are the indexes of all measured points. To identify the clusters of data we use the Ward minimum variance criterion (Ward, 1963), in which pairs of clusters are first identified considering the smallest values of d ij and then progressively merged minimizing the within cluster variance after merging. This procedure produces a dendrogram with decreasing number of branches upward. The length of the branches represents the difference in d ij required for separated clusters to be merged into one cluster. As the number of clusters is unknown a priori we apply 24 different methods (R package Nbclust; Charrad et al., 2014) that evaluate the potential number of clusters in our dataset and select the number that receive the highest number of votes (five clusters in our case).

RESULTS
The crystals are oscillatory zoned, in most cases show hourglass sector zoning and are generally characterized by an outer rim, which appears relatively bright in Back Scattered Electron images (BSE; Figure 2). Most of the crystals show multiple (2-4) regions of resorption, commonly in either the inner portion of the crystals or immediately before the outer rim (Figure 2a). In a few cases, the outer resorbed rim is also very bright (Figure 2b). Oscillatory zoning is never conspicuous inside the inner resorption zone but, when present, it is continuous through the hourglass sector, which is, in turn, truncated by the outer resorption zone and the outer rim of the crystals. We also observed cpx that partially includes olivine that shows normal zoning and a diffused pattern toward the outer portion of the crystal (Figure 2a). Plagioclase is also included in cpx but is not generally closer to the core than the innermost resorption zone (Figure 2a). Plotting the analyses in Harker diagrams with gray-scale reflecting the month of eruption of the sample, show no chemical trends of the cpx with time (Figure 3), which is consistent with previous observations (Halldórsson et al., 2018). Additionally, the variation of chemistry from the core to the rim of single cpx can span the entire range of chemical variability observed throughout the eruption. A general inspection of these diagrams highlights that no two elements can be used to identify groups of cpx with distinct compositions (Figure 3). In these situations an unsupervised learning approach can be useful (Bergen et al., 2019) to explore more than two chemical dimensions, without any need for an a priori assumption on the mechanisms responsible for the presence of compositional clusters. However, we stress that while this approach always identifies compositional clusters, their validity must be cross-checked with petrography and geochemistry.
HC was applied to all cpx major element analyses (Supplementary Table 1), with a modal average of 5 clusters present in our dataset ( Figure 4A). Considering all elements together, the different clusters are discernible by their chemistry (Figure 5). There is no combination of two elements, however, that could be used to unequivocally distinguish all five clusters, justifying our use of the clustering approach. Considering the dendrogram (Figure 4b), Clusters 1 and 2 are chemically similar, while Cluster 3 and 5 exhibit the largest compositional difference with respect to the other clusters, especially in Al 2 O 3 , TiO 2 , MgO, CaO, and Na 2 O (Figure 5). Cluster 4 differs from 1 to 2 for its lower Al 2 O 3 , CaO, and to a minor extent Na 2 O content, with higher FeO and MnO (Figure 5).
While no particular temporal trends were observed in cpx chemistry (Figure 3), the compositional clusters occupy specific textural positions, which indicate that most pyroxenes grew following the same sequence of events. Cluster 3 corresponds to the composition of the {−111} hourglass sector of the cpx crystals, which is visible in the majority of crystals (Figures 6a,b,g). Analyses collected in the innermost portions of the cpx crystals belong either to Cluster 3 (dark in BSE and representing portions of the {−111} hourglass sector) or to Cluster 1. Cluster 2 is generally observed in the vicinity of a partly resorbed rim (Figures 6a,c-f). The oscillatory-zoned mantle of the crystals is constituted by an alternation of Cluster 1 and 2. Cluster 4 generally includes analyses of the outer rims of the cpx, although in January and February 2015 it also appears in more internal portions of the crystals (Figure 6g). To test if the appearance of Cluster 4 in the inner portions of crystals was associated with cluster mis-assignment we checked the composition of the analyses closest to the core assigned to Cluster 4 and find that these composition are indeed typical of Cluster 4 as they occupy the mean compositions of Cluster 4 for all elements. Cluster 5 represents only 2% of the analyses (Figure 4c) and corresponds exclusively to the bright and partially resorbed zone just before the outer rims of the cpx. The large chemical variance of this cluster (Figure 5), reflects a thickness for these zones that is at the limit or smaller than the EPMA spatial resolution. Therefore these analyses are of the lowest quality and at least some of the analyses represent mixtures of Cluster 5 composition and the chemistry of the neighboring zones.
The same approach we used for major elements was applied to analyses (no. 128) of selected trace elements present in sufficiently high concentration in cpx (Sc, Ti, V, Cr, Mn, Co, Ni). Including elements in low concentrations only adds noise to the data and complicates the clustering results (Supplementary Table 2). The number of clusters that received the highest amount of votes with the different tests are 2 and 5 (Figure 7k). Harker diagrams of the dataset suggest a better visual partitioning in five clusters with respect to two clusters (Figures 7f-j). The comparison between the textural positions of the clusters identified using major and trace elements (we identify the clusters identified by trace element adding "TE" after the cluster number) reveals some correspondence between the two. As Cluster 3, also Cluster 3TE corresponds to the compositions of the {−111} hourglass sector (Figures 7b,c). Cluster 5TE is identified in the outer portions of the crystals, which corresponds to the textural location of Cluster 4. Cluster 2TE, 4TE, and 1TE, constitute cores and mantles of the crystals and occupy progressively more external positions within the crystals (Figures 7a-e).

DISCUSSION
The application of HC to the major element chemical analyses of cpx of the 2014-2015 Holuhraun-Bárdarbunga eruption allowed the identification of five compositional clusters and their relative  position within the crystals. Cluster 1 and 2 constitute most of the inner portions and mantles of the crystals, respectively; Cluster 5 represents a sparsely present bright (in BSE) outer portion of the cpx; the outer rims are part of Cluster 4, and analyses belonging to Cluster 3 are from the {−111} hourglass sectors (Figure 6). In the following we discuss the processes that could be responsible for the chemical variability observed in the cpx crystals and discuss potential relationships between mineral chemistry and the dynamics of magma migration during the eruption. The cpx crystals we analyzed have core to rim distances between 250 and 500 µm, which, considering experimental growth rates under kinetically controlled growth at temperature typical for basaltic magmas of 5 × 10 −10 to 1.5 × 10 −9 m/s , imply that they grew over a period of 2-9 days. Using liquid-cpx thermobarometry (Putirka, 2008;Neave and Putirka, 2017) we calculated cpx crystallization pressures to estimate the structure of the plumbing system. We first selected only cpx crystals that were in equilibrium with the melt, following the procedure described in Neave et al. (2019). Considering the aphyric nature of the erupted magmas (Halldórsson et al., 2018) we use whole rock chemistry with H 2 O content of 0.4wt.% (based on Bali et al., 2018) as an equivalent for the melt in equilibrium with the cpx crystals. Analyses in Clusters 1, 2, 4 passed the equilibrium test while no analyses from Cluster 5 were found to be in equilibrium, and only two analyses from the {−111} hourglass sector (i.e., Cluster 3) passed the equilibrium screening (Figure 8).
Similarly to previous determinations, the pressure estimates vary between about 100 and 400 MPa with the highest density of estimates around 200 MPa (Neave et al., 2019; Figure 8). The calculated storage pressure progressively decreases from Cluster 1, to 2 and 4, suggesting concurrent magma rise and crystallization. However, as these variations are within the uncertainty of the geobarometer we do not discuss further their bearing on the pre-eruptive crystallization depth ( Figure 8G). Nevertheless, the distribution of pressures obtained with cpxliquid barometry indicate that over the 2-9 days before eruption, magma crystallized mainly at pressures around 200 MPa (i.e., 5-6 km depth), which is close to the depth range over which magma within the dyke has been interpreted to be transported laterally to the to the eruption site (Woods et al., 2019).
Prior to using cpx chemistry to retrieve the conditions of magma storage or the dynamics of magma migration to the surface, it is important to understand whether the chemical variations we have identified are due to variations in pressure (e.g., Nimis, 1995;Nimis and Ulmer, 1998), temperature (e.g., Lindsley, 1983), melt chemistry (e.g., Neave and Putirka, 2017), or growth kinetics Lanzafame et al., 2013;Mollo et al., 2013Mollo et al., , 2012Mollo et al., , 2011Ubide et al., 2019). Additionally, experiments show that cpx zoning resulting from kinetic effects can be produced even under nominally equilibrium conditions (Neave et al., 2019), and crystal can grow rapidly generating skeletal/dendritic structures, which are subsequently filled during growth at lower rate Welsch et al., 2016). All these effects should be taken in consideration before using cpx chemistry to identify storage conditions or processes occurring in the volcanic plumbing system.
The chemical composition of Cluster 3 is in agreement with previously measured cpx {−111} hourglass sectors (Ubide et al., 2019; Figure 6). Cluster 5 and 1 shows enrichment of Al, Ti and depletion of Na and lower Si, Ca, Mg (more extreme for Cluster 5) with respect to Cluster 2 (Figure 5). The chemistry of Cluster 1 and 5, could reflect either growth from magmas of different compositions, or as shown experimentally, be the consequence of more rapid crystal growth at relatively high undercooling with respect to Cluster 2 . To test these two hypotheses, we selected cpx crystals with the highest textural complexity and with chemical composition covering all clusters identified with major elements.
Trace element profiles consistently show that Cr decreases from core to rim in a stepwise manner while Zr, Sc, V, Ti, covary (Figures 9d,e; we report only Cr, Zr and Sc in the figure for clarity). The profiles show that Cluster 1 is generally characterized by higher concentrations of Zr and Sc and moderate to no decrease of Cr with respect to Cluster 2 (Figures 9d,e). Zr is incompatible in a bulk crystallizing assembly constituted by olivine, cpx and plagioclase, as is Sc, at least until the relative fraction of cpx in the crystallizing assembly does not reach about 50 wt.% and the bulk partition coefficient becomes larger than 1; partition coefficients from Earthref.org). As Cr is compatible in the crystallizing assemblage, with crystallization the concentrations of Zr and Sc in the residual melt phase (and therefore in the cpx) increases while Cr should decrease. Zr and Cr are inversely correlated in some portions of the profile, as it would be expected if the growth of the cpx zones occurred from magmas at different degrees of crystallinity, while they covary in other regions (Figures 9d,e), as it would be expected if the crystal growth occurred at increasing cooling rates (or undercooling; Mollo et al., 2013). No clear relationship is observable between the clusters identified by major element chemistry and the Cr-Zr variations (Figures 9d,e), nor between Cr-Zr and the trace element clusters (Figure 10). This implies that while variations of undercooling were occurring during crystal growth, they were not the only processes modulating the major element composition of the cpx. Therefore, Cluster 1 zones either (1) crystallized from a magma of similar initial composition but slightly more evolved and potentially cooler than Cluster 2, or (2) these two clusters were formed from magmas with distinct chemistry. Regarding the lack of relationships between Zr-Cr and the trace element clusters, this suggest that the ratio between these two elements is not the most significant factor controlling the trace element clusters. To test if any evidence exist of the interaction between magmas of distinct chemistry, we performed HC on the trace element profile dataset (Supplementary Tables 3-5) and include the LA-ICP-MS spots (Figure 10; clusters recognized on the base of trace element chemistry are identified by TE). Including all data, the result is similar to that obtained considering only the trace element spots data, however, the few analyses associated with Cluster 5TE are not identified as a separated cluster when including all data. Cluster 1 and 2 correspond to Cluster 4TE and Cluster 2TE, respectively and their trace element chemistry suggest that Cluster 1 cpx could crystallized at lower temperature (i.e., in a more fractionated magma) with respect to cpx from Cluster 2.
Thus, altogether major and trace element data suggest that Cluster 1 and 2 zones were formed by crystallization from the same magma at lower and higher temperature, respectively ( Figures 5, 9, 10). The higher temperature of zones belonging to Cluster 2 is also confirmed by their common growth over resorption surfaces (Figures 6a,c,d-f).
The trace element chemical evolution of cpx crystals from core to rim, can be plausibly explained with their crystallization from the same parental magma but at different temperatures and crystallinities (Figure 10). Few points located in the outer rim (Cluster 5TE), show higher content of incompatible elements such as Zr, Ti, and V and lower content of Sc FIGURE 11 | Evolution of the fraction of Cluster 1 relative to Cluster 2. Note that we did not have samples from November 2014. (Figure 10). These few analyses suggest that some of the outer portions of the cpx crystals crystallize from a distinct melt. This magma does not need to be derived by a different parental magma with respect to the Bárdarbunga-Holuhraun 2014-2015 magma, but could be representing a highly crystalline lens of cold magma left over from previous eruptions. This scenario is in agreement with the major elements chemistry (Figure 5) of the rare bright outer zone (Cluster 5) of some cpx crystals (Figures 7c, 8a), which could indicate crystal growth under rapid cooling conditions .
In summary, major and trace elements suggest that Cluster 1 and 2 rims grew from magmas with the same or similar composition but at lower (more chemically evolved melt) and higher temperature (lower concentration of incompatible elements; Figure 9), respectively. As Cluster 1 and 2 constitute the oscillatory zoned inner portions of the crystals, such temperature variations must have occurred repeatedly during the 2-9 days of cpx growth. The entrainment of cpx crystals in a higher temperature magma propagating through the dyke would lead to partial resorption of crystals followed by growth at lower degrees of undercooling and the formation of Cluster 2 growth zones. These distinct variations in temperature could be due to a number of factors, including: (1) thermal variability within the dyke due to non-uniform thermal exchange of magma with the wall rocks (e.g., complex geometry), which would result in high degrees of undercooling and crystallization of cpx marked by higher contents of Al, Ti and Na as observed for Cluster 1 and 5 compositions (Figure 5; Mollo et al., 2013); or (2) the temperature of the magma feeding the eruption was heterogeneous and mingling and mixing occurred within the dyke during transport resulting in the formation of Cluster 1 and 2 cpx zones. The fraction of Cluster 1 progressively decreases with respect to Cluster 2 from October 2014 to the end of the eruption (Figure 11). In the first scenario, this decrease would reflect the thermal maturation of the feeder dyke leading to a decreases of cpx crystallization at high degrees of undercooling. Within the second hypothesis, such decrease of the fraction of Cluster 1 cpx crystals would result from a decrease of the fraction of lower temperature magma feeding the dyke during the eruption (Figure 12).
Major elements suggest that Cluster 5 was formed at the highest degrees of undercooling with respect to the other compositional clusters (Figure 5; Mollo et al., 2013). Because of the limited spatial resolution, the trace element composition of Cluster 5 is difficult to determine with high precision, however, trace element profiles reveal that the presence of the bright rim in the cpx crystals (Figure 9a) likely grew from a magma with distinct chemistry (Figure 10). On the basis of these data we tentatively suggest that Cluster 5 compositions were formed when the magma came into contact with chemically different pocket/s of relatively cold magma while in the proximity of the eruption site (Figure 12), an hypothesis proposed also by Hartley et al. (2018). Finally, the rims (Cluster 4) were formed at shallow depths, but possibly still within the dyke, as this compositional cluster can also be observed in the inner portions of the crystals from January 2015 (Figure 5g).

CONCLUSION
The Holuhraun-Bárdarbunga eruption lasted for approximately 6 months during which the bulk magma chemistry remained substantially the same (Halldórsson et al., 2018). However, mineral chemistry can gives more resolution on the chemical processes occurring within magmatic systems (Wallace and Bergantz, 2002;Perugini et al., 2005;Kahl et al., 2011;Cheng et al., 2017;Forni et al., 2018;Probst et al., 2018). Using HC we have identified five groups of cpx compositions in major and trace element space that constitute different zones of the crystals. This, together with the petrographic context, allow us to interpret the evolution of the crystallization conditions during the eruption and to identify the different processes responsible for the chemical variability measured in the cpx crystals. Major and trace element analyses suggest that the chemical variability in the inner portions of the cpx crystals is mainly modulated by growth at different temperatures and resulting kinetic effects (Figures 5,  9, 10; Mollo et al., 2013;Ubide et al., 2019). The chemistry of the outer portions of the crystals reflect copious crystallization at shallow depths (Figures 5, 6, 8, 9). The chemistry of Zr-rich rims present in some of the crystals suggest the potential interaction of magma with pockets of cold and crystal-rich magma during the propagation of the magma to the eruption site (Figure 12). Compositional data, after removal of kinetic effects (Neave et al., 2019) suggest pre-eruptive pressure and temperature distributions in agreement with previous estimates (Figure 8; Halldórsson et al., 2018).
Our analysis suggests that the oscillatory zoning in the cpx mantles could be the result of either the interaction between magma propagating to the eruption site and pockets of magma that were cooling within the dyke, or by the mingling-mixing of magmas released from the reservoirs at different temperatures (Figure 12). The fraction of cpx rims generated at higher degrees of undercooling [Cluster 1/(Cluster1 + Cluster2)] progressively decline during the eruption (Figure 11), highlighting either thermal maturation of the dyke during the eruption, or the decrease of the fraction of the cooler magma entering the dyke. We propose that the approach presented here can be used to trace the evolution of long-lasting eruptions and can be useful to link monitoring parameters to magmatic processes occurring at depth.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
LC conceived the study, led the fieldwork, and wrote the first draft of the manuscript. GS and LP assisted with the sample collection and data interpretation. MP collected the LA-ICP-MS data and contributed to their interpretation. TS assisted with the unsupervised learning approach and helped with the data interpretation. EB helped to interpreting the geochemical data and shared her expertise on the Holuhraun-Bárdarbunga 2014-2015 eruption. All authors contributed to the final version of the manuscript.

ACKNOWLEDGMENTS
We are grateful to Jean-Marie Boccard for the preparation of excellent rock thin sections.