FPLUME-1.0: An integral volcanic plume model accounting for ash aggregation

Abstract. Eruption Source Parameters (ESP) characterizing volcanic eruption plumes are crucial inputs for atmospheric tephra dispersal models, used for hazard assessment and risk mitigation. We present FPLUME-1.0, a steady-state 1D cross-section averaged eruption column model based on the Buoyant Plume Theory (BPT). The model accounts for plume bending by wind, entrainment of ambient moisture, effects of water phase changes, particle fallout and re-entrainment, a new param5 eterization for the air entrainment coefficients and a model for wet aggregation of ash particles in presence of liquid water or ice. In the occurrence of wet aggregation, the model predicts an “effective” grain size distribution depleted in fines with respect to that erupted at the vent. Given a wind profile, the model can be used to determine the column height from the eruption mass flow rate or vice-versa. The ultimate goal is to improve ash cloud dispersal forecasts by better constraining the 10 ESP (column height, eruption rate and vertical distribution of mass) and the “effective” particle grain size distribution resulting from eventual wet aggregation within the plume. As test cases we apply the model to the eruptive phase-B of the 4 April 1982 El Chichón volcano eruption (México) and the 6 May 2010 Eyjafjallajökull eruption phase (Iceland). The modular structure of the code facilitates the implementation in the future code versions of more quantitative ash aggregation parameterization as 15 further observations and experiments data will be available for better constraining ash aggregation processes.

sized blocks to micron-sized fine ash (diameter ≤ 63µm). Sustained volcanic plumes present a negatively buoyant basal thrust region where the mixture rises due to its momentum. As ambient air is entrained by turbulent mixing, it heats and expands, thereby reducing the average density of the mixture. It leads to a transition to the convective region, in which positive buoyancy drives the mix-25 ture up to the so-called Neutral Buoyancy Level (NBL), where the mixture density equals that of the surrounding atmosphere. Excess of momentum above the NBL (overshooting) can drive the mixture higher forming the umbrella region, where tephra disperses horizontally first as a gravity current (e.g. Costa et al., 2013;Carazzo et al., 2014) and then under passive wind advection forming a volcanic cloud (see Fig. 1). 30 Quantitative observations and models of volcanic plumes are essential to provide realistic source terms to atmospheric dispersal models, aimed at simulating atmospheric tephra transport and/or the resulting fallout deposit (e.g. Folch, 2012). Plume models range in complexity from 1D integral models built upon the Buoyant Plume Theory (BPT) of Morton et al. (1956) to sophisticated multiphase Computational Fluid Dynamics (CFD) models (e.g. Suzuki et al., 2005;Esposti Ongaro 35 et al., 2007;Suzuki and Koyaguchi, 2009;Herzog and Graf, 2010;Suzuki and Koyaguchi, 2013).
The latter group of models are valuable to understand physical phenomena and the role of different parameters but, given their high computational cost, coupling with atmospheric dispersal models at an operational level is still unpractical. Moreover even sophisticated 3D multiphase models can have serious problems to accurately describe the physical processes related to e.g. closure equations, 40 computational spatial resolution, etc. For this reason, simpler 1D cross-section averaged models or even empirical relationships between plume height and eruption rate (e.g. Mastin et al., 2009;Degruyter and Bonadonna, 2012) are used in practice to furnish Eruption Source Parameters (ESP) to atmospheric transport models, the results of which strongly depend on the source term quantification (i.e. determination of plume height, eruption rate, vertical distribution of mass and particle grain size 45 distribution).
Many plume models based on the BPT have been proposed after the seminal studies of Wilson (1976) and Sparks (1986) to address different aspects of plume dynamics. For example, Woods (1993) proposed a model to include the latent heat associated with condensation of water vapour and quantify its effects upon the eruption column. Ernst et al. (1996) presented a model considering 50 particle sedimentation and re-entrainment from plume margins. Bursik (2001) analyzed how the interaction with wind enhances entrainment of air, plume bending, and decrease of the total plume height for a given eruption rate. Several other plume models exist (see Costa et al., 2015, and references therein). considering different modelling approaches, simplifying assumptions and model parameterizations. It is well recognized that the values of the air entrainment coefficients have a large 55 influence on the results of the plume models. On the other hand, volcanic ash aggregation (e.g. Brown et al., 2012) can occur within the eruption column or, under certain circumstances, downstream within the ash cloud . In any case, the formation of ash aggregates (with typical sizes around few hundreds of µm and less dense than the primary particles) dramatically impacts particle transport dynamics thereby reducing the atmospheric residence time of aggregating particles 60 and promoting the premature fallout of fine ash. As a result, atmospheric transport models neglecting aggregation tend to overestimate far-range ash cloud concentrations, leading to an overestimation of the risk posed by ash clouds on civil aviation and an underestimation of ash loading in the near field. So far, no plume model tries to predict the formation of ash aggregates in the eruptive column and how it affects the particle grain size distribution erupted at the vent. This can be explained in often identified with terminal settling velocity classes although, strictly, a particle settling velocity class is defined not only by particle size but also by its density and shape. We propose a model 95 for volcanic plumes as a multiphase homogeneous mixture of water (in any phase), entrained air, and n particle classes, including a parameterization for the air entrainment coefficients and a wet aggregation model. Because the governing equations based upon the BPT are not adequate above NBL, we also propose a new semi-empirical model to describe such a region.

100
The steady-state cross-section averaged governing equations for axisymmetric plume motion in a turbulent wind (see Fig. 1) are the following (for the meaining of the used symbols see Tables 1 and   2): whereM = πr 2ρû is the total mass flow rate,P =Mû is the total axial (stream-wise) momentum flow rate, θ is the plume bent over angle with respect to the horizontal (i.e. θ = 90 • for a plume raising vertically),Ê =M (Ĥ +gz + 1 2û 2 ) is the total energy flow rate,Ĥ is the enthalpy flow rate of 125 the mixture,T =T (Ĥ) is the mixture temperature,M a is the mass flow rate of dry air,M w =Mx w is the mass flow rate of volatiles (including water vapour, liquid and ice), h wa is the enthalpy per unit mass of the water in the atmosphere,M i =Mx p f i is the mass flow rate of particles of class i(i = 1 : n), x and y are the horizontal coordinates, z is height, and s is the distance along the plume axis (see Tables 1 and 2 for the definition of all symbols and variables appearing in the manuscript).

130
The equations above derive from conservation principles assuming axial (stream-wise) symmetry and considering bulk quantities integrated over a plume cross-section using a top-hat profile in which a generic quantity φ has a constant valueφ(s) at a given plume cross-section and vanishes outside (here we refer to section-averaged quantities as bulk quantities, denoted by a hat). We have derived these equations by combining formulations from different previous plume models (Netter-135 ville, 1990;Woods, 1993;Ernst et al., 1996;Bursik, 2001;Costa et al., 2006;Woodhouse et al., 2013) in order to include in a single model effects from plume bending by wind, particle fallout and re-entrainment at plume margins, transport of volatiles (water) accounting also for ingestion of ambient moisture, phase changes (water vapour condensation and deposition) and particle aggregation. Equation (2a) expresses the conservation of total mass, accounting in the Right Hand Side 140 (RHS) for the mass of air entrained through the plume margins and the loss/gain of mass by particle fallout/re-entrainment. Equations (2b) and (2c) express the conservation of axial (stream-wise) and radial momentum respectively, accounting on the RHS for contributions from buoyancy (first term), entrainment of air, and particle fallout/re-entrainment. Note that generally the buoyancy term, acting only along the vertical direction z, represents a sink of momentum in the basal gas-thrust jet region 145 (whereρ > ρ a ) and a source of momentum where the plume is positively buoyant (ρ < ρ a ). Equation (2d) express the conservation of energy, accounting on the RHS for gain of energy (enthalpy, potential and kinetic) by ambient air entrainment (first term), loss/gain by particle fallout/re-entrainment (second term), and gain of energy by conversion of water vapour into liquid (condensation) or into ice (deposition). Equations (2e), (2f) and (2g) express, respectively, the conservation of mass of dry 150 air, water (vapour, liquid and ice) and solid particles. The latter set of equations, one for each particle class, account on the RHS for particle re-entrainment (first term), particle fallout (second term) and particle aggregation. Here we have included to terms (A + i and A − i ) that account for the creation of mass from smaller particles aggregating into particle class i and for the destruction of mass resulting from particles of class i contributing to the formation of larger-size aggregates. Finally, Eqs. (2h) to 155 (2j) determine the 3D plume trajectory as a function of the length parameter s. All these equations constitute a set of 9 + n first order ordinary differential equations in s for 9 + n unknowns:M ,P , θ, Ê ,M a ,M w ,M i (for each particle class), x, y and z. Note that, using the definitions ofM -P -Ê, the equations can also be expressed in terms ofû-r-T given the bulk density.
Assuming an homogeneous mixture, the bulk densityρ of the mixture is: wherex p ,x l andx s are, respectively, the mass fractions of particles, liquid water and ice, ρ p is the class-weighted average density of particles (pyroclasts), ρ l and ρ s are liquid water and ice densities, and ρ g is the gas phase (i.e. dry air plus water vapour) density. Under the assumption of mechanical equilibrium (i.e. assuming the same bulk velocityû for all phases and components) is holds that: The enthalpy flow rate of the mixture is a non-decreasing function of the temperatureT , given by: where h v , h l , and h s are, respectively, the enthalpy per unit mass of water vapour, liquid, and ice: where c s = 2108 J K −1 kg −1 is the specific heat of ice, T 0 is a reference temperature, h l0 = 3.337× 175 10 5 J kg −1 is the enthalpy of the liquid water at the reference temperature, c l = 4187 J K −1 kg −1 is the specific heat of liquid water, h v0 = 2.501×10 6 J kg −1 is the enthalpy of vapour water at the reference temperature and c v = 1996 J K −1 kg −1 is the specific heat of vapour water. For convenience, the reference temperature T 0 is taken equal to the temperature of triple point of the water (T 0 = 273.15 K).
The energy and the enthalpy flow rate are related by: For the integration of Eq.(2d) and for evaluating the aggregation rate terms in Eq.(2g), the temper-atureT and the mass fractions of ice (x s ), liquid water (x l ) and vapour (x v ) need to be evaluated.
These quantities are obtained by the direct inversion of Eq.(5), with the use of eqs.(2d) and (7) and by assuming that the pressure inside the plume P is equal to the atmospheric pressure at the same 185 altitude (z).
The model uses a pseudo-gas assumption considering that the mixture of air and water vapour behaves as an ideal gas: where P v and P a are, respectively the partial pressures of the water vapour and of the air in the plume, n v and n a are the molar fractions of vapour and air in the gas phase (n v + n a = 1) and m v = 0.018 kg/mole and m a = 0.029 kg/mole are the molar weights of vapour and air. Following Woods (1993) and Woodhouse et al. (2013), we consider that, if the air-water mixture becomes 195 saturated in water vapour, condensation or deposition occurs and the plume remains just saturated.
This assumption implies that the partial pressure of water vapour P v equals the saturation pressure of vapour over liquid (e l ) or over ice (e s ) at the bulk temperature, and the saturation pressures over liquid and ice are given (in hPa) by (Murphy and Koop, 2005): Equation (9) holds forT ≥ T f and Eq. (10) is valid forT ≤ T f , where T f is the temperature of the triple point of the water (here set at P f = 611.2 Pa, T f = 273.16 K). Therefore, ifT > T f and P v < e l the plume is undersaturated and there is no water vapour condensation (i.e.x v =x w and 205x l =x s = 0). In contrast, if P v ≥ e l , the vapour in excess is immediately converted into liquid and: The vapour and air mass fractions x v and x a are evaluated by combining Eq.(11) and (8b). On the other hand, ifT ≤ T f and P v < e s the plume is undersaturated and there is no water vapour deposition. In contrast, if P v ≥ e s , the vapour in excess is immediately converted into ice and:

215
(P − e s ) n v = e s n â x l = 0 Again, the vapour and air mass fractions x v and x a are evaluated by combining Eq.(12) and (8b).
For the particle re-entrainment parameter f we adopt the fit proposed by Ernst et al. (1996) using data for plumes not affected by wind: where P o = r 2 oû 2 o and F o = r 2 oûoĤo are the specific momentum and thermal fluxes at the vent (s = 0), andĤ o is the enthalpy per unit mass of the mixture at the vent. This expression may overestimate 225 re-entrainment for bent over plumes (Bursik, 2001). Finally, particle terminal settling velocity u si is parameterized as (Costa et al., 2006;Folch et al., 2009): where d i is the class particle diameter and C d is a drag coefficient that depends on the Reynolds number Re = d i u siρ /μ. Several empirical fits exist for drag coefficients of spherical and non-spherical 230 particles (e.g. Wilson and Huang, 1979;Arastoopour et al., 1982;Ganser, 1993;Dellino et al., 2005).
In particular, Ganser (1993) gives a fit valid over a wide range of particle sizes and shapes covering the spectrum of volcanic particles considered in volcanic column models (lapilli and ash): where K 1 and K 2 are two shape factors depending on particle sphericity, Ψ, and particle orientation.

235
Given that the C d depends on Re (i.e. on u s ), Eq. 14 is solved iteratively using a bisection algorithm.
Given a closure equation for the turbulent air entrainment velocity u e , and an aggregation model (defining the mass aggregation coefficients A + i and A − i ), Eqs. (2a) to (2i) can be integrated along the plume axis from the inlet (volcanic vent) up to the neutral buoyancy level. Inflow (boundary) conditions are required at the vent (s = 0) for, e.g., total mass flow rateM o , bent over angle θ o = 240 90 • , temperatureT o , exit velocityû o , fraction of waterx wo , null air mass flow rateM a = 0, vent coordinates (x o ,y o and z o ), and mass flow rate for each particle classM io . The latter is obtained from the total mass flow rate at inflow given the particle grain size distribution at the vent: where f io is the mass fraction of class i at the vent.

Entrainment coefficients
Turbulent entrainment of ambient air plays a key role on the dynamics of jets and buoyant plumes. In the basal region of volcanic columns, the rate of entrainment dictates whether the volcanic jet enters into a collapse regime by exhaustion of momentum before the mixture becomes positively buoyant or whether it evolves into a convective regime reaching much higher altitudes. Early laboratory experiments (e.g. Hewett et al., 1971) already indicated that the velocity of entrainment of ambient air is proportional to velocity differences parallel and normal to the plume axis (see inset in Fig. 1): where α s and α v are dimensionless coefficients that control the entrainment along the stream-wise (shear) and cross-flow (vortex) directions respectively. Note that, in absence of wind (i.e. u a = 0), 255 the equation above reduces to u e = α sû and the classical expression for entrainment velocity of Morton et al. (1956) is recovered. In contrast, under a wind field, both an along-plume (proportional to the relative velocity differences parallel to the plume) and a cross-flow (proportional to the wind normal component) contributions appear. However it is worth noting that Eq. (17) Devenish, 2013). However, more recent experimental (Kaminski et al., 2005) and sensitivity analysis numerical studies (Charpentier and Espíndola, 2005) concluded that pice-wise constant functions are valid only as a first approach, implying that 1D integral models assuming constant entrainment coefficients do not always provide satisfactory results. This has 270 also been corroborated by 3D numerical simulations of volcanic plumes (Suzuki and Koyaguchi, 2013), which indicate that 1D integral models overestimate the effects of wind on turbulent mixing efficiency (i.e. the value of α v ) and, consequently, underestimate plume heights under strong wind fields. For example, recent 3D numerical simulation results for small-scale eruptions under strong wind fields suggest lower values of α v , in the range 0.1 − 0.3 (Suzuki and Koyaguchi, 2015). For 275 this reason, besides the option of constant entrainment coefficients, FPLUME allows considering also a parameterization of α s and α v based on the local Richardson number. In particular, we use the empirical parameterization of Kaminski et al. (2005) and Carazzo et al. (2006Carazzo et al. ( , 2008a) that describes α s for jets and plumes as a function of the local Richardson number as: where A(z s ) is an entrainment function depending on the dimensionless length z s = z/2r o (r o is the vent radius) and Ri = g(ρ a −ρ)r/ρ aû 2 is the Richardson number. Beside the local Richardson number, the entrainment coefficient α s depends on plume orientation (e.g. Lee and Cheung, 1990;Bemporad, 1994), therefore we modify Eq. 18 as: Moreover in order to use a compact analytical expression and extend it to values of z s ≤ 10 we fitted the experimental data of Carazzo et al. (2006Carazzo et al. ( , 2008b considering the following empirical function: and in order to extrapolate to low z s we multiply A(z s ) for the following function h(z s ) that affects the behavior only for small values of z s : where c i are dimensionless fitting constants. Best-fit results and entrainment functions resulting from fitting Eqs. (20a)-(20c) are shown in Table 3 whereû o is the mixture velocity at the vent andū a is the average wind velocity. For illustrative purposes, Fig. 3 shows the entrainment coefficients α s and α v predicted by Eqs. (19) and (21) for weak and strong plume cases under a prescribed wind profile. It is important stressing that air entrainment rates play a first-order role on eruptive plume dynamics and a simple description in terms of entrainment coefficients, both assuming them as empirical constants or describing them like in 305 (18), represents an over-simplication of the real physics characterizing the processes. A better quantification of entrainment rates is one of the current main challenges of the volcanological community (see Costa et al., 2015, and references therein).

Modelling of the Umbrella Region
The umbrella region is defined as the upper region of the plume, from about the NBL to the top of the 310 column. This region can be dominated by fountaining processes of the eruptive mixture that reaches the top of the column, dissipating the excess of momentum at the NBL, and then collapsing as a gravity current (e.g. Woods and Kienle, 1994;Costa et al., 2013). The 1-D BPT should not be extended to this region because it assumes that the mixture still entrains air with the same mechanisms as below NBL and, moreover, predicts that the radius goes to infinity towards the top of the column.
For these reasons, we describe the umbrella region adopting a simple semi-empirical approximation.
In the umbrella region (from the NBL to the top of the column), we neglect air entrainment and assume that the mixture is homogeneous, i.e. the content of air, water vapour, liquid water, ice, and total mass of particles do not vary with z. Pressure P (z) is considered equal to the atmospheric pressure P a (z) evaluated at the same level, whereas temperature decreases with z due to the adiabatic 320 cooling: As a consequence, the density of the mixture varies accordingly. The total height of the volcanic plume H t , above the vent, is approximated as (e.g. Sparks, 1986): where c H is a dimensionless parameter (typically c H = 1.32), H b is the height of the Neutral Buoyancy Level (above the vent) and r o the radius at the vent. Between H b and H t , the coordinates x and y of the position of the plume centre and the plume radius r are parameterized as a function of The position of the plume centre is assumed to vary linearly with the same slope at the NBL, whereas the effective plume radius is assumed to decrease as a Gaussian 330 function: where x b , y b , r b are, respectively, the coordinates x and y of the center of the plume and the plume 335 radius at the NBL, and Finally, assuming that the kinetic energy of the mixture is converted to potential energy, the vertical velocity is approximated to decrease as the square root of the distance from the NBL: where u zb is the vertical velocity of the plume at the NBL.

340
Although the proposed empirical parameterization of the region above the NBL is qualitatively consistent with the trends predicted by 3D numerical models (Costa et al., 2015), a more rigorous description requires further research.

Plume Wet Aggregation Model
Particle aggregation can occur inside the column or in the ash cloud during subsequent atmospheric 345 dispersion (e.g. Carey and Sigurdsson, 1982;, thereby affecting the sedimentation dynamics and deposition of volcanic ash. Our model explicitly accounts for aggregation in the plume by adding source (A + i ) and sink (A − i ) terms for aggregates and aggregated particles in their respective particle mass balance Eqs. (2g) and by modifying the settling velocity of aggregates.
Given the complexity of aggregation phenomena, not yet fully understood, we consider only the oc-350 currence of wet aggregation and neglect dry aggregation mechanisms driven by electrostatic forces or disaggregation processes resulting from particle collisions that can break and decompose aggregates. Costa et al. (2010) and Folch et al. (2010) proposed a simplified wet aggregation model in which particles aggregate on a single effective aggregated class characterized by a diameter d A (i.e. aggregation only involves particle classes having an effective diameter smaller than d A , typically in 355 the range 100-300 µm). Obviously the assumption that all particles aggregate into a single particle class is simplistic and considering a range of aggregating classes would be more realistic. However, there are no quantitative data available for such a calibration. Hence, considering this assumption it follows that: where k is the (given) index of the aggregated class and the sum over j spans all particle classes having diameters lower than d A . The mass of particles of class i (d i < d A ) that aggregate per unit of time and length in a given plume cross-section is: whereṅ i is the number of particles of class i that aggregate per unit volume and time, estimated as: In the expression above, N i is the number of particles of diameter d i in an aggregate of diameter d A , andṅ tot is the total particle decay per unit volume and time. Costa et al. (2010) considered that N i is given by a semi-empirical fractal relationship (e.g. Jullien and Botet, 1987;Frenklach, 2002;Xiong and Friedlander, 2001): where k f is a fractal pre-factor and D f is the fractal exponent. Costa et al. (2010);Folch et al. (2010) assumed constant values for k f and D f that were calibrated by best-fitting tephra deposits from 18 May 1980 Mount St. Helens and 17-18 September 1992 Crater Peak eruptions. However, for the granulometric data from these deposits they used a cut-off considering only particles larger than 375 about 10µm, for which the gravitational aggregation kernel dominates. This poses a problem if one wants to extend the granulometric distribution to include micrometric and sub-micrometric particles, for which the Brownian kernel is the dominant one (it is known that Brownian particle-particle interaction has typical values of D f ≈ 2, with values ranging between 1.5 and 2.5, e.g. Xiong and Friedlander (2001)). Actually, preliminary model tests involving micrometric and sub-micrometric 380 particle classes considering constant values for D f and k f have revealed a strong dependency of results (fraction of aggregated mass) on both granulometric cut-off and bin width (particle grain size discretization). In order to overcome this problem, we assume a size-dependent fractal exponent as: where D (2002): Figure 4 shows the values of D f (d) and k f (d) predicted by Eqs. (32) and (33) for a range of D f o .

390
We have performed different tests to verify that, in this way, the results of the aggregation model become much more robust independently of the distribution cut-off (Φ min = 8, 10, 12) and bin width (∆Φ = 1, 0.5, 0.25), with maximum differences in the aggregated mass laying always below 10%.
The total particle decay per unit volume and timeṅ tot is given by: where α m is a mean (class-averaged) sticking efficiency, φ is the solid volume fraction, n tot is the total number of particles per unit of volume that can potentially aggregate andf is a correction factor that accounts for conversion from gaussian to top-hat formalism (see Appendix A for details). The expression above comes from integrating the collection kernel over all particle sizes, and involves the product of the (averaged) sticking efficiency times the collision frequency function accounting where k b is the Boltzmann constant andμ is the mixture dynamic viscosity (≈ air viscosity at the 405 bulk temperatureT ). The term A T I derives from the collision kernel due to turbulence as result of inertial effects β T I,ij (e.g. Pruppacher and Klett, 1996;Jacobson, 2005): whereν is the mixture kinematic viscosity and is the dissipation rate of turbulent kinetic energy, computed assuming the Smagorinsky-Lilly model: where k s ≈ 0.1 − 0.2 is the constant of Smagorinsky. The term A S derives from the collision kernel due to laminar and turbulent fluid shear β S,ij : where Γ S is the fluid shear, computed as: Finally, the term A DS derives from the differential sedimentation collision kernel β DS,ij (e.g. Costa et al., 2010): where u si denotes the settling velocity of particle class i. Note that, with respect to the original 420 formulation of Costa et al. (2010), using the same approach and approximation, we have included the additional term A T I due to the turbulent inertial kernel that, thanks to the similarity between Eqs. (40) and (36), can be easily derived. Once these kernels are integrated, expressions for the terms in Eq. (34) yield: is the diameter to volume fractal relationship and v j is the particle volume.
Note that for spherical particles in the Euclidean space (D f = 3) v j = πd 3 j /6 and ξ = (6/π) 1/3 . The total number of particles per unit of volume available for aggregation is related to particle class mass concentration at each section of the plumeĈ j and can be estimated as (see Appendix B): where d aj and d bj are the particle diameters of the limits of the interval j and: Finally, the class-averaged sticking efficiency α m appearing in (34) is computed as: where f k is the particle class mass fraction, and α ij is the sticking efficiency between the classes i and j. In presence of a pure ice phase we assume that ash particles stick as ice particles (α m = 0.09).
In contrast, in presence of a liquid phase, the aggregation model considers: where St cr = 1.3 is the critical Stokes number, q = 0.8 is a constant, and St ij is the Stokes number 445 based on the binder liquid (water) viscosity: where Obviously, our aggregation model requires the presence of water either in liquid or solid phases, i.e.   3. Estimate the total number of particles per unit of volume available for aggregation n tot depending onĈ j using Eq. (42).
5. Compute the total particle decay per unit volume and timeṅ tot using Eq. (34) depending also on the solid volume fraction. Eq. (31) assuming size-dependent fractal exponent D f and pre-factor k f . 7. Compute class particle decayṅ i using Eq. (30).
8. Finally, compute the mass sink term for each aggregating class A − i using Eq. (29) and the mass source term A + i for the aggregated class using Eq. (28) to introduce these terms in the 500 particle class mass balance equations Eqs. (2g).

Test Cases
As we mentioned above, here we apply FPLUME to two eruptions relatively well characterized by previous studies. In particular we consider the strong plume formed during 4 April 1982 by El Chichón 1982 eruption (e.g. Sigurdsson et al., 1984;Bonasia et al., 2012) and the weak plume 505 formed during the 6 May 2010 Eyjafjallajökull eruption (e.g. Bonadonna et al., 2011;Folch, 2012).

Phase-B El Chichón 1982 eruption
El Chichón volcano reawakened in 1982 with three significant Plinian episodes occurring during March 29th (phase A) and April 4th (phases B and C). Here we focus on the second major event, starting at 01:35 UTC on April 4th and lasting nearly 4.5h (Sigurdsson et al., 1984). Bonasia et al. 510 (2012) used analytical (HAZMAP) and numerical (FALL3D) tephra transport models to reconstruct ground deposit observations for the three main eruption fallout units. Deposit best-fit inversion results for phase-B suggested column heights between 28 and 32 km (above vent level, a.v.l.) and a total erupted mass ranging between 2.2×10 12 and 3.7×10 12 kg. Considering a duration of 4.5h, the resulting averaged mass eruption rates are between 1 × 10 8 and 2.3 × 10 8 kg/s. TGSD of phases B 515 and C were estimated by Rose and Durant (2009) weighting by mass, by isopach volume and using the Voronoi method. Bonasia et al. (2012) found that the reconstruction of the deposits is reasonably achieved taking into account the empirical Cornell aggregation parameterization (Cornell et al., 1983). In this simplistic approach, 50% of the 63-44 µm ash, 75% of the 44-31 µm ash and 100% of the less than 31 µm ash are assumed to aggregate as particles with a diameter of 200 µm and density 520 of 200 kgm −3 . Note that here, as in previous studies , we use a modified version of Cornell et al. (1983) parameterization that assumes that 90% and not 100% of the particle smaller than 31 µm fall as aggregates.
We use this test case to verify whether FPLUME can reproduce results from these previous stud-  Figure 5 shows the wind profile and the FPLUME results for bulk velocity 530 and plume radius. The model predicts a total plume height of 28 km (a.s.l.), a mass eruption rate of 2.7 × 10 8 kg/s, and a total erupted mass of 4.4 × 10 12 kg. These values are consistent but slightly higher than those from previous studies (Bonasia et al., 2012). Regarding the aggregation model, we did several sensitivity runs to look into the impact of the fractal exponent D f o on the fraction of aggregates, ranging this parameter between 2.85 and 3.0 at 0.01 steps values (see Fig. 6). As antici-535 pated in the original formulation Folch et al., 2010), the results of the aggregation model are sensitive to this parameter. Values of D f o = 2.96 fit very well the total mass fraction of aggregates predicted by Cornell but not the fraction of the aggregating classes (Fig. 7 b). In contrast, we find a more reasonable fit with D f o = 2.92, although in this case the relative differences for the total mass fraction of aggregates are of about 15%, with our model under-predicting with respect to 540 Cornell (Fig. 7 a).
A clear advantage of a physical aggregation model of ash particles inside the eruption column, with respect an empirical parameterization like (Cornell et al., 1983), is that it allows estimating the fraction of very fine ash that escapes to aggregation processes and is transported distally within the cloud. As we mentioned above, based on the features of the observed deposits, Cornell et al. (1983) 545 proposed that 100% of particles smaller than 31 µm fall as aggregates that is quite reasonable as most of fine ash falls prematurely. However assessing the small mass fraction of fine ash that escapes to aggregation processes is crucial for aviation risk mitigation and for comparing model simulations with satellite observations. For example, in the case of El Chichón 1982 eruption, for D f o = 2.92, the model predicts that ≈ 10% of fine ash between 20 and 2 µm in diameter escapes to aggregation 550 processes. This value is an order of magnitude larger than that estimated by Schneider et al. (1999) using TOMS and AVHRR but we need to consider that we do not account for dry aggregation that can be dominant for very fine particles.

6 May 2010 Eyjafjallajökull eruption phase
The infamous April-May 2010 Eyjafjallajökull eruption, that disrupted the European North Atlantic 555 region airspace (e.g. Folch, 2012), was characterized by a very pulsating behavior, resulting on nearly continuous production of weak plumes that oscillated on height between 2 and 10 km (a.s.l.) during the 39 day-long eruption (e.g. Gudmundsson et al., 2012). During 4-8 May, Bonadonna et al. (2011) performed in-situ observations of tephra accumulation rates and PLUDIX Doppler radar measurements of settling velocities at different locations which then used to determine erupted mass, mass 560 eruption rates and grain size distributions. The authors estimated a TGSD representative of 30 min of eruption by combining ground-based grain-size observations (using a Voronoi tessellation technique) and ash mass retrievals (7-9Φ particles) from MSG-SEVIRI satellite imagery for 6 May between 11:00 and 11:30 UTC. On the other hand, they also report the in-situ observation of sedimentation of dry and wet aggregates falling as particle clusters and poorly structured and liquid accretionary 565 pellets (AP1 and AP3 according to Brown et al. (2012) nomenclature). Bonadonna et al. (2011) did also grain-size analyses of collected aggregates using scanning electron microscope (SEM) images.
The combination of all these data allowed them to determine how the original TGSD was modified by the formation of different types of aggregates (see Fig. 8). The total mass fraction of aggregates was estimated to be about 25% with aggregate sizes ranging between 1Φ (500 µm) and 4Φ (62.5 570 µm). These results constitute a rare and valuable dataset to test the aggregation model implemented in FPLUME. However, several challenges can be anticipated. First, our model assumes a single ag-gregated class and, as a consequence, we may expect to reproduce only the total mass fraction of aggregates, but not to match the resulting mass fraction distribution class by class. Second, the proportion of dry versus wet aggregates is unknown and, wet aggregation could have occurred within 575 the plume but also by local rain showers that scavenged coarse particles (Bonadonna et al., 2011), moreover the presence of meteoritic water in the plume (not considered here) could significantly enhance aggregation. For these reasons, we aim to capture the correct order of magnitude of total mass fraction of ash that went into aggregates.
Preliminary simulations using time-averaged plume heights of 3.5-4.5 km (a.v.l.) did not result in This may suggest that wet aggregates could have formed within the plume not continuously but during sporadic higher-intensity column pulses. In order to check this possibility, we performed a parametric study to compute the total mass fraction of wet aggregates as function of mass flow rate that controls the value of the column height.

590
Input values for FPLUME are summarized in Table 5. The wind profile (see Fig. 9) was extracted from the ERA-Interim re-analysis dataset interpolating values at the vent coordinates. As shown in Fig. 10, 10% in mass of wet aggregates is predicted by our model for column heights ranging between 6 and 7 km (a.v.l.) and 20% for column heights from 7.2 to 8.3 km (a.v.l.). For the considered input parameters and ambient conditions (wind and moisture profile), we observed the formation of 595 a window in the plume containing liquid water only for column heights above 5.3 km (a.v.l.). For illustrative purposes, Fig.11 shows the resulting grain size distribution for a column height of 6.5 km (a.v.l.) and two different values of the fractal exponent D f . As anticipated, the model can predict the total mass fraction of aggregates, but an error (< 10%) exists for some particular classes. However, it should be kept in mind that mass fraction of aggregates is not controlled only by the eruptive col-600 umn height but depends on several variables such as particle concentration (that is a function of the mass flow rate), presence of liquid water (that can form above a given level depending on the local meteorological conditions), etc.

Conclusions
We presented FPLUME, a 1D cross-section averaged volcanic plume model based on the BPT 605 that accounts for plume bending by wind, entrainment of ambient moisture, effects of water phase changes, particle fallout and re-entrainment, a new parameterization for the air entrainment coeffi-cients and an ash wet aggregation model based on Costa et al. (2010). Given conditions at the vent (mixture exit velocity, temperature and magmatic water content) and a wind profile, the model can solve for plume height given the eruption rate or vice-versa. FPLUME can also be extended above

Code availability
The code FPLUME-1.0 is available under request for research purposes.

625
Appendix A: Correction factorf for mass distribution for top-hat versus Gaussian formalism Denoting with R the top-hat radius of the plume and with b the Gaussian length scale the relationship between them can be written as (e.g. Davidson, 1986): Assuming a Gaussian profile for the concentration, C(r), the mean value between r = 0 (where the 630 concentration is maximum) and r = R is: that impliesĈ = 0.25C 0 . Following similar calculations we have also: Therefore if we use average (top-hat) variables in Eq. (34) we need to keep in mind that concentration appears in the nonlinear terms and therefore we should use the following correction factors: and so on ( · denotes the average using the top-hat filter, e.g.Ĉ = C ). Because terms in Eq. (34) scale with concentration with a power of two we need to account for a correction factorf =f 2 . The factorf can be also used to correct underestimation of Eulerian time scale with respect Lagrangian time scale (e.g. Dosio et al., 2005).