PLUME-MoM-TSM 1.0.0: A volcanic columns and umbrella cloud spreading model

In this paper a new version of the integral model for volcanic plumes PLUME-MoM is presented. The model describes the steady-state dynamics of a plume in a 3-D coordinate system, and the two-size moment (TSM) method is adopted to describe changes in grain-size distribution along the plume, associated with particle loss from plume margins and with particle aggregation. For this reason, the new version is named PLUME-MoM-TSM. For the first time in a plume model, the full 5 Smoluchowski coagulation equation is solved, allowing to quantify the formation of aggregates during the rise of the plume. In addition, PLUME-MOM-TSM allows to model the phase change of water, which can be either magmatic, added at the vent as liquid from external sources, or incorporated through ingestion of moist atmospheric air. Finally, the code includes the possibility to simulate the initial spreading of the umbrella cloud intruding from the volcanic column into the atmosphere. A transient shallow water system of equations models the intrusive gravity current, allowing to compute the upwind spreading. 10 The new model is applied first to the eruption of Calbuco volcano in southern Chile in April 2015, and then to a sensitivity analysis of the upwind spreading of the umbrella cloud to mass flow rate and meteorological conditions (wind speed and humidity). This analysis provides an analytical relationship between the upwind spreading and some characteristic of the volcanic column (height of the neutral buoyancy level and plume bending), which can be used to better link plume models and volcanic-ash transport and dispersion models. 15


Introduction
In the last years, the simulation of ash dispersal in the atmosphere has become an ordinary activity for volcanic observatories and meteorological offices. Volcanic ash represents one of the major hazard associated with explosive eruptions, affecting both the zones proximal to the volcanic vent, where it can damage vegetation, infrastructures and pose a health risk to people, and 20 the surrounding atmosphere, requiring effective and prompt actions to regulate flights near to the ash plume (Bonadonna et al., 2020). spreading can lead to an important radial expansion, including upwind transport, which is sometime omitted in numerical simulations by the simple coupling column models with VATD models. In Costa et al. (2013), an analytical model describing the spreading of the umbrella cloud as a gravity current was coupled with an advection-diffusion-sedimentation model, and 60 the conditions under which gravity-driven transport were analyzed. Baines (2013) studied the dynamics of this gravitational spreading into a density-stratified crossflow, deriving a semi-analytic model able to produce a good description of the behaviour of the ash cloud within the vicinity of intrusion from the column. A different approach to model the spreading of the umbrella cloud is presented in Johnson et al. (2015) where, following the work of Ungarish and Huppert (2002), the flow is modeled as a transient shallow-layer gravity current subject to gravitational forces and drag between the current and the surrounding not possible to retrieve the amount of mass associated to a size interval. This represents a problem when the output of a plume model is used as input of a model for transport and deposition of ash in the atmosphere.
To cope with this problem, we describe here a hybrid method combining the advantages of the sectional and moment methods. With this hybrid approach, the size spectrum is still partitioned in sections (or bins), and then the transport equations 105 for a limited number of moments for each section are solved. The number density function (NDF) is then represented by using a piecewise reconstruction, continuous on each section. Among these hybrid methods, in several enginnering context, like spray modeling, the two-size moment (TSM) method has proven to produce good results for size reduction simulations (i.e. evaporation) and for coalescence, by using in each section a second-order accurate reconstruction of the NDF from the moments (Nguyen et al., 2016), leading also to a good representation of the NDF with a small number of sections. For this 110 reasons, TSM is a good candidate for the modeling of the aggregation processes occuring in the volcanic column. In this paper, we introduce a modified version of the TSM method presented in Nguyen et al. (2016). While in the original formulation the internal variable is the volume (or the radius) of the particles and in each section the NDF is reconstructed with a linear function, here the NDF is defined as a function of particle mass. This greatly simplify the formulation when aggregation is considered (two particles of mass M 1 and M 2 respectively aggregate in a new particle of mass M = M 1 + M 2 ). In addition, 115 while in Nguyen et al. (2016) for some section the support of the linear reconstruction (the subset of the section containing the elements which are not mapped to zero) can be smaller than the section width (which means that the existing particles size spectrum do not cover entirely the respective section), here in those section we switch from a linear reconstruction to a power law, with zero value at one of the boundaries of the section. The advantages of this different approach will be explained in details in the following of the section. 120 The starting point of the TSM method is the partition of the size spectrum in n sections (or bins). Here, we use a uniform partition φ i = φ 0 + i∆φ, i = 0, . . . , n s , in the non-dimensional Krumbein φ 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). Given the particles density ρ = ρ(φ) and their sphericity, it is possible to associate to each size of 125 the original partition φ i the corresponding mass particle m i , and thus to obtain a partition of the particles mass spectrum in n sections corresponding to the n bins in the Krumbein scale. In this way, the internal variable characterizing the variability of the particle population becomes the mass, and its variability can be characterized by the ash particles number density function, here denoted by η (m, x, t). (N umber Density F unction, N DF ) (2) 130 This function represents the number of ash particles per unit volume of the gas-particle mixture in the infinitesimal mass interval [m; m + dm] at the location x in space and at the time t (according to this definition, the number density function η(m, x, t) has units L −3 M −1 ). With this choice, when we integrate the number density function in a fixed control volume V and over a mass interval [m 1 ; m 2 ], we obtain the number of particles in V at time t with mass between m 1 and m 2 .
One important property of the NDF adopted is that the moments of the distribution correspond to quantities that have meaningful physical interpretations and are therefore, at least in principle, directly measurable in the field, as well as in experiments.
For example, we observe that the moment of order zero, m (0) , represents the total number of ash particles per unit volume of the mixture, while the moment of order one, m (1) , represents the total mass of ash particles per unit volume of the gas-particle 140 mixture, i.e. the bulk density of the ash particles, here denoted with ρ B p . From a physical point of view, these moments are commonly the main variables of interest, since when there is no aggregation the global moment of order 0 is conserved, while with aggregation the global moment of order 1 is conserved. In addition, dividing the first moment by the average ash particle density, we obtain the volume fraction of particles with respect to the volume of the gas-particle mixture.
As previously stated, the mass spectrum is partitioned in n sections, and we can define the local moments of order zero and 145 one on the k−th section defined by the mass interval [m k−1 , m k ]: Also in this case, the local moments have an important physical interpretation. In fact, the local moments of order one M k , k = 1, . . . , n, represent the mass within each mass section per unit volume of the mixture and thus, when normalized, they give the total grain size distribution. Solving for appropriate transport equation of the moments M k allows then to describe the 150 evolution of the total grain size distribution within the plume. We observe that, differently from the global moment of order one, the local moments of order one are not conserved by aggregation, because mass moves from sections corresponding to finer particles to sections corresponding to coarser ones.
The local moments of other quantities ψ(m) which depend on the mass of the particles, as for example the density, or the settling velocity, can be defined in terms of the continuous distribution of mass fraction η(m) as This term has the same units of ψ(m) and represents its integral average, weighted differently according to the order of the moment. For example ψ (0) k represents the average of ψ in the mass interval [m k−1 , m k ] weighted with the number of ash particles of different masses, while ψ (1) k represents a mass-weighted average in the same interval. When the moments of the mass distribution η are known for each section, but not the actual expression of η, the integral in 160 Eq. (5) must be approximated with a quadrature formula, and to do this we need the "weights" and the "nodes" (or "abscissae") of the quadrature. Here, we use an approach based on the Gauss-Legendre quadrature, which can be written as 6 https://doi.org/10.5194/gmd-2020-227 Preprint. Discussion started: 2 October 2020 c Author(s) 2020. CC BY 4.0 License.
where the i−th node x i is the i−th root of the Legendre polynomial P n (x) and the weights ω i are given by the Abramowitz-Stegun 1972 formula (Abramowitz and Stegun, 1965). This n−points quadrature formula is exact for polynomials of degree 165 2n − 1 or less. From these quadrature nodes and weights defined for the standard interval [−1; 1], it is possible to compute the nodes and weights for the integral over the interval [m k−1 , m k ]: wherẽ With this approach, when the integrand of Eq. (5) changes (for example because the mass distribution changes in space or time), there is no need to reevaluate the quadrature nodes and weights, which depends only on the interval and can be computed in advance. The only thing to update is the value of the integrand at the quadrature points. In the TSM approach, as presented in Nguyen et al. (2016), mass distribution η(m) is approximated in each section by a linear function. Here, non-linear approximations ≈η k (m) of the mass distribution are considered. On each section k, a continuous reconstruction is assumed: where three possible continuous reconstructions are used for eachη k (m), as illustrated in Fig. 2: The parameters α k , β k (representing the values of the reconstructed function at m k−1 and m k ) and γ k are computed in order to have In this way, the domain of each functionη k (m) is always the whole k−th section, and there is no need to update the quadrature nodes.
So far we presented the TSM approach for a single family of particles. When more than one particle family is considered, for example if we consider aggregation and we want to distinguish between aggregated and non-aggregated particles, the subscript 185 i p is used to distinguish the mass distribution functions η ip (m, x, t), i p = 1, . . . , n p , where n p is the number of particles families. Similarly, the notation m (j) ip is introduced for the j−th moment of the i p −th distribution, and the notations N ip,k and M ip,k for the number and mass of particles of the i p −th family in the k−th section, respectively.

Plume equations
In this section we present the set of conservation equations describing the steady rise (d/dt = 0) of the volcanic mixture 190 through the atmosphere. As in Bursik (2001), the model assumes a homogeneous mixture of different phases with thermal and mechanical equilibrium between all of them. In this version of the model, in addition to air and particles, we also explicitly model the presence of water, which can be in gas, liquid or solid state, and the presence of additional volcanic gases (transported as passive components). The equation set for the plume rise model is solved in a 3-D coordinate system (x,y,z) centered horizontally at the plume base location, and vertically at the sea level. The plume is assumed with a circular section with radius 195 r, parallel to the horizontal plane (x,y) and orthogonal to the vertical coordinate z (see Fig. 3).
In the (x,y,z) coordinate system, we denote with u, v and w the three components of plume velocity, and with U sc = √ u 2 + v 2 + w 2 its magnitude. Furthermore, we define φ as the angle between the axial direction and the horizontal plane, for which it holds the following relationship We also consider the presence of an horizontal atmospheric wind, with velocity components u a , v a and magnitude U atm .
Because of the wind, the location of the centers of the horizontal section, which defines the plume center-line, can be bent.
The rising multiphase mixture is formed by different families of solid particles (differing in shape, density, thermal properties), dry air, water (which can be in vapour, liquid and ice form) and other volcanic gases, such as CO 2 and SO 2 . When  more than one particle family is considered, they are denoted with the subscript i p = 1, . . . , n p ; in addition, if aggregation is 205 considered, we use the first n p − 1 indexes for the original particles and the last one for the aggregated particles. With this notation, for the particles with mass m of the family i p , we can write the following equation for the vertical mass flux through an horizontal cross-section of the plume (see Fig. 3): where both the plume radius r and the vertical component of the plume axial velocity w vary with the vertical coordinate z.

210
Equation (12) describes the variation with the vertical coordinate z of the mass flux of particles of size m, due to fallout of pyroclasts (first term on the right-hand side) and aggregation (second term on the right hand side). Here, w ip (·) is the particle settling velocity (see Appendix A2) and the coefficient p, which is a function of the radial entrainment coefficient α, expresses the probability that an individual particle will fall out of the plume: The last term on the right hand side of Eq. (12) represents the aggregation term which, following Von Smoluchowski, can be written as: where β(z, m, m ) is the aggregation (or coagulation) kernel describing the rate at which a particle of mass m aggregate with a particle of mass m . The equations state that the aggregation of particles of mass m and m , forming a particle of size m + m , is proportional to the number density of particles of such masses (η(m) and η(m )) and to the aggregation kernel. The negative term on the right hand side of the first equation represents the loss of ash particles of the original families due to aggregation.

225
The positive term on the right hand side of second equation represents the increase of particles of mass m at the vertical coordinate z due to aggregation of smaller particles. The aggregation kernel can be defined as the product of two quantities: the "collisions frequency" between particles of mass m and m and the "efficiency of aggregation" (i.e., the probability of particles of mass m coalescing with particles of mass m when colliding). In PLUME-MoM-TSM several aggregation kernel are implemented, from simple constant kernels to more complex formulations describing wet aggregation, where the collisions 230 frequency depends on particles properties (size, density) and the efficiency of aggregation depends on the amount of water in the mixture (see Appendix A3).
Let us consider now a mass interval Ω, which can be the whole ash particle mass spectrum or a smaller subset, as the k−th section for the family i p , defined by [m ip,k−1 ; m ip,k ]. Multiplying Eq. (12) for m i and integrating over Ω, we obtain the following differential equations for the moments: where the last term is defined by: − Ω m i η np (m, z) We remark that so far we have written a general form of the moment equations, which are valid independently from the particular approach used to solve them (QMOM or TSM). As previously stated, when the TSM method is adopted, the integral 250 terms on the right-hand side of the two previous equations (w (·) ip andŜ (·) ip,k ) are calculated first by reconstructing the continuous mass distributions in each section from the known moments (number and mass of ash particles in the k−th section per unit volume of the mixture), and then approximating the integrals with the Gauss-Legendre quadrature formulas.
We also observe that when the equations for the mass of particles are added together (summed over the particles families and over the sections), the aggregation terms cancel out, and the following transport equation is obtained: which gives the total mass of ash particles per unit volume of the mixture lost because of particle sedimentation from the plume.
The right-hand side of the previous expression can be used in the mass conservation equation for the volcanic mixture accounting for atmosphere entrainment and particle sedimentation, which writes as: where ρ atm is the density of the atmosphere and U is the entrainment velocity defined as in Hewett et al. (1971): In the previous equation, α |U sc − U atm cos φ| is entrainment by radial inflow minus the amount swept tangentially along the plume margin by the wind, and γ |U atm sin φ| is entrainment from wind. Equation (22) states that the variation of mass flux 265 along the plume axis (left-hand side term) is due to ingestion of atmospheric air into the column (first right-hand side term) and loss of particles from the column margins (second right-hand side term).
The horizontal and vertical components of the momentum balance are derived from Newton's second law and the variation of mass flux as: and On the horizontal directions (Eqs. 23 and 24), the momentum of the eruptive mixture varies due to the wind (first right-hand side term) and the loss of particles (second right-hand side term), while on the vertical axis (Eq. 25) the momentum is affected 275 by the gravitational acceleration term and the fall-out of particles.
Mixture temperature T mix varies when plume rises and in this version of PLUME-MoM temperature, in a manner similar to that described in Folch et al. (2016), is computed from the total specific energy of the mixture E, defined as the sum of mixture specific entalphy, potential energy and kinetic energy as:

280
The specific enthalpy of the mixture H is written as the mass weighted average of the enthalpies of its components: where the terms on the right-hand side express the contributions to mixture enthalpy of dry air, solid particles, water vapour, liquid water, ice and other volcanic gases, respectively. All the mass fractions are denoted with the notation x (·) and they refer to the mass fraction of the component (·) with respect to the whole mixture. Similarly, the notation C (·) is used to represent the 285 specific heat capacities of the different components. In Eq. (27), h wv0 and h lw0 are the enthalpies at a reference temperature (T ref = 273.15 K) per unit mass of water vapour and liquid water, respectively.
By denoting with w atm the specific humidity of the atmosphere (mass of water vapour in a unit mass of moist air, expressed here as kilograms of vapour per kilogram of air), the following conservation equation for the total energy of the eruptive mixture can be written: On the right-hand side of Eq. (28) the terms responsible for the variation of mixture energy are listed. These terms include the interaction of the eruptive mixture with the atmosphere (which is assumed to contain water) and the loss of thermal and kinetic energy due to the sedimentation of solid particles.
Eq. (29) states that the entrainment of dry air into the column is responsible for the variation of dry air mass flux. Analogously, the transport of water is affected by the ingestion into the the column of atmospheric water as: where x w is the mass fraction of water (which can be in either vapour, liquid and ice form) in the mixture.

300
Finally, PLUME-MoM-TSM accounts for the presence of additional volcanic gases: From Eq. (31) we see that the volcanic gases mass flux is assumed to be constant along the plume. We are thus assuming that all volcanic gases reach the top of the plume (or the neutral buoyancy level) where they are released into the atmosphere.
The system of Eqs. (19-31), describing the steady dynamics of volcanic plume rise in the atmosphere, is solved numerically 305 by PLUME-MoM-TSM. A 4-5th order Dormand-Prince Runge-Kutta method (Dormand and Prince, 1980), which automatically adapt the integration step to control the numerical error, is implemented in a Fortran 90 code for the numerical integration of the system of ordinary differential equations. The computational time requested for a run is of the order of a few tenths of seconds for a simulation without aggregation, and slightly longer for simulations modeling also the aggregation process.

Phase changes of water 310
The present version of PLUME-MoM accounts for phase change of water as a function of the pressure and temperature conditions. Water (in vapour, liquid or ice phase) in the mixture can be either magmatic, added at the vent as liquid, or incorporated through ingestion of moist atmospheric air. Within the eruptive column, phase changes lead the water to be partitioned into vapour, liquid and ice: 315 where x w is the water mass fraction in the eruptive mixture, while x wv , x lw and x i are the mass fractions of water vapour, liquid water and ice phases, respectively.
In the following, we will use the subscript da, wv and vg to refer to dry air, water vapour and volcanic gases, respectively.
With this notation, assuming that the mixture of dry air, water vapour and volcanic gases behaves as an ideal gas, we can write the following relationships between atmospheric pressure P atm , phases partial pressures P (·) and molar fractions n (·) : 320 P atm = P da + P wv + P vg , P da = n da P atm , P wv = n wv P atm , Molar fractions of water vapour (n wv ), volcanic gases (n vg ) and dry air (n da ) in the gas phase of the volcanic mixture (for which it holds n wv + n da + n vg = 1) are given by: where mw wv = 0.018 kg mole −1 and mw da = 0.029 kg mole −1 are the molar weights of water vapour and dry air. The molar weight of the additional mixture of volcanic gases mw vg is: where x k and mw k are the mass fraction with respect to the mixture and the molar weight of the k−th volcanic gas. For CO 2 330 and SO 2 (most abundant volcanic gases after water), mw k are 0.044 kg mole −1 and 0.064 kg mole −1 respectively.
Once the conservation equation for the total energy (Eq. 28) is integrated, giving a new value of the total specific energy E, the specific enthalpy of the eruptive mixture H is computed from Eq. (26). Then, mixture temperature (T mix ) and water mass fractions (x wv , x lw , x i ) are computed from the known value of H by finding the values which satisfy Eq. (27) and Eqs. (32-36).

Addition of external liquid water
The addition of external liquid water to the volcanic mixture is one of the new features of PLUME-MoM-TSM. As in Koyaguchi and Woods (1996), we assume that magma-water mixing takes place at the vent location and that the thermal equilibrium between the eruptive mixture and the external water is reached before the interaction of the mixture with the atmosphere.
where x s and x wv are the mass fractions of particles and magmatic water vapour in the erupted magma and T mix0 is its temperature.
If external liquid water is added at the vent, the eruptive mixture formed by magma and liquid water is initialized by computing a new enthalpy (H ventnew ): where x lwext is the mass fraction of external liquid water at temperature T lwext , and x erupt = 1 − x lwext is the magmatic mass fraction. From the value of H ventnew , we update magma temperature and water partitions by following the procedure described in Appendix A.
To test the consistency of the results produced by PLUME-MoM-TSM in the presence of external water, in the Supplemen-350 tary Materials we reproduced some of the figures presented by Koyaguchi and Woods (1996).

Umbrella spreading model
As previously stated, PLUME-MoM-TSM includes the possibility to simulate the initial spreading of the umbrella cloud associated with a volcanic column. When this option is activated, the column equations are solved up to the neutral buoyancy level, where the density of the mixture equals that of the atmosphere. Above this height, the initial transport of the ash is 355 dominated by a radial spreading, due to the fact that the plume still has a vertical momentum, which produce an overshoot in the atmosphere and then an intrusion at the neutral buoyancy level as a gravity current. This dynamics is described by a different set of equations from those of the column, and here we use a system of "shallow-water" equations similar to that developed by Baines (2013); Johnson et al. (2015).
The intruding fluid is assumed to spread horizontally at the neutral buoyancy level as a shallow layer (with respect to the 360 atmosphere thickness) so that the vertical accelerations are negligible and the pressure distribution within the umbrella is hydrostatic. In addition, we assume that, due to turbulent mixing, the density of fluid is uniform within the umbrella, and this implies that variations in the thickness of the layer generate horizontal pressure gradients driving the motion, which is resisted by turbulent drag. Furthermore, for the umbrella cloud model, we neglect the loss of particles from its base, and thus the density can be considered constant. Differently from the model of Johnson et al. (2015), where the volumetric flux intruding from the 365 column into the atmosphere is provided by imposing a radial axysymmetric flux at a prescribed distance from the plume centerline (with fixed initial thickness), here the flux is provided, in both mass and momentum equations, through appropriate source terms derived from the neutral buoyancy level horizontal section of the column. In this way, it is possible to consider also the initial horizontal velocity of the umbrella cloud resulting from the horizontal momentum of the volcanic plume.
velocity, the equations for the conservation of volume and horizontal momentum write as: On the left-hand sides of the momentum equations, N is the Brunt-Väisälä frequency, which represents and upper limit for the frequency of internal gravity waves in the atmosphere. Its value is computed at the neutral buoyancy level from the atmospheric density profile as follows: On the right-hand size of all the equations in system (40),w is the volumetric flux source derived from the vertical velocity at the vent: and γ x and γ y are two functions which take into account the vertical gradient of the plume radius at the neutral buoyancy level: In this way, the volume and the momentum flux derived from the plume model are correctly passed to the umbrella model.
The system of partial differential equations is solved with a transient finite-volume code based on the numerical solver Following Reckziegel et al. (2019), the total grain size distribution released at the vent for both the phases has been initialized with of a bi-modal distribution with 11 particle classes from −1φ to 9φ, with densities increasing linearly with φ from 2300 410 to 2700 kg m −3 . Particle sphericity, mostly affecting the loss of particles through the terminal velocity equation (Textor et al., 2006b), has been fixed to 0.8. Following Van Eaton et al. (2016), we assumed 5 wt% of magmatic water, plus an additional 5 wt% of cooler surface water, possibly associated with the presence of a summit glacier prior to the eruption. The initial temperature of the magmatic mixture has been fixed to 1173 K and, after the mixing with the external water, a temperature of 1006 K has been found for the base of the plume. The initial plume velocity for both the phases has been fixed to 250 m s −1 , 415 and the radius has been changed for the two phases in order to match the height of the intrusion of the umbrella cloud in the atmosphere and its main spreading direction. This results in a mass eruption rate of 9.54×10 6 kg s −1 for the first phase and 2.02×10 7 kg s −1 for the second one. Varying the initial velocity, while keeping fixed the height of the neutral buoyancy level, produced very similar results in terms of mass eruption rates.
In level is relatively small, ranging from a maximum of 0.24 for φ = −1 to fractions smaller that 0.01 for φ ≥ 3. This means that almost all the particles erupted reach the neutral buoyancy level, and that the initial grain size distribution differs very little 430 from that feeding the umbrella cloud (bottom-left panel). In any case, almost immediately above the vent, due to the large entrainment, the mass of particles, as that of water, become rapidly a small fraction of the total mass of the mixture, with the dry air being the more abundant phase (middle-right panel). At the neutral buoyancy level, reached at 12.7 km above sea level, the plume radius is 3935 m (see bottom-right panel) and the mass flow rate injected in the umbrella cloud is 6. atmosphere observed during the eruption. Three values were tested (0.5, 1 and 2) by carrying on the simulations for 1.5 hour and 2 hours for the first and second phase, respectively. At each time step we computed: (i) the upwind spreading, defined 445 here as the maximum horizontal distance from the vent in the opposite direction to that of the wind at the neutral buoyancy level; (ii) the radius of a circle with area equivalent to that of the modeled umbrella cloud. These values are plotted, at intervals of 300 s, in the left panels of Fig. 7, with red lines for the upwind distances and blue lines for the equivalent radius. First of all, we can see from the plots that for all values of C D a steady upwind distance is reached at the end of the simulation, and that by changing (increasing or decreasing) the drag coefficient of a factor 2 with respect to the value C D = 1, the upwind C D = 1 seems to better match the observation, with a small overestimate at the beginning and a small underestimate at the end of the simulation. For this reason, in the following of the paper, we used this value as reference value.

455
In the right panels of Fig. 7  star in the right panels of Fig. 7) are produced also for subplinian eruptions. For the simulation presented in this work, almost 90% of the maximum upwind spreading is reached in 30 minutes after the onset of the intrusion at the neutral buoyancy level.
In the analysis presented so far, we have not accounted for the effect of aggregation. To better understand if aggregation 465 could affect the spreading of the umbrella for this scenario, we repeated the plume rise simulation for the first phase of the April 2015 Calbuco eruption by enabling aggregation between particles and analyzing the changes in the height of the neutral buoyancy level and the grain size distribution of particles still present in the plume and released from its margin at that height.
For this set of simulations, two coarser bins with an initial negligible amount of particles have been added to the total grain size distribution. This allow for the formation of aggregates larger than the particles present in the original size spectrum. Density 470 of the aggregates has been fixed to 1500 kg m −3 , resulting in lower settling velocity than the original particles, the size begin equal.
https://doi.org/10.5194/gmd-2020-227 Preprint. Discussion started: 2 October 2020 c Author(s) 2020. CC BY 4.0 License. it is limited to the coarser particles. This is mostly due to the fact that the dominant collisional mechanism is the differential settling velocity, which is larger for coarser particles. It is important to note that the neutral buoyancy level differs from that obtained without aggregation by less than 0.1%, and also the total amount of particles lost from the plume does not differ 490 significantly.
In order to better understand the effect of aggregation, we performed additional simulations with a constant kernel, starting from a value β = 10 −15 m 3 s −1 , up to a value β = 10 −13 m 3 s −1 . The simulation with β = 10 −15 m 3 s −1 shows small differences with that without aggregation, with small amounts of aggregated particles produced (represented in panel B1 of Fig. 8 with the red bars). Because of the adoption of a constant kernel, aggregates form in each bin of the initial grain size distribution.

495
By increasing β of an order of magnitude to a value of 10 −14 m 3 s −1 , the amount of aggregates become significant (more than 50% for some bins), and also the relative amounts of particles released from the plume margins at the neutral buoyancy level differs from those obtained without aggregation. Despite that, we observe that the height of the neutral buoyancy level of the plume does not change, and also the overall bi-modal shape of the total (aggregated and non-aggregated particles) distribution does not change substantially. This is because, with the binary aggregation and the constant aggregation kernel considered 500 here, most of the aggregation occurs "intra-bin", i.e. produce an aggregate with a volume, or mass, which belongs to the same bin of particles forming the new aggregate. In fact, for bin sizes of 1 phi, and a uniform mass distribution within each bin, the number of smallest particles within a bin is almost one order of magnitude smaller than the number of largest particles.
Thus, with a constant kernel, aggregation of the smallest particles occurs much frequently, and the aggregate produced from a binary aggregation of these particles produce a new particle with an equivalent diameter (diameter calculated from the mass) 505 still belonging to the same bin. In order to produce a larger change in the particle size distribution in the plume at the neutral buoyancy level, we have to increase the constant kernel to β = 10 −13 m 3 s −1 , for which the fraction of aggregates is close to 97% of the total mass of particles in the plume (Fig. 8, panels E1 and E2). It is interesting to observe that even with such a large amount of aggregates, the total amount of particles lost during the rise does not change with respect to the simulation without aggregation, and also the height of the neutral buoyancy level differs from that obtained without aggregation only by 1 m. This 510 result allows us to disregard aggregation for the sensitivity analysis presented in the next section, where the main outputs of the plume module used to initialize the umbrella spreading module are the volumetric flow rate and the height of the gas/particles mixture injected by the plume into the atmosphere at the neutral buoyancy level.

Umbrella cloud upwind spreading
In this section we present a sensitivity analysis of the upwind spreading of the umbrella cloud with respect to mass flow rate and 515 to some parameters characterizing the meteorological profile above the vent. This study builds on the ideas of Carey and Sparks (1986), where the ratio of upwind spreading to the amount of axial displacement at the neutral buoyancy level was analyzed, with results based on theoretical and empirical information on column behaviour. The simulations presented here employ pressure, temperature and density atmospheric profiles modified from the International Standard Atmosphere (ISA) model.
This model provides a common reference for temperature and pressure at various altitudes, and it is based on the subdivision 520 of the atmosphere in layers where the absolute temperature T atm against geopotential altitude h is assumed to vary linearly.
Geopotential and orthometric heights differ little for low terrestrial atmosphere (Daidzic, 2019), and for this reason we can assume a linear dependence of temperature also against orthometric height. In the following, as common in meteorology, we refer to the negative of the rate of temperature change with altitude with the term "lapse rate", and we use for it the notation Γ.
Once the bottom values and the lapse rate are known in an atmospheric layer, pressure, temperature and density profiles within 525 the layer can be calculated by solving simultaneously the equation for vertical pressure gradient resulting from hydrostatic balance, the equation for temperature and the ideal gas law: When the option for a modified standard atmosphere is selected in PLUME-MoM-TSM, the first two equations are added to the set of ordinary differential equations solved numerically by the model, and the last one is employed to calculate explicitly 530 the right-hand term of the first equation.
ISA lapse rates do not account for humidity effects and air is assumed to be dry. For this reason, when humidity is considered, appropriate corrections to the atmospheric profiles given by the ISA model need to be introduced. Moist unsaturated air cools with height at a slightly lower rate than dry air (Dutton, 2002) and in this work we adopt the following lapse rate correction to account for the presence of water vapour, as proposed in Daidzic (2019): where Γ atm and Γ d are the moist and dry atmospheric lapse rates, respectively. According to Ahrens (2012), the average specific humidity in the lower troposphere ranges from 0.004 kg kg −1 at 60 degrees (north or south) to 0.018 kg kg −1 at the equator, and its value exponentially decrease to values lower than 10 −5 kg kg −1 above the tropopause. It follows that the sensitivity of lapse rates to the specific humidity value is minimal in the lower troposphere, and negligible above it.

540
Density of moist atmosphere also varies with water content, and the relationship between dry and moist air densities can be written as: where R wv and R air are the gas constant of water vapour and air, respectively.
With these equations, ISA lapse rates Γ d and density ρ d are corrected for the presence of humidity and it is possible to 545 solve the system of Eqs. (43) for the modified pressure and temperature profiles. The specific humidity profile used in the simulations presented here has been derived from the data presented in Anderson et al. (1986), providing the values at the boundaries between the atmospheric layers, while the troposhpere bottom value has been varied in the range 0.004-0.018 kg kg −1 . Inside each layer, the logarithm of specific humidity is assumed to vary linearly.
In this work we also assume the wind speed being a continuous function of height which varies linearly in several layers.

550
Observation indicate that on average wind speed profile has local minima at altitudes of about 20 km and about 90 km, where air masses are relatively stationary, and local maxima at altitudes of about 11 km and 70 km (Modica et al., 2007;Brasefield, 1954). While a quite wide range of values is observed for the first maximum at the tropopause, the value of the wind speed at 70 km varies less, and a speed of 65 m s −1 is assumed in this work. For the wind speed minimum at 20 km altitude, we fix a value of 10 m s −1 . With these assumptions, the relevant atmospheric profiles obtained for a moist atmosphere and with a wind 555 speed of 20 m s −1 at 11 km are plotted in Fig. 9, where the red markers in the last two panels denote the values we varied for the sensitivity analysis.  Figure 9. Atmospheric profiles for a dry atmosphere and with a 20 m s −1 wind speed at the tropopause (11 km asl). The red markers denote the values that can be changed in the modified standard atmosphere adopted in the model. The first three profiles (pressure, temperature and density) also change when the specific humidity at the sea level is changed.
Maximum upwind spreading of the umbrella cloud with respect to vent location represents an important information that coupled plume and umbrella models, or fully three-dimensional models (Cerminara et al., 2015;Suzuki and Koyaguchi, 2009) can provide. It has been noted that using only the output of integral plume models as source of VATD models generally 560 underestimate the upwind depositions patterns of volcanic ash, and corrections accounting for the initial spreading of the umbrella as a gravity current result in a better agreement of ash distribution close to the eruptive source (Costa et al., 2013;Webster et al., 2020).
Most of the work done so far to couple the spreading of the umbrella cloud with advection-diffusion-sedimentation models focused on large events as the 1991 Pinatubo eruption, with a mass flow rate of the order of 10 9 kg s −1 . Less work has 565 been devoted to smaller scenarios and for this reason here we focus on lower values of mass flow rate, in order to gain a better understanding of the relationship between mass flow rate, wind speed and upwind spreading. It is well known that upwind spreading increases with increasing mass flow rate and decreasing wind, but a quantitative relationship between them is unknown. Here, we varied three input parameters of the model: 1) the decadic logarithm of mass flow rate in the range 6 − 8 (i.e.ṁ = 10 6 − 10 8 kgs −1 ); 2) the wind speed at the tropopause in the range 35 − 80 m s −1 ; 3) the specific humidity at the 570 sea level in the range 0.004-0.018 kg kg −1 . A uniform probability has been assigned to all the input variables and the space of parameters has been sampled with a Latin Hypercube Sampling generating 700 elements. The other main input parameters used for this analysis are listed in Table 1. As previously stated, aggregation is disregarded here, because the effect on the input parameters of the umbrella module is negligible. In addition, as shown in de' Michieli Vitturi et al. (2015), also the initial total grain size distribution has a minimal control on plume height, and consequently on upwind umbrella spreading. For this reason, it has been kept fixed for all the simulations. For each simulation, the horizontal upwind distance from the vent of the umbrella cloud has been calculated from the ouptut of umbrella spreading module (when a steady value is reached), and correlation plots with the three input parameters are presented in the first three panels of Fig. 10, where each red dot represents a different simulation. In general, we observe increasing upwind spreading with increasing mass eruption rate and decreasing wind speed at the tropopause, while the relative 580 humidity does not seems to affect it. For a fixed value of the input parameters on the x-axis, the vertical spreading of the red dot cloud is associated with variations in the other input parameters. Each black line represents the average of the upwind spreading values at a given input value. In the first three panels each red dot represents a different simulation, whereas the black lines represent the mean of the output parameter at a given input value. For a fixed value of the input parameters on the x-axis of a plot, the vertical spreading of the results is associated with variations in the other input parameters.
The global sensitivity of the upwind spreading with respect to the three input parameters has been computed by using a variance based method (Saltelli et al., 2010;Scollo et al., 2008;de' Michieli Vitturi et al., 2016). If we denote with Y the 585 output parameter of interest and with x i the input parameters considered, the main sensitivity indexes S i , also called the Sobol indices, express the fraction of the variability in the output Y , that is reduced when the value of the input x i is fixed. These indices are calculated by: The value of the Sobol indices are shown in the bar plot presented in the right panel of Fig. 10. The largest value is obtained 590 for the sensitivity with respect to the mass flow rate (yellow bar), with a smaller sensitivity with respect to the wind speed at the tropopause (red bar) and a negligible effect of the specific humidity at sea level (blue bar). This highlight the fact the upwind spreading is mostly controlled by the first two parameters. For this reason, we plotted in Fig. 11 the maximum upwind spreading as a function of the logarithm of the mass flow rate and the tropopause wind speed only, highlighting an exponential dependence of the spreading on a linear combination of the two 595 input variables. In fact, if we denote with d up the upwind spreading and withṁ and v the mass flow rate and the tropopause wind speed, respectively, we can write: A least squares minimization procedure gives as optimal values for the fitting parameters a = 0.1495 m, b = −1.601 and θ = −1.066 • , with a coefficient of determination R 2 = 0.997 and a square root of the variance of the residuals of 173.1 meters.

600
With these values, the relationship expressed by Eq. (47) gives an estimation of the upwind spreading once the mass flow rate and the tropopause wind speed are known. While from a modeling perspective this represents an interesting result, it is not always applicable, because the input values needed to compute the spreading cannot be observed directly. For this reason, we plotted in the right panel of Fig. 11 the upwind umbrella distance as a function of the height of the NBL and the downwind horizontal distance from the vent of the centerline, both extracted from the output of PlumeMoM-TSM plume module. As 605 expected, upwind spreading decreases with wincreasing bending of the column (expressed by the downwind distance of the centerline from the wind), because the bending is associated with wind intensity. On the contrary, upwind spreading increases with increasing height of the intrusion in the atmosphere. Also in this case there is an exponential dependence of the upwind spreading on a linear combination of the two input variables and, if we denote with h and d down the NBL height and the downwind distance, respectively, we can write: In this case, the least squares minimization procedure gives as optimal values for the fitting parameters This second relationship, from an operational point of view, results more useful. On one hand, the two input values are produced by all the plume integral model, and in this way a correction can be applied to obtain the upwind spreading without the need of an umbrella cloud spreading model. On the other hand, the input values are also observable data than can be estimated with monitoring techniques, allowing to have an immediate initial estimation of the upwind spreading, and thus of 620 the upwind area potentially affected by the presence of ash, once the volcanic column is observed.
The results presented here are consistent with those discussed in Carey and Sparks (1986), who showed that for very weak plumes no true umbrella cloud forms, and a direct coupling of column and dispersal models is appropriate, and that for column able to reach the tropopause upwind spreading is significant and an umbrella model is required. Finally, it is worth to note that in this analysis we have not varied some of the parameters which can have a large control on the eruptive column, such as 625 the initial water content and temperature (de' Michieli Vitturi et al., 2016;Costa et al., 2016). Further studies will address this point, by considering the uncertainty associated to additional input parameters of the models. of turbulent volcanic plumes. Despite that, integral models are able to predict column heights and neutral buoyancy levels consistent with those calculated by 3D models . This result, coupled with the fact the integral models require few seconds to run, pushes toward a further development of these models. In this work we have presented PLUME-MoM-TSM, a new version of the volcanic plume model PLUME-MoM described in de' Michieli Vitturi et al. (2015) and based on the method of moments for the description of the particles distribution. With respect to the previous version, where a finite 635 set of moments of the total grain size distribution were tracked, here the adoption of the two-size moment allowed to model the evolution of the mass and the number of particles associated to each single bin of the distribution. A first important consequence is that, with such approach, it is possible to use PLUME-MoM-TSM outputs as input of VATD models. In addition, the method used in this work to describe the grain size distribution allows for the solution of the full Smoluchowski coagulation equation, and thus to model accurately particle aggregation. We applied here the two-size moment method to the steady-state integral 640 equations for the plume, but the procedure described is general and the code is based on modular structure, which make easy to port the procedure to other Eulerian models for transport of volcanic ash in the atmosphere. PLUME-MoM-TSM also accounts for phase change of water, resulting in the formation of a liquid water or ice inside the plume, and it includes a module for the spreading of the umbrella cloud as a gravity current. These new features have been tested by applying the model to the April 2015 eruption of Calbuco volcano in Chile, allowing to constrain model parameters 645 and to quantify the effect of wet aggregation on model results. The analysis shows that aggregation has a minimal control on plume characteristics and on the loss of particles from its margins.
Finally, a sensitivity analysis provided a relationship between mass flow rate, wind speed at the tropopause and maximum upwind spreading of the umbrella cloud, and also between observable parameters quantifying plume bending (height of the neutral buoyancy level and downwind displacement of the plume axis at the same height) and the umbrella upwind spreading.
These results are in agreement and extend those presented in Carey and Sparks (1986), and provide important information for both modelers and observatories, even if further studies are necessaries to extend the analysis to more general initial vent conditions. 1. For T mix ≥ T ref , water can be in vapour and liquid form only and T mix is given by: At equilibrium condition, when both liquid water and vapour are present, it must hold: where both the terms can be written as functions of x wv , Eqs. (A1), is ≥ T ref . If this is not the case, T mix must be less than T ref .
2. For T mix < T ref − 40, water can be in vapour and ice form only and the temperature of the eruptive mixture can be computed as: We follow the procedure adopted for T mix ≥ T At equilibrium, when both vapour and ice are present, it must hold: where P wv and e s are functions of x wv only, Eqs.
The right-hand side term contains n wv which, at this stage, is unknown. However, at equilibrium conditions, n wv is a function of T mix and it can be computed by inverting Eq. (33) as: where e is vapour partial pressure which, depending on mixture temperature, can be evaluated from Eq. (A2) for T mix = As done in Mastin (2007), we assume that liquid water mass fraction is a function of T mix as: The specific enthalpy of the eruptive mixture at equilibrium (H eq ) is known by solving the energy equation, and, for a generic mixture temperatureT , it must hold: where H| Tmix=T is the mixture enthalpy at the generic temperatureT .

A2 Settling velocity
Several models are implemented in PLUME-MoM-TSM to compute particle settling velocity (w ip ) as a function of particle diameter (D), particle density (ρ p ) and shape factor (ψ).
The second model (Ganser, 1993) defines the settling velocity for the Stokes' regime (Re << 0.005) as: The drag coefficient (C D ) in Eq.A13 is calculated with the following expression: where Re is the Reynolds number, K 1 is the Stokes' shape factor and K 2 is the Newton's shape factor: Finally, settling velocity can be computed following Pfeiffer et al. (2005) as: where m is the particle mass and A cs the cross-sectional area of the particle. Depending on the Reynolds number, the drag coefficient used in Eq. (A18) is: Re ψ −0.828 + 2 √ 1.07 − ψ, Re ≤ 100.
(A19) 740 We remark here that the Reynolds number is a function of the settling velocity, and thus for the formulations where the drag coefficient is a function of the Reynolds number a nonlinear equation has to be solved. This is done in PLUME-MoM-TSM with an iterative procedure.
in PLUME-MoM-TSM. In their work, a general model describing the aggregation processes was presented, but a simplified version computationally cheaper was then adopted. Here, the full original version is retained and implemented in the plume model, according with the framework of the TSM method of moments.
The kernel aggregation is defined as the product of a collision frequency function δ ij and a sticking efficiency function λ ij , with the collision frequency function given by the sum of three contributions: a Brownian frequency function δ B ; a 750 fluid shear frequency function δ S ; and a differential sedimentation frequency function δ DS . Considering two particles of the family i p with masses m i and m j and diameters in meters d i and d j , we can write where k b is the Boltzmann constant, µ is the dynamic viscosity of air, and Γ S is the fluid shear.
Following Costa et al. (2010), sticking efficiency λ ij in presence of ice only (no liquid water) takes a constant value: When liquid water is present in the plume (without ice), sticking is a function of the viscous Stokes number St ij (Costa et al., 2010;Liu and Litster, 2002), and the following parametrization is adopted: where St cr = 1.3 is a critical Stokes number above which there is no sticking, and q = 0.8 is an empirical parameter.

760
These two expressions for sticking efficiency do not take into account the amount of ice/liquid water available in the plume as limiting factor for aggregation. Here, we modified the formulation presented in Costa et al. (2010), by taking the sticking efficiency as a weighted linear combination of the two expressions above: where α lw and α i are the volumetric fraction of liquid water and ice in the mixture, respectively.
C4 of Mastin (2007)) and those of Koyaguchi and Woods (1996) (green lines digitized from Fig. 5). It can be seen a general agreement between our results and those of the other models, with a maximum difference of 20% between the maximum column height found by PLUME-MoM-TSM and the one found by Plumeria for the 10% external water case. For the 0% case, PLUME-MoM-TSM predicts column collapse for a slightly lower mass eruption rate than the other models, while a good fit between PLUME-MoM-TSM and Koyaguchi and Woods (1996) results is observed for the other cases.  Koyaguchi and Woods (1996) and Mastin (2007), input parameters are: magma temperature = 1000 K, water vapour mass fraction = 0.03 wt and relative humidity in the atmosphere = 100%. The exit velocity for the dry scenario is 100 m s −1 , while an adjustment on exit velocities was done for the runs with external water to match the mass flux of magma and gas of the dry case. For this reason, exit velocities are 300, 329, 289 and 236 m s −1 for the runs involving 10%, 20%, 30% and 40% of external water, respectively (Mastin, 2007). RIGHT PANEL: plume height as a function of mass fraction of external water added at the vent. Orange lines show the PLUME-MoM-TSM results, while light blue lines are digitized from Fig. 4 of Koyaguchi and Woods (1996). Curves are shown for vent radii of 10, 50, 100, 500 and 100 m. Exit velocity is equal to 100 m s −1 , water vapour content to 0.03 wt and mixture temperature to 1000 K.
The right panel of Fig. B1 reproduces Fig. 4 of Koyaguchi and Woods (1996), where it is shown the height of the eruptive column as a function of the mass of external water added at the vent. Each curve is drawn for a constant mass flow rate of 2.24 x 10 5 , 5.6 x 10 6 , 2.24 x 10 7 , 5.6 x 10 8 and 2.24 x 10 9 kg s −1 . These values correspond to vent radii of 10, 50, 100, 500 and 1000 m, respectively. The PLUME-MoM-TSM results (orange lines) appear to be similar to those of (Koyaguchi and Woods, 1996) (light blue lines), although PLUME-MoM-TSM predicts column collapse for slightly lower values of external water than 780 Koyaguchi and Woods (1996). Fig. B2 and B3 present the variation of mixture temperature, velocity, relative density, gas mass fraction and water mass fraction (liquid and ice) as functions of column height for two eruptive scenarios. Each panel shows the results for amounts of external water of 3, 10, 20 and 30 wt%. The two figures reproduce Fig. 3 of Koyaguchi and Woods (1996).  Figure B2. Variation with column height of mixture temperature, velocity, relative density, gas mass fraction and liquid+ice mass fraction.
As done in Koyaguchi and Woods (1996), each curve accounts for a different amount of external water added to the magmatic mixture: 3, 10, 20, 30 wt%. The dry case is set for a vent radius of 50 m and an exit velocity of 100 m s −1 , resulting in a mass flow rate of 5.7 x 10 6 kg s −1 . The wet cases are obtained by keeping constant vent radius and mass flow rate, while varying exit velocity. Magma is assumed to contain 0.03 wt fraction of volatiles (only water vapour) and to have an initial temperature of 1000 K. Following Koyaguchi and Woods (1996), Fig. B4 shows the variation of initial mixture density and temperature as the mass 785 of external liquid water increases. The figure shows how the addition of a certain amount of external water alters the properties of the magmatic mixture at the vent. The three curves are done for different initial temperatures of the eruptive mixture (1000 and 1400 K) and volatile contents (0.01 and 0.05 wt fraction, only water vapour is considered). The temperature of the external water is 273 K.  Figure B4. Variation of initial mixture density (ρmix) and temperature (Tmix) as the mass of water added to the magma is increased.
Continuous lines show the results of PLUME-MoM, while dotted lines are digitized from Koyaguchi and Woods (1996). Curves are shown for initial magma temperatures of 1000 and 1400 K, while magma volatile content ranges from of 0.01 to 0.05 wt fraction.
Appendix C: List of model variables 790