Audrey Michaud-Dubuy
email: michauddubuy@ipgp.fr

Guillaume Carazzo

Edouard Kaminski

Frédéric Girault

Audrey Michaud

A revisit of the role of gas entrapment on the stability conditions of explosive volcanic columns

Keywords: Plinian eruptions, particle sedimentation, grain-size distribution, fragmentation, open porosity, pyroclastic density currents

Highlights

• We study the collapse conditions using a 1D model with grain-size distribution (GSD).

• Combined effects of gas entrapment and variable entrainment favor column collapse.

• This drastic effect is reduced when accounting for open porosity in the model.

• We predict GSD in pyroclastic flows as a function of GSD and MDR at the vent.

• We test the model against the historical 186 CE Taupo and 79 CE Vesuvius eruptions.

Introduction

Explosive volcanic eruptions stand as one of the most powerful and dangerous natural phenomena on Earth. During these extreme events, the magma ascending from depth is fragmented in the conduit and expelled at the vent as a dense turbulent mixture of hot gas and pyroclasts. The mass discharge rate of these high-velocity turbulent jets usually ranges between 10 6 and 10 9 kg s -1 [START_REF] Carey | The intensity of Plinian eruptions[END_REF]. Depending on its mass discharge rate, a volcanic jet can follow remarkably different dynamical evolutions during its rise in the atmosphere. In the buoyant regime (also called "Plinian" regime), the volcanic mixture forms a vertical column that rises up to tens of kilometers before spreading out laterally to form an horizontal umbrella cloud [START_REF] Sparks | The dimensions and dynamics of volcanic eruption columns[END_REF]. In the fountain regime (or "collapse" regime), the turbulent jet collapses to the ground and produces pyroclastic density currents (PDC) rushing down the volcano flanks. When eruption conditions are close to those of column collapse, the regime is transitional: an unstable buoyant column still rises to high altitudes but occasionally generates PDC. The two eruptive regimes, which can occur one after another and even alternate during the same eruption, mainly determine the associated hazards. Whereas a rain of ash and pumices produced during the Plinian regime may cause infrastructure damages (e.g., [START_REF] Wilson | Volcanic hazard impacts to critical infrastructure: A review[END_REF], major perturbations of air traffic (e.g., [START_REF] Miller | Volcanic ash hazards to aviation[END_REF][START_REF] Schmidt | Assessing hazards to aviation from sulfur dioxide emitted by explosive Icelandic eruptions[END_REF] and breathing difficulties (e.g., [START_REF] Horwell | The respiratory health hazards of volcanic ash: a review for volcanic risk mitigation[END_REF][START_REF] Horwell | Grain-size analysis of volcanic ash for the rapid assessment of respiratory health hazard[END_REF][START_REF] Horwell | Physicochemical and toxicological profiling of ash from the 2010 and 2011 eruptions of Eyjafjallajökull and Grímsvötn volcanoes, Iceland using a rapid respiratory hazard assessment protocol[END_REF], PDC most commonly lead to massive human and material losses (e.g., [START_REF] Spence | Building vulnerability and human casualty estimation for a pyroclastic flow: a model and its application to Vesuvius[END_REF][START_REF] Wilson | Volcanic hazard impacts to critical infrastructure: A review[END_REF]. The prediction of column behavior is therefore fundamental to assess the impact of explosive volcanic eruptions. The understanding and prediction of the source (and environmental) conditions leading to either a buoyant plume or a collapsing fountain remains a major goal in physical volcanology.

The dual behavior described above rises from the evolution of the volcanic column bulk density during its ascent in the atmosphere. At the vent, the hot mixture of gas and pyroclasts has a bulk density greater than the ambient air, and its initial momentum drives the plume ascent. Thereafter, owing to turbulent mixing, cold atmospheric air is entrained into the flow and heated by the hot pyroclasts, leading to a rapid expansion of the gas and an associated decrease of the bulk density of the jet [START_REF] Sparks | A model for the formation of ignimbrite by gravitational column collapse[END_REF][START_REF] Woods | The fluid dynamics and thermodynamics of eruption columns[END_REF][START_REF] Woods | The dynamics of explosive volcanic eruptions[END_REF]. In the meantime, the jet momentum decreases with altitude due to its negative buoyancy. A stable Plinian column forms when the bulk density of the volcanic mixture becomes lower than that of the atmospheric air before complete exhaustion of its initial momentum. Then, the volcanic plume rises by natural convection until it reaches a neutral buoyancy level and spreads out laterally. In the fountain regime, the jet consumes its initial momentum before becoming buoyant and thus collapses to the ground producing PDC.

Since the 1970's, various methods have been used to study the stability of volcanic columns produced by explosive eruptions (or equivalently the conditions leading to PDC production). A first generation of theoretical models were developed based on a simplified 1D approach [START_REF] Wilson | Explosive volcanic eruptions: III. Plinian eruption columns[END_REF][START_REF] Sparks | The dimensions and dynamics of volcanic eruption columns[END_REF][START_REF] Woods | The fluid dynamics and thermodynamics of eruption columns[END_REF] stemming from the widely used "tophat" formalism of [START_REF] Morton | Turbulent gravitational convection from maintained and instantaneous source[END_REF]. These models rely on the "dusty-gas" hypothesis where particles are considered small enough to remain in thermal and mechanical equilibrium with the gas. In that case, the volcanic mixture can be represented as a single "equivalent" gaseous phase. This convenient formalism has been used to quantitatively study the relationship between the maximum height of the plume and the source mass discharge rate [START_REF] Settle | Volcanic eruption clouds and the thermal power output of explosive eruptions[END_REF][START_REF] Wilson | Explosive volcanic eruptions-IV. The control of magma properties and conduit geometry on eruption column behaviour[END_REF][START_REF] Carey | The intensity of Plinian eruptions[END_REF] and to determine the conditions leading to column collapse. [START_REF] Wilson | Explosive volcanic eruptions-IV. The control of magma properties and conduit geometry on eruption column behaviour[END_REF] showed that the mass fraction of free gas in the eruptive mixture at the vent (gas + pyroclastic fragments) and mass discharge rate (MDR) strongly control the transition between the stable Plinian plume and the collapsing fountain regimes, a conclusion later confirmed by laboratory experiments [START_REF] Woods | A laboratory study of explosive volcanic eruptions[END_REF]. Other effects have been considered in order to determine quantitatively this regime boundary. [START_REF] Woods | Particle fallout, thermal disequilibrium and volcanic plumes[END_REF] demonstrated that particle sedimentation has a minor effect on the column dynamics, but thermal disequilibrium significantly changes the column behavior and can induce a column collapse. [START_REF] Woods | The decompression of volcanic jets in a crater during explosive volcanic eruptions[END_REF] and [START_REF] Koyaguchi | Effects of the crater on eruption column dynamics[END_REF] calculated the conditions for which jet decompression in a crater can yield subsonic velocities and column collapse. [START_REF] Degruyter | Impact of wind on the condition for column collapse of volcanic plumes[END_REF] showed that high velocity atmospheric winds can significantly increase the amount of air engulfed in the volcanic column hence making it more stable.

More sophisticated 2D axisymetric and 3D models, based on the time-dependent solution of the Navier-Stokes and energy conservation equations, have been developed to describe the fluid dynamics of the eruptive mixture and the surrounding atmosphere [START_REF] Valentine | Numerical models of Plinian eruption columns and pyroclastic flows[END_REF][START_REF] Neri | Influence of eruption parameters on the thermofluid dynamics of collapsing volcanic columns[END_REF][START_REF] Suzuki | A numerical study of turbulent mixing in eruption clouds using a three-dimensional fluid dynamics model[END_REF][START_REF] Esposti Ongaro | Transient 3D numerical simulations of column collapse and pyroclastic density current scenarios at Vesuvius[END_REF]. These models improved our understanding of the dynamics of explosive volcanic columns by exploring ranges of parameters that are beyond the limitations imposed by 1D models. However, a recent intercomparison exercise revealed that the results of 2D and 3D models diverge in some of their predictions depending on the assumptions made to solve the governing equations (Costa et al., 2016a). On the other hand, 1D models predictions show similarities with those of 3D models, suggesting that 1D models can be used to adequately describe the general behavior of volcanic columns.

Among the different studies on the stability of volcanic plumes, two - Kaminski and Jaupart (2001) and Carazzo et al. (2008a) -have compared the model predictions of collapse with field constraints in a systematic way, and have shown that previous models tended to significantly favor the buoyant regime in comparison to the natural cases. Kaminski and base of the column and promote column collapse. Carazzo et al. (2008a) studied the effects of reduced turbulent entrainment due to negative buoyancy at the base of the column and showed that it also promotes column collapse. Gas entrapment and reduced entrainment act in the same direction and, when taken independently, yield model predictions consistent with field data. However, these two phenomena have not been considered together yet, and it can be argued that their combined effect could be to favor too much the collapse regime compared to the natural cases, hence reducing their performances in determining accurately the conditions of collapse. Here we propose to combine the two effects, which requires to take into account the evolution of the total grain size distribution (TGSD) due to sedimentation [START_REF] Girault | The effect of total grainsize distribution on the dynamics of turbulent volcanic plumes[END_REF].

Our paper is organized as follows. In section 2, we describe our physical 1D model of explosive volcanic columns. In section 3, we analyse the results of the model taking into account variable entrainment as a function of buoyancy, sedimentation, gas entrapment, and open porosity, and we investigate the control of these parameters and phenomena on the conditions leading to column collapse. We further quantify the fountain height and characterize the PDC produced in the collapse regime. In section 4, we compare our theoretical predictions with field data from two well-known historical eruptions (Taupo and Vesuvius). In section 5, we discuss additional effects that may influence the transition between stable and collapsing eruptive columns, and we conclude in section 6.

Physical model of explosive volcanic columns

Conservation equations and constitutive laws

Our model relies on a 1D steady-state "top-hat" formalism for a conical jet in which all the dynamical variables are considered constant inside the jet and zero outside [START_REF] Morton | Turbulent gravitational convection from maintained and instantaneous source[END_REF]. The horizontal rate of entrainment of the surrounding fluid is assumed to be proportional to the local vertical ascension rate of the plume, through a constant entrainment coefficient 𝛼 ! . Our work is based on the improved version of the [START_REF] Woods | The fluid dynamics and thermodynamics of eruption columns[END_REF] model, that explicitely considers the conservation of energy and its effect on the evolution of the bulk density of the flow. As in [START_REF] Woods | The fluid dynamics and thermodynamics of eruption columns[END_REF], we consider thermal and mechanical equilibrium between the volcanic gas and the particles.

For a particle-laden volcanic jet and a calm stratified atmosphere, the three macroscopic conservation equations of mass, momentum and energy flow rates in steady-state are written as [START_REF] Woods | The fluid dynamics and thermodynamics of eruption columns[END_REF][START_REF] Woods | Particle fallout, thermal disequilibrium and volcanic plumes[END_REF][START_REF] Bursik | Effect of wind on the rise height of volcanic plumes[END_REF][START_REF] Costa | A three-dimensional Eulerian model for transport and deposition of volcanic ashes[END_REF][START_REF] Girault | The effect of total grainsize distribution on the dynamics of turbulent volcanic plumes[END_REF]:

𝑑 𝑑𝑧 𝜌𝑈𝑅 ! = 2𝜌 ! 𝑈 ! 𝑅 + 𝑑𝑄 ! 𝑑𝑧 ! ! !!! , (1)
𝑑 𝑑𝑧 𝜌𝑈 ! 𝑅 ! = 𝑔 𝜌 ! -𝜌 𝑅 ! + 𝑈 𝑑𝑄 ! 𝑑𝑧 ! ! !!! , (2)
𝑑 𝑑𝑧 𝜌𝑈𝑅 ! 𝑐 ! 𝑇 = 2𝜌 ! 𝑈 ! 𝑅𝑐 ! 𝑇 ! -𝜌 ! 𝑔𝑈𝑅 ! + 𝑐 ! 𝑇 𝑑𝑄 ! 𝑑𝑧 ! ! !!! , (3)
where R(z) is the column radius, U(z) is the vertical velocity, g is the acceleration of gravity, 𝜌, 𝑐 ! and T are the density, the specific heat and the temperature at constant pressure of the bulk mixture, respectively, 𝜌 ! , 𝑐 ! and 𝑇 ! are those of the atmosphere (all variables are defined in the Notation section). 𝑄 ! = 𝑥 ! 𝜌𝑈𝑅 ! , where 𝑥 ! is the mass proportion of 𝜙-sized particles in the GSD, is the discharge rate of 𝜙-sized particles (kg s -1 ) distributed within 20 classes of grain sizes ranging from 𝑑 !"# = 10 𝜙 (1 𝜇m) to 𝑑 !"# = -9 𝜙 (0.5 m) with one 𝜙 intervals.

Magma temperature is taken as the average of andesitic magma (T 0 = 1200 K). As illustrated in Fig. 1, 𝑈 ! is the entrainment rate at the edge of the plume, defined as 𝑈 ! = 𝛼 ! 𝑈 [START_REF] Morton | Turbulent gravitational convection from maintained and instantaneous source[END_REF], which is here function of the column buoyancy relative to the ambient air, and is expressed as [START_REF] Kaminski | Turbulent entrainment in jets with arbitrary buoyancy[END_REF]:

𝛼 ! = 𝐶 2 + 1 - 1 𝐴 𝑅𝑖, (4)
where Ri = g ρ a -ρ ( )R / ρ a U 2 is the Richardson number in the plume, and A and C are dimensionless parameters depending on the flow structure. C is taken as a constant (≈ 0.135)

whereas A evolves as a function of the downstream distance from the source and buoyancy

(see Appendix A).
To follow the evolution of the density in the plume as a function of entrainement and temperature variations with height, we use the same constitutive equations as in [START_REF] Woods | The fluid dynamics and thermodynamics of eruption columns[END_REF]:

1 𝜌 = 1 -𝑥 ! 𝜌 ! + 𝑥 ! 𝑅 ! 𝑇 𝑃 ! , (5)
𝑥 ! = 1 + 𝑥 ! ! -1 𝜌 ! 𝑈 ! 𝑅 ! ! 𝜌𝑈𝑅 ! , (6)
𝑅 ! = 𝑅 ! + 𝑅 ! ! -𝑅 ! 1 -𝑥 ! 𝑥 ! 𝑥 ! ! 1 -𝑥 ! ! , (7)
𝑐 ! = 𝑐 ! + 𝑐 ! ! -𝑐 ! 1 -𝑥 ! 1 -𝑥 ! ! , (8)
where 𝜌 ! = 2000 kg m -3 is the average density of the particles, 𝑥 ! (𝑧) is the effective gas mass fraction, 𝑅 ! ! = 461 J K -1 kg -1 and 𝑅 ! = 287 J K -1 kg -1 are the bulk column and the air gas constants, respectively, 𝑃 ! (𝑧) is the atmospheric pressure, and the subscript 0 denotes a value at the vent. 𝑃 ! 𝑧 , 𝑇 ! (𝑧) and 𝜌 ! (𝑧) depend on the atmospheric conditions which we take as typical polar, mid-latitude or tropical [START_REF] Glaze | Sensitivity of buoyant plume heights to ambient atmospheric conditions: implications for volcanic eruption columns[END_REF]Carazzo et al., 2008b).

Particle sedimentation

To account for the mass loss of 𝜙-sized particles from the edges of the column, we consider it to be proportional to the mass discharge rate of particles 𝑄 ! and to the terminal fall velocity 𝑉 ! , such as [START_REF] Woods | Particle fallout, thermal disequilibrium and volcanic plumes[END_REF][START_REF] Ernst | Sedimentation from turbulent jets and plumes[END_REF][START_REF] Girault | The effect of total grainsize distribution on the dynamics of turbulent volcanic plumes[END_REF]:

𝑑𝑄 ! 𝑑𝑧 = -𝑝 ! 𝑄 ! 𝑅 𝑉 ! 𝑈 , (9)
where 𝑝 ! is a probability of sedimentation experimentally determined and taken equal to 0.27 ± 0.01 [START_REF] Ernst | Sedimentation from turbulent jets and plumes[END_REF][START_REF] Girault | The effect of total grainsize distribution on the dynamics of turbulent volcanic plumes[END_REF][START_REF] Sigurdsson | and Sigurdsson[END_REF]. For a given particle size, the fallout velocity 𝑉 ! is calculated using the formulae of [START_REF] Bonadonna | Thickness variations and volume estimates of tephra fall deposits: the importance of particle Reynolds number[END_REF]:

where 𝑑 ! is the particle diameter of a 𝜙-sized particle, 𝜇(𝑧) is the dynamic viscosity of air [START_REF] Sutherland | The viscosity of gases and molecular forces[END_REF], and 𝑅𝑒 ! = 𝜌𝑑 ! 𝑉 ! /𝜇 is the particle Reynolds number.

V φ = 3.1d φ g ρ p -ρ ( ) ρ a d φ 4g 2 ρ p -ρ ( ) 2 225µρ a ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ 1/3 d φ 2 g ρ p -ρ ( ) 18µ ⎧ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ for , for , (10)
for ,

Grain-size distribution and amount of gas at the vent

In explosive eruptions, GSD of pyroclasts results from a fragmentation sequence in the conduit before the eruption. The rapid decompression of magma during its ascent causes a "primary" fragmentation [START_REF] Alibidirov | Magma fragmentation by rapid decompression[END_REF], which disintegrates bubbly magma into fragments, then followed by a "secondary" fragmentation of larger fragments into finer ash (Kaminski and Jaupart, 1998). During this sequence, magmatic gas separates into two phases: an entrapped one contained in bubbles within the clasts, and a continuous one carrying fragments and ashes. This latter phase corresponds to the effective amount of "free"

gas that has to be considered for the modeling of the turbulent flow. Within this framework, clast size plays a key role in setting the amount of gas released at fragmentation: large fragments (pumices) retain a larger amount of gas than small fragments (ashes) do. Field data (Kaminski and Jaupart, 1998) and fragmentation experiments (e.g., [START_REF] Kueppers | Explosive energy" during vol-canic eruptions from fractal analysis of pyroclasts[END_REF] have shown that volcanic rocks fragment according to a power-law distribution:

𝑁 𝑅 ≥ 𝑟 = 𝜆𝑟 !! , 11
where 𝑁 𝑅 ≥ 𝑟 is the number of fragments of size R larger than or equal to 𝑟, 𝜆 is a normalization constant and 𝐷 is the power-law exponent.

The two main parameters that control the stability of a volcanic plume are the mass discharge rate and the momentum flow rate at the vent, or, equivalently, the mass discharge rate and exit velocity [START_REF] Wilson | The control of volcanic column heights by eruption energetics and dynamics[END_REF][START_REF] Woods | The fluid dynamics and thermodynamics of eruption columns[END_REF][START_REF] Sparks | Volcanic Plumes[END_REF]. Because the exit velocity is mainly controled by the amount of gas in the volcanic mixture, 𝑥 ! ! , (e.g., [START_REF] Wilson | Explosive volcanic eruptions-IV. The control of magma properties and conduit geometry on eruption column behaviour[END_REF][START_REF] Koyaguchi | Effects of the crater on eruption column dynamics[END_REF]see Appendix B) the transition between the Plinian and Fountain regimes is often given as a threshold mass discharge rate for a given gas mass fraction. Here, we calculate 𝑥 ! ! as the total fraction of exsolved gas minus the amount of exsolved gas trapped in large particles, using the model of Kaminski and Jaupart (1998;[START_REF] Bursik | Effect of wind on the rise height of volcanic plumes[END_REF].

If bubbles inside a fragment are not connected to the exterior, the fraction of gas released by each fragment is given by

! !"# ! !"# = !!
! for particle with a radius 𝑟 larger than 3𝑏, and

! !"# ! !"#
= 1 for particles smaller than 3𝑏, where 𝑏 is the mean bubble radius in the magma. These relationships can then be used to calculate the total gas released at fragmentation (Fig. 2a)

provided that the TGSD -hence D -is known, as well as the smallest and largest particle sizes 𝑟 !"# and 𝑟 !"# , and the mean bubble size. As discussed by Kaminski and Jaupart (1998), the exact values of 𝑟 !"# and 𝑟 !"# do not significantly affect the calculations, and uncertainties on their determination change the result by ± 5% only. Indeed, where fine particles dominate (i.e., D > 3), the fragments do not significantly entrap gas since all their bubbles are connected to the exterior. The exact value of 𝑟 !"# is therefore not critical as long as it is much smaller than 3𝑏. On the other hand, where large particles dominate (i.e., D < 3) the average fragment size is much larger than the bubble size, hence fragments do not significantly release gas. The results are also weakly sensitive to the exact value of 𝑏 or to more complex bubble size distribution in the magma. Fig. 2a gives the fraction of gas released at fragmentation as a function of D for three different values of the mean bubble size (from 10 𝜇m to 1 mm) and

shows that this does not change the results by more than 5%. The main parameter controlling the mass fraction released at fragmentation is therefore the power law exponent D. For the rest of the paper we set 𝑟 !"# = 0.5 𝜇m, 𝑟 !"# = 0.25 m (see section 2.1), and 𝑏 = 100 𝜇m, and we take a power-law exponent D ranging from 2.5 (coarsest distribution) to 3.3 (finest distribution), a value beyond which the results are not longer affected by the precise value of

D.

Several studies have shown that inside pumices, some of the bubbles are connected to the exterior (e.g., [START_REF] Toramaru | Formation of propagation pattern in two-phase flow systems with application to volcanic eruptions[END_REF][START_REF] Klug | Permeability development in vesiculating magmas : implications for fragmentation[END_REF] and contribute to an additional release of gas through this "open porosity". Measurements made on trachytic pumices from the Vesuvius 79 CE [START_REF] Shea | Transitions between fall phases and pyroclastic density currents during the AD79 eruption at Vesuvius: building a transient conduit model from the textural and volatile record[END_REF], on dacitic pumices from the Novarupta 1912 [START_REF] Nguyen | Explosive to effusive transition during the largest volcanic eruption of the 20th century (Novarupta 1912, Alaska)[END_REF], on andesitic pumices from Soufriere Hills 1997 [START_REF] Formenti | Vesicle connectivity in pyroclasts and implications for the fluidization of fountain-collapse pyroclastic flows, Montserrat (West Indies)[END_REF],

Taranaki 1655 CE [START_REF] Platz | Transition from effusive to explosive phases in andesite eruptions -a case-study from the AD 1655 eruption of Mt. Taranaki, New Zealand[END_REF], and Lascar 1993 [START_REF] Formenti | Vesicle connectivity in pyroclasts and implications for the fluidization of fountain-collapse pyroclastic flows, Montserrat (West Indies)[END_REF], and on rhyolitic pumice from Kos Plateau Tuff 161,000 BP (Bouvet de [START_REF] Bouvet De Maisonneuve | Characterization of juvenile pyroclasts from the Kos Plateau Tuff (Aegean Arc): insights into the eruptive of a large rhyolitic eruption[END_REF],

and Mount Mazama 7700 BP eruptions [START_REF] Klug | Structure and physical characteristics of pumice from the climactic eruption of Mount Mazama (Crater Lake), Oregon[END_REF], lead to open porosity between 60 to 70% for explosive eruptions regardless of the magma composition.

To take into account open porosity in our model, we relax the gas entrapment hypothesis of Kaminski and Jaupart (1998) by introducing a new parameter 𝜉 =

! !"#$%&&'( ! !"#$$%&
where 𝑉 !"#$%&&'( is the volume of gas initially entraped in the fragment but now released

through the connected open porosity, and 𝑉 !"#$$%& is the volume of gas initially trapped into the fragments at fragmentation. For particles with a radius larger than 3𝑏, the definition of 𝜉 leads to

! !"# ! !"# = !!(!!!)!! ! !
, and we keep

! !"# ! !"#
= 1 for particles smaller than 3𝑏. Thus, in the extreme case where no open porosity develops, the volume of gas initially entraped but released is null (𝜉 = 0) and

! !"# ! !"# = !!
! . On the other hand, if the open porosity reaches 100%, all the gas initially entraped is released (𝜉 = 1) and

! !"# ! !"#
= 1 for all particles. Fig. 2b gives the effective fraction of gas released after fragmentation and the development of an open porosity as a function of D for four different values of 𝜉 (0, 20, 65 and 80%).

Exit velocity at the base of the eruptive column

In explosive eruptions, the volcanic mixture generally exits the vent at a sonic velocity (𝑈 ! ) and with a pressure (𝑃 ! ) larger than the atmospheric pressure (𝑃 ! ). After its rapid decompression to the atmospheric pressure, the velocity of the mixture 𝑈 ! is supersonic and can be expressed as [START_REF] Woods | The decompression of volcanic jets in a crater during explosive volcanic eruptions[END_REF]:

𝑈 ! = 𝑈 ! + 𝑆 𝑄 𝑃 ! -𝑃 ! ( 12
)
where 𝑆 is the cross-sectional area of the conduit, and 𝑄 is the mass discharge rate feeding the eruption. In the case of a free decompression (i.e., not controled by the shape of the crater)

and for a cylindrical conduit, it is possible to obtain 𝑈 ! as a function of the amount of gas in the mixture for a given mass discharge rate (see Appendix B). We will use these conditions in the next parts of the article and we will discuss in section 5.1 the consequences of the presence of a crater.

For a given set of source conditions, the dynamical parameters of the plume, i.e. its velocity, radius, density and temperature are then calculated at each step of altitude z in a midlatitude atmosphere using equations ( 1)-( 12).

Results

Prediction of column collapse

To compare the model predictions with the previous studies of Kaminski and Jaupart (2001) and Carazzo et al. (2008a), we define the column regimes as a function of source mass discharge rate and total gas content, which can both be retrieved from field data.

Fig. 3 shows the predictions of the model for column collapse when accounting for reduced entrainment and particle sedimentation without gas entrapment (i.e., 𝑥 ! ! = 𝑥 !"! ). The transition curve is highly influenced by the value of D: the critical mass discharge rate before collapse is shifted by about an order of magnitude between the coarsest (D = 2.5) and finest (D = 3.3) population. For mass discharge rates larger than 10 7 kg s -1 , low values of D tend to increase the critical mass discharge rate at which collapse occurs. In this case, the loss of particles by sedimentation decreases significantly the column mass discharge rate during its rise, but is not large enough to drain out the thermal reservoir available to heat up the entrained cold atmospheric air, which helps the generation of a buoyant plume. We thus obtain that when sedimentation only is taken into account, coarse distributions (i.e. low values of D) promote the formation of stable plumes, a result apparently at odds with the conclusions of Kaminski and Jaupart (2001). To settle these conflicting results, we study the net effect of sedimentation when gas entrapment is also considered.

Fig. 4 shows the combined effect of particle sedimentation and gas entrapment compared to the case considering gas entrapment but no sedimentation. We consider here that the exponent D controls the effective "free" gas content (Fig. 2a) and that there is no open porosity (𝜉 = 0). We find that gas entrapment, hence the characteristics of the population of particles produced by fragmentation, has in general a dominant role on column collapse compared to particle sedimentation. However, for D values smaller than 2.8, sedimentation and gas entrapment are of equal importance and act together to hamper and even make impossible the formation of stable plumes. This result is consistent with the observations of Kaminski and Jaupart (1998) who compiled values of D systematically larger than 3.0 in all the pumice fallout deposits they considered. The shift between the two transition curves calculated at D = 2.8 (Fig. 4) suggests that, for small values of D, gas entrapment results in the decrease of the vertical velocity near the vent, which further enhances the effect of particle sedimentation.

The model accounting for gas entrapment and a constant entrainment presented by Kaminski and Jaupart (1998), and the model described by Carazzo et al., (2008a) accounting for a variable entrainment but no gas entrapment, both make the formation of a buoyant column less likely. The combined effects of sedimentation, gas entrapment and reduced entrainment lead to the drastic effect on the transition described here (Fig. 4) and casts some doubt on their ability to reproduce natural data previously explained by each model considered separately. These two models are however end-members and the introduction of open porosity may yield a more balanced conclusion.

Predictions for the dynamics of collapsing fountains

Girault et al. (2014) showed that the power-law exponent of the TGSD at the vent reduces the maximum height reached by a stable plume by 30% for mass discharge rates larger than 10 7 kg s -1 . To investigate a similar effect in the case of collapsing volcanic fountains, we performed calculations for a power-law exponent D ranging from 2.5 to 3.3 and for two different values of total exsolved gas content of 2 and 4 wt% (Fig. 6). For each value of gas content, we tested two different initial MDR corresponding to conditions near the plume/fountain transition (i.e., 3.10 8 kg s -1 for 𝑥 ! ! = 2 wt%, and 10 10 kg s -1 for 𝑥 ! ! = 4 wt%) and far from it (i.e., 10 10 kg s -1 for 𝑥 ! ! = 2 wt%, and 10 11 kg s -1 for 𝑥 ! ! = 4 wt%). We further consider an open porosity of 65%. We compare the maximum height reached by the fountain when accounting for particle sedimentation (Sed) and/or gas entrapment (GE) and/or open porosity (OP) or none of these effects. Fig. 6 shows that when mainly composed of fine particles (i.e., high D value), the fountain will reach a relatively high altitude allowing the resulting PDC to cover larger distances on ground. On the other hand, when the fountain has a coarser distribution of fragments (i.e., low D value), it will reach a lower maximum height, and one can expect the resulting PDC to be restricted to closer distances from the vent due to the smaller potential energy available. [START_REF] Bursik | The dynamics and thermodynamics of large ash flows[END_REF] showed that varying the MDR by an order of magnitude changed the run-out distance reached by the PDC by a factor of ~2.4 (Fig. 12a of [START_REF] Bursik | The dynamics and thermodynamics of large ash flows[END_REF], an increase that can be explained by doubling the height of the column in the plume model [START_REF] Woods | The fluid dynamics and thermodynamics of eruption columns[END_REF]. In our model, a change in fountain height of a factor 1.5 is obtained when varying D from 2.5 to 3.3 (red curves in Fig. 6a, b, c, d),

suggesting in turn that the run-out distance will increase by a factor of ≈ 2 for this change of D.

Our model can also be used to infer the grain size distribution in the PDC generated by a collapsing column. The results are presented in Fig. 7 based on four particle classes: bombs (-9 to -7 𝜙), lapilli (-6 to -1 𝜙), coarse ash (0 to 6 𝜙), and fine ash (7 to 10 𝜙). For low values of D, the collapsing fountain is more concentrated and enriched in coarse fragments (bombs and lapilli) at its maximum height. These large fragments shall form dense PDC loaded with bombs and lapilli. Conversely, for high values of D, the collapsing fountain is mainly made of fine ash particles at its maximum height, which will produce dilute PDC mostly made of coarse and fine ash particles. Varying the initial gas content and the MDR feeding the fountain does not affect these conclusions (Fig. 7), suggesting that D has a major control on the population of particles in the PDC produced during column collapse.

Comparison with natural cases

We now compare our model predictions with historical eruptions providing well-constrained values of gas contents, MDR and power-law exponent D. There are however only a few examples of such well-documented events, and here we choose two famous historical ones:

the ≈186 CE Taupo and the 79 CE Vesuvius eruptions. All the available geological data are summarized in Table 1.

The ≈186 CE Taupo eruption

The Taupo eruptive episode consists in a series of Plinian and phreato-Plinian eruptions [START_REF] Walker | The Taupo pumice: product of the most powerful known (ultraplinian) eruption[END_REF][START_REF] Froggatt | Stratigraphy and nature of Taupo Pumice Formation[END_REF][START_REF] Talbot | Dilute gravity current and rain-flushed ash deposits in the 1.8 ka Hatepe Plinian deposit[END_REF][START_REF] Wilson | The Taupo eruption, New Zealand II. The Taupo ignimbrite[END_REF][START_REF] Wilson | Stratigraphy, chronology, styles and dynamics of late Quaternary eruptions from Taupo volcano, New Zealand[END_REF][START_REF] Wilson | The Taupo eruption, New Zealand. I General aspects[END_REF][START_REF] Houghton | The 1800a Taupo eruption: "III wind" blows the ultraplinian type event down to Plinian[END_REF]. The stratigraphic sequence starts with a phreato-Plinian ash layer (Y1) overlain by a Plinian fallout (Y2 or Hatepe), and two phreato-Plinian ash layers (Y3 and Y4 or Rotongaio). The main sequence corresponds to a Plinian fallout deposit (Y5 or Taupo Plinian pumice) interbedded with intraplinian ignimbrites (Y6 or early ignimbrite flow units) and covered by PDC deposits (Y7 or Taupo ignimbrite). Here, we focus on the Y2, Y5-6 and Y7 phases as they represent the transition between Plinian, transitional, and collapse regimes, respectively. During these phases, the MDR increases from 1.8×10 ! kg s -1 during Y2 [START_REF] Carey | The intensity of Plinian eruptions[END_REF] to 2×10 ! kg s -1 during Y5/6 [START_REF] Carey | The intensity of Plinian eruptions[END_REF][START_REF] Houghton | The 1800a Taupo eruption: "III wind" blows the ultraplinian type event down to Plinian[END_REF][START_REF] Carazzo | The timing and intensity of column collapse during explosive volcanic eruptions[END_REF], and to 1.4×10 !" kg s -1 during Y7 [START_REF] Carey | The intensity of Plinian eruptions[END_REF][START_REF] Bursik | The dynamics and thermodynamics of large ash flows[END_REF]. These values have been debated in the literature.

Here we retain average values for Y2, Y5/6 and Y7, and we use the minimum and maximum MDR proposed in the literature to estimate an error bar in Fig. 8 (see Table 1).

Dunbar and Kyle (1993) estimated the dissolved water content to 4.3 ± 0.5 wt% for the Y2 and Y5-6 phases and to 3.6 ± 0.5 wt% for the last collapsing phase Y7. We correct these values for the presence of crystals and lithic fragments, which do not contain volatiles (Kaminski and Jaupart, 2001). Considering a proportion of crystals and lithics contents of 20% (Dunbar et al. 1989) and assuming complete degassing, we deduce a total gas mass fraction in the mixture of 3.44 ± 0.4 wt% for both the Y2 and Y5-6 phases and 2.88 ± 0.4 wt% for the Y7 phase. The power-law exponent D is taken to be 3.2, as calculated by Kaminski and Jaupart (1998), the magma temperature is set to 1133 K, as measured by Shane (1998), and the exit velocity is set to 306 ± 10 m s -1 , as determined from the isopleth maps of the transitional phase Y5-6 [START_REF] Walker | The Taupo pumice: product of the most powerful known (ultraplinian) eruption[END_REF]. We note that using the power-law exponent D =

and an open porosity of 60-70% gives

! ! ! ! !"!
≈ 0.94 (Fig. 2b) hence a small net effect of gas entrapment. Together with the total gas mass fraction of 3.44 wt%, this value gives a calculated exit velocity of 275 ± 5 m s -1 for the Y5-6 phase (Appendix B), which is consistent with the exit velocity constrained from the distribution of lithic fragments.

Fig. 8 compares the eruptive parameters inferred for the Taupo eruption (i.e., mass discharge rate and total gas content) and the transition between the stable plumes and the collapsing fountains predicted when considering particle sedimentation, gas entrapment and an open porosity of 60-70%. The good agreement between the model predictions and the field data illustrates the global consistency of the model.

The 79 CE Vesuvius eruption

The 79 CE Vesuvius eruption is one of the most well-documented historical Plinian event [START_REF] Lirer | Two Plinian pumice-fall deposits from Somma-Vesuvius, Italy[END_REF][START_REF] Sigurdsson | The eruption of Vesuvius in AD 79: reconstruction from historical and volcanological evidence[END_REF][START_REF] Wilson | The Taupo eruption, New Zealand II. The Taupo ignimbrite[END_REF][START_REF] Sigurdsson | Influence of magma withdrawal on compositional gradients during the AD 79 Vesuvius eruption[END_REF][START_REF] Carey | Temporal variations in column height and magma discharge rate during the 79 A.D. eruption of Vesuvius[END_REF][START_REF] Cioni | Volatile content and degassing processes in the AD 79 magma chamber at Vesuvius (Italy)[END_REF][START_REF] Cioni | Compositional layering and syn-eruptive mixing of a periodically refilled shallow magma chamber: the AD79 Plinian eruption of Vesuvius[END_REF]2008;[START_REF] Balcone-Boissard | Textural and geochemical constraints on eruptive style of the 79AD eruption at Vesuvius[END_REF]Shea et al., 2012 and references therein) and allows a thorough comparison with our theoretical predictions. The eruption started with a short phreatomagmatic phase and then pursued with a paroxysmal

Plinian phase for about 17 hours. First, a buoyant stable column rose up to ≈27 km high and produced a massive layer of white pumice (WP) fallout deposits that can be decomposed in 4 chronostratigraphic levels (levels 0.25W, 0.5W, 0.75W and 1W of the white fallout sequence; [START_REF] Carey | Temporal variations in column height and magma discharge rate during the 79 A.D. eruption of Vesuvius[END_REF]. This phase is also called EU2 in the literature [START_REF] Cioni | Compositional layering and syn-eruptive mixing of a periodically refilled shallow magma chamber: the AD79 Plinian eruption of Vesuvius[END_REF]2008;[START_REF] Balcone-Boissard | Textural and geochemical constraints on eruptive style of the 79AD eruption at Vesuvius[END_REF][START_REF] Shea | Transitions between fall phases and pyroclastic density currents during the AD79 eruption at Vesuvius: building a transient conduit model from the textural and volatile record[END_REF]. After a first partial collapse, the magma composition changed and yielded more unstable conditions. The column reached a maximum height of ≈33 km and generated several grey pumice (GP) fallout deposits interbedded with surges (levels W/G-S1, S1-S2, S2-S3, S3-S4 of the grey fallout sequence; [START_REF] Carey | Temporal variations in column height and magma discharge rate during the 79 A.D. eruption of Vesuvius[END_REF]. This transitional phase is also called EU3 in the literature [START_REF] Cioni | Compositional layering and syn-eruptive mixing of a periodically refilled shallow magma chamber: the AD79 Plinian eruption of Vesuvius[END_REF]2008;[START_REF] Balcone-Boissard | Textural and geochemical constraints on eruptive style of the 79AD eruption at Vesuvius[END_REF][START_REF] Shea | Transitions between fall phases and pyroclastic density currents during the AD79 eruption at Vesuvius: building a transient conduit model from the textural and volatile record[END_REF]. A first total collapse ended this magmatic phase. The eruption pursued for another few hours with a last short-lived plume before a second total collapse marking the end of the eruption [START_REF] Balcone-Boissard | Textural and geochemical constraints on eruptive style of the 79AD eruption at Vesuvius[END_REF].

The mass eruption rates of successive phases were calculated, from the column heights inferred from the isopleth maps and the temperate discharge curves of [START_REF] Sparks | The dimensions and dynamics of volcanic eruption columns[END_REF]: the MDR increased from 6.1×10 ! to 7.1×10 ! kg s -1 for the WP fall sequence and then up to 1.5×10 ! kg s -1 for the first partial collapse involving GP [START_REF] Carey | Temporal variations in column height and magma discharge rate during the 79 A.D. eruption of Vesuvius[END_REF]Kaminski and Jaupart, 2001). We retain an uncertainty of a factor of 2 for the MDR, which corresponds to the change in MDR required to reach the same column height when using the model of [START_REF] Sparks | The dimensions and dynamics of volcanic eruption columns[END_REF] or a more recent one [START_REF] Girault | The effect of total grainsize distribution on the dynamics of turbulent volcanic plumes[END_REF]. The total dissolved gas content of the melt is estimated to 6.00 ± 0.44 wt% at the beginning of the Plinian stable phase, and then decreased to 4.95 ± 0.20 wt% at the WP/GP boundary [START_REF] Cioni | Volatile content and degassing processes in the AD 79 magma chamber at Vesuvius (Italy)[END_REF]. Using the same method as for the Taupo eruption, we correct these values for the presence of crystals and lithic fragments. The crystals content increased from 20% during the WP phase to 28% during the GP phase [START_REF] Shea | Leucite crystals: surviving witnesses of magmatic processes preceding the 79 AD eruption at Vesuvius, Italy[END_REF], and the lithics content increased from 12%

during the WP phase to 20% during the GP phase [START_REF] Sigurdsson | The eruption of Vesuvius in AD 79: reconstruction from historical and volcanological evidence[END_REF]. Assuming complete degassing, we find that the total gas mass fraction in the magma decreased from 4.08 ± 0.30 wt% during the WP phase to 2.57 ± 0.10 wt% during the GP phase. The power law exponent D is taken to be 3.0 as measured by Costa et al. (2016b), the magma temperature is set to 1323 K, as measured by [START_REF] Gurioli | Complex changes in eruption dynamics during the 79 AD eruption of Vesuvius[END_REF] for the GP phase, and the exit velocity is set to 260 ± 10 m s -1 for the GP phase from the isopleth maps of [START_REF] Carey | Temporal variations in column height and magma discharge rate during the 79 A.D. eruption of Vesuvius[END_REF]. We note that using the power-law exponent D = 3.0 and an open porosity of 60-70% gives

! ! ! ! !"!
≈ 0.84 (Fig. 2b). Together with the total gas mass fraction of 2.57 wt%, this value gives a calculated exit velocity of 230 ± 5 m s -1 for the W/G-S1 phase (Appendix B), which is consistent with the exit velocity constrained from the distribution of lithic fragments.

Fig. 9 compares the eruptive parameters inferred for the Vesuvius eruption (i.e., mass discharge rate and total gas content) and the transition between the stable plumes and the collapsing fountains predicted when considering particle sedimentation, gas entrapment and an open porosity of 60-70%. The model fails to reproduce the conditions of the transition between the Plinian stable regime of the WP fallout sequence and the collapse regime at the beginning of the GP fallout sequence. As all the model parameters are set by the field constraints, they cannot be tuned to improve the model predictions. However, we note that field deposits from this eruption are characterized by a power-law exponent D = 3.0, a D value that defines the transition between a population dominated by coarser fragments (i.e., D

< 3) and a population dominated by ash (i.e., D > 3). We suggest that the hypothesis of thermal equilibrium has to be reconsidered for D = 3.0. For D > 3, the ash particles dominating the population are small enough to remain in thermal equilibrium [START_REF] Woods | Particle fallout, thermal disequilibrium and volcanic plumes[END_REF]. For D < 3, coarse pumices rapidly settle to the ground leaving fine and coarse ash particles in thermal equilibrium in the column. For D = 3.0, there is no dominant population: the average size of the fragments is intermediate between fine ash and lapilli fragments, and is too large to insure perfect thermal equilibrium. In that case column collapse in favored and the critical mass discharge rate before collapse can be reduced by up to an order of magnitude (Fig. 5 of [START_REF] Woods | Particle fallout, thermal disequilibrium and volcanic plumes[END_REF], which would shift the theoretical predictions towards a better agreement with field data from the 79 CE Vesuvius eruption.

Discussion

The effect of crater shape on exit velocity

Gas entrapment -more or less modulated by the open porosity -reduces the amount of gas in the mixture and, in turn, the vertical velocity at the base of the plume. The presence of a crater with a specific shape can similarly reduce the velocity at the base of the plume [START_REF] Woods | The decompression of volcanic jets in a crater during explosive volcanic eruptions[END_REF][START_REF] Koyaguchi | Effects of the crater on eruption column dynamics[END_REF], and it is always possible to define a crater shape that will perfectly mimic the effect of gas entrapment. As field data provide only constraints on the exit velocity, it is not possible to decide which model should be chosen. Overall, both decompression in a crater and gas entrapment can occur during an eruption and act together to set the conditions for collapse. In the case of the Taupo eruption, the global consistency between our model predictions with gas entrapment and open porosity and field data can be taken as a strong argument in favor of this framework. Adding a crater constraining decompression in the case of Taupo will actually decrease the agreement between model predictions and the data. In the case of Vesuvius, our model does not predict correctly the transition, but the predictions would not be changed if a crater is added as the exit velocity will remain constrained by the lithics.

Effect of wind on column collapse

The presence of wind can affect the plume dynamics and stability [START_REF] Bursik | Effect of wind on the rise height of volcanic plumes[END_REF][START_REF] Degruyter | Impact of wind on the condition for column collapse of volcanic plumes[END_REF][START_REF] Woodhouse | Interaction between volcanic plumes and wind during the 2010 Eyjafjallajökull eruption, Iceland[END_REF][START_REF] Suzuki | 3D numerical simulation of volcanic eruption clouds during the 2011 Shinmoe-dake eruptions[END_REF]2015;[START_REF] Mastin | Testing the accuracy of a 1-D volcanic plume model in estimating mass eruption rate[END_REF][START_REF] Girault | Combined effects of total grain-size distribution and crosswind on the rise of eruptive volcanic columns[END_REF]. Thereotical and numerical studies show that for high MDR and/or low wind velocity, the volcanic column forms a strong plume not affected by the wind field [START_REF] Bonadonna | Sedimentation from strong volcanic plumes[END_REF]. For low MDR and/or high wind velocity, the volcanic plume trajectory is strongly controlled by the wind strength and direction, and the volcanic column forms a weak bended plume [START_REF] Bursik | Effect of wind on the rise height of volcanic plumes[END_REF]. The critical MDR before collapse is also strongly increased by a crosswind relative to no wind conditions, thus promoting the formation of a stable plume [START_REF] Degruyter | Impact of wind on the condition for column collapse of volcanic plumes[END_REF][START_REF] Girault | Combined effects of total grain-size distribution and crosswind on the rise of eruptive volcanic columns[END_REF]. [START_REF] Girault | Combined effects of total grain-size distribution and crosswind on the rise of eruptive volcanic columns[END_REF] showed that this effect is more pronounced for small MDR when D < 3 than for large MDR when D > 3, and that it strongly depends on the shape of the wind profile. The incorporation of crosswind in our model would thus requires to exactly know the shape of the wind profile, a parameter that is unknown for historical eruptions, and to parameterize the turbulent entrainment coefficient due to wind (commonly named 𝛽). The latter is currently unconstrained from laboratory experiments (Costa et al., 2016a), and its relation with the turbulent entrainment coefficient 𝛼 still remains unclear (Aubry et al., 2017a;2017b). We therefore consider that a systematic study of this effect remains outside the scope of the present study.

We however check the validity of our no wind model when applied to the two historical eruptions considered here, using the method of [START_REF] Carazzo | Laboratory experiments of forced plumes in a density-stratified crossflow and implications for volcanic plumes[END_REF] to estimate the impact of wind on the plume dynamics. According to the isopleths maps determined from field data, the maximum wind velocity reached about 28 m s -1 during the WP phase of the 79 CE Vesuvius eruption, 31 m s -1 during the GP phase [START_REF] Carey | Temporal variations in column height and magma discharge rate during the 79 A.D. eruption of Vesuvius[END_REF], and 27 m s -1 during the ≈186 CE Taupo eruption [START_REF] Carey | Quantitative models of fallout and dispersal of tephra from volcanic eruption columns[END_REF]. Together with the source conditions given in Table 1, we calculate the ratio of the dimensionless wind velocity to the dimensionless plume velocity 𝑊 * /𝑈 * which is used as a proxy to infer the plume regime [START_REF] Carazzo | Laboratory experiments of forced plumes in a density-stratified crossflow and implications for volcanic plumes[END_REF]. The 79 CE Vesuvius column is found to be a slightly distorted plume at the beginning of the eruption (𝑊 * /𝑈 * ≈ 0.2) and a strong plume before the column collapses (𝑊 * /𝑈 * ≈ 0.1). The ≈186 CE Taupo column is found to be a strong plume during the entire eruption (𝑊 * /𝑈 * ≈ 0.06). For these two historical eruptions, atmospheric winds did not affect significantly the plume dynamics, and can thus be neglected when investigating the mechanisms leading to column collapse.

Conclusions

We have tested the combined effects of gas entrapment and reduced entrainment using a 1D turbulent plume model accounting for particle sedimentation. Compared with the predictions of previous studies accounting for gas entrapment or reduced entrainment only, this model favors too much column collapse compared with natural cases. We thus take into account open porosity, i.e. the fraction of bubbles preserved inside a fragment after fragmentation but connected to the exterior, which helps to generate stable plumes. We then predict various grain-size distributions in PDC depending on the TGSD (hence D value) at the volcanic vent.

Low D values tend to promote the formation of relatively small fountains producing concentrated PDC with large particles reaching relatively low run-out distances. High D values tend to promote the formation of higher fountains producing dilute PDC containing mostly fine particles and covering larger distances. Finally, we compare the model predictions to the two well-known eruptions of Taupo (≈186 CE) and Vesuvius (79 CE). We find a good agreement between theoretical predictions and field data for the Taupo eruption. For the Vesuvius eruption characterized by a power-law exponent D = 3.0, we suggest that thermal disequilibrium could explain the discrepancy between our predictions and field data. where 𝑓 is the frictional forces acting on the flow, P the flow pressure in the conduit, and 𝑔 the gravitational acceleration.

Before exsolution, we consider that the liquid can be modelled as Poiseuille flow,

where the mass discharge rate 𝑄 is linked to the conduit radius 𝑎 by

𝑄 = 𝜌 ! 𝜋𝑎 ! Δ𝜌𝑔 8𝜇 ! , ( 𝐵3
)
where 𝜌 ! and 𝜇 ! are the liquid density and viscosity taken as 2500 kg m -3 and 10 5 Pa s, respectively, and Δ𝜌 is the density difference between the liquid and the country rock, which we set at 100 kg m -3 . Combining equations (B1) and (B3) gives

𝑆 𝑄 = 8𝜋𝜇 ! 𝑄𝜌 ! Δ𝜌𝑔 !/! , (𝐵4)
which can be replaced in equation (12).

We now seek for an expression for 𝑈 ! and 𝑃 ! in equation (12). Above the fragmentation level, the particle-gas mixture density depends on the effective amount of free gas 𝑥 ! in the volcanic mixture, such as

𝜌 = !!! ! ! ! + ! ! ! ! ! ! !! , (𝐵5)
where we assume that 𝜌 ! and T remain constant in the conduit [START_REF] Woods | The decompression of volcanic jets in a crater during explosive volcanic eruptions[END_REF].

By replacing equations (B1) and (B5) into equation (B2), one can write

𝑑𝑃 𝑑𝑧 1 -𝑢 ! 𝑥 ! 𝑅 ! 𝑇 𝜌 ! 𝑃 ! = -𝜌𝑔 -𝑓, (𝐵6)
which can be used to find the velocity at the vent 𝑈 ! by setting the coefficient on the left-hand side to zero, which yields

𝑈 ! = 𝑃 ! 𝜌 ! 𝑥 ! 𝑅 ! 𝑇 !!/! , (𝐵7) 𝑃 ! = 𝑄 𝑆 𝑥 ! 𝑅 ! 𝑇 !/! , (𝐵8)
which can be by using equation (B4) as:

𝑃 ! = 𝑄𝑥 ! 𝑅 ! 𝑇𝜌 ! Δ𝜌𝑔 8𝜋𝜇 ! !/! . ( 𝐵9
)
From equations (B4), (B5), (B7), (B9) together with equation ( 12) one can then calculate 𝑈 ! , the vertical velocity of the volcanic column after decompression. [START_REF] Carey | The intensity of Plinian eruptions[END_REF] and MDR calculated from the maximum height proposed by [START_REF] Houghton | The 1800a Taupo eruption: "III wind" blows the ultraplinian type event down to Plinian[END_REF]; the MDR feeding the Y6 phase is calculated from the MDR of Y5 and from masses of deposits for Y5 and Y6 given in [START_REF] Carazzo | The timing and intensity of column collapse during explosive volcanic eruptions[END_REF]; the MDR feeding the Y5-6 phase is the addition of the two precedent values.

Tables

** Average value from MDR proposed by [START_REF] Carey | The intensity of Plinian eruptions[END_REF] and MDR proposed by [START_REF] Bursik | The dynamics and thermodynamics of large ash flows[END_REF]. The curve calculated without sedimentation displayed in Fig. 4 of Carazzo et al. (2008a) is the same as the one calculated with D = 3.3. We consider a mid-latitude atmosphere, the magma temperature is taken as the average of andesitic magma (T 0 = 1200K), and only sedimentation is introduced in the model compared to Carazzo et al. (2008a). and dashed curves represent that from our model considering gas entrapment only (as in Kaminski and Jaupart, 2001). We consider a mid-latitude atmosphere, and the magma temperature is taken as the average of andesitic magma (T 0 = 1200K). ) and total gas content (in wt%). The curves give the threshold mass discharge rate for column collapse when accounting for sedimentation and gas entrapment. Calculations are made for a power-law exponent D = 3.2 calculated from field deposits (Kaminski and Jaupart, 1998), an initial magma temperature of 1133 K [START_REF] Shane | Correlation of rhyolitic pyroclastic eruptive units from the Taupo volcanic zone by Fe-Ti oxide compositional data[END_REF], and an exit velocity of 306 m s -1 [START_REF] Walker | The Taupo pumice: product of the most powerful known (ultraplinian) eruption[END_REF]. The open circle stands for the stable Plinian phase Y2, the grey square for the two stages of transitional regime Y5-Y6, and the black square for the final total collapse Y7 (Table 1). an initial magma temperature of 1323 K [START_REF] Gurioli | Complex changes in eruption dynamics during the 79 AD eruption of Vesuvius[END_REF], and an exit velocity of 260 m s - 1 [START_REF] Carey | Temporal variations in column height and magma discharge rate during the 79 A.D. eruption of Vesuvius[END_REF]. The open circles stand for stable Plinian phases (white pumices, WP or EU2) and the grey square for the first stage of transitional regime (grey pumice, GP or EU3) (Table 1).

Fig. 5

 5 Fig.5shows the combined effect of particle sedimentation, gas entrapment and post-

  Fig.6shows that particle sedimentation alone has a negligible (yet positive) effect on the maximum fountain height. When gas entrapment alone or gas entrapment plus sedimentation are taken into account, the maximum fountain height strongly decreases when D decreases, in agreement with the results of[START_REF] Girault | The effect of total grainsize distribution on the dynamics of turbulent volcanic plumes[END_REF], which emphasizes the dominant role of gas entrapment. The change in fountain height can reach up to a factor of 5 for conditions far from the transition (Fig.6a,b), to a factor of 15 for conditions near the transition (Fig.6 c,d). Finally, when an open porosity of 65% is taken into account (Sed+GE+OP), the effect of gas entrapment is largely reduced and intermediate fountain heights are obtained (red curves in Fig.6).

C

  shear stress parameter, 0.135 𝑐 ! specific heat at constant pressure of the atmosphere, 998 J K -1 𝑐 ! specific heat at constant pressure of the particles, 1100 J K -1 𝐷 exponent of the power-law distribution of particles, dimensionless 𝑑 ! particle diameter, m f frictional forces per unit of volume acting on the flow, N m -3 𝑔 gravitational acceleration, 9.81 m s -2 n total volatile content of the melt, wt% 𝑁 ! number of classes of particle size, 20 P flow pressure in the conduit, Pa 𝑃 ! atmospheric pressure, Pa 𝑃 ! flow pressure at the vent, Pa 𝑝 ! probability of sedimentation, 0.27 ± 0.01 Q mass discharge rate feeding the eruption, kg s -1 𝑄 ! mass discharge rate of particles, kg s -1 r particle radius, m 𝑟 !"# maximum particle size considered in this study, m 𝑟 !"# minimum particle size considered in this study, m 𝑅 column radius, m 𝑅 ! column radius at the vent, m 𝑅 ! gas constant of the atmospheric air, 287 J K -1 kg -1 𝑅𝑒 ! Reynolds number of the particles, dimensionless 𝑅 ! bulk constant of the mixture in the column, 𝑅 ! ! = 461 J K -1 kg -1 𝑅𝑖 local Richardson number, dimensionless S cross-sectional area of the conduit, m 2 𝑇 flow temperature, K 𝑇 ! ambient atmospheric temperature for a mid-latitude atmosphere, 273 K 𝑇 ! average andesitic magma temperature, 1200 K u velocity of the magma-volatiles mixture, m s -1 𝑈 average vertical velocity in the column, m s -1 𝑈 ! column average vertical velocity after decompression, m s -1 𝑈 ! entrainment rate at the edge of the plume, m s -1 𝑈 ! average vertical velocity at the vent, m s -1 𝑉 !"#$%&&'( , volume of gas released due to the development of open porosity, m 3 𝑉 !"#$$%& , volume of gas initially trapped at fragmentation, m 3 𝑉 ! particle fallout velocity, m s -1 𝑥 ! effective gas mass fraction, dimensionless 𝑥 !"! total gas content, dimensionless 𝑧 vertical axis, m z*, dimensionless height, z* = 𝑧/(2𝑅 ! ) 𝛼 ! variable entrainment coefficient, dimensionless 𝜆 normalization constant 𝜇 ! dynamic viscosity of air at the vent, 𝜇 ! = 1.832 ×10 !! Pa s 𝜇 ! liquid viscosity, 10 5 Pa s 𝜙 particle size notation, 𝑑 ! = 2 !! 𝜌 bulk density of the volcanic mixture, kg m -3 𝜌 ! atmospheric density, kg m -3 𝜌 ! liquid density, 2500 kg m -3 𝜌 ! particle density, 2000 kg m -3 𝜌 ! flow density at the vent, kg m -3 Δ𝜌 density difference between the liquid and the country rock, 100 kg m -3 𝜉 fraction of gas initially trapped and released by open porosity, dimensionless 𝜌𝑢

1 FiguresFig. 1 .

 11 Figures

Fig. 2 .

 2 Fig. 2. (a) Fraction of gas released at fragmentation as a function of the power law exponent

Fig. 3 .

 3 Fig. 3. Threshold mass discharge rate (in kg s -1 ) at the transition between buoyant and

Fig. 4 .

 4 Fig. 4. Threshold mass discharge rate (in kg s -1 ) at the transition between buoyant and

Fig. 5 .

 5 Fig. 5. Threshold mass discharge rate (in kg s -1 ) at the transition between buoyant and

Fig. 9 .

 9 Fig. 9. Eruptive conditions of the 79 CE Vesuvius eruption in terms of source mass discharge rate (in kg s -1 ) and gas content (in wt%). The curves give the threshold mass discharge rate for

Table 1 .

 1 Eruption conditions of the 79 CE Vesuvius and ≈186 CE Taupo events. The MDR feeding the Y5 phase is an average value from MDR proposed by

	Estimated	minimal vent	velocity	(calculated	from isopleth	maps in	literature, m s -1 )	140 +/-10 (10)	or 306 +/-10	(15)	/	/	/
			Plume	maximum	height (km)			31 -37 (10),	50 (15)	/	/	/
			Magma	temperature,	T0 (K)
				Power-law	exponent D
				Mass Discharge Rate,	MDR (kg s -1 )				3.3 x10 8 -1.1 x10 9 * (9,	10)	5.8 x10 8 -1.9 x10 9 * (9,	10, 11)	9 x10 8 -3 x10 9 * (9, 10,	11)	7.7 x10 9 ** (9, -2 x10 10	12)
		Total crystal	content and	lithics content	(%)
		Total gas content,	xtot (corrected for	phenocrysts and	lithics, wt%)					2.9 +/-0.4
				Dissolved gas	content, n0 (wt%)							3.6 +/-0.5 (7)
					Eruption (deposit)					Taupo 186 CE	(Y5)	Taupo 186 CE	(Y6)	Taupo 186 CE	(Y5-6)	Taupo 186 CE	(Y7)	*

Acknowledgments

The authors warmly thank the reviewers, T. Koyaguchi and W. Degruyter, and the editor, J.

Marti for their constructive comments which greatly improved the original version of this manuscript. We are grateful to H. Balcone-Boissard for fruitful discussions on porosity of pumices from the Vesuvius eruption. This work is a part of A. Michaud-Dubuy's PhD thesis (Institut de Physique du Globe de Paris). Support was provided by RAVEX (ANR contract ANR-16-CE03-0002). This is IPGP contribution number 3944.

Appendix A

The model of turbulent entrainment of [START_REF] Kaminski | Turbulent entrainment in jets with arbitrary buoyancy[END_REF] introduced two dimensionless variables: C, which gives the fraction of kinetic energy available for turbulent entrainment, and A which depends on the half width of velocity and density profiles in the jet. [START_REF] Carazzo | The route to self-similarity in turbulent jets and plumes[END_REF]2008a;2008b) further showed that A is not a constant but evolves as a function of the distance from the source and buoyancy. They propose the following empirical formula:

where A j and A p are the values of A for a pure jet (Ri=0), and a pure plume (no mass discharge rate at the source), respectively. Lm is the Fisher length scale, defined as:

where M 0 and F 0 are the momentum and buoyancy flow rates at the source, respectively.

These parameters were estimated from laboratory measurements [START_REF] Carazzo | The route to self-similarity in turbulent jets and plumes[END_REF]) and can be fitted as a function of the dimensionless height = z /(2R 0 ) by:

and

This set of equations is injected into Eq. ( 4) to calculate the entrainment coefficient 𝛼 ! at each altitude step z.

Appendix B

We derive an expression for the exit velocity after decompression following [START_REF] Woods | The decompression of volcanic jets in a crater during explosive volcanic eruptions[END_REF]. For an impermeable conduit, the mass discharge rate Q remains constant with depth, so that the mass conservation equation gives:

where 𝜌 and 𝑢 are the density and the velocity of the magma-volatiles mixture, respectively, and 𝑆 = 𝜋𝑎 ! is the cross-sectional area of the conduit with 𝑎 its radius. For simplicity, we assume that the conduit has a constant cross-section.

The momentum conservation equation in the conduit can be written as [START_REF] Wilson | Explosive volcanic eruptions-IV. The control of magma properties and conduit geometry on eruption column behaviour[END_REF][START_REF] Woods | The decompression of volcanic jets in a crater during explosive volcanic eruptions[END_REF][START_REF] Koyaguchi | An analytical study for 1-dimensional steady flow in volcanic conduits[END_REF] sedimentation, gas entrapment and an open porosity of (a) 𝜉 = 20%, (b) 𝜉 = 65%, and (c) 𝜉 = 80%. We consider a mid-latitude atmosphere, and the magma temperature is taken as the average of andesitic magma (T 0 = 1200K). effects. The green curve corresponds to the model of Kaminski and Jaupart (1998;[START_REF] Bursik | Effect of wind on the rise height of volcanic plumes[END_REF], the dashed curve to the model of Carazzo et al. (2008a) and the blue curve to the model of [START_REF] Girault | The effect of total grainsize distribution on the dynamics of turbulent volcanic plumes[END_REF]. Calculations are made at conditions far above the transition curves with an initial exsolved gas content and an initial mass discharge rate of (a) 2 wt% and 10 10 kg s -1 , and (b) 4 wt% and 10 11 kg s -1 , and just above the transition curves for initial values of (c) 2 wt% and 3.10 8 kg s -1 , and (d) 4 wt% and 10 10 kg s -1 . We consider a mid-latitude atmosphere, and the magma temperature is taken as the average of andesitic magma (T 0 = 1200K) in all calculations.


Andrea D'ariano
email: dariano@ing.uniroma3.it

Francesco Corman

Taku Fujiyama

Lingyun Meng

Paola Pellegrini

Andrea D ' Ariano

Simulation and Optimization for Railway Operations Management

In the forthcoming decades, railway transport is expected to face a significant growth of traffic flows which will mostly have to be accommodated over the existing infrastructures. An increase of capacity utilization is needed to avoid reduction of reliability and punctuality of transport services. Furthermore, energy-efficient operations and multimodality opportunities are also topics of growing interest.

Intelligent railway operations management requires accurate modelling and simulation of train and pedestrian traffic flows and optimal management of key decisions at strategic, tactical, and operational levels. This special issue focuses on advanced mathematical modelling and optimal control in the railway domain and related multimodal transport networks. We aim to identify new methods for improving the effectiveness and efficiency of railway operations, including the development of advanced algorithmic techniques for timetabling, capacity management, infrastructure management, and traffic and passengers flow management.

This special issue has selected a compendium of research papers addressing recent theoretical and practical advances on railway operations management. 32 papers were submitted to this special issue, 14 of which were accepted for publication. As the guest editors of this special issue, we next summarize the 14 accepted papers.

In the paper titled "A Short Turning Strategy for Train Scheduling Optimization in an Urban Rail Transit Line: The Case of Beijing Subway Line 4" by M. Zhang et al., the authors propose a mixed-integer nonlinear programming model for the train scheduling with consideration of short turning strategy and train circulation plan. Mixed-integer linear programming approach is used to solve this optimization problem. Two case studies are carried out based on the data of Beijing subway line 4. The simulation results show that operation pattern with short turning train services can acquire a better train schedule to meet passenger demand and a better train circulation plan.

In the paper titled "Integrated Optimization on Train Control and Timetable to Minimize Net Energy Consumption of Metro Lines," Y. Zhou et al. present an integrated optimization model on train control and timetable to minimize the net energy consumption with consideration of utilizing regenerative energy. An improved model and algorithms on train control are proposed to attain energy-efficient speed profiles. Case studies on Beijing Metro Line 5 illustrate that the improved train control approach can save traction energy consumption by 20% in comparison with the commonly adopted train control sequence in timetable optimization.

In the paper titled "Application of Data Clustering to Railway Delay Pattern Recognition" by F. Cerreto et al., the authors employ K-means clustering to identify recurrent delay patterns on a high traffic railway line north of Copenhagen, Denmark. The clusters identify behavioral patterns in the very large ("big data") data sets generated automatically and continuously by the railway signal system. The results reveal where corrective actions are necessary, showing where recurrent delay patterns take place. The demonstrated methodology is scalable and can be potentially transferred to any system of transport.

In the paper titled "Defining Reserve Times for Metro Systems: An Analytical Approach" by L. D' Acierno et al., the authors provide an analytical approach for determining operational parameters for metro systems so as to support the planning and implementation of energy-saving strategies. They develop a suitable methodology for estimating reserve times that represent the main rate of extra time needed to put eco-driving strategies in place. The approach proposed by the authors is applied in Line 1 of the Naples metro system, whose service frequency was duly taken into account, to analyze operation configurations and to quantify the amount of saved energy.

In the paper titled "Using Smart Card Data Trimmed by Train Schedule to Analyze Metro Passenger Route Choice with Synchronous Clustering," W. Li et al. analyze smart card data in Shanghai metro systems to understand mobility patterns. To do that, they cluster travelers and smart card data transactions by looking at the pure travel time. This is then appropriately converted to a model of route choice throughout the network. They found out that those steps can improve substantially the amount of insight into the travelers' behavior.

In "A Simulation Platform for Combined Rail/Road Transport in Multiyards Intermodal Terminals," X. Chen et al. tackle multimodal freight terminals where the combination between road and rail transport is investigated by means of simulation in terminals featuring multiple rail yards. The key features of the proposed approach based on Time Petri Nets are an increased realism of train operations and container movements in Qianchang railway terminal, including train routing dispatching rules. The validation on historical data allows evaluating influence of design parameters on container operations.

In the paper titled "On Individual Repositioning Distance along Platform during Train Waiting," F. Leurent and X. Xie propose a stochastic model in which user's journey is decomposed into phases of, successively, walking in the access station, platform positioning, waiting for boarding, train riding, and walking in the egress station. Egress times and exit instants are random variables that are characterized by distribution and mass probability functions under closed form for both single and distributed walking speed. Maximum likelihood estimation is proposed and applied to a case study of commuter rail line in Paris, France, with satisfactory numerical results.

In the paper titled "Efficiency Assessment of Transit-Oriented Development by Data Envelopment Analysis: Case Study on the Den-en Toshi Line in Japan," J. Guo et al. assess the efficiency of Transit-Oriented Development (TOD), that is, an urban planning approach that encourages a modal shift from private to public transportation, by applying the Data Envelopment Analysis (DEA) method. The ridership of public transportation is considered as the direct output characteristic of TOD efficiency, and nine indicators of ridership are selected as inputs of TOD. The Tokyo Den-en Toshi Line in Japan is investigated as a typical case of TOD.

In the paper titled "Stop Plan of Express and Local Train for Regional Rail Transit Line" by Q. Luo et al., the Logit model is used to analyze the behavior of passengers choosing trains by considering the sensitivity of travel time and travel distance. Based on the composition of passenger travel time, an integer programming optimization model for train stop scheme is proposed, which aims at minimizing the total passenger travel time for a certain regional rail line in Shenzhen. Genetic Algorithms are used to solve the problem. The simulation result shows the feasibility of the proposed model and the efficiency of the proposed algorithm.

In the paper titled "Fuzzy Approach in Rail Track Degradation Prediction" by M. Karimpour et al., an adaptive network-based fuzzy inference system (ANFIS) model is proposed to estimate rail track degradation for the curves and straight sections of Melbourne tram track system. A fuzzy approach is proposed due to the nonlinear and noisy nature of the data according to the data that were available on the Melbourne tram network. Experimental results demonstrate that the developed model is capable of estimating the long-term behavior of rail tracks and predicting the gauge values with an R2 of 0.6 and 0.78 for curves and straights, respectively.

In the paper titled "PULSim: User-Based Adaptable Simulation Tool for Railway Planning and Operations," Y. Cui et al. introduce a user-based, customizable platform to provide the ability of defining sophisticated workflows for users. As the preconditions of the platform, the design aspects for modelling the components of a German railway system and building the workflow of railway simulation are elaborated. Based on the model and the workflow, an integrated simulation platform with open interfaces is developed. Users and researchers gain the ability to rapidly develop their own algorithms, supported by the tailored simulation process in a flexible manner.

In the paper titled "The Planners' Perspective on Train Timetable Errors in Sweden" by C.-W. Palmqvist et al., typical errors in train timetables of railways, relevant reasons, and potential benefits of new tools and processes are investigated and reported through reviewing the state of practice and the state of the art in timetable planning, studying the research literature and railway management documents from several European countries, and conducting interviews with timetable planners in Southern Sweden.

In the paper titled "Calculation Method for Load Capacity of Urban Rail Transit Station Considering Cascading Failure," J. Huang et al. quantify the capacity of urban rail transit stations in terms of number of passengers. Their paper relies on an association network of facilities, set up based on the analysis of passenger service chain in station. The analysis of a case study of Lujiabang Road Station in Shanghai shows that this algorithm can search for capacity bottlenecks and help trace the load variation of facilities in different scenarios, providing support for passenger flow organization.

In the paper titled "Optimal Operation of High-Speed Trains Using Hybrid Model Predictive Control" by Y. Yang et al., a control framework is proposed for the controller

design of an automatic train operation (ATO) system. Firstly, a piecewise linear system is constructed through piecewise linearization of the Chinese high-speed train's nonlinear dynamics. Secondly, the piecewise linear system is transformed into a mixed logical dynamical system. For the latter system, a hybrid model predictive controller is designed to realize the precise control. Simulation results validate the effectiveness of the framework.

Taken together, these 14 papers point to interesting research directions in mathematical modeling, simulation, and particularly optimization to tackle practically oriented railway operations management challenges in several countries worldwide, including China, Denmark, Italy, France, Japan, Australia, Germany, and Sweden.

Andrea D'Ariano Francesco Corman

Taku Fujiyama Lingyun Meng Paola Pellegrini