the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
PLUMEMoMTSM 1.0.0: a volcanic column and umbrella cloud spreading model
Mattia de' Michieli Vitturi
Federica Pardini
In this paper, we present a new version of PLUMEMoM, a 1D integral volcanic plume model based on the method of moments for the description of the polydispersity in solid particles. The model describes the steadystate dynamics of a plume in a 3D coordinate system, and a modification of the twosize 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 PLUMEMoMTSM. For the first 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, PLUMEMOMTSM 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 shallowwater system of equations models the intrusive gravity current, allowing computation of the upwind spreading.
The new model is applied first 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 flow 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 volcanicash transport and dispersion models.
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 volcanicash 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: onedimensional (1D) 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 threedimensional (3D) models (Oberhuber et al., 1998; Suzuki and Koyaguchi, 2009; Cerminara et al., 2016; Ongaro et al., 2007). A recent intercomparison exercise (Costa et al., 2016) has shown that 1D and 3D 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, 1D models 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 (Herzog et al., 1998). 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 watersaturated 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 advection–diffusion–sedimentation model, and the conditions under which gravitydriven transport were analyzed. Baines (2013) studied the dynamics of this gravitational spreading into a densitystratified crossflow, deriving a semianalytic 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 shallowlayer 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. (2013, 2016), where the results of a “shallowwater” model are compared with satellite observations in order to estimate growth rates and mass fluxes. More recently, Folch et al. (2016) adopted a simple semiempirical model to describe the umbrella region, and their results are qualitatively consistent with the results obtained with more complex 3D 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 PLUMEMoM 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.
In this section, we present the equations solved by PLUMEMoMTSM. 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.
2.1 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 PLUMEMoM (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 twosize 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 secondorder 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 “twosize 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}_{\mathrm{1}}+{M}_{\mathrm{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 ${\mathit{\varphi}}_{i}={\mathit{\varphi}}_{\mathrm{0}}+i\mathrm{\Delta}\mathit{\varphi},\phantom{\rule{0.33em}{0ex}}i=\mathrm{0},\mathrm{\dots},{n}_{\mathrm{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 $\mathit{\eta}(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+\mathrm{d}m]$ at the location x in space and at time t (according to this definition, the number density function $\mathit{\eta}(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 jth 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 ${\mathit{\rho}}_{\mathrm{p}}^{B}$. 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\mathrm{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, ${\mathit{\psi}}_{k}^{\left(\mathrm{0}\right)}$ represents the average of ψ in the mass interval $[{m}_{k\mathrm{1}},{m}_{k}]$ weighted by the number of solid particles, while ${\mathit{\psi}}_{k}^{\left(\mathrm{1}\right)}$ represents a massweighted 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 npoints quadrature formula is exact for polynomials of degree 2n−1 or less. From these quadrature nodes and weights defined for the standard interval $[\mathrm{1},\mathrm{1}]$, it is possible to compute the nodes and weights for the integral over the interval $[{m}_{k\mathrm{1}},{m}_{k}]$:
where
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 reevaluate 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, nonlinear approximations $\approx {\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}_{k}\left(m\right)$ 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 ${\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}_{k}\left(m\right)$, 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 ${\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}_{k}\left(m\right)$ 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 nonaggregated particles), the subscript i_{p} is used to distinguish the different number density functions ${\mathit{\eta}}_{{i}_{\mathrm{p}}}\left(m\right),{i}_{\mathrm{p}}=\mathrm{1},\mathrm{\dots},{n}_{\mathrm{p}}$, where n_{p} is the number of particles families. Similarly, the notation ${m}_{{i}_{\mathrm{p}}}^{\left(j\right)}$ is introduced for the jth moment of the i_{p}th number density function, and the notations ${N}_{{i}_{\mathrm{p}},k}$ and ${M}_{{i}_{\mathrm{p}},k}$ for the number and mass of particles of the i_{p}th family in the kth section, respectively.
2.2 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 1D integral model where plume properties are averaged over crosssections. 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 3D 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}_{\mathrm{sc}}=\sqrt{{u}^{\mathrm{2}}+{v}^{\mathrm{2}}+{w}^{\mathrm{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}_{\mathrm{p}}=\mathrm{1},\mathrm{\dots},{n}_{\mathrm{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 crosssection 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 righthand side) and aggregation (second term on the righthand side, ${S}_{{i}_{\mathrm{p}}}$). Here, ${w}_{{i}_{\mathrm{p}}}(\cdot )$ 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 righthand side of Eq. (12) represents the aggregation term which, following Smoluchowski (1917), can be written as
where $\mathit{\beta}(z,m,{m}^{\prime})$ 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}^{\prime}$, are proportional to the number density of particles of such masses (η(m) and η(m^{′})) and to the aggregation kernel. The negative term on the righthand 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 PLUMEMoMTSM, 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}_{\mathrm{p}},k\mathrm{1}},{m}_{{i}_{\mathrm{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 by
with $\stackrel{\mathrm{\u203e}}{\mathrm{\Omega}}=\left\{(m,{m}^{\prime})>\mathrm{0}\left\right(m+{m}^{\prime})\in \mathrm{\Omega}\right\}$. The set $\stackrel{\mathrm{\u203e}}{\mathrm{\Omega}}$ 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}_{\mathrm{p}},k\mathrm{1}},{m}_{{i}_{\mathrm{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 righthand side of the two previous equations (${w}_{{i}_{\mathrm{p}}}^{(\cdot )}$ and ${\widehat{S}}_{{i}_{\mathrm{p}},k}^{(\cdot )}$) 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 righthand 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, ${\mathit{\alpha}}_{\mathit{\u03f5}}{U}_{\mathrm{sc}}{U}_{\mathrm{atm}}\mathrm{cos}\mathit{\zeta}$ is entrainment by radial inflow minus the amount swept tangentially along the plume margin by the wind, and ${\mathit{\gamma}}_{\mathit{\u03f5}}\left{U}_{\mathrm{atm}}\mathrm{sin}\mathit{\zeta}\right$ is entrainment from wind. Equation (22) states that the variation of mass flux along the plume axis (lefthandside term) is due to ingestion of atmospheric air into the column (first righthandside term) and loss of particles from the column margins (second righthandside 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
Along the horizontal directions (Eqs. 23 and 24), the momentum of the eruptive mixture varies due to the wind (first righthandside term) and the loss of particles (second righthandside 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 PLUMEMoM 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 righthand 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 temperature range $[{T}_{\mathrm{ref}}\mathrm{40},{T}_{\mathrm{ref}}]$ 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 righthand 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, PLUMEMoMTSM 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 PLUMEMoMTSM. A fifthorder 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 fifthorder 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.
2.3 Phase changes of water
The present version of PLUMEMoM 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}_{\mathrm{wv}}+{n}_{\mathrm{da}}+{n}_{\mathrm{vg}}=\mathrm{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.
2.3.1 Addition of external liquid water
The addition of external liquid water to the volcanic mixture is one of the new features of PLUMEMoMTSM. 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}_{{\mathrm{vent}}_{\mathrm{new}}}$):
where ${x}_{{\mathrm{lw}}_{\mathrm{ext}}}$ is the mass fraction of external liquid water at temperature ${T}_{{\mathrm{lw}}_{\mathrm{ext}}}$, and ${x}_{\mathrm{erupt}}=\mathrm{1}{x}_{{\mathrm{lw}}_{\mathrm{ext}}}$ is the magmatic mass fraction. From the value of ${H}_{{\mathrm{vent}}_{\mathrm{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 PLUMEMoMTSM in the presence of external water, in the Supplement, we reproduced some of the figures presented by Koyaguchi and Woods (1996).
2.4 Umbrella spreading model
As previously stated, PLUMEMoMTSM 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 “shallowwater” 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 shallowwater 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, $\mathit{u}=(u,v)$ the horizontal components of the velocity, and ${\mathit{u}}_{\mathrm{a}}=({u}_{\mathrm{a}},{v}_{\mathrm{a}})$ the atmospheric wind at the neutral buoyancy level, the equations for the conservation of volume and horizontal momentum are written:
On the lefthand 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 righthand 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 righthand side of all the equations in system (Eq. 40), $\stackrel{\mathrm{\u0303}}{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 finitevolume 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 lefthand 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.
In this section, we present a suite of simulations to highlight the novelties introduced in this version of PLUMEMoM 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.
3.1 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 subPlinian column with maximum plume heights exceeding 15 km above the vent (Pardini et al., 2018). Preceded by an hourlong period of volcanotectonic 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 bimodal 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 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 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 ECMWFERA5 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 (topleft 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 topright 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 (middleleft 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 $\mathit{\varphi}=\mathrm{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 (middleright panel). At the neutral buoyancy level, reached at 12.7 km above sea level, the plume radius is 3935 m (see bottomright 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 bottomright panel presents a 3D 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 PLUMEMoMTSM (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 spreading 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 shallowwater 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 subPlinian 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 GOES13 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 depthaveraged umbrella cloud model uses a constant wind (in time and space), extracted at the vent coordinates and at the neutral buoyancy level from the ECMWFERA5 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 cloudtop IR brightness temperatures as retrieved by the NOAA GOES13 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 crosswind 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 ($\mathit{\varphi}=\mathrm{2},\mathrm{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 $\mathit{\beta}={\mathrm{10}}^{\mathrm{15}}$ m^{3} s^{−1} up to a value $\mathit{\beta}={\mathrm{10}}^{\mathrm{13}}$ m^{3} s^{−1}. The simulation with $\mathit{\beta}={\mathrm{10}}^{\mathrm{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 bimodal shape of the total (aggregated and nonaggregated 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 “intrabin”; i.e., it produces 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 ϕ, 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 $\mathit{\beta}={\mathrm{10}}^{\mathrm{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.
3.2 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 ECMWFERA5 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 PLUMEMoMTSM, 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 righthand 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., $\dot{m}={\mathrm{10}}^{\mathrm{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 sampling 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 variancebased method (Saltelli et al., 2010; Scollo et al., 2008; de' Michieli Vitturi et al., 2016). 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 and the tropopause wind speed only, highlighting a potential powerlaw dependence of the spreading on the two input variables. In fact, if we denote with d_{up} the upwind spreading and with $\dot{m}$ and v the mass flow rate and the tropopause wind speed, respectively, we can write
A leastsquares minimization procedure gives $a=\mathrm{1.22}\times {\mathrm{10}}^{\mathrm{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 PLUMEMoMTSM 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 powerlaw 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 leastsquares minimization procedure gives $A=\mathrm{9.055}\times {\mathrm{10}}^{\mathrm{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 (de' Michieli Vitturi et al., 2016; Costa et al., 2016).
PLUMEMoMTSM, like all the integral models of volcanic plumes, is based on important simplifications, from the steadiness of the plume to the 1D formulation. For this reason, it cannot describe the details of the complex transient and 3D dynamics 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 (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 PLUMEMoMTSM, a new version of the volcanic plume model PLUMEMoM 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 twosize 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 PLUMEMoMTSM 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 twosize moment method to the steadystate 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.
PLUMEMoMTSM 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 depthaveraged 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 shakerpan 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.
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 ${\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}_{k}\left(m\right)$ from the two moments (N_{k} and M_{k}) in the mass interval $[{m}_{k\mathrm{1}},{m}_{k}]$.
The linear reconstruction ${\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}_{k}\left(m\right)$ defined by Eq. (10), CASE 2, satisfies the conditions ${\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}_{k}^{\left(\mathrm{0}\right)}={N}_{k}$ and ${\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}_{k}^{\left(\mathrm{1}\right)}={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\mathrm{1}},{m}_{k}]$, and this can be verified by calculating the values of ${\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}_{k}\left(m\right)$ 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 ${\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}_{k}\left({m}_{k\mathrm{1}}\right)<\mathrm{0}$, then we switch to Eq. (10), CASE 1. In this case, the conditions ${\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}_{k}^{\left(\mathrm{0}\right)}={N}_{k}$ and ${\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}_{k}^{\left(\mathrm{1}\right)}={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 ${\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}_{k}\left({m}_{k}\right)<\mathrm{0}$, then we switch to Eq. (10), CASE 3. In this case, the conditions ${\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}_{k}^{\left(\mathrm{0}\right)}={N}_{k}$ and ${\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}_{k}^{\left(\mathrm{1}\right)}={M}_{k}$ are satisfied for the following values of β_{k} and γ_{k}:
We observe that we cannot have both ${\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}_{k}\left({m}_{k\mathrm{1}}\right)<\mathrm{0}$ and ${\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}_{k}\left({m}_{k}\right)<\mathrm{0}$ because of the positivity of the moments M_{k} and N_{k}.
B1 Waterphase transition
For the partition of water in the gas, liquid, and solid phases, three different cases are considered: T_{mix}≥T_{ref}, ${T}_{\mathrm{ref}}\mathrm{40}<{T}_{\mathrm{mix}}<{T}_{\mathrm{ref}}$, and ${T}_{\mathrm{mix}}\le {T}_{\mathrm{ref}}\mathrm{40}$, where T_{ref} is the reference temperature of 273.15 K. Starting from the first case, we solve for the unknown variables (mass fraction of the different phases) and we check if their values are consistent with the hypothesis of thermal and mechanical equilibrium between all phases. If the equilibrium is not possible, we skip to a different case.

For T_{mix}≥T_{ref}, water can be in vapor and liquid form only and T_{mix} is given by
$$\begin{array}{}\text{(B1)}& {T}_{\mathrm{mix}}={\displaystyle \frac{H{x}_{\mathrm{lw}}\left({h}_{\mathrm{lw}\mathrm{0}}{C}_{\mathrm{lw}}{T}_{\mathrm{ref}}\right){x}_{\mathrm{wv}}\left({h}_{\mathrm{wv}\mathrm{0}}{C}_{\mathrm{wv}}{T}_{\mathrm{ref}}\right)}{{x}_{\mathrm{da}}{C}_{\mathrm{atm}}+{x}_{\mathrm{s}}{C}_{\mathrm{s}}+{x}_{\mathrm{lw}}{C}_{\mathrm{lw}}+{x}_{\mathrm{wv}}{C}_{\mathrm{wv}}+{x}_{\mathrm{vg}}{C}_{\mathrm{vg}}}}.\end{array}$$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:
$$\begin{array}{}\text{(B2)}& {e}_{\mathrm{l}}\left({T}_{\mathrm{mix}}\right)=\mathrm{611.2}\mathrm{exp}\left(\mathrm{17.67}{\displaystyle \frac{\left({T}_{\mathrm{mix}}\mathrm{273.16}\right)}{\left({T}_{\mathrm{mix}}\mathrm{29.65}\right)}}\right).\end{array}$$At equilibrium condition, when both liquid water and vapor are present, it must hold:
$$\begin{array}{}\text{(B3)}& {\mathrm{\Delta}}_{\mathrm{wv}}\equiv {P}_{\mathrm{wv}}{e}_{l}=\mathrm{0},\end{array}$$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 ${\mathrm{\Delta}}_{\mathrm{wv}}{}_{{x}_{\mathrm{wv}}={x}_{\mathrm{w}}}<\mathrm{0}$ and ${T}_{\mathrm{mix}}{}_{{x}_{\mathrm{wv}}={x}_{\mathrm{w}}}\ge {T}_{\mathrm{ref}}$, the plume is undersaturated and there is no water vapor condensation (x_{w}=x_{wv}). If ${\mathrm{\Delta}}_{\mathrm{wv}}{}_{{x}_{\mathrm{wv}}=\mathrm{0}}>\mathrm{0}$ and ${T}_{\mathrm{mix}}{}_{{x}_{\mathrm{wv}}=\mathrm{0}}>{T}_{\mathrm{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}.

For ${T}_{\mathrm{mix}}\le {T}_{\mathrm{ref}}\mathrm{40}$, water can be in vapor and ice form only and the temperature of the eruptive mixture can be computed as
$$\begin{array}{}\text{(B4)}& {T}_{\mathrm{mix}}={\displaystyle \frac{H{x}_{\mathrm{wv}}\left({h}_{\mathrm{wv}\mathrm{0}}{C}_{\mathrm{wv}}{T}_{\mathrm{ref}}\right)}{{x}_{\mathrm{da}}{C}_{\mathrm{atm}}+{x}_{\mathrm{s}}{C}_{\mathrm{s}}+{x}_{\mathrm{wv}}{C}_{\mathrm{wv}}+{x}_{\mathrm{i}}{C}_{\mathrm{i}}+{x}_{\mathrm{vg}}{C}_{\mathrm{vg}}}}.\end{array}$$We follow the procedure adopted for T_{mix}≥T_{ref} with the difference of considering the equations holding for a mixture containing water as vapor and ice. From Folch et al. (2016), the saturation pressure of vapor over ice (e_{s}) is a function of mixture temperature T_{mix} as
$$\begin{array}{}\text{(B5)}& \begin{array}{rl}& {e}_{\mathrm{s}}\left({T}_{\mathrm{mix}}\right)\\ & =\mathrm{611.22}\times {\mathrm{10}}^{\mathrm{9.097}\left({\scriptscriptstyle \frac{\mathrm{273.16}}{{T}_{\mathrm{mix}}}}\mathrm{1}\right)\mathrm{3.566}{\mathrm{log}}_{\mathrm{10}}\left({\scriptscriptstyle \frac{\mathrm{273.16}}{{T}_{\mathrm{mix}}}}\right)+\mathrm{0.876}\left(\mathrm{1}{\scriptscriptstyle \frac{{T}_{\mathrm{mix}}}{\mathrm{273.16}}}\right)}.\end{array}\end{array}$$At equilibrium, when both vapor and ice are present, it must hold:
$$\begin{array}{}\text{(B6)}& {\mathrm{\Delta}}_{\mathrm{wv}}\equiv {P}_{\mathrm{wv}}{e}_{\mathrm{s}}=\mathrm{0},\end{array}$$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 ${\mathrm{\Delta}}_{\mathrm{wv}}{}_{{x}_{\mathrm{wv}}={x}_{\mathrm{w}}}<\mathrm{0}$ and ${T}_{\mathrm{mix}}{}_{{x}_{\mathrm{wv}}={x}_{\mathrm{w}}}<{T}_{\mathrm{ref}}\mathrm{40}$, the plume is undersaturated and water is in vapor form only (x_{w}=x_{wv}). On the contrary, in the case of ${\mathrm{\Delta}}_{\mathrm{wv}}{}_{{x}_{\mathrm{wv}}=\mathrm{0}}>\mathrm{0}$ and ${T}_{\mathrm{mix}}{}_{{x}_{\mathrm{wv}}=\mathrm{0}}<{T}_{\mathrm{ref}}\mathrm{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,x_{w}] and we check if the corresponding temperature given by Eq. (B4) is $<{T}_{\mathrm{ref}}\mathrm{40}$.

For ${T}_{\mathrm{ref}}\mathrm{40}<{T}_{\mathrm{mix}}<{T}_{\mathrm{ref}}$, water can be in vapor, liquid, and ice form. Following Mastin (2007) and references therein, we assume that in this range of temperatures liquid and ice coexist but not in equilibrium, with liquid water present in a subcooled state to temperatures well below freezing. In this case, specific mixture enthalpy can be computed from Eq. (27) as a function of mixture temperature only. Indeed, we show in the following that water mass fractions can be written as functions of T_{mix}. Water vapor mass fraction is computed by inverting Eq. (34) as
$$\begin{array}{}\text{(B7)}& {x}_{\mathrm{wv}}\left({T}_{\mathrm{mix}}\right)={\displaystyle \frac{\left(\frac{{x}_{\mathrm{da}}}{m{w}_{\mathrm{da}}}+\frac{{x}_{\mathrm{vg}}}{m{w}_{\mathrm{vg}}}\right){n}_{\mathrm{wv}}m{w}_{\mathrm{wv}}}{(m{w}_{\mathrm{wv}}\mathrm{1})}}.\end{array}$$The righthandside 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
$$\begin{array}{}\text{(B8)}& {n}_{\mathrm{wv}}={\displaystyle \frac{e{}_{T={T}_{\mathrm{mix}}}}{{P}_{\mathrm{atm}}}},\end{array}$$where e is vapor partial pressure which, depending on mixture temperature, can be evaluated from Eq. (B2) for T_{mix}=T_{ref} or from Eq. (B5) for ${T}_{\mathrm{mix}}={T}_{\mathrm{ref}}\mathrm{40}$.
As done in Mastin (2007), we assume that liquid water mass fraction is a function of T_{mix} as
$$\begin{array}{}\text{(B9)}& \phantom{\rule{0.33em}{0ex}}{x}_{\mathrm{lw}}\left({T}_{\mathrm{mix}}\right)={x}_{\mathrm{lw}}{}_{{T}_{\mathrm{mix}}={T}_{\mathrm{ref}}}\left[{\displaystyle \frac{{T}_{\mathrm{mix}}({T}_{\mathrm{ref}}\mathrm{40})}{\mathrm{40}}}\right],\end{array}$$where ${x}_{\mathrm{lw}}{}_{{T}_{\mathrm{mix}}={T}_{\mathrm{ref}}}$ can be computed from ${x}_{\mathrm{wv}}{}_{{T}_{\mathrm{mix}}={T}_{\mathrm{ref}}}$ assuming that at T_{ref} no ice is present. Finally, ice mass fraction is
$$\begin{array}{}\text{(B10)}& \phantom{\rule{0.33em}{0ex}}{x}_{\mathrm{i}}\left({T}_{\mathrm{mix}}\right)={x}_{\mathrm{w}}{x}_{\mathrm{wv}}{x}_{\mathrm{lw}}.\end{array}$$The specific enthalpy of the eruptive mixture at equilibrium (H_{eq}) is known by solving the energy equation, and, for a generic mixture temperature $\stackrel{\mathrm{\u203e}}{T}$, it must hold:
$$\begin{array}{}\text{(B11)}& {\mathrm{\Delta}}_{H}\equiv {H}_{\mathrm{eq}}H{}_{{T}_{\mathrm{mix}}=\stackrel{\mathrm{\u203e}}{T}}=\mathrm{0},\end{array}$$where $H{}_{{T}_{\mathrm{mix}}=\stackrel{\mathrm{\u203e}}{T}}$ is the mixture enthalpy at the generic temperature $\stackrel{\mathrm{\u203e}}{T}$. Before solving Eq. (B11), we check if equilibrium is possible in the interval $[{T}_{\mathrm{ref}}\mathrm{40},{T}_{\mathrm{ref}}]$ by computing ${\mathrm{\Delta}}_{H}{}_{{T}_{\mathrm{mix}}={T}_{\mathrm{ref}}\mathrm{40}}$ and ${\mathrm{\Delta}}_{H}{}_{{T}_{\mathrm{mix}}={T}_{\mathrm{ref}}}$. If ${\mathrm{\Delta}}_{H}{}_{{T}_{\mathrm{mix}}={T}_{\mathrm{ref}}\mathrm{40}}>\mathrm{0}$ and ${\mathrm{\Delta}}_{H}{}_{{T}_{\mathrm{mix}}={T}_{\mathrm{ref}}}>\mathrm{0}$ or ${\mathrm{\Delta}}_{H}{}_{{T}_{\mathrm{mix}}={T}_{\mathrm{ref}}\mathrm{40}}<\mathrm{0}$ and ${\mathrm{\Delta}}_{H}{}_{{T}_{\mathrm{mix}}={T}_{\mathrm{ref}}}<\mathrm{0}$, it means that equilibrium cannot be reached in the range $[{T}_{\mathrm{ref}}\mathrm{40},{T}_{\mathrm{ref}}]$ and we examine the case ${T}_{\mathrm{mix}}<{T}_{\mathrm{ref}}\mathrm{40}$. Otherwise, we apply a bisection procedure to solve for T_{mix} and water partitions in the interval $[{T}_{\mathrm{ref}}\mathrm{40},{T}_{\mathrm{ref}}]$. 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 PLUMEMoMTSM to compute particle settling velocity (${w}_{{i}_{\mathrm{p}}}$) as a function of particle diameter (D), particle density (ρ_{p}), and shape factor (ψ).
The first model is taken from Textor et al. (2006b) and references therein:
where ${k}_{\mathrm{1}}=\mathrm{1.19}\times {\mathrm{10}}^{\mathrm{5}}$ m^{2} kg^{−1} s^{−1}, k_{2}=8 m^{3} kg^{−1} s^{−1}, k_{3}= 4.833 m^{2} kg${}^{\mathrm{1}/\mathrm{2}}$ s^{−1}, and C_{D} is the drag coefficient that we set equal to 0.75.
The second model (Ganser, 1993) defines the settling velocity for the Stokes regime (Re≪0.005) as
The drag coefficient (C_{D}) in Eq. (B13) 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 shape factor:
Finally, settling velocity can be computed following Pfeiffer et al. (2005) as
where m is the particle mass and A_{cs} the crosssectional area of the particle. Depending on the Reynolds number, the drag coefficient used in Eq. (B18) is
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 PLUMEMoMTSM 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 PLUMEMoMTSM. 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.
To test the consistency of our results with those presented by Koyaguchi and Woods (1996) and later by Mastin (2007) for the addition of external water at the base of the plume, we applied PLUMEMoMTSM to reproduce some of the figures presented in their works. The left panel of Fig. C1 shows the effects of external water on plume height as a function of mass flow rate. PLUMEMoMTSM results (orange lines) are compared with those of Plumeria (light blue lines digitized from Fig. C4 of Mastin (2007)) and those of Koyaguchi and Woods (1996) (green lines digitized from Fig. 5). A general agreement can be seen between our results and those of the other models, with a maximum difference of 20 % between the maximum column height found by PLUMEMoMTSM and the one found by Plumeria for the 10 % external water case. For the 0 % case, PLUMEMoMTSM predicts column collapse for a slightly lower mass eruption rate than the other models, while a good fit between PLUMEMoMTSM and Koyaguchi and Woods (1996) results is observed for the other cases.
The right panel of Fig. C1 reproduces Fig. 4 of Koyaguchi and Woods (1996), where the height of the eruptive column as a function of the mass of external water added at the vent is shown. Each curve is drawn for a constant mass flow rate of 2.24 ×10^{5}, 5.6 ×10^{6}, 2.24 ×10^{7}, 5.6 ×10^{8}, and 2.24 ×10^{9} kg s^{−1}. These values correspond to vent radii of 10, 50, 100, 500, and 1000 m, respectively. The PLUMEMoMTSM results (orange lines) appear to be similar to those of Koyaguchi and Woods (1996) (light blue lines), although PLUMEMoMTSM predicts column collapse for slightly lower values of external water than Koyaguchi and Woods (1996). Figures C2 and C3 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 (5 in Fig. C3), 10, 20, and 30 wt %. The two figures reproduce Fig. 3 of Koyaguchi and Woods (1996).
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.
The numerical code, benchmark tests, and documentation are available at https://github.com/demichie/PLUMEMoMTSM (last access: 23 June 2020). Postprocessing scripts (to plot the solution variables and to create animations) are also available. Furthermore, each example has a page description on the model wiki (https://github.com/demichie/PLUMEMoMTSM/wiki (last access: 23 June 2020), where detailed information on how to run the simulations is given. The digital object identifier (DOI) for the version of the code documented in this paper is https://doi.org/10.5281/zenodo.3904379 (de' Michieli Vitturi and Pardini, 2020).
No data sets were used in this article.
Both authors equally contributed to the model development, the numerical simulations, the analysis of the results, and the preparation of the manuscript.
The authors declare that they have no conflict of interest.
This work has been supported by the Italian Department of Civil Protection, INGVDPC agreement A 2020. Federica Pardini's research activities were partially funded by the Italian MIUR project Premiale AshRESILIENCE, the MIUR FISR project “Sale Operative Integrate e Reti di Monitoraggio del Futuro”, and the European project EUROVOLC (grant no. 731070). The authors also thank Alessandro Tadini, Kyle Mohr, and Silvia Giansante for testing the code and for the feedback provided, and Giovanni Biagioli for the help in the definition of the modified standard atmosphere. Moreover, the authors thank the reviewers (Larry Mastin, Chris Johnson, and Julien Savre) for their thoughtful comments and suggestions which helped to improve the paper.
This research has been supported by the Italian Department of Civil Protection (grant no. INGVDPC agreement A 2020), the Italian MIUR (grant no. project Premiale AshRESILIENCE), the MIUR (FISR project “Sale Operative Integrate e Reti di Monitoraggio del Futuro”), and the EU (EUROVOLC (grant no. 731070)).
This paper was edited by Thomas Poulet and reviewed by Larry Mastin, Chris Johnson, and Julien Savre.
Abraham, G., Karelse, M., and Van Os, A.: On the magnitude of interfacial shear of subcritical stratified flows in relation with interfacial stability, J. Hydraul. Res., 17, 273–287, 1979. a, b
Abramowitz, M. and Stegun, I. A.: Handbook of mathematical functions with formulas, graphs, and mathematical table, US Department of Commerce, National Bureau of Standards Applied Mathematics , Washington DC, USA, 1965. a
Ahrens, C. D.: Meteorology today: an introduction to weather, climate, and the environment, Cengage Learning, Boston, MA, 2012. a
Anderson, G. P., Clough, S. A., Kneizys, F., Chetwynd, J. H., and Shettle, E. P.: AFGL atmospheric constituent profiles (0.120 km), Technical Report, Air Force Geophysics Lab, Hanscom, MA, USA, 46 pp., 1986. a
Ascher, U. M., Ruuth, S. J., and Spiteri, R. J.: Implicitexplicit RungeKutta methods for timedependent partial differential equations, Appl. Numer. Math., 25, 151–167, 1997. a
Baines, P. G.: The dynamics of intrusions into a densitystratified crossflow, Phys. Fluids, 25, 076601, https://doi.org/10.1063/1.4811850, 2013. a, b, c, d
Bonadonna, C., Biass, S., Menoni, S., and Gregg, C. E.: Assessment of risk associated with tephrarelated hazards, in: Forecasting and Planning for Volcanic Hazards, Risks, and Disasters, Volume 2, edited by: Papale, P., Elsevier, Amsterdam, Netherlands, 329–378, https://doi.org/10.1016/B9780128180822.000081, 2021. a
Brasefield, C.: Winds at altitudes up to 80 kilometers, J. Geophys. Res., 59, 233–237, 1954. a
Brown, R., Bonadonna, C., and Durant, A.: A review of volcanic ash aggregation, Phys. Chem. Earth, 45, 65–78, 2012. a, b
Bursik, M.: Effect of wind on the rise height of volcanic plumes, Geophys. Res. Lett., 18, 3621–3624, 2001. a, b, c
Bursik, M., Sparks, R., Gilbert, J., and Carey, S.: Sedimentation of tephra by volcanic plumes: I. Theory and its comparison with a study of the Fogo A plinian deposit, Sao Miguel (Azores), B. Volcanol., 54, 329–344, 1992. a
Carey, S. and Sparks, R.: Quantitative models of the fallout and dispersal of tephra from volcanic eruption columns, B. Volcanol., 48, 109–125, 1986. a, b, c
Castruccio, A., Clavero, J., Segura, A., Samaniego, P., Roche, O., Le Pennec, J.L., and Droguett, B.: Eruptive parameters and dynamics of the April 2015 subPlinian eruptions of Calbuco volcano (southern Chile), B. Volcanol., 78, 1–19, https://doi.org/10.1007/s0044501610588, 2016. a
Cerminara, M., Esposti Ongaro, T., and Berselli, L. C.: ASHEE1.0: a compressible, equilibrium–Eulerian model for volcanic ash plumes, Geosci. Model Dev., 9, 697–730, https://doi.org/10.5194/gmd96972016, 2016. a, b
Colucci, S., Landi, P., and de' Michieli Vitturi, M.: CrystalMom: a new model for the evolution of crystal size distributions in magmas with the quadraturebased method of moments, Contrib. Mineral. Petr., 172, 100, https://doi.org/10.1007/s0041001714216, 2017. a
Costa, A., Folch, A., and Macedonio, G.: A model for wet aggregation of ash particles in volcanic plumes and clouds: 1. Theoretical formulation, J. Geophys. Res.Sol. Ea., 115, 1978–2012, 2010. a, b, c, d, e, f, g, h, i
Costa, A., Folch, A., and Macedonio, G.: Densitydriven transport in the umbrella region of volcanic clouds: Implications for tephra dispersion models, Geophys. Res. Lett., 40, 4823–4827, https://doi.org/10.1002/grl.50942, 2013. a, b, c
Costa, A., Suzuki, Y. J., Cerminara, M., Devenish, B. J., Esposti Ongaro, T., Herzog, M., Van Eaton, A. R., Denby, L. C., Bursik, M., de' Michieli Vitturi, M., Engwell, S., Neri, A., Barsotti, S., Folch, A., Macedonio, G., Girault, F., Carazzo, G., Tait, S., Kaminski, E., Mastin, L. G., Woodhouse, M. J., Phillips, J. C., Hogg, A. J., Degruyter, W., and Bonadonna, C.: Results of the eruptive column model intercomparison study, J. Volcanol. Geoth. Res., 326, 2–25, https://doi.org/10.1016/j.jvolgeores.2016.01.017, 2016. a, b, c, d
Courant, R., Friedrichs, K., and Lewy, H.: On the partial difference equations of mathematical physics, IBM J. Res. Dev., 11, 215–234, 1967. a
Daidzic, N. E.: On atmospheric lapse rates, Int. J. Aviat. Aeronaut. Aerosp., 6, https://doi.org/10.15394/ijaaa.2019.1374, 2019. a, b
Degruyter, W. and Bonadonna, C.: Improving on mass flow rate estimates of volcanic eruptions, Geophys. Res. Lett., 39, L16308, https://doi.org/10.1029/2012GL052566, 2012. a
de' Michieli Vitturi, M. and Pardini, F.: PLUMEMoMTSM v. 1.0.0 (Version 1.0.0), Zenodo, https://doi.org/10.5281/zenodo.3904379, 2020. a
de' Michieli Vitturi, M., Neri, A., and Barsotti, S.: PLUMEMoM 1.0: A new integral model of volcanic plumes based on the method of moments, Geosci. Model Dev., 8, 2447–2463, https://doi.org/10.5194/gmd824472015, 2015. a, b, c, d
de' Michieli Vitturi, M., Engwell, S., Neri, A., and Barsotti, S.: Uncertainty quantification and sensitivity analysis of volcanic columns models: Results from the integral model PLUMEMoM, J. Volcanol. Geoth. Res., 326, 77–91, https://doi.org/10.1016/j.jvolgeores.2016.03.014, 2016. a, b
de' Michieli Vitturi, M., Esposti Ongaro, T., Lari, G., and Aravena, A.: IMEX_SfloW2D 1.0: a depthaveraged numerical flow model for pyroclastic avalanches, Geosci. Model Dev., 12, 581–595, https://doi.org/10.5194/gmd125812019, 2019. a
Devenish, B.: Using simple plume models to refine the source mass flux of volcanic eruptions according to atmospheric conditions, J. Volcanol. Geoth. Res., 256, 118–127, 2013. a
Dormand, J. R. and Prince, P. J.: A family of embedded RungeKutta formulae, J. Comput. Appl. Math., 6, 19–26, 1980. a
Dutton, J. A.: The ceaseless wind: An introduction to the theory of atmospheric motion, Courier Corporation, McGrawHill, New York, 2002. a
Folch, A., Costa, A., and Macedonio, G.: FPLUME1.0: An integral volcanic plume model accounting for ash aggregation, Geosci. Model Dev., 9, 431–450, https://doi.org/10.5194/gmd94312016, 2016. a, b, c, d, e, f
Ganser, G. H.: A rational approach to drag prediction of spherical and nonspherical particles, Powder Technol., 77, 143–152, 1993. a
Graf, H.F., Herzog, M., Oberhuber, J. M., and Textor, C.: Effect of environmental conditions on volcanic plume rise, J. Geophys. Res.Atmos., 104, 24309–24320, 1999. a
Herzog, M., Graf, H.F., Textor, C., and Oberhuber, J. M.: The effect of phase changes of water on the development of volcanic plumes, J. Volcanol. Geoth. Res., 87, 55–74, 1998. a
Hewett, T., Fay, J., and Hoult, D.: Laboratory experiments of smokestack plumes in a stable atmosphere, Atmos. Environ., 5, 767–789, 1971. a
Johnson, C. G., Hogg, A. J., Huppert, H. E., Sparks, R. S. J., Phillips, J. C., Slim, A. C., and Woodhouse, M. J.: Modelling intrusions through quiescent and moving ambients, J. Fluid Mech., 771, 370–406, 2015. a, b, c, d, e, f, g
Koyaguchi, T. and Woods, A. W.: On the formation of eruption columns following explosive mixing of magma and surfacewater, J. Geophys. Res.Sol. Ea., 101, 5561–5574, 1996. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o
Künzli, P., Falcone, J.L., Rossi, E., Albuquerque, P., and Chopard, B.: HPC multiscale simulation of transport and aggregation of volcanic particles, in: Proceedings of the 17th International Symposium on Parallel and Distributed Computing (ISPDC), Proceedings of 2018 17th International Symposium on Parallel and Distributed Computing (ISPDC), 25–28 June 2018, Geneva, Switzerland, 25–32, 2018. a
Liu, L. and Litster, J.: Population balance modelling of granulation with a physically based coalescence kernel, Chem. Eng. Sci., 57, 2183–2191, 2002. a
Marchisio, D. L. and Fox, R. O.: Computational models for polydisperse particulate and multiphase systems, Cambridge University Press, Cambridge, UK, 2013. a
Marchisio, D. L., Vigil, R. D., and Fox, R. O.: Quadrature method of moments for aggregation–breakage processes, J. Colloid Interf. Sci., 258, 322–334, 2003. a
Mastin, L. G.: A userfriendly onedimensional model for wet volcanic plumes, Geochem. Geophy. Geosy., 8, Q03014, https://doi.org/10.1029/2006GC001455, 2007. a, b, c, d, e, f, g, h, i
Modica, G. D., Nehrkorn, T., and Myers, T.: An Investigation of Stratospheric Winds in Support of the High Altitude Airship, available at: https://ams.confex.com/ams/pdfpapers/128256.pdf (last access: 23 February 2021), 2007. a
Morton, B., Taylor, G., and Turner, J.: Turbulent gravitational convection from maintained and instantaneous sources, in: P. Roy. Soc. AMath. Phy., 234, 1–23, 956. a, b
Nguyen, T. T., Laurent, F., Fox, R. O., and Massot, M.: Solution of population balance equations in applications with fine particles: mathematical modeling and numerical schemes, J. Comput. Phys., 325, 129–156, 2016. a, b, c, d
Oberhuber, J. M., Herzog, M., Graf, H.F., and Schwanke, K.: Volcanic plume simulation on large scales, J. Volcanol. Geoth. Res., 87, 29–53, 1998. a
Ongaro, T. E., Cavazzoni, C., Erbacci, G., Neri, A., and Salvetti, M.V.: A parallel multiphase flow code for the 3D simulation of explosive volcanic eruptions, Parallel Comput., 33, 541–560, 2007. a
Pardini, F., Burton, M., Arzilli, F., La Spina, G., and Polacci, M.: SO_{2} emissions, plume heights and magmatic processes inferred from satellite data: The 2015 Calbuco eruptions, J. Volcanol. Geoth. Res., 361, 12–24, 2018. a
Pfeiffer, T., Costa, A., and Macedonio, G.: A model for the numerical simulation of tephra fall deposits, J. Volcanol. Geoth. Res., 140, 273–294, 2005. a
Pouget, S., Bursik, M., Sparks, R., Hogg, A., Johnson, C., Singh, T., and Pavolonis, M.: Gravity current model of the volumetric growth of volcanic clouds: remote assessment with satellite imagery and estimation of mass eruption rate, American Geophysical Union, Fall Meeting 2013, San Francisco (USA), December 2013, abstract id. V23C2861, 2013. a
Pouget, S., Bursik, M., Johnson, C. G., Hogg, A. J., Phillips, J. C., and Sparks, R. S. J.: Interpretation of umbrella cloud growth and morphology: implications for flow regimes of shortlived and longlived eruptions, B. Volcanol., 78, 1, https://doi.org/10.1007/s0044501509930, 2016. a, b, c
Reckziegel, F., Folch, A., and Viramonte, J.: ATLAS1.0: Atmospheric Lagrangian dispersion model for tephra transport and deposition, Comput. Geosci., 131, 41–51, 2019. a, b
Romero, J. E., Morgavi, D., Arzilli, F., Daga, R., Caselli, A., Reckziegel, F., Viramonte, J., DíazAlvarado, J., Polacci, M., Burton, M., and Perugini, D.: Eruption dynamics of the 22–23 April 2015 Calbuco Volcano (Southern Chile): Analyses of tephra fall deposits, J. Volcanol. Geoth. Res., 317, 15–29, 2016. a, b
Rossi, E., Pollastri, S., and Bonadonna, C.: Ash aggregation in volcanic plumes: when, where, how. A new theoretical perspective, EGUGA, 12854, 20th EGU General Assembly, EGU2018, Proceedings from the conference held 4–13 April, 2018 in Vienna, Austria, p. 12854, 2018. a
Saltelli, A., Annoni, P., Azzini, I., Campolongo, F., Ratto, M., and Tarantola, S.: Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index, Comput. Phys. Commun., 181, 259–270, https://doi.org/10.1016/j.cpc.2009.09.018, 2010. a
Scollo, S., Tarantola, S., Bonadonna, C., Coltelli, M., and Saltelli, A.: Sensitivity analysis and uncertainty estimation for tephra dispersal models, J. Geophys. Res.Sol. Ea., 113, 1–17, https://doi.org/10.1029/2006JB004864, 2008. a
Self, S.: Largescale phreatomagmatic silicic volcanism: a case study from New Zealand, J. Volcanol. Geoth. Res., 17, 433–469, 1983. a
Sisson, T.: Blast ashfall deposit of May 18, 1980 at Mount St. Helens, Washington, J. Volcanol. Geoth. Res., 66, 203–216, 1995. a
Smoluchowski, M.: Attempt at a mathematical theory of coagulation kinetics of colloidal solutions, Z. Phys. Chem., 92, 129–168, 1917. a
Sparks, R.: The dimensions and dynamics of volcanic eruption columns, B. Volcanol., 48, 3–15, 1986. a, b
Sparks, R. S. J., Bursik, M., Carey, S., Gilbert, J., Glaze, L., Sigurdsson, H., and Woods, A.: Volcanic plumes, Wiley, United States, 1997. a
Suzuki, Y. and Koyaguchi, T.: A threedimensional numerical simulation of spreading umbrella clouds, J. Geophys. Res.Sol. Ea., 114, B03209, https://doi.org/10.1029/2007JB005369, 2009. a, b
Textor, C., Graf, H., Herzog, M., Oberhuber, J., Rose, W. I., and Ernst, G.: Volcanic particle aggregation in explosive eruption columns. Part II: Numerical experiments, J. Volcanol. Geoth. Res., 150, 378–394, 2006a. a
Textor, C., Graf, H.F., Herzog, M., Oberhuber, J. M., Rose, W. I., and Ernst, G. G.: Volcanic particle aggregation in explosive eruption columns. Part I: Parameterization of the microphysics of hydrometeors and ash, J. Volcanol. Geoth. Res., 150, 359–377, 2006b. a, b, c, d
Þórarinsson, S.: Surtsey: the new island in the North Atlantic, Viking Press, New York, USA, 1967. a
Ungarish, M. and Huppert, H. E.: On gravity currents propagating at the base of a stratified ambient, J. Fluid Mech., 458, 283–301, 2002. a
Van Eaton, A. R., Muirhead, J. D., Wilson, C. J., and Cimarelli, C.: Growth of volcanic ash aggregates in the presence of liquid water and ice: an experimental approach, B. Volcanol., 74, 1963–1984, 2012. a
Van Eaton, A. R., Amigo, Á., Bertin, D., Mastin, L. G., Giacosa, R. E., González, J., Valderrama, O., Fontijn, K., and Behnke, S. A.: Volcanic lightning and plume behavior reveal evolving hazards during the April 2015 eruption of Calbuco volcano, Chile, Geophys. Res. Lett., 43, 3563–3571, 2016. a, b, c, d, e, f
Veitch, G. and Woods, A. W.: Particle aggregation in volcanic eruption columns, J. Geophys. Res.Sol. Ea., 106, 26425–26441, 2001. a
Wallace, K. L., Schaefer, J. R., and Coombs, M. L.: Character, mass, distribution, and origin of tephrafall deposits from the 2009 eruption of Redoubt Volcano, Alaska – Highlighting the significance of particle aggregation, J. Volcanol. Geoth. Res., 259, 145–169, 2013. a
Webster, H. N., Devenish, B. J., Mastin, L. G., Thomson, D. J., and Van Eaton, A. R.: Operational modelling of umbrella cloud growth in a lagrangian volcanic ash transport and dispersion model, Atmosphere, 11, 200, https://doi.org/10.3390/atmos11020200, 2020. a, b
Woodhouse, M., Hogg, A., Phillips, J., and Sparks, R.: Interaction between volcanic plumes and wind during the 2010 Eyjafjallajökull eruption, Iceland, J. Geophys. Res.Sol. Ea., 118, 92–109, 2013. a, b
Woods, A. W.: The fluid dynamics and thermodynamics of eruption columns, B. Volcanol., 50, 169–193, 1988. a
Woods, A. W.: Moist convection and the injection of volcanic ash into the atmosphere, J. Geophys. Res.Sol. Ea., 98, 17627–17636, 1993. a, b
 Abstract
 Introduction
 Model equations
 Applications
 Conclusions
 Appendix A: Reconstruction from moments
 Appendix B: Constitutive equations
 Appendix C: Addition of external water
 Appendix D: List of model variables
 Code availability
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Model equations
 Applications
 Conclusions
 Appendix A: Reconstruction from moments
 Appendix B: Constitutive equations
 Appendix C: Addition of external water
 Appendix D: List of model variables
 Code availability
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References