Modelling the size distribution of aggregated volcanic ash and implications for operational atmospheric dispersion modelling

We have developed an aggregation scheme for use with the Lagrangian atmospheric transport and dispersion model NAME, which is used by the London Volcanic Ash Advisory Centre (VAAC) to provide advice and guidance on the location of volcanic ash clouds to the aviation industry. The aggregation scheme uses the fixed pivot technique to solve the Smoluchowski coagula5 tion equations to simulate aggregation processes in an eruption column. This represents the first attempt at modelling explicitly the change in the grain size distribution (GSD) of the ash due to aggregation in a model which is used for operational response. To understand the sensitivity of the output aggregated grain size distribution (AGSD) to the model parameters we conducted a simple parametric study and scaling analysis. We find that the modelled AGSD is sensitive to the density distribution and grain size distribution assigned to the non-aggregated ash at the source. Our ability to accurately forecast the long-range transport 10 of volcanic ash clouds is, therefore, still limited by real-time information on the physical characteristics of the ash. We assess the impact of using the AGSD on model simulations of the Eyjafjallajökull 2010 ash cloud, and consider the implications for operational forecasting. Using the time-evolving AGSD at the top of the eruption column to initialise dispersion model simulations had little impact on the modelled extent and mass loadings in the distal ash cloud. Our aggregation scheme does not account for the density of the aggregates; however, if we assume that the aggregates have the same density of single grains 15 of equivalent size the modelled extent of the Eyjafjallajökull ash cloud with high concentrations of ash, significant for aviation, is reduced by ∼3%. If we assume that the aggregates have a lower density (500 kg m−3) than the single grains of which they are composed and make-up 75% of the mass in the ash cloud the extent is 1.2 times larger. 1 https://doi.org/10.5194/acp-2021-254 Preprint. Discussion started: 31 May 2021 c © Author(s) 2021. CC BY 4.0 License.

tion equations to simulate aggregation processes in an eruption column. This represents the first attempt at modelling explicitly the change in the grain size distribution (GSD) of the ash due to aggregation in a model which is used for operational response.
To understand the sensitivity of the output aggregated grain size distribution (AGSD) to the model parameters we conducted a simple parametric study and scaling analysis. We find that the modelled AGSD is sensitive to the density distribution and grain size distribution assigned to the non-aggregated ash at the source. Our ability to accurately forecast the long-range transport 10 of volcanic ash clouds is, therefore, still limited by real-time information on the physical characteristics of the ash. We assess the impact of using the AGSD on model simulations of the Eyjafjallajökull 2010 ash cloud, and consider the implications for operational forecasting. Using the time-evolving AGSD at the top of the eruption column to initialise dispersion model simulations had little impact on the modelled extent and mass loadings in the distal ash cloud. Our aggregation scheme does not account for the density of the aggregates; however, if we assume that the aggregates have the same density of single grains 15 of equivalent size the modelled extent of the Eyjafjallajökull ash cloud with high concentrations of ash, significant for aviation, is reduced by ∼3%. If we assume that the aggregates have a lower density (500 kg m −3 ) than the single grains of which they are composed and make-up 75% of the mass in the ash cloud the extent is 1.2 times larger.

20
In volcanic plumes ash can aggregate together, bound by hydro-bonds and electrostatic forces. Aggregates typically have diameters > 63 µm (Brown et al., 2012) and their fall velocity differs to that of the single grains of which they are composed (Lane et al., 1993;James et al., 2003;Taddeucci et al., 2011;Bagheri et al., 2016). Neglecting aggregation in atmospheric dispersion models could, therefore, lead to errors when modelling the rate of removal of ash from the atmosphere, and consequently inaccurate forecasts of the concentration and extent of volcanic ash clouds used by civil aviation for hazard assessment (e.g. 25 Folch et al., 2010;Mastin et al., 2013;Beckett et al., 2015;Mastin et al., 2016).
The theoretical description of aggregation is still far from fully understood, mostly due to the complexity of particle-particle interactions within a highly turbulent fluid. There have been several attempts to provide an empirical description of the aggregated grain size distribution (AGSD) by assigning a specific cluster settling velocity to fine ash  or fitting the distribution used in dispersion models to observations of tephra deposits retrospectively (e.g. Cornell et al.,30 1983; Bonadonna et al., 2002;Mastin et al., 2013Mastin et al., , 2016. Cornell et al. (1983) found that by replacing a fraction of the fine ash with aggregates which had a diameter of 200 µm they were able to reproduce the observed dispersal of the Campanian Y-5 ash. Bonadonna et al. (2002) found that the ash deposition from co-pyroclastic density currents and the plume associated with both dome collapses and Vulcanian explosions of the 1995-1999 eruption of Soufrière Hills Volcano (Montserrat) were better described by considering variation in the aggregate size and in the grain-size distribution within individual aggregates. 35 Mastin et al. (2016) determined optimal values for the mean and standard deviation of input AGSDs for ash from the eruptions of Mount St Helen's, Crater Peak (Mount Spurr), Ruapehu and Mount Redoubt using the Ash3d model. They assumed that the aggregates had a Gaussian size distribution and found that for all the eruptions the optimal mean aggregate size was 150-200 µm.
There have been only a few attempts to model the process of aggregation explicitly. Veitch and Woods (2001) were the first 40 to represent aggregation in the presence of liquid water in an eruption column using the Smoluchowski Coagulation Equations (SCE) (Smoluchowski, 1916). Textor et al. (2006a, b) introduced a more sophisticated aggregation scheme to the Active Tracer High-resolution Atmospheric Model (ATHAM), also designed to model eruption columns, which included a more robust representation of microphysical processes and simulated the interaction of hydrometeors with volcanic ash. They suggest that wet rather than icy ash has the greatest sticking efficiency and that aggregation is fastest within the eruption column where 45 ash concentrations are high and regions of liquid water exist. More recently microphysical based aggregation schemes which represent multiple collision mechanisms have been introduced to atmospheric dispersion models FALL3D Folch et al., 2010) and WRF-Chem (Egan et al., 2019), and an eruption column model, FPLUME (Folch et al., 2016). They all use an approximate solution of the SCE which assumes that aggregates can be described by a fractal geometry and particles aggregate onto a single effective aggregate class defined by a prescribed diameter. 50 Here we introduce an aggregation scheme coupled to a one-dimensional steady state buoyant plume model which uses a discrete solution of the SCE based on the Fixed Pivot Technique (FPT) (Kumar and Ramkrishna, 1996). As such we are able to model explicitly the evolution of the AGSD with time in the eruption column. We have integrated our aggregation scheme into 2010, and consider the implications of using an AGSD for operational forecasting. We discuss the results in Section 6, before the conclusions are presented in Section 7.

The Aggregation Scheme
We use a one-dimensional steady state buoyancy model, where mass, momentum and total energy are derived for a control 70 volume, and time variations are assumed to be negligible (Devenish, 2013(Devenish, , 2016. It combines the effects of moisture (liquid water and water vapour) and the ambient wind, and includes the effects of humidity and phase changes of water, on the growth of the plume. The governing equations are given by: where s is the distance along the plume axis, M z = Q m w p is the vertical momentum flux, M i = (u pi − U i ) Q m is the horizontal momentum flux relative to the environment, H = c pp T Q m is the enthalpy flux, Q t = Q m n t is the total moisture flux within the plume, and Q m = ρ p πb 2 v p is the mass flux. The meaning of the symbols used throughout are given in Tables 1 and 2. The 85 entrainment rate depends on the ambient and plume densities, and when the plume is rising buoyantly is given by: where ρ p is the plume density: and u e is the entrainment velocity: 90 u e = (k s |∆u s |) 1.5 + (k n |∆u n |) Here two entrainment mechanisms are considered, one due to velocity differences parallel to the plume axis (u s ) and one due to the velocity differences perpendicular (u n ) to the plume axis, k s and k n are the entrainment coefficients associated with each respective entrainment mechanism (note k s is given the symbol α and k n the symbol β in Devenish (2013)).
As aggregation is controlled by the amount of available water it is essential that we adequately consider the entrainment of 95 water vapour, its condensation threshold, and phase changes from water vapour to ice and liquid water, and vice versa. As such, we have modified the scheme presented by Devenish (2013) to introduce an ice phase. Ice is produced whenever T < 255K, the critical temperature in the presence of volcanic ash, following Durant et al. (2008); Costa et al. (2010);Folch et al. (2016).
It is assumed that there is no source liquid water or ice flux, given the high temperatures, and that there is no entrainment of ambient liquid water (only water vapour). Liquid water condensate and ice are formed whenever the water vapour mixing ratio 100 (r v ) is larger than the saturation mixing ratio (r s ), which is determined using the Clausius-Clapeyron equation: where ε = 0.62 is the molecular mass of water vapour to dry air, p d is the dry ambient pressure and e s is the saturation vapour pressure, which for liquid is given by a modification of Tetens' empirical formula (Emanuel, 1994, pg. The mass fractions of water (n l ) and ice (n ice ) can then be expressed as: where n t is the total moisture fraction (n t = n v + n l,ice ) and n d is the dry gas fraction. It is assumed that any liquid condensate and ice that forms remains in the plume and thus the total water content is conserved.
The SCE are solved using the FPT, which transforms a continuous domain of masses (whilst conserving mass) into a discrete space of sections, each identified by the central mass of the bin, i.e. the pivot. The growth of the aggregates is described by the 115 sticking efficiency between the particles and their collision frequencies. The approach is computationally efficient but can be affected by numerical diffusion if the number of bins is too coarse compared to the population under analysis. The coupling of the FPT with the one-dimensional buoyant plume model is applied at the level of the mass flux conservation equations. The mass flux is modified such that the mass fractions of the dry gas (n d ), total moisture (n t , which is the mass fraction of vapour (n v ) only, as neither liquid water or ice are entrained) and solid phases (n s ) are treated separately: where n d + n t + n v = 1. As there is no entrainment of solids, Eq.14 can be expressed as: We assume a discretized GSD composed of N bins , where the mass fractions of a given size (x i ) are divided across a set of bins, such that N bins i=1 x i = 1. Assuming each size shares an amount of mass flux that is proportional to x i , Eq.16 becomes: where we used the linearity of the sum with respect to the derivative operator. This is the continuity equation for solid mass flux in the case of a steady-state process. The continuity equation can be seen as a set of N bins equations, one for each i-th 130 section, where aggregation is then taken into account by introducing source (B i ) and sink (D i ) terms in the right-hand side (rhs) of Eq.17. The continuity equation for the i-th bin then becomes: In the FPT the source term B i states that a given particle of the i-th section can be created when the sum of the masses m sum of two interacting particles k and j is between the pivots [i-1, i] and [i, i+1]. A fraction of m sum is then proportionally attributed 135 to the i-th pivot according to how close the mass m sum is to m i . The redistribution of m sum among the bins is done in such a way that the mass is conserved by definition. The sink term D i , on the other hand, is just related to the number of collisions where N i is the number of particles of a given mass per unit volume: K k,j is the aggregation kernel between particles belonging to bins k and j respectively, and δ kj is the Kronecker delta function.

145
As such the SCE have been transformed into a set of Ordinary Differential Equations (ODEs) which are solved for each bin representing the i-th mass. The process of aggregation between two particles of mass m k and m j , at a given location s along the central axis of the plume, depends on the aggregation kernel (K k,j ), which can be expressed in terms of the sticking efficiency (α k,j ) and the collision rate (β k,j ) of particles: 150 where α k,j is a dimensionless number between 0 and 1 which quantifies the probability of the particles successfully sticking together after a collision. β k,j describes the average volumetric flow of particles (m 3 s −1 ) involved in the collision between particles k and j. We consider five different mechanisms, following Costa et al. (2010): Brownian motion (β B k,j ), interactions due to the differential settling velocities between the particles (β DS k,j ), and the interaction of particles due to turbulence: the inertial turbulent kernel (β T I k,j ) and the fluid shear, both laminar β LS k,j and turbulent β T S k,j : where d k and d j are the diameters and V k and V j are the sedimentation velocities of the colliding particles: where C D is the drag coefficient and Re is the Reynolds number.: 170 and the sedimentation velocity is evaluated using an iterative scheme following Arastoopour et al. (1982). The laminar fluid shear is taken to be Γ = |dw p /dz|. The dissipation rate of turbulent kinetic energy per unit mass ( ) is constrained by the parameters controlling the large-scale flow, the magnitude of velocity fluctuations (about 10% of the axial plume velocity) and the size of the largest eddies, which we take to be the plume radius (Textor and Ernst, 2004): The total contribution from collisions due to each of the different mechanisms is represented by a linear superposition of each of the kernels (taking the maximum of the Shear Laminar and Shear Turbulent kernels): The different collision mechanisms are evaluated at each position s along the central axis of the plume.
We assume that ash can stick together due to the presence of a layer of liquid water on the ash, following Costa et al. (2010).

180
In this framework the energy involved in the collision of particles k and j, identified from the relative kinetic energy of the bodies (i.e. rotations are not taken into account), can be parametrized in terms of the collision Stokes number (St v ): which is a function of the average density of the two colliding particles (ρ), the liquid viscosity (µ l ) and the relative velocities between the colliding particles (U r ), here approximated as: where h is the thickness of the liquid layer and h a is the surface asperity (Liu et al., 2000;Liu and Litster, 2002). Unfortunately this information is poorly constrained for volcanic ash. Costa et al. (2010) propose the following parameterization for the sticking efficiency: using the experimental data of Gilbert and Lane (1994), which considered particles with diameters between 10 and 100 µm, and set St cr = 1.3 and q = 0.8 (see Figure 12 in Gilbert and Lane (1994) and Figure 1 in Costa et al. (2010)).
The influence of the ambient conditions, such as the relative humidity, on liquid bonding of ash aggregates still remains poorly constrained. Moreover, when trying to derive environmental conditions from one-dimensional plume models, it should be remembered that this description of a 3-dimensional turbulent flow simply represents an average of the flow conditions, 200 and lacks details on local 'pockets' of liquid water due to clustering of the gas mixture (Cerminara et al., 2016b). In these local regions the concentration of water vapour can be high enough to reach the saturation condition and trigger the formation of liquid water. Further, aggregation can occur even when the bulk value of the relative humidity is relatively low (Telling and Dufek, 2012;Telling et al., 2013;Mueller et al., 2016). As such we allow sticking to occur in regions where the relative humidity is < 100% and liquid water is not yet present in the one-dimensional description of the plume, and we scale the 205 sticking efficiency (α k,j ) by the relative humidity: In the presence of ice we assume that the sticking efficiency is constant and α k,j = 0.09, following Costa et al. (2010) and Field et al. (2006).  Mass fraction of water vapour -Default n 0 v = 0.00 n g Mass fraction of gas n g = n d + n v n s Mass fraction of solids -n t Mass fraction of total moisture content -  To consider the influence of uncertainty on the source and internal model parameters on the simulated AGSD we have conducted a simple sensitivity study whereby the input parameters are varied one at a time. As such we assess the difference between the simulated output using the set of default parameters (the control case) from a perturbed case. This approach assumes model variables are independent when considering the effects of each on model predictions.  Table 3, along with the range of values for each parameter considered in the sensitivity study. Values are based 220 on the existing literature and the sources used are also listed in Table 3. The scheme is initialized with a GSD with diameters ranging from 1 µm to 16 mm. Bins are defined on the Phi Scale, where the Phi diameter is calculated as the negative logarithm to the base 2 of the particle diameter in millimeters (Krumbein, 1938). The mass is distributed uniformly across the bins such that 50% is on grains with diameter ≤ 125 µm and 36% of the mass is on grains with diameter ≤ 32 µm. The output AGSD at the top of the plume, defined as the point at which W p < 0, is assessed. Given the nature of the SCE, the aggregation scheme 225 does not track explicitly the mass fraction of aggregates versus single grains within a given size bin. Instead we consider how the mode of the output AGSD varies, and compare the mass fraction on particles with diameter ≤ 32 µm (m 32 ), which predominantly lose mass to larger aggregates, for each sensitivity run.
First, we consider how the AGSD varies as conditions within the plume change over time given the local meteorological and eruption conditions (plume height). Figure 1 shows the relative humidity (RH), temperature (T ) and the fractions of liquid 230 water (n l ), water vapour (n v ) and ice (n i ) with height along the plume axis at different times during the eruption. Note that the maximum height of the modelled plume axis, when the plume is bent over as in this case, is the maximum observed plume height minus the plume radius (Mastin, 2014;Devenish, 2016). At 19:00 UTC on the 04/05/2010 the maximum observed plume height is 7000 m asl, liquid water starts to form at 4122 m asl, but no ice forms in the plume. At 12:00 UTC on the 05/05/2010 the observed maximum plume height is lower, reaching just 5500 m asl, liquid water is present from 3684 m asl 235 and again no ice is formed. However, at 13:00 UTC on the 06/05/2010 no liquid water forms in the plume, only ice, which is present from 5882 m asl and the maximum observed plume height is 10,000 m asl. At 12:00 UTC on the 07/05/2010 the maximum observed plume height is 5500 m asl and there is neither liquid water or ice in the plume, only water vapour. Figure 2 shows the output AGSDs for each of these times, compared to the input GSD. We find that in all the cases considered the mode of the AGSD is always the same, with most of the mass now residing in the 125 -250 µm bin. When ash spends more time in The mode and m 32 of the simulated AGSD for each sensitivity run output at 19:00 04/05/2010 are listed in Table 4. Using 245 the control parameters the mode of the AGSD lies at 125 -250 µm and m 32 is 32% at this time (c.f. Figure 2a). The aggregation scheme is sensitive to the values assigned to the sticking parameters (St cr and q) and the parameters which define the particle characteristics, the input GSD and the particle density (note that here all model particles are assigned the same density, as suchρ = ρ s ). Figure 3 shows how the cumulative distribution of the AGSD changes as these parameters are varied within their known ranges. The parameters used to set the sticking efficiency between the particles (St cr and q) are currently poorly is 24% and the modal bin is 500 -1000 µm (Figure 3c). We find that using a relatively coarse input GSD (from the eruption of Hekla 1991) there is very little aggregation; there is no change in m 32 or the modal grain size (Figure 3d) from the input  Woodhouse et al. (2016) perpendicular ( The AGSD shows little sensitivity to the model values assigned to define the plume conditions within the ranges investigated: the entrainment parameters (k s and k n ), the initial mass fraction of dry gas and water vapour (n 0 d ,n 0 v ), the plume temperature 260 at the source (T 0 ) or the source mass flux (Q m ) (see Table 4 and Supplementary Figures S4 -S7). When we consider that the mass flux may have an order of magnitude uncertainty, and vary the input mass flux to the aggregation scheme by a factor of 10 m 32 varies by just 1%.    (c.f. Figure 2a). Note that the blue lines represent simulations using the control values. and β k,j in turn. Starting with the collision rate we can write Eq. 32 as where are taken to be constant (including ρ s ). Here we have assumed that the particles settle with Stokes' terminal velocity (i.e. we 270 neglect the second term on the rhs of Eq. 29); this will lead to quantitative discrepancies with the collision kernel calculated in Section 3 for larger diameters but the qualitative behaviour will be correct. We have also assumed that β T S i,j > β LS i,j ; this assumption does not affect our conclusions below. Since aggregation is associated with the presence of liquid water or ice As in Section 3 we restrict attention to diameters in the range [1, 10 4 ] µm. Figure 4 shows the variation of β k,j given by Eq. 39 with d k for three fixed values of d j . The difference between assuming Stokes' terminal velocity and using the terminal velocity as given by Eqs (28)-(30) becomes clear for large diameters. Note that β k,j is symmetric in the indices j and k.

280
In the special case that d k = d j Eq. 39 becomes For d j ∼ 1 µm the first term dominates. The second term dominates for all values of d j 10 µm: for d j ∼ 10 µm we get β j,j ∼ 10 −14 m 3 s −1 ; for d j ∼ 100 µm we get β j,j ∼ 10 −11 m 3 s −1 ; for d j ∼ 1000 µm we get β j,j ∼ 10 −8 m 3 s −1 .
In the case that d k d j Eq. 39 becomes Scale analysis (using the values above) shows that the third term can be neglected and the second term is only comparable with the last term when d j ∼ 0.1 µm which is outside the range of interest. Noting that the smallest values of d j , d k ∈ [1, 10 4 ] µm that satisfy d k d j are d j ∼ 10 µm and d k ∼ 1 µm we see that, for all d j 10 µm, the fourth term will dominate and so β k,j is effectively constant (since we are considering fixed d j ). Thus, for d j ∼ 10 µm we get β k,j ∼ 10 −13 m 3 s −1 ; for d j ∼ 100 290 µm we get β k,j ∼ 10 −9 m 3 s −1 ; for d j ∼ 1000 µm we get β k,j ∼ 10 −5 m 3 s −1 . These are consistent with what is observed in  and that this difference increases as d j increases in magnitude. This explains the kink in Figure 4 when d k = d j and why it becomes sharper as d j increases.
For d k d j scale analysis shows that For d j ∈ [1, 10 4 ] µm and d k d j , the second term dominates. Thus, to leading order β k,j ∝ d 4 k for d k d j which is consistent with what is observed in Figure 4 when β k,j is computed with Stokes' terminal velocity or for d k not too large when β k,j is computed with the terminal velocity given by Eqs (28)-(30).
We now turn to the sticking efficiency. On making use of Eq. 33 and Eq. 34 we can write St v as 300 are assumed to be constant and we have also assumed that the two colliding particles have the same density (as in Section 3).
Using the same values of the parameters as above (with constant ρ s ), the constants have the following orders of magnitude: Note that St v and hence α k,j are symmetric in the indices k and j.
Consider first the special case d k = d j . Then Eq. 45 becomes The first term dominates for d j 1 µm; the second term dominates for d j 10 µm.
In the case d k d j Eq. 45 becomes   d k d j that d k ∼ 10 4 µm and so V d 2 k d j 1 is always satisfied. For these values we would expect α < 1. As d j increases the range of d k -values for which V d 2 k d j 1 and α < 1 increases (e.g. if d j ∼ 100 µm then d k 10 µm for V d 2 k d j 1 to hold). Figure 5 shows the variation of α k,j with St cr , q and ρ s computed with the terminal velocity given by Eqs (28) Raising V d 2 k d j > 1 to the power q > 1 will enhance its value whereas raising it to the power q < 1 will diminish its value; similarly for V d 2 j d k . Thus, for example, as shown in Figure 5d when d j = 10 µm and d k d j , V d 2 k d j > 1 for d k 10 −4 m 340 and so, for q > 1, α k,j is smaller than it would be for q = 1 whereas for q < 1 is is larger. These patterns hold true in Figures 5e and f for d k d j though with diminishing values of α k,j for increasing values of d j . Similarly, we see in Figure 5f, for example, that, for d k d j , V d 2 j d k > 1 for all values of d k ∈ [1, 100] µm and so α k,j is closer to unity for q < 1 and vice-versa for q > 1.
The assumption of constant ρ s in the analysis above can be relaxed. We then see that β k,j depends linearly on ρ s except 345 when d k = d j (when β k,j is independent of ρ s ). The sticking efficiency also depends on ρ s : to leading order St v varies like ρ 2 s and so α k,j decreases with increasing ρ s . Figures 5g-i show the variation of α k,j over the range of ρ s -values considered in Section 3. There is more variation of α k,j with ρ s compared with the variation of α k,j with St cr (for fixed q and ρ s ) over the range of St cr -values shown in Figures 5a-c. Figure 6 shows the variation of the collision kernel, K k,j = α k,j β k,j (Eq. 22), with St cr , q and ρ s computed with the terminal 350 velocity given by Eqs (28)-(30). It is immediately clear that while the sticking efficiency tends to be largest for the particles with the smallest diameters this is negated by the relatively small values of the collision rate for particles of the same size.
The net effect is that the largest values of the collision kernel tend to be found for the particles with the largest diameters. The largest range of values occurs for the smallest value of d j and vice-versa. This reflects the dominance of differential settling in the collision kernel.

355
It should be noted that α k,j computed with Stokes' terminal velocity shows a larger variation than that shown in Figure 5.
Since Stokes' terminal velocity is larger than that calculated from Eqs (28)-(30) for large diameters and α k,j is dominated by differential settling, then for large diameters α k,j becomes smaller than the values shown in Figure 5. Conversely Figure 4 shows that β k,j is larger when using Stokes' terminal velocity. The net effect is that the collision kernel computed with Stokes' terminal velocity is similar in magnitude to that shown in Figure 6.  Figure 6 shows that the variation of K k,j with St cr (over the range of values shown in Figures 6a-c) is smaller than the variation of K k,j with ρ s (Figures 6g-i) and that both of these are smaller than the variation of K k,j with q (Figures 6d-f). For a given value of d k , the value of K k,j increases with increasing St cr but decreases with increasing ρ s . Figures 6d-f show that the variation with q is more complicated but the largest values of K k,j occur for the smallest values of q (for a given value of d k ). This explains why the mode of the AGSD in Figure 3 shifts to larger diameters with increasing St cr , decreasing q or ρ s .

365
The behaviour of α k,j , β k,j and its product, K k,j with respect to changes in St cr , q and ρ s explains why there is less variation of the AGSD with St cr compared with q and ρ s . However, this behaviour cannot explain all the variation of the AGSD with ρ s which is much larger than that with either St cr or q. The additional factor is explained by Eq. 21 which shows that the particle number density for a given size bin, N i , increases with decreasing ρ s since m i is proportional to ρ s .

370
We investigate the impact of representing aggregation on dispersion model simulations of the distal ash cloud from the eruption of Eyjafjallajökull volcano in 2010. We consider the period between the 04/05/2010 -08/05/2010, as we have measurements of the GSD and density of the non-aggregated grains sampled from deposits for this time (Bonadonna et al., 2011). The aggregation scheme is coupled to NAME, such that NAME uses the output AGSD at the top of the plume and Mass Eruption Rate (MER) calculated from the buoyant plume scheme initialised with the observed plume heights, at every time step. The aggre-375 gation scheme is initialized with the measured GSD of the non-aggregated ash (see Figure 3d). However, in the initialization of NAME we consider only grains which have a diameter of ≤ 125 µm and take just 5% of the total MER to represent the mass in the distal ash cloud, following the approach of the London VAAC . Figure 7 shows an example of the output AGSD normalized to 125 µm. It is compared to the input size distribution of the non-aggregated grains. Mass has been lost from grains with diameter ≤ 16 µm, however the mode of the normalized distributions are the same, lying at 64 -125 380 µm. The density distribution of the non-aggregated grains is also shown; densities range from 2039 -2738 kg m −3 for this size range (Bonadonna and Phillips, 2003;Bonadonna et al., 2011). Model particles are released with a uniform distribution over the depth of the modelled (bent-over) plume (Devenish, 2013(Devenish, , 2016. The setup of the NAME runs is given in Table 5 and we use the control internal model parameters in the aggregation scheme (Table 3). Figure 8a shows the modelled 1-hour averaged total column mass loadings in the ash cloud at 00:00 on the 05/05/2010, 24 385 hours after the release start, using the measured GSD (normalized to 125 µm) and density distribution of the non-aggregated Eyjafjallajökull ash. In comparison Figure 8b shows the modelled plume using the normalized time-varying AGSD. As the density of the Eyjafjallajökull aggregates is not known the measured density distribution of the single grains is applied. Current regulations in Europe state that airlines must have a safety case accepted in order to operate in ash concentrations greater than 2 x 10 −3 g m −3 . We assume a cloud depth of 1 km and consider the area of the ash cloud with mass loadings > 2 g m −2 390 to compare the differences in the modelled areas which are significant for aircraft operations. Using the aggregated GSD the extent of the ash cloud is only slightly smaller, it is reduced by just ∼3%, reflecting the slight increase in the fraction of larger (aggregated) grains in the ash cloud which have a greater fall velocity and hence shorter residence time in the atmosphere.
However, it is expected that porous aggregates, specifically cored clusters (PC3 , see Brown et al. (2012) and Bagheri et al. (2016) for topology definitions) may have lower densities than single grains of ash of equivalent size (Bagheri et al., 2016;395 Gabellini et al., 2020;Rossi et al., 2021). Figure 9 shows the modelled ash cloud when we assume that the aggregates have densities of 1000 and 500 kg m −3 (Taddeucci et al., 2011;Gabellini et al., 2020;Rossi et al., 2021). As the aggregation scheme does not track explicitly the mass fraction of aggregates versus single grains we must also make an assumption about how much of the mass released is represented by aggregates with the lower density. Here we consider the case where 25, 50 and 75% of the total mass on the ash ≤ 125 µm is represented by aggregates. Assigning a lower density to the aggregates reduces 400 their fall velocity and the extent of the simulated ash cloud increases: if we assume that 75% of the mass of ash ≤ 125 µm is represented by aggregates, when they are assigned a density of 1000 kg m −3 the simulated ash cloud with mass loadings > 2 g m −2 is 32,042 km 2 , this increases to 34,697 km 2 when they are assigned a density of 500 kg m −3 . Figure 10 shows the relative increase in the area of the ash cloud with concentrations > 2 g m −2 as a function of the mass fraction of aggregates in the ash cloud and their density. The circle with a diameter of 1 represents the extent of the modelled cloud when aggregation 405 is not considered (area 29,278 km 2 ). The largest modelled ash cloud is ∼1.2 times bigger. This is achieved when we use the AGSD, assign the aggregates a density of 500 kg m −3 , and assume that aggregates constitute 75% of the total mass released in NAME (ash ≤ 125 µm). We have integrated an aggregation scheme into the atmospheric dispersion model NAME. The scheme is coupled to a one 410 dimensional buoyant plume model and uses the FPT to solve the SCE to simulate aggregation processes in an eruption column.
The time-evolving AGSD at the top of the plume is provided to NAME as part of the source conditions. This represents the first attempt at modelling explicitly the change in the GSD of the ash due to aggregation in a model which is used for operational response, as opposed to assuming a single aggregate class (Cornell et al., 1983;Bonadonna et al., 2002;Costa et al., 2010).
Our scheme predicts that mass is preferentially removed from bins representing the smallest ash (≤ 64 µm). This agrees well 415 with field and laboratory experiments which have also observed that aggregates mainly consist of particles <63 µm in diameter (Bonadonna et al., 2011;James et al., 2002James et al., , 2003. This suggests aggregation will be more prevalent when large quantities of fine ash are generated by the eruption. Previous sensitivity studies of dispersion model simulations of volcanic ash clouds have highlighted the importance of constraining the GSD of ash for operational forecasts, as this parameter strongly influences its residence time in the atmosphere 420 (Scollo et al., 2008;Beckett et al., 2015;Durant, 2015;Poret et al., 2017;Osman et al., 2020;Poulidis and Iguchi, 2020). Here we show that the modelled AGSD is also sensitive to the GSD, and the density, of the non-aggregated ash at the source. When the scheme is initialized with a coarse GSD there are fewer particles per unit volume (lower number concentrations) within the plume and aggregation is reduced. Whereas when particle densities are low, for the same mass flux, there are higher number concentrations and hence more aggregation.

425
Dispersion model simulations are influenced by the interplay between the size and density distributions of the aggregated ash. Aggregates can have higher fall velocities than the smaller single grains of which they are composed, and therefore act to reduce the extent and concentration of ash in the atmosphere (Rossi et al., 2021). However, porous aggregates can also have lower densities than the single grains, and this can act to 'raft' ash to much greater distances (Bagheri et al., 2016;Rossi et al., 2021). In our case study of the Eyjafjallajökull 2010 eruption we found that although mass was lost from bins representing 430 smaller grain sizes the mode of the AGSD did not differ from the source GSD of the erupted non-aggregated ash, for example the output AGSD at 19:00 04/05/2010 has lost mass from ash ≤ 16 µm but the mode remains at 64-125 µm ( Figure 7). As such, we found that using the time-varying AGSD to initialise our dispersion model, rather than the size distribution of the single grains, had little impact on the simulated ash cloud. When we considered that aggregates may have (lower) densities of 1000 and 500 kg m −3 and make up 25 -75% of the total mass of the simulated AGSD we found that the area of the ash cloud 435 with concentrations significant for aircraft operations (> 2 g m −2 ) varied by a factor of just ∼1.2.
Previous studies which have considered the sensitivity of dispersion model forecasts of volcanic ash clouds to the density distribution of the ash have also suggested that simulations are relatively insensitive to this parameter (Scollo et al., 2008;Beckett et al., 2015). In fact, in this case, the modelled ash cloud is more sensitive to the input GSD of the non-aggregated ash at the source, than due to any change to the GSD or density of the ash due to aggregation. Osman et al. (2020) compared NAME 440 simulations initialised with the default GSD used by the London VAAC (which is relatively fine) and the published GSD of ash from the 1991 eruption of Hekla which is much coarser. They found that simulations of the extent of the Eyjafjallajökull 2010 ash cloud with concentrations > 2 g m −2 varied by a factor of ∼2.5.
It should be remembered that operational forecasts are also sensitive to other eruption source parameters needed to initialize dispersion model simulations. Dioguardi et al. (2020) found that given the uncertainty on the MER, forecasts of the area of 445 the Eyjafjallajökull ash cloud with concentrations > 2 x 10 −3 g m −3 varied by a factor of 5. When generating operational forecasts, uncertainty on the plume height, vertical distribution, MER and GSD of the non-aggregated ash at the source could therefore outweigh any error associated with not representing aggregation processes.
In the operational application of atmospheric dispersion models at VAACs it is often assumed that a large fraction of the total erupted mass is deposited close to the source. Only ash with diameters ≤ 100 µm are considered by the London VAAC in their 450 forecasts of the long-range transport of ash clouds . Bonadonna et al. (2011) constrained the GSD of the ash from the 2010 eruption of Eyjafjallajökull volcano using ground sampling and satellite retrievals. They suggest that 41% of the total mass of erupted tephra was represented by grains with a diameter of ≤ 125 µm. Ash of this size can travel significant distances; > 700 km given the plume heights and meteorological conditions during the Eyjafjallajökull eruption . However, measured ash concentrations over the UK, North Sea and North Atlantic made by research aircraft, ground-455 based lidar and satellite retrievals suggest that only ∼5% of the total erupted mass was actually present in the distal ash cloud (Devenish et al., 2012a, b;Dacre et al., 2011). Satellite measurements of the mass loading in ash clouds from El Chichón, Láscar, and Hudson volcanoes (Rose et al., 2000) also suggest that the distal fine ash fraction is just ∼5%. It has been assumed that aggregation, which can enhance the removal of fine ash, might explain the discrepancy between the observed mass on small grains at the source and the relatively small fraction which makes it into the distal ash cloud. Aggregates, up to 600 µm in 460 diameter, composed of grains < 63 µm, were observed in Iceland during the Eyjafjallajökull eruption (Bonadonna et al., 2011).
However, our simulated AGSD still has ∼30% of the total mass on tephra (aggregates and single grains) with diameters ≤ 125 µm, which, given their size and density, would travel further than Iceland before depositing due to sedimentation alone. Below we consider the limitations of our aggregation scheme. Alternatively, other near-source processes which act to prematurely remove ash from the ash cloud, such as gravitational instabilities (Carazzo and Jellinek, 2012;Durant, 2015;Manzella et al., 465 2015), hydrometeor formation , particle-particle interactions (Eychenne et al., 2015;Del Bello et al., 2017) and topography induced perturbations to local wind fields which can enhance sedimentation on the lee-side of mountains (Watt et al., 2015;Eychenne et al., 2017;Poulidis et al., 2017), could be playing a more dominant role. Or, the mass loading of fine ash in the distal ash cloud could be severely under-estimated by the current approaches, as postulated by Cashman and Rust (2019). Further work is needed to constrain the uncertainty on measured mass loadings 470 and GSDs, associated with the methods used for both ground sampling (e.g. Bonadonna et al., 2015) and satellite retrievals (e.g. Stevenson et al., 2015). We would also benefit from further in situ measurements of the GSDs and densities of falling aggregates.

Limitations
To be considerate of computational costs for operational systems we have limited aggregation processes to the eruption column 475 only. However, it is likely that, while ash concentrations remain high, aggregation will continue in the dispersing ash cloud.
As we do not represent electric fields in our scheme we are also unable to explicitly simulate aggregation through electrostatic attraction (Pollastri et al., 2021). Further work is needed to consider this contribution and the implications for the long-range transport of the ash cloud. Our approach may therefore underestimate the amount of aggregation, which could further shift the mode of the aggregated GSD to larger grain sizes. We also disregard disaggregation due to collisions with other aggregates 480 and ash grains (Del Bello et al., 2015;Mueller et al., 2017). This process has received little attention and remains relatively under-constrained, and as such has also been neglected here.
Volcanic plumes are highly turbulent flows characterized by a wide range of interacting length and timescales. The length scale of the largest eddies (the integral scale) is the plume radius (e.g. Cerminara et al., 2016a). Whereas the smallest eddies are at the Kolmogorov scale, the point at which viscosity dominates and the turbulent kinetic energy is dissipated into heat. In the 485 treatment of the collision kernels in our scheme we have assumed that the Saffman Turner limit is satisfied: that the particles are smaller than the smallest turbulent scale and as such are completely coupled with the flow. However, larger particles lie outside this limit and, if sufficiently large, will be uncorrelated with the flow. Further work is needed to consider the treatment of large uncorrelated particles, for example the application of the Abrahamson limit in the treatment of the collision kernels could be explored (Textor and Ernst, 2004). 490 We consider that particle sticking can occur due to viscous dissipation in the surface liquid layer on the ash (Liu et al., 2000;Liu and Litster, 2002). This is based on the assumption that large amounts of water (magmatic, ground water, and atmospheric) will be available, and the assumption that this mechanism will play a dominant role over other possible sticking mechanisms e.g. electrostatic forces . Using scaling analysis (Section 4) we show that the modelled AGSDs are particularly sensitive to the parameters used in the aggregation scheme to control the sticking efficiency of two colliding 495 particles, the critical Stokes number (St cr ) and parameter q (an exponent). Varying these parameters is in some sense equivalent to changing the amount of viscous dissipation acting on the surface of the particles, which is in turn related to the thickness of the surface water layers. Both of these parameters are poorly constrained and would benefit from further calibration with field and laboratory studies. In particular the depth of the liquid layers on ash grains needs to be better understood and applied here.
The sticking efficiency also depends on the relative velocities between the colliding particles. In Eq. 34 we have neglected any 500 effect of the particle inertia induced by the background turbulent flow which represents a further source of uncertainty.
Our 1-dimensional treatment of the SCE does not allow us to represent the change in density of the simulated aggregates, or track explicitly the mass fraction of aggregates versus single grains within a given size bin. Using their aggregation scheme in FPLUME Folch et al. (2016) predicted 10% of the erupted mass from Eyjafjallajökull 2010 was represented by aggregates when column heights ranged from 6 to 7 km (agl) and 20% for column heights between 7.2 and 8.3 km (agl). However, it 505 should be remembered that these values are sensitive to the treatment of ambient conditions within their model (they assume aggregation only occurs in the presence of liquid water) and the fractal exponent they applied to describe the geometry of the aggregates. Our scheme could be significantly improved by using a multi-dimensional description which represents the fluctuation in the density of the growing aggregates and retains information on the mass fraction of aggregated particles.
However, this would also require a better understanding of the structure (porosity) of aggregates. 510 We only consider aggregation due to the presence of water layers on the surface of the particles, following Costa et al. (2010). This approach could neglect the presence of particle clusters (PC), which usually require less water to form, and might therefore underestimate the impact of aggregation: it is suggested that 46% of the particles < 10 µm fell as particle clusters during a field study of the Eyjafjallajökull ash cloud (Bonadonna et al., 2011). We assume that aggregation can occur for RH > 0%, and scale the sticking efficiency by the RH. As the residence time in the presence of liquid water (RH = 1) increases 515 more aggregation occurs, represented by a decrease in mass on smaller grains (m 32 is reduced). If we were to consider that aggregation can only take place when liquid water is present in the modelled eruption column then, for the period considered in this study, aggregation would only occur in the top ∼1 km, and in some instances no liquid water was formed (e.g. 13:00 06/05/2010 and 12:00 07/05/2010, Figure 1). Folch et al. (2010) found, using their 1-dimensional plume model, that there was only a 30 s window for ash to aggregate in the presence of liquid water in the initial phase of the eruption at Mount St Helen's 520 in 1980, which generated a plume which rose 32 km. For less vigorous plumes, like that generated from the eruption of Crater Peak 1992, despite the lower plume heights, they found that there was a longer, albeit still limited, time period (∼ 45 s) for which liquid water was present. Further work is needed to better constrain the influence of the ambient conditions, such as the relative humidity, on liquid bonding of ash aggregates to improve simulations of their formation in volcanic ash clouds.

525
We have integrated an aggregation scheme into the atmospheric dispersion model NAME. The scheme uses an iterative buoyant plume model to simulate the eruption column dynamics and the Smoluchowski Coagulation Equations are solved with a sectional technique which allows us to simulate the AGSD in discrete bins. The modelled AGSD at the top of the eruption column is then used to represent the time-varying source conditions in the dispersion model simulations. Our scheme is based on the assumption that particle sticking is due to viscous dissipation of surface liquid layers on the ash, and scale analysis 530 indicates that our output AGSD is strongly controlled by under-constrained parameters which attempt to represent these liquid layers. The modelled AGSD is also sensitive to the the physical characteristics assigned to the particles in the scheme: the initial GSD and density distribution. Our ability to accurately forecast the long range transport of volcanic ash clouds is, therefore, still limited by real-time information on the physical characteristics of the ash. We found that using the time-evolving AGSD in dispersion model simulations of the Eyjafjallajökull 2010 eruption had very little impact on the modelled extent of the distal 535 ash cloud with mass loadings significant for aviation. However, our scheme does not represent all the possible mechanisms by which ash may aggregate (i.e. electrostatic forces), nor does it distinguish the density of the aggregated grains. Our results indicate the need for more field and laboratory experiments to further constrain the binding mechanisms and composition of aggregates; their size distribution and density.
https://doi.org/10.5281/zenodo.4607421. The NAME code is available under license from the Met Office.