PLUME-MoM-TSM 1.0.0: a volcanic column and umbrella cloud spreading model

. In this paper, we present a new version of PLUME-MoM, a 1-D integral volcanic plume model based on the method of moments for the description of the polydispersity in solid particles. The model describes the steady-state dynamics of a plume in a 3-D coordinate system, and a modiﬁcation of 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 ﬁrst time in a plume model, the full Smoluchowski coagulation equation is solved, allowing us to quantify the formation of aggregates during the rise of the plume. In addition, PLUME-MOM-TSM allows us 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 computation of the upwind spreading. The new model is applied ﬁrst to the eruption of the Calbuco volcano in southern Chile in April 2015 and then to a sensitivity analysis of the upwind spreading of the umbrella cloud to mass ﬂow rate and meteorological conditions (wind speed and humidity). This analysis provides an analytical relationship between the upwind spreading and some observable 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.


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 hazards associated with explosive eruptions, affecting both the proximity of the volcanic vent, where it can damage vegetation and infrastructures, and pose a health risk to people, and the surrounding atmosphere, requiring effective and prompt actions to regulate flights near the ash plume (Bonadonna et al., 2021).
In order to properly model the transport of ash and its associated deposition on the ground, it is important to provide the correct inputs to volcanic-ash transport and dispersion (VATD) models, and so far the best approach has been to use the output of numerical models simulating the rise of volcanic columns up to the neutral buoyancy level, where ash clouds transported by the wind spread in the atmosphere. Column models can be broadly categorized in two classes: one-dimensional (1-D) integral models (Bursik, 2001;Mastin, 2007;Degruyter and Bonadonna, 2012;Devenish, 2013;Woodhouse et al., 2013;de' Michieli Vitturi et al., 2015;Folch et al., 2016) and three-dimensional (3-D) models Suzuki and Koyaguchi, 2009;Cerminara et al., 2016;Ongaro et al., 2007). A recent intercomparison exercise (Costa et al., 2016) has shown that 1-D and 3-D models, notwithstanding that the formulations are different, produce consistent predictions of some of the more relevant quantities describing volcanic plumes, as, for example, the neutral buoyancy level. For this reason, in operative contexts where multiple scenarios have to be considered and the simulation time is an important constraint, 1-D mod-els are still preferred. This class of models is based on the theory of turbulent buoyant plumes, as introduced by Morton et al. (1956), and it has been first applied in the volcanic context by Sparks (1986) and Woods (1988), who modeled the rising plume as a homogeneous mixture of gas and ash particles. Later works incorporated additional processes, such as the effects of moisture (Woods, 1993;Mastin, 2007) and ambient wind (Bursik, 2001).
Despite the relatively high number of integral models developed in the past couple of decades, some important features of volcanic plumes have been neglected so far, or they have not been considered together in a single model. Most notably, only a few models account for the thermodynamic phase relations for water . As pointed out in several works (Koyaguchi and Woods, 1996;Woods, 1993;Sparks et al., 1997;Mastin, 2007), the presence of water, coming both from entrainment during some Surtseyan eruptions (Þórarinsson, 1967) or from entrainment of water-saturated air (Graf et al., 1999), and the thermodynamics of condensation and evaporation, can strongly affect the rise and collapse of strong volcanic plumes. Conversely, in Woodhouse et al. (2013), it is shown that for small weak plumes, the effect on height of external humidity entrained in the eruption column is almost negligible when compared to the role played by the uncertainty in parameters like source buoyancy flux, atmospheric stratification, and wind.
Another process that is frequently disregarded in the modeling of volcanic plumes is solid particle aggregation (Brown et al., 2012). So far, the number of works quantitatively describing aggregation processes in an eruption column is limited (Veitch and Woods, 2001;Textor et al., 2006b, a;Brown et al., 2012;Folch et al., 2016;Künzli et al., 2018;Rossi et al., 2018). This is because this process is quite complex, from both physical and modeling perspectives, requiring a careful treatment of the evolving grain size distribution. In fact, aggregation leads to the formation, from fine ash, of coarser particles, resulting in larger settling velocities than the primary particles and thus with an increased fallout of fine ash from the plume margins. A proper treatment of aggregation is thus important to predict the effective grain size distribution of solid particles reaching the neutral buoyancy level and then advected by the wind into the atmosphere.
For sustained eruptions, when neutral buoyancy level is reached, the eruptive mixture still has a vertical momentum, allowing the column to further rise and then to collapse as a gravity current spreading into the atmosphere (see Fig. 1). During the eruption of Mount St. Helens on 18 May 1980, this phenomenon was clear, and satellite data were used by Sparks (1986) to validate a simple model of umbrella spreading based on continuity of mass flow rate. For large mass eruption rates, this spreading can lead to an important radial expansion, including upwind transport, which is sometimes omitted in numerical simulations making use of coupled column and 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 advectiondiffusion-sedimentation model, and 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 behavior of the ash cloud in the proximity 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 atmosphere. The results of Johnson et al. (2015) show that the current spreads in a different way from what was suggested by simple scaling arguments in previous works, such as that of Costa et al. (2013). The same approach is used in Pouget et al. (2013Pouget et al. ( , 2016, where the results of a "shallow-water" model are compared with satellite observations in order to estimate growth rates and mass fluxes. More recently, Folch et al. (2016) adopted a simple semi-empirical model to describe the umbrella region, and their results are qualitatively consistent with the results obtained with more complex 3-D numerical codes (Costa et al., 2016). A simple approach is also presented in Webster et al. (2020), where a parameterization is used to model the lateral spread of the umbrella and to correct the initial conditions for transport and dispersion in the Lagrangian model NAME, resulting in substantial improvements to ash cloud predictions for large explosive eruptions.
In this work, we extend the original PLUME-MoM model, which is the first volcanic plume model based on the method of moments for the description of the polydispersity in solid particles, by adding all the important features described above: phase transition of water, particle aggregation, and umbrella cloud spreading. The model is first applied to the April 2015 Calbuco eruption in Sect. 3.1, by comparing simulation results with observations to constrain some model parameters. Then, in Sect. 3.2, a sensitivity analysis is performed with the newly developed model in order to obtain a relationship first between eruption parameters and the upwind spreading of the umbrella cloud first, and then between the latter and other parameters characterizing the plume.

Model equations
In this section, we present the equations solved by PLUME-MoM-TSM. First of all, the new formulation for the description of the particle size distribution, based on the method of moments, is presented in Sect. 2.1. Then, in Sect. 2.2, the steady transport equations for the plume are derived from conservation principles (mass, momentum, and energy). The equations for the phase transition of water and the shallowwater equations describing the spreading of the umbrella cloud are finally presented in Sect. 2.3 and 2.4, respectively.

Solid particle number density function
As previously stated, particle aggregation is a complex process and its modeling requires an accurate description of the particle size distribution within a column and of its changes during the rise of the volcanic mixture. While particle fallout from the plume margin changes the size distribution but does not modify the size spectrum, aggregation leads to the formation of coarser particles, whose size was not present in the initial size distribution erupted at the vent. Furthermore, when a discretization in bins is employed to describe the variability in size, computational problems arise in the modeling of aggregation processes, because particles belonging to two bins of the original distribution could form an aggregate whose size is not represented in the chosen discretization. In this work, we try to solve these problems in the description and modeling of polydispersity in size of solid particles by adopting the method of moments.
In the first version of PLUME-MoM (de ' Michieli Vitturi et al., 2015), the quadrature method of moments (QMOM) was used to model the polydispersity of tephra. According to the QMOM (Marchisio et al., 2003;Marchisio and Fox, 2013;Colucci et al., 2017), the particle size distribution is approximated from the moments as a sum of Dirac delta functions, which is particularly convenient for evaluating the integrals appearing in the closure terms of the moments transport equations. When aggregation processes are considered, QMOM provides an effective alternative to the sectional (or "discrete") method, where the size spectrum is partitioned in discrete bins. Its advantages are a limited number of variables and the dynamic calculation of the more representative sizes within the distribution. The main disadvantage is that, given the moments of a distribution, it is not possible to retrieve the amount of mass associated with 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 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 engineering contexts, 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 first two moments (Nguyen et al., 2016), leading also to a good representation of the NDF with a small number of sections. For this reason, TSM is a good candidate for the modeling of the aggregation processes occurring 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 size (described by the volume or the radius) of the particles and in each section two moments of the size distribution are used to reconstruct the NDF with a linear function of the size (hence the name "two-size moment"), here the NDF is defined as a function of particle mass. This greatly simplifies 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, 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), we switch here from a linear reconstruction to a power law in these sections, with zero value at one of the boundaries of the section. The advantages of this approach will be explained in detail in the following.
The starting point of the TSM method is the partition of the size spectrum into n sections (or bins). Here, given a maximum size φ 0 and a size increment φ, we define a uniform partition φ i = φ 0 + i φ, i = 0, . . ., n s , in the nondimensional 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 the original partition φ i the corresponding particle mass m i and thus to obtain a partition of the particles' mass spectrum in n s sections corresponding to the n s 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 solid particles' number density function (NDF), here denoted by η(m, x, t). This function, which is formally a mass distribution, represents the number of solid 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 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 . We introduce the following notation for the j th moment of the number density function η: We observe that in the notation of the moment on the lefthand side of the equation we dropped the explicit dependence from x and t. Also in the rest of the paper, for the sake of simplicity, we will use without ambiguity the notation without x and t to denote the moments, even when they vary in space and time.
One important property of the NDF that is adopted is that its moments 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 0, m (0) , represents the total number of solid particles per unit volume of the mixture, while the moment of order 1, m (1) , represents the total mass of solid particles per unit volume of the gas-particle mixture, i.e., the bulk density of the solid particles, here denoted by ρ 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 solid 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 orders 0 and 1 on the kth 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 1 M k 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 appropriate transport equation for the moments M k allows one to describe the evolution of the total grain size distribution within the plume. We observe that, conversely to the global moment of order 1, the local moments of order 1 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 particle mass, as, for example, the density, or the settling velocity, can be defined in terms of the continuous number density function η(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, ψ k represents the average of ψ in the mass interval [m k−1 , m k ] weighted by the number of solid particles, while ψ (1) k represents a mass-weighted average in the same interval.
When the moments of the number density function η are known for each section but not the actual expression of η, the integral in Eq. (4) must be approximated with a quadrature formula. 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 where the ith node x i is the ith 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 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 wherẽ With this approach, when the integrand of Eq. (4) changes (for example, because the number density function, i.e., the mass distribution, changes in space or time), there is no need to re-evaluate the quadrature nodes and weights, which depend 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), the number density function η(m) is approximated in each section by a linear function. Here, non-linear approximations ≈η k (m) of the number density function 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 havẽ In this way, the domain of each functionη k (m) is always the whole kth section, and there is no need to update the quadrature nodes. Further details on the selection of the case and the computation of the parameters α k and β k are given in Appendix A. 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 i p is used to distinguish the different number density functions η i p (m), i p = 1, . . ., n p , where n p is the number of particles families. Similarly, the notation m (j ) i p is introduced for the j th moment of the i p th number density function, and the notations N i p ,k and M i p ,k for the number and mass of particles of the i p th family in the kth 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 through the atmosphere. The model we present, based on the buoyant plume theory of Morton et al. (1956), is a 1-D integral model where plume properties are averaged over cross-sections. As in Bursik (2001), the model assumes a homogeneous mixture of different phases with thermal and mechanical equilibrium among all of them. In this version of the model, in addition to air and solid 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, such as CO 2 and SO 2 (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 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 u, v, and w as the three components of plume velocity, and 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 the following relationship holds: We also consider the presence of a 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 sections, which defines the plume centerline, can be bent, and the following equations hold for the centerline coordinates x s , y s , z s : When more than one particle family is considered, they are denoted by the subscript i p = 1, . . ., n p ; in addition, if aggregation is 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 a 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. Equation (12) describes the variation with the vertical coordinate z of the mass flux of particles of size m, due to fallout of tephra (first term on the right-hand side) and aggregation (second term on the right-hand side, S i p ). Here, w i p (·) is the particle settling velocity (see Appendix B2) 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 (Bursik et al., 1992):  The last term on the right-hand side of Eq. (12) represents the aggregation term which, following Smoluchowski (1917), 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 mass m+m , are 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 solid particles of the original families due to aggregation. The positive term on the righthand side of the second equation represents the increased number 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 "collision 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 kernels are implemented, from simple constant kernels to more complex formulations describing wet aggregation, where the collision frequency depends on particles properties (size, density) and the efficiency of aggregation depends on the amount of water in the mixture (see Appendix B3). Let us consider now a mass interval , which can be the whole solid particle mass spectrum or a smaller subset, as the kth section for the family i p , defined by [m i p ,k−1 , m i p ,k ]. Multiplying Eq. (12) by m i and integrating over , we obtain the following differential equations for the moments: where the last term is defined bŷ The set represents the set of couples of particles which form, when aggregating, particles with mass in the interval . Let us consider now the kth section for the family i p , defined by [m i p ,k−1 , m i p ,k ]. For i = 0 and i = 1, Eq. (16) gives the equations governing the evolution of the number of particles and mass of particles of the family i p in the kth mass interval per unit volume of the mixture, respectively: 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 terms on the right-hand side of the two previous equations (w (·) i p andŜ (·) i p ,k ) are calculated first by reconstructing the continuous number density functions in each section from the known moments (number and mass of solid particles in the kth 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 solid 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 air 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 along the plume axis (left-hand-side term) is due to ingestion of atmospheric air into the column (first righthand-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 Along 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 righthand-side term), while along the vertical axis (Eq. 25) momentum is affected by the gravitational acceleration term and the fallout of particles. Mixture temperature T mix varies when the 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 enthalpy, potential energy, and kinetic energy as The specific enthalpy of the mixture H is written as the massweighted 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 vapor, liquid water, ice, and other volcanic gases, respectively. All the mass fractions are denoted by 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 specific heat capacities at constant pressure 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 vapor and liquid water, respectively. The values of specific heat capacities and enthalpies at a reference temperature are listed in Table D1 in Appendix D. The freezing temperature is here fixed at T ref = 273.15 K, and we assumed that in the temper- vapor, liquid, and ice forms may coexist. Further details on the relationships between mixture temperature and mass fractions of water in gas, liquid, and solid states are given in Appendix B1.
Denoting the specific humidity of the atmosphere (mass of water vapor in a unit mass of moist air, expressed here as kilograms of vapor per kilogram of air) by sh atm , 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.
The transport of dry air through the column is modeled as Equation (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 column of atmospheric water as where x w is the total mass fraction of water (which can be in either vapor, liquid, and ice form) in the mixture. The partitioning of total water in the different phases is detailed in Sect. 2.3 and in Appendix B1.
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. (11) and (19-31), describing the steady dynamics of volcanic plume rise in the atmosphere, is solved numerically by PLUME-MoM-TSM. A fifth-order sevenstage Dormand-Prince Runge-Kutta method (Dormand and Prince, 1980) is implemented in a Fortran 90 code for the numerical integration of the system of ordinary differential equations. This method is based on a fifth-order method used to advance the solution, which is compared with a fourthorder method to estimate the integration error and to automatically reduce or increase the integration step. The computational time requested for a run is of the order of a few tenths of a second for a simulation without aggregation and slightly longer for simulations also modeling the aggregation process.

Phase changes of water
The present version of PLUME-MoM accounts for phase change of water as a function of the pressure and temperature conditions. Water (in vapor, liquid, or ice phases) 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 vapor, liquid, and ice: 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 vapor, liquid water, and ice phases in the mixture, respectively. In the following, we will use the subscripts da, wv, and vg to refer to dry air, water vapor, and volcanic gases, respectively. With this notation, assuming that the mixture of dry air, water vapor, 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 (·) : Molar fractions of water vapor (n wv ), volcanic gases (n vg ), and dry air (n da ) in the gas phase of the volcanic mixture (for which 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 vapor 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 fractions with respect to the mixture and the molar weight of the kth volcanic gas. For CO 2 and SO 2 (most abundant volcanic gases after water), mw k values are 0.044 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 Eqs. (27) and (32-36). The details of the procedure employed to compute these values are given in Appendix B1.

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.
Before addition of external water, specific enthalpy at the vent is where x s and x wv are the mass fractions of particles and magmatic water vapor 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 vent new ): where x lw ext is the mass fraction of external liquid water at temperature T lw ext , and x erupt = 1 − x lw ext is the magmatic mass fraction. From the value of H vent new , we update magma temperature and water partitions by following the procedure described in Appendix B1.
To test the consistency of the results produced by PLUME-MoM-TSM in the presence of external water, in the Supplement, 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 tephra is 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 intrusive 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 atmosphere thickness) so that vertical acceleration is 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, only opposed by turbulent drag, associated with the formation of shear layers between the ambient and umbrella cloud (Abraham et al., 1979). Furthermore, for the umbrella cloud model, we neglect the loss of solid 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 column into the atmosphere is provided by imposing a radial axisymmetric 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. It is important to observe that the transitional region between the rising plume and umbrella cloud, characterized by the plume overshoot and then the descent to the neutral buoyancy level, has an aspect ratio of the flow not small in either direction and is likely to be highly turbulent. For these reasons, this transition region is not properly described by either the rising plume or the shallow-water model for the umbrella cloud, and the derivation of the inflow of the umbrella cloud model directly from the output of the plume model at the neutral buoyancy level represents a simplification of the real dynamics.
With these assumptions, if we denote h as the umbrella cloud thickness, u = (u, v) the horizontal components of the velocity, and u a = (u a , v a ) the atmospheric wind at the neutral buoyancy level, the equations for the conservation of volume and horizontal momentum are written: On the left-hand sides of the momentum equations, N is the Brunt-Väisälä frequency, which represents an 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: where ρ nbl is the atmospheric density at the neutral buoyancy level. The first term on the right-hand side of both momentum equations is the drag force, which is proportional to the square of the relative velocity through the dimensionless drag coefficient C D . On the right-hand side of all the equations in system (Eq. 40),w is the volumetric flux source derived from the plume vertical velocity w nbl at the neutral buoyancy level and defined as a function of the horizontal coordinates (x, y), with the origin of the coordinate system fixed at the center of the plume at the neutral buoyancy level: The variables u nbl and v nbl are the horizontal components of plume velocity at the neutral buoyancy level, while γ x and γ y are two functions which take into account the vertical gradient of the plume radius at the neutral buoyancy level: where χ is the indicator function. 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 on a uniform grid in the (x, y) horizontal plane at the umbrellacloud height with a transient finite-volume code based on the numerical solver presented in de' Michieli Vitturi et al. (2019), where an implicit-explicit Runge-Kutta scheme (Ascher et al., 1997) is adopted to compute explicitly the flux terms on the left-hand side of the equations and the source terms, and implicitly the drag terms in the momentum equations. The solution is advanced in time until a steady upwind spreading is reached, with the integration time step controlled by a CFL condition, i.e., by a constrain on the ratio between the distance traveled in one time step by the fastest waves in the solution and the size of the cells of the computational grid (Courant et al., 1967). An example of the output of the coupled plume-umbrella models is presented in Fig. 4.

Applications
In this section, we present a suite of simulations to highlight the novelties introduced in this version of PLUME-MoM and to constrain model parameters. First, we model the spreading of the umbrella cloud for the eruption of the Calbuco volcano in southern Chile in April 2015 by comparing simulation results with rates of umbrella cloud expansion and upwind spreading derived from satellite observations. This allows us to constrain the coefficient C D in Eq. (40), controlling the spreading of the umbrella. Then, by using a modified standard atmosphere, an analysis of the effects of wind, atmospheric humidity, and mass flow rate on the upwind spreading of the umbrella cloud is presented.

April 2015 Calbuco eruption
We present here model results obtained for the two phases of the eruption of the Calbuco volcano in April 2015 (Romero et al., 2016;Castruccio et al., 2016;Van Eaton et al., 2016). Calbuco is one of the most active volcanoes of the southern Chilean Andes, and on the evening of 22 April 2015, for the first time in 43 years, Calbuco erupted, producing a sub-Plinian column with maximum plume heights exceeding 15 km above the vent (Pardini et al., 2018). Preceded by an hour-long period of volcano-tectonic events, the first phase started at 21:04 UTC and lasted approximately 1.5 h, producing tephra dispersed mostly in north and northeast directions and causing a 20 km exclusion zone to be declared. After a 5.5 h break, at 03:54 UTC on 23 April 2015, a second, larger explosive phase begun. This second phase lasted approximately 6 h, with a drop in activity of 1 h approximately 2 h after its onset (Van Eaton et al., 2016). After the second phase, discrete small explosions caused further ash emission up to 2 km height, with a major pulse which resulted in a 4 km high column on 30 April (Romero et al., 2016). Ash from the different phases of the eruption reached Brazil, Chile, Argentina, and Uruguay, causing the canceling of international flights into and out of several major cities in South America.
Following Reckziegel et al. (2019), for both phases, the total grain size distribution released at the vent (here fixed at an elevation of 2003 m) has been initialized with of a bi-modal distribution with 11 particle classes from −1φ to 9φ (see Table 1), with densities increasing linearly with φ from 2300 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 wa-  ter, 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 phases has been fixed to 250 m s −1 , 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. The atmospheric profiles used to run the simulations derive from the ECMWF-ERA5 reanalysis dataset. For the two eruptive phases, we extracted geopotential height, atmospheric density, pressure, temperature, specific humidity, and horizontal wind velocity components at the vent location and eruption starting time.
In Fig. 5, some of the plume model outputs obtained for the first phase are reported. In each panel, values below the neutral buoyancy level are plotted with a solid line, while dashed lines represent plume model results above it. Plume velocity (top-left panel) initially decelerates very rapidly in the gas thrust region and then decelerates more gently and almost linearly with height, up to the neutral buoyancy level. Mixture density, because of the large amount of air entrained from the column margin and heated in the gas thrust region, rapidly decreases below atmospheric density, making the plume positively buoyant (see the plot of the relative density in the top-right panel). The water vapor initially present in the rising mixture, with the additional contribution of the moist air entrained in the plume, stays in the gaseous state up to a height above sea level of approximately 8.5 km, above which ice starts to form. At the neutral buoyancy level, almost all the water present in the plume is constituted by ice (middle-left panel). Particles are lost during the rise, but the fraction of the initial amount lost up to the neutral buoyancy 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 that erupted reach the neutral buoyancy level, and that the initial grain size distribution differs very little from that feeding the umbrella cloud (bottomleft panel). In any case, almost immediately above the vent, due to the large entrainment, the mass of particles, as that of water, rapidly becomes 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.64 × 10 8 kg s −1 (corresponding to a volume flow rate of 2.21 × 10 9 m 3 s −1 ). The bottom-right panel presents a 3-D view of the plume, showing the weak bending of the plume axis due to both the moderate wind and the magnitude of the eruption.
The simulation performed for the second phase produced similar results, and the main output variables are plotted in Fig. 6. For this simulation, the height of the neutral buoyancy level increased to approximately 13.6 km, at which the plume radius is 3814 m, with a mass flow rate injected in the umbrella cloud of 8.38×10 8 kg s −1 (corresponding to a volume flow rate of 3.22 × 10 9 m 3 s −1 ).
For both phases, the outputs at the neutral buoyancy level of the plume model of PLUME-MoM-TSM (volumetric flow rate, radius, and horizontal velocities) were used as input of the umbrella cloud module. For this application, we varied the drag coefficient C D of the umbrella cloud model in order to find the value which best reproduced the spread-ing in the atmosphere observed during the eruption. From the experiments of Abraham et al. (1979), Baines (2013), and Johnson et al. (2015), both inferred coefficients in the range C D = 0.001 to 0.01, while Pouget et al. (2016), by comparing the results of the shallow-water intrusion model developed by Johnson et al. (2015) with satellite data from seven eruptions, better reproduced observations of umbrella cloud structure and morphological evolution with a value C D = 0.1, which was the largest value of the investigated range. Here, three values were tested (C D = 0.01, 0.1, and 1) by carrying on the simulations for 1.5 and 2 h for the first and second phases, respectively. At each time step, we computed (i) the upwind spreading, defined 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 blue lines for the equivalent radius and red lines for the upwind distances. In the right panels of Fig. 7, the edges of the umbrella cloud at the end of the simulations are plotted, with different lines for the different values of C D . From both the left and right panels, we can observe the expected larger spreading of the umbrella cloud for smaller values of C D . These results highlight the fact that intrusive gravity current dominates in the initial stages of the dispersion of tephra at the neutral buoyancy level and large upwind and crosswind spreadings of the umbrella cloud with respect to the vent location (denoted by a yellow star in the right panels of Fig. 7) are produced also for sub-Plinian eruptions. In order to constrain the value of the drag coefficient, the blue lines in the left panels of Fig. 7 are compared with the series of values for the umbrella radius reported in Van Eaton et al. (2016), obtained by detecting the edge of the umbrella from GOES-13 satellite images and plotted here with blue markers. From the results plotted in the figure, we can see that the values obtained with C D = 1 seem to better match the observations, with a small overestimate at the beginning and a small underestimate at the end of the simulation.
A different result is obtained by analyzing the modeled upwind spreading of the umbrella cloud. First of all, we can see from the red lines in the left panels of Fig. 7 that for values of C D greater than 0.1, a steady upwind distance is approached at the end of the simulation, and that by changing (increasing or decreasing) the drag coefficient by an order of magnitude, the upwind spreading at the final time changes approximately by a factor of 2. In this case, when model results are compared with the upwind spreading derived from the processing of the observations presented in Van Eaton et al. (2016), represented by the red crosses in the left panels of Fig. 7, we see that the drag coefficient value C D = 0.1 produces the best results. The two different values of C D obtained by comparing the equivalent radius and the upwind spreading with observational data can be due to the several approximations in the numerical simulations. In particular, the depth-averaged umbrella cloud model uses a constant wind (in time and space), extracted at the vent coordinates and at the neutral buoyancy level from the ECMWF-ERA5 reanalysis data. Thus, while upwind spreading, which is measured close to the vent, is slightly affected by this approximation, the spreading of the whole umbrella cloud (controlling the equivalent radius) could be affected by larger errors, because downwind wind far from the vent can be different from the assumed constant value. Furthermore, we notice that the value C D = 0.1 better agrees with the results presented in Pouget et al. (2016). For these reasons, and because in the rest of the paper we are mostly interested in quantifying the upwind spreading of the umbrella cloud, we use the value C D = 0.1 as reference value. For this value of the drag coefficient, simulated umbrella cloud thicknesses 1.5 h after the onset of each phase are compared in Fig. 8 with cloud-top IR brightness temperatures as retrieved by the NOAA GOES-13 geostationary satellite. First of all, it is important to remark that the images on the left have been cropped from a larger satellite image and have not been reprojected in the same Universal Transverse Mercator (UTM) projection zone used for the right panels (WGS 84/UTM zone 18S), so the two areas represented in the left and right panels do not correspond. In addition, as already stated, model simulations use a constant wind (in time and space), extracted at the vent coordinates and at the neutral buoyancy level, and thus downwind spreading can differ substantially from the real one. In any case, a qualitative comparison shows that the cross-wind spreading of the two phases of the eruption matches well with that predicted by the model. In addition, we observe from the contour plots in the right panels that the larger volumetric flow rate injected at the neutral buoyancy level for the second phase resulted in a thicker cloud, with a total height of column and umbrella of approximately 15.5 and 17.5 km above sea level for the first and second phases, respectively. Also, these values compare well with those reported in Van Eaton et al. (2016).
In the analysis presented so far, we have not accounted for the effect of aggregation. To better understand if aggregation 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 (φ = −2, −3) with an initial negligible number of particles have been added to the total grain size distribution. This allows for the formation of aggregates larger than the particles present in the original size spectrum. Because we are considering wet aggregation, the density of the aggregates has been fixed to 1500 kg m −3 (Costa et al., 2010), resulting in lower settling velocities than the original particles, with the size being equal.
The results of the simulation without aggregation are plotted in the top panels of Fig. 9. The grain size distribution within the plume at the neutral buoyancy level (NBL) is plotted in panel A1 of Fig. 9, the grain size distribution of the particles lost from the plume margins at the same height is plotted in panel A2, and the fraction of the initial total solid flux lost with height (which indicates the number of particles lost from the plume edges during the rise) is shown in panel A3. The grain size distribution at the neutral buoyancy level does not differ significantly from that at the vent, and this is due to the fact that the number of particles lost during the rise is very small (approximately 4 %).
As a first aggregation test, we used a modified version of the wet aggregation kernel introduced in Costa et al. (2010), defined as the product of a collision frequency function and a sticking efficiency function (see Appendix B3). According to this model, to aggregate, a collision must occur between the particles and when that collision occurs the particles must adhere to one another. Collisions between particles are associated with different mechanisms, as Brownian motion, laminar and turbulent fluid shear, and differential settling velocity. Following Costa et al. (2010), all these contributions are considered in the calculation of the collision frequency function. In our formulation, the sticking efficiency depends on the amount of liquid water and ice present in the volcanic plume. For the simulations of the Calbuco eruption, as shown in Fig. 5, only ice forms above a height of 8 km a.s.l. For this reason, aggregation is quite limited and the grain size distributions of the particles within the plume and those particles released from the plume margin at the NBL show very small differences with the simulation without aggregation, and 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 number of particles lost from the plume does not differ significantly.
In order to better understand the effect of aggregation, we performed additional simulations with a constant kernel, starting from a value of β = 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 compared to that without aggregation, with small numbers of aggregated particles produced (represented in panel B1 of Fig. 9 with the red bars). Because of the adoption of a constant kernel, aggregates form in each bin of the initial grain size distribution. By increasing β by an order of magnitude to a value of 10 −14 m 3 s −1 , the number of aggregates becomes significant (more than 50 % for some bins), and also the relative number of particles released from the plume margins at the neutral buoyancy level differs from that 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 here, most of the aggregation occurs "intra-bin"; i.e., it produces an aggregate with a volume, or mass, which belongs to the same bin of particles forming the new aggre-gate. In fact, for bin sizes of 1 φ, and a uniform mass distribution within each bin, the number of smallest particles within a bin is almost 1 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 produces a new particle with an equivalent diameter (diameter calculated from the mass) 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. 9, panels E1 and E2). It is interesting to observe that even with such a large number of aggregates, the total number 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 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-particle 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 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 behavior. With respect to the previous section, where atmospheric profiles used to run the simulations were derived from the ECMWF-ERA5 reanalysis data, the simulations presented here employ atmospheric profiles of pressure, temperature, and density 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 of the atmosphere in layers where the absolute temperature T atm against geopotential altitude h (i.e., the vertical coordinate referenced to Earth's mean sea level) 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 (i.e., the distance of a point along the plumb line to the geoid). In the following, as common in meteorology, we refer to the negative value of the rate of temperature change with altitude using the term "lapse rate", and we use for it the notation of . Once the bottom values and the lapse rate are known in an atmospheric layer, pressure, temperature, and density profiles within 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: where R is the specific gas constant of the atmosphere, which can vary with humidity. 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 the right-hand term of the first equation. ISA lapse rates (see Table 1) do not account for humidity effects and air is assumed to be dry, and the specific gas constant used in the ideal gas law is that of dry air (R = R air ). 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 vapor, 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 • (north or south) to 0.018 kg kg −1 at the Equator, and its value exponentially decreases 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. 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 constants of water vapor 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 solve the system of Eq. (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 values at the bottom of the troposphere have 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 to be a continuous function of height which varies linearly in several layers. Observations indicate that on average the 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 and 70 km (Modica et al., 2007;Brasefield, 1954). While quite a 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 speed of 20 m s −1 at 11 km are plotted in Fig. 10, where the red markers in the last two panels denote the values we varied for the sensitivity analysis.
Maximum upwind spreading of the umbrella cloud with respect to vent location represents an important information that coupled plume and umbrella models or fully threedimensional models (Cerminara et al., 2016;Suzuki and Koyaguchi, 2009) can provide. It has been noted that using only the output of integral plume models as the source of VATD models generally underestimates the upwind deposition 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 such as the 1991 Pinatubo eruption, with a mass flow rate of the order of 10 9 kg s −1 . Less work has been devoted to smaller scenarios, and for this reason we focus here 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 kg s −1 ), adjusted by changing vent diameter; (2) the wind speed at the tropopause in the range 35-80 m s −1 ; and (3) the specific humidity at 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 Latin hypercube sam- pling generating 700 elements. The other main input parameters used for this analysis are listed in Table 3. As previously stated, aggregation is disregarded. 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. As for the previous application, settling of particles from the plume edges is considered, and the settling velocity model from Textor et al. (2006b) is adopted (see Appendix B2 for more details). For each simulation, the horizontal upwind distance from the vent of the umbrella cloud has been calculated from the output 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. 11, 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 humidity does not seem 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.
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;. If we denote the output parameter of interest with Y and the input parameters considered with x i , the main sensitivity index S i , also called the Sobol index, expresses the fraction of the variability in the output Y that is reduced when the value of the input x i is fixed. It is calculated by The values of the Sobol index are shown in the bar plot presented in the right panel of Fig. 11. The largest value is obtained 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 highlights the fact the upwind spreading is mostly controlled by the first two parameters.
For this reason, we plotted in Fig. 12 the maximum upwind spreading as a function of the logarithm of the mass flow rate Figure 10. Atmospheric profiles for a dry atmosphere with a 20 m s −1 wind speed at the tropopause (11 km a.s.l.). 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. Figure 11. Upwind spreading of the umbrella cloud as a function of input parameters (a-c) and main sensitivity indices (d). 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. and the tropopause wind speed only, highlighting a potential power-law dependence of the spreading on the two 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 a = 1.22 × 10 −3 , b = 10.75, and c = 1.453 as optimal values for the fitting parameters, with a coefficient of determination of R 2 = 0.987 and a square root of the variance of the residuals of 680.8 m. 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. 12 the upwind umbrella distance as a function of the height of the NBL and the downwind horizontal distance of the centerline at the NBL from the vent, both extracted from the output of PLUME-MoM-TSM plume module. As expected, upwind spreading decreases with increasing bending of the column (expressed by the downwind distance of the centerline from the vent), 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 a power-law dependence of the upwind spreading d up on the two input variables and, if we denote with h and d down the NBL height and the downwind distance of the centerline from the vent (measured at the NBL), respectively, we can write In this case, the least-squares minimization procedure gives A = 9.055×10 −6 , B = 2.993, and C = 0.806 as optimal values for the fitting parameters, with a coefficient of determination of R 2 = 0.985 and a square root of the variance of the residuals of 719.5 m. The results of this fitting and the residuals are shown in Fig. 13, showing that the largest absolute error is obtained for the larger values of the upwind spreading. In any case, these errors of the fitting result in an overestimation of the modeled upwind spreading of less than 5 %. It is important to remark that the fitting parameters we obtained are valid for the range of tropopause wind speeds (35-80 ms −1 ) and mass flow rates (10 6 -10 8 kg s −1 ) investigated in this analysis, and that extrapolations outside these ranges should be avoided.
This second relationship, from an operational point of view, may be more useful. On one hand, the two input values are produced by all plume integral models, 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 us to have an immediate initial estimation of the upwind spreading, and thus of 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, while, for columns able to reach the tropopause, upwind spreading is significant and an umbrella model is required. Further studies will address this point by considering the uncertainty associated with additional input parameters of the models. We also remark that the values of the upwind spreading presented here are strongly dependent on the drag coefficient C D , in this analysis fixed at 0.1, and that the lower values suggested by Baines (2013) and Johnson et al. (2015) would produce a larger upwind spreading of the umbrella cloud. Additional simulations we performed (not shown here) suggest that a decrease of the drag coefficient of 1 order of magnitude results approximately in doubling the upwind spreading distance. Finally, it is worth noting that in this analysis we have not varied some of the parameters which can have a large control on the eruptive column, such as the initial water content and temperature Costa et al., 2016). 4 Conclusions PLUME-MoM-TSM, like all the integral models of volcanic plumes, is based on important simplifications, from the steadiness of the plume to the 1-D formulation. For this reason, it cannot describe the details of the complex transient and 3-D dynamics of turbulent volcanic plumes. Despite that, integral models are able to predict column heights and neutral buoyancy levels consistent with those calculated by 3-D models (Costa et al., 2016). This result, coupled with the fact that integral models require few seconds to run, pushes toward 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 particle distributions. With respect to the previous version, where a finite set of moments of the total grain size distribution was tracked, here the adoption of the two-size moment allowed us to model the evolution of the mass and the number of particles associated with each single bin of the distribution. A first important consequence is that, with such an approach, it is possible to use PLUME-MoM-TSM outputs as the 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 equations for the plume, but the procedure described is general and the code has a modular structure, which make it easy to port the procedure to other Eulerian models for transport of volcanic ash in the atmosphere. PLUME-MoM-TSM also accounts for phase changes of water, resulting in the formation of 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 the Calbuco volcano in Chile, allowing us to calibrate the drag coefficient of the umbrella cloud depth-averaged model 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. We observe that dry aggregation could produce different results, because the aggregates would have lower densities and thus lower settling velocity, strongly affecting the deposition pattern from the rising plume. In addition, we remark that most aggregates that have been mapped in proximal fallout (Self, 1983;Sisson, 1995;Wallace et al., 2013) and produced in shaker-pan experiments (Van Eaton et al., 2012) have a size distribution narrower than that produced by the aggregation kernel we adopted. This suggests that in the future the model could be updated with new collision and sticking kernels, informed by laboratory experiments and data coming from fallout deposits.
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 necessary to extend the analysis to more general initial vent conditions.

Appendix A: Reconstruction from moments
The TSM method is based on the reconstruction of the number density function from the moments, which are updated by the solution of the plume equations. We provide here the details of the TSM reconstruction of the number density functionη k (m) from the two moments (N k and M k ) in the mass interval [m k−1 , m k ].
The linear reconstructionη k (m) defined by Eq. (10), CASE 2, satisfies the conditionsη (1) k = M k for the following values of α k and β k : We observe now that the linear reconstruction is a physical approximation of a number density function only if its values are positive in the interval [m k−1 , m k ], and this can be verified by calculating the values ofη k (m) at m k−1 and m k .
If for the linear reconstruction given by CASE 2 and by the values of α k and β k defined above it holdsη k (m k−1 ) < 0, then we switch to Eq. (10), CASE 1. In this case, the conditionsη (1) k = M k are satisfied for the following values of α k and γ k : Conversely, if for the linear reconstruction given by CASE 2 and by the values of α k and β k defined above it holds η k (m k ) < 0, then we switch to Eq. (10), CASE 3. In this case, the conditionsη (0) k = N k andη (1) k = M k are satisfied for the following values of β k and γ k : We observe that we cannot have bothη k (m k−1 ) < 0 and η k (m k ) < 0 because of the positivity of the moments M k and N k .
Following Folch et al. (2016), for T ≥ 29.65 K, the saturation pressure of vapor over liquid (e l ) depends on mixture temperature in the following way: At equilibrium condition, when both liquid water and vapor are present, it must hold: where both terms can be written as functions of x wv , by appropriately combining Eqs. (33), (34), (B1), and (B3).
Before trying to solve Eq. (B3), we consider the two extreme values -x wv = 0 (liquid only) and x wv = x w (vapor only) -in order to check if they are compatible with the condition T mix ≥ T ref . If wv | x wv =x w < 0 and T mix | x wv =x w ≥ T ref , the plume is undersaturated and there is no water vapor condensation (x w = x wv ). If wv | x wv =0 > 0 and T mix | x wv =0 > T ref , all the water is in liquid form (x w = x lw ). Otherwise, we apply a bisection iterative method to solve Eq. (B3) for x wv in the interval [0, x w ]. Once the solution is obtained, we check if the corresponding temperature, as given by Eq. (B1), 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 vapor and ice form only and the temperature of the eruptive mixture can be computed as We follow the procedure adopted for T mix ≥ At equilibrium, when both vapor and ice are present, it must hold: where P wv and e s can be written as functions of x wv only, by appropriately combining Eqs. (33), (34), (B5), and (B4). Again, before solving Eq. (B6), we examine the two extreme cases -x wv = 0 (ice only) and x wv = x w (vapor only) -and we check if the corresponding temperature as evaluated from Eq. (B4) is less than T ref − 40. In the case of wv | x wv =x w < 0 and T mix | x wv =x w < T ref −40, the plume is undersaturated and water is in vapor form only (x w = x wv ). On the contrary, in the case of wv | x wv =0 > 0 and T mix | x wv =0 < T ref −40, all the water is present as ice (x w = x i ).
If equilibrium conditions can not be obtained for the two extreme cases, we solve Eq. (B6) for x w by applying a bisection procedure in the interval [0, .
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 vapor partial pressure which, depending on mixture temperature, can be evaluated from Eq. (B2) As done in Mastin (2007), we assume that liquid water mass fraction is a function of T mix as where x lw | T mix =T ref can be computed from x wv | T mix =T ref assuming that at T ref no ice is present. Finally, ice mass fraction is The specific enthalpy of the eruptive mixture at equilibrium (H eq ) is known by solving the energy equation, and, for a generic mixture temperature T , it must hold: where H | T mix =T is the mixture enthalpy at the generic temperature T . Before solving Eq. (B11), we check if equilibrium is possible in the inter- . During the procedure, we check that water partitions are positive values and their sum is equal to water mass fraction. In the event that these conditions do not hold, we skip to a different case.

B2 Settling velocity
Several models are implemented in PLUME-MoM-TSM to compute particle settling velocity (w i p ) as a function of particle diameter (D), particle density (ρ p ), and 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. (B18) is Re ψ −0.828 + 2 √ 1.07 − ψ, Re ≤ 100.
C DRe 100 + (Re−100) 1000−100 (C DRe 1000 − C DRe 100 ), 100 ≤ Re < 1000 1, Re > 1000. (B19) 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.

B3 Aggregation kernel
Following Costa et al. (2010), a model for wet aggregation based on the classical Smoluchowski equation is implemented in PLUME-MoM-TSM. In their work, a general model describing the aggregation processes was presented, but a simplified version that was computationally cheaper was then adopted. Here, the full original version is retained and implemented in the plume model, in accordance 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 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 (here constant and equal to 0.0045 s −1 , following Costa et al. (2010), Table 1). Following Costa et al. (2010), the 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. These two expressions for sticking efficiency do not take into account the amount of ice and/or liquid water available in the plume as the 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. Figure C1. (a) Effect of addition of external water on column height as a function of mass flow rate. As in Koyaguchi and Woods (1996) and Mastin (2007), input parameters are as follows: magma temperature is 1000 K, water vapor mass fraction is 0.03, and relative humidity in the atmosphere is 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). (b) 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 vapor mass fraction to 0.03, and mixture temperature to 1000 K.
column height for two eruptive scenarios. Each panel shows the results for amounts of external water of 3 (5 in Fig. C3), 10, 20, and 30 wt %. The two figures reproduce Fig. 3 of Koyaguchi and Woods (1996). Figure C2. Variation with column height of mixture temperature (a), velocity (b), relative density (c), gas mass fraction (d), liquid mass fraction (e), and ice mass fraction (f). 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, and 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 ×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 a 0.03 wt fraction of volatiles (only water vapor) and to have an initial temperature of 1000 K. Figure C3. As in Fig. C3, but the model is initialized with a vent radius of 1000 m. Figure C4. Variation of initial mixture density (ρ mix ) and temperature (T mix ) as the mass of water added to the magma increases. 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 0.01 to 0.05 wt fraction.
Following Koyaguchi and Woods (1996), Fig. C4 shows the variation of initial mixture density and temperature as the mass 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 vapor is considered). The temperature of the external water is 273 K. Mass of particles of the i p th family in the kth section kg m −3