PLUME-MoM 1.0: A new integral model of volcanic plumes based on the method of moments

Abstract. In this paper a new integral mathematical model for volcanic plumes, named PlumeMoM, is presented. The model describes the steady-state dynamics of the plume in a 3D coordinate system, accounting for continuous variability in particle size distribution of the pyroclastic mixture ejected at the vent. Volcanic plumes are composed of pyroclastic particles of many different sizes ranging 5 from a few microns up to several centimeters and more. Proper description of such a multiparticle nature is crucial when quantifying changes in grain-size distribution along the plume and, therefore, for better characterization of source conditions of ash dispersal models. The new model is based on the method of moments, which allows description of the pyroclastic mixture dynamics not only in the spatial domain but also in the space of parameters of the continuous size-distribution of the 10 particles. This is achieved by formulation of fundamental transport equations for the multiparticle mixture with respect to the different moments of the grain-size distribution. Different formulations, in terms of the distribution of the particle number, as well as of the mass distribution expressed in terms of the Krumbein log scale, are also derived. Comparison between the new moments-based formulation and the classical approach, based on the discretization of the mixture in N discrete phases, 15 shows that the new model allows the same results to be obtained with a significantly lower computational cost (particularly when a large number of discrete phases is adopted). Application of the new model, coupled with uncertainty quantification and global sensitivity analyses, enables investigation of the response of four key output variables (mean and standard deviation of the grain-size distribution at the top of the plume, plume height and amount of mass lost by the plume during the ascent) to 20 changes in the main input parameters (mean and standard deviation) characterizing the pyroclastic mixture at the base of the plume. Results show that, for the range of parameters investigated and without considering interparticle processes such as aggregation or comminution, the, the grain-size

the second and third moment assume a spherical shape. Or if this is not the case, what do they assume? And what does the superscript i on Di mean? Line 20: it would be helpful to explain briefly what the Sauter mean diameter means physically. I thought this was a diameter based on number concentration of particles but in Fig. 2 it shows the Sauter mean diameter as being larger than the mass-based mean, which one would not expect from a numbers-based diameter.
A. The term on the right-hand side changing for different moments is the exponent i in D i . The text has been changed clarifying the assumpion of spherical particles for the expressions of the moments. The Sauter diameter is not based on number concentration but it is defined as the ratio between total volume and total surface area. To clarify we added this physical interpretation also in the text.
R2. Page 3750: line 14: add such before as settling velocity. line 15: change function to functions.
A. We changed the text as suggested by the referee. : what is the physical meaning of the different moments of density? Line 8: What are D* and rho*? Lines 11-12: you define volumetric averaged density as the mass of particles per unit volume. Does per unit volume mean per unit volume in the jet or plume? Per unit volume of each particle? If per particle (as implied later), perhaps say "average mass per unit volume of particles" A. Eq. 3 has been removed to avoid confusion and the text has been changed in the following way: "density of pumices ρ s,pum (D) with diameter < 2 mm is assumed to decrease and to reach the lithic density value when the fragment diameter decreases below 8 µm." The physical meaning of the different moments of density is explained in the text: "Otherwise, there is no reason, e.g., for ρ (1) s,j and ρ (3) s,j to be the same, as they represent length and volume weighted density averages, respectively." This means that the moment of order 1 is an average with a weight given by the particle length, the third moment is an average with the weight given by the volumes. In the same way, for example, the moment of order zero will just integrate the densities and divide by the total number of particles. D * and rho * are the constant values assumed for a monodisperse distribution. This has been now clarified in the text.
The referee is right and the definition has been corrected writing "average mass per unit volume of particles".
R2. Page 3752: line 10: Put Pfeiffer et al. (2005); Textor et al. (2006a, b) in parentheses. Equation 8: what is the physical meaning of the moments of settling velocity? Maybe this could be addressed by just moving the statement (line 3 of the next page) that they represent surface and volume-weighted averages to the line immediately above the equation.
A. The parentheses have been added. The physical meaning of the dif-ferent moments is simply a different way to evaluate the weighted average, i.e. weighting the settling velocities of particles of different sizes with a quantity proportional do the particle surface (D 2 ) or to the particle volume (D 3 ). R2. Page 3753: Line 5: change particles to particles, to make it possessive.
A. Done! R2. Page 3756: Line 13-14: bulk density means mass of particles per unit volume of par-ticles? Mass of particles per unit volume of plume mixture? If the latter, maybe say mass of particles per unit volume of plume mixture, or something similar.
A. The definition has been added to the text: "bulk density of the particles of the j-th family (i.e. the mass of particles of the j-th family per unit volume of gas-particles mixture)".
R2. Page 3757: Equations 22, 23: If this is a 3-D coordinate system, shouldnt there be three momentum conservation equations? Also, in the equation for horizontal momentum (22), Ive generally thought of the the change in momentum of the gaseous phase (first term on the right-hand side) as being equal to the horizontal momentum contained in the entrained air (2 * r * rho a tm * U e * w). Your formulation is a little different. Perhaps you could add a sentence explaining your formulation.
A. As stated in the introduction, in this work we present an extension of the Eulerian steady-state volcanic plume model presented in Barsotti et al. (2008) (derived from Bursik, 2001). Eqs. 22 and 23 use the same formulation of Eqs. 4 and 5 in Barsotti et al. (2008), with the only difference that the loss of particles is computed in terms of the moments. The reference to the original work has now been added to the text.
R2. Page 3758: Equation 24 (energy equation): Do you mention anywhere that you are not considering phase changes of water? It looks from this equation like you are not considering it, but I dont see that you mention this point in the anyplace in the text. Line 18: What do the superscript Bs represent? Also, defining the rho terms on this line as bulk densities seems misleading (at least to me). In order for the denominators on the top and bottom of eq. 26 to be consistent, the rhos must be the mass (or air, water vapor etc.) per unit volume plume mixture, not per unit volume of air or water vapor. Perhaps refer to them as the mass of each component per unit volume of plume fluid. (maybe I missed it).
A. A paragraph has been added at the beginning of the section where the main assumptions, as the fact that we are not considering phase change of water, are listed. The superscript B is used for the bulk properties and now it is defined in the text when the bulk density is introduced. As stated before, the bulk density refers to the mass of a phase per unit volume of the mixture, and thus the definition given by Eq.26 is consistent with the interpretation given by the referee.
R2. Page 3760: Line 15: change particles number to particle numbers A. The text has been changed. R2. Page 3762: Equation 40: If these are ODEs, the LHS of eq. (40) should be dy/ds rather than partial(y)/partial(s). Equation 41: perhaps add a multiplication symbol on the RHS between ds and f().
A. The partial derivatives have been changed to dy/ds. R2. Page 3764: Line 14 and elsewhere: change particles family to particle family Line 23: change particles size distribution to particle size distribution A. The text has been changed. R2. Page: 3765: Line 9: change "then writes as" to "is then written as" Equations 44: does log2 mean the log base 10 of 2? Maybe write as log 1 0(2). Or, if you mean the natural log, then write ln(2). Pages 3765-3768: this is quite a slog through this section. Its not clear exactly where youre going. Adding a sentence at the beginning of Section 4.2.1 stating your objective in deriving these formulas might help readers keep their interest.
A. "then writes as" has been changed to "is then written as". All the occurrences of "log" have been changed to "ln". Furthermore, a sentence has been added before the subsection 4.2.1 to better introduce what is the objective in deriving the formulas.
R2. Page 3769: Line 5: the term "weak plume without wind" makes no sense to me, since a weak plume is generally defined as one that is bent over by wind [Sparks et al., 1997, p. 279]. "Low-flux plume without wind" may be more accurate. Line 16: you use a normal distribution? Not lognormal? Line 13: can you provide a reference for the standard atmosphere cited here? Line 17: change "expresses" to "expressed". Line 18: Are you describing three different model runs, or a single model run with the output portrayed in three different ways?
A. The text has been changed using the term suggested by the referee "Low-flux plume without wind". The normal distribution is referred to the phi scale, it is a lognormal when the diameter is expressed in meters. A reference for the standard atmosphere has been added: "Champion, K.; Cole, A. & Kantor, A. Standard and reference atmospheres Handbook of Geophysics and the Space Environment, Air Force Geophysics Laboratory, 1985, 14" In the figure we present three different model runs with the same initial condition, but with different modeling approaches used to describe the GSD: discretization in bins, moments of the number of particles as a function of diameter expressed in meters, moments of the mass fraction of particles expressed as a function of the diamter expressed in the phi scale.
R2. Page 3772: a nicely written and illuminating summary of the Latin Hypercube and gPCE alternatives to MCMC modeling.
A. Thank you! R2. Page 3773: Still illuminating. Im not a specialist in the mathematics of these tech-niques and cant critique them. Its a little unclear to me what the form of equation 55 would be if fully expanded. For example, if zeta were a vector of 3 variables, would P1 represent three families of polynomial coefficients; one for each variable?
A. Thank you again for the positive comment. Some references have been added where the technique is explained in more details. If zeta were a vector of 3 variables, P1 would be a polynomial of the 3 variables. The polynomials P1,P2,...Pn have to be orthogonal, i.e. the integral of the product of two polynomials with index i and j (with i = j) has to be equal to zero.
R2. Page 3774: the values contoured in the lower panels of Fig. 7 were not initially very clear to me. You call them response functions in the caption, and on page 3773 (line 17). Are these the values of gamma(zeta) in eq. 55? Perhaps referring to gamma(zeta) as a response function would clarify. Also, it would help to mention that the values contoured in the lower panels are the same as those plotted in the upper panels (e.g. top mean phi for panel 1). Lines 22-30: I would say that the points you make in these few lines are the most significant of the paper, for readers interested in geologically relevant findings.
A. As suggested by the referee, we refer now to eq.55 when presenting the response functions. We also added that the values contoured in the lower panels are the same as those plotted in the upper panels: "The variables contoured in the lower panels are the same as those on the horizontal axes in the upper panels".
R2. Page 3775: Line 4: change reduce to reduces. Line 5: change relvance to rele-vance. Line 6: change entrained aid to entrained air. lines 13-14: The mean and the SD of the TGDS at the top of the eruptive column clearly reflects the corresponding values at the bottom, with a small effect on the mean size at the top of larger values of the bottom SD. What does this mean?
A. All the typos have been corrected. The sentence was really confusing and it has been rephrased in the following way: " a small effect of the bottom standard deviation on the mean size at the top, resulting in an increase in the average grain size with increasing values of the initial standard deviation".
R2. Page 3776: Line 10: Im a little confused about which of the sensitivity indices (Si or Ti) is displayed in Fig. 9.
A. The indices plotted are the main indices S i . This is now written in the caption of figure 9.
R2. Figure 2 caption, line 4: change forth to fourth. Figure 5: the light gray curves on this figure are hard to see on my computer screen. Darkening them should make them more visible. Figure 6 caption: change Two parameters Latin Hypercube to Two-parameter Latin Hypercube A. The caption of figure 2 has been corrected. In Figure 5 now the lines are more visible. The caption of figure 6 has been corrected.

Abstract.
In this paper a new integral mathematical model for volcanic plumes, named PlumeMoM, is presented. The model describes the steady-state dynamics of the plume in a 3D coordinate system, accounting for continuous variability in particle size distribution of the pyroclastic mixture ejected at the vent. Volcanic plumes are composed of pyroclastic particles of many different sizes ranging 5 from a few microns up to several centimeters and more. Proper description of such a multiparticle nature is crucial when quantifying changes in grain-size distribution along the plume and, therefore, for better characterization of source conditions of ash dispersal models. The new model is based on the method of moments, which allows description of the pyroclastic mixture dynamics not only in the spatial domain but also in the space of parameters of the continuous size-distribution of the 10 particles. This is achieved by formulation of fundamental transport equations for the multiparticle mixture with respect to the different moments of the grain-size distribution. Different formulations, in terms of the distribution of the particle number, as well as of the mass distribution expressed in terms of the Krumbein log scale, are also derived. Comparison between the new moments-based formulation and the classical approach, based on the discretization of the mixture in N discrete phases, 15 shows that the new model allows the same results to be obtained with a significantly lower computational cost (particularly when a large number of discrete phases is adopted). Application of the new model, coupled with uncertainty quantification and global sensitivity analyses, enables investigation of the response of four key output variables (mean and standard deviation of the grain-size distribution at the top of the plume, plume height and amount of mass lost by the plume during the ascent) to 20 changes in the main input parameters (mean and standard deviation) characterizing the pyroclastic mixture at the base of the plume. Results show that, for the range of parameters investigated and without considering interparticle processes such as aggregation or comminution, the, the grain-size distribution at the top of the plume is remarkably similar to that at the base and that the plume height is only weakly affected by the parameters of the grain distribution. The adopted approach can 25 be potentially extended to the consideration of key particle-particle effects occurring in the plume including particle aggregation and fragmentation.

Introduction
In the past decades, numerical simulation of volcanic eruptions has greatly advanced and models are now often able to deal with the multiphase nature of volcanic flows. This is the case, for ex-30 ample, of models describing the dynamics of pyroclastic particles in a volcanic plume, or that of bubbles and crystals dispersed in the magma rising in a volcanic conduit. Despite this, in numerical models, the polydispersity associated with the multiphase nature of volcanic flows is often ignored or largely simplified (Valentine and Wohletz, 1989;Neri at al., 2003;Dartevelle, 2003;Dufek and Bergantz, 2007;Esposti Ongaro et al., 2007;de' Michieli Vitturi et al., 2010). For instance, in 35 most of the existing conduit models, crystals and bubbles are treated as simple flow components and described by volume fractions only, while in plume dynamics and ash dispersal models the grain size distribution of pyroclasts is discretized in a finite number of classes (i.e. phases). Both approaches make proper treatment of the continuous variability of the dimension of pyroclastic particles and gas bubbles difficult. Literature results (Llewellin et al., 2002;Pal, 2003;Costa et al., 2010) clearly 40 show that this variability can largely affect relevant physical/chemical processes that occur during the transport of the dispersed phase such as, for example, the nucleation and growth of bubbles and the coalescence/breakage of bubbles and crystals in the conduit or the aggregation of pyroclastic particles in a volcanic plume.
A theoretical framework and the corresponding computational models, namely the method of 45 moments for disperse multiphase flows, have been developed in the past decades, mostly in the chemical engineering community (Hulburt and Katz, 1964;Marchisio et al., 2003), to track the evolution of these systems not only in the physical space, but also in the space of properties of the dispersed phase (called internal coordinates). According to this method, a population balance equation is formulated as a continuity statement written in terms of a density function. From the 50 density function some integral quantities of interest (namely the moments, i.e. specific quantitative measures of the shape of the density function) are then derived and their transport equations are formulated.
In this work we present an extension of the Eulerian steady-state volcanic plume model presented in Barsotti et al. (2008) (derived from Bursik (2001)) obtained by adopting the method of moments.

55
In contrast to the original works where pyroclastic particles are partitioned into a finite number of classes with different size and properties, the new model is able to consider a continuous size distribution function of pyroclasts, f (D), representing the number or the mass fraction of particles (per unit volume) with diameter between D and D + dD. Accordingly, conservation equations of the plume are expressed in terms of the transport equations for the moments of the ash particles size 60 distribution. In particular in the following we present the new multiphase model formulation based on the implementation of the quadrature method of moments McGraw (2006) and we investigate the sensitivity of the model to uncertain or variable input parameters such as those describing the grainsize distribution of the mixture. To quantify and incorporate this epistemic uncertainty affecting the input parameters (characterizing lack-of-knowledge) into our application of the model we tested two 65 different approaches, a modification of the Monte Carlo method based on Latin hypercube sampling and a stochastic approach, namely the generalized Polynomial Chaos Expansion method. This paper is organized as follows: in Section 2 we present the method of moments applied to two different descriptions of particles distribution. In Section 3 the equations of the model for the two formulations are described. Section 4 is devoted to the numerical discretization of the model In contrast to previous works, where the solid particles are partitioned in a finite number of classes with different size (Barsotti et al., 2008), we introduce here a continuous size distribution function representing the number (or mass) concentration of particles (per unit volume) as a function of the particles diameter. In general, this particle size distribution (PSD) is a function of time t, of the 80 spatial coordinate and of the diameter of the particles.
First, we present the method of moments for a particle size distribution f (D), representing the number concentration of particles (particles per unit volume) with diameter between D and D +dD, where D is expressed in meters. When more than one family of particles are present, for example lithics and pumice, we will use the subscript j to distinguish among them. Consequently, the function 85 f j (D) will denote the number concentration of particles of the j-th family.
Given a particle size distribution f j (D), we observe that its "shape" can be quantified through the moments M (i) j (Hazewinkel, 2001), defined by (1) The particular definition of f j (D) we adopt, expressing the number concentration of particles of 90 size D, allows the following physical interpretation of the first four moments: is total number of particles of the j-th family(per unit volume); -M (1) j is sum of the particles diameter of the j-th family(per unit volume); -M (2) j is total surface area of particles of the j-th family(per unit volume); j : total volume of particles of the j-th family (per unit volume) or the local volume fraction of the j-th dispersed phase, also denoted with α s,j . The multiplying factor π 6 is obtained assuming spherical particles. For particles with different shape, if volume scales with the third power of length, we can still relate the particle volume V with the particle length D through a volumetric shape factor k v such as V = k v L 3 .
We also note that the central moments (i.e., those taken about the mean) can be expressed as func-100 tion of the raw moments (i.e., those taken about zero as in Eq. (1)), and in this way it is possible to relate the moments of the distribution with the mean, variance, skewness, kurtosis. Furthermore, a mean particle size can be defined as the ratio of the moments M j . Similarly, it is possible to define 105 the mean particle length averaged with respect to particle number density L j,10 = M the sum of the lengths of particles (per unit volume) divided by the number of particles (per unit volume), and the mean particle length averaged with respect to particle volume-fraction L j,43 = The motivation for the introduction of the moments is to minimize computational costs by avoid-110 ing the discretization of the size distribution in several classes, and nevertheless to capture the polydispersity of the flow through the correct description of the evolution of the moments (Carneiro, 2011). The moments approach also allows to treat interparticles processess such as particle aggregation and fragmentation that strongly depend on and affect the GSD of the mxture (Marchisio et al., 2003). The moments and the corresponding transport velocities appear naturally in the mathematical 115 formulation as a direct consequence of the integration of the Eulerian particle equations over the diameter spectrum, as will be shown in the next section.

Moments of other quantities
In the plume model, several quantities characteristic of the particles, such as settling velocity, density and specific heat capacity, are also defined as functions of the particle diameter, and thus we can 120 define their moments as done for the particle size distribution f j (D). In general, for a quantity ψ j that is a function of the diameter D, we define its moments as As a first example, we consider here the moments of particles density ρ s . In particular, following Bonadonna and Phillips (2003), density of lithics is assumed to be constant, whereas density 125 of pumices ρ s,pum (D) with diameter D < D 2 (here equal to 2 mm) is assumed to decrease and to reach the lithic density value when the fragment diameter decreases below D 1 (here equal to 8 µm).
Substituting the expression for the particles density of the j-th particle family in Eq.
(2), we obtain the moments of the density as: (3) 130 We remark that moments of different order are generally different, they will only be equal (ρ s,j , l = m) in two limiting cases: for a monodisperse distribution with diameter D * and density ρ * s,j to be the same, as they represent length and volume weighted den-135 sity averages, respectively. For our application, we are interested mostly in the volumetric averaged density ρ (3) s,j , i.e. the average mass per unit volume of particles from now on denoted withρ s,j . The moments defined by Eq. (3) can also be used to define other properties of the gas-particles mixture. For example, it follows from the definition of the moments that if we have a mixture of a gas with density ρ g and a family of polydisperse distributions of particles with density ρ s,j = ρ s,j (D), 140 the mixture density is given by: and consequently the mass fraction of the j-th solid phase with respect to the gas-particles mixture is given by: 145 We also remark that here the gas phase is a mixture of atmospheric air, entrained in the plume during the rise in the atmosphere, and a volcanic gas component, generally water vapour. In the following, we will use the subscript atm to denote the atmospheric air and wv for the volcanic water vapour.
Differently from the approach used in Barsotti et al. (2008), where a constant settling velocity for each class is provided by the user, here several models have implemented in the code (Pfeiffer 150 et al., 2005;Textor et al., 2006a, b). For the application presented in this work, the settling velocity is defined as a function of the particle diameter and density as in Textor et al. (2006a): where k 1 = 1.19 × 10 5 m 2 kg −1 s −1 , k 2 = 8 m 3 kg −1 s −1 and k 3 = 4.833 m 2 kg −1/2 s −1 . The drag coefficient C D is a parameter accounting for the particles surface roughness, and for this work we 155 used a value of 0.75 as in Carey and Sparks (1986).
Proceeding as done for the particle density, it is possible to evaluate the moments w (i) s,j of the settling velocity w s,j (D), defined as and representing weighted integrals of the settling velocity over the size spectrum. Again, moments 160 of different order are generally different. There is no reason, e.g., for w (2) s,j and w s,j to be the same, as they represent surface and volume weighted averages, respectively.
Finally, it is possible to define the moments C (i) s,j of the particles' specific heat capacity C s,j : We observe that for the specific heat capacity generally we are not interested in a volumetric average 165 but in the mass average, denoted here with the notationC s,j and given by the following expression:

Mass fraction distribution
While in chemical engineering, where the method of moments is commonly used, the particle number distribution f j (D) is generally preferred to describe the polidispersity of the particles, in volcanology it is more common to use a mass fraction distribution γ j (φ), defined as a function of the Krumbein phi (φ) scale: where D is the diameter of the particle expressed in meters, and D 0 is a reference diameter, equal to 1 mm (to make the equation dimensionally consistent).

170
In this case, the distribution γ j (φ) represents the mass fraction of particles (mass per unit mass of the gas-particles mixture) of the j-th family with diameter between φ and φ + dφ. Again, the shape of the distribution γ j (φ) can be characterized by its moments Π i j , defined by Also in this case the particular definition of γ j (φ) allows a physical interpretation of the moments: 175 for example, the moment Π (0) j is the mass fraction of the j−th solid phase x s,j with respect to the gas-particles mixture. As done with the particle number distribution, it is possible to define a mean particle size in terms of the moments of the mass fraction distribution as Π j ; this ratio, for i = 0, gives the mass averaged diameter, corresponding to the volume averaged diameter Again, it is possible to define the moments of other quantities ψ j (φ) in terms of the continuous distribution of mass fraction γ j (φ) as For example, when the mass fraction distribution γ j (φ) is used, the mass averaged heat capacitȳ C s,j is given by the following expression: and the volumetric averaged density, i.e. the mass of particles per unit volume, can be evaluated from 3 Plume Model

190
In this section we describe the assumption and the equations of the model. As in Bursik (2001), the model assumes an homogeneous mixture of particles and gases with thermal and mechanical equilibrium between all phases. Aggregation and breakage effects are not considered and consequently density does not change with time. Finally, the model does not consider effects of humidity and water phase changes.

195
The equation set for the plume rise model is solved in a 3-D coordinate system (s, η, θ) by considering the bulk properties of the eruptive mixture (see Figure 1). The plume is assumed with a circular section in the plane normal to the centerline trajectory with curvilinear coordinate s, a top-hat profile of the velocity along the certerline, an inclination on the ground defined by an angle η between the axial direction and the horizon, and an angle θ in the horizontal plane (x, y) with respect to the 200 x−axis. These angles are needed to describe the evolution of weak explosive eruptions which are strongly affected by atmospheric conditions.
Following Bursik et al. (1992) and Ernst et al. (1996), the conservation of flux of particles with size D of the j−th family is given by: where r is characteristic plume radius, U sc represents the velocity of the plume cross section along its centerline (a top-hat profile is assumed) and p is the probability that an individual particle will fall out of the plume, defined as a function of an entrainment coefficient α as Equation (14) states that the number of particles of the j−th family with size D lost from the plume 210 is proportional to the number of particles at the plume margin, given by f j (D) · 2πr, to the settling velocity w s,j (D) and to the probability factor p.
Now, multiplying both the sides of equation (14) for D i and then integrating over the size spectrum [0, +∞], we obtain the following conservation equations for the moments M

215
If we compare our formulation with that presented in Barsotti et al. (2008), where the effects of a polydisperse solid phase are taken into account partitioning the size spectrum in a finite number N of solid classes, the set of equations (16) replaces the N mass conservation equations for the N particulate classes.
From equation (14), if we multiply both the terms by the mass of the particles of size D, given by 220 π 6 D 3 ρ s,j (D), we obtain the additional equation: and, integrating over the size spectrum: where on the left hand-side the term π 6 M (3) s,j represents the volume average bulk density of the 225 particles of the j-th family (i.e. the mass of particles of the j-th family per unit volume of gas-particles mixture, denoted with the superscript B, ρ B s,j ), while on the right-hand side the term [w s,j ρ s,j ] (3) represents the mass averaged settling velocity of the particles of the j-th family multiplied by the volume averaged particles density. Equation (18)

235
First of all, we derive the conservation equation for the mixture mass. As in the plume theory, we assume that the entrainment, due to both turbulence in the rising buoyant jet and to the crosswind field, is parameterized through the use of two entrainment coefficients, α and γ . The theory assumes that the efficiency of mixing with ambient air is proportional to the product of a reference velocity (the vertical plume velocity in one case and the wind field component along the plume cen-240 terline in the other), by α and γ (Morton, 1959;Briggs, 1975;Wright, 1984;Weil, 1988). Thus, following Hewett et al. (1971) and Bursik (2001), we define the entrainment velocity U as a function of windspeed, U atm , as well as axial plume speed, U sc : where α |U sc − U atm cos φ| is entrainment by radial inflow minus the amount swept tangentially 245 along the plume margin by the wind, and γ |U atm sin φ| is entrainment from wind. With this notation, the total mass conservation equation solved by the model becomes stating that the variation of mass flux (left-hand side term) is due to air entrainment (first right-hand side term) and loss of solid particles (second right-hand side term), as obtained from Eq. (18).

250
From Newton's second law and the variation of mass flux, we can derive also the horizontal and vertical components of the momentum balance solved by the model as: and

255
where the two components of plume velocity along the horizontal and vertical axes are u and w, respectively, and they are linked by the relation U sc = √ u 2 + w 2 . In the right-hand side of equation (21) the terms related to the exchange of momentum due to the wind (Barsotti et al., 2008) and to momentum loss from the fall of solid particles appear. Similar contributions are evident in the right-hand side term of equation (22) where the vertical momentum is changed by the gravitational 260 acceleration term and the fall-out of particles. Now, following the notation adopted above and denoting with T the mixture temperature, the equation for conservation of thermal energy solved by the model writes as The first term on right-hand side describes the cooling of the plume due to ambient air entrainment, 265 the second one takes into account atmospheric thermal stratification, and the third term allows for heat loss due to loss of solid particles. Again, this last term is obtained writing the heat loss for the particles of size D, and then integrating over the size spectrum. A thermal equilibrium between solid and gaseous phases is assumed. In Eq. (23) C atm and C mix are the heat capacity of the entrained atmospheric air and of the mixture, respectively, the latter being defined as: or, in terms of the bulk densities ρ B atm = x atm ρ mix , ρ B wv = x wv ρ mix and ρ B s,j = π 6 M (3) jρ s,j , as From this expression, if we multiply all the terms at the numerator and the denominator of the righthand side by U sc r 2 and we differentiate with respect to s, we obtain after some cancellation and 275 algebra manipulations the following equation for the variation of the mixture specific heat with s: Now, substituting the expressions for the derivatives appearing in each term on the right-hand side, we obtain the equation for the variation rate of mixture specific heat in terms of the moments: Similarly, a gas constant R g can be defined as a weighted average of the gas constant for the entrained atmospheric air R atm and the gas constant of the volcanic water vapour R wv and a conservation equation can be derived, knowing that the variation of gaseous mass fraction with height is solely due to entrained air: This formulation reduces, for particular cases, to the expressions of Woods (1988) and Glaze and Baloga (1996). Equations (27) and (29) are needed in order to close the system of equations and recover the new values of the temperature and the mixture density once the system of ordinary differential equations is integrated. Otherwise, without the solutions of Equations (27) and (29), we 290 should use the old values of ρ mix and C mix at s to obtain the values of the temperature at s + ds from the lumped term (ρ mix U sc r 2 C mix T ) obtained integrating Eq. (23).

Mass fraction distribution
Similarly as done for the distribution of particle number f j (D) and the moments M In this case, the conservation of mass flux of particles with size φ of the j−th family write as: Multiplying both sides of the equation by φ i and integrating over the size spectrum [−∞, +∞], we obtain the following conservation equations for the moments of the continuous distributions γ j (φ): For i = 0, the equations of conservation of the moments give: expressing the loss of mass flux of the particles of the j−th family and thus we can write the total mass conservation equation as From the variation of mass flux, as done for the distribution of particle number f j (D) and the 310 moments M (i) j , we derive the horizontal and vertical components of the momentum balance: The equation for conservation of thermal energy is and the equation for the variation rate of mixture specific heat in terms of the moments of the mass fraction distribution write as: The formulation of the equations for the gas constant R g and the coordinates of the (x, y, z) remain 320 unchanged.

Numerical scheme
The plume rise equations are solved with a predictor-corrector Heun's scheme (Petzold and Ascher, 1998) that guarantees a second-order accuracy, keeping the execution time in the order of seconds.
If we rewrite the system of ordinary differential equations with the following compact notation: where y is the vector of the quantities in the left-hand sides of the conservation equations presented in the previous section, then the procedure for calculating the numerical solution by way of Heun's method (Süli and Mayers, 2003) is to first calculate the intermediate valuesỹ i+1 and then the solution y i+1 at the next integration point 330ỹ i+1 = y i + dsf (s i , y i ), predictor step

Quadrature method of moments
We observe that to calculate the right-hand side for both the predictor and corrector step we need not only the moments M (i) , but also the additional moments [w s ] i , [w s ρ s ] (i) and [w s ρ s C s ] (i) . As in Marchisio and Fox (2013), the integral in the definition of these moments is replaced by a quadrature 335 formula and the moments, for a generic variable ψ = ψ(D), are approximated as: Here ω l and D l are known as "weights" and "nodes" (or "abscissae") of the quadrature, respectively, and the accuracy of a quadrature formula is quantified by its degree. The degree of accuracy is equal to d if the interpolation formula is exact when the integrand is a polynomial of order less 340 than or equal to d and there exists at least one polynomial of order d + 1 that makes the interpolation formula inexact. In particular, an N −point Gaussian quadrature rule, is a quadrature rule constructed to yield an exact result for polynomials of degree 2N − 1 or less by a suitable choice of the nodes D l and weights ω l for l = 1, . . . , N (Golub and Welsch, 1969).
The Wheeler algorithm, as presented in Marchisio and Fox (2013), provides an efficient O(N 2 ) 345 algorithm for finding a full set of weights and abscissas for a realizable moment set. The resulting nodes D l are always within the support (and therefore represent realizable values of the particle size), and the weights ω l are always positive, ensuring that, when the quadrature is used, accurate results are obtained Marchisio and Fox (2013). Nevertheless, these properties are respected only if the moment set is realizable, meaning that there exists a particle size distribution resulting in that 350 specific set of moments.
A strategy that might overcome the problem of moment corruption (i.e. the transformation during the integration of the moment-transport equations of a realizable set of moments into an unrealizable one) is replacing unrealizable moment sets as soon as they appear. An algorithm of this kind was developed by McGraw (McGraw, 2006). The algorithm first verifies whether the moment set is 355 realizable (by looking at the second-order difference vector or by looking at the Hankel-Hadamard determinants (Gautschi, 2004)). If the moment set is unrealizable it proceeds with the correction.
In the numerical model here presented, the implementation of the correction algorithm of Wright is derived from the version presented in Marchisio and Fox (2013).
Thus, in both the predictor and corrector step, the following algorithm is used:

360
the nodes D j,l and weights ω j,l are calculated with the Wheeler algorithm for l = 1, . . . , N ; the quadrature formula (41)  the right-hand side of the ODE's system (39) is evaluated explicitly; the solution is advanced with the predictor (or the corrector) step of the Heun's scheme; for each particle family j, the moments M We observe that if the j−th family of particles is monodisperse with diameterd j , the Wheeler algorithm fed with the first two moments only gives as result a single quadrature node D j,1 =d j with weight ω j,1 = 1. This allows us also to use the model for the simplified case where the solid particle distribution is partitioned in a finite number of classes with constant size, assigning to each 370 class a monodisperse distribution.

Initial condition
Initial conditions at the vent include the initial plume radius (r 0 ), mixture velocity (U sc,0 ) and temperature (T 0 ), gas mass fraction (n 0 ) and the particles size distribution through the initial moments M (i) 0 . In the next section we derive analytically the moments of a specific initial distribution (a nor-375 mal distribution in the Krumbein scale) for both the formulations based on the number of particles as a function of the particles diameter expressed in meters and the formulation based on the mass concentration expressed as a function of the phi scale.

Lognormal distribution
For the application presented in this work, the initial distribution f (D) at the base of the plume is 380 defined as a function of the particles diameter expressed in meters (m), in order to give a corresponding normal distribution with parameters µ and σ for the mass concentration expressed as a function of the Krumbein phi (φ) scale (when all the particles have the same density): where K 0 is a parameter that has to be chosen in order to satisfy the initial condition on the solid 385 mass fraction.
Given the parameters µ and σ, the initial distribution f (D) is then written as: where C 0 , analogously to K 0 , is a parameter that has to be fixed in order to satisfy the initial condition prescribed for the mass (or volume) fraction of particles. 390 We observe that if we introduce the following re-scaled variables for the diameter, the mean and the variance: then it is possible to rewrite the particle distribution f (D) in terms of a lognormal distribution in the variableD with parametersμ andσ: Consequently, we can evaluate the moments M (i) of f (D) analytically from the moments of the lognormal distribution as: and we obtain, for the third moment: where α 0 s is the initial volume fraction of the particles in the solid-gas mixture.
From the expressions of the moments it follows also that, if the mass concentration expressed as a function of the Krumbein scale has a normal distribution, the Sauter mean diameter D A expressed in meters can be evaluated as or, if expressed in φ, as Processes involving the mutual interaction between particles and the interaction between the particles and the carrier fluid (friction and cohesion between the particles; viscous drag; chemical reac-410 tions between fluid and solid components) operate at the surface of the particles. For this reason the Sauter mean diameter, based on the specific area of the particles, is a convenient descriptor and it is important to remark that it differs from the mean µ of the lognormal distribution by a factor proportional to the variance σ 2 . For numerical models describing the multiphase (particulate) nature of the matter and which approximate the particle size distribution with an average size, it is hence more ap-

415
propriate to use, as particle size representative of a lognormal distribution, the Sauter mean diameter than the mean diameter µ. The difference between the two approximations is smaller the narrower the particle size distribution. We must also remark that, while for particles in the inertial-dominated regime (e.g. Re p > 2000) Loth et al. (2004) showed that the Sauter mean diameter is the effective diameter, regardless of particle shape, particle size distribution, particle density distribution or net 420 volume fraction, for particles in the creeping flow regime (Re p << 1) the effective mean diameter is the volume-width diameter.
When the Sauter mean diameter is used, also the variance and the standard deviation SD should be based on the specific surface area (Rietema, 1991). Hence: 425 or, expressed as a function of the moments: Finally, we note that if the particle density is constant and the mass concentration expressed as a function of the Krumbein scale has a lognormal distribution and both the Sauter mean diameter L 32 = M (3) /M (2) and the mean particle length averaged with respect to particle number density 430 L 10 = M (1) /M (0) (or if the first 4 moments) are known then we can solve for the re-scaled mean and varianceμ andσ the following system: Once the re-scaled mean and variance are known, we can obtain µ and σ in the Krumbein φ scale.
When the initial distribution is expressed for the mass fractions instead of the particle number, 435 and the mass fraction written as a function of the Krumbein scale has a normal distribution with mean µ and variance σ 2 , then the continuous distribution is given by Eq. (42). We observe that this expression of the distribution is not based on the assumption of constant density for the particles of different size.
In this case, the moments Π (i) are given by the following expression 440 where the symbols and !! denote the integer part and the double factorial (n!! = m k=0 (n − 2k), where m = n/2 − 1), respectively. Now, as the 0−th moment is equal to the mass fraction of particles, we obtain K 0 = x s . Furthermore, we observe that the mass fraction averaged diameter in the φ scale is given by the ratio 445 Π (1) /Π (0) , while the variance of the mass fraction distribution can be evaluated as Π (2) Π (0) − (Π (1) ) 2 /(Π (0) ) 2 .
These two quantities correspond to the parameters (µ, σ 2 ) generally used to describe the mass fraction when a normal distribution in the φ scale is assumed. For this reason, when we want to track the changes of the mass fraction averaged diameter and its standard deviation (or variance) in the φ scale during the plume rise, it is preferred to use a formulation based on the moments Π (i) than the 450 moments M (i) .

Simulation inputs
We applied the model to three different test cases with different vent and atmospheric conditions: -Test Case 1 -low-flux plume without wind;
The parameters used for the different test cases are listed in Table 1, while the atmospheric conditions are plotted in Fig. 3. For the low-flux plumes a mass flow rate of 1.5 × 10 6 kg/s has been fixed, while for the strong plume the value is 1.5 × 10 9 kg/s. The temperature pressure and density profiles 460 used for the test case without wind (Test Case 1) are those defined by the International Organization for Standardization for the International Standard Atmosphere (Champion et al., 1985), while the profiles for the other two test cases come from reanalysis data. For all the runs presented here, a single family of particles has been used, with a normal distribution (with parameters µ and σ) for the mass concentration as a function of the diameter expressed in 465 the φ scale and with density varying with the particle diameter.
We first present a comparison of the plume profiles obtained with the 3 different descriptions presented in the previous sections and highlighted in the three colored boxes of Fig. 2 for the Test Case 2: method of moments for the particle number being function of the size expressed in meters; method of moments for the particle mass fraction being function of the size expressed in the φ scale; 470 discretization in uniform bins in the φ scale. For this comparison, the mass flow rate at the vent is 1.5 × 10 6 kg/s and a rotating wind is present, as shown in Fig. 3, while the mean and the standard deviation of the initial total grain size distribution are respectively 2 and 1.5, expressed in the φ scale.
The results of the numerical simulations obtained with the three different formulations are presented in Fig. 4 and they perfectly match, showing that the method of moments (dotted lines), both applied 475 to the continuous distribution of the particle number (red) or to the mass distribution (green), gives .5-2.5 0.5-2.5 0.5-2.5 Table 1. Input parameters used for the numerical simulations. Vent height is the elevation of the base of the column above sea level. The values ρ1,2 and D1,2 are used to compute the density of the particles as a function of the diameter, according to the formulation of Bonadonna and Phillips (2003). The values reported for µ and σ define the range used for the uncertainty quantification and sensitivity analysis.  In blue the profiles obtained using 13 bins, in red the profiles obtained using a continuous distribution of the particle number density and in green using a continuous distribution of the mass fraction.

Simulation results
In this section we want to study the variation during the ascent of solid mass flux (due to particles settling) and of the mean and the variance of the mass distribution along the column. As shown in the In Fig. 5 we present the results relative to the Test Case 2 for an initial particle size distribution with mean diameter 2 and standard deviation 1.5, expressed in the φ scale. In the left and middle panels the mean, the variance and the skew of the mass fraction distribution are shown respectively, while in the right panel the cumulative loss of solid mass flux is plotted as a percentage of the initial 495 value. We observe a decrease in the mean size of the particles, due to the different settling velocities of particles of different sizes. A decrease in the variance of the size distribution with height is also observed from the second plot. We remark that the particles have a normal distribution only at the base of the column (resulting in a null skewness), and the negative skew at the top of the column indicates that the tail on the left side of the grain size distribution is longer than the tail on the 500 right side, i.e. the mass is more concentrated on the right of the spectrum of particle sizes (finer particles). For this reason we do not have to look at the mean and the variance plotted in Fig. 5 as the parameters of a normal (and symmetric) distribution. Nonetheless changes in the mean, the variance and the skewness are observed, we remark that these changes are quite small and for this reason the parameters of the total grain size distribution at the top of the eruptive column are a good 505 approximation of the parameters at the base of the column, and vice versa. However, this is true for the specific input condition of this test case and not in general. For this reason, it is important to quantify the uncertainty of this assumption for different initial total grain size distributions and different atmospheric conditions.

Uncertainty and sensitivity analysis 510
When dealing with volcanic processes and volcanic hazards, our understanding of the physical system is limited, and vent parameters (volatile contents, temperature, grain size distribution, etc.) are often not well-constrained or are constrained with significant uncertainty. These factors mean that it is difficult to predict the characteristic of the ash cloud released from the volcanic column with certainty. An alternative is to quantify the probability of the outcomes (for example the grain size 515 distribution at the top of the column) by coupling deterministic numerical codes with stochastic approaches. It is our goal in this work also to assess the ability to systematically quantify the uncertainty and the sensitivity of the plume model outcomes to uncertain or variable input parameters, in particular to those characterizing the grain size distribution at the base of the eruptive column.
Uncertainty quantification (UQ) or nondeterministic analysis is the process of characterizing input An alternative approach to uncertainty quantification is the so-called generalized Polynomial utilizing the notions of projection, orthogonality, and weak convergence (Ghanem and Red-Horse, 1999). PCE was developed by Norbert Wiener in 1938 and it soon become widely used because of its efficiency when compared to Monte Carlo simulations. The term "Chaos" here simply refers to the uncertainties in input, while the word "Polynomial" is used because the propagation of uncertainties 550 is described by polynomials. If ζ is the vector of uncertain input variables, the aim of the gPCE is to express the response function Y in the form of a polynomial ξ as follows: where P 1 , . . . , P m are polynomials which form an orthogonal basis. The choice of the polynomials basis depends on the probability distribution of the input variables. In particular, for a uniform 555 distribution, the basis of the expansion is given by the Lagrange polynomials. For the application presented in this work the coefficients of the expansion have been evaluated using a spectral projection where the computation of the required multi-dimensional integrals is based on the tensor product of one-dimensional Gaussian quadrature rules. In order to compute the quadrature points, the grid used in our work is the Clenshaw-Curtis grid (Fig. 6, right), representing a good solution 560 for a multi-dimensional Gaussian quadrature with a small number of variables (Eldred and Burkardt, 2009).
We present here the results of several tests performed coupling the plume model with the Dakota toolkit (Adams et al., 2013)  Monte Carlo simulations, with a number of runs required to produce the same accuracy reduced by a factor 10 (81 simulations vs 1000 simulations). If more parameters were varied, the computational cost would increase for both gPCE and LHS, although the advantage of gPCE would be reduced.
As mentioned previously, the aim of the gPCE is to express the output of the models as polyno-585 mials and these polynomials can be used to obtain response surfaces for the output parameters as functions of the unknown input parameters through the polynomials befined by Eq. (54). In the four bottom panels of Fig. 7  the plume height for this test case shows a non-linear dependency but at the same time a small sensitivity to the initial grain size distribution, with changes, for the specific conditions here considered, 595 smaller than 1% of the average height. This can be explained by the fact that a large amount of air is entrained in the column during the ascent and the contribution of the solid fraction to the overall dynamics becomes small compared to that exerted by the gas. Finally, we observe that the loss of particles is mostly controlled by mean size of the TGSD.
In Figure 8 the same contour plots are shown for the polynomial expansion computed for Test 600 Case 2 (top) and Test Case 3 (bottom) with 81 quadrature points. The results show again that the total grain size distribution at the base of the vent represents a reasonable approximation of that at the top of the column. For these test cases, both the column height and the solid mass lost appear to be mostly controlled by the mean size of the TGSD at the base of the column, with a small sensitivity of the height to the initial grain size distribution. We also observe that the maximum percentage 605 of loss in the solid mass flux is about 15% for the strong plume simulations, and it is attained for larger mean sizes and smaller variance of the initial TGSD. This value is noticeably smaller than that obtained for the weak bent test case (≈ 40%) and for the weak test case without wind (≈ 60%).
Despite the loss of particles, in both the cases the range of variation of the column height is quite small and, as previously mentioned, this is due to the large amount of air entrained in the volcanic 610 column that reduces the contribution of the solid fraction to the overall dynamics. As an example to understand the relevance of the entrained air, for a simulation performed for the low-flux plume without wind and with µ = 2 and σ = 1.5 in the φ scale, the mass flow rate at the top of the column is 1.2 × 10 8 kg/s, compared to the value at the base of 1.5 × 10 6 kg/s.

Sensitivity analysis
With the polymomial chaos expansion it is also possible to easily obtain the variance-based sensitivity indices (Saltelli et al., 2008) with no additional computational cost. In contrast with some instances, where the term sensitivity is used in a local sense to denote the computation of response    V ar xi [E(Y |x i )] against the total variance V ar(Y ). Formulas for the indices are: and where Y = f (x) and x −i = (x 1 , . . . , x i−1 , x i+1 , . . . , x m ). Similarly, it is also possible to define the sensitivity indices for higher order interactions such as the two-way interaction S i,j . The calculation of S i and T i requires the evaluation of m-dimensional integrals which are typically approximated by Monte-Carlo sampling. However, in stochastic expansion methods, it is possible to approximate the 635 sensitivity indices as analytic functions of the coefficients in the stochastic expansion.
The results of the sensitivity analysis for the four outputs and the three test cases investigated are presented in the bar plot of Figure 9. For each of the four groups (one for each of the different output functions) the three bars represent the main sensitivity indices for the three test cases (test 1 on the left, test 2 in the middle and test 3 on the right) while the different colors are for the 640 sensitivity indices with respect to different variables (blue is for the mean of the initial TGSD, green for the standard variation of the initial TGSD and brown for the 2nd order coupled interaction).
Again, the sensitivity analysis confirms that the mean and the standard deviation of the grain size distributions at the top of the eruptive column are controlled primarily by the respective parameters for the weak test case without wind the dispersion of the distribution and second-order interaction also play a major role in controlling plume height variability. However, as already observed with the uncertainty quantification analysis, we remark that the variability in the plume height, when the mean and the standard deviation of the TGSD vary in the investigated ranges, is extremely small for 650 all the test cases (less than 1% with respect to the average values) and thus the investigation of how the variability in model output can be apportioned to the variability in individual input variables is less relevant for the plume height than for the other output parameters.

Conclusions
In this work we have presented an extension, based on the method of moments, of the Eulerian 655 steady-state volcanic plume model presented in Barsotti et al. (2008) (derived from Morton (1959); Ernst et al. (1996); Bursik (2001)). Two different formulations, one based on a continuous distribution of the number of particles as a function of the size and a second based on the continuous distribution of the mass fraction, have been presented. The tracking of the moments of mass distribution, defined as a function of the Krumbein phi scale, has the advantage that with the first three 660 moments only we are able to recover the mean and the standard deviation of the total grain size distribution. The results of a comparison between the two formulations based on the method of moments and the classical formulation based on the discretization of the mass distribution in bins show that the different approaches produce the same results, with an advantage of the method of moments in terms of computational costs. Furthermore, a formulation based on continuous description 665 of particle size, is better suited to properly describe complex inter-particle processes such as particle aggregation and fragmentation that are likely to play an important role in the plume evolution. In particular, the method of moments has already been successfully applied to model aggregation and breakage precesses in particulate systems (Marchisio et al., 2003).
An uncertainty quantification analysis has also been applied to the formulation based on the mo-670 ments of the mass distribution. The results show, for the range of conditions here investigated and neglecting likely relevant interparticle processes such as particle aggregation and communition, a small change of the mean and variance of the particle mass distribution along the column, indicating that the total grain size distribution at the base of the vent represents a reasonable approximation of that at the top of the column. Furthermore, based on the plume model assumptions and outcomes, we 675 observe a small sensitivity of the plume height to the initial grain size distribution, with variations of the order of tens of meters for a plume rising to several kilometers.
For the application presented in this work, involving only two parameters, the comparison between the latin hypercube sampling technique and the generalized polynomial chaos expansion method shows that the latter only requires 81 simulations to produce the same results, in terms of cumulative