Reader TTS

ABSTRACT Although rare, basaltic Plinian eruptions represent a considerable volcanic hazard. The low viscosity of crystal-poor basaltic magma inhibits magma fragmentation; however, Las

Sierras-Masaya volcano, Nicaragua, has produced multiple basaltic Plinian eruptions. Here, we quantify the geochemistry and volatile concentrations of melt inclusions in samples of the

Fontana Lapilli and Masaya Triple Layer eruptions to constrain pre-eruptive conditions. Combining thermometry and geochemical modelling, we show that magma cooled to ~1000 °C prior to

eruption, crystallising a mush that was erupted and preserved in scoriae. We use these data in a numerical conduit model, which finds that conditions most conducive to Plinian eruptions are

a pre-eruptive temperature <1100 °C and a total crystal content >30 vol.%. Cooling, crystal-rich, large-volume basaltic magma bodies may be hazardous due to their potential to erupt

with Plinian magnitude. Rapid ascent rates mean there may only be some minutes between eruption triggering and Plinian activity at Masaya. SIMILAR CONTENT BEING VIEWED BY OTHERS CRYPTIC

EVOLVED MELTS BENEATH MONOTONOUS BASALTIC SHIELD VOLCANOES IN THE GALÁPAGOS ARCHIPELAGO Article Open access 28 July 2020 THE FORMATION OF TONALITIC AND GRANODIORITIC MELT FROM VENUSIAN

BASALT Article Open access 31 January 2022 THE MAGMATIC ARCHITECTURE AND EVOLUTION OF THE CHANG’E-5 LUNAR BASALTS Article 20 March 2023 INTRODUCTION Plinian eruptions represent the extreme

endmember of explosive basaltic volcanism. The low viscosity of crystal-poor basaltic magma typically precludes magma fragmentation, but basaltic Plinian eruptions have occurred, such as the

60 ka Fontana Lapilli (FL) and 2.1 ka Masaya Triple Layer (MTL) eruptions of Las Sierras-Masaya volcano, Nicaragua; the 122 BCE Etna eruption, Italy; the 1886 CE Tarawera eruption, New

Zealand, and the 2017–2018 CE eruption of Ambae, Vanuatu1,2,3,4,5,6,7. Basaltic Plinian eruptions may eject >1 km3 of tephra and have high mass eruption rates reaching up to 108 kg s−1,

sustaining an eruption column over 25 km in height2,4,5, with substantial impacts for communities and the environment. Examples of sub-Plinian basaltic eruptions include the 1085 CE eruption

of Sunset Crater, USA6. The rapid ascent rates of basaltic magma, estimated between 4–75 m s−1 for the 2017–2018 CE Ambae eruption7, may lead to large volumes of tephra being ejected within

a few minutes of eruption onset, leaving minimal time to implement mitigation procedures. Several mechanisms have been proposed to explain basaltic Plinian eruptions, including syn-eruptive

crystallisation, bubble nucleation during rapid ascent, and high initial volatile concentrations8,9,10,11,12,13,14,15,16,17. Las Sierras-Masaya volcanic system is located approximately 25 

km south from Nicaragua’s capital city, Managua18 (Fig. 1). The present-day Masaya volcano is very active, producing low-intensity, localised explosive activity and is a large source of

magmatic degassing19 in the Central American Volcanic Arc. However, Masaya caldera has produced several basaltic Plinian eruptions: the 6 ka San Antonio Tephra (SAT), the 2.1 ka MTL and 1.8 

ka Masaya Tuff-Ticuantepe Lapilli eruptions (MT-TIL)13,20. The 60 ka FL eruption has been associated with the formation of the older Las Sierras caldera, after which the Masaya shield

volcano was constructed18,20. Deposits from all these eruptions are found at the present-day location of Managua2,3,13,21. The last major eruption of juvenile material at Masaya was the 1772

CE lava flow20. Since 1772, approximately 47 km3 of non-erupted magma is estimated to have been degassed20. Despite the diverse eruptive activity at Las Sierras-Masaya volcanic system,

whole-rock compositions of erupted products appear to show chemical homogeneity over 60,000 years, suggesting highly explosive activity is not controlled by variations in magma

composition20,22,23. This absence of geochemical cyclicity makes the timescale of a future Plinian eruption difficult to forecast as the transition to an explosive system may not be preceded

by an identifiable geochemical change, increasing the volcanic risk posed to the ~2 million people who live within the area24. The aim of this study is to understand the magmatic conditions

that could drive a highly explosive eruption of this hazardous volcano. We present a combined petrological and numerical investigation (detailed in Methods) into the conditions which

promote Plinian activity at this system, providing new major, trace and volatile element data for the FL and MTL eruptions. We demonstrate that the FL magma experienced a complex

pre-eruptive history, stagnating at shallow depth, enabling cooling and crystallisation of a mush within the magma reservoir. We incorporate petrological data on the MTL eruption from this

study and Bamber et al.16 into a numerical conduit model25,26,27 to provide simulations of the conduit dynamics of a basaltic Plinian eruption of Masaya. We correlate the results of our

sensitivity analysis to the diverse eruptive history of the volcano. RESULTS AND DISCUSSION THE FONTANA LAPILLI ERUPTION: TEXTURAL OBSERVATIONS AND GEOCHEMICAL MODELLING FL samples are

basaltic andesitic (~53 wt.% SiO2; ~4 wt.% total alkalis), with plagioclase (7 vol.%), olivine (1 vol.%) clinopyroxene (<1 vol.%) phenocrysts and trace amounts of magnetite. Microlites

(<100 µm) often have hopper to swallowtail habits. Microlite phases are plagioclase, clinopyroxene, olivine and magnetite. Scoriae are texturally heterogeneous, with microlite-poor and

microlite-rich regions (Fig. 2a), where crystallinity increases from 5 to 63 vol.% and glass varies between 95 and 37 vol.%. Vesicularity varies between 14 and 26 vol.% for the two

endmembers. This microtextural heterogeneity is comparable to MTL samples (Fig. 2b–d). FL samples contain enclaves with a mineralogical assemblage similar to the surrounding matrix (Fig. 

2e–h). We define the matrix as the assemblage of phenocrysts, microlites and glass described above. Enclaves have a higher plagioclase and clinopyroxene crystal fraction than the matrix

(Supplementary Table 1) and a different habit than microlite matrix crystals despite their comparable size. Enclave crystals may have a crystallisation history distinct from matrix

phenocrysts and microlites. Vesicles within enclaves are smaller and more spherical than matrix vesicles. We define two categories of enclaves based on microtextural characteristics,

referred to as Type A and Type B. Type A enclaves have features indicative of plastic deformation28, such as irregular contacts with the matrix, subparallel alignment of the mineral

assemblage with matrix microlites, and a large interstitial glass component (41–51 vol.%) (Fig. 2e–g). Spherical to sub-spherical vesicles in interstitial glass suggest that Type A enclaves

contained dissolved volatiles, vesiculating during magma ascent29. Type B enclaves (Fig. 2h) have regular shaped contacts with the matrix, crystal alignment discordant to that of matrix

microlites, and less interstitial glass (11–16 vol.%) with very few vesicles. The area fraction of enclaves ranges between 0.04 and 0.26 and the area of enclaves ranges between 0.3-0.6 mm2.

Plagioclase and clinopyroxene (phenocrysts and microlites) compositions within the matrix (An51-90; Wo28-46En35-51Fs8-23) and enclaves (An22-69; Wo30-44En34-51Fs12-23) in FL samples overlap

(Supplementary Fig. 1). Enclave hosted olivine crystals are more fayalitic (Fo45-68) than matrix microlites (Fo68-73). Enclave and matrix crystal compositional homogeneity suggests a

possible common origin, as phases accumulating within a crystal-rich mush zone of the magma reservoir30. Enclave phases tend to show more evolved compositions, suggesting that they may have

crystallised at lower temperatures within the reservoir. Preferred orientation of crystals resembles the sedimentation structures of cumulates (Fig. 2g), trapping interstitial liquid.

Preferred orientation of crystals may also indicate flow alignment, but lamination and randomly orientated, sub-spherical vesicles in interstitial glass may be more suggestive of crystal

settling, as vesicles are not elongated in the direction of preferred orientation31. We interpret Type A enclaves as representative of the mush zone within the storage reservoir32. Type B

enclaves have less interstitial glass and occasional small, sub-spherical vesicles and are more likely sourced from the adjacent rigid crust of the mush32. The distribution of Type B

enclaves throughout the FL eruption phases suggests continued disaggregation of the mush, as opposed to entrained plutonic fragments which may be concentrated in the opening phase. The mush

and its interstitial liquid were preserved within scoriae on eruption and quenching29. Most FL melt inclusions (MI) are more evolved than the whole-rock composition, showing a decrease in

MgO and an enrichment in SiO2 and alkalis (Fig. 3a–c). MI have been corrected for post-entrapment crystallisation (PEC). Following PEC corrections, most MI are still not in equilibrium with

both their host crystals and more anorthite-rich plagioclase crystals (An73-90) (Supplementary Data 1). FL microlite-poor glass and MI compositions generally overlap. However, some

microlite-poor glasses are trachy-andesitic, overlapping with microlite-rich glass. Type A enclave interstitial glasses are generally more evolved than FL matrix glass but overlap with

trachy-andesitic, microlite-poor glass. Type B enclave interstitial glasses have the most evolved and varied compositions. We used Rhyolite-MELTS v.1.233 to examine the expected liquid

compositions that would be generated from the average FL MI composition with progressive cooling and crystallisation. FL compositions measured using EPMA plot upon a common modelled liquid

line of descent (LLD) (Fig. 3b, c), indicating melt evolution controlled by plagioclase, clinopyroxene and olivine crystallisation, reaching a total crystallinity of 60 vol.% at 960 °C and

40 MPa (Supplementary Data 2). The modelled crystal fraction at 40 MPa and 980–990 °C agrees with measured Type A enclave crystal fractions (Supplementary Table 1). Figure 4 shows the

distribution of trace elements in glasses and MI (Supplementary Data 3). For FL samples, Sr decreases with increasing Zr, as Sr is compatible in plagioclase, the dominant crystallising

phase. Incompatible elements Ba and La increase with increasing Zr. A fractional crystallisation trend was modelled using Petrograph software34 (Fig. 4). FL MI and glasses fit well with the

modelled fractional crystallisation trend, indicating crystallisation-controlled melt evolution. According to the fractional crystallisation model, the observed trace element distribution is

produced by 50–60 vol.% crystallisation, in agreement with Type A enclaves that have crystal contents of 50–60 vol.%. The microlite-poor, trachy-andesitic melt enriched in incompatible

elements but depleted in Sr could originate from expulsion of interstitial melt from the crystal mush on its compaction35, or upon volatile saturation31. Melt differentiation models for both

FL major and trace element chemistry confirm crystallisation-controlled melt evolution and that it is unlikely that the observed trends were produced by magma mixing36. PEC-corrected FL MI

contain 0.91–2.27 wt. % H2O and 98–1103 ppm CO2 (Fig. 5). For FL MI, H2O and CO2 concentrations generally decrease with decreasing MgO content. In general, more primitive MI are found in

more anorthitic phenocrysts. Calculated FL MI entrapment pressures range between 26 and 189 MPa, assuming volatile saturation at the time of entrapment. However, the CO2 budget of MI-hosted

bubbles has not been measured. Formation of bubbles during MI cooling and PEC can partition CO2 into the bubble from the glass phase37. MI–hosted bubbles may host up 40–99% of the total

inclusion CO237,38,39,40; therefore, the MI glass phase provides only a minimum estimate of the pre-eruptive CO2 content. In turn, our entrapment pressures represent minimum estimates of the

storage region depth. This entrapment pressure range is consistent with Rhyolite-MELTS simulations, which find that the FL mineralogical assemblage cannot be produced at pressures exceeding

200 MPa. At higher pressures, clinopyroxene is the dominant crystallising phase (Supplementary Data 2), but clinopyroxene phenocrysts are rarely observed in our samples. The LLD modelled at

a fixed pressure of 200 MPa does not fit our measured glass compositions (Fig. 3). Modelling a degassing trend using MagmaSat41,42 finds that FL MI plot upon a closed-system trend with an

initial exsolved vapour content of 5 wt.% (Fig. 5). However, analysis of bubble-bearing MI demonstrates that a degassing trend may be produced by measuring only the MI glass phase and not

accounting for the bubble CO2 content37. For our FL samples, only one MI has a bubble; therefore, our modelled trend most probably represents a true degassing trend rather than being

produced by CO2 loss to MI-hosted bubbles. Our MI may have been affected by H+ diffusion, which may re-equilibrate H2O concentrations. We selected MI from small, rapidly quenched scoriae 1–3

 cm in size to minimise the effect of post-fragmentation H+ diffusion43. H+ diffusion can occur prior to fragmentation at pre- and syn-eruptive conditions. The smallest MI closest to the rim

of the smallest phenocryst hosts will be the most susceptible to diffusive loss of H2O43,44. Indeed, MI closest to the phenocryst rim show the lowest H2O contents. However, these lower H2O

contents show some overlap with larger MI hosted in phenocryst cores, which would be less prone to diffusive H2O loss. The short period of storage and rapid magma ascent may have limited the

timescale for diffusive re-equilibration. FL MI likely experienced only minimal diffusive H2O re-equilibration within the storage region, preserving the observed degassing trend. This

exsolving vapour phase may have driven the ascent dynamics within the conduit, imparting increased buoyancy and gas pressure. THE FONTANA LAPILLI PRE-ERUPTIVE CONDITION As mineral-melt

thermometry is sensitive to melt H2O concentration, we estimate pre-eruptive temperatures using the maximum H2O content and the average H2O content of 1.59 wt. % calculated from FL MI in

this study and other works11,45. Using the maximum H2O content, we find matrix FL plagioclase phenocrysts and microlites crystallised within the temperature range 1050–1062 °C (Supplementary

Fig. 2). The standard error of the Putirka46 thermometer (Eq. 23) is ±43 °C. Pyroxene phenocrysts crystallised within the range 1063–1075 °C and microlites between 1051–1103 °C. The

standard error of the Putirka46 thermometer (Eq. 33) is ±45 °C. The majority of matrix less anorthitic plagioclase and pyroxene phenocrysts and microlites (_n_plg = 88/118; _n_pyx = 11/11)

are in equilibrium with the average FL MI composition, using the equilibrium criterion of each thermometer. Our temperature range agrees with experimentally determined pre-eruptive

conditions for the FL eruption of 40–80 MPa and 1010–1060 °C11. Our estimated pressure range from Rhyolite-MELTS simulations (30–40 MPa) agrees with the lower end of this experimentally

determined pressure range11. Costantini et al.10 find higher pre-eruptive temperatures of 1108–1130 °C for clinopyroxene crystallisation but used a lower input H2O content for their

thermometry calculations. Type A enclave plagioclase crystallised over a lower temperature range of 996–1006 ± 43 °C, and pyroxene between 988–1035 ± 45 °C. Type B enclave phases

crystallised between 966–1036 °C. Enclave phases are in equilibrium with the average Type A and Type B enclave glass compositions respectively. Enclave crystals are small and euhedral and

most likely crystallised towards the higher end of this range at low undercooling. Using a lower H2O content produces a higher calculated temperature for plagioclase (1076–1089 °C) and

pyroxene (1060–1112 °C) phenocrysts and microlites. Type A enclave phases crystallise between 995–1044 °C and Type B phases between 973–1045 °C. However, calculated temperatures are still

mostly <1100 °C. Using the simulated output crystal fraction, measurements of phenocryst fraction (Φplg = 7 vol.%; Φol = 1 vol.%; Φcpx = <1 vol.%) and the results from mineral-melt

thermometry16,36, we find the FL magma that crystallised the matrix phases was last stored at pressures of 30–40 MPa and temperatures of 1050–1060 °C (Supplementary Data 2). The measured

enclave crystal fraction can be reproduced at pressures of 30–40 MPa, but between 990 and 1000 °C. Textural observations, mineral-melt thermometry and geochemical modelling of melt evolution

suggest that enclaves crystallised from a melt genetically related to that which crystallised the matrix phases but represent a more prolonged stage of cooling and crystallisation. Enclave

phases formed within the same pre-eruptive storage region as matrix phases but at an earlier point in time, crystallising in contact with their interstitial liquid and may be considered

antecrystic47. Unzoned, euhedral phenocrysts and microlites such as low-An plagioclase and pyroxene are in equilibrium with the average MI composition and can be considered autocrystic,

forming from the last batch of melt prior to eruption47. High-An plagioclase crystals are not in equilibrium with the average MI composition and may originate from deeper within the plumbing

system as part of a second antecrystic crystal cargo injected into the magma chamber prior to eruption. However, it can be difficult to distinguish between autocrysts and antecrysts in

mush-bearing systems due to the remobilisation of crystal cargo47. At depth, magma chambers likely do not exist as a single, continuous body of melt, but consist of crystalline mush and

transient, eruptible melt lenses32. Cooling rates within the magma chamber would be greatest at the margin, where there is greater thermal contrast with country rock. Crystallisation would

proceed first along the chamber walls and floor where magma temperatures are lower48. Settling of crystals would occur towards the chamber floor48. Numerical simulations of a shallow,

differentiating, basaltic magma reservoir find convection transfers mass and heat, producing compositionally different liquids and temperature contrasts of several hundred degrees early

within the reservoir’s cooling history48. Crystallisation, assembly and mobilisation of crystal mushes can occur over days to months at basaltic systems49,50. Crystal Size Distribution (CSD)

calculations show enclave plagioclase could have crystallised within several days (Supplementary Fig. 3; Supplementary Table 2). CSD calculations provide a minimum residence time, as

settling of mush-forming phases in the chamber is not considered. Crystallisation and volatile exsolution could provide sufficient overpressure within the magma chamber to trigger an

eruption36,51,52. We suggest that a batch of magma entered the reservoir, stagnating and crystallising a mush at the edge of the reservoir at temperatures of ~1000–1040 °C and 30–40 MPa over

several days. A fresh batch of volatile-rich, higher-temperature magma (>1050 °C) then entered the reservoir, where it was stored for only a short period of time at 30–40 MPa,

crystallising a small plagioclase-dominated assemblage (Fig. 6). The addition of volatiles or higher-temperature melt into the intra-crystalline mush network likely induced melting and/or

internal overpressure, unlocking the mush and enabling its eruption35,49. The closed-system degassing trend modelled for the FL eruption suggests fresh magma intruding the reservoir was

volatile-rich, providing overpressure within the sealed storage region that could trigger an eruption. As plagioclase-melt equilibrium is sensitive to melt H2O content, disequilibrium

between MI and their anorthitic hosts suggests that plagioclase hosts crystallised from a melt with a higher H2O content within the deeper plumbing system, similar to the SAT, MTL and MT-TIL

eruptions13,23. FL MI which record entrapment pressures >40 MPa may represent the volatile-rich magma injection originating from greater depth. These MI may have had insufficient time to

re-equilibrate during the short storage time within the shallow reservoir, preserving higher volatile concentrations and entrapment pressures. COMPARISON WITH OTHER ERUPTIONS OF LAS

SIERRAS-MASAYA VOLCANIC SYSTEM We compare our FL pre-eruptive scenario with the SAT and MTL eruptions. SAT MI have low H2O contents of 1.03–2.22 wt.%13, comparable to FL MI contents.

Similarly, SAT samples show both basaltic-basaltic andesitic and andesitic glass. Pérez et al.13 calculated that the SAT andesitic matrix glass (Fig. 3a) can be produced by crystallisation

of 25% plagioclase, 15% clinopyroxene, 3% olivine and 5% magnetite, comparable to our measured enclave crystallinity (Supplementary Table 1). The SAT andesitic glass is interpreted as a

residual melt from crystallisation and differentiation of a batch of SAT basaltic magma at the reservoir walls, consistent with our FL pre-eruptive condition. Although cumulates were not

observed in SAT samples, the andesitic magma indicates cooling, crystallisation and melt differentiation within the storage reservoir, producing a glass composition comparable to our FL

enclave glass13. MTL samples also show microtextural heterogeneity16 (Fig. 2b–d). However, there is minimal evolution concerning MTL major (Fig. 3a) and trace element geochemistry (Fig. 4).

Compared to FL MI, a narrower range of 0.85–1.25 wt. % H2O and 129–457 ppm CO2 is measured in MTL PEC-corrected MI (Fig. 5). Our results for the MTL eruption are comparable to Hlinka et

al23. but are lower than the maximum H2O content of 2 wt. % measured for MTL MI using Raman spectroscopy16. The small size of some MTL MI may indicate that they were susceptible to diffusive

H2O loss, despite the small size of scoriae43. Diffusive re-equilibration may have produced the low MTL MI H2O contents. Higher H2O contents of 5–6 wt. % H2O have been suggested for the

SAT, MTL and MT-TIL eruptions13,23. MI then experienced re-equilibration of H2O during a short storage period within the shallow reservoir13,23. Thus, the homogeneous H2O contents measured

in MTL MI most likely represent the magmatic conditions prior to eruption. Calculated MTL MI entrapment pressures range between 33–83 MPa, but most cluster at lower values. As we measure

only the MTL MI glass phase, reported CO2 contents and calculated entrapment pressures represent minima. We observe only a small variation comparing the glass CO2 content of bubble-bearing

(54% of MTL MI) and bubble-free MI, suggesting negligible CO2 loss to vapour bubbles, consistent with previous work23. Our maximum measured CO2 content of 457 ppm is consistent with previous

work which finds that most CO2 contents do not exceed 400 ppm for the MTL, MT-TIL and recent eruptions23. FL MI cluster at the lower end of H2O and CO2 contents, with entrapment pressures

<100 MPa. However, FL MI record higher entrapment pressures than MTL MI within the most primitive MI (Supplementary Fig. 4). The preservation of a degassing trend and a wide range in

volatile contents in FL MI suggests that the FL magma was stored for less time than the MTL magma, limiting the extent of degassing and homogenisation of melt compositions. A more prolonged

period of storage and pre-eruptive degassing for the MTL eruption may have produced the invariant, lower MI volatile concentrations and homogeneous melt chemistry, obscuring any original

degassing trend (Fig. 6). Storage within the shallow reservoir, degassing and fractionation at Masaya has produced geochemically homogeneous erupted products over time20,22. The FL and SAT

basaltic-andesitic and andesitic glass suggests that these magmas were stored for shorter durations at shallow depth compared to the MTL and MT-TIL magmas. Evidence suggests that the Las

Sierras-Masaya magmas were stored within a reservoir at shallow depth (<4 km) prior to eruption, degassing and crystallising the mineral assemblage13,16,20,22,23. Magma stagnation and

degassing at shallow crustal depth is suggested as a consistent pre-eruptive stage at Masaya regardless of eruptive style, producing the low volatile contents of Masaya MI23. Our

pre-eruptive scenario for the FL eruption is consistent with this model. The FL eruption had the lowest pre-eruptive temperature of the Las Sierras-Masaya Plinian eruptions. The more recent

Plinian eruptions show higher pre-eruptive temperatures of 1075–1161 °C (SAT), 1080–1100 °C (MTL) and 1080 °C (MT-TIL)13,16. Cooling within the magma reservoir can lead to crystallisation of

a mush, as observed for FL samples. Back-scattered electron (BSE) images of MTL enclaves suggest inclusions of mush may have been erupted (Supplementary Fig. 5). The transition from

steady-state open conduit activity to a Plinian eruption is suggested to be modulated by sealing of the conduit, producing a closed system and build-up of overpressure23. We suggest that

cooling and crystallisation of a mush at shallow crustal depth could provide a mechanism for sealing the system. The mush, together with the injection of volatile-rich magmas into the

storage region, may contribute to the increase in overpressure of the system and assist in unlocking the mush35 (Fig. 6). Therefore, the eruptive style at Las Sierras-Masaya volcanic system

is closely related to the initial, pre-eruptive conditions. Although brittle magma fragmentation is most likely driven by syn-eruptive microlite crystallisation and restricted outgassing

within the conduit15,16, the initial, pre-eruptive conditions may promote this scenario by producing a lower temperature, higher viscosity magma more susceptible to fragmentation within the

conduit. The relative importance of these initial conditions and their control on eruptive style can be investigated using numerical modelling. CONDUIT DYNAMICS OF A BASALTIC PLINIAN

ERUPTION We used the 1D, multiphase, steady-state conduit model of La Spina et al.25,26,27 to simulate a basaltic Plinian eruption at Masaya caldera. The governing equations driving the

multiphase flow are those illustrated in La Spina et al.27 and for completeness have been reported in the Supplementary Methods. The constitutive equations needed to characterize the magma

ascent dynamics at Masaya caldera are described in the Methods section. The initial conditions (Supplementary Table 3) are informed by the results of this study and Bamber et al.16. The

maximum measured MI CO2 content from this study is used as the initial dissolved CO2 concentration. Initial pressure, temperature, crystal content and dissolved H2O concentration values are

from Bamber et al.16. We use the strain-rate criterion of Papale53 to simulate fragmentation, where failure occurs when the Deborah number exceeds 0.0153. We set the characteristic time of

crystallisation to 10 s, consistent with previous work on the Etna 122 BCE Plinian eruption15. Our reference simulations focus on characterising the MTL conduit dynamics. The MTL

pre-eruptive conditions likely represent the current system and future possible Plinian activity. The MTL vent location is thought to lie within the present-day Masaya caldera, whereas the

FL vent location is not well known but thought to lie outside the modern caldera2. The FL crystal cargo, with fragments of mush, makes it difficult to assign an initial crystal content and

to describe the syn-eruptive crystallisation process. Instead, the FL pre and syn-eruptive conditions constrained by this study are explored through a sensitivity analysis. We address the

rheological behaviour of the heterogeneous magma by defining two endmembers from natural observations: a microlite-poor and microlite-rich case (Fig. 2b–d). Here, crystal-poor and

crystal-rich simulations refer to the microlite-poor and microlite-rich endmembers defined from observations, as we do not consider enclaves within our simulations. We vary the initial

conditions until the output reproduces the observed crystal content in our samples, which we take as representative of the crystal content at the vent. The initial crystal fraction (_φ_) and

the initial pre-eruptive temperature (_T_) differ between the crystal-poor (_T_ = 1090 °C; _φ_ = 10 vol. %) and crystal-rich (_T_ = 1070 °C; _φ_ = 17 vol. %) simulations. The initial

temperatures agree with the estimated MTL pre-eruptive temperature range, accounting for the standard error of mineral-melt thermometry16. The initial crystal fraction of the crystal-poor

endmember agrees with the observed MTL phenocryst fraction. Although the crystal-rich endmember is produced using a higher initial crystal fraction, the agreement between the output crystal

fraction and microtextural observations indicates a higher crystal fraction may have been needed to produce the crystal-rich endmember. The crystal-poor simulation (Fig. 7, solid lines)

produces a total exit crystal fraction of 0.24, and the crystal-rich simulation (Fig. 7, dashed lines) a total fraction of 0.51, comparable with the measured crystal fractions for the

microlite-poor and microlite-rich regions in MTL samples of Plinian phase TLL216. The crystal-poor simulation produces an effusive eruption, as the critical Deborah number is not reached

during magma ascent (Fig. 7a) and there is no fragmentation within the conduit. For the crystal-rich simulation, fragmentation occurs ~280 m below the vent (Fig. 7a). Syn-eruptive

crystallisation in the crystal-rich simulation drives a considerable increase in mixture viscosity (Fig. 7b), reaching the critical Deborah number for fragmentation. The lower viscosity,

crystal-poor endmember (Fig. 7b) enters the conduit at a higher velocity on the order of 10 m s−1 (Fig. 7c), producing a higher mass flow rate of 1.1 × 108 kg s−1 (Fig. 7d). The higher

viscosity crystal-rich endmember enters the conduit at lower velocity on the order of 1 m s−1, resulting in a lower mass flow rate of 6 × 106 kg s−1. Gas-melt coupling is maintained during

ascent for the crystal-poor endmember and prior to fragmentation for the crystal-rich endmember (Fig. 7c). Both mass eruption rates lie within the range in mass eruption rate of 106−108 kg 

s−1 indicated by the MTL deposit3. The average decompression rate is 0.91 MPa s−1 and 0.05 MPa s−1 for the crystal-poor and crystal-rich simulations respectively. Our simulated decompression

rates are comparable to estimates of the FL decompression rate from experiments and calculations using bubble number densities, which are on the orders of 0.1 and 1 MPa s−1

respectively12,54. The rapid ascent velocity results in a short ascent time, between 1–13 min considering both endmember simulations. The calculated timescale is consistent with CSD analysis

of the MTL samples, which suggests rapid syn-eruptive crystallisation within minutes during magma ascent16. Experiments demonstrate that rapid syn-eruptive crystallisation can occur within

minutes at high undercooling, driving basaltic magma fragmentation during ascent15. The lower initial temperature of 1070 °C for the crystal-rich simulation provided a higher undercooling,

facilitating rapid crystallisation and magma fragmentation. The temperature of the crystal-rich simulation increases during magma ascent due to the release of latent heat of crystallisation.

The amount of latent heat produced by syn-eruptive crystallisation exceeds the decrease in temperature produced by adiabatic cooling during ascent, which is instead observed for the

crystal-poor simulation. Heterogeneous textures can reflect portions of the same initial magma that experienced different extents of microlite crystallisation55. Microlite-rich domains are

expected to experience a more extended period of crystallisation, ascending at lower velocities and in cooler zones of the conduit55. A 40 °C decrease in the temperature of crystallisation

may have produced crystalline tachylite glass as opposed to sideromelane during the 2001 CE Etna flank eruption56. Our simulations demonstrate the two endmembers defined from natural

observations can be produced by decreasing the initial temperature by 20 °C and increasing the initial crystal content by 7 vol.%. With a small change in temperature and initial crystal

content (Fig. 7e, f), a large micro-textural difference can be produced during magma ascent. This result suggests that the two textural endmembers can be produced from a small difference in

initial conditions, which then experience different degrees of syn-eruptive crystallisation during ascent. This difference in crystallisation produces a large viscosity contrast between the

endmembers of 104−107 Pa s. We use the magma ascent model to estimate the magma reservoir volume that supplied the MTL eruption. We find calculations using a chamber volume of 10 km3 are

consistent with the duration (~1 hr) and the minimum erupted volume (0.038 km3 DRE) estimated from land outcrops of MTL Plinian phase II3, where the duration is estimated from the mass flux

and erupted mass of the deposit (Fig. 8). The total MTL erupted volume is ~2 km3 (DRE) considering offshore deposits13. The rate of magma flux into the present-day Masaya reservoir has been

estimated at ~0.19 km3 yr−1 20. Using our chamber volume estimate of 10 km3, the volume of the Plinian eruptions will have accumulated within a minimum of ~50 years. Although mush may

crystallise rapidly at Masaya, a critical crystal fraction may need to be reached within the reservoir for incoming, volatile-rich fresh magma to generate sufficient overpressure to trigger

the eruption. As crystallisation likely proceeds from the reservoir margin, it may take years to assemble a magma volume capable of producing a Plinian eruption. SENSITIVITY ANALYSIS We

perform a sensitivity analysis to identify the most important parameters driving explosive and effusive activity at Las Sierras-Masaya volcanic system and their effect on output parameters

such as mass eruption rate, exit velocity, total crystal content and mixture viscosity. We selected a range of values for the input parameters which incorporate the pre-eruptive conditions

determined for the eruptions of Las Sierras-Masaya volcanic system. We use our constraints on initial temperature, crystal content, and dissolved H2O and CO2 from this study for the FL

eruption. Our sensitivity analysis uses the average MTL MI composition, which we assume is representative of the melts being delivered to Masaya’s shallow magma chamber due to the

compositional homogeneity of erupted products through time20. Explosive solutions (Fig. 9a–j, square symbols) occur when the critical Deborah number is reached, resulting in brittle magma

fragmentation. In our sensitivity analysis, <20% of solutions are explosive, consistent with the rarity of highly explosive eruptions at this system. Non-explosive solutions (Fig. 9a–j,

diamond symbols) represent a wider range in eruptive style27, which we define using the average exit velocity (represented by the colour of the diamond symbols). Model solutions for which

the fragmentation threshold is not reached, but exit the vent at a high ascent velocity, likely represent lava fountaining activity27. Lava fountain activity is expected to occur when exit

velocities exceed 1 m s−1 for low fountain activity (jet height is 0.1–50 m) and 30 m s−1 for high fountain (jet height >50 m) activity27. Non-explosive solutions with a lower exit

velocity likely represent effusive eruptions such as lava flows. As we adopt a small characteristic time of crystallisation, we may overestimate the number of explosive solutions15. The

initial temperature has a strong control (Fig. 9a) on the frequency of explosive vs. non-explosive solutions (Supplementary Fig. 6a). As initial temperature decreases from 1100 °C, the

frequency of explosive solutions increases, reaching a maximum of 80% at temperatures of 1060 °C, due to the increase in melt viscosity with decreasing temperature and the control of

temperature on syn-eruptive crystallisation (Supplementary Fig. 6b). Following this trend, we can infer that at lower temperatures the percentage of explosive simulations would exceed 80%.

There is a positive correlation between increasing temperature and the frequency of non-explosive solutions. Only non-explosive solutions are observed at temperatures exceeding ~1100 °C.

This result agrees with the FL and MTL pre-eruptive temperatures of <1100 °C constrained by mineral-melt thermometry16. Combining these results, we suggest that a low pre-eruptive

temperature may have been a driving cause of Plinian activity at this volcanic system. The decrease in exit velocity with decreasing temperature indicates that a transition from

high-velocity to low-velocity effusive activity may signify a future shift to explosive activity. The conduit radius has minimal influence on the relative frequency of

explosive/non-explosive solutions (Fig. 9c) but exerts a control on the mass flow rate (Supplementary Fig. 6c). As a high mass eruption rate is required to produce a sub-Plinian to Plinian

eruption (106–108 kg s−1), we expect most explosive solutions of Plinian magnitude to occur at larger values of the conduit radius. The increased frequency of Plinian activity at Masaya

following the SAT eruption may be attributed to the SAT eruption establishing a larger conduit radius. Following the MT-TIL eruption, activity has been largely controlled by the fault

structures remaining from caldera formation and the extensional stress exerted by the Managua graben. The 1772 CE lava flow was erupted from a fracture linked to the Cofradías fault57. The

sensitivity analysis shows a small increase in the frequency of explosive solutions between H2O contents of 1 and 1.5 wt.% as volatile-saturated conditions are approached (Fig. 9d). The

frequency of non-explosive eruptions increases in water under-saturated conditions. Once water-saturated conditions are exceeded (~1.5 wt.% H2O), the frequency of explosive eruptions

slightly decreases with increasing H2O content. These results suggest that high H2O concentrations may not be required to produce a Plinian eruption at this system, but only sufficient H2O

to establish volatile-saturated conditions. Magmas stored at shallow crustal depth may be able to erupt explosively if volatile-saturated conditions are reached. The maximum measured MTL

pre-eruptive H2O content16 agrees with these results. High initial CO2 concentrations have been implicated as a driving mechanism of explosive basaltic activity58,59. Even if our

measurements of MI CO2 represent minima, the sensitivity analysis results show that the frequency of explosive solutions is not particularly sensitive to the initial CO2 content, even up to

10,000 ppm (Fig. 9e). Our measurements and numerical simulations suggest high initial volatile concentrations were not a critical condition promoting Plinian activity at this volcanic

system. The frequency of explosive solutions increases with increasing initial crystal content (Fig. 9f), due to the increase in mixture viscosity. We can assume that enclaves would increase

the likelihood of an explosive eruption. Numerical results show explosive solutions start to occur when the total crystal content exceeds 30 vol.%, as a higher crystal content will increase

mixture viscosity (Fig. 9g). FL, MTL and MT-TIL samples show high crystal contents, between 48 and 64 vol.%10,16,23. Instead, present-day lava lake ejecta shows low microlite contents of 5

vol.%23. The exit velocity and the degree of gas-melt coupling show a relationship with eruptive style (Fig. 10). Numerical results show that at high values of Reynolds number (>102) and

low viscosity (<104 Pa s) magma is not able to fragment within the conduit and exits the vent at high speed (on the order of 10–102 m s−1), producing lava fountaining activity27. At

higher viscosity (between 105 and 108 Pa s), the exit velocity decreases and both explosive and effusive activities become more likely. Results show that a strong coupling between gas and

melt is not required to produce an explosive eruption, in agreement with La Spina et al.27. As the viscosity increases even further (>108 Pa s), magma ascent rate continues to decrease,

and gas can decouple efficiently from the melt, preventing the formation of high strain rates, and thus producing only effusive eruptions. Although we account for the heterogeneity of the

MTL magma by defining two endmembers in our simulations, we recognise that there could be a complex, nonlinear behaviour during magma ascent. The effect of the two intermingling endmembers

and their relative contribution to bulk mixture viscosity is difficult to constrain. In a magmatic system, the mass eruption rate and viscosity ratio between the two endmembers can control

the development of instabilities and physical mixing60,61. The viscosity ratio between the endmembers is expected to be large due to their differences in temperature and crystallinity16,

where viscosity increases from 102 to 106 Pa s. Future work investigating the impact of the physical mixing of two endmembers with different crystallinities on bulk mixture viscosity is

needed. For our simulations, the large viscosity ratio between the endmembers should produce turbulent flow, enabling efficient syn-eruptive mingling in the conduit61 but not sufficient to

homogenise the two endmembers, preserving the heterogeneous textures observed in the samples. Low-viscosity basaltic systems such as Kīlauea rarely produce explosive eruptions without

external water27. A numerical investigation found that the range of pre-eruptive conditions at Kīlauea, with temperatures ~1150 °C, crystal contents <5 vol.% and H2O contents <1 wt.%,

are not favourable to explosive eruptions, instead producing fountaining activity27. The Las Sierras-Masaya magmas are lower temperature, more crystalline and have a higher volatile content.

They are therefore more likely to be higher viscosity. Although pre-eruptive temperatures for recent and Plinian activity at Masaya can overlap, there is a tendency for samples of recent

activity to reach pre-eruptive temperatures which exceed 1100 °C23. CONCLUSION: PRE- AND SYN-ERUPTIVE CONDITIONS OF BASALTIC PLINIAN ERUPTIONS To draw a complete picture of the conditions

promoting basaltic Plinian volcanism, we compare our constraints on pre- and syn-eruptive conditions from Las Sierras-Masaya volcanic system with other basaltic Plinian eruptions

(Supplementary Table 4). The Las Sierras-Masaya, Etna and Tarawera Plinian eruptions show pre-eruptive temperatures <1100 °C. Low pre-eruptive temperature was likely a fundamental

condition to achieve high undercooling during magma ascent, rapid crystallisation, and brittle magma fragmentation15, evidenced by the high microlite crystal fractions observed in the Etna

122 BCE and Tarawera samples8,9. From our sensitivity analysis, we find that a critical condition for an explosive Plinian eruption is a pre-eruptive temperature <1100 °C and total

crystal content >30 vol.%, in agreement with Arzilli et al.’s15 predictions for the Etna 122 BCE Plinian eruption. Cooling, crystal-rich basaltic magma may accumulate in mush-bearing

systems such as Etna, which show a wide range in explosive activity62. The Las Sierras-Masaya and Etna 122 BCE Plinian eruptions shows similar volatile concentrations, suggesting that lower

initial volatile concentrations are sufficient to drive a basaltic Plinian eruption. However, a free fluid phase may have also driven the ascent dynamics of the FL and Etna 122 BCE

eruptions63. Similar to the FL eruption, crystallisation, closed system degassing and accumulation of volatiles within the uppermost storage region is expected to have occurred prior to the

Etna 122 BCE eruption63. However, CO2 fluxing and open system degassing also characterised ascent of the Etna 122 BCE magma from the deeper plumbing system63. Sub-Plinian to Plinian

eruptions show fast decompression rates (0.1–0.2 MPa s−1) and rapid ascent rates (>10 m s−1) within the conduit, indicating that rapid decompression and ascent provided little time for

gas-melt separation (Fig. 7c), even within short conduits such as those of Masaya and Sunset Crater12,64. Therefore, there may be limited time between the onset of potential precursory

activity and eruption. Although the last Plinian eruption of Las Sierras-Masaya volcanic system occurred ~1800 years ago, multiple episodes of Plinian activity suggests that under the right

conditions, there could be a sudden transformation to a highly explosive system. Changes in the pre-eruptive system could be identified through examining micro-textures of erupted products,

which may show crystalline enclaves indicative of processes occurring at depth such as formation of a mushy magma; monitoring changes in gas geochemistry to identify the arrival of fresh,

volatile-rich magma; or monitoring deformation and seismic signals to identify possible inflation or deflation following changes in the storage region. As large magma volumes can assemble at

Masaya within 50 years, the explosive transition could occur rapidly. As these conditions are shared with other basaltic Plinian eruptions, cool storage of crystal-rich basaltic magma at

shallow depth may be a crucial pre-eruptive condition driving highly explosive and hazardous eruptions at basaltic volcanoes. METHODS SAMPLING Samples of the FL eruption are 1–3 cm sized

scoriae, taken from units A/B to lower G (LG) 2. We use the same nomenclature for our samples, where FonD refers to unit D of the FL deposit. Samples are taken from proximal sites

12°03’18.4°N, 86°14’38.2°W and 12°01’36.9°N, 86°15’00.1°W (Fig. 1). Samples of the MTL eruption are from the same location as Bamber et al.16. MAJOR AND TRACE ELEMENT ANALYSIS Homogenised

glass chips were prepared for whole-rock compositional analysis at the Department of Physics and Geology, Universitá degli studi di Perugia, Italy, due to the presence of enclaves. Scoriae

of units A/B and E of the FL eruption were crushed to a fine powder and heated within a crucible until melted. This melt was then rapidly quenched in air on removal from the furnace. Glass

chips were homogenous and crystal-free. Glasses were mounted in epoxy resin and polished. For major and trace element analyses of glasses and minerals, scoriae (A/B to LG for the FL samples

and MTL_L2 for the MTL samples) were mounted as 300 µm thick wafers on glass slides and polished. Electron probe microanalysis (EPMA) was conducted using the JEOL JXA-8530F FEG electron

microprobe at the Photon Science Institute, University of Manchester, to determine major element chemistry of seven FL samples. BSE images of sample textures were also acquired using the

same instrument. BSE images were used to inspect MI prior to EPMA analysis and visibly crystallised inclusions were not analysed. We used a dual condition set-up to measure both crystals and

glasses whilst minimising potential devolatilisation of the glass. The operating conditions for analysis of crystals were: an acceleration voltage of 15 kV, a beam current of 7–10 nA and

beam diameter of 1–5 µm. For the analysis of glasses, the beam current was reduced to 2–7 nA and the beam diameter increased to 5–15 µm. Peak count times for all elements were 20–40 s,

except Mn, where peak count times were 20–50 s (Supplementary Table 5). Na and K were measured first to minimise loss via devolatilisation. Calibration standards were orthoclase (K, Si and

Al), wollastonite (Ca and Si), fayalite (Fe), periclase (Mg), jadeite (Na), tephroite (Mn), rutile (Ti), apatite (P) Cr2O3 (Cr), NiO (Ni) and pyrite and anhydrite (S). The basaltic glass

standard A9965 was used as a secondary standard. Data processing was accomplished using the Probe for EPMA software (v.12.9.5). Analytical precision (1σ as calculated by Probe for EPMA)

ranges between <1 and 3% for major elements with concentrations greater than 1%. For major element concentrations <1% analytical precision ranges between 12 and 14%. For element

concentrations <0.1%, analytical precision ranges between 41 and 52%. Whole-rock composition of FL samples was determined by measuring four glass chips using EPMA. The reported whole-rock

compositions are an average of a number of points (_n_ = 30–86) measured from each glass chip. Full compositional data are provided in Supplementary Data 1. Trace element concentrations

were determined using Laser Ablation Inductively Coupled Plasma Mass Spectrometry (LA-ICP-MS) on four samples at the Department of Physics and Geology, Università degli Studi di Perugia,

Italy. A Teledyne Photon Machine G2 laser ablation system was used, equipped with a Two-Volume ANU HelEx 2 cell and coupled with a Thermo Fisher Scientific quadrupole-based iCAP Qc ICP-MS

spectrometer. We used a circular laser beam, and a beam diameter of either 10–15 µm or 25 µm depending on MI size or microlite density of the glass. The external calibration material was

NIST-SRM610. The internal standard was Si. Precision and accuracy were monitored through repeated analyses of the USGS BCR2-G reference material66, analysed as an unknown together with FL

and MTL samples over three different analytical sessions. At 25 microns, precision and accuracy are always better than 10% with the only exceptions of Nb and Er, characterized by an accuracy

of ~13.6 and ~12%, respectively. At 10–15 microns, precision is always better than 12% with the only exception of Cs (~16%) and Tm (~20%). The accuracy is always better than ~12% with the

only exception of Yb (~14%). Complete trace element data are provided in Supplementary Data 3. MICRO-TEXTURAL ANALYSIS BSE images of MTL samples were also acquired using a FEI Quanta 650

FEG-SEM at the Department of Earth and Environmental Sciences, University of Manchester, applying an acceleration voltage of 15 kV and beam current of 10 nA. We used ImageJ software67 for

quantitative textural analysis. Crystal fractions were calculated as follows: $$\phi =\frac{{A}_{x}}{A}$$ (1) where _A__x_ represents the area of each crystal phase (_x_) and _A_, the

vesicle-free area of the sample68. The crystal fractions presented in this study are all calculated on a vesicle-free basis. The area fraction of enclaves was calculated using the total

area, including vesicles. We used _CSDCorrections_ (v.1.6)69 to produce CSDs of plagioclase crystals. We used _CSDSlice_70 to perform the stereological correction to provide the 3D crystal

habit. The CSD graph was used to calculate the crystallisation time (_τ_) of phenocrysts, microlites and crystals within enclaves, using the linear regression of the slope from the CSD

analysis and the following equation69. $${slope}=\,\frac{-\,1}{{Gt}\,\tau }$$ (2) where _Gt_ represents a range of experimentally derived growth rates. ANALYSIS OF VOLATILE CONCENTRATIONS

Volatile concentrations in FL and MTL samples were measured by Secondary Ion Mass Spectrometry (SIMS) at the NERC Ion Microprobe Facility, University of Edinburgh, UK, using a Cameca IMS

7f-Geo instrument. Individual scoria of the FL eruption (2) and MTL eruption (5) hosting phenocrysts were bulk mounted in epoxy resin. MI were only observed in plagioclase. Scoriae were then

polished to produce a flat surface suitable for analysis and to expose crystal-hosted MI at the surface. MI were visually inspected prior to analysis using BSE images, and crystallised

inclusions were not analysed. Some MI contained fluid bubbles, but most MI analysed are bubble-free. Samples were gold coated prior to analysis and placed under a high vacuum to outgas. Two

separate protocols were used to measure CO2 and H2O. The first protocol measured 12C+, using a primary O− ion beam with an accelerating voltage of 13 keV, a secondary accelerating voltage of

5 keV minus a 50 eV offset, a beam current of 5 nA and a beam diameter of 10 µm. The analysed area was pre-sputtered using the primary beam for 4 min prior to data acquisition. Count times

(seconds) are indicated in brackets. The masses 30Si+(2), 40Ca2+(2) and 24Mg2+ (5) were also measured alongside 12C+(10). The interference by 24Mg2+ on the 12C+ peak was resolved by using a

high mass resolution configuration of 2997. Masses 30Si+ and 40Ca2+ were acquired to allow comparison with major element data. For the carbon analysis, 15 cycles were acquired, but the first

7 cycles were discarded to minimise the effects from surface contamination. A second protocol was used to measure 1H+, using a primary O- ion beam with an accelerating voltage of 13 keV and

secondary accelerating voltage of 5 keV minus 75 eV offset. The analysed area was pre-sputtered for 1 minute prior to data acquisition. The masses 30Si+(2), 44Ca+(5) and 26Mg+(2) were also

measured alongside 1H+(5), over 10 cycles for each analysis with the beam centred in the pit made by the previous analysis. An extra peak at 0.7(1) was also measured to determine backgrounds

for each cycle. Background counts were monitored by analysing the CO2 and H2O-free standard N72, and mean background counts were subtracted from acquired data. Calibration curves for H2O

and CO2 (Supplementary Fig. 7) were constructed using the measured background-corrected counts per second of a suite of tholeiitic basaltic glass standards (M40, N72, M36, M21, M5, M10 and

M47)71, considering the SiO2 content. H2O and CO2 contents of unknowns were then calculated using their background-corrected counts per second and the measured SiO2 concentration of the MI,

as determined by EPMA. The EPMA analysis was performed after the SIMS analysis, to avoid surface contamination from carbon coating the samples. Precision (1σ) was calculated through repeat

measurements of the standard array and is ±3% for CO2 and ±2% for H2O. Accuracy (1σ) is ±13% for CO2 and ±5% for H2O across the standard array, calculated by comparing measured compositions

with published values. Images of MI and phenocryst hosts can be found in the Supplementary Information (Supplementary Figs. 10–21). Phenocryst, MI and bubble sizes can be found in

Supplementary Data 1. POST-ENTRAPMENT CRYSTALLISATION CORRECTION OF MELT INCLUSIONS Plagioclase-hosted MI were corrected for post-entrapment crystallisation (PEC) following Neave et al.72.

The plagioclase-melt equilibrium model of Putirka46, which accounts for the effect of H2O on equilibration, was used to calculate the plagioclase composition in equilibrium with the melt.

This equilibrium plagioclase composition was then added incrementally to the measured MI composition. For FL MI, the PEC correction is achieved when the corrected composition matches the

SiO2-FeO systematics predicted by a linear regression through glass data from this study and Costantini et al.10 on the FL eruption. For MTL MI, we matched the PEC-corrected compositions to

the MgO-Al2O3 systematics predicted by a linear regression through glass data of the MTL eruption16. The mean PEC correction for FL MI is 3% (range 0–8%). The mean PEC correction for MTL MI

is 3% (range 0–14%). This correction was then extended to volatile elements, which are assumed to be incompatible. The correction was also extended to trace element compositions, which were

corrected based on their compatibility in plagioclase, evaluated using experimental partition coefficients on plagioclase-melt partitioning for basaltic melts73,74,75. All trace elements

were corrected as incompatible except Sr, Cs and Pb. Raw and corrected MI compositions are presented in Supplementary Data 1 and 3. GEOCHEMICAL MODELLING Rhyolite-MELTS simulations were

performed using the average PEC-corrected FL MI composition, with fO2 fixed at NNO-0.65, consistent with the results of Bamber et al.16, to determine the pressure-temperature window for the

pre-eruptive storage condition. We performed simulations between fixed pressures of 100–10 MPa, with a pressure step of 10 MPa (Supplementary Fig. 8). The Moore et al.76 model was used as an

initial test to calculate the input H2O content (Supplementary Fig. 9). Isobaric simulations were performed under volatile-saturated conditions, using the MagmaSat solubility model

incorporated in Rhyolite-MELTS41. The temperature investigated ranged between 1150 and 950 °C, with a temperature step of 1 °C. For the liquid line of descent (LLD) model, isobaric

simulations were performed at 40 MPa and 200 MPa with initial H2O fixed at 2.3 wt.% and 5.2 wt.% respectively, to maintain H2O-saturated conditions. The initial CO2 was 0.12 wt.%.

Simulations were performed between temperatures of 1080–960 °C with a step of 5 °C. Isobars were produced using the MagmaSat software41. The average PEC-corrected MI composition was used to

produce H2O-CO2 isobars for the FL samples. Temperature was fixed at 1060 °C, consistent with mineral-melt thermometry results from this study and the maximum experimentally determined

pre-eruptive temperature11. The average MTL MI composition16 was used to construct the MTL H2O-CO2 isobars, instead with temperature fixed at 1100 °C, consistent with the maximum temperature

defined for the MTL pre-eruptive condition from petrological analysis16. Saturation pressure calculations were performed using MagmaSat, with an input temperature of 1060 °C for FL MI and

an input temperature of 1100 °C for MTL MI. All MagmaSat calculations were performed using the software VESIcal42. The fractional crystallisation model was produced using Petrograph

software34 and uses a _D_Zr partition coefficient of 0.001. _D_Sr is 1.714, _D_Ba is 0.327 and _D_La is 0.06173. The partition coefficients used were obtained from the experiments of

Aigner-Torres et al.73 on trace element partitioning between plagioclase and basaltic melt. These experiments were selected as the whole-rock composition used and the experimental conditions

investigated (1180 °C, NNO-0.75) can be considered representative of the FL pre-eruptive conditions. As plagioclase is the dominant crystallising phase, we assume that plagioclase partition

coefficients are representative of the bulk partition coefficient for the assemblage. As a linear trend can only be defined for the FL, a fractional crystallisation trend was modelled only

for this eruption, using the most primitive MI as the starting composition. MINERAL-MELT THERMOMETRY We used equation 23 of Putirka46 and an input water content of 2.27 wt.% to determine the

temperature of plagioclase crystallisation for the FL eruption. We also performed calculations using a lower water content of 1.59 wt.%, calculated as an average H2O content from MI of this

study, Goepfert and Gardner11 and Wehrmann et al.45. The standard error of the thermometer is ±43 °C. The _K_D (Ab-An) equilibrium criterion (0.28 ± 0.11 for T > 1050 °C and 0.1 ± 0.05

for T < 1050 °C) was used to test for plagioclase-melt equilibrium. We used an input pressure of 80 MPa, the maximum experimentally determined pre-eruptive storage pressure11. We used

equation 33 of Putirka46 and the same input pressure conditions to determine the temperature of pyroxene crystallisation. The standard error of the thermometer is ±45 °C. The _K_D (Fe-Mg)

equilibrium criterion (0.28 ± 0.08) was used to test for clinopyroxene-melt equilibrium. CONSTITUTIVE EQUATIONS OF THE 1D STEADY-STATE CONDUIT MODEL To simulate basaltic Plinian eruptions at

Las Sierras-Masaya Volcanic system, we use the 1D steady-state model of La Spina et al.25,26,27, which simulates magma ascent within a cylindrical conduit. The governing equations of the

model are provided in the Supplementary Methods. The model is presented in a general form which can be applied to a variety of volcanoes and eruptive styles25,26,27. To simulate the basaltic

Plinian eruptions of Las Sierras-Masaya volcanic system, a series of constitutive equations must be defined to describe the magma ascent dynamics for this particular volcanic system and

regime, including rheological, crystallisation, volatile solubility, outgassing, wall-friction, and fragmentation models. The rheology of the magmatic mixture has been modelled following La

Spina et al.64. The mixture viscosity is calculated as a function of the crystal-free, bubble-free melt viscosity, multiplied by two factors which account for crystals and bubbles. The

viscosity of the melt phase is calculated following Giordano et al.77. We used the average MTL MI composition (Supplementary Table 6) as the melt composition16. We account for the effect of

bubbles on the mixture viscosity following Llewellin et al.78, Mader et al.79 and La Spina et al.64. The effect on mixture viscosity is controlled by the capillary number (_Ca_), a

dimensionless number that describes the ratio between the relaxation time of the bubble in response to deformation and the timescale of deformation in the fluid79. At low capillary numbers,

bubbles tend to be spherical, whereas at high capillary numbers they are highly elongated. As the MTL samples dominantly show spherical bubbles, for our simulations we assume _Ca_ <<1.

The crystallisation model considers plagioclase, clinopyroxene and olivine. It is assumed that no crystals are removed during ascent. Similarly, to La Spina et al.25, the equilibrium mass

fraction \({x}_{{{{{{{\rm{c}}}}}}}_{j}}^{{{{{{\rm{l}}}}}},{{{{{\rm{eq}}}}}}}\)of each crystal phase _j_ is calculated using the polynomial function:

$${x}_{{{{{{{\rm{c}}}}}}}_{j}}^{{{{{{\rm{l}}}}}},{{{{{\rm{eq}}}}}}}\left({T}^{* },{x}_{{{{{{\rm{d}}}}}}}^{* }\right)=\,{\zeta }_{j,1}\cdot {\left({T}^{* }\right)}^{2}+\,{\zeta }_{j,2}\cdot

{\left({x}_{{{{{{\rm{d}}}}}}}^{* }\right)}^{2}+{\zeta }_{j,3}\cdot {T}^{* }\cdot {x}_{{{{{{\rm{d}}}}}}}^{* }+{\zeta }_{j,4}\cdot {T}^{* }+\,{\zeta }_{j,5}\cdot {x}_{{{{{{\rm{d}}}}}}}^{*

}+\,{\zeta }_{j,6}$$ (3) where \({T}^{* }\) is the temperature of the liquid in  °C and \({x}_{{{{{{\rm{d}}}}}}}^{* }\) is the water content dissolved in the liquid in weight percent. With

respect to La Spina et al.25, we neglect all terms related to pressure in the polynomial function, since we only consider water-saturated conditions. The coefficients _ζ__i,j_ are provided

in Supplementary Table 7. We also account for disequilibrium crystallisation, as it is likely that there is rapid syn-eruptive crystallisation15,16. The characteristic time of

crystallisation _τ_(c) is set to 10 s15. The solubility model for water and carbon dioxide has been described following Henry’s Law for gas mixtures15. The solubility parameters controlling

the equilibrium dissolved mass fraction of gas follow Arzilli et al.15. Outgassing has been described following La Spina et al.26. The degree of decoupling between gas and melt is controlled

by the relaxation parameter τ(f), which is defined to consider a permeability model based on Forchheimer’s law below the fragmentation level, and Stokes’ law elsewhere26. The parameters

controlling τ(f) follow La Spina et al.27. To simulate fragmentation during magma ascent, we adopt the strain-rate criterion of Papale53. Failure of the magma occurs when the Deborah number

(the ratio between the Maxwell relaxation timescale and the timescale of deformation) exceeds 0.0127. The pre-eruptive storage chamber for the Plinian eruptions of Masaya caldera is

estimated at a depth of ~2 km13,16. Geophysical investigation of the present Masaya plumbing system resolves a magma body with a depth-to-centre of ~3 km24. Due to the similarity of the

estimated reservoir depth for both Plinian and present-day activity, we use an initial pressure estimate of 42 MPa for our simulations16. Results can be found in Supplementary Data 4.

SENSITIVITY ANALYSIS The sensitivity analysis was performed using the open-source software DAKOTA toolkit (Design Analysis Kit for Optimization and Terascale Applications)80 developed at

Sandia National Laboratories. We selected a range of values for the input parameters which incorporate the pre-eruptive conditions determined for the eruptions of Las Sierras-Masaya volcanic

system: FL, SAT, MTL and MT-TIL3,10,11,13,16,21 and recent activity at Masaya caldera20,23 (Supplementary Table 8). The sensitivity analysis was performed using 10,000 simulations with

different combinations of the input parameters, varied assuming a uniform distribution within the aforementioned ranges. Results of the sensitivity analysis can be found in Supplementary

Data 5. CHAMBER VOLUME CALCULATION We estimate the magma chamber volume following the method of La Spina et al.64, coupling the steady-state conduit flow model with a quasi steady-state

model of the chamber. The following equation describes the variation in pressure through time as a function of the chamber volume and mass eruption rate81,82,83. $$\frac{\partial

{p}_{{ch}}}{\partial t}=\,\frac{4{EK}}{\rho {V}_{{ch}}(3K+4E)}\,\left({Q}_{{in}}-{Q}_{{out}}\right)$$ (4) The pressure of the magma chamber is indicated by _p__ch_, the volume by _V__ch_,

_E_ is the elastic modulus of the country rock, _K_ is the bulk modulus of the magma and _ρ_ is the average magma density. We use an elastic modulus of _E_ = 2 × 1010 Pa, a bulk modulus of

_K_ = 1010 Pa and an average magma density of _ρ_ = 2700 kg m−3 51. The mass flow rate is represented by _Q_, however, we assume negligible input into the magma chamber (_Q__in_), thereby

setting _Q__in_ to 0. The magma flux out of the chamber is represented by _Q__ou_t. We use the 1D steady-state conduit model to perform multiple simulations with the inlet pressure

decreasing accordingly with Eq. (4) and use the computed mass flow rates at each pressure step in our calculation. We stop the simulations when the pressure within the conduit is too low to

maintain conduit stability and collapse of the conduit walls is expected to occur, which we assume represents the end of the eruption. The collapse of the conduit wall is predicted following

the Mohr-Coulomb criteria64,84. Simulations were performed using both endmembers of the MTL eruption. As a first order approximation, we determine the eruption duration for three estimates

of chamber volume: 1 km3, 10 km3 and 100 km3. We assume that the most probable order of magnitude estimate of the magma chamber volume is that which shows the best agreement with the

estimated eruption duration and erupted volume, taken as a range which accounts for both endmembers. DATA AVAILABILITY All analytical data supporting the findings of this study are provided

within the article and Supplementary Data. This data and the data produced by numerical simulations are available at https://doi.org/10.6084/m9.figshare.20764858. CODE AVAILABILITY The magma

ascent model used in this study has been adapted from the MAMMA model available on GitHub: https://github.com/demichie/MAMMA. The version of the MAMMA model used in this work is available

at https://doi.org/10.6084/m9.figshare.20764858. REFERENCES * Williams, S. N. Plinian airfall deposits of basaltic composition. _Geology_ 11, 211–214 (1983). Article  CAS  Google Scholar  *

Costantini, L., Bonadonna, C., Houghton, B. & Wehrmann, H. New physical characterization of the Fontana Lapilli basaltic Plinian eruption, Nicaragua. _Bull. Volcanol._ 71, 337–355

(2009). Article  Google Scholar  * Pérez, W., Freundt, A., Kutterolf, S. & Schmincke, H. The Masaya Triple Layer: A 2100 year old basaltic multi-episodic Plinian eruption from the Masaya

Caldera Complex (Nicaragua). _J. Volcanol. Geotherm. Res._ 179, 191–205 (2009). Article  Google Scholar  * Coltelli, M., Del Carlo, P. & Vezzoli, L. Discovery of a Plinian basaltic

eruption of Roman age at Etna volcano, Italy. _Geology_ 26, 1095–1098 (1998). Article  CAS  Google Scholar  * Walker, G. P. L., Self, S. & Wilson, L. Tarawera 1886, New Zealand – a

basaltic Plinian fissure eruption. _J. Volcanol. Geotherm. Res._ 21, 61–78 (1984). Article  Google Scholar  * Alfano, F. et al. Subplinian monogenetic basaltic eruption of Sunset Crater,

Arizona, USA. _Geol. Soc. Am. Bull._ 131, 661–674 (2018). Article  Google Scholar  * Moussallam, Y. et al. Fast ascent rate during the 2017-2018 Plinian eruption of Ambae (Aoba) volcano: a

petrological investigation. _Contrib. Mineral. Petrol._ 174, 1–24 (2019). Article  CAS  Google Scholar  * Sable, J. E., Houghton, B. F., Del Carlo, P. & Coltelli, M. Changing conditions

of magma ascent and fragmentation during the Etna 122 BC basaltic Plinian eruption: Evidence from clast microtextures. _J. Volcanol. Geotherm. Res._ 158, 333–354 (2006). Article  CAS  Google

Scholar  * Sable, J., Houghton, B., Wilson, C. & Carey, R. _Eruption mechanisms during the climax of the Tarawera 1886 basaltic Plinian eruption inferred from microtextural

characteristics of the deposits_. (Stud. Volcanology: Leg. George Walk. Spec. Publ. IAVCEI 2, Geological Society, 2009). * Costantini, L., Houghton, B. & Bonadonna, C. Constraints on

eruption dynamics of basaltic explosive activity derived from chemical and microtextural study: The example of the Fontana Lapilli Plinian eruption, Nicaragua. _J. Volcanol. Geotherm. Res._

189, 207–224 (2010). Article  CAS  Google Scholar  * Goepfert, K. & Gardner, J. Influence of pre-eruptive storage conditions and volatile contents on explosive Plinian style eruptions of

basic magma. _Bull. Volcanol._ 72, 511–521 (2010). Article  Google Scholar  * Szramek, L. Mafic Plinian eruptions: is fast ascent required? _J. Geophys. Res. Solid Earth_ 121, 7119–7136

(2016). Article  Google Scholar  * Pérez, W., Freundt, A. & Kutterolf, S. The basaltic Plinian eruption of the ~6 ka San Antonio Tephra and formation of the Masaya caldera, Nicaragua.

_J. Volcanol. Geotherm. Res._ 401, 106975 (2020). Article  Google Scholar  * Moitra, P., Gonnermman, H. M., Houghton, B. F. & Tiwary, C. S. Fragmentation and Plinian eruption of

crystallizing basaltic magma. _Earth Planet. Sci. Lett._ 500, 97–104 (2018). Article  CAS  Google Scholar  * Arzilli, F. et al. Magma fragmentation in highly explosive basaltic eruptions

induced by rapid crystallization. _Nat. Geosci._ 12, 1023–1028 (2019a). Article  CAS  Google Scholar  * Bamber, E. C. et al. Pre-and syn-eruptive conditions of a basaltic Plinian eruption at

Masaya Volcano, Nicaragua. _J. Volcanol. Geotherm. Res._ 392, 106761 (2020). Article  CAS  Google Scholar  * Di Genova, D. et al. In situ observation of nanolite growth in volcanic melt: A

driving force for explosive eruptions. _Sci. Adv_ 6, eabb0413 (2020). Article  Google Scholar  * Girard, G. & van Wyk de Vries, B. The Managua Graben and Las Sierras-Masaya volcanic

complex (Nicaragua); pull-apart localization by an intrusive complex: results from analogue modelling. _J. Volcanol. Geotherm. Res._ 144, 37–57 (2005). Article  CAS  Google Scholar  *

Burton, M. R., Oppenheimer, C., Horrocks, L. A. & Francis, P. W. Remote sensing of CO2 and H2O emission rates from Masaya volcano, Nicaragua. _Geology_ 28, 915–918 (2000). Article  CAS 

Google Scholar  * Zurek, J., Moune, S., Williams-Jones, G., Vigouroux, N. & Gauthier, P.-J. Melt inclusion evidence for long term steady-state volcanism at Las Sierras-Masaya volcano,

Nicaragua. _J. Volcanol. Geotherm. Res._ 378, 16–28 (2019). Article  CAS  Google Scholar  * Pérez, W. Basaltic Plinian and violent Surtseyan eruptions from the Masaya Caldera Complex,

Nicaragua. _PhD thesis_, University of Kiel, 195 (2007). * Walker, J. A., Williams, S. N., Kalamarides, R. I. & Feigenson, M. D. Shallow open-system evolution of basaltic magma beneath a

subduction zone volcano: the Masaya Caldera Complex, Nicaragua. _J. Volcanol. Geotherm. Res._ 56, 379–400 (1993). Article  CAS  Google Scholar  * Hlinka, L., Longpré, M.-A., Pérez, W.,

Kutterolf, S. & Monteleone, B. Top-down control on eruptive style at Masaya volcano inferred from melt composition. _Earth Planet. Sci. Lett._ 572, 117138 (2021). Article  CAS  Google

Scholar  * Stephens, K. J. & Wauthier, C. Satellite geodesy captures offset magma supply associated with lava lake appearance at Masaya volcano, Nicaragua. _Geophys. Res. Lett._ 45,

2669–2678 (2018). Article  Google Scholar  * La Spina, G., Burton, M. & de’ Michieli Vitturi, M. Temperature evolution during magma ascent in basaltic effusive eruptions: a numerical

application to Stromboli volcano. _Earth Planet. Sci. Lett._ 426, 89–100 (2015). Article  Google Scholar  * La Spina, G., Polacci, M., Burton, M. & de’ Michieli Vitturi, M. Numerical

investigation of permeability models for low viscosity magmas: Application to the 2007 Stromboli effusive eruption. _Earth Planet. Sci. Lett._ 473, 279–290 (2017). Article  Google Scholar  *

La Spina, G. et al. Explosivity of basaltic lava fountains is controlled by magma rheology, ascent rate and outgassing. _Earth Planet. Sci. Lett._ 553, 116658 (2021). Article  Google

Scholar  * Freundt, A. & Schmincke, H.-U. Mixing of rhyolite, trachyte and basalt magma erupted from a vertically and laterally zoned reservoir, composite flow P1, Gran Canaria.

_Contrib. Mineral. Petrol._ 112, 1–19 (1992). Article  CAS  Google Scholar  * Hermes, O. D. & Cornell, W. C. Quenched crystal mush and associated magma compositions as indicated by

intercumulus glasses from Mt. Vesuvius, Italy. _J. Volcanol. Geotherm. Res._ 9, 133–149 (1981). Article  Google Scholar  * Ellis, B. S., Bachmann, O. & Wolff, J. A. Cumulate fragments in

silicic ignimbrites: the case of the Snake River Plain. _Geology_ 42, 431–434 (2014). Article  Google Scholar  * Pistone, M. et al. Gas-driven filter pressing in magmas: Insights into

in-situ melt segregation from crystal mushes. _Geology_ 43, 699–702 (2015). Article  CAS  Google Scholar  * Marsh, B. D. _Magma chambers_. (Encycl. of Volcan. San Diego, California, Academic

Press, 2015). * Gualda, G. A. R., Ghiorso, M. S., Lemons, R. V. & Carley, T. Rhyolite-MELTS: a modified calibration of MELTS optimized for silica-rich, fluid-bearing magmatic systems.

_J. Petrol._ 53, 875–890 (2012). Article  CAS  Google Scholar  * Petrelli, M., Poli, G., Perugini, D. & Peccerillo, A. PetroGraph: a new software to visualize, model, and present

geochemical data in igneous petrology. _Geochem. Geophys. Geosyst_. 6, Q07011 (2005). * Lissenberg, C. J., MacLeod, C. J. & Bennett, E. N. Consequences of a crystal mush-dominated magma

plumbing system: a mid-ocean ridge perspective. _Philos. Trans. R. Soc. A_, 377, 20180014 (2018). * Arzilli, F. et al. The unexpected explosive sub-Plinian eruption of Calbuco volcano (22–23

April 2015; southern Chile): triggering mechanism implications. _J. Volcanol. Geotherm. Res._ 378, 35–50 (2019b). Article  CAS  Google Scholar  * Moore, L. R. et al. Bubbles matter: an

assessment of the contribution of vapor bubbles to melt inclusion volatile budgets. _Am. Mineral._ 100, 806–823 (2015). Article  Google Scholar  * Wieser, P. E. et al. Reconstructing magma

storage depths for the 2018 Kīlauean eruption from melt inclusion CO2 contents: the importance of vapor bubbles. _Geochem. Geophys. Geosyst._ 22, e2020GC009364 (2021). Article  CAS  Google

Scholar  * Rasmussen, D. J., Plank, T. A., Wallace, P. J., Newcombe, M. E. & Lowenstern, J. B. Vapor-bubble growth in olivine-hosted melt inclusions. _Am. Mineral._ 105, 1898–1919

(2020). Article  Google Scholar  * Hartley, M. E., Maclennan, J., Edmonds, M. & Thordarson, T. Reconstructing the deep CO2 degassing behaviour of large basaltic fissure eruptions. _Earth

Planet. Sci. Lett._ 393, 120–131 (2014). Article  CAS  Google Scholar  * Ghiorso, M. S. & Gualda, G. A. R. An H2O-CO2 mixed fluid saturation model compatible with rhyolite-MELTS.

_Contrib. Mineral. Petrol._ 169, 1–30 (2015). Article  CAS  Google Scholar  * Iacovino, K., Matthews, S., Wieser, P. E., Moore, G. M. & Bégue, F. VESIcal Part I: an open-source

thermodynamic model engine for mixed volatile (H2O-CO2) solubility in silicate melts. _Earth Space Sci._ 8, e2020EA001584 (2021). Article  Google Scholar  * Lloyd, A. S., Plank, T.,

Ruprecht, P., Hauri, E. H. & Rose, W. Volatile loss from melt inclusions in pyroclasts of differing sizes. _Contrib. Mineral. Petrol._ 165, 129–153 (2013). Article  CAS  Google Scholar 

* Qin, Z., Lu, F. & Anderson, A. T. Diffusive re-equilibration of melt and fluid inclusions. _Am. Mineral._ 77, 565–576 (1992). CAS  Google Scholar  * Wehrmann, H., Hoernle, K.,

Portnyagin, M., Wiedenbeck, M. & Heydolph, K. Volcanic CO2 output at the Central American subduction zone inferred from melt inclusions in olivine crystals from mafic tephras. _Geochem.

Geophys. Geosyst_. 12, Q06003 (2011). * Putirka, K. D. Thermometers and barometers for volcanic systems. _Rev. Mineral. Geochem._ 69, 61–120 (2008). Article  CAS  Google Scholar  * Zellmer,

G. F. Gaining acuity on crystal terminology in volcanic rocks. _Bull. Volcanol_. 83, 78 (2021). * Gutiérrez, F. & Parada, M. A. Numerical modeling of time-dependent fluid dynamics and

differentiation of a shallow basaltic magma chamber. _J. Petrol._ 51, 731–762 (2010). Article  Google Scholar  * Cooper, K. M. Time scales and temperatures of crystal storage in magma

reservoirs: implications for magma reservoir dynamics. _Philos. Trans. R. Soc. A_, 377, 20180009 (2019). * Albert, H. et al. Magma interactions, crystal mush formation, timescales and unrest

during caldera collapse and lateral eruption at ocean island basaltic volcanoes (Piton de la Fournaise, La Réunion). _Earth Planet. Sci. Lett._ 515, 187–199 (2019). Article  CAS  Google

Scholar  * Tait, S., Jaupart, C. & Vergniolle, S. Pressure, gas content and eruption periodicity of a shallow, crystallising magma chamber. _Earth Planet. Sci. Lett._ 92, 107–123 (1989).

Article  CAS  Google Scholar  * Tramontano, S., Gualda, G. A. R. & Ghiorso, M. S. Internal triggering of volcanic eruptions: tracking overpressure regimes for giant magma bodies. _Earth

Planet. Sci. Lett._ 472, 142–151 (2017). Article  CAS  Google Scholar  * Papale, P. Strain-induced magma fragmentation in explosive eruptions. _Nature_ 397, 425–428 (1999). Article  CAS 

Google Scholar  * Shea, T. Bubble nucleation in magmas: a dominantly heterogeneous process? _J. Volcanol. Geotherm. Res._ 343, 155–170 (2017). Article  CAS  Google Scholar  * Cimarelli, C.,

Di Traglia, F. & Taddeucci, J. Basaltic scoria textures from a zoned conduit as precursors to violent Strombolian activity. _Geology_ 38, 439–442 (2011). Article  Google Scholar  *

Taddeucci, J., Pompilio, M. & Scarlato, P. Conduit processes during the July–August 2001 explosive activity of Mt. Etna (Italy): inferences from glass chemistry and crystal size

distribution of ash particles. _J. Volcanol. Geotherm. Res._ 137, 33–54 (2004). Article  CAS  Google Scholar  * Caravantes González, G. et al. Structures controlling volcanic activity within

Masaya caldera, Nicaragua. _Volcanica_ 2, 25–44 (2019). Article  Google Scholar  * Allison, C., Roggensack, K. & Clarke, A. B. Highly explosive basaltic eruptions driven by CO2

exsolution. _Nat. Commun._ 12, 217 (2021). Article  CAS  Google Scholar  * Allard, P. A CO2-rich gas trigger of explosive paroxysms at Stromboli basaltic volcano, Italy. _J. Volcanol.

Geotherm. Res._ 189, 363–374 (2010). Article  CAS  Google Scholar  * Freundt, A. & Tait, S. R. The entrainment of high-viscosity magma into low-viscosity magma in eruption conduits.

_Bull. Volcanol._ 48, 325–339 (1986). Article  Google Scholar  * Blake, S. & Campbell, I. H. The dynamics of magma-mixing during flow in volcanic conduits. _Contrib. Mineral. Petrol._

94, 72–81 (1986). Article  Google Scholar  * Ubide, T. & Kamber, B. S. Volcanic crystals as time capsules of eruption history. _Nat. Commun_ 9, 326 (2018). Article  Google Scholar  *

Moretti, R. et al. Degassing vs eruptive styles at Mt. Etna volcano (Sicily, Italy). Part I: Volatile stocking, gas fluxing, and the shift from low-energy to highly explosive basaltic

eruptions. _Chem. Geol._ 482, 1–17 (2018). Article  CAS  Google Scholar  * La Spina, G. et al. Conduit dynamics of highly explosive basaltic eruptions: the 1085 CE Sunset Crater sub-Plinian

events. _J. Volcanol. Geotherm. Res._ 387, 106658 (2019). Article  Google Scholar  * Jarosewich, E., Nelen, J. A. & Norberg, J. A. Reference samples for electron microprobe analysis.

_Geostandards Newsletter_ 4, 43–47 (1980). Article  Google Scholar  * Rocholl, A. Major and trace element composition and homogeneity of microbeam reference material: Basalt glass USGS

BCR‐2G. _Geostan. Geoanaly. Res_. 22, 33–45 (1998). Article  CAS  Google Scholar  * Abramoff, M. D., Magalhães, P. J. & Ram, S. J. Image processing with ImageJ. _Biophoton. Int._ 11,

36–42 (2004). Google Scholar  * Hammer, J. E., Cashman, K. V., Hoblitt, R. P. & Newman, S. Degassing and microlite crystallization during pre-climatic events of the 1991 eruption of Mt.

Pinatubo, Philippines. _Bull. Volcanol._ 60, 875–890 (1999). Article  Google Scholar  * Higgins, M. D. Measurements of crystal size distributions. _Am. Mineral._ 85, 1105–1116 (2000).

Article  CAS  Google Scholar  * Morgan, D. & Jerram, D. On estimating crystal shape for crystal size distribution analysis. _J. Volcanol. Geotherm. Res._ 154, 1–7 (2006). Article  CAS 

Google Scholar  * Shishkina, T. A., Botcharnikov, R. E., Holtz, F., Almeev, R. R. & Portnyagin, M. V. Solubility of H2O- and CO2-bearing fluids in theoliitic basalts at pressures up to

500 MPa. _Chem. Geol._ 277, 115–125 (2010). Article  CAS  Google Scholar  * Neave, D. A., Hartley, M. E., Maclennan, J., Edmonds, M. & Thordarson, T. Volatile and light lithophile

elements in high-anorthite plagioclase-hosted melt inclusions from Iceland. _Geochim. Cosmochim. Acta_. 205, 100–118 (2017). Article  CAS  Google Scholar  * Aigner-Torres, M., Blundy, J.,

Ulmer, P. & Pettke, T. Laser Ablation ICPMS study of trace element partitioning between plagioclase and basaltic melts: an experimental approach. _Contrib. Mineral. Petrol._ 153, 647–667

(2007). Article  CAS  Google Scholar  * Bindemann, I. N., Davis, A. M. & Drake, M. J. Ion microprobe study of plagioclase-basalt partition experiments at natural concentration levels of

trace elements. _Geochim. Cosmochim. Acta_. 62, 1175–1193 (1998). Article  Google Scholar  * Laubier, M., Grove, T. L. & Langmuir, C. H. Trace element mineral/melt partitioning for

basaltic and basaltic andesitic melts: an experimental and laser ICP-MS study with application to the oxidation state of mantle source regions. _Earth Planet. Sci. Lett._ 392, 265–278

(2014). Article  CAS  Google Scholar  * Moore, G., Vennemann, T. & Carmichael, I. S. E. An empirical model for solubility of H2O in magmas to 3 kilobars. _Am. Mineral._ 83, 36–42 (1998).

Article  CAS  Google Scholar  * Giordano, D., Russell, J. K. & Dingwell, D. B. Viscosity of magmatic liquids: a model. _Earth Planet. Sci. Lett._ 271, 123–134 (2008). Article  CAS 

Google Scholar  * Llewellin, E. W., Mader, H. M. & Wilson, S. D. R. The constitutive equation and flow dynamics of bubbly magmas. _Geophys. Res. Lett_. 29, 2170 (2002). * Mader, H. M.,

Llewellin, E. W. & Mueller, S. P. The rheology of two-phase magmas: a review and analysis. _J. Volcanol. Geotherm. Res._ 257, 135–158 (2013). Article  CAS  Google Scholar  * Adams, B. M.

et al. Multilevel parallel object-oriented framework for design optimization, parameter estimation, uncertainty quantification, and sensitivity analysis version 6.6. _User’s Manual_.

Technical report, SAND2014-4633 (Sandia National Laboratories, 2017). * Woods, A. W. & Koyaguchi, T. Transitions between explosive and effusive eruptions of silicic magmas. _Nature_ 370,

641–644 (1994). Article  CAS  Google Scholar  * Melnik, O. Dynamics of two-phase conduit flow of high-viscosity gas-saturated magma: large variations of sustained explosive eruption

intensity. _Bull. Volcanol._ 62, 153–170 (2000). Article  Google Scholar  * Melnik, O. & Sparks, R. S. J. Controls on conduit magma flow dynamics during lava dome building eruptions. _J.

Geophys. Res._ 110, B02209 (2005). Google Scholar  * Aravena, Á., de’ Michieli Vitturi, M., Cioni, R. & Neri, A. Stability of volcanic conduits during explosive eruptions. _J. Volcanol.

Geotherm. Res._ 339, 52–62 (2017). Article  CAS  Google Scholar  * Ryan, W. B. F. et al. Global Multi-Resolution Topography synthesis. _Geochem. Geophys. Geosyst_. 10, Q03014 (2009).

Download references ACKNOWLEDGEMENTS The work contained in this publication contains work conducted during a PhD study supported by the Natural Environment Research Council (NERC) EAO

Doctoral Training Partnership and is fully funded by NERC whose support is gratefully acknowledged (Grant Number NE/L002469/1). We acknowledge the RCUK NERC-NSFGeo DisEqm project

(NE/N018575/1), and support from INGV Pisa from a CASE studentship. F.A. thanks the KINETICVOLC project funded by the European Union, NextGenerationEU, Ministry of University and Research,

Promotion and Development Fund DM737/2021. We thank Instituto Nicaragüense de Estudios Territoriales (INETER), J.A. Saballos and D. Chavarría for their extensive and valuable assistance

during fieldwork in Nicaragua. Access to the Edinburgh Ion Microprobe Facility was supported by NERC grant IMF686/0519 to MRB, and we thank C. Talavera and J. Craven for their assistance

with the SIMS analyses. We thank the University of Perugia for use of their LA-ICP-MS facility, and D. Perugini and D. Morgavi for insightful discussions, laboratory access and assistance

with sample preparation. We thank D.A. Neave for assistance with post-entrapment crystallisation corrections. Sensitivity analyses were performed using the ARCHER2 National Supercomputing

Service. AUTHOR INFORMATION AUTHORS AND AFFILIATIONS * Department of Earth and Environmental Sciences, The University of Manchester, Oxford Road, Manchester, M13 9PL, UK Emily C. Bamber, 

Giuseppe La Spina, Fabio Arzilli, Margherita Polacci, Margaret E. Hartley, Jonathan Fellowes & Mike Burton * Istituto Nazionale di Geofisica e Vulcanologia - Osservatorio Etneo, Sezione

di Catania, Piazza Roma, 2, 95125, Catania, Italy Giuseppe La Spina * School of Science and Technology, Geology Division, University of Camerino, Camerino, Italy Fabio Arzilli * Department

of Geology, University at Buffalo, Buffalo, NY, 14260, USA Mattia de’ Michieli Vitturi * Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Pisa, Via Cesare Battisti, 53, 56125,

Pisa, Italy Mattia de’ Michieli Vitturi * Department of Physics and Geology, Università degli Studi di Perugia, Piazza dell’Universitá, 06123, Perugia, Italy Maurizio Petrelli Authors *

Emily C. Bamber View author publications You can also search for this author inPubMed Google Scholar * Giuseppe La Spina View author publications You can also search for this author inPubMed

 Google Scholar * Fabio Arzilli View author publications You can also search for this author inPubMed Google Scholar * Mattia de’ Michieli Vitturi View author publications You can also

search for this author inPubMed Google Scholar * Margherita Polacci View author publications You can also search for this author inPubMed Google Scholar * Margaret E. Hartley View author

publications You can also search for this author inPubMed Google Scholar * Maurizio Petrelli View author publications You can also search for this author inPubMed Google Scholar * Jonathan

Fellowes View author publications You can also search for this author inPubMed Google Scholar * Mike Burton View author publications You can also search for this author inPubMed Google

Scholar CONTRIBUTIONS E.C.B. conceived the research project with G.L.S., F.A., M.Po., M.B., M.E.H. and M.M.V. E.C.B., M.Po and F.A. collected the Masaya samples. E.C.B. performed the

numerical simulations with G.L.S. and M.M.V. E.C.B. performed the chemical analyses with F.A., M.E.H., and J.F. M.Pe performed the LA-ICP-MS analysis with E.C.B. All authors contributed to

production of the manuscript. CORRESPONDING AUTHOR Correspondence to Emily C. Bamber. ETHICS DECLARATIONS COMPETING INTERESTS The authors declare no competing interests. PEER REVIEW PEER

REVIEW INFORMATION _Communications Earth & Environment_ thanks Penny Wieser, Ery Hughes and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Primary Handling Editors: Emma Liu, Joe Aslin. ADDITIONAL INFORMATION PUBLISHER’S NOTE Springer Nature remains neutral with regard to jurisdictional claims in published maps and

institutional affiliations. SUPPLEMENTARY INFORMATION DESCRIPTION OF ADDITIONAL SUPPLEMENTARY FILES SUPPLEMENTARY INFORMATION SUPPLEMENTARY DATA 1 SUPPLEMENTARY DATA 2 SUPPLEMENTARY DATA 3

SUPPLEMENTARY DATA 4 SUPPLEMENTARY DATA 5 RIGHTS AND PERMISSIONS OPEN ACCESS This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use,

sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative

Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated

otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds

the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. Reprints and

permissions ABOUT THIS ARTICLE CITE THIS ARTICLE Bamber, E.C., La Spina, G., Arzilli, F. _et al._ Basaltic Plinian eruptions at Las Sierras-Masaya volcano driven by cool storage of

crystal-rich magmas. _Commun Earth Environ_ 3, 253 (2022). https://doi.org/10.1038/s43247-022-00585-5 Download citation * Received: 11 October 2021 * Accepted: 12 October 2022 * Published:

29 October 2022 * DOI: https://doi.org/10.1038/s43247-022-00585-5 SHARE THIS ARTICLE Anyone you share the following link with will be able to read this content: Get shareable link Sorry, a

shareable link is not currently available for this article. Copy to clipboard Provided by the Springer Nature SharedIt content-sharing initiative