A Crystal Mush Perspective Explains Magma Variability at La Fossa Volcano (Vulcano, Italy)

: The eruptive products of the last 1000 years at La Fossa volcano on the island of Vulcano (Italy) are characterized by abrupt changes of chemical composition that span from latite to rhyolite. The wide variety of textural features of these products has given rise to several petrological models dealing with the mingling/mixing processes involving mafic-intermediate and rhyolitic magmas. In this paper, we use published whole-rock data for the erupted products of La Fossa and combine them in geochemical and thermodynamic modelling in order to provide new constrains for the in-terpretations of the dynamics of the active magmatic system. The obtained results allow us to pic-ture a polybaric plumbing system characterized by multiple magma reservoirs and related crystal mushes, formed from time to time during the differentiation of shoshonitic magmas, to produce latites, trachytes and rhyolites. The residing crystal mushes are periodically perturbated by new, fresh magma injections that, on one hand, induce the partial melting of the mush and, on the other hand, favor the extraction of highly differentiated interstitial melts. The subsequent mixing and mingling of mush-derived melts ultimately determine the formation of magmas erupted at La Fossa, whose textural and chemical features are otherwise not explained by simple assimilation and fractional crystallization models. In such a system, the compositional variability of the erupted products reflects the complexity of the physical and chemical interactions among recharging magmas and the crystal mushes.


Introduction
Understanding the conditions of magma storage and their geochemical evolution in shallow crustal magmatic systems is one of the key points of igneous petrology, as it provides critical elements to link the architecture and mechanisms acting in a plumbing system to the assessment of future eruptive scenarios [1,2]. It is now widely recognized that magmatic reservoirs can remain for a long time in a low-melt state, with the development of crystal mushes [3][4][5][6]. Under these low-melt conditions, differentiated melts formed in situ during the solidification of the crystal mush may be periodically remobilized by hotter mafic recharging magmas [7]. In this scenario, batches of crystal-poor silicic magmas can be extracted from crystal mushes of intermediate composition through: (i) compaction-induced segregation of the interstitial melt [8,9], (ii) gravitational instability of the crystalline network [10][11][12], (iii) filter-press induced by volatile exsolution [13]. For all these processes, thermomechanical models indicate an optimal window for melt extraction in the range of crystal content of 50-70 vol.% in the mush [8,14,15]. The efficiency of crystal-melt separation has been documented for large volcanic systems, feeding voluminous ignimbrites [16][17][18]. Nonetheless, in the case of relatively small arc volcanic systems, crystal-melt separation mechanisms and their geochemical footprint on crystal-poor silicic magmas remain poorly documented.
The recent (i.e., in the last 1000 years) volcanic sequence of La Fossa volcano on the island of Vulcano, Italy represents a study case suitable to investigate the relationships between a mush-dominated crustal plumbing system and the geochemical variability of the erupted magmas. Volcanic products are in fact characterized by remarkable chemical heterogeneities, either in the form of mingled crystal-rich, latitic-trachytic magmatic enclaves in rhyolitic lavas and pyroclasts, or as mingled and partially mixed lati-trachytic and rhyolitic domains in pyroclasts and lavas [19][20][21][22]. Notably, these features are observed both at the scale of the eruptive sequence and within single eruptive units [19,22] and they have been attributed to mixing processes between intermediate mafic and rhyolitic melts, according to different petrogenetic models. [23] suggested a mingling process between an uprising rhyolitic magma and a shallower less evolved plug, followed by dispersion of latitic enclaves in the rhyolite to explain the textural features of the Pietre Cotte rhyolitic lava flow (La Fossa, AD 1739). [24,25] proposed instead that the crystal-rich enclaves occurring in the same eruption can be reconciled with the injection of the latitic magma into the more evolved rhyolitic magma. [26], and more recently [21,27], suggested that crystal-rich enclaves occurring in the Pietre Cotte lava are the result of remobilization of crystal-rich cumulates, derived from a crystal mush. Following this latter interpretation, the association between crystal-poor rhyolites and (heterogeneous) crystal-rich enclaves at La Fossa can provide the ultimate evidence of the twofold role of crystal mushes in the origin of differentiated silicic magmas and their eventual segregation to feed volcanic eruptions in small arc systems. In addition, the occurrence of K-Ba-rich trachytic magma (i.e., Pal D eruption) [21,28,29] suggests that a complex interplay between diverse mush-derived melts must be considered.
In this paper, we took in exam the magmatic products erupted during the recent (i.e., in the last 1000 years) activity of La Fossa volcano with the aim to determine the control of crystal mushes in the petrochemical evolution of intermediate and evolved magmas. At this purpose, we combine geochemical and thermodynamic modelling, using literature data of major and trace elements compositions of whole-rock samples encompassing the entire compositional range of the recent erupted magmas. Our results suggest that the more evolved trachytic and rhyolitic magma compositions cannot be fully explained through fractional crystallization plus crustal assimilation processes alone, but rather include the interaction among melts derived from distinct crystal mush regions. In such a system, the crystal mush behaves like a filter, where the interactions between residual (intercrystalline) melts with those derived from partial melting of distinct mush domains (induced by hot mafic recharging magmas) may hide the actual composition of magmas triggering and/or feeding small volume eruptions.

Geological Background and Stratigraphy of the Last 1000 Years
La Fossa is the active volcanic center of the island of Vulcano in the Aeolian archipelago, a continental volcanic arc located in southern Tyrrhenian Sea (Figure 1a) [20]. The volcanic activity of the island started at ~130 ka with high-K calcalkaline (mostly pre-50 ka) and shoshonitic (post-50 ka) magmas [30]. La Fossa volcano developed during the last ~6 ka in the middle of La Fossa Caldera, a volcano-tectonic depression located in the northern sector of Vulcano [31,32]. Together with Vulcanello, a volcanic center developed in the northern border of La Fossa Caldera, the La Fossa volcano represents the most recent volcanic structure on the island (Figure 1a) [28,33]. La Fossa erupted both with explosive (subplinian, long lasting vulcanian eruptions and ash emissions) and effusive styles, resulting in the emplacement of pyroclastic fallout deposits (ash layers, blocks and breadcrust bombs), subordinated pyroclastic density current deposits, lava flows and breccia deposits from phreatic/phreatomagmatic explosions ( [34] and ref. therein). Since the last eruption, that occurred in AD 1888-1890, the volcano is in a state of rest, characterized by an intense fumarolic activity and moderate seismicity ( [34] and ref. therein). [28] summarized the eruptive history of the last 1000 years in two main eruptive clusters: the Palizzi-Commenda eruptive cluster (PCEC), developed during the 12-13th centuries, and the Gran Cratere eruptive cluster (GCEC, 15th century-1890 AD) (Figure 1b). The stratigraphic sequence displays a large variety of eruptive products and a wide spectrum of magma compositions, ranging from latite to trachyte and rhyolite (Figures 1c and 2). The PCEC includes the products belonging to the Palizzi eruptive unit (PEU) and of the Commenda eruptive unit (CEU). The GCEC includes the products of the Lower and Upper Pietre Cotte eruptive units and the products of the last eruption (AD 1888-1890) ( Figure  1b).

Chemical and Petrographic Features of the Eruptive Products
The products emplaced at La Fossa and Vulcanello in the last 1000 years belong to the shoshonitic series and range from shoshonite to alkali rhyolite in composition ( Figure  2 and Supplementary Tables S1 and S2). The alkali content of the eruptive products in the trachytic field varies from ~9 to ~12 wt.% (whole rock), resulting in two distinct evolutionary trends: one directed towards the rhyolitic field with constant alkali content (9-10 wt.% Na2O + K2O) (Figure 2), and the other showing an enrichment in K2O (up to 7.5 wt.% K2O, with total alkali content being 11-12.5 wt.%) with nearly constant SiO2 ( Figure 2). The high-K2O products correspond to the trachyte of Pal D (cf. K-rich trachyte in the figures), a pyroclastic fallout deposit, belonging to the PEU, which products show the highest K2O content of the entire volcanic succession of the island [21]. The trend showing constant alkali, highlighted with a line in Figure 2, is represented by the pyroclastic products of Upper Pietre Cotte and of the AD 1888-1890 eruption (cf. Trachyte Hybrids in the figures) [19,22,42,43]. With regards to trace elements compositions, the rhyolitic products show markedly lower contents of Sr and Ba and higher Rb compared to latitic and trachytic magmas ( Figure 3). Interestingly, the Ba content of the K-rich trachyte is higher with respect to the other trachytic, latitic and shoshonitic products, up to ~2200 ppm ( Figure 3). Compatible Cr decreases as the degree of evolution increases, while La increases although with more scattering. The Eu anomaly, expressed by the Eu/Eu* ratio (with Eu/Eu* = EuN/(SmNxGdN) 1/2 ), becomes markedly negative with the increasing degree of evolution from latite to rhyolite, with the exception of the K-rich trachyte that shows higher values around 0.8 ( Figure 3). Literature data from [19][20][21][22]27,[35][36][37][38][39][40][41]. The straight lines in diagrams visually indicate the mixing trend between latitic and rhyolitic magmas, as suggested by [19,22]. The products showing evident textural mixing features and constant alkali (Trachyte Hybrids) plot along this line. See the Supplementary Table S1 for the complete data set. Major elements are recalculated to 100 wt.% on anhydrous basis.
In general, all the rhyolitic products of La Fossa exhibit evidence of mixing and mingling with less differentiated magmas, in the form of crystal-rich (with porphyritic index >35%) magmatic enclaves with quenched groundmass, phenocrysts in chemical disequilibrium with the host rock demonstrated by resorbed and reaction rims, bimodal mineral compositions and lati-trachytic variably porphyritic glassy domains (i.e., banding) in lavas and pyroclastic products (Figure 4) [20][21][22][23]27,43]. When the textural, petrographic and geochemical features of these products are considered (Figures 2, 3 and 4), it is evident that their compositions are determined by mixing processes between latitic and rhyolitic end members as suggested by [19,22,42,43]. Interestingly, if the products originated by mixing are not considered, there is a marked compositional gap in the bulk rock composition between the erupted latites and rhyolites.

The Plumbing System of La Fossa: Crystallization Conditions and Magma Differentiation
The temperature and pressure distribution inferred by several petrological studies is consistent with the polybaric nature of the plumbing system. [22] indicates crystallization temperatures of 1080 ± 10 and 1050 ± 10 °C for the latite and trachyte, respectively, involved in the AD 1888-1890 eruption. Pre-eruptive temperature of latitic and rhyolitic magmas involved in the mixing process leading to the extrusion of the 1739 AD Pietre Cotte lava flow are about 1010 °C for the latite and 950 °C for the rhyolite [25]. Temperature and pressure estimates obtained with mineral-melt thermometry by [21], suggest crystallization temperature of 1004 ± 14 °C and pressure of 160 ± 54 MPa for the K-rich trachyte of Pal D, 1007 ± 9 °C and 208 ±30 MPa for other trachytes (Palizzi lava), 1027 ± 5 °C and 307 ± 47 MPa for latites (Pietre Cotte enclaves) and 955 ± 8 °C for the rhyolitic magma of Pal B. The pressure of the uppermost rhyolitic reservoir of 30-60 MPa has been suggested by [22,44] based on the study of equilibration pressure of fluid inclusions found in metamorphic xenoliths. Based on plagioclase-liquid hygrometer, [21,45] estimated H2O contents ranging from ~2 to ~3.5 wt.% for the latitic and trachytic magmas respectively, and of 1.7-2.2 wt.% for the Pal B rhyolite. These values are in accordance with H2O contents in melt inclusions hosted in the K-rich magma (1.0-2.5 wt.% H2O) measured by [29], and higher with respect to the values between 0.3-0.6 wt.% and between 0.8-1.9 wt.% in the melt inclusions of shoshonitic and latitic compositions, respectively [46].
The great variability of petrographic and chemical characteristics of the products of La Fossa volcano, the widespread evidence of mixing processes and the relatively small volumes of the erupted magmas, overall confirm the vertical distribution of the plumbing system (as inferred from thermobarometry) and the mutual interactions among physically and chemically distinct magma bodies [20][21][22][23]25,27,37,47]. In terms of magma differentiation, [19,48] explain the geochemical and isotopic characteristics of the products emplaced at La Fossa and Vulcanello with a two-stage model. In the first stage, the shoshonitic basaltic magmas differentiate up to shoshonitic and latitic compositions with RAFC process (assimilation plus fractional crystallization in a continuously refilled and tapped reservoir) in the lower crust (~21 km) [44]. The second stage is dominated by AFC and mixing between evolved and intermediate magmas in several small reservoirs in the upper crust (up to 1-2 km) [19,48]. This view implies that all the magma bodies feeding La Fossa and Vulcanello eruptions reasonably originated by a similar parental magma and then evolved independently in various batches at several depths ( [20] and ref. therein). The least evolved magma emplaced in the last eruptive period is represented by the shoshonitic products of Vulcanello, although this is not considered the primary magma of the system based on its isotopic composition ( 87 Sr/ 86 Sr ~0.7046), SiO2 content (52-54 wt.%), and relatively low Mg# (45)(46)(47)(48)(49)(50) [19,48]. In fact, [35] suggest that the Vulcanello shoshonites differentiated from the parental shoshonitic basalts, represented by the melt inclusions hosted in olivine (Fo90) in the La Fossa products, as also reported in [46]. The transition from silica-undersaturated shoshonites to silica-saturated latites of La Fossa and Vulcanello can be explained by considering low rates of crustal assimilation (~2%) [35]. The assimilation of crustal material is also suggested by the isotopic signatures of rhyolitic products. Indeed, the mafic to intermediate rocks of La Fossa and Vulcanello have 87 Sr/ 86 Sr ~0.7046-0.7048, while the rhyolites reach values of 87 Sr/ 86 Sr ~0.7052-0.7056 [19,48]. The isotopic variations are thus correlated with the degree of evolution of the rocks with the most contaminated composition occurring in rhyolites, which are consequently interpreted as the result of AFC processes starting from a shoshonitic magma. The oscillatory variations of the isotopic composition with time testify to frequent arrivals of mafic magmas, AFC processes and syn-eruptive mixing [19]. The crustal contaminant is likely the upper crust material ascribable to the metamorphic rocks of the Calabrian Arc [19,49]. Based on different lines of evidence, the deepest reservoir (basaltic to shoshonitic) is located between ~21 and 18 km [44], a shoshonitic to latitic storage level has been recognized between ~17 and 12 km [20,37]. The shallow crustal reservoirs hosting latitic to trachytic magmas are located at shallower depths comprised between ~8 and ~3 km [21,37,44,50,51]. The uppermost reservoir hosting rhyolites is located at a depth of about 1-2 km [22,44,51].
In addition to rhyolite-MELTS, we used the energy-constrained thermodynamic model magma chamber simulator (MCS) [54,55]. The MCS code simulates the evolution of magmas that thermally and reactively interacts with contaminant rocks (through assimilation and partial melting) during magma cooling and crystallization. In order to model the AFC process of La Fossa magmas with MCS, we used the same conditions used for rhyolite-MELTS. MCS runs started at Tliq and terminated once the thermal equilibrium between the cooling magma and assimilant was achieved (Tend in Figure 5). In the simulations performed at 200 MPa and 4 wt.% H2O, MCS indicates 24% crystallization and a mineral assemblage composed by 14% olivine, 4% Fe-Ti oxide, 51% clinopyroxene and 31% plagioclase. At 200 MPa and 2 wt.% H2O the crystallization is 30% and the mineral assemblage is composed by 12% olivine, 43% clinopyroxene, 3% Fe-Ti oxide and 42% plagioclase. At 200 MPa and 0.5 wt.% H2O, the crystallization is 43% and the mineral assemblage is composed by 31% clinopyroxene, 11% olivine, 55% plagioclase and 3% Fe-Ti oxide. At 500 MPa and 4 wt.% H2O, a modal assemblage consisting of 79% clinopyroxene, 3% Fe-Ti oxide, 2% orthopyroxene and 16% garnet is reproduced with 24% crystallization. At 500 MPa and 2 wt.% H2O, the crystallization is 29% and the predicted mineral assemblage consists of 68% clinopyroxene, 2% Fe-Ti oxide, 5% orthopyroxene, 24% plagioclase and 1% garnet. Finally, at 500 MPa and 0.5 wt.% H2O, 58% clinopyroxene, 40% plagioclase, 1% Fe-Ti oxide and 1% orthopyroxene are reproduced after 36% of crystallization. The composition of both magma and wallrock for the runs at 200 MPa and 500 MPa are reported in the Supplementary Tables S4 and S5, respectively. The results obtained using MCS are comparable with those obtained using rhyolite-MELTS, although these indicate slightly lower percentage of crystallization ( Figure 5) and the formation of orthopyroxene and garnet (runs performed at 500 MPa) that are never observed in La Fossa magmas.

Trace Element Geochemical Modelling
Several previous works, on the basis of isotopic signatures and trace elements compositions of erupted products, invoked open-system crystallization processes with contamination of continental crust to explain the formation of rhyolitic magmas at La Fossa, starting from the less differentiated shoshonites [19,22,23,27,48,49]. The very low contents of Ba and Sr, and the marked negative Eu anomaly of La Fossa rhyolites (Figure 3), testify to processes of fractional crystallization with a significant involvement of feldspars. We thus examined the evolution of La Fossa magmas through trace element modelling using the FC-AFC equations proposed by [56]. We used bulk partition coefficients (D) estimated using modal mineral assemblages of the relative proportions of crystallizing phases obtained through mass balance calculations by [21], that better replicate the observed trend of trachyte-rhyolite differentiation (see next section), and a series of mineral/melt partition coefficients selected among those available in the literature (Supplementary Table S6) specific for latitic, trachytic and rhyolitic magmas, including those estimated for La Fossa products by [23,39]. Thus, two steps of magmatic differentiation (shoshonites-trachytes and trachytes-rhyolites) have been modelled using a minimum (low D) and a maximum (high D) bulk D, via pure FC and AFC with r between 0.2 and 0.5 (r, assimilated mass/crystallized mass), consistent with previous AFC models reported in the literature [19,23,48,49], in order to cover a series of initial trace elements concentration ( Figure 6). Additionally, we performed energy-constrained AFC trace elements modelling with MCS [57] using the phase equilibria reproduced at 200 and 500 MPa with 0.5 wt.% and 4 wt.% H2O and the same mineral/melt partition coefficients (Supplementary Table S6) already used in the FC-AFC equations of [56], thus obtaining a low-and high-bulk D. The results are reported in the Supplementary Table S7 and in Figure 6.
Starting from the same shoshonitic composition used for rhyolite-MELTS and MCS simulations, the models show that the trachytic products (including the K-rich trachyte) are not entirely reproduced by AFC processes. Indeed, after ~40% of crystal fractionation, the modelled residual liquids plot in a region where none of the intermediate-evolved samples of Vulcano occur, whereas the K-rich trachytic magma composition is completely out of this trend. The AFC process for the genesis of the rhyolitic magma has been modeled starting from the ending points of the shoshonite-trachyte fractionation path, thereby assuming a trachytic magma composition as the parental intermediate liquid from which the rhyolitic magma originated ( Figure 6). By these points, in accordance to previous works [19][20][21][22][23]27], rhyolites are possibly explained through an additional 40-50 % AFC with r between 0.2 and 0.5. The high fraction of sanidine involved in the fractionated assemblage [21] drives the liquid evolution toward depletion in Sr and Ba ( Figure 6).  [56] and the MCS code (in middle panels). The used bulk D have been estimated using modal mineral assemblages of the relative proportions of crystallizing phases obtained through the mass balance calculations [21] and selected mineral/melt partition coefficients and by the phase equilibria produced by MCS (see the text for further explanations and Supplementary Tables S6 and S7). The numbers at the ends of the differentiation paths refer to the crystal fractions. Literature data sources as in Figure 2. Only MCS runs performed at 0.5 wt.% H2O are showed for the sake of simplicity (see Supplementary Table S7 for the runs performed at 4 wt.% H2O).

Divergence between Natural Trends and Geochemical Models
The liquid lines of descent modeled by rhyolite-MELTS and MCS do not fully reproduce the whole compositional trend of natural products ( Figure 5). In particular, it is worth noting that latitic and trachytic magma compositions (including the K-rich trachyte) with K2O > 4.5 wt.%, Al2O3 > 17 wt.% and rhyolitic magma compositions with SiO2 > 70 wt.% are not achieved using rhyolite-MELTS. Similarly, the melt differentiation trends predicted by MCS do not fully reproduce the evolved melt erupted at La Fossa.
The divergence between thermodynamic models and natural trends is possibly explained by the fact that the thermodynamic modelling underestimates the amount of crystallization, resulting in lower values than those required for the genesis of rhyolitic magmas. Indeed, by using mass balance calculations, [21] suggest that the differentiation of highly differentiated rhyolitic melts from a lati-trachytic magma composition would require additional ~26-39% of sanidine fractionation in the solid assemblage, which is never predicted by thermodynamic models. These findings are in accordance with the thermodynamic simulations carried out by [58] that fail to reproduce the trachytic and rhyolitic melt compositions erupted at Vulcano between 50 and 8 ka, after about 50% of fractionated solid. Fractionation of these melts, in turn, would require additional increase of sanidine crystallization of 59-84% [58].
In a similar fashion, the AFC trace element modelling, including MCS, accounting for assimilation of continental crust, does not completely reproduce the intermediate latitic-trachytic products nor the K-rich trachyte, even considering a higher plagioclase/sanidine ratio in the modal assemblage with respect to other trachytes. Indeed, as suggested from mass balance calculations and temperature gradient experiments by [21], the K-rich magma would require a higher plagioclase/sanidine ratio in the fractionated assemblage, thus leading to an increase in K2O, Ba and Rb compared to the other latitic and trachytic products [21] (Figure 6). Moreover, the Ba-enrichment showed by the K-rich trachyte is not reproduced by energy-constrained AFC processes modelled by MCS, although this latter model seems to better reproduce the major element composition of this magma compared to rhyolite-MELTS ( Figure 5).
Compared to mass balance calculations of [21] and thermodynamic modelling, the AFC model performed with the equations of [56], would require a slightly lower amount of crystal removal. This divergence could be due to the fact that, to model the evolution from trachyte to rhyolite, the actual partition coefficients may vary during differentiation or other mechanisms of liquid fractionation may occur. Furthermore, the crystal-poor nature of the rhyolitic magmas erupted at La Fossa in the last 1000 years remains an apparent paradox considered the result of this modeling. Bearing in mind these considerations and considering also the fact that the lati-trachytic magmas (particularly the K-rich trachyte) are not entirely explained by trace elements and thermodynamic models, we explored other mechanisms able to explain the formation of intermediate magmas and crystal-poor differentiated magmas at La Fossa.

Combined Mush Crystallization-Melt Extraction Modelling
The solidification of intermediate magmas and the consequent development of a crystal mush is at the base of the production of the rhyolitic melt. In turn, this may also account for the compositional gap observed at SiO2 contents between ~62-69 wt.% ( Figure  2), corresponding to an area in the binary compositional diagrams occupied only by trachytic products showing evident mixing and mingling features (Figures 2, 5 and 6). In this framework, the formation of crystal-poor rhyolites has been explained by segregation of the interstitial melt from latitic-trachytic crystal mushes below La Fossa volcano, likely favored by the second boiling of the mushy magma [21], whilst mixed/mingled products could represent the result of late interaction between the segregating rhyolitic melt, remobilized portions of the crystal mush and the mafic magma responsible for remobilization (see below).
We simulated the extraction of rhyolitic melt from crystal mushes of intermediate compositions with the help of the geochemical model proposed by [15] (Figure 7). Through this model, consisting of equations for the conservation of mass, it is possible to simulate the evolution of a trace element in a magma mush environment simultaneously undergoing crystallization and extraction of interstitial melt. The model input includes the initial trace elements concentration and constant bulk partition coefficients; melt segregation rate is determined by a Gaussian curve centered on intermediate mush crystallinities (50-70 crystallinity peak) [15]. We used the high bulk D already used in the AFC models. The starting points of melt extraction are represented by the ending points of the shoshonite-trachyte differentiation trend shown in Figure 6. An additional starting point of melt extraction, plotting in the middle of the gap region outlined by the shoshonitetrachyte differentiation step has been also tested (Figure 7). Results of the melt extraction model are reported in the Supplementary Table S8. Although these starting compositions do not usually occur in the erupted products of La Fossa volcano, the presence of felsic crystal mushes from which rhyolite may be extracted can be envisaged, following the occurrence of sanidine-rich enclaves (sanidine-rich mush rock, classified as alkali-syenite) composed by idiomorphic medium-grained sanidine with minor plagioclase, femic minerals and interstitial glass [22]. These types of magmatic enclaves have been found in the rhyolitic products of the AD 1888-1890 eruption and interpreted as fragments of the crystallizing magma chamber at the boundary with the wall-rock [22]. In trace elements variation diagrams, this sample plots in the region of compositional gap highlighted by magmatic differentiation with its composition being not so dissimilar from the starting points of melt extraction (Figures 6 and 7). The extraction of interstitial melt occurs when the system reaches a mid-crystallinity peak of 50-70 vol.%, in agreement with the crystal fraction predicted by the trachyte to rhyolite differentiation modelled through simple AFC modelling ( Figure 6) and also in accordance with the models of interstitial melt extraction, suggesting that the most favorable window for crystal-melt separation occurs at a crystallinity of ~50-70 vol.% [14]. The extracted liquids match with the rhyolite composition erupted at La Fossa if low amounts, 10% to 30% of melt are extracted from the initial composition (Figure 7 and Supplementary Material S8). The rhyolitic melts at La Fossa, that are enriched in Rb, Zr and depleted in Eu, Sr, and Ti, can be thus explained by extraction of interstitial melts from crystal mushes of intermediate composition largely dominated by feldspars [14]. Similar conclusions have been proposed by [58] for the older (50-8 ka) differentiated trachytic and rhyolitic products of Vulcano, which can be reproduced by in situ crystallization starting from intermediate magmas.
The volume of rhyolitic magma emplaced by the Pietre Cotte lava flow has been estimated to ~0.002 km 3 [25,59]. Assuming that this volume of extracted magma constituted 10% of the interstitial melt present in the mush, it can be estimated that a volume of rhyolitic melt of ~0.02 km 3 was present in the shallower region of the plumbing system. Hence, considering a crystal mush with 50-70 vol.% crystals, the volume of the mush reservoir could be estimated in the order of ~0.04-0.1 km 3 .  Table S8 for results of the model). It is worth noting that the cumulate compositions left over by melt extraction compare well with the latitic crystal-rich enclaves occurring in the La Fossa rhyolites with respect to the cumulates compositions calculated with the models of [56] (see text). Literature data sources as in Figure 2.

The Crystal Mush Factory: A Petrologic Model for Intermediate Magmas
The widespread occurrence of magma mixing/mingling features in the rhyolitic products erupted at La Fossa (Figure 4) testifies to a pervasive interaction of the rhyolitic melts extracted from the trachytic crystal mush with other, compositionally distinct, mushy regions of the plumbing system. This observation is corroborated by [26], who interpreted the magmatic enclaves hosted in the rhyolitic lava flow of Rocche Rosse in the similar context of the nearby island of Lipari, as remobilized and thermally modified portions of the crystal mush leftover by the extraction of rhyolitic melt following a mafic recharge. Remobilization and chemical modifications of crystal mush may be caused by the rising of new batches of hotter mafic magmas able to reduce the crystallinity of the mush [18,60,61]. In the case of La Fossa system, new magma inputs may have either interacted directly (and mingled) with extracted rhyolitic melts (Figure 4d,e) or, indirectly, induced remelting within the crystal-rich mushy region [21]. In both cases, partial melting of the crystal mush could reduce the overall crystallinity and promote remobilization of both interstitial melt and crystal-rich fragments of the mush. Thus, the previously extracted rhyolitic melt, the remobilized crystal rich enclaves and the newly arrived recharging magma erupt together [26].
In order to test this hypothesis, we modelled the partial melting of the cumulates left over by the extracted rhyolitic melt, using the equation for aggregated non-modal fractional melting [62] (Figure 8 and Supplementary Table S9). It is worth stressing that the composition of cumulate leftover by rhyolite extraction calculated with the model of [15] compares much better with the composition of latitic magmas, including the crystal-rich enclaves, with respect to the calculated cumulate composition by simple AFC modelling (Figures 7 and 8). The thermal modification of the mush could be favored by periodical ascent of relatively mafic magma from depth, which in the case of last 1000 years history of the magmatic system may be represented by shoshonitic magma that erupted directly only at Vulcanello [20,35,36,63]. Calculated mixing lines between Vulcanello shoshonites and melts originated by the partial melting of the crystal mush show that all the latitic and trachytic compositions are possibly originated by variable degree of mixing between mafic melts and liquids produced by the partial melting of the crystal mush and cumulate crystals (Figures 8 and 9). This kind of process is able to explain the textural characteristics of the magmatic enclaves hosted in rhyolites, such as the sanidine overgrowth of high-Sr plagioclase [23] (Figure 4). According to [26], in order to form such high-Sr plagioclase, the Sr content in the melt must be much higher than that measured in Lipari and Vulcano magmas, yet consistent with the partial melting of a feldspar-rich crystal mush.  [62], of cumulates leftover from rhyolitic melt extraction (model of Figure 7) and calculated mixing lines between mush-melts and recharging magmas represented by Vulcanello shoshonites (see Supplementary Table S9). The light blue cumulates represent the bulk compositions used for the melting simulations. F (melting curve) refer to the liquid fraction. Numbers along mixing line refer to proportion of mush melts. Mixing in variable proportion between recharging magmas, partially melted crystal mushes and the extracted rhyolites accounts for the compositional and textural features of lati-trachytic products that show evidence of mixing and for the K-Ba-rich trachyte. The ideal compositional field of these mixing processes is highlighted in green. Literature data sources as in Figure 2.
This process of combined partial melting of feldspar-rich crystal mush and consequent assimilation of resulting melt by shoshonitic magmas is also supported by thermodynamic modelling carried out with rhyolite-MELTS ( Figure 9). Thermodynamic simulations have been performed in the same pressure, H2O-melt content and temperature paths (as those reported in Figure 5) by adding 0.2 g of sanidine per iteration as contaminant to a shoshonitic magma simultaneously undergoing AFC (Supplementary Table S10). The assimilation of a sanidine-rich mush may account for the genesis of intermediate lati-trachytic products, including the K-rich trachyte. Indeed, the K-rich trachyte genesis is influenced by strong interactions between shoshonitic evolving magma and a large amount of mush derived-melts resulting from the thermal modification of the crystal mush. This process is able to explain the enrichment not only in K2O, but also in Al2O3, TiO2, Ba, Sr and Eu (e.g., [64][65][66]) that characterize this unique composition in the eruptive record of La Fossa (Figures 2, 3, 8 and 9).  Table S10). Numbers at the end of each MELTS simulation refer to the amount of crystallization. The ideal mixing field between mush melts, latites and extracted rhyolite (same as Figure 8) is highlighted in green. Literature data sources as in Figure 2.
Mixing processes in variable proportion between shoshonitic magmas and melts derived from both the partial melting of the crystal mush and extraction of highly evolved interstitial melt (rhyolites) explain also the compositional and textural features of trachytic products that show evidence of mixing (lying in the straight line in major and trace elements variation diagrams and connecting the rhyolites with the lati-trachytes in Figures 2 and 3 and green fields in Figures 8 and 9). In this frame, the presence of a vertically extended mush column below La Fossa volcano could explain why in the eruptive record of the last 1000 years, scarcely differentiated magma compositions have been erupted predominantly at Vulcanello. At La Fossa volcano, the feeding basaltic-shoshonitic melts that periodically contribute to the thermal rejuvenation and growth of the vertically developed crystal mush system are modified by the crystal mush itself, through the interaction with the melts derived from its melting and/or solidification. The relatively short timescales of crystal residence (18-120 months) at La Fossa obtained through Sr diffusion in plagioclase by [37] are consistent with this scenario. In fact, these short storage times are apparently in contrast with the time gap occurring during eruptive events of the last 1000 years (which are on the order of 10 1 -10 2 years) [28,34], but are reconcilable with a fast remobilization of the mushy magma following mafic recharges. At Vulcanello, instead, shoshonitic magmas directly fed the eruptions in response to a less-developed mush system, possibly with the help of the structural influence of the border fault system of the La Fossa Caldera [35,36,63].

Conclusions
In this paper we examined the complex association of the latitic, trachytic and rhyolitic magmas erupted at La Fossa volcano in the last 1000 years and their relationship with the vertically-extended crystal mush that characterizes the shallow plumbing system. Shoshonitic magmas erupted at Vulcanello, but not at La Fossa cone, testify to a common parental magma from which latitic, trachytic and rhyolitic magmas eventually originated after a complex interplay of magma crystallization and interaction with mush-derived melts. While the most differentiated rhyolitic magmas periodically originate within latitictrachytic crystal mushes and segregated into crystal-poor reservoirs, hybrid intermediate compositions can originate from the mixing of rhyolitic melts and melts derived from the partial melting of the crystal mush, periodically perturbated by the arrival of fresh shoshonitic magmas from the deeper portion of the system. The shoshonitic magma thus contributes to the thermal rejuvenation of the crystal mush and its eventual growth at depth, but not actively in the eruption dynamics during the intermediate-evolved eruptions of La Fossa volcano. In such a system, the textural and geochemical variability of the erupted products testifies to complex histories that involve mush remobilization, partial melting of mushy regions, interactions between recharging magmas, mush-derived melts and extracted rhyolitic magmas.