Earth and Planetary Science Letters Explosivity of basaltic lava fountains is controlled by magma rheology, ascent rate and outgassing

The dichotomy between explosive volcanic eruptions, which produce pyroclasts, and effusive eruptions, which produce lava, is deﬁned by the presence or absence of fragmentation during magma ascent. For lava fountains the distinction is unclear, since the liquid phase in the rising magma may remain continuous to the vent, fragment in the fountain, then re-weld on deposition to feed rheomorphic lava ﬂows. Here we use a numerical model to constrain the controls on basaltic eruption style, using Kilauea and Etna as case studies. Based on our results, we propose that lava fountaining is a distinct style, separate from effusive and explosive eruption styles, that is produced when magma ascends rapidly and fragments above the vent, rather than within the conduit. Sensitivity analyses of Kilauea and Etna case studies show that high lava fountains ( > 50 m high) occur when the Reynolds number of the bubbly magma is greater than ∼ 0.1, the bulk viscosity is less than 10 6 Pas, and the gas is well-coupled to the melt. Explosive eruptions (Plinian and sub-Plinian) are predicted over a wide region of parameter space for higher viscosity basalts, typical of Etna, but over a much narrower region of parameter space for lower viscosity basalts, typical of Kilauea. Numerical results show also that the magma that feeds high lava fountains ascends more rapidly than the magma that feeds explosive eruptions, owing to its lower viscosity. For the Kilauea case study, waning ascent velocity is predicted to produce a progressive evolution from high to weak fountaining, to ultimate effusion; whereas for the Etna case study, small changes in parameter values lead to transitions to and from explosive activity, suggesting that eruption transitions may occur with little warning.


Introduction
The style of a volcanic eruption strongly influences the resultant hazard (Blong, 2013;Cassidy et al., 2018). Eruption style varies over a broad spectrum, but the most fundamental distinction is between effusive and explosive eruptions (Gonnermann and Manga, 2013). Effusive eruptions produce lava flows and lava domes and usually pose only a local hazard (Blong, 2013). Explosive eruptions produce pyroclastic material, which is ejected from the vent as a jet or plume (Walker, 1973), and may pose a local, regional, or global hazard. For silicic eruptions there is a clear dichotomy between effusive and explosive activity (Woods and Koyaguchi, 1994); even for hybrid eruptions, which exhibit simultaneous effusive and explosive activity, the component styles are clearly distinct and spatially localized (Castro and Dingwell, 2009;Schipper et al., 2013). At basaltic volcanoes, the distinction may be less clear.
Basaltic Plinian and sub-Plinian eruptions are unambiguously explosive, producing abundant fine ash ( 60 μm) that feeds an eruptive column up to tens of kilometres in height (Coltelli et al., 1998;Houghton et al., 2004;Gonnermann and Manga, 2013). Basaltic eruptions producing only lava flows are unambiguously effusive. However, basaltic lava fountains, whilst generally described as explosive (Mueller et al., 2019;Polacci et al., 2019), form a spectrum from ash-dominated to fluid-dominated ( Fig. 1) and so defy simple classification. Ash-dominated fountains are composed mainly of solid clasts from hundreds of microns up to a few millimetres in size, produced by brittle failure of magma, and share many characteristics with sub-Plinian eruptions (Parfitt, 1998;Houghton and Gonnermann, 2008). Fluid-dominated fountains are composed mainly of fluidal clasts from a millimetre scale (Pele's tears, Pele's hair) up to centimetre or decimetre scale (spatter and bombs) (Gonnermann and Manga, 2013;Cannata et al., 2019), produced by ductile fragmentation in the fountain. The larger clasts may coalesce on deposition to feed rheomorphic lava flows. At the fluid extreme, a dome fountain (Takahashi and Griggs, 1987) may produce almost no clastic material, and so must be considered effusive.
The degree to which a lava fountain is explosive or effusive is more than just a question of terminology. The nature and extent of the hazard a lava fountain poses depends strongly on the abundance, size, and other physical characteristics of the clasts. Fine ash generated from the fountain may loft into a buoyant plume that poses a downwind respiratory or aviation hazard, whereas fluidal clasts that coalesce on deposition may feed a lava flow that poses a hazard to downslope communities and infrastructure. Furthermore, the different clast characteristics suggest that different physical processes are responsible for their formation; it is unlikely that the same fragmentation process can generate such different pyroclastic materials as spindle bombs, scoria, Pele's hair/tears, and fine ash. These important distinctions are obfuscated by the current terminology.
Here we investigate our hypotheses by integrating all three fragmentation criteria (brittle, inertial, fluid-dynamic) into a 1D numerical conduit model for magma ascent, focussing our attention on two case studies: Kilauea and Etna. We vary key input parameters over a range of values that are consistent with each case study to discover general trends and thresholds. We then analyse how some initial conditions affect the style of activity. Finally, we provide some constraints on the style of activity as function of magma viscosity, ascent rates and outgassing.

1D steady-state model for magma ascent
Here we summarize the 1D steady-state model for magma ascent in a cylindrical conduit described by La Spina et al. (2015Spina et al. ( , 2016Spina et al. ( , 2017a, which explicitly accounts for the complex and nonlinear coupling of temperature changes, viscosity evolution, nonideal gas behaviour, outgassing and both disequilibrium crystallisation and exsolution. In this model, the flow of the magmatic two-phase multicomponent mixture is treated as a continuum, and each of the two phases and each component is characterized by the volume fraction (α k ), mass density (ρ k ), mass fraction (x k ), velocity (u k ), specific internal energy (e k ), specific entropy (s k ), pressure ( P k ) and temperature (T k ). Below the fragmentation level, magma is a mixture of a liquid phase (denoted by the index l), composed of melt, crystals and various dissolved gas components, and a gas phase (denoted by the index g), comprising bubbles of various exsolved gas components. Above the fragmentation level, the mixture consists of a dispersed particle phase (still denoted by index l) and a continuous gas phase composed of various exsolved gases. Since the gas mixture is always referred to as the exsolved gas phase, we can use, without ambiguity, the subscripts g 1 and g 2 to refer to the gas species considered here: H 2 O and CO 2 . The subscripts d 1 and d 2 are used to refer to those gas species when they are dissolved in the melt. The subscript m is used for the melt, and c for the total crystal content. Within each phase, we assume that all components have the same pressure, temperature and velocity (i.e. P m = P c = P d i = P l and P g i = P g , and analogously for the temperature and velocity). Furthermore, the saturation constraints of α l + α g = 1 and x l + x g = 1 hold all along the conduit. Using the notations described above, we can define the mixture parameters as follows: e = x l e l + x g e g ; s = x l s l + x g s g . (1) Relative motion between gas and liquid (or particle) phases is permitted, and pressures may differ between the gas and the liquid (or particle) phases. The temperature of the two phases is assumed to be in equilibrium. This assumption is supported by Oppenheimer et al. (2018) which shows that the temperature difference between the gas phase and the melt phase is negligible for bubble radii smaller than 10 cm, which is typical of Strombolian activity and lava fountains at Etna (Polacci et al., 2006).
The conservation equations for the mixture mass and momentum are: where g is the gravitational acceleration, r is the conduit radius, which is assumed to be constant, and φ f is a non-dimensional variable, which is 0 below the fragmentation level and 1 otherwise. Furthermore, f D l and f D g are the Darcy-Weisbach friction factors, which are functions of Reynolds number and conduit wall roughness as per the empirically-derived Moody diagram (Weisbach, 1845;Brown, 2003). The mixture energy equation has been derived following the formulation in Zein et al. (2010) and in Rodio and Abgrall (2015), including, also, heat loss due to gravitational force: The variation in the liquid volume fraction (and thus also in the gas volume fraction) along the conduit is defined by the following balance equation: Here, the relaxation parameter τ (p) (m 2 s −1 ) defines the disequilibrium between the gas and liquid pressures. This parameter is expected to be a function of the bulk viscosity, both volumetric fractions, and mixture density in the bubbly magma, and is assumed to be a very small constant (<< 1) in the gas/particle zone (La Spina et al., 2017b). In the bubbly magma, higher bulk viscosity leads to higher gas overpressure with respect to the surrounding melt. For a crystallized andesitic magma, La Spina et al. (2017b) showed that the gas overpressure is negligible throughout most of the conduit, except for a few meters below the fragmentation level, where the overpressure can reach up to 10 MPa. For a basalt, however, the viscosity is several orders of magnitude lower than a crystallized andesitic magma, which results in an overpressure of several orders of magnitude lower than that found in La Spina et al. (2017b). Since in this work we focus on basaltic magmas, it is reasonable to assume pressure equilibrium between the two phases, obtained by setting τ (p) = 10 −5 m 2 s −1 .
The following balance equations determine the exsolution of each gas component and the corresponding dissolved contents: where x md d i is the mass fraction of the dissolved gas phase i with respect to the liquid crystal-free phase, while x md,eq d i is the value at equilibrium. The exsolution rate (and thus how close the dissolved gas mass fraction is to the equilibrium value) is controlled by the relaxation parameter τ (e) (s). This characteristic time is likely to be controlled by pressure, temperature, dissolved water content and bulk composition of the melt. Since experimental constraints on this characteristic time are lacking, following La Spina et al. (2016), we consider (as a first order approximation) disequilibrium exsolution assuming a constant characteristic time τ (e) = 1 s.
The variations of the volume fractions of each crystal component are controlled by the following governing equation: Here β j is the volume fraction of the crystal component j while β eq j is the same physical quantity at equilibrium. The disequilibrium crystallization process is controlled by the characteristic time Finally, the following differential equation governs the relative motion of the gas with respect to the liquid phase: The degree of decoupling between the gas and the liquid phase is controlled by a relaxation parameter τ ( f ) (kg −1 m 3 s). When this relaxation parameter is very small, the relative velocity approaches zero, whereas large values of τ ( f ) result in strong decoupling. The velocity relaxation parameter τ ( f ) adopted here is defined according to the outgassing model described in the Supplementary Materials.

Constitutive equations
The model presented above is expressed in a general formulation, and can be used to describe magma ascent dynamics for a wide range of volcanoes (La Spina et al., 2015Spina et al., , 2016Spina et al., , 2019. In order to simulate the magma ascent dynamics at Kilauea and Etna we adopted specific constitutive equations (such as rheological, solubility, crystallization, wall friction, and outgassing models), which are reported in the Supplementary Materials. Extending previous works (Aravena et al., 2017;Arzilli et al., 2019;La Spina et al., 2019), here we implemented three different fragmentation criteria: the brittle (strain-rate) criterion (Papale, 1999), the inertial criterion (Namiki and Manga, 2008), and the fluid-dynamic breakup criterion (Jones et al., 2019).
For silicic magmas, magmatic failure is thought to occur when the Deborah number (De), which is the ratio between the Maxwell relaxation timescale (λ r ) and the timescale of deformation (λ d ), is greater than 0.01 (Papale, 1999): where μ is the magma viscosity, γ is the strain rate and G ∞ is the unrelaxed shear modulus at infinite frequency of the magma. This fragmentation mechanism is considered to be responsible for the formation of ash during Plinian and sub-Plinian eruptions, and is referred to as brittle fragmentation. The low viscosity of basaltic magmas, however, makes the formation of high strain rates necessary to trigger fragmentation difficult to achieve. Namiki and Manga (2008) proposed a new mechanism for basaltic magma fragmentation, based on the inertia of the magma generated by bubble expansion. They performed several decompression experiments in a shock tube apparatus with an analogue material covering the viscosity range of basaltic magmas (from 1 to 3 × 10 3 Pa s). With these experiments they found that fragmentation occurs when the expansion Reynolds number (Re exp ) is in the order of unity (Namiki and Manga, 2008): where ρ is the density of magma, α l is the magma volume fraction, v exp is the expansion velocity of the magmatic mixture, and L is the characteristic length of expansion. The characteristic length is calculated as the distance from where the exsolution starts to the position of interest in the conduit. The expansion velocity is also calculated following Namiki and Manga (2008): where P h is the pressure where the exsolution starts, and P l is the pressure at the position of interest in the conduit. Pele's hairs and tears have been suggested to be examples of textures related to this type of fragmentation (Namiki and Manga, 2008;Gonnermann, 2015), which is also considered to occur within the uppermost hundreds of meters of the conduit. Recently, Jones et al. (2019) presented a new criterion for fluid dynamic break-up of thinning filaments of low-viscosity magmas. They introduced a modified Deborah number (De * ), which is defined as the ratio of the dominant fluid instability timescale (λ inst ) and the extensional timescale (λ ext ). Their experiments indicate that fluid dynamic break-up occurs when where λ cap and λ vis are capillary and viscous instability timescales, respectively. The extensional timescale is the reciprocal of the extension rate of the filament γ ext , which can be computed aṡ where d the thickness of the filament. The dominant timescale of fluid instability is defined as the greater of the capillary instability timescale λ cap = ρ l d 3 /σ 1/2 and the viscous instability timescale λ visc = μ l d/σ . Jones et al. (2019) suggest that this mechanism may occur during magma ascent within the conduit or even above it, producing Pele's hairs and tears.

Initial conditions for the numerical simulations
We conduct simulations based on case study basaltic volcanoes from contrasting tectonic settings: Kilauea is a typical basaltic hot spot volcano erupting hot (∼1150 • C), crystal-poor magma; Etna is a typical arc basalt volcano erupting cooler (∼1100 • C), crystal-rich magma.
We perform two suites of numerical simulations. First, we perform runs for each of the two case study volcanoes using a set of reference parameter values (Table 1) to investigate the behaviour under each fragmentation criterion in turn. For numerical simulations of the Kilauea test case, we considered a cylindrical conduit of 2 km depth (Ferguson et al., 2016) and 20 m radius. The magma chamber pressure and temperature are, respectively, 55 MPa and 1150 • C (Garcia et al., 2000;Ferguson et al., 2016), while total water and carbon dioxide contents are, respectively, 0.7 wt.% and 400 ppm (Sides et al., 2014;Ferguson et al., 2016). The initial phenocryst content has been set to 0. For numerical simulations of the Etna test case we considered a cylindrical conduit of 9 km depth (Métrich et al., 2004;Corsaro et al., 2007;La Spina et al., 2016) and 25 m radius. The magma chamber pressure and temperature are, respectively, 250 MPa and 1090 • C (La Spina et al., 2016), while total water and carbon dioxide contents are, respectively, 3.4 wt.% and 4000 ppm (Métrich et al., 2004). The initial crystal content has been set to 0. Melt compositions used in the numerical simulations for Kilauea and Etna test cases are illustrated in Table 2. For Kilauea, the melt composition is obtained averaging the melt inclusion compositions for the 2008 and 2010 transient explosive eruptions at Kilauea presented in Sides et al. (2014). For Etna the composition is obtained by averaging the melt inclusion compositions for the 2001 Etna eruption (LV2) presented in Métrich et al. (2004).
We then perform a sensitivity analysis that sweep over ranges of parameter values appropriate for each case study (Table 1) to investigate the sensitivity of predicted eruptive behaviour to each parameter. For the Kilauea test case, we considered the following range of variables: 55-65 MPa for magma chamber pressure; 1060-1160 • C for inlet temperature; 1-30 m for conduit radius; 0.5-6 wt.% for total water content; 200-10000 ppm of total carbon dioxide, 0-20 vol.% of crystals. For the Etna test case we considered the following range of variables: 230-270 MPa for magma chamber pressure; 1050-1150 • C for inlet temperature; 1-30 m for conduit radius; 2-5 wt.% for total water content; 2000-20000 ppm of total carbon dioxide, 0-20 vol.% of crystals.

Numerical simulations of lava fountaining at Kilauea and Etna
Simulations assuming only the brittle criterion produce numerical results that satisfies the boundary conditions at the vent for both case studies (Fig. 2), but the fragmentation threshold (Deborah number >0.01) is not met within the conduit. This scenario would conventionally be considered effusive; however, the simula-tions reach a choked flow condition at the vent (i.e. ascent velocity equal to the speed of sound of the liquid/gas mixture) so that the pressure at the vent is many times higher than atmospheric pressure (Fig. 2a) -a condition generally associated with explosive eruption. This choked condition arises because the melt velocity is high throughout the conduit, reaching (at the vent of the conduit) ∼60 m s −1 for Kilauea, and ∼75 m s −1 for Etna (Fig. 2c). These high ascent rates result in very short magma ascent times, i.e.
∼3 min for Kilauea and ∼40 min for Etna. Furthermore, both temperature profiles along the conduit show a continuous cooling from depth up to the surface (Fig. 2b) (Ferguson et al., 2016). For the Kilauea test case, numerical results indicate a stronger coupling between gas and melt compared to the Etna test case. Indeed, the gas-slip velocity (i.e. the difference in velocity between gas and melt) for Kilauea reaches a maximum of ∼3 m s −1 at the vent, whereas for Etna we find ∼200 m s −1 (Fig. 2d). Notwithstanding the higher initial volatile contents of typical Etnean magmas Métrich et al., 2004;Ferguson et al., 2016 that results in a greater volume fractions of bubbles (Fig. 2e), the melt velocity within the conduit for the Etna test case is mostly lower than for the Kilauea test case, mainly due to stronger friction with the wall of the conduit caused by the higher viscosities of typical Etnean magmas (Giordano and Dingwell, 2003;Vona et al., 2011;Kolzenburg et al., 2018). Indeed, numerical results show a magmatic viscosity for the Kilauea test case of about 10 2 Pa s, whereas for Etna it reaches ∼10 5 Pa s close to the vent (Fig. 2f), consistent with typical viscosities of their magmas (Moune et al., 2007;Giordano and Dingwell, 2003;Andronico et al., 2009;Vona et al., 2011;Kolzenburg et al., 2018). Finally, strain-rates calculated for both test cases show small values for most of the ascent of the conduit (between 10 −5 and 10 −2 s −1 ), with a sudden increase in the last few hundreds of meters (up ∼10 s −1 ), but still not enough to trigger brittle fragmentation (Fig. 2g,h).
Simulations assuming the inertial criterion or fluid-dynamic breakup criterion produce similar outcomes: for both canonical case studies, there is no numerical solution that satisfies the boundary conditions at the vent (either atmospheric pressure or choked flow condition). This is because, for both criteria and both case studies, the fragmentation threshold is exceeded as soon as the gas volume fraction in the system rises to >0.01, and remains above this threshold value for the rest of the magma's ascent of the conduit. This implies that the gas would escape from deep in the conduit, leaving the melt behind and not resulting in an eruption.
We infer, therefore, that neither the inertial nor fluid-dynamic criteria can apply during ascent of the conduit, otherwise gas would escape from deep in the conduit, suppressing any eruption. Furthermore, for our reference cases the conduit model indicates that no brittle fragmentation criterion is likely to be met within the conduit, hence these are not explosive eruptions in the conventional sense (i.e. sub-Plinian or Plinian). However, magma arrives at the vent with a significant overpressure and with high velocity, hence the solutions do not represent conventional effusive eruptions either.
We therefore suggest that numerical solutions in which magma does not fragment in the conduit, but has a high exit velocity, represent a different style of activity with respect to the canonical explosive and effusive eruptions, and this is lava fountaining. This interpretation is consistent with the height calculated for the erupted jets in the reference case studies under a simple ballistic assumption (Head and Wilson, 1987). We calculate a height of ∼180 m for Kilauea and ∼290 m for Etna, both comparable with observed heights of lava fountains for these volcanoes (Andronico et al., 2008(Andronico et al., , 2009Mueller et al., 2019).
Whilst the simulations imply that fragmentation does not occur within the conduit, that does not preclude fragmentation once the magma is ejected from the vent. Fig. 2h shows that the Deborah number increases very rapidly towards the vent for both reference case studies, approaching the critical value for the Etna case. When the magma exits the vent, it is both overpressured (Fig. 2a) and highly vesicular (Fig. 2e); consequently, the sudden decompression to atmospheric pressure will drive very rapid bubble expansion. For cases where the magma is already near the critical Deborah number at the vent, we hypothesise that this expansion may be enough to trigger brittle fragmentation generating fine ash. We further hypothesise that brittle fragmentation is most likely to occur at the margins of the jet, where lava can expand into the surrounding atmosphere unimpeded. Both hypotheses are consistent with observations of lava fountains at Etna, which often comprise a coherent central jet that is mantled with fine ash (Fig. 1b).
For cases where the magma is far from meeting the brittle criterion at the vent, fragmentation by the inertial or fluid-dynamic breakup criteria may still occur above the vent. Using Eq. (6) in Jones et al. (2019) we can estimate the conditions under which the fluid-dynamic breakup criterion may be met for a thinning, stretching filament of fluid lava in the jet. Assuming a fountain containing spherical clasts of 5 cm in diameter (Jones et al., 2019), we estimate that, for the Kilauea case, fluid-dynamic breakup should occur once filaments thin to ∼100 μm, producing small, fluidal pyroclasts consistent with typical Pele's hairs dimensions (Moune et al., 2007). For the Etna case, fluid dynamic break-up should occur at ∼2 μm, too small to be consistent with Pele's hairs or tears. Neither the brittle nor fluid-dynamic breakup criterion is able to explain the formation of larger fluidal spatter clasts. The inertial criterion, which is based on expansion at the bubble scale, is also unlikely to be responsible for this larger-scale disruption. This leads us to conclude that an additional fragmentation mechanism is required to explain the formation of spatter.

Sensitivity analysis of controls on eruption style
To determine how eruptive style changes with pre-eruptive system parameters, we perform a sensitivity analysis. This was done using the DAKOTA (Design Analysis Kit for Optimization and Terascale Applications) toolkit (Adams et al., 2017), an open-source software developed at Sandia National Laboratories that provides a flexible and extensible interface between analysis codes and iterative systems analysis methods, such as uncertainty quantification, sensitivity analysis, optimization and parameter estimation. For each test case, we performed 10,000 conduit model simulations (assuming the brittle criterion for fragmentation), for which reservoir pressure, conduit radius, initial temperature, initial total (dissolved + exsolved) H 2 O and CO 2 content, and initial crystal content at depth are uniformly varied within the ranges illustrated in Table 1.
Simulation results are categorised into four eruption styles. If the brittle criterion is reached in the conduit, the eruption is "explosive" (i.e. representative of sub-Plinian or Plinian activity). If the brittle criterion is not reached, the ballistic jet height is calculated: "high fountain" denotes jet height is >50 m; "weak fountain" denotes jet height is 0.1-50 m; "effusive" denotes jet height is <0.1 m. The boundary between high and weak fountains is taken from Sides et al. (2014), whilst the boundary with effusive simulations is taken arbitrarily at 0.1 m. These boundaries, however, are not intended to be considered as hard borders between eruptive styles, but as a qualitative way to distinguish between different styles and establish general trends.
The sensitivity analysis shows that, among the input parameters investigated here, the initial temperature, crystal content, and H 2 O content are the dominant controls on the style of activity for both Kilauea and Etna (Fig. 3a-c, Fig. 4a-c). High lava fountains are favoured at higher initial temperature and lower initial crystal content. This is because higher temperature and lower crystal content result in lower magma viscosity (Fig. 3d, Fig. 4d), allowing faster ascent from depth and higher melt exit velocities (Fig. 3e, Fig. 4e). As a result, mass eruption rates for high fountains show higher values compared to the other styles, including explosive eruptions (Fig. 3f, Fig. 4f). Conversely, at lower initial temperature and higher initial crystal content, viscosity is higher, and weak fountains and effusive eruptions become increasingly common. At lower initial temperature and higher initial crystal content we also generate explosive solutions, but only when initial H 2 O content is also high (Fig. 3a-c, Fig. 4a-c). This is because high H 2 O content drives rapid bubble expansion; when combined with higher viscosity melt, this leads to higher Deborah number, promoting brittle fragmentation (Arzilli et al., 2019). Fig. 3. Numerical results of the sensitivity analysis for Kilauea test case. Frequency of effusive (blue dots), weak fountain (green circles), high fountain (yellow diamonds), and explosive (red triangles) solutions for Kilauea test case resulting from the 10,000 runs performed in the sensitivity analysis. Specifically, we plot the number of solutions as function of: (a) the inlet temperature; (b) the initial crystal content; (c) the H 2 O content; (d) the viscosity at the fragmentation for the explosive simulations, otherwise at the vent; (e) the melt exit velocity; (f) the mass eruption rate; (g) the cooling rate; (h) the decompression rate. Since inlet temperature, initial crystal content and H 2 O content are inputs of the sensitivity analysis, results are normalized as a percentage such that the sum of frequencies is always 100 for (a-c); for all the other quantities, the sum of the frequencies is variable, and is expressed as raw number (i.e. out of 10,000 runs).
Numerical results show also that cooling rates for high fountains are higher than those of the other styles of activity (Fig. 3g,  Fig. 4g). Cooling rates within the conduit for high fountains range between 10 −3 and 1 • C s −1 , whereas for weak fountains they are between 10 −4 and 10 −2 • C s −1 , for explosive eruptions they are between 10 −5 and 10 −2 • C s −1 , and for effusive eruptions they are between 10 −8 and 10 −3 • C s −1 . The higher cooling rates for high fountains are again a consequence of the high ascent rates associated with the style of activity, because the lower ascent time inhibits syn-eruptive crystallisation and thus reduces the release of latent heat in the system (La Spina et al., 2015). This also explains the different decompression rates found in our sensitivity studies (Fig. 3h, Fig. 4h). Numerical results show decompression rates between 10 −2 and 1 MPa s −1 for high fountains, 10 −3 and 10 −1 MPa s −1 for weak fountains, 10 −4 and 10 −1 MPa s −1 for explosive eruptions, and 10 −6 and 10 −2 MPa s −1 for effusive eruptions.
Within the parameter space investigated here for both test cases (Table 1), only a small fraction of the solutions (∼1%) for the Kilauea case study are explosive (Fig. 3); all the remaining solutions represent either effusive activity or lava fountaining. In contrast, ∼15% of the solutions for the Etna case study are explosive (Fig. 4). These findings are consistent with observations of the principal styles of volcanic activity at Kilauea and Etna: explosive eruptions that do not involve groundwater are very rare on Kilauea (Garcia et al., 2000;Sides et al., 2014;Ferguson et al., 2016;Mueller et al., 2019), but Etna has produced at least 24 sub-Plinian eruptions and one Plinian eruption during the Holocene (Coltelli et al., 1998;Houghton et al., 2004). Specifically, for Kilauea, explosive eruptions are predicted for temperatures <1080 • C, initial crystal content 10 vol.%, and H 2 O content >3 wt.%. These pre-eruptive conditions fall far from typical values estimated for Kilauea (i.e. a temperature of ∼1150 • C, a crystal content <5 vol.%, and a H 2 O content 1 wt.%) (Garcia et al., 2000;Sides et al., 2014;Ferguson et al., 2016), suggesting that a high lava fountain is very unlikely to transition into an explosive sub-Plinian or Plinian eruption. For Etna, explosive eruptions require an initial temperature 1100 • C and a H 2 O content >2 wt.%. These pre-eruptive conditions are consistent with those estimated numerically for the 122 BCE Plinian eruption at Etna (Arzilli et al., 2019).
Analysis of the numerical results shows that magma viscosity at the vent or fragmentation front is a powerful discriminant of eruption style (Fig. 3d, Fig. 4d). We plot the Reynolds number of the liquid phase against viscosity (Fig. 5a,c) to explore the role of inertia on eruption style. Parameters are taken at the fragmentation depth if the brittle criterion is met within the conduit, otherwise at the vent. The different eruption styles are grouped into clear fields, with limited overlap, showing that style is strongly influenced by Reynolds number. High fountains are expected at higher Reynolds number and lower magma viscosities (Fig. 5a,c). Specifically, within the parameter space investigated for both Kilauea and Etna, high lava fountaining is predicted for bulk viscosity <10 6 Pa s and Reynolds number 0.1 (Fig. 5a,c). We also plot the gas-slip velocity (i.e. gas velocity minus melt velocity) normalized to the melt ascent velocity, both taken at the fragmentation depth if the brittle criterion is met within the conduit, otherwise at the vent, against viscosity (Fig. 5b,d) to explore the role of outgassing on eruption style. Again, results show that style is strongly affected by the degree of decoupling between gas and melt (Fig. 5b,d). For both Kilauea and Etna, high fountains are produced when the gasslip velocity is less than 10 times melt velocity, and this condition is typically met when viscosity <10 5 Pa s. For higher melt viscos- Fig. 4. Numerical results of the sensitivity analysis for Etna test case. Frequency of effusive (blue dots), weak fountain (green circles), high fountain (yellow diamonds), and explosive (red triangles) solutions for Etna test case resulting from the 10,000 runs performed in the sensitivity analysis. Specifically, we plot the number of solutions as function of: (a) the inlet temperature; (b) the initial crystal content; (c) the H 2 O content; (d) the viscosity at the fragmentation for the explosive simulations, otherwise at the vent; (e) the melt exit velocity; (f) the mass eruption rate; (g) the cooling rate; (h) the decompression rate. Since inlet temperature, initial crystal content and H 2 O content are inputs of the sensitivity analysis, results are normalized as a percentage such that the sum of frequencies is always 100 for (a-c); for all the other quantities, the sum of the frequencies is variable, and is expressed as raw number (i.e. out of 10,000 runs). ity, numerical results show gas-slip velocity is 10-100 times the melt velocity (Fig. 5b,d) and eruptions display different styles of activity (explosive, effusive, or weak fountains), depending on the initial conditions of the magmatic system. Specifically, these results suggest, counterintuitively, that explosive eruptions do not always require strong coupling between gas and melt, and that it is possible to trigger brittle fragmentation even when the gas-slip velocity is 100 times the ascent velocity (Fig. 5b,d).
The region where fluid-dynamic breakup may occur above the vent is indicated in Fig. 5 for both case studies, following Jones et al. (2019). For magma viscosities typical of Kilauea (grey region in Fig. 5a,b), our model predicts that fluid-dynamic breakup will occur in most high fountains. This is consistent with Pele's hairs and tears being typical pyroclastic products of lava fountaining at Kilauea (Moune et al., 2007). For magmatic viscosities typical of Etna (Fig. 5c,d grey region), fluid-dynamic breakup is predicted for only a small percentage of fountains, consistent with Pele's hairs and tears being produced only rarely during lava fountaining at Etna. However, Pele's hairs and tears may still be produced if the magma viscosity is sufficiently low; this was the case of the Etna 2002-2003 eruption, in which Pele's hairs were erupted together with crystal-poor sideromelane clasts, whereas ash was produced together with crystal-rich tachylite clasts (Andronico et al., 2009).
Our steady-state simulations cannot directly inform on temporal changes in eruption style. Nonetheless, we can use the mapping between system parameters and eruption style encapsulated in Figs. 3-5 to infer how eruption style is likely to evolve in response to typical temporal trends in basaltic eruption parameters. For example, during the course of a high-fountaining eruption, the flux wanes with time, resulting in a decrease in Reynolds number. The Kilauea case study suggests that this causes a shift to weak fountaining or effusive activity (Fig. 5a), consistent with observations (Garcia et al., 2000), whereas it is unlikely that high fountaining would evolve into a highly explosive sub-Plinian or Plinian eruption. The Etna case study also shows that fountaining activities may evolve into effusive eruptions. However, if the magmatic system is in the right thermodynamic conditions, lava fountaining could evolve instead into a highly explosive eruption (Fig. 5c,d). This indicates that a modest shift in eruption rate or outgassing efficiency could cause a dramatic shift in eruption style. These shifts could be enhanced by feedback in the system. For example, at higher viscosities ( 10 6 Pa s), a decrease in the magma ascent rate as reservoir pressure decays gives more time for crystals to grow during ascent. A small increase in crystal content can lead to a dramatic increase in viscosity (Vona et al., 2011), which may trigger brittle fragmentation within the conduit (Arzilli et al., 2019).

Conclusions
In this work we investigate the parameters that control basaltic eruption styles using a numerical steady state model for magma ascent, and considering Kilauea and Etna as case studies. We implement three fragmentation criteria (brittle, inertial, and fluiddynamics breakup criterion) in order to understand which conditions lead to fragmentation of basaltic magmas.
Our results indicate that the inertial criterion and fluid-dynamic breakup criterion do not apply during ascent in the conduit because they predict fragmentation at very low gas volume fraction (<0.01), leading to gas escape at depth, leaving melt behind, and not resulting in an eruption. Our numerical results also highlight that, for some initial conditions, magma can arrive at the vent with a significant overpressure and with high velocity, but without fragmenting inside the conduit. We thus suggest that this type of activity represents a different style with respect to the canonical explosive and effusive eruptions, and we identify it with lava fountaining. We speculate that fragmentation may occur once the magma is ejected from the vent, by any of the three criteria, depending on conditions. For cases where magma is already near the critical Deborah number at the vent, we hypothesise that the additional expansion as the magma is ejected into the atmosphere may be enough to trigger brittle fragmentation generating fine ash. For cases where the magma is far from meeting the brittle criterion at the vent, fragmentation by the inertial or fluid-dynamic breakup criteria may still occur above the vent, producing small, fluidal pyroclasts consistent with Pele's hairs and tears. The inertial and fluid-dynamic breakup criteria, however, are unlikely to be responsible for larger-scale disruption, such as the formation of larger fluidal spatter clasts. This leads us to conclude that an additional fragmentation mechanism is likely to be responsible for the formation of spatter.
We also perform a sensitivity study in order to investigate how eruptive style changes with pre-eruptive system parameters. Our results on Kilauea and Etna case studies show that high lava fountains are favoured at higher temperatures and lower initial crystal content, hence at lower viscosity. Specifically, high lava fountains occur when the Reynolds number of ascent is greater than ∼0.1, the magma viscosity is less than 10 6 Pa s, and the difference between the gas and the melt velocity is lower than the melt velocity. Furthermore, owing to lower magma viscosities, high lava fountains show higher ascent rates and lower ascent times than explosive eruptions, resulting in higher cooling and decompression rates. Finally, we find that lower viscosity basalts, typical of Kilauea, rarely produce explosive eruptions, whereas higher viscosity basalts, typical of Etna, may produce explosive, effusive, or fountaining activity over a small region of parameter space. Thus our results suggest that, unless there is a significant cooling of the magmatic system and a strong increase in volatile content, it is unlikely that lava fountains at lower viscosity basaltic volcanoes (similar to Kilauea) would evolve into a highly explosive eruption, whereas at higher viscosity basaltic volcanoes (such as Etna) transitions to and from explosive activity occur readily and potentially with little warning.
The simulations and analyses presented here demonstrate that basaltic eruptions show complex dependence on a range of fragmentation mechanisms, and that these mechanisms are not completely understood. Better understanding of basaltic explosive eruptions, their fragmentation mechanisms, and the transitions between different eruptive styles, is essential for mitigating the risk associated with these usually quiet and gentle, but still potentially very dangerous, eruptions.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.