Stability of volcanic conduits during explosive eruptions

Accepted Manuscript Stability of volcanic conduits during explosive eruptions

Álvaro Aravena, Mattia de'Michieli Vitturi, Raffaello Cioni, Augusto Neri PII: DOI: Reference:

S0377-0273(17)30142-7 doi: 10.1016/j.jvolgeores.2017.05.003 VOLGEO 6088

To appear in:

Journal of Volcanology and Geothermal Research

Received date: Revised date: Accepted date:

16 March 2017 4 May 2017 4 May 2017

Please cite this article as: Álvaro Aravena, Mattia de'Michieli Vitturi, Raffaello Cioni, Augusto Neri , Stability of volcanic conduits during explosive eruptions, Journal of Volcanology and Geothermal Research (2017), doi: 10.1016/j.jvolgeores.2017.05.003

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Stability of volcanic conduits during explosive eruptions Aravena, Álvaro1; de'Michieli Vitturi, Mattia2; Cioni, Raffaello1; Neri, Augusto2. 1 Dipartimento di Scienze della Terra, Università di Firenze, Via Giorgio La Pira 4, I-50121,

T

Firenze, Italy.

IP

2 Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Pisa, Via della Faggiola 32, I-56126,

CR

Pisa, Italy.

US

Corresponding author: Á. Aravena.

AC

CE

PT

ED

M

AN

E-mail: [email protected]

ACCEPTED MANUSCRIPT Abstract Geological evidences of volcanic conduit widening are common in most pyroclastic deposits (e.g. presence of lithic fragments from different depths), suggesting a continuous modification of the conduit geometry during volcanic eruptions. However, the controlling factors of the

T

mechanisms driving conduit enlargement (e.g. erosion, local collapse) are still partially unclear,

IP

as well as the influence of conduit geometry on the eruptive dynamics. Although numerical

CR

models have been systematically employed to study volcanic conduits, their mechanical stability and the eruptive dynamics related to non-cylindrical conduits have been poorly addressed.

US

We present here a 1D steady-state model which includes the main processes experimented by

AN

ascending magmas (i.e. crystallization, rheological changes, fragmentation, drag forces, outgassing and degassing), and the application of two mechanical stability criteria (Mohr –

M

Coulomb and Mogi – Coulomb), in order to study the collapse conditions of volcanic conduits

ED

during a representative explosive rhyolitic eruption. It emerges that mechanical stability of volcanic conduits is mainly controlled by its radial dimension, and a minimum radius for

PT

reaching stable conditions can be computed, as a function of water content and inlet

CE

overpressure. Additionally, for a set of input parameters thought typical of explosive rhyolitic volcanism, we estimated a minimum magma flux for developing a mechanically stable conduit kg/s). Results are consistent with the unsteady character usually observed in

AC

(

sub-Plinian eruptions, opposite to mainly stationary Plinian eruptions, commonly characterized by higher magma discharge rates. We suggest that cylindrical conduits represent a mechanically stable configuration only for large radii. Because the instability conditions are not uniform along the conduit, the widening processes probably lead to conduit geometries with depth-varying width. Consequently, as our model is able to consider volcanic conduits with depth-dependent

ACCEPTED MANUSCRIPT radius, two plausible and previously untested geometries have been studied, evidencing major and complex modifications in some eruptive parameters (particularly, exit pressure and mass discharge rate), and suggesting that the geometry acquired by the conduit as it is widened influences the eruptive dynamics.

T

Keywords: Volcanology. Conduit stability. Rhyolitic magmas. Numerical modelling. Explosive

AC

CE

PT

ED

M

AN

US

CR

IP

volcanism.

ACCEPTED MANUSCRIPT 1. Introduction Volcanic conduit processes have proved to be important controlling factors of eruptive dynamics [Houghton et al., 2004; Kennedy et al., 2005; Pioli et al., 2009]. In explosive eruptions, the eruptive dynamics seem to be mainly controlled by the relation between magma supply and

T

discharge rate [Scandone and Malone, 1985; Bursik, 1993], which are functions of magma

IP

chamber conditions (e.g. overpressure, volatile content, temperature) and volcanic conduit

CR

geometry. Several studies have addressed the relation between eruptive dynamics and the nature of the magma source [Civetta et al., 1997; Gurioli et al., 2005; La Spina et al., 2015], whereas

US

the relative importance of volcanic conduit geometry on the eruptive dynamics remains poorly

AN

understood [Costa et al., 2007; de'Michieli Vitturi et al., 2008; Koyaguchi et al., 2010]. Geological evidences of changes in volcanic conduit geometry (i.e. erosive processes) are

M

common in the volcanic record, as revealed by the occurrence of lithic fragments in most

ED

pyroclastic deposits, reaching frequently mass concentrations higher than 20% [Calder et al., 2000; Maeno and Taniguchi, 2009]. Lithic fragments can proceed from different depths because

PT

of the occurrence of craterization and/or conduit enlargement. These erosion mechanisms are

CE

two different processes that often act together during volcanic eruptions, and the knowledge of subsurface stratigraphy of a given volcano can help to distinguish them [Barberi et al., 1989].

AC

However, the controlling factors of the erosive mechanisms (and thus, incorporation of lithic fragments during magma ascent) are still partially unclear, and quantitative studies for analyzing them are scarce. Macedonio et al. [1994] identified four main erosive mechanisms: (1) impact of pyroclasts, (2) fluid shear stress, (3) conduit wall collapse and (4) volcanic tremor; and indicated the zones where each one is probably dominant. In particular, in order to study conduit wall collapse, Macedonio et al. [1994] used a criterion based on the difference between conduit and

ACCEPTED MANUSCRIPT lithostatic pressures, establishing a maximum value for a stable conduit, as a function of the tensile strength of the country rocks. Two additional remarkable works have addressed some aspects of conduit erosion during effusive eruptions [Dragoni and Santini, 2007; Piombo et al., 2016], showing that shear traction would make the conduit section closer to a circular shape and

T

that conduit enlargement is needed to explain the typical evolution of effusion rate during

IP

basaltic eruptions. Nonetheless, quantifications of the stability of volcanic conduits have never

CR

been addressed. The main limitation for studying numerically this issue is our insufficient knowledge about conduit stability criteria and conduit geometry evolution during volcanic

US

eruptions.

AN

On the other hand, wall collapse conditions along a vertically extended cavity have been systematically analyzed in the geothermal research field and petroleum industry for the analysis

M

of borehole stability [Rahman et al., 2000; Aadnoy and Edland, 2001]. For studying rock failure

ED

phenomena, the most frequently used procedure is the application of the Mohr – Coulomb criterion. It considers the maximum and minimum principal stresses and some parameters of the

PT

country rocks, ignoring the intermediate principal stress, and therefore it is expected to be a

CE

conservative criterion [Al-Ajmi and Zimmerman, 2006]. Because numerous cases have been reported where the Mohr – Coulomb criterion predicted deficient results and because it has been

AC

suggested a significant effect of the intermediate principal stress on rocks behavior [Vernik and Zoback, 1992; Single et al., 1998], Al-Ajmi and Zimmerman [2006] developed the Mogi – Coulomb collapse criterion, whose formulation includes the intermediate principal stress, so predicting a larger resistance of rocks. Although these collapse criteria can be applied to the analysis of mechanical stability of volcanic conduits, they have never been used in volcanological studies.

ACCEPTED MANUSCRIPT In this study, we address numerically the mechanical stability of volcanic conduits. For that, we analyzed the collapse conditions of cylindrical conduits, using a modified version of the 1D steady-state model presented by de'Michieli Vitturi et al. [2011] and applying the collapse criteria described by Al-Ajmi and Zimmerman [2006]. Additionally, we test the effects of non-

T

cylindrical geometries in the eruptive dynamics, by analyzing conduits with depth-varying width.

IP

This study consists of five parts: (1) we present a general description of the 1D steady-state

CR

model used in this work; (2) we describe a procedure for applying the collapse criteria presented by Al-Ajmi and Zimmerman [2006] (Mohr – Coulomb and Mogi – Coulomb) and the set of input

US

and geometric parameters that we imposed in this study, in order to model representative

AN

explosive rhyolitic eruptions; (3) we present the results related to cylindrical conduit simulations, including a global sensitivity analysis of our results; (4) we describe simulations of volcanic

M

conduits with depth-dependent radius (here referred as non-cylindrical geometries), focusing on

ED

the influence of conduit geometry on the eruptive dynamics (e.g. exit pressure, mass discharge rate); and (5) we discuss our results, trying to provide some lights about the applicability of these

2. Methodology

CE

explosive eruptions.

PT

collapse criteria in volcanic conduits, and the significance of conduit stability in the dynamics of

AC

2.1 Steady-state conduit model In this work, we employed a modified version of the 1D multiphase steady-state model presented in de'Michieli Vitturi et al. [2011] (http://demichie.github.io/MAMMA). The system of conservation equations is derived from the theory of thermodynamically compatible systems [Romenski et al., 2010], considering the effects of the main processes that magmas experiment during ascent, such as crystallization, rheological changes, fragmentation, physical interaction

ACCEPTED MANUSCRIPT with the conduit walls, outgassing and degassing. The system is described as a mixture of two phases (

), each one characterized by a volumetric fraction ( ), density ( ), velocity ( )

and specific entropy ( ). Below the fragmentation level, phase 1 is a mixture of crystals, dissolved water and melt (continuous phase); while phase 2 is composed by the exsolved gas

T

bubbles (dispersed phase). On the other hand, above the fragmentation level, phase 1 is

IP

constituted by magma fragments (dispersed phase) and phase 2 is the exsolved gas mixture

CR

(continuous phase). Following Sparks [1978] and Starostin et al. [2005], magma fragmentation

US

is assumed to occur when the exsolved gas volumetric fraction reaches a critical value ( ).

AN

Each component of the system is characterized by an equation of state, and the specific energy is given as a function of density and entropy, obtained by a linearized form of the Mie-Grüneisen

M

equations. Additionally, the pressure ( ) and temperature ( ) of both phases are derived from

CE

PT

ED

the internal energy ( ):

AC

The system of equations is composed by the conservation laws of total mass (equation 3), momentum (equation 4), energy (equation 5) and mass of phase 1 (equation 6).

acceleration of gravity, of phase ,

is the mixture density,

is the mixture viscosity,

is the mixture temperature,

is the mixture velocity,

IP

is the vertical coordinate,

is the conduit radius,

is the

is the mass fraction

CR

where

T

ACCEPTED MANUSCRIPT

is the relaxation parameter which controls the gas

is the mass fraction of the dissolved gas,

same parameter,

is the volume fraction of crystals and

US

exsolution rate,

is the equilibrium value of the

AN

is the density of crystals.

where

ED

M

Additionally, the phase 1 volume fraction is governed by the following equation:

is the relaxation parameter which controls the pressure difference between both

PT

phases.

AC

CE

Furthermore, we include the following equation for the relative velocity between the phases:

where

is a drag factor, is the continuous phase index ( or , function of the relative position

of the fragmentation level), while

controls the sign of the viscous term contribution ( or

,

it is also function of the continuous phase index). Finally, we include the mass conservation laws of crystals (equation 9) and dissolved gas (equation 10).

ACCEPTED MANUSCRIPT

where

is the relaxation parameter of the crystallization process,

is the equilibrium

T

volume fraction of crystals and the subscript refers to the liquid component of the system.

IP

For the numerical solution of the model equations we employ a numerical shooting technique:

CR

for a given inlet pressure, the model adjusts the inlet flow rate until the appropriate boundary

US

condition (choked flow or atmospheric pressure) is reached. For the spatial integration of the equations a well-established proportional-integral step-size control technique is adopted, with the

AN

relaxation terms treated implicitly to guarantee the stability of the numerical scheme. 2.2 Constitutive equations

M

In order to apply the described 1D steady-state model to a representative case of explosive

ED

rhyolitic volcanism, we employed appropriate constitutive equations to describe the rheological

PT

evolution, crystallization, gas exsolution, outgassing and the equations of state of the components of the system.

CE

2.2.1 Magma rheology model

) is evaluated using the following

AC

According to Hess and Dingwell [1996], melt viscosity ( empirical relationship:

where

is the dissolved

concentration in weight percent,

resulting viscosity value is measured in

.

is expressed in Kelvin and the

ACCEPTED MANUSCRIPT On the other hand, below the fragmentation level, phase 1 viscosity is computed using the following equation:

is a factor which considers the effect of the crystals [Caricchi et al., 2007]. It is function ) and it was imposed following Costa [2005]:

,

and

are fitting parameters (supplementary Table S1).

AN

where ,

US

CR

IP

of the volumetric fraction of crystals (

T

where

M

Our model includes an additional term for controlling the effect of bubbles on rheology.

ED

However, despite it has been suggested a non-negligible effect of bubbles on rheology [Llewellin and Manga, 2005], since there are no systematic agreements about the dependence between

PT

bubbles and magma viscosity (at least, hardly applicable to general cases) [Mader et al., 2013;

CE

Truby et al., 2015] and its influence is expected to be sensibly lower than other processes such as crystallization and exsolution of volatiles [Manga and Loewenberg, 2001; Truby et al., 2015],

AC

here we have ignored its effect. 2.2.2 Solubility model Water is considered the only dissolved volatile component. We adopted the following solubility model:

ACCEPTED MANUSCRIPT where

is the pressure of the gas component,

is the solubility coefficient and

is the

solubility exponent. Based on Zhang [1999], we set the following values:

and

.

T

2.2.3 Crystallization model

is the maximum crystallinity,

is the initial volume fraction of crystals and

US

where

CR

IP

We adopted the crystallization model presented by de'Michieli Vitturi et al. [2010] (equation 16).

AN

is the minimum function (supplementary Table S1). 2.2.4 Outgassing model

M

Because a non-linear relationship between pressure gradient and gas flow rate has been

ED

recognized, the outgassing process is described using the Forchheimer’s law [Rust and Cashman,

where

AC

CE

forces (quadratic term):

PT

2004; Degruyter et al., 2012], which includes the influence of viscous (linear term) and inertial

is the velocity difference between both phases, the subscript

gas phase, while

and

refers to the exsolved

are the Darcian and inertial permeabilities, respectively, which are

computed using the following equations:

ACCEPTED MANUSCRIPT

is the throat-bubble size ratio, and

is the bubble density

are fitting parameters. In this work, we employed the following values:

,

,

and

T

number and

is the average bubble size,

.

IP

where

CR

2.2.5 Equations of state

US

In order to define the specific internal energy of the melt, crystals and dissolved water, we adopted a linearized version of the Mie-Grüneisen equations of state [Le Métayer et al., 2005; La

while

is the specific heat capacity at constant volume,

are the density and sound speed at a reference state,

PT

and

represents the formation energy,

is the pressure at a reference state. Here the subscript

water or crystals.

CE

where

ED

M

AN

Spina and de'Michieli Vitturi, 2012]:

is the adiabatic exponent, refers to the melt, dissolved

AC

For the exsolved water, we adopted the Van der Waals equations of state:

ACCEPTED MANUSCRIPT where

and

are the critical pressure and temperature of water, respectively.

The parameters that we used in the equations of state are shown in supplementary Table S2, and Table 4 exhibits a notation summary.

T

2.3 Conduit collapse criteria

IP

In order to evaluate the stability conditions of cylindrical conduits, we compared the minimum , predicted by Mohr - Coulomb and Mogi

CR

pressure needed to avoid conduit collapse (

- Coulomb collapse criteria, Table 1) [Al-Ajmi and Zimmerman, 2006] and the conduit pressure , computed with the 1D steady-state model). In particular, we used the pressure

US

profile (

)

AN

profile of phase 1, but in any case, the pressure relaxation parameter here adopted (

produces negligible differences between both pressure profiles (i.e. continuous and dispersed

M

phases). The Mohr - Coulomb and Mogi - Coulomb collapse criteria have been developed for

ED

studying boreholes stability, where an appropriate methodology for calculating the critical mud pressure in order to avoid borehole collapse is a crucial research topic [Gholami et al., 2014;

PT

Curry et al., 2016]. Indeed, it is required for performing drilling activities and lowering the

CE

operating costs [Al-Ajmi and Zimmerman, 2006]. These stability criteria have been widely applied in the petroleum industry during the last decade [Gholami et al., 2014] and, based on the

AC

intrinsic similarities between these systems and volcanic conduits, we assume that these collapse criteria can be applied to the analysis of mechanical stability of volcanic conduits, as shown here. Note that here we study exclusively the collapse conditions of volcanic conduits, whereas the fracture criteria presented by Al-Ajmi and Zimmerman [2006] and thus failure processes due overpressure conditions were not addressed in this work. The key parameters of the country rocks for applying the Mohr – Coulomb and Mogi – Coulomb collapse criteria are: cohesion ( ), angle of internal friction ( ) and the horizontal and vertical

ACCEPTED MANUSCRIPT stress gradients (

,

and

). Additionally, the presence of fluids in the country rock

pores (e.g. presence of aquifers) is expected to produce a destabilizing effect, whereas a strong difference between the two horizontal stress gradients tend to increase the country rock strength. In this study, we used representative values of natural conditions for the country rock mechanical

T

parameters (Table 2), and we analyzed three different measures for addressing the instability

IP

character of the conduit. It is important to highlight that the definition of the country rock

CR

mechanical parameters presents a strong uncertainty, mainly derived from the vertical variability

US

of lithologies and structures commonly observed in volcanic fields, and also the eventual depthdependence of other parameters, such as cohesion and friction angle. However, in any case, the

AN

expected cohesion variations because of the lithostatic pressure in the depth ranges here studied are much lower than the intrinsic cohesion variability derived from lithological variations and

ED

values of cohesion and friction angle.

M

geological structures and thus, following Al-Ajmi and Zimmerman [2006], we assumed constant

2.4 Input and geometric parameters

PT

The input parameters that we employed in the numerical simulations can be considered

CE

representative of typical conditions of explosive rhyolitic volcanism. We varied three input parameters: inlet pressure, water content and conduit dimensions (Table 3). We imposed a

AC

conduit length of 5 km, which implies a lithostatic pressure at the conduit bottom of 127.4 MPa. Since radius variations do not limit the applicability of the 1D steady-state model that we adopted, in order to test the effect of different conduit geometries, three cases were considered and compared: (1) cylindrical conduit (geometry C), (2) cylindrical conduit with a local enlargement along a vertical distance

(geometry NC1), and (3) conduit with two coaxial

cylindrical portions connected by a linear enlarging transitional zone of length

(geometry

ACCEPTED MANUSCRIPT NC2) (Figure 1). For geometries with depth-dependent diameter (i.e. NC1 and NC2), we fixed the middle level of the variable-radius zone ( , Figure 1) around the fragmentation level computed in the equivalent cylindrical case (i.e.

using the same

input parameters of inlet pressure and water content), for different values of

and

T

(Table 3). This assumption and the choice of these conduit geometries are based on the expected

IP

nature of the gravitational collapse process, because of the typical location of the instability

CR

zones, as discussed in the following sections. On the other hand, although the model is

US

potentially able to deal with conduits characterized by elliptical cross-sections, this would make more complex the application of the stability criteria here adopted, so we preferred to restrict this

AN

work to vertical conduits with circular cross sections. 3. Results

M

3.1 Cylindrical conduits

) of an unstable conduit, where the minimum

ED

Figure 2 shows a typical pressure profile (

) are estimated using both collapse

PT

pressures needed to avoid conduit collapse (

criteria. Three measures of the instability character of the conduit have been evaluated:

CE

(1) The integral along the conduit of

.

(3)

AC

(2) The length of the unstable domain of the conduit. .

The collapse conditions are preferably reached near the fragmentation level, in correspondence with a sharp gradient in the upward decrease of

; and Mohr – Coulomb criterion predicts

more unstable conduits than Mogi – Coulomb criterion, as expected. In the following sections, we only use the third instability measure (hereafter, instability index) because it provides useful

ACCEPTED MANUSCRIPT data even for stable conduits and, in any case, similar conclusions can be obtained with the other two instability measures. Figure 3 illustrates the relationships between some key variables of the eruption and the instability index when inlet pressure, water content and geometry are varied. Conduit stability is

T

highly conditioned by the resulting fragmentation depth (Figure 3a), which is explained because

IP

deep fragmentation produces low-pressure fields in domains characterized by high lithostatic

CR

pressures, leading to favorable conditions for conduit collapse. As shallow fragmentation levels are favored by wide conduits, a monotonic tendency can be observed between instability index

US

and conduit radius (supplementary Figure S1). Because of the direct relation between conduit

AN

radius and mass discharge rate, this fact is clearly reflected in the relation between the instability index and magma flux (Figure 3b). An interesting result is that a minimum mass discharge rate kg/s and

M

for a stable conduit can be estimated (

kg/s, depending on the stability

ED

criteria), for the set of constitutive equations here adopted (representative of rhyolitic explosive eruptions with a 5000 m-length conduit).

PT

For given values of inlet pressure and water content, the instability index is a decreasing function

CE

of conduit radius (supplementary Figure S1), and a minimum radius for a stable conduit can be inferred (critical radius, hereafter). Figure 4 shows the critical radius as a function of inlet

AC

pressure and water content, using both stability criteria. Despite the obtained values of critical radius are sensitive to mechanical parameters of the country rocks (supplementary Figure S2), there are significant differences between both collapse criteria, which confirm an important effect of the intermediate principal stress on the country rock strength. Critical radius results to be mainly controlled by the inlet pressure, while the water content only exhibits a weak influence.

ACCEPTED MANUSCRIPT Additionally, in order to test the influence of the conduit length, we developed some additional simulations using values between 3,000 m and 9,000 m. Stability conditions are favored by short conduits (supplementary Figures S3 and S4), which is partially explained by the deeper magma fragmentations obtained when long conduits are considered. Nonetheless, it does not imply a

T

significant control of the conduit length on the critical radius results, as it is shown in the

IP

following section.

CR

3.1.1 Sensitivity analysis of the critical radius

US

We developed a global sensitivity analysis of the critical radius based on a Monte Carlo sampling scheme. We assumed uniform probability distributions for three input parameters: (a) ), between -12.4 and 7.6 MPa with respect to the lithostatic value (these

AN

inlet overpressure (

(b) water content (

M

values have been set in order to be consistent with the pressure range employed in section 3.1); ), ranging from 4.5 to 6.5 wt. %; and (c) conduit length (

ED

and 6000 m.

), between 4000

PT

For each set of input parameters, we performed eight simulations with different conduit radii (random values between 35 m and 115 m), which allow us to estimate the critical radius ( ),

CE

using both stability criteria (

for Mohr – Coulomb,

for Mogi - Coulomb).

AC

As a sensitivity parameter, we calculated the first order Sobol sensitivity index ( ) [Sobol, 2001], which is defined following equation 25, and represents the fraction of the output variability ( ) only attributed to the variability assigned to the input parameter

.

ACCEPTED MANUSCRIPT On the other hand, the expression for calculating the total effect sensitivity index (

) is shown

in equation 26, and measures the contribution to the output variance ( ) of the input

,

T

including all interactions with the other parameters.

and

, for Mohr - Coulomb and

CR

critical radius is the inlet overpressure (

IP

Sensitivity analysis results (Figure 5) confirm that the most important factor controlling the

US

Mogi - Coulomb collapse criteria, respectively), while water content exhibits a weak influence, higher when the Mohr – Coulomb criterion is employed (

and

, for Mohr -

AN

Coulomb and Mogi - Coulomb collapse criteria, respectively). Conduit length has a minor effect

and

M

on the critical radius, more significant when the Mogi – Coulomb criterion is employed ( , for Mohr - Coulomb and Mogi - Coulomb collapse criteria, respectively).

ED

Because the sum of the first order sensitivity indexes is only slightly lower than 1 and the

PT

differences between both sensitivity indexes are relatively small, the output variance is mainly controlled by the first order contributions, and the effect of interactions between input parameters

CE

plays a secondary role.

AC

3.2 Conduits with depth-dependent radius A viable evolution of a volcanic conduit passes through the initial formation of a vertical dyke in which the magma ascending path rapidly focuses into more cylindrical zones of lower friction [Bruce and Huppert, 1989]. We can then assume that these zones are progressively widened up to reach the minimum radius for stable conditions, as a consequence of variations in both inlet pressure and water content. However, collapse conditions may vary along the conduit (i.e. they are preferably reached near and above the fragmentation level, Figure 2), and therefore,

ACCEPTED MANUSCRIPT preservation of a cylindrical geometry is very unlikely in natural conditions. In this section, we explore the eruptive dynamics of conduits with depth-dependent radius (here referred as noncylindrical conduits), in order to study their changes with plausible conduit geometries produced by wall instabilities. Additional simulations were performed, introducing small geometric

T

perturbations in the conduit (10% and 20% of relative difference between maximum and

IP

minimum radius), using the non-cylindrical geometries previously described (section 2.4,

CR

Figures 1 and 6).

US

3.2.1 Geometry NC1

Figure 6a-b exhibits the general results related to geometry NC1. The pressure profile along the

AN

conduit exhibits a low-pressure zone coincident with the conduit widening, as expected

M

(supplementary Figure S5). This geometry produces deeper fragmentation levels than the equivalent cylindrical case (differences up to ~200 m, Figure 6a), which is maximized for high and

(supplementary Figure S6). Additionally, this geometry generates

ED

values of

PT

higher magma discharge rates than the equivalent cylindrical case (mass discharge rate increment between ~3% and ~13%, Figure 6b). On the other hand, counterintuitively, exit pressures are up

CE

to ~14% higher than the values obtained with the equivalent cylindrical case, while the

AC

consequent increase in exit velocity (at choked conditions) is lower than 1% (Figure 6a-b). However, during the subsequent decompression to atmospheric conditions, an important acceleration of the mixture is expected, which would increase its velocity near the vent and would produce higher convective plumes [Carcano et al., 2013]. All the tendencies previously described are more pronounced for water-rich magmas. 3.2.2 Geometry NC2

ACCEPTED MANUSCRIPT Figure 6c-d shows the results related to geometry NC2. The pressure profile shows a slightly larger pressure drop at the fragmentation level and lower pressures above it, in relation to the equivalent cylindrical case (supplementary Figure S5). This geometry produces a more complex behavior:

T

(a) Wide conduits with high water contents: they produce deeper fragmentation levels than the

IP

equivalent cylindrical case (up to ~100 m of depth difference, Figure 6c) and slightly higher

CR

mass fluxes (up to 3% of the mass flux with respect to the equivalent cylindrical case, Figure 6d).

US

(b) Narrow conduits with low water contents: they produce shallower magma fragmentations

AN

than the equivalent cylindrical case (40 m or so of depth difference, Figure 6c) and lower mass discharge rates (about 90% of the mass flux of the equivalent cylindrical case, Figure 6d).

M

Both cases generate lower exit velocities and pressures than the equivalent cylindrical case, with

ED

more important differences for narrow conduits. These changes in the eruptive dynamics can have important effects on decreasing the buoyancy of eruptive plumes.

PT

4. Discussion

CE

Despite stability conditions of eruptive conduits are highly sensitive to mechanical parameters of the country rocks and the applied collapse criterion, for eruptions of high viscosity magmas, we

AC

suggest that cylindrical conduits are mechanically stable only for large radii, and therefore, the occurrence of widening processes is very likely. This is not a surprising result because it is consistent with the commonly observed presence of lithic fragments in pyroclastic deposits [Barberi et al., 1989; Varekamp, 1993; Calder et al., 2000; Maeno and Taniguchi, 2009], whose concentration is determined by the relative importance of the erosive rate with respect to the magma discharge rate through the conduit. Consequently, as conduit widening produces higher

ACCEPTED MANUSCRIPT mass discharge rates, it is not necessarily manifested in lithic-rich horizons in pyroclastic deposits or in an increase of the relative amount of lithic clasts. Conduit collapse conditions predicted with both stability criteria exhibit significant differences, indicating that the intermediate principal stress could have an important role in volcanic conduit

T

stability. For a given set of input parameters (thought typical of explosive rhyolitic volcanism)

IP

and using the Mogi – Coulomb criterion, we determined that rhyolitic magmas require conduit

CR

radii larger than 50 m for reaching stability conditions. Critical radius (i.e. minimum radius required to avoid conduit collapse) is mainly controlled by the inlet overpressure (Figures 4 and

US

5), while the water content has a minor influence. So, under the assumption that magma chamber

AN

conditions during a typical volcanic eruption follow a depressurizing trend, a continuous widening process is expected in the conduit. This process could explain the pervasive and

M

continuous presence of lithic fragments in most pyroclastic deposits, even with stationary

ED

properties and conditions of the magma source (e.g. water content, temperature, composition). Moreover, the expected conduit widening could increase the magma discharge rate during

PT

explosive events, favoring magma chamber evacuation and counterbalancing the effect of a

CE

depressurizing magma reservoir (indeed, for a given conduit radius and water content, a decrease

2005].

AC

in magma chamber overpressure is manifested in lower mass discharge rates) [Macedonio et al.,

Based on the presence of a minimum mass discharge rate required to support stable conduits (between

and

kg/s, function of the stability criterion), we suggest that the

unsteady character commonly observed in events with relatively low mass fluxes (e.g. subPlinian eruptions) can be related to unstable conditions of the volcanic conduit, manifested in episodic collapse events and intermittent discharge rates [Bursik, 1993; Cioni et al., 2003].

ACCEPTED MANUSCRIPT Conversely, the sustained and quasi-stationary Plinian eruptions, characterized by higher mass discharge rates [Cioni et al., 2000], would reflect a prolonged conduit stability, and the erosive processes could be preferentially controlled by fluid shear stress and the impact of pyroclasts [Macedonio et al., 1994]. These unsteady phenomena have been observed in some case studies,

T

such as the 1630 AD eruption of Furnas volcano [Cole et al., 1995], the Mangaone subgroup

IP

eruptions [Jurado-Chichay and Walker, 2001] and the 1886 Tarawera eruption [Carey et al.,

CR

2007], where the variability of the eruptive dynamics was interpreted as a consequence of variable processes of vent wall collapse. Indeed, the unsteadiness of volcanic eruptions and its

US

quantification have been recently studied, and these phenomena are observed in a wide range of

AN

viscosities [Dominguez et al., 2016], whose comprehension can be relevant for studying explosivity mechanisms. In particular, our results suggest that conduit stability could be a

M

dominant process in the unsteady character of Subplinian eruptions.

ED

On the other hand, because instability conditions are not uniform along the conduit, we suggest that, when collapse conditions are reached, they trigger localized widening processes. As the

PT

other erosive mechanisms also act preferentially around (fluid shear stress) and above the

CE

fragmentation level (volcanic tremor and impact of pyroclasts) [Macedonio et al., 1994], it seems unlikely to preserve cylindrical geometries in volcanic conduits. Consequently, we considered

AC

two additional conduit geometries (Figures 1 and 6). A geometry with a localized enlargement (NC1) produces a more abrupt pressure decrease near the fragmentation level and higher pressure values above it. This could inhibit the collapsing character at shallower levels of the conduit and, consequently, it seems a stable geometry above the wider domain. On the other hand, because a deepening tendency of the fragmentation level is observed (comparing with the equivalent cylindrical case) and considering the importance of the

ACCEPTED MANUSCRIPT fragmentation level in the stability conditions, a downward-propagation of the erosive mechanisms could be expected (i.e. extension of the wider zone), with relevant consequences on the eruptive dynamics, such as higher exit pressures and mass discharge rates. On the other hand, a conduit geometry characterized by a wider upper part (NC2) produces an

T

abrupt pressure decrease near the fragmentation level and lower pressure values above it, which

).

US

changes of

-increasing tendency with no major

CR

collapse mechanisms and the impact of pyroclasts (i.e.

IP

could lead to a continuous erosive process above the fragmentation level, both because of wall

The tested non-cylindrical geometries (i.e. conduits with depth-dependent width) do not

AN

necessarily tend to stable conditions, which is probably reflected in a continuous erosive process and a pervasive presence of lithic fragments in pyroclastic deposits, even with a steady magma

M

source (e.g. water content and inlet overpressure). On the other hand, considering that a uniform

ED

widening process of volcanic conduits (i.e. permanent cylindrical geometry) would produce shallower fragmentation levels and an increase in mass discharge rate, exit velocity and exit

PT

pressure (supplementary Figure S5), it is important to note that non-cylindrical geometries do not and

.

CE

exhibit an intermediate behavior between the cylindrical cases using

Consequently, the sequence of geometries that the conduit acquires as it is eroded would be

AC

manifested in different variation trends of the most important eruptive parameters, suggesting that conduit evolution could play a primary role for controlling syn-eruptive changes of the eruptive dynamics. Finally, despite the difficulty to predict the most probable erosive sequence in a specific event, some parameters, such as the presence of different lithologies in the country rocks, could represent key information to understand the geometric evolution of a volcanic conduit and its

ACCEPTED MANUSCRIPT implications on the eruptive dynamics. Nonetheless, because the scarce available information about country rocks properties in volcanic fields (e.g. cohesion, friction angle), there are not enough data for developing further analysis and corroborating the expected relation between the dynamics of explosive eruptions and the properties of country rocks.

T

5. Concluding remarks

IP

(a) In rhyolitic explosive eruptions, cylindrical conduits represent a mechanically stable

CR

geometry only for large radii.

(b) For a given set of magma properties, it is possible to calculate a minimum mass discharge

US

rate for a mechanically stable conduit, suggesting that the unsteady character usually

AN

observed in low-mass flux events (e.g. sub-Plinian eruptions) could be produced by episodic collapse events of the volcanic conduit, opposite to the mainly stationary high-

M

mass flux events (e.g. Plinian eruptions).

ED

(c) The minimum conduit radius needed to avoid conduit collapse is mainly controlled by the inlet overpressure. If an eruption proceeds following a typical depressurizing trend of the

PT

magma chamber, we expect a continuous widening process of the conduit, even with a

CE

steady magma source (e.g. water content, temperature, composition), so justifying the ubiquitous presence of lithic fragments observed in pyroclastic deposits.

AC

(d) Geometric features that conduits acquire during volcanic eruptions represent a first order factor in the eruptive dynamics, modifying parameters as the magma discharge rate, exit pressure and fragmentation level. (e) Transient models coupled with the inclusion of an additional solid phase (lithic fragments) are required for studying the temporal evolution of conduits geometry.

Acknowledgments

ACCEPTED MANUSCRIPT This research has been partially framed in the project V1 “Valutazione della pericolosità vulcanica in termini probabilistici", financed by the Italian Civil Protection Department (INGVDPC 2012-2015). During this research, Alvaro Aravena was supported by the grant CONICYT

AC

CE

PT

ED

M

AN

US

CR

IP

T

Nº72160016.

ACCEPTED MANUSCRIPT References

AC

CE

PT

ED

M

AN

US

CR

IP

T

Aadnoy, B. S., and C. Edland (2001), Borehole stability of multilateral junctions, Journal of Petroleum Science and Engineering, 30(3), 245-255. Al-Ajmi, A. M., and R. W. Zimmerman (2006), Stability analysis of vertical boreholes using the Mogi– Coulomb failure criterion, International Journal of Rock Mechanics and Mining Sciences, 43(8), 12001211. Barberi, F., R. Cioni, M. Rosi, R. Santacroce, A. Sbrana, and R. Vecci (1989), Magmatic and phreatomagmatic phases in explosive eruptions of Vesuvius as deduced by grain-size and component analysis of the pyroclastic deposits, Journal of volcanology and geothermal research, 38(3-4), 287-307. Bruce, P. M., and H. E. Huppert (1989), Thermal control of basaltic fissure eruptions, Nature, 342(6250), 665-667. Bursik, M. (1993), Subplinian eruption mechanisms inferred from volatile and clast dispersal data, Journal of Volcanology and Geothermal Research, 57(1), 57-70. Calder, E., R. Sparks, and M. Gardeweg (2000), Erosion, transport and segregation of pumice and lithic clasts in pyroclastic flows inferred from ignimbrite at Lascar Volcano, Chile, Journal of Volcanology and Geothermal Research, 104(1), 201-235. Carcano, S., L. Bonaventura, T. Esposti Ongaro, and A. Neri (2013), A semi-implicit, second-orderaccurate numerical model for multiphase underexpanded volcanic jets, Geoscientific Model Development, 6(6), 1905-1924. Carey, R., B. Houghton, J. Sable, and C. Wilson (2007), Contrasting grain size and componentry in complex proximal deposits of the 1886 Tarawera basaltic Plinian eruption, Bulletin of Volcanology, 69(8), 903-926. Caricchi, L., L. Burlini, P. Ulmer, T. Gerya, M. Vassalli, and P. Papale (2007), Non-Newtonian rheology of crystal-bearing magmas and implications for magma ascent dynamics, Earth and Planetary Science Letters, 264(3), 402-419. Cioni, R., R. Sulpizio, and N. Garruccio (2003), Variability of the eruption dynamics during a subplinian event: the Greenish Pumice eruption of Somma–Vesuvius (Italy), Journal of volcanology and geothermal research, 124(1), 89-114. Cioni, R., P. Marianelli, R. Santacroce, and A. Sbrana (2000), Plinian and subplinian eruptions, Encyclopedia of volcanoes. Academic, San Diego, 477-494. Civetta, L., G. Orsi, L. Pappalardo, R. Fisher, G. Heiken, and M. Ort (1997), Geochemical zoning, mingling, eruptive dynamics and depositional processes—the Campanian Ignimbrite, Campi Flegrei caldera, Italy, Journal of Volcanology and Geothermal Research, 75(3), 183-219. Cole, P., G. Queiroz, N. Wallenstein, J. Gaspar, A. Duncan, and J. Guest (1995), An historic subplinian/phreatomagmatic eruption: the 1630 AD eruption of Furnas volcano, Sa˜ o Miguel, Azores, Journal of volcanology and geothermal research, 69(1), 117-135. Costa, A. (2005), Viscosity of high crystal content melts: dependence on solid fraction, Geophysical Research Letters, 32(22). Costa, A., O. Melnik, and R. Sparks (2007), Controls of conduit geometry and wallrock elasticity on lava dome eruptions, Earth and Planetary Science Letters, 260(1), 137-151. Curry, D., A. M. Lourenco, L. W. Ledgerwood III, R. Maharidge, S. A. B. d. Fontoura, and N. Inoue (2016), The Effect of Borehole Pressure on the Drilling Process in Salt, SPE Drilling & Completion. de'Michieli Vitturi, M., A. Clarke, A. Neri, and B. Voight (2008), Effects of conduit geometry on magma ascent dynamics in dome-forming eruptions, Earth and Planetary Science Letters, 272(3), 567-578. de'Michieli Vitturi, M., A. Clarke, A. Neri, and B. Voight (2010), Transient effects of magma ascent dynamics along a geometrically variable dome-feeding conduit, Earth and Planetary Science Letters, 295(3), 541-553.

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN

US

CR

IP

T

de'Michieli Vitturi, M., A. Clarke, A. Neri, and B. Voight (2011), Assessing the influence of disequilibrium crystallization and degassing during magma ascent in effusive and explosive eruptions, paper presented at AGU Fall Meeting Abstracts. Degruyter, W., O. Bachmann, A. Burgisser, and M. Manga (2012), The effects of outgassing on the transition between effusive and explosive silicic eruptions, Earth and Planetary Science Letters, 349, 161170. Dominguez, L., L. Pioli, C. Bonadonna, C. Connor, D. Andronico, A. Harris, and M. Ripepe (2016), Quantifying unsteadiness and dynamics of pulsatory volcanic activity, Earth and Planetary Science Letters, 444, 160-168. Dragoni, M., and S. Santini (2007), Lava flow in tubes with elliptical cross sections, Journal of volcanology and geothermal research, 160(3), 239-248. Gholami, R., A. Moradzadeh, V. Rasouli, and J. Hanachi (2014), Practical application of failure criteria in determining safe mud weight windows in drilling operations, Journal of Rock Mechanics and Geotechnical Engineering, 6(1), 13-25. Gurioli, L., B. F. Houghton, K. V. Cashman, and R. Cioni (2005), Complex changes in eruption dynamics during the 79 AD eruption of Vesuvius, Bulletin of Volcanology, 67(2), 144-159. Hess, K., and D. Dingwell (1996), Viscosities of hydrous leucogranitic melts: A non-Arrhenian model, American Mineralogist, 81(9-10), 1297-1300. Hoek, E., and E. Brown (1997), Practical estimates of rock mass strength, International Journal of Rock Mechanics and Mining Sciences, 34(8), 1165-1186. Houghton, B., C. Wilson, P. Del Carlo, M. Coltelli, J. Sable, and R. Carey (2004), The influence of conduit processes on changes in style of basaltic Plinian eruptions: Tarawera 1886 and Etna 122 BC, Journal of Volcanology and Geothermal Research, 137(1), 1-14. Jurado-Chichay, Z., and G. Walker (2001), The intensity and magnitude of the Mangaone subgroup plinian eruptions from Okataina Volcanic Centre, New Zealand, Journal of volcanology and geothermal research, 111(1), 219-237. Kennedy, B., O. Spieler, B. Scheu, U. Kueppers, J. Taddeucci, and D. B. Dingwell (2005), Conduit implosion during Vulcanian eruptions, Geology, 33(7), 581-584. Koyaguchi, T., Y. J. Suzuki, and T. Kozono (2010), Effects of the crater on eruption column dynamics, Journal of Geophysical Research: Solid Earth, 115(B7). La Spina, G., and M. de'Michieli Vitturi (2012), High-resolution finite volume central schemes for a compressible two-phase model, SIAM Journal on Scientific Computing, 34(6), B861-B880. La Spina, G., M. Burton, and M. d. M. Vitturi (2015), Temperature evolution during magma ascent in basaltic effusive eruptions: A numerical application to Stromboli volcano, Earth and Planetary Science Letters, 426, 89-100. Le Métayer, O., J. Massoni, and R. Saurel (2005), Modelling evaporation fronts with reactive Riemann solvers, Journal of Computational Physics, 205(2), 567-610. Llewellin, E., and M. Manga (2005), Bubble suspension rheology and implications for conduit flow, Journal of Volcanology and Geothermal Research, 143(1), 205-217. Macedonio, G., F. Dobran, and A. Neri (1994), Erosion processes in volcanic conduits and application to the AD 79 eruption of Vesuvius, Earth and planetary science letters, 121(1), 137-152. Macedonio, G., A. Neri, J. Martì, and A. Folch (2005), Temporal evolution of flow conditions in sustained magmatic explosive eruptions, Journal of volcanology and geothermal research, 143(1), 153-172. Mader, H., E. Llewellin, and S. Mueller (2013), The rheology of two-phase magmas: A review and analysis, Journal of Volcanology and Geothermal Research, 257, 135-158. Maeno, F., and H. Taniguchi (2009), Sedimentation and welding processes of dilute pyroclastic density currents and fallout during a large-scale silicic eruption, Kikai caldera, Japan, Sedimentary Geology, 220(3), 227-242.

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN

US

CR

IP

T

Manga, M., and M. Loewenberg (2001), Viscosity of magmas containing highly deformable bubbles, Journal of Volcanology and Geothermal Research, 105(1), 19-24. Pioli, L., B. Azzopardi, and K. Cashman (2009), Controls on the explosivity of scoria cone eruptions: Magma segregation at conduit junctions, Journal of Volcanology and Geothermal Research, 186(3), 407415. Piombo, A., A. Tallarico, and M. Dragoni (2016), Role of mechanical erosion in controlling the effusion rate of basaltic eruptions, Geophysical Research Letters, 43(17), 8970-8977. Rahman, M., D. Naseby, and S. Rahman (2000), Borehole collapse analysis incorporating timedependent pore pressure due to mud penetration in shales, Journal of Petroleum Science and Engineering, 28(1), 13-31. Romenski, E., D. Drikakis, and E. Toro (2010), Conservative models and numerical methods for compressible two-phase flow, Journal of Scientific Computing, 42(1), 68-95. Rust, A., and K. V. Cashman (2004), Permeability of vesicular silicic magma: inertial and hysteresis effects, Earth and Planetary Science Letters, 228(1), 93-107. Scandone, R., and S. D. Malone (1985), Magma supply, magma discharge and readjustment of the feeding system of Mount St. Helens during 1980, Journal of Volcanology and Geothermal Research, 23(3-4), 239-262. Single, B., R. Goel, V. Mehrotra, S. Garg, and M. Allu (1998), Effect of intermediate principal stress on strength of anisotropic rock mass, Tunnelling and Underground Space Technology, 13(1), 71-79. Sobol, I. M. (2001), Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates, Mathematics and computers in simulation, 55(1), 271-280. Sparks, R. S. J. (1978), The dynamics of bubble formation and growth in magmas: a review and analysis, Journal of Volcanology and Geothermal Research, 3(1-2), 1-37. Starostin, A., A. Barmin, and O. Melnik (2005), A transient model for explosive and phreatomagmatic eruptions, Journal of Volcanology and Geothermal Research, 143(1), 133-151. Truby, J., S. Mueller, E. Llewellin, and H. Mader (2015), The rheology of three-phase suspensions at low bubble capillary number, paper presented at Proc. R. Soc. A, The Royal Society. Varekamp, J. C. (1993), Some remarks on volcanic vent evolution during plinian eruptions, Journal of volcanology and geothermal research, 54(3-4), 309-318. Vernik, L., and M. D. Zoback (1992), Estimation of maximum horizontal principal stress magnitude from stress‐induced well bore breakouts in the Cajon Pass Scientific Research borehole, Journal of Geophysical Research: Solid Earth, 97(B4), 5109-5119. Zhang, Y. (1999), H2O in rhyolitic glasses and melts: measurement, speciation, solubility, and diffusion, Reviews of Geophysics, 37(4), 493-516.

ACCEPTED MANUSCRIPT Tables with captions Table 1. Mohr - Coulomb and Mogi - Coulomb criteria for collapse pressure [Al-Ajmi and Zimmerman, 2006]. : Minimum pressure needed to avoid conduit collapse. : Radial stress.

: Maximum horizontal stress.

: Minimum horizontal stress.

: Vertical stress. : Poisson ratio.

US

.

IP

. : Cohesion. : Angle of internal friction.

CR

. .

: Tangential stress.

T

: Vertical stress.

.

AN

.

: Pore pressure.

M

.

.

ED

Mohr – Coulomb Collapse occurs if

Case

PT

1 2

1 2 3

AC

Case

CE

3

Mogi – Coulomb Collapse occurs if

.

ACCEPTED MANUSCRIPT Table 2. Parameters related to the application of stability criteria, representative of typical conditions of country rocks in nature [Hoek and Brown, 1997]. As a reference, reported values of cohesion ranges between 1 MPa and 10 MPa, whereas the friction angle varies between and 45º. Value

IP

T

Parameter (Symbol)

5 MPa

Rock cohesion ( )

,

)

AC

CE

PT

ED

M

AN

Both horizontal stress gradients (

38º

26 kPa/m

)

US

Vertical stress gradient (

CR

Angle of friction ( )

18 kPa/m

30º

ACCEPTED MANUSCRIPT Table 3. Input and geometric parameters used in the numerical simulations. (1) In order to test the influence of conduit length, we performed some additional simulations varying the values of conduit length and inlet pressure. They are reported in supplementary Figures S3 and S4, and section 3.1.1. (2) Some additional simulations were performed in order to improve the quality of

AC

CE

PT

ED

-

M

-

IP

Geometry NC2 4.5 – 6.5 wt. % 115 – 135 MPa 850ºC 5000 m 35 m – 75 m

1.1 – 1.2

1.1 – 1.2

Fragmentation level of the cylindrical case

Fragmentation level of the cylindrical case

100 m – 1000 m

100 m – 1000 m

AN

-

Geometry NC1 4.5 – 6.5 wt. % 115 – 135 MPa 850ºC 5000 m 35 m – 75 m

CR

Geometry C 4.5 – 6.5 wt. % 115 – 135 MPa 850ºC 5000 m 25 – 115 m -

US

Parameter Water content Inlet pressure (1) Temperature Conduit length (1) Fixed radius (2) Minimum radius Maximum radius/Minimum radius Depth of the central point of the variableradius domain Length of the variableradius zone of the conduit

T

the interpolations of critical radius (defined in section 3.1).

ACCEPTED MANUSCRIPT Table 4. Notation summary. Units

AC

CE

PT

ED

M

IP CR US

AN

Name Volumetric fraction Adiabatic exponent Drag factor Solubility exponent Parameter for controlling the crystal influence on viscosity Viscosity Density Solubility coefficient Relaxation parameter for controlling the crystallization rate Relaxation parameter for controlling the gas exsolution rate Relaxation parameter for controlling the pressure disequilibrium Parameter for controlling the direction of the viscous term contribution Relative velocity between phases Auxiliary parameter in equations of state Auxiliary parameter in equations of state Specific heat capacity Fitting parameter of rheological model Fitting parameter of rheological model Fitting parameter of rheological model Internal energy Formation energy Fitting parameter of outgassing model Throat-bubble size ratio Acceleration of gravity Darcian permeability Inertial permeability Fitting parameter of outgassing model Pressure Average bubble size Specific entropy Velocity Dissolved water concentration Mass fraction Vertical coordinate Fitting parameter of rheological model Sound speed Auxiliary parameter in rheological model Bubble density number Conduit radius Temperature Global subscripts and superscripts Referred to Crystals Critical value Dissolved gas Equilibrium value Exsolved gas Phases 1 and 2 Continuous phase Components of phase 1 Liquid component of phase 1 Maximum value At constant volume Initial or reference value Phase 1 Phase 2 No subscript Mixture or global parameter

T

Symbol

ACCEPTED MANUSCRIPT

US

CR

IP

T

Figures with captions

AN

Figure 1. (a) Cylindrical conduit (geometry C). (b) Cylindrical conduit with a local enlargement

M

along a vertical distance ( ) (geometry NC1). (c) Conduit with two coaxial cylindrical portions

AC

CE

PT

ED

connected by a linear enlarging transitional zone of length

(geometry NC2).

US

CR

IP

T

ACCEPTED MANUSCRIPT

AN

Figure 2. Depth versus pressure in the conduit and minimum pressures needed to avoid conduit collapse, using both stability criteria. These results were obtained in a specific simulation (radius

M

of 35 m, water content of 6.5% and inlet pressure of 125 MPa), and they represent a typical case

AC

CE

PT

ED

of an unstable conduit.

US

CR

IP

T

ACCEPTED MANUSCRIPT

Figure 3. (a) Instability index versus fragmentation level, for different radii in cylindrical

AN

conduits. (b) Instability index versus mass discharge rate, for different radii in cylindrical

M

conduits. We present here the results related to simulations with variable values of inlet pressure

AC

CE

PT

ED

(115 – 135 MPa) and water content (4.5 – 6.5 wt. %).

US

CR

IP

T

ACCEPTED MANUSCRIPT

Figure 4. Critical radius (i.e. minimum radius for a stable conduit) for different conditions of

AC

CE

PT

ED

M

AN

inlet pressure and water content. (a) Mohr – Coulomb criterion. (b) Mogi – Coulomb criterion.

M

AN

US

CR

IP

T

ACCEPTED MANUSCRIPT

Figure 5. Summary of results related to the sensitivity analysis of critical radius. (a) Critical

ED

radius (Mohr - Coulomb) versus inlet overpressure. (b) Critical radius (Mohr - Coulomb) versus

PT

water content. (c) Critical radius (Mohr - Coulomb) versus conduit length. (d) Critical radius (Mogi - Coulomb) versus inlet overpressure. (e) Critical radius (Mogi - Coulomb) versus water

CE

content. (f) Critical radius (Mogi - Coulomb) versus conduit length. (g) First order sensitivity

AC

index. (h) Total effect sensitivity index. We used 587 different sets of inlet parameters (4,696 simulations).

CR

IP

T

ACCEPTED MANUSCRIPT

US

Figure 6. (a) Deepening of the fragmentation level versus ratio of exit pressures for geometry

AN

NC1 and the equivalent cylindrical case. (b) Ratio of mass discharge rates versus ratio of final velocities for geometry NC1 and the equivalent cylindrical case. (c) Deepening of the

M

fragmentation level versus ratio of exit pressures for geometry NC2 and the equivalent

ED

cylindrical case. (d) Ratio of mass discharge rates versus ratio of final velocities for geometry NC2 and the equivalent cylindrical case. We present here results related to a set of simulations

AC

CE

PT

with variable values of inlet pressure (115 MPa – 135 MPa) and

(35 – 75 m).

ACCEPTED MANUSCRIPT Highlights: Numerical modelling of cylindrical and non-cylindrical volcanic conduits.



Application of collapse criteria for studying volcanic conduits stability.



Mechanical stability of a conduit is mainly controlled by its radial dimensions.



Geometry acquired by the conduit as it is widened influences the eruptive dynamics.

AC

CE

PT

ED

M

AN

US

CR

IP

T

