Magma recharge in persistently active basaltic–andesite systems and its geohazards implications: the case of Villarrica volcano, Chile

We report whole-rock chemistry, mineral chemistry, and volatile content from Villarrica volcano’s major recent paroxysms and background activity. Composition of the volcanic products are basalt to basaltic andesite with whole-rock SiO 2 content between 50 and 56 wt%, and a mineralogy dominated by olivine (Fo 71-80 ), clinopyroxene (Mg# ~ 50) and plagioclase (An 60–80 ). Volatile contents in melt inclusions are up to 1.5 wt% H 2 O, 500 ppm CO 2 , 1230 ppm sulfur and 580 ppm chlorine. Regardless of the type of activity, there are no substantial variations in whole-rock composition or the volatile content when the activity switches from background activity to a major paroxysm, strongly suggesting that this shift does not just depend on the arrival of new magma in the shallow magmatic system. Geothermobarometry constrains crystallization of the major mineral phases at various depths between 3 and 12.7 km, suggesting that degassing of a volatile-rich recharge magma occurs deeper than 12 km, producing efficient mixing throughout the whole system, and sustaining the lava lake activity in Villarrica’s summit crater. The occurrence of a permanent lava lake also suggests that the magma recharge must be close to continuous and therefore sudden changes between background and paroxysmal volcanic activity are likely controlled by relatively small changes in the rate of recharge and/or the volatile release rate in the magmatic system. This has important implications for the understanding of eruption triggers and the forecasting of volcanic eruptions.


Introduction
Magma recharge and mixing in long-lived volcanic systems are considered the key processes that can trigger volcanic eruptions (e.g., Sparks et al. 1977;Murphy et al. 1998;Eichelberger and Izbekov 2000;Richer et al. 2004;Morgavi et al. 2017), although not always do these processes drive a volcanic system into an eruption (e.g., Popocatépetl: Mangler et al. 2022;Krafla: Blake and Cortés 2018).
In intermediate and evolved magmatic systems, in which the compositional contrast between a resident magma and a primitive recharge is substantial, the arrival of new magma (e.g., by dyke injection) usually results in an eruption triggered by the remobilization of a resident crystal mush driven by heat and/or mass transfer (e.g., Richer et al. 2004;Bergantz et al. 2015).In these systems it is usually, although not always (e.g., Ruprecht et al. 2020), relatively straightforward to determine the end-member compositions involved in magma recharge and mixing, to estimate the conditions associated with magma recharge (e.g., relative volumes, pressure/depth and temperature), and to constrain the overall process (Sparks and Marshall 1986;Humphreys et al. 2013), although it has been shown that there is often a cryptic dispersal of mafic-magma-derived components in the more evolved host magma (e.g., Humphreys et al. 2013).In terms of numerical modelling of this process and its overall dynamics, it is generally accepted that mushy, crystal-rich resident magmas (e.g., Cashman et al. 2017;Magee et al. 2018;Lissenberg et al. 2019) can be remobilized, and efficiently mixed with a comparatively crystal-poor recharge magma (Bachmann and Bergantz 2006;Bergantz et al. 2015).This means that the ultimate eruptive behavior of such systems is directly controlled by the volume and frequency of magma recharge, the parameters of which can be measured or estimated by current monitoring techniques.
At the other end of the compositional spectrum lie the primitive magmatic systems (i.e., SiO 2 ~ 50-53 wt%), in which the differences between a resident magma and a recharge magma are subtle.In these systems, efficient mixing is promoted by similar compositions, thermal conditions and viscosities, and the mechanical stirring from volatile degassing, leading to the generation of a fully hybridized basaltic-andesite magma (e.g., Stromboli: Cortés et al. 2005;Arenal: Streck et al. 2002, 2005), while efficient heat transfer promotes the development and sustainability of a lava lake (e.g., Harris et al. 1999;Moussallam et al. 2016).Another key aspect on these systems is that the arrival of a new mafic recharge in the system does not necessarily culminate in an eruption (Blake and Cortés 2018).Therefore, what dictates the onset of an eruptive episode is the efficiency of the mixing in conjunction with any geometrical constraints imposed by the magma plumbing system (e.g., Cashman et al. 2017;Magee et al. 2018) and its interplay with regional tectonics (e.g., Cembrano and Lara 2009).As such, a general model of magma recharge and mixing would be complex, without a clear relationship between the arrival of a recharge and an eruptive event.
Here, we are interested on investigating such primitive magmatic systems and specifically, what controls their behavior when new magma arrives in the shallower volcano plumbing.We focus on Villarrica volcano, Chile (39.42° S, 71.93° W; Fig. 1), a well-known basaltic andesite volcano located in the Southern Volcanic Zone of the Andes (SVZ, Stern 2004), characterized by several relevant historical eruptions (e.g., 1963, 1971, 1984 and 2015) and steady-state activity manifested by a lava lake that has been active since at least 1985 (Calder et al. 2004;Witter et al. 2004;Moussallam et al. 2016).We compare geochemical and textural data between eruptive products of the major historical eruptions and the lavas emplaced during steady state or background activity and examine the link between its eruptive activity and the influence of the local and regional tectonic regime (Cembrano and Lara 2009;Schonwalder-Angel et al. 2018;Simmons et al. 2020a).

Basaltic-andesite volcanism in the SVZ
Andean magmatism is a direct consequence of slab dehydration and partial melting of the overlying mantle wedge (López-Escobar et al. 1977;Thorpe 1984).Depending on the crustal thickness, variable degrees of continental crustal contamination may occur (Davidson et al. 1990(Davidson et al. , 1991;;Wörner et al. 1994).In the SVZ, where the crust is relatively thin (~ 35 km;Cembrano and Lara 2009) and composed predominantly of post-Cretaceous intermediate volcanic rocks, a crustal signature is negligible or difficult to recognize (e.g., Thorpe 1984).Magma residence times at the nearby Puyehue-Cordón Caulle complex (Sigmarsson et al. 2002;Jicha et al. 2007) and Llaima volcano (Reubi et al. 2011) are relatively short (< 1000 yr) as calculated by U-Th-Ra disequilibria in mineral separates, limiting the timescales over which crustal contamination can occur.
More evolved magma compositions and longer residence times may be a function of the regional stress regime (Cembrano and Lara 2009;Schonwalder-Angel et al. 2018;Simmons et al. 2020a).At volcanoes in the compressive tectonic regime of the northern SVZ, prolonged magma storage has resulted in generally more evolved or differentiated (andesitic to dacitic) magmas.Volcanoes located in the central SVZ section (between 37° S and 41.5° S) are influenced by the interplay of dextral-transpressional tectonics in the volcanic arc, associated with the Liquiñe-Ofqui fault zone and compressional stress imposed by the back arc (Göllner et al. 2021).This interaction generates the geometry and the appropriate kinematics for the development of second-order structures such as tension cracks, shear fractures and volcanic fissures, facilitating relatively unimpeded magma ascent from the source, with relatively little time for magmatic differentiation to occur (Cembrano and Lara 2009).The specific tectonic setting hosts the distinctive high alumina basaltic to basaltic-andesite stratovolcanoes of the region, including 13 centers with historical activity (Petit-Breuilh 2006) and several others with post-glacial (< 12.5 ka) activity (Stern et al. 2007).

Villarrica volcano
With more than 60 eruptions VEI > 2 since 1553, Villarrica (Fig. 1) is by far the most active volcano in the SVZ of the Andes (Petit-Breuilh 2006;Global Volcanism Programme 2023).Radiogenic and oxygen isotopic compositions suggest that the magma source for its volcanic products is the subcontinental lithospheric mantle enriched by slab-derived fluids or melts, possibly with a minor component from the lower crust (Tormey et al. 1991;Sigmarsson et al. 2002).Probabilistic analyses of eruption frequency, developed by Dzierma and Wehrmann (2010) and Wehrmann and Dzierma (2011) The presence of an active lava lake (Moussallam et al. 2016) and persistent degassing (Palma et al. 2009(Palma et al. , 2011) ) are evidence of sustained or nearly continuous recharging of the shallow magma reservoir.Punctuating this behavior are regular paroxysmal eruptions (Barberi et al. 1993), the last of which occurred in February 2015 with little seismic warning (Johnson and Palma 2015;Aiuppa et al. 2017).The event culminated in violent Strombolian fountaining (Valentine and Gregg 2008) that lasted ~ 3 h, later returning to normal levels of background activity characterized by persistent degassing and reappearance of an active lava lake near the vent of the volcanic edifice.Historical explosive eruptions (VEI 2-3) occurred in 1922, 1948, 1963 and 1971, generating > 10 km high eruption columns, and extensive lava flows.During the 1964 eruption, lahars generated by melting of the summit glaciers, caused loss of lives and infrastructure in the town of Coñaripe (Moreno and Clavero 2006;Petit-Breuilh 2006), highlighting the importance of a better understanding of the factors that control the onset of an eruptive event.
Volcanism at Villarrica has been divided into three main stages (Fig. 2), separated by two major caldera-forming events (Moreno and Clavero 2006).Stage I (> 90-14 ka) involved the construction of a large stratocone, which terminated in the caldera-forming Licán Ignimbrite eruption that deposited ~ 10 km 3 (DRE) of basaltic andesite ignimbrite (Clavero 1996;Moreno and Clavero 2006;Lohmar et al. 2007Lohmar et al. , 2012)).Stage II (14-3.7 ka) involved the construction of a second stratocone on the north-western margin of the first caldera and ended in a smaller caldera-forming eruption associated with a ~ 3.3-5.0km 3 (DRE) basaltic andesite ignimbrite known as the Pucón Ignimbrite (Clavero 1996;Moreno and Clavero 2006;Silva-Parejas et al. 2010).Both major caldera collapses are inferred to have been triggered by the interaction of magma with external water (Lohmar et al. 2007(Lohmar et al. , 2012)).Stage III (< 3.7 ka) represents the construction of the currently active cone, again on the north-western margin of the two previous calderas (Fig. 2), characterized by basaltic andesite lavas and pyroclastic successions of violent Strombolian to sub-Plinian fallout and surge deposits (Moreno and Clavero 2006;Costantini et al. 2011;Pioli et al. 2015).The volcanic system has also developed two distinctive monogenetic volcanic fields during different stages of activity in the Holocene (Moreno and Clavero 2006), located to the east (Los Nevados) and south of the main volcanic edifice (indicated as "a" and "b" in Fig. 2).Moreno and Clavero (2006) and Petit-Breuilh (2006) summarized the historical activity at Villarrica since 1553, which has displayed a variety of eruptive styles (e.g., strombolian to violent strombolian with fountaining, lava flows and lahar development) in 1922, 1948, 1964 and 1971, involved significant explosive activity (VEI 2 to 3), and produced extensive lava flows from flank fissures.By contrast, the 1984-1985 eruption (VEI 2) generated only minor lava flows from summit-crater overflow.Current Villarrica products regardless of eruptive style span a narrow compositional range from basalt to basaltic andesite (Hickey-Vargas et al. 1989;Witter et al. 2004;Constantini et al. 2011).

Sample suite
This study focuses on 30 samples that are representative of Villarrica's different eruptive styles as described in the previous section (Supplementary Table 1).The 1999 and 2000 background activity samples were collected on, or just outside the crater rim during a period when explosive activity in the crater was elevated and the free surface of the lava lake within the crater was relatively high (Ortiz et al. 2003;Calder et al. 2004).The 2004 samples were also collected on the crater rim, but during a time of background degassing, only mild explosive activity, and a withdrawn free lava lake surface (Gurioli et al. 2008).Brief descriptions, sample locations and modal mineralogy of the complete sample set, including samples from the Stage III Chaimilla violent Strombolian fallout deposit (2.6-1.0 ka; Constantini et al. 2011), and from the Holocene Pucón (Lohmar et al. 2007(Lohmar et al. , 2012) ) and Licán ignimbrites (Silva- Parejas et al. 2010) are presented in Supplementary Table 1.Additional geochemical data were compiled from previous studies, including those of López-Escobar et al. (1977), Deruelle (1982), Harmon et al. (1984), Hickey-Vargas et al. (1989), Sigmarsson et al. (2002) and Witter et al. (2004).

Methods
Seventy-one lineaments of Holocene monogenetic volcanic structures located to the east (Los Nevados satellite cones; area "a" in Fig. 2) and to the south (Chaillupén satellite cones; area "b" in Fig. 2) of Villarrica were measured using Google Earth images (Supplementary Table 3).Fissures were measured directly from the images, while the orientations of the feeder dykes were estimated for each cone following the morphometric approach of Tibaldi (1995).Images were imported and measured in FIJI (Schindelin et al. 2012) equipped with the Strike_Results.ijmplugin (Cortés 2019).Lineament orientations were then imported into Stereonet (Allmendinger et al. 2012;Cardozo and Allmendinger 2013) to generate rose diagrams (Fig. 2a, b).
Whole-rock major and trace-element compositions were determined by X-ray fluorescence spectroscopy at The Open University (UK) and Keele University (UK) using ARL 8420 + dual goniometer wavelength dispersive XRF spectrometers.The estimated reproducibility and accuracy for these analyses, monitored using several in-house rock standards, are better than 1% relative.Loss on ignition (LOI) was determined by heating ~ 1 g of sample material to 900-1050 °C for 1.5-2 h (Supplementary Table 2) and is the total effect of dehydration (loss of weight) and oxidation of ferrous iron (gain of weight).Oxidation often outweighs the effect of dehydration in the analyzed Villarrica samples.
Mineral major element compositions, and major element, S and Cl concentrations in olivine-, clinopyroxeneand plagioclase-hosted melt inclusions and groundmass glasses from scoria from the 1999-2004 period and the 1971 eruption were determined using a CAMECA SX100 electron microprobe at The Open University (UK) with internal PAP correction (Pouchou and Pichoir 1991).Major element glass analyses were obtained using a 20 kV accelerating voltage, 10 nA beam current, a defocused beam (10 µm diameter) and peak counting times per element between 10 and 30 s chosen to minimize the total counting times per analysis.Concentrations of S and Cl were obtained using a 20 kV accelerating voltage, 30 nA beam current, a defocused beam (10-20 µm diameter) and extended peak counting times.
Water and CO 2 concentrations in melt inclusions were analyzed in five olivine-hosted inclusions from the 1971 tephra and in one clinopyroxene-hosted from the background using a Thermo Nicolet Nexus FTIR spectrometer coupled with a Continuμm IR microscope at The Open University.For all spectra, a standard EverGlo mid-infrared source optics, Ge-on-KBr beamsplitter and a liquid nitrogen-cooled MCT-A* detector (11,700-750 cm −1 ) were used.Concentrations of H 2 O were calculated from the height of the total water (i.e., molecular H 2 O + OH − ) peak at 3550 cm −1 , using the Beer-Lambert law: H 2 O (wt%) = 10 × (MA/ρdε), where M is the molecular weight of H 2 O (18.02),A is the height of the absorption peak, ρ is the sample density (g/l), d is the sample thickness (cm) and ε is the molar absorption coefficient (l mol −1 cm −1 ).Sample thicknesses (± 3 μm) were measured using a Mitutoyo Digimatic Indicator.The calculated glass density (2844 kg/m 3 ) is that for the average Villarrica melt inclusion composition measured, assuming a nonlinear temperature dependence of melt volume.For all calculations, a molar absorption coefficient at 3550 cm −1 of 63 l/ mol-cm was assumed.
Fifty-nine melt inclusions hosted in olivine crystals separated from samples from the major eruptive events (1971, 1984 and 2015 eruptions) and from the background activity were measured using the CAMECA ims 4f ion probe (SIMS) at The University of Edinburgh.The instrument is equipped with a Cs microbeam ion source, a He cryo-pump, liquid nitrogen cold trap, and a sample airlock than can hold, pre-pump at working pressures of ~ 3e −6 Pa and bake eight samples prior to analysis.Individual crystals were first mounted in EpoxiCure resin to prevent beam damage and a high background, and then the epoxy blocks were gold-coated prior deploying in the instrument for analysis.
The modal mineralogy and petrography of selected Villarrica samples (Supplementary Table 1) were determined using a combination of optical petrography and high-resolution backscattered electron (BSE) images taken on a Hitachi S-4000 Scanning Electron Microscope (SEM).Crystal size distributions (hereafter CSD; Cashman and Marsh 1988;Marsh 1988Marsh , 1998;;Higgins 2006) of plagioclase was determined using BSE images of representative samples of the different eruptive products (Supplementary Table 1).Images were imported into FIJI (Schindelin et al. 2012), in which best-fit ellipses were applied to determine long and short axes of the different crystals; a minimum of 200 measurements were performed for each mineral phase to have a representative sample of each CSD (Mock and Jerram 2005;Morgan and Jerram 2006).Mean crystal-aspects' ratios were determined using the CSDSlice methodology of Morgan and Jerram (2006).The fabric was considered massive (i.e., with no foliation) and measurements were not corrected for crystal roundness.Resulting data were exported using the CSD output plugin into CSDCorrections 1.4 (Higgins 2006) to calculate the CSD distribution of each mineral phase.The total slide area, volumetric phase abundance and vesicularity of the samples were determined using FIJI's built in features.

Tectonic constrains inferred from monogenetic cone lineaments
A main, N-S prevalent orientation was found in the Chaillupén south monogenetic field (40 measurements), while two prevalent orientations, N-S and NE-SW were found in the Los Nevados, east monogenetic field (31 measurements).The measures are consistent with the orientation of the main trace of the LOFZ and the proposed orientation for the maximum compressive stress associated with the dextral transpression of the structure (Göllner et al. 2021;Simmons et al. 2020a) but parallel to the orientation of the fissures that opened in the main edifice during the 1971 eruption (Castruccio and Contreras 2016).

Whole-rock major and trace element geochemistry
The eruptive products of Villarrica volcano range from calcalkaline basalt to basaltic andesite (50-56 wt% SiO 2 ; Fig. 3, Supplementary Table 2).Juvenile clasts from the Licán and Pucón ignimbrites comprise the more evolved sample cluster (at SiO 2 contents of 54-57 wt%).Subtle major element variations can be observed between the most mafic and more evolved whole-rock compositions (Fig. 4).Fig. 4 Major element-variation diagrams in (wt%) vs SiO 2 (wt%) for Villarrica whole-rock samples.All analyses are normalized to 100 wt% on a volatile free basis.Symbols and source of the data as in Fig. 3 in the mafic rock types with less than 53 wt% SiO 2 and decreasing abundances in more SiO 2 -rich compositions.The high field strength elements (HFSE), such as Zr, Nb behave as incompatible elements throughout the entire Villarrica suite.Among the transition metals, Ni and Cr loosely correlate negatively with SiO 2 .Rb/Zr ratio shows little variation (between 0.14 and 0.22) across all samples and eruptive styles, while the ratio Ni/Zr varies between 0.2 and 0.7 averaging ~ 0.6 for all the products except for the ignimbrites (Fig. 6).

Petrography and mineral chemistry
Texturally, the Villarrica eruptive products range from porphyritic (15-20% phenocrysts) to less porphyritic flow textures (10% of phenocrysts) and are variably vesicular (from 5% in a dense, glassy fragment from the 1984 ice-chilled lava, to more than 56% in a cauliflower bomb from the Pucón Ignimbrite and reticulites from the background activity).
Olivine occurs as a conspicuous phenocryst phase in most of the lava flow samples while it is less common in the samples representing background activity.These phenocrysts are euhedral to subhedral, often with the presence of resorbed rims.For all samples (recent, historical, and Holocene eruptions), Fig. 5 Selected trace element (ppm) variation diagrams vs SiO 2 (wt%).Symbols and source of the data as in Fig. 3 olivine ranges from Fo 69 to Fo 86 with an average of Fo 76 mol% (Figs.7a-d, 8a; Supplementary Table 4a).
Small amounts of clinopyroxene occur as euhedral and subhedral, typically unzoned phenocrysts, mainly detected in the background activity scoria (e.g., VILL-99-3, VILL-99-6, VILL-99-7, and VILL-9910).Representative core-mineral chemistry is reported in Fig. 8b and Supplementary Table 4b.The speciation of iron in clinopyroxene was estimated using charge balance following the approach of Cameron and Papike (1981), while the intrinsic oxygen fugacity was estimated using the relation proposed by Cortés et al. (2006) based on such speciation.Clinopyroxene phenocrysts are mostly confined to the augite field in the pyroxene quadrilateral (Morimoto 1988; Fig. 8a), with two analyses (from sample VILL-99-10) being diopside.
For plagioclase, irrespective of whether the analyses were performed in either core or rim of the crystals, two chemically and texturally distinct populations of phenocrysts are observed in the 1984 samples (Fig. 9).The first population is composed of subhedral, slightly reversely zoned phenocrysts of bytownitic compositions (An 77-86 ), sometimes with ferromagnesian mineral inclusions and melt inclusions that are all indicative of disequilibrium textures.The other population is formed by euhedral to subhedral phenocrysts, which lack both inclusions and zoning, and analyses close to labradoritic compositions (An 60-69 ).
Plagioclase phenocrysts from samples from the 1971 eruption and the Pucón ignimbrite have similar labradoritic compositions (An 58-68 ) and are also texturally similar to those found in the 1984 lava.Tephra from the background activity has phenocrysts with core compositions showing both bytownitic and labradoritic compositions, while the rims, regardless the core composition are close to a labradoritic composition, similar to plagioclases found in the 1971 lava samples (Supplementary Table 4c).
Oxide phenocrysts are euhedral to subhedral with no ilmenite exsolution, consistent with observations by Witter et al. (2004).Oxides in the recent tephra are typically chromian spinel (Cr# ~ 42-62).Following the cation allocation proposed by Sack and Ghiorso (1991), in terms of the molar proportions of end-members, the average oxides found in the samples have a composition of 34% spinel, 30% magnetite, 8% ulvöspinel, and 28% chromite (Supplementary Table 4d).
All plagioclase CSDs lie within a distinct band and range from a steep concave-upward slope to a moderately kinked, straight slope (Fig. 10).Each CSD pattern can be divided into two populations defined by size (> 0.5 mm and < 0.1-0.5 mm), representing two different magma batches.Slope values (based on the straight line calculated with a least square fit) vary between − 1.1 and − 4.7 [mm −1 ]  Based on an average growth rate of ~ 10 -8 mm/sec, consistent with decompression during ascent (Moschini et al. 2023) residence time can be estimated at ~ 100 days for the phenocrysts and ~ 7 days for the smaller population.
Not considering the smaller population, the CSDs of the plagioclase phenocrysts for both recent products and historical lava flows are similar (Fig. 10).

Major element, S and Cl compositions of melt inclusions and groundmass glass
Melt inclusions in the suite of recent tephra occur in olivine, clinopyroxene and plagioclase phenocrysts.They show variable shapes and sizes, typically consisting of light brown glass and, occasionally, show one or more shrinkage bubbles (Wallace et al. 2015).The melt inclusions trapped in olivine phenocrysts tend to be round, triangular, or irregular in shape, with sizes ranging from < 10 µm to a maximum of ~ 100 µm (Fig. 7e).In rare cases, the melt inclusions display evidence of post-entrapment crystallization (Steele-Macinnis et al. 2011) along the inclusion-host mineral interface or in the form of crystals that extend from the wall into the inclusions.Such melt inclusions were avoided during the analyses.Overall, the studied melt inclusions are typically glassy occasionally containing one, or in rare cases two, small shrinkage bubbles, and appear not to be affected by secondary processes after entrapment occurred, such as Plagioclase-hosted melt inclusions tend to be elongated (< 100 µm) and irregular in shape (Fig. 7f).They often appear to be concentrated in the cores of the plagioclase phenocrysts and consist of brown glass, occasionally affected to a minor degree by host mineral crystallization along their boundaries.
A total of 50 melt inclusions in olivine, clinopyroxene and plagioclase and fifty groundmass glasses were analyzed for major elements, sulfur and chlorine (Supplementary Tables 4e, f).Analyses in olivine and clinopyroxene were corrected for post-entrapment crystallization following Putirka (2008) and Simmons et al. (2020a), respectively.After correcting for post-entrapment crystallisation (PEC), reconstructed melt inclusion compositions are basalt to basaltic andesite and show only small variations in major element concentrations (Fig. 11: 52.6-57.9wt% SiO 2 ).They are typically more evolved than the corresponding wholerock compositions, following and extending the trend displayed by the Villarrica whole-rock analyses toward more evolved compositions.The corresponding groundmass glass compositions in the different units vary between 53.6 and 56.1 wt% SiO 2 and are within the same range of values measured in the melt inclusions trapped in the mineral phases within the same sample when compared on an anhydrous basis (Fig. 11).
The concentration of sulfur in melt inclusions varies considerably from ~ 1230 ppm to values at around and below the microprobe detection limit for sulfur (3σ limit ~ 78 ppm).Only ~ 25% of the inclusions were found preserving sulfur concentrations in excess of 250 ppm (Supplementary Table 4f, Fig. 12a).Chlorine concentrations in melt inclusions vary from ~ 130 to ~ 600 ppm, although most inclusions cluster around 200 to 300 ppm of chlorine.Our results also show that in general, inclusions with relatively high concentrations of sulfur also preserve relatively high chlorine contents (Fig. 12a).Groundmass glass analyses of scoria yielded sulfur and chlorine concentrations that are significantly lower than those of the melt inclusions, ranging from below the detection limit (see above) to 190 ppm and from ~ 140 to 440 ppm, respectively.

Water and CO 2 concentrations in melt inclusions and groundmass glass
FTIR results (Supplementary Table 5) suggest water contents between 0.3 and 0.8 wt% in the samples from the background activity, consistent with our SIMS results, and a previously published analysis (a single result of 0.29 wt% using SIMS; Witter et al. 2004).Dissolved CO 2 could not be detected using FTIR in any of the analyzed melt least square fits to the populations that are smaller and bigger than 0.5 mm for the CSD measured in the 1984 lava sample.Regression line equations and the coefficient of determination (r 2 ) are also given in the plot inclusions (detection limit ~ 100 ppm CO 2 ).Of the 59 melt inclusions analyzed using SIMS, 15 are from the background activity, 8 from the 1971 tephra, 5 from the 1984 lava and 16 from the 2015 tephra.Regardless the eruption, the melt inclusions contain on average between ~ 100 and 150 ppm of CO 2 and between 0.4 and 0.9 wt% of H 2 O, with CO 2 up to 560 ppm and H 2 O up to 1.55 wt% (Supplementary Table 5; Fig. 12b).Although in the considered eruptions the range of H 2 O values has a wider spread than in the background activity, there is no relationship between the magnitude of the eruptive style and CO 2 /H 2 O ratios, and high CO 2 content has been measured in some olivine-hosted melt inclusions collected from the background activity (Fig. 12b).
In contrast with similar systems such as Stromboli (Métrich et al. 2001), there is no significant relationship between CO 2 and H 2 O content, with melt inclusions containing high CO 2 concentrations and very low-H 2 O content (Fig. 12b).As most MIs have a large range of H 2 O content at constant and low CO 2 content, i.e., there is no clear degassing path from higher CO 2 concentrations showing a vertical trend at higher H 2 O (Fig. 12b), CO 2 degassing is something likely happening deep in the system, and only the shallower pathway of degassing in which H 2 O gets reduced at low but constant CO 2 has been recorded in the MIs.There are few MIs containing relatively high amounts of CO 2 in relation with H 2 O content, which could be a consequence of: H + diffusion in originally enriched MIs (Gennaro et al. 2019); disequilibrium degassing (Pichavant et al. 2013); diffusive fractionation due to vesiculation and different magma-ascent velocities (Yoshimura 2015); and "CO 2 fluxing", due to the flushing of CO 2 through the shallow system (Yoshimura and Nakamura 2013) or magma mixing and convection (Witham 2011).

Equilibrium between plagioclase and melt
For feldspars, Putirka (2008) introduced a partition coefficient for Na-Ca partitioning between plagioclase and melt.This coefficient is not constant, but temperature dependent (K D ~ 0.25 for T < 1050 °C, K D ~ 0.1 otherwise).The criterion suggests that in general, in the Villarrica magmas, plagioclase phenocrysts are far from equilibrium with Villarrica glass compositions.Only a few calculations in samples VILL-99-7 and VILL-99-1 satisfied the equilibrium criterion, yielding temperatures between 1135 and 1194 °C (Putirka 2008, Eq. 23), and water contents clustering at 0.5-0.9wt% and 1.3-1.8wt% H 2 O (Putirka 2008, Eq. 25a).The values are consistent with some of the measured values using SIMS and FTIR (Supplementary Table 4), in which the measured water content can be considered a minimum value due to possible water loss by diffusion from the melt inclusions.

The Villarrica system
To understand the mechanisms that control the Villarrica system, we first evaluate the overall implications of Villarrica's location in relation to the LOFZ (Fig. 1) and its exerted stress field on the volcano (e.g., Simmons et al. 2020a).Our lineament analysis based on the two monogenetic volcanic fields associated with the main volcanic edifice strongly supports that at least since the Holocene, the orientation of the maximum compressive stress at Villarrica is NE-SW, which is consistent with the overall geometry of the LOFZ main branch running N-S, and with its dextral, strike-slip movement associated with regional transpression (Göllner et al. 2021).Our findings are also consistent with the tectonic model proposed at a regional scale by Cembrano and Lara (2009) or locally, at Quetrupillán, by Simmons et al. (2020a).The main implication is a prevalent orientation for the development of tensile stress exerted on the volcanic edifice and the development of fissures, both at the central cone (e.g., 1971 eruption) and at its flanks.As it has been shown at Llaima volcano (Schonwalder-Angel et al. 2018), the tensile stress orientation promotes the unimpeded ascent of magmas through the volcanic system, consistent with the lack of evolved eruptive products at Villarrica (Fig. 3) and their relatively short residence time in the magmatic system.
Following Edmonds et al. (2022), this local tensional regime is also the cause for long term, persistent, pervasive, and efficient deep degassing at Villarrica.Geobarometry suggests that the depth for a second boiling, i.e., degassing associated with crystallization, is at depths of > 12 km (~ 350 MPa), explaining the lack of CO 2 -rich melt inclusions (their actual concentrations are consistent with depths of less than 1 km; Fig. 12b), which would contain ~ 2500 ppm at depths of ~ 12 km if they were not degassed (e.g., Pan et al. 1991).Melt inclusions hosted in olivine are consistently degassed, regardless of whether they occur in products from the background activity or the major eruptions.This deep and open-system degassing, in conjunction with the low viscosity of the magmas, is likely to be the main mechanism to produce pervasive and efficient mixing of resident and recharging magma (Moussallam et al. 2016;Risso 2018;Boschetty et al. 2022;Edmonds et al. 2022), consistent with some CO 2 -rich/H 2 O-poor melt inclusions (Witham 2011).The mechanical effect of bubble nucleation and growth is likely the agent to promote efficient mass and heat transfer within the Villarrica system, as the buoyant bubbles help to entrain the recharge melt into the resident magma (Wiesmaier et al. 2015).The mechanism is also required to sustain the presence of a permanent lava lake at the vent (Moussallam et al. 2016).The permanent activity of the lava lake and the persistent degassing strongly suggest that the arrival of a volatile-rich recharge is a nearly continuous process in the Villarrica system.
Results from geothermobarometry confirm the presence of two depths in the transcrustal system, where magmas pond and form crystal mush zones located at 9-12.7 km and at around 3-3.5 km, consistent with previous estimations (Morgado et al. 2015;Boschetty et al. 2022).Measurements of plagioclase CSD constrain the timescale for the continuous pervasive mixing between the reservoirs to days/months.

Background activity vs major eruptions
The above results satisfy the dynamics at Villarrica during background activity, i.e., permanent recharge and degassing, pervasive and continuous mixing sustaining the permanent lava lake at the summit, but they do not explain the sudden shift between background activity and major eruptive events.The crystal size distribution of the plagioclase phase and overall, the petrographic textures among similar eruptive products (e.g., tephra form the background activity and from the 1971 eruption) strongly suggest that, compositionally and chemically, there are no differences between products of the background activity and those products from a major eruption.On the other hand, degassing rates measured prior and during the 2015 eruption show a sudden increase in volatile release rates, which has been related to the separate ascent of over-pressured deeply sourced bubbles (Aiuppa et al. 2017;Edmonds et al. 2022).This interpretation is consistent with the low-H 2 O and -CO 2 contents we found in melt inclusions from the 2015 eruption (Supplementary Table 4).In the case of this eruption, or other major historic eruptions, a comparatively more volatile-rich and hotter (but compositionally similar) recharge has been invoked as eruption trigger (Pizarro et al. 2019).
While a substantial change in magma recharge might explain the occurrence of major, paroxysmal events, it fails as a mechanism to explain the background activity, which is permanently sustained between major eruptive events, involving persistent open-system degassing and permanent recharge.An additional issue to consider here is that irrespective of eruption style, the eruptive products are always degassed basaltic andesite.As partial melting of the mantle produces basaltic melt, which evolves while ascending through the transcrustal plumbing system, this constant hybrid composition of the eruptive products further supports a model of continuous injection of evolving melts that subsequently mix efficiently by the mechanical action of the exsolved volatiles (Wiesmaier et al. 2015).Monotonous basaltic andesite compositions at Arenal volcano eruptive products have been explained similarly (e.g., Streck et al. 2002Streck et al. , 2005)).
Therefore, unlike it has been shown at Stromboli (Petrone et al. 2022), we propose that the shift between background activity and a major eruption is not due to the sudden arrival of a substantially more voluminous or different recharge, but rather due to a relatively small increase in the mixing rate and/or the volatile exsolution of a crystal-free recharge into the crystal-mush resident system located at depths of ~ 12-12.7 km.Based on the CSD of the plagioclase phase, such recharge is occurring permanently likely on a timeframe of days to months under the dynamic but unstable equilibrium conditions that sustain the lava lake.Thus, small changes in the recharge dynamics would have dramatic effects on the resulting eruptive style.As modeled by Bergantz et al. (2015), continuous but low-rate magma injection into a resident crystal mush penetrates and spreads through this mush, while high-rate magma injections would fluidise it, which can then be fully transported into the shallower system and eventually evacuating it.A faster injection would also have the effect of the sudden release of volatiles (e.g., Aiuppa et al. 2017), which would also contribute and enhance the evacuation of the shallower magma system.
Our model explains why there are no fundamental changes between the eruptive products from background and paroxysmal activity, and why melt inclusions in both cases contain similar volatile concentrations.It is also supported by the observation that there are several crystal populations in the eruptive products regardless of the type of activity, as reported by Boschetty et al. (2022).Our model also supports the idea that this switch in eruptive activity can be sudden, with little or no warnings (e.g., Stromboli's 2019 paroxysmal activity, Andronico et al. 2021;Petrone et al. 2022), and with similar timescales of the recharge during background activity.

Conclusions
Villarrica is a basaltic andesite volcanic system with shallow and deep reservoirs (~ 3 km and ~ 9-12.7 km, respectively), in which the continuous arrival of a volatile-rich recharge followed by its volatile exsolution and degassing, promotes the continuous and pervasive mixing of their resident magmas.The crystal mush in the lower system is continually but slowly mixed between reservoirs by this process, in which the efficient heat transfer is promoted by the mechanical action of the exsolved volatiles sustaining both the shallow system and the lava lake at the surface during background activity.As magma recharge is a continuous process at Villarrica, its arrival into the system is not what triggers major eruptions, but rather, the small changes in the recharge injection rate and/or degassing rate of the volatiles, which could abruptly mobilize fully the resident crystal mush into the shallow system producing its sudden evacuation switching the eruptive behavior into a paroxysmal event.
Because no substantial changes are required in the magma recharge mechanisms, switching activities might occur in the same timeframe of background activity (i.e., days to few months), posing dramatic implications at Villarrica and similar volcanic systems when forecasting paroxysmic eruptions and their associated hazards.

Fig. 1
Fig. 1 Location of Villarrica volcano (red filled triangle), other volcanic centers in the region (white filled triangles), main settlements (white filled circles) and the two main branches of the Liquiñe-Ofqui Fault zone (red solid lines).The figure shows the Chile-Argentina border with a white dashed line and the inset shows the coverage of the image in relation to the Southern Volcanic Zone of the Andes (see text for more details).Modified from Simmons et al. (2020b); satellite image from Google Earth

Fig. 6
Fig. 6 Rb/Zr versus Ni/Zr diagram of the studied samples and samples from the literature.Symbols and source of the data as in Fig. 3

Fig. 7
Fig. 7 BSE images of the Villarrica volcanic products.a Sample VILL-99-10 (background activity) with plagioclase microlites (plag) and olivines (ol), in a glassy groundmass with 30% of vesicles (vs).b 1947 Lava with phenocrysts of plagioclase and olivine in a groundmass of plagioclase, olivine, and oxides.c Sample VILL-11-05 (1971 Lava) with phenocrysts of plagioclase and olivine with resorbed rims, in a groundmass of olivine, plagioclase, and oxides.d Sample VILL-

Fig. 10
Fig. 10 Natural logarithm of the population density [mm −4 ] versus size [mm] diagram of the CSDs of plagioclases in samples from the background activity (solid blue circles), 1984 lava (solid green triangles) and 1971 tephra and lava (red asterisks).Dashed lines represent

Fig. 11
Fig. 11 TAS diagram of melt inclusions hosted in olivine (circles) and plagioclase (triangles), and of groundmass glasses for background tephra (blue crosses) and for the 1971 tephra (red crosses).Colors are red for the 1971 products, green for the 1984 products, and

Fig. 12 a
Fig. 12 a Plot of chlorine (ppm) vs sulfur (ppm) for melt inclusions hosted in olivine (ol-hosted; solid circle) and plagioclase (plaghosted; solid triangle) and in groundmass glass (asterisk).Colors are red for the 1971 products, green for the 1984 products, and blue for samples from the background (b.g.) activity.b Plot of CO 2 (ppm) , indicate a probability of a relevant eruption in the range of 60-95% within the next 20 years.Villarrica's latest activity includes a large violent Strombolian explosion on March 3, 2015, intermittent Strombolian activity between September 2018 and August 2019, and large explosions until February 2021.