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 (Numerical Atmospheric Dispersion modelling Environment), 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 coagulation 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 GSD to the model parameters, we conducted a simple parametric study and scaling analysis. We find that the modelled aggregated GSD is sensitive to the density distribution and grain size distribution assigned to the non-aggregated particles at the source. 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 assess the impact of using the aggregated GSD on model simulations of the 2010 Eyjafjallajökull ash cloud and consider the implications for operational forecasting. Using the time-evolving aggregated GSD at the top of the eruption column to initialize 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 of equivalent size, the modelled area of the Eyjafjallajökull ash cloud with high concentrations of ash, significant for aviation, is reduced by ∼ 2 %, 24 h after the start of the release. If we assume that the aggregates have a lower density (500 kg m⁻³) than the single grains of which they are composed and make up 75 % of the mass in the ash cloud, the extent is 1.1 times larger.

1 Introduction

In volcanic plumes ash can aggregate, bound by hydro-bonds and electrostatic forces. Aggregates typically have diameters > 63 µm (Brown et al., 2012), and their fall velocity differs from 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. Folch et al., 2010; Mastin et al., 2013, 2016; Beckett et al., 2015).

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 (GSD) by assigning a specific cluster settling velocity to fine ash (Carey and Sigurdsson, 1983) or fitting the distribution used in dispersion models to observations of tephra deposits retrospectively (e.g. Cornell et al., 1983; Bonadonna et al., 2002; Mastin et al., 2013, 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 the Soufrière Hills volcano (Montserrat) were better described by considering the variation in the aggregate size and in the grain size distribution within individual aggregates. Mastin et al. (2016) determined optimal values for the mean and standard deviation of input aggregated GSDs for ash from the eruptions of Mount St Helens, Crater Peak (Mount Spurr), Mount 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 to represent aggregation in the presence of liquid water in an eruption column using the Smoluchowski coagulation equations (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 the 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 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 (Costa et al., 2010; Folch et al., 2010), WRF-Chem (Egan et al., 2020), and an eruption column model, FPLUME (Folch et al., 2016). They all use an approximate solution of the Smoluchowski coagulation equations, which assumes that aggregates can be described by a fractal geometry and particles aggregating onto a single effective aggregate class defined by a prescribed diameter.

Here we introduce an aggregation scheme coupled to a one-dimensional steady-state buoyant plume model, which uses a discrete solution of the Smoluchowski coagulation equations based on the fixed pivot technique (Kumar and Ramkrishna, 1996). As such, we are able to model explicitly the evolution of the aggregated GSD with time in the eruption column. We have integrated our aggregation scheme into the Lagrangian atmospheric dispersion model NAME (Numerical Atmospheric Dispersion modelling Environment). NAME is used operationally by the London Volcanic Ash Advisory Centre (VAAC) to provide real-time forecasts of the expected location and mass loading of ash in the atmosphere (Beckett et al., 2020). In our approach, the aggregated GSD at the top of the plume is supplied to NAME to provide a time-varying estimate of the source conditions. This means that aggregation is considered as being a key process inside the buoyant plume above the vent but neglected in the atmospheric transport. This choice ensures aggregation is represented where ash concentrations are highest (and aggregation most likely), while also respecting the need for reasonable computation times for an operational system. The paper is organized as follows. In Sect. 2, we present the aggregation scheme. In Sect. 3, we perform a parametric study to investigate the sensitivity of the modelled aggregated GSD to the internal model parameters. We show that the modelled size distribution of the aggregates is sensitive to the sticking parameters and the initial erupted GSD and density of the non-aggregated particles. In Sect. 3.1 we present a scale analysis to understand the dependency of the collision kernel on these parameters. In Sect. 4, we assess the impact of using the modelled aggregated GSD on the simulated extent and mass loading of ash in the distal volcanic ash cloud from the eruption of Eyjafjallajökull volcano in 2010 and consider the implications of using an aggregated GSD for operational forecasting. We discuss the results in Sect. 5, before the conclusions are presented in Sect. 6.

2 The aggregation scheme

We use a one-dimensional steady-state integral plume model, where mass, momentum, and total energy are derived for a control volume (Devenish, 2013, 2016). It combines the effects of moisture and the ambient wind and includes the effects of humidity and phase changes of water on the growth of the plume. As aggregation is controlled by the amount of available water, it is essential that we adequately consider the entrainment of 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. The governing equations are given by the following:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle \frac{\text{d}{M}_{\text{z}}}{\text{d}s}}=\left({\mathit{\rho }}_{\text{a}}-{\mathit{\rho }}_{\text{p}}\right)g\mathit{\pi }{b}^{\mathrm{2}}\text{(2)}& {\displaystyle }{\displaystyle \frac{\text{d}{M}_{x,y}}{\text{d}s}}=-Q_{\text{m}}{\displaystyle \frac{\text{d}{U}_{x,y}}{\text{d}s}}\text{(3)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\text{d}H}{\text{d}s}}=& \left(\left(\mathrm{1}-q_{\text{v}}^{\text{a}}\right)c_{\text{pd}}+q_{\text{v}}^{\text{a}}{c}_{\text{pv}}\right)T_{\text{a}}{\displaystyle \frac{\text{d}{Q}_{\text{m}}}{\text{d}s}}-g{Q}_{\text{m}}{\displaystyle \frac{{\mathit{\rho }}_{\text{a}}}{{\mathit{\rho }}_{\text{p}}}}{\displaystyle \frac{w_{\text{p}}}{v_{\text{p}}}}\\ +& [L_j-\mathrm{273}\left(c_{\text{pv}}-c_{\text{p}j}\right)]{\displaystyle \frac{\text{d}{Q}_j}{\text{d}s}}\end{array}\text{(4)}& {\displaystyle }{\displaystyle \frac{\text{d}{Q}_{\text{t}}}{\text{d}s}}=E{q}_{\text{v}}^{\text{a}}\text{(5)}& {\displaystyle }{\displaystyle \frac{\text{d}{Q}_{\text{m}}}{\text{d}s}}=E,\end{array}$$

where s is the distance along the plume axis, j=l or ice, such that the two phases cannot coexist, M_z=Q_mw_p is the vertical momentum flux, $M_{x,y}=\left(u_{\text{p}x,y}-U_{x,y}\right)Q_{\text{m}}$ is the horizontal momentum flux relative to the environment, H=c_ppTQ_m is the enthalpy flux, Q_t=Q_mn_t is the total moisture flux within the plume, and Q_m=ρ_pπb²v_p is the mass flux. The bulk specific heat capacity is given by the following:

$$\begin{array}{}\text{(6)}& c_{\text{pp}}=n_{\text{d}}{c}_{\text{pd}}+n_{\text{v}}{c}_{\text{pv}}+n_j{c}_{\text{p}j}+n_{\text{s}}{c}_{\text{ps}},\end{array}$$

where j=l or ice. The meaning of the symbols used throughout are given in Tables 1 and 2. The entrainment rate depends on the ambient and plume densities, and when the plume is rising buoyantly is given as follows:

$$\begin{array}{}\text{(7)}& E=\mathrm{2}\mathit{\pi }b\sqrt{{\mathit{\rho }}_{\text{a}}{\mathit{\rho }}_{\text{p}}}{u}_{\text{e}},\end{array}$$

where ρ_p is the plume density, as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\text{p}}}}={\displaystyle \frac{n_{\text{g}}}{{\mathit{\rho }}_{\text{g}}}}+{\displaystyle \frac{n_{\text{s}}}{{\mathit{\rho }}_{\text{s}}}}+{\displaystyle \frac{n_j}{{\mathit{\rho }}_j}},\end{array}$$

where j=l or ice, and u_e is the entrainment velocity, as follows:

$$\begin{array}{}\text{(9)}& u_{\text{e}}={\left({\left(k_{\text{s}}\left|\mathrm{\Delta }{\mathit{u}}_{\text{s}}\right|\right)}^f+{\left(k_{\text{n}}\left|\mathrm{\Delta }{\mathit{u}}_{\text{n}}\right|\right)}^f\right)}^{\mathrm{1}/f}.\end{array}$$

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 that k_s is given the symbol α and k_n the symbol β in Devenish, 2013). The radial and cross-flow entrainment terms are raised to an exponent, f, which controls the relative importance of these two terms. Devenish et al. (2010) found that f=1.5 gave the best agreement with large eddy simulations of buoyant plumes in a crosswind and field observations, and we adopt this here.

Ice is produced whenever T<255 K, 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 (r_v) is larger than the saturation mixing ratio (r_s), which is determined using the Clausius–Clapeyron equation as follows:

$$\begin{array}{}\text{(10)}& r_{\text{s}}={\displaystyle \frac{\mathit{\epsilon }{e}_{\text{s}}}{p_{\text{d}}}},\end{array}$$

where ε=0.62 is the ratio of 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 as follows (Emanuel, 1994, p. 117):

$$\begin{array}{}\text{(11)}& e_{\text{s,l}}=\mathrm{6.112}\text{exp}\left({\displaystyle \frac{\mathrm{17.65}(T-\mathrm{273.5})}{T-\mathrm{29.65}}}\right),\end{array}$$

and, for ice, is given by the following (Murphy and Koop, 2005, p. 1558):

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& \text{log}{e}_{\text{s,ice}}=-\mathrm{9.09718}\left({\displaystyle \frac{\mathrm{273.16}}{T}}-\mathrm{1}\right)-\mathrm{3.56654}\\ & \mathrm{\text{log}}\left({\displaystyle \frac{\mathrm{273.16}}{T}}\right)+\mathrm{0.876793}\left(\mathrm{1}-{\displaystyle \frac{T}{\mathrm{273.16}}}\right)\\ & +\text{log}\left(\mathrm{610.71}\right).\end{array}\end{array}$$

The mass fractions of water (n_l) and ice (n_ice) can then be expressed as follows:

$$\begin{array}{}\text{(13)}& {\displaystyle }{n}_{\text{l}}=\text{Max}(\mathrm{0},n_{\text{t},T>\mathrm{255}\text{K}}-n_{\text{d}}{r}_{\text{s,l}})\text{(14)}& {\displaystyle }{n}_{\text{ice}}=\text{Max}(\mathrm{0},n_{\text{t},T<\mathrm{255}\text{K}}-n_{\text{d}}{r}_{\text{s,ice}}),\end{array}$$

where n_t is the total moisture fraction ($n_{\text{t}}=n_{\text{v}}+n_{\text{l, ice}}$), n_{l, ice} is the mass fraction of either liquid water or 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 Smoluchowski coagulation equations are solved using the fixed pivot technique, which transforms a continuous domain of masses (while 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 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 fixed pivot technique 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, as follows:

$$\begin{array}{}\text{(15)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\text{d}{Q}_m}{\text{d}s}}=& {\displaystyle \frac{\text{d}}{\text{d}s}}\left[\left({\mathit{\rho }}_{\text{p}}\mathit{\pi }{b}^{\mathrm{2}}{v}_{\text{p}}\right)n_{\text{d}}\right]+{\displaystyle \frac{\text{d}}{\text{d}s}}\left[\left({\mathit{\rho }}_{\text{p}}\mathit{\pi }{b}^{\mathrm{2}}{v}_{\text{p}}\right)n_{\text{t}}\right]\\ +& {\displaystyle \frac{\text{d}}{\text{d}s}}\left[\left({\mathit{\rho }}_{\text{p}}\mathit{\pi }{b}^{\mathrm{2}}{v}_{\text{p}}\right)n_{\text{s}}\right]=E,\end{array}\end{array}$$

where $n_{\text{d}}+n_{\text{t}}+n_{\text{v}}=\mathrm{1}$. As there is no entrainment or fallout of solids, Eq. (15) can be expressed as follows:

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\text{d}}{\text{d}s}}\left[\left({\mathit{\rho }}_{\text{p}}\mathit{\pi }{b}^{\mathrm{2}}{v}_{\text{p}}\right)n_{\text{d}}\right]+{\displaystyle \frac{\text{d}}{\text{d}s}}\left[\left({\mathit{\rho }}_{\text{p}}\mathit{\pi }{b}^{\mathrm{2}}{v}_{\text{p}}\right)n_{\text{t}}\right]=E\text{(17)}& {\displaystyle \frac{\text{d}}{\text{d}s}}\left[\left({\mathit{\rho }}_{\text{p}}\mathit{\pi }{b}^{\mathrm{2}}{v}_{\text{p}}\right)n_{\text{s}}\right]=\mathrm{0}.\end{array}$$

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 ${\sum }_{i=\mathrm{1}}^{N_{\text{bins}}}{x}_i=\mathrm{1}$. Assuming each size shares an amount of mass flux that is proportional to x_i, Eq. (17) becomes the following:

$$\begin{array}{}\text{(18)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\text{d}}{\text{d}s}}\left[\left({\mathit{\rho }}_{\text{p}}\mathit{\pi }{b}^{\mathrm{2}}{v}_{\text{p}}\right)n_{\text{s}}\right]\sum _{i=\mathrm{1}}^{N_{\text{bins}}}{x}_i=\mathrm{0}\to \\ & \sum _{i=\mathrm{1}}^{N_{\text{bins}}}{\displaystyle \frac{\text{d}}{\text{d}s}}\left[\left({\mathit{\rho }}_{\text{p}}\mathit{\pi }{b}^{\mathrm{2}}{v}_{\text{p}}\right)n_{\text{s}}{x}_i\right]=\mathrm{0},\end{array}\end{array}$$

where we used the linearity of the sum with respect to the derivative operator. This is the continuity equation for the mass flux of solids in the case of a steady-state process. The continuity equation can be seen as being a set of N_bins equations, one for each ith section, where aggregation is then taken into account by introducing source (B_i) and sink (D_i) terms in the right-hand side of Eq. (18). The continuity equation for the ith bin then becomes the following:

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\text{d}}{\text{d}s}}\left[\left({\mathit{\rho }}_{\text{p}}\mathit{\pi }{b}^{\mathrm{2}}{v}_{\text{p}}\right)n_{\text{s}}{x}_i\right]=b^{\mathrm{2}}{m}_i\left[B_i-D_i\right].\end{array}$$

In the fixed pivot technique, the source term B_i states that a given particle of the ith 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 to the ith 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 and the sticking processes of the ith particles with all the other pivots available, and there is no need to redistribute mass. The fixed pivot technique applied to Eq. (19) then becomes the following:

$$\begin{array}{}\text{(20)}& {\displaystyle }\begin{array}{rl}& B_i=\sum _{m_i\le \left(m_k+mj\right)<m_{i+\mathrm{1}}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\delta }}_{\text{kj}}\right)\left({\displaystyle \frac{m_{i+\mathrm{1}}-m_{\text{sum}}}{m_{i+\mathrm{1}}-m_i}}\right)K_{k,j}{N}_k{N}_j\\ & +\sum _{m_{i-\mathrm{1}}\le \left(m_k+mj\right)<m_i}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\delta }}_{\text{kj}}\right)\left({\displaystyle \frac{m_{\text{sum}}-m_{i-\mathrm{1}}}{m_i-m_{i-\mathrm{1}}}}\right)K_{k,j}{N}_k{N}_j\end{array}\end{array}$$

$$\begin{array}{}\text{(21)}& D_i=\sum _{j=\mathrm{1}}^{N_{\text{bins}}}{K}_{i,j}{N}_i{N}_j,\end{array}$$

where N_i is the number of particles of a given mass per unit volume, as follows:

$$\begin{array}{}\text{(22)}& N_i={\displaystyle \frac{{\mathit{\rho }}_{\text{p}}{n}_{\text{s}}{x}_i}{m_i}}.\end{array}$$

K_k,j is the aggregation kernel between particles belonging to bins k and j, respectively, and δ_kj is the Kronecker delta function. As such, the Smoluchowski coagulation equations have been transformed into a set of ordinary differential equations which are solved for each bin representing the ith 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 the particles, as follows:

$$\begin{array}{}\text{(23)}& K_{k,j}={\mathit{\alpha }}_{k,j}{\mathit{\beta }}_{k,j},\end{array}$$

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 (cubic metres per second; hereafter m³ s⁻¹) involved in the collision between particles k and j. We consider the following five different mechanisms (after Pruppacher and Klett, 1996, Costa et al., 2010, and Folch et al., 2016): Brownian motion (${\mathit{\beta }}_{k,j}^{\text{B}}$), interactions due to the differential settling velocities between the particles (${\mathit{\beta }}_{k,j}^{\text{DS}}$), turbulent inertia (${\mathit{\beta }}_{k,j}^{\text{TI}}$), and both laminar (${\mathit{\beta }}_{k,j}^{\text{LS}}$) and turbulent fluid shear ${\mathit{\beta }}_{k,j}^{\text{TS}}$).

$$\begin{array}{}\text{(24)}& {\displaystyle }{\mathit{\beta }}_{k,j}^{\text{B}}={\displaystyle \frac{\mathrm{2}{k}_{\text{B}}T}{\mathrm{3}{\mathit{\mu }}_{\text{a}}}}{\displaystyle \frac{{\left(d_k+d_j\right)}^{\mathrm{2}}}{d_k{d}_j}}\text{(25)}& {\displaystyle }{\mathit{\beta }}_{k,j}^{\text{DS}}={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}{\left(d_k+d_j\right)}^{\mathrm{2}}|V_k-V_j|\text{(26)}& {\displaystyle }{\mathit{\beta }}_{k,j}^{\text{TI}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\displaystyle \frac{\mathit{\pi }{\mathit{ϵ}}^{\mathrm{3}/\mathrm{4}}}{g{\mathit{\nu }}_{\text{a}}^{\mathrm{1}/\mathrm{4}}}}{\left(d_k+d_j\right)}^{\mathrm{2}}|V_k-V_j|\text{(27)}& {\displaystyle }{\mathit{\beta }}_{k,j}^{\text{LS}}={\displaystyle \frac{\mathrm{\Gamma }}{\mathrm{6}}}{\left(d_k+d_j\right)}^{\mathrm{3}}\text{(28)}& {\displaystyle }{\mathit{\beta }}_{k,j}^{\text{TS}}={\left({\displaystyle \frac{\mathrm{1.7}\mathit{ϵ}}{{\mathit{\nu }}_{\text{a}}}}\right)}^{\mathrm{1}/\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{\mathrm{8}}}{\left(d_k+d_j\right)}^{\mathrm{3}},\end{array}$$

where d_k and d_j are the diameters, and V_k and V_j are the sedimentation velocities of the colliding particles, as follows:

$$\begin{array}{}\text{(29)}& V_{k,j}=\sqrt{\left({\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}{\displaystyle \frac{d_{k,j}}{C_{\text{d}}}}g{\displaystyle \frac{{\mathit{\rho }}_{\text{s}}-{\mathit{\rho }}_{\text{a}}}{{\mathit{\rho }}_{\text{a}}}}\right)},\end{array}$$

where C_d is the drag coefficient, and Re is the Reynolds number, as follows:

$$\begin{array}{}\text{(30)}& C_{\text{d}}={\displaystyle \frac{\mathrm{24}}{Re}}\left(\mathrm{1}+\mathrm{0.15}R{e}^{\mathrm{0.687}}\right)\end{array}$$

$$\begin{array}{}\text{(31)}& Re={\displaystyle \frac{V_{k,j}{d}_{k,j}}{{\mathit{\nu }}_{\text{a}}}},\end{array}$$

and the sedimentation velocity is evaluated using an iterative scheme following Arastoopour et al. (1982). The laminar fluid shear is taken to be $\mathrm{\Gamma }=|\mathrm{d}{w}_{\text{p}}/\mathrm{d}z|$. 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, as follows (Textor and Ernst, 2004):

$$\begin{array}{}\text{(32)}& \mathit{ϵ}={\displaystyle \frac{(\mathrm{0.1}{v}_{\text{p}}{)}^{\mathrm{3}}}{b}}.\end{array}$$

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 laminar shear and turbulent shear kernels):

$$\begin{array}{}\text{(33)}& {\mathit{\beta }}_{k,j}={\mathit{\beta }}_{k,j}^{\text{B}}+\text{Max}({\mathit{\beta }}_{k,j}^{\text{LS}},{\mathit{\beta }}_{k,j}^{\text{TS}})+{\mathit{\beta }}_{k,j}^{\text{TI}}+{\mathit{\beta }}_{k,j}^{\text{DS}}.\end{array}$$

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). 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 parameterized in terms of the collision Stokes number (St_v) as follows:

$$\begin{array}{}\text{(34)}& S{t}_{\text{v}}={\displaystyle \frac{\mathrm{8}\widehat{\mathit{\rho }}{U}_r}{\mathrm{9}{\mathit{\mu }}_{\text{l}}}}{\displaystyle \frac{d_k{d}_j}{d_k+d_j}},\end{array}$$

which is a function of the average density of the two colliding particles ($\widehat{\mathit{\rho }}$), the liquid viscosity (μ_l), and the relative velocities between the colliding particles (U_r), here approximated as follows:

$$\begin{array}{}\text{(35)}& {\displaystyle }{U}_r={\displaystyle \frac{\mathrm{8}{k}_{\text{B}}T}{\mathrm{3}\mathit{\pi }{\mathit{\mu }}_{\text{a}}{d}_k{d}_j}}+|V_k-V_j|+{\displaystyle \frac{\mathrm{4}}{\mathit{\pi }}}{\mathrm{\Gamma }}_{\text{max}}(d_k+d_j)\text{(36)}& {\displaystyle }{\mathrm{\Gamma }}_{\text{max}}=\text{max}\left({\displaystyle \frac{\mathrm{\Gamma }}{\mathrm{6}}},{\displaystyle \frac{\mathrm{1}}{\mathrm{8}}}{\left({\displaystyle \frac{\mathrm{1.7}\mathit{ϵ}}{{\mathit{\nu }}_{\text{a}}}}\right)}^{\mathrm{1}/\mathrm{2}}\right).\end{array}$$

Following a collision, particles stick together if the relative kinetic energy of the colliding particles is completely depleted by viscous dissipation in the surface liquid layer on the particles (Liu et al., 2000). The condition for this to occur is given by the following:

$$\begin{array}{}\text{(37)}& S{t}_{\text{v}}<S{t}_{\text{cr}}=\ln \left({\displaystyle \frac{h}{h_{\text{a}}}}\right),\end{array}$$

where h is the thickness of the liquid layer, and h_a is the surface asperity or surface roughness (Liu et al., 2000; Liu and Lister, 2002). Unfortunately, this information is poorly constrained for volcanic ash. Instead, Costa et al. (2010) propose the following parameterization for the sticking efficiency:

$$\begin{array}{}\text{(38)}& {\mathit{\alpha }}_{k,j}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\left(\frac{S{t}_{\text{v}}}{S{t}_{\text{cr}}}\right)}^q}},\end{array}$$

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 Fig. 12 in Gilbert and Lane, 1994, and Fig. 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 three-dimensional turbulent flow simply represents an average of the flow conditions 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. Furthermore, 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 sticking efficiency (α_k,j) by the relative humidity as follows:

$$\begin{array}{}\text{(39)}& {\mathit{\alpha }}_{k,j}={\mathit{\alpha }}_{k,j}\cdot RH.\end{array}$$

In the presence of ice, we assume that the sticking efficiency is constant, and ${\mathit{\alpha }}_{k,j}=\mathrm{0.09}$, following Costa et al. (2010) and Field et al. (2006).

Table 1 List of Latin symbols. Quantities with a superscript of 0 indicate values at the source.

Table 2 List of Greek symbols.

3 Aggregation model sensitivities

To consider the influence of uncertainty on the source and internal model parameters on the simulated aggregated GSD, 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.

For our case study, we consider the 2010 eruption of the Eyjafjallajökull volcano, Iceland (location 63.63^∘ lat, −19.62^∘ long; summit height 1666 m a.s.l. – above sea level), between 4 and 8 May 2010. We use the time profile of plume heights given in Webster et al. (2012), which are based on radar data, pilot reports, and Icelandic coastguard observations. Meteorological data, used by the aggregation scheme and NAME simulations, are from the global configuration of the Unified Model (UM) which, for this period, had a horizontal resolution of ∼ 25 km (at mid-latitudes) and a temporal resolution of 3 h. Figure 1 shows the relative humidity (RH), temperature (T), and mixing ratios of liquid water (n_l/n_d), water vapour (n_v/n_d), and ice (n_i/n_d) 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 04 May 2010, the maximum observed plume height is 7000 m a.s.l., liquid water starts to form at 4122 m a.s.l., but no ice forms in the plume. At 12:00 UTC on 5 May 2010, the observed maximum plume height is lower, reaching just 5500 m a.s.l., liquid water is present from 3684 m a.s.l., and again no ice is formed. However, at 13:00 UTC on 6 May 2010, no liquid water forms in the plume, only ice, which is present from 5582 m a.s.l., and the maximum observed plume height is 10 000 m a.s.l. At 12:00 UTC on 7 May 2010, the maximum observed plume height is 5500 m a.s.l., and there is neither liquid water nor ice in the plume – only water vapour.

The control values used for the source and internal model parameters in the aggregation scheme are given in Table 3, along with the range of values for each parameter considered in the sensitivity study. Values are based on the existing literature, and the sources used are also listed in Table 3. The scheme is initialized with a GSD with a uniform distribution of mass across 14 bins, representing particles 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 millimetres (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 (Fig. 2). The output-aggregated GSD at the top of the plume, defined as the point at which W_p=0, is assessed. Given the nature of the Smoluchowski coagulation equations, the aggregation scheme does not track explicitly the mass fraction of aggregates versus single grains within a given bin size. Instead, we consider how the mode of the output-aggregated GSD varies and compare the mass fraction on particles with diameter ≤ 32 µm (m₃₂), which predominantly lose mass to larger aggregates, for each sensitivity run. First, we consider how the aggregated GSD varies as conditions within the plume change over time given the local meteorological and eruption conditions (plume height). Figure 2 shows the output aggregated GSDs, for the same times as the plume conditions shown in Fig. 1, compared to the input GSD. We find that, in all the cases considered, the mode of the aggregated GSD is always the same, with most of the mass now residing in the 125–250 µm bin. When particles spend more time in the presence of liquid water, m₃₂ decreases slightly; m₃₂=32 % at 19:00 UTC on 4 May 2010 when liquid water is present from 4122 m a.s.l., but m₃₂=33 % at 12:00 UTC on 5 May 2010 when liquid water is only available over a more limited depth (Fig. 2a, b). Aggregation still occurs when there is only ice present and no liquid water (6 May 2010 at 13:00 UTC; Fig. 2c) and when there is no ice or liquid water present (12:00 UTC on 7 May 2010; Fig. 2d).

Figure 1 Mixing ratios of water vapour (n_v/n_d), liquid water (n_l/n_d), and ice (n_ice/n_d) with height along the buoyant plume axis, for the eruption of Eyjafjallajökull volcano between 4 and 7 May 2010. Also shown are the relative humidity and temperature. Note the variation in axis scales. The maximum height of the modelled plume axis, when the plume is bent over as in this case, is the maximum observed plume height (provided in the figure titles) minus the plume radius.

Figure 2 Modelled aggregated GSDs corresponding to the times and phase conditions shown in Fig. 1. The aggregation scheme is initialized with a GSD with a uniform distribution of mass, as indicated by the dark blue bars.

Woodhouse et al. (2016)Aubry et al. (2017)Costa et al. (2016)Woodhouse et al. (2016)Devenish (2013)Woods (1988)Devenish (2013)Costa et al. (2016)Costa et al. (2016)Costa et al. (2010)Gilbert and Lane (1994)Costa et al. (2010)Gilbert and Lane (1994)Bonadonna and Phillips (2003)Bonadonna et al. (2011)Gudnason et al. (2017)

Table 3 Model variables used in the aggregation scheme to represent the eruption conditions. The control values listed for each parameter are based on the defaults used in the existing literature. The range of parameter values considered in the sensitivity study are also given.

Table 4 Properties of the simulated aggregated GSD from the model sensitivity runs. The output is for 19:00 UTC on 4 May 2010. Using control values (Table 3), the mode is at 125–250 µm and m₃₂ 32 %.

The mode and m₃₂ of the simulated aggregated GSD for each sensitivity run output at 19:00 UTC on 4 May 2010 are listed in Table 4. Using the control parameters, the mode of the aggregated GSD lies at 125–250 µm, and m₃₂ is 32 % at this time (cf. Fig. 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, $\widehat{\mathit{\rho }}={\mathit{\rho }}_{\text{s}}$). Figure 3 shows how the cumulative distribution of the aggregated GSD 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 understood and, therefore, under-constrained. Figure 3a–b show the aggregated GSD when St_cr and q are varied by a factor of 2. When q=0.4, m₃₂ is 30 %, and the mode of the aggregated GSD moves to 500–1000 µm; when q=1.6, the mode lies at 64–125 µm, and for St_cr in the range 0.65–2.6, m₃₂ varies from 31 %–33 %. When particles have a (low) density of 500 kg m⁻³, m₃₂ is 24 %, and the modal bin is 500–1000 µm (Fig. 3c). We find that, when using a relatively coarse input GSD (from the eruption of Hekla in 1991), there is very little aggregation, and there is no change in m₃₂ or the modal grain size from the input GSD (Fig. 3d). Whereas, when using the Eyjafjallajökull (2010) GSD, which is much finer, the mode of the aggregated GSD is shifted to larger sizes. Output from the sensitivity runs for other times during the eruption (corresponding to those in Fig. 2) are provided in the Supplement and show the same behaviour (Figs. S1–S3).

The aggregated GSD shows little sensitivity to the model values assigned to define the plume conditions within the ranges investigated, namely the entrainment parameters (k_s and k_n), the initial mass fraction of dry gas and water vapour ($n_{\text{d}}^{\mathrm{0}}$,$n_{\text{v}}^{\mathrm{0}}$), the plume temperature at the source (T₀), or the source mass flux (Q_m; see Table 4 and Figs. 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₃₂ varies by just 1 %.

Figure 3 Sensitivity of the output aggregated GSD to the sticking efficiency parameters, (a) St_cr and (b) q, and the physical characteristics assigned to the particles, (c) particle density ρ_s, and (d) input GSD. Output is for 19:00 UTC on 4 May 2010, and the plume height is 7000 m a.s.l. (cf. Fig. 2a). Note that the blue lines represent simulations using the control values, and AGSD is the aggregated GSD.

3.1 Scale analysis of the collision kernel

To gain insight into the dependence of the aggregation kernel (K_k,j) on the critical Stokes number (St_cr), the parameter (q), the size of the particles (d_k,j), and their density (ρ_s), we performed a scale analysis of the collision rate (β_k,j) and sticking efficiency (α_k,j), the details of which are provided in Appendix A. All of the parameters, other than the one being varied, are kept fixed at the control (default) values listed in Table 3.

First, we consider how the behaviour of the collision rate (β_k,j) changes as the particle size and density varies. Figure 4 shows the variation in the collision rate between two particles of the same fixed density, where the diameter of one of the colliding particles (d_j) is kept fixed (diameters of 10, 100, and 1000 µm are considered) and the diameter of the second particle (d_k) is allowed to vary between 1 and 10 000 µm, consistent with the GSD of the tephra considered in this study (Table 4). For typical values of each of the parameters that occur in the different kernels considered in the collision rate equation (Eq. 33), and assuming Stokes drag, the scale analysis in Appendix A shows that the collision rate is dominated by differential settling (${\mathit{\beta }}_{k,j}^{\text{DS}}$) when d_k≫d_j or d_k≪d_j. As it is the larger particle which determines the collision rate in these limiting cases, then for d_k≪d_j the collision rate is effectively constant; the scale analysis gives the correct order of magnitude for β_k,j in these cases. When d_k≫d_j, then, to leading order, the collision rate increases to the fourth power of the diameter of the colliding particle (${\mathit{\beta }}_{k,j}\propto d_k^{\mathrm{4}}$) and is independent of d_j. However, Fig. 4 shows that when d_k≳10² µm, β_k,j departs from this power law when the Reynolds-number-dependent terminal velocity is used (Eqs. 29–31). When d_k=d_j, then the collision rate is dominated by shear, except for the very smallest particles (of the order of 1 µm) when it is dominated by Brownian motion. The scale analysis gives the correct order of magnitude for these cases and explains the kinks seen in Fig. 4, when d_k=d_j, and why they become sharper as d_j increases. When the assumption of constant ρ_s is relaxed, it is easily seen that β_k,j depends linearly on ρ_s (through ${\mathit{\beta }}_{k,j}^{\text{TI}}$ and ${\mathit{\beta }}_{k,j}^{\text{DS}}$) for d_k≠d_j; when d_k=d_j, the collision rate is independent of ρ_s.

Figure 4 The variation in the collision rate (β_k,j) between particles as the size of the particle d_k varies for three values of d_j. The solid lines are calculated with Stokes terminal velocity, and the dashed lines are calculated with the terminal velocity given by Eqs. (29)–(31). The black line is proportional to $d_k^{\mathrm{4}}$.

We now turn to the sensitivity of the sticking efficiency (α_k,j) to the critical Stokes number (St_cr), the sticking parameter (q), and the density of the particles (ρ_s). It follows immediately from Eq. (37) that increasing St_cr increases the range of values of the collision Stokes number (St_v) for which coalescence can occur. When St_v≪St_cr, ${\mathit{\alpha }}_{k,j}\approx \mathrm{1}$, and coalescence is almost certain, q has the effect of enhancing or reducing the effect of $S{t}_{\text{v}}/S{t}_{\text{cr}}$, with q>1 reducing the effect in this limit and so increasing the sticking efficiency further and vice versa for q<1. When St_v≫St_cr, ${\mathit{\alpha }}_{k,j}\ll \mathrm{1}$, and there is effectively no coalescence; q>1 has the effect of increasing the value of $S{t}_{\text{v}}/S{t}_{\text{cr}}$ relative to its value with q=1 and, hence, reducing the sticking efficiency still further, whereas the converse applies when q<1. How the sticking efficiency depends on diameter is determined by the dependence of St_v on d_j and d_k, and this is given by the scale analysis in Appendix A. When d_k=d_j, St_v (via U_r, as given by Eq. 35) is dominated by shear ($S{t}_{\text{v}}\propto d_j^{\mathrm{2}}$), except for the smallest particles (of the order of 1 µm) when it is dominated by Brownian motion ($S{t}_{\text{v}}\propto \mathrm{1}/d_j$). When d_k≠d_j, St_v is dominated by differential settling; for d_k≪d_j, we find that $S{t}_{\text{v}}\propto d_j^{\mathrm{2}}{d}_k$, whereas, for d_k≫d_j, we have $S{t}_{\text{v}}\propto d_k^{\mathrm{2}}{d}_j$. It is the size of this term relative to St_cr which determines whether α_k,j is close to one or not.

Figure 5 shows the sensitivity of the sticking efficiency (α_k,j) to the critical Stokes number (St_cr), the sticking parameter (q), and the density of the particles (ρ_s) for three fixed values of d_j (10, 100, and 1000 µm). Figure 5 clearly shows the asymmetry in the dependence of α_k,j on d_j and d_k, as highlighted above. Figure 5a–c show the sensitivity of α_k,j to the critical Stokes number (St_cr) as follows: as St_cr increases, α_k,j increases towards one for fixed d_j and d_k indicating, as expected, a greater propensity for the particles to coalesce. As d_j increases, α_k,j tends to decrease, indicating the increasing importance of the ratio $S{t}_{\text{v}}/S{t}_{\text{cr}}$ in the evaluation of α_k,j. Figure 5d–f show the sensitivity of α_k,j to the parameter q which acts to alter the shape of the sticking matrix. There is more variation with q than with St_cr; because q appears in α_k,j as an exponent, a change in the value of q is not simply a multiplicative change as it is with a change in the value of St_cr. Figure 5g–i show the sensitivity of α_k,j to ρ_s. This occurs because to leading order $S{t}_{\text{v}}\propto {\mathit{\rho }}_{\text{s}}^{\mathrm{2}}$ and so α_k,j decreases with increasing ρ_s.

Figure 6 shows the variation in the collision kernel, $K_{k,j}={\mathit{\alpha }}_{k,j}{\mathit{\beta }}_{k,j}$ (Eq. 23), with St_cr, q, and ρ_s. It is immediately clear that, while the sticking efficiency tends to increase when particles are small (Fig. 5), this effect is negated by the reduction in the collision rate (Fig. 4) 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 (Fig. 6). 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.

The collision kernel inherits its sensitivity to St_cr and q from the sticking efficiency, α_k,j. Figure 6 shows that the variation in K_k,j with St_cr (over the range of values shown in Fig. 6a–c) is smaller than the variation in K_k,j with ρ_s (Fig. 6g–i) and that both of these are smaller than the variation in K_k,j with q (Fig. 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. Figure 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 aggregated GSD in Fig. 3 shifts to larger diameters with increasing St_cr, decreasing q or ρ_s (see also Table 4). The behaviour of the sticking efficiency (α_k,j), collision rate (β_k,j) and its product, the collision kernel (K_k,j), with respect to changes in St_cr, q, and ρ_s explains why there is less variation in the aggregated GSD with St_cr compared with q and ρ_s. However, this behaviour cannot explain all the variation in the aggregated GSD with ρ_s, which is much larger than that with either St_cr or q. The additional factor is explained by the fact that the particle number density for a given size bin, N_i, increases with decreasing ρ_s, since m_i is proportional to ρ_s (Eq. 22).

Figure 5 The variation in sticking efficiency (α_k,j), with St_cr (a–c), q (d–f), and ρ_s (g–i) for three fixed values of d_j, i.e. d_j=10 µm (a, d, g), d_j=100 µm (b, e, h), and d_j=1000 µm (c, f, i). The terminal velocity is calculated using Eqs. (29)–(31). The diameter d_k0=1 µm.

Figure 6 The variation in the collision kernel (K_k,j) with St_cr (a–c), q (d–f), and ρ_s (g–i) for three fixed values of d_j, i.e. d_j=10 µm (a, d, g), d_j=100 µm (b, e, h), and d_j=1000 µm (c, f, i). The terminal velocity is calculated using Eqs. (29)–(31). The diameter d_k0=1 µm.

4 Dispersion modelling: a case study of the eruption of Eyjafjallajökull in 2010

We now investigate the impact of representing aggregation on dispersion model simulations of the distal ash cloud from the eruption of the Eyjafjallajökull volcano in 2010. We consider the period between 4 and 8 May 2010, as we have measurements of the GSD and density of the non-aggregated grains for this time (Bonadonna et al., 2011). The aggregation scheme is initialized with the measured GSD of the non-aggregated tephra, with diameters between 1 µm and 8 mm, provided in Bonadonna et al. (2011), which is based on both deposit and satellite measurements. Figure 7 shows the output aggregated GSD at 19:00 UTC on 4 May 2010. Following aggregation, there are fewer grains with diameters ≤16 µm, but the mode of the distribution remains at 500–1000 µm. Furthermore, the total fraction of mass on ash with diameters ≤125 µm has changed very little; prior to aggregation, 41 % of the total mass is represented by ash with diameters ≤125 µm. This is reduced to just 39 % following aggregation. The density distribution of the non-aggregated Eyjafjallajökull grains is also shown; densities range from 2738 to 990 kg m⁻³ for this size range (Bonadonna and Phillips, 2003; Bonadonna et al., 2011). The aggregation scheme is coupled to NAME, such that it uses the output-aggregated GSD at the top of the plume and the mass eruption rate (MER) calculated from the buoyant plume scheme, which is initialized with the observed plume heights, at every time step. When NAME is used by the London VAAC, it is initialized with small particles which are expected to remain in the atmosphere and contribute most to the distal ash cloud (Beckett et al., 2020; Osman et al., 2020). We follow this approach in our simulations; we use the aggregated GSD up to 125 µm, and the MER is scaled to represent the mass on these grains only. For example, at 19:00 UTC on 4 May 2010, 39 % of the total mass erupted is released over the seven bins representing ash with diameters ≤125 µm (Fig. 7). The exact diameter of each model particle is allocated such that the log of the diameter is uniformly distributed within each size bin. These model particles are then released with a uniform distribution over the depth of the modelled (bent-over) plume. See Devenish (2013, 2016) for details of how the plume radius (depth) is constrained. The set-up 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 h averaged total column mass loadings in the ash cloud at 00:00 UTC on 5 May 2010, 24 h after the release start, using the measured GSD and density distribution of the non-aggregated Eyjafjallajökull particles. In comparison, Fig. 8b shows the modelled plume using the time-varying aggregated GSD. 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 to operate in ash concentrations greater than $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ g m⁻³. We assume a cloud depth of 1 km and consider the area of the ash cloud with mass loadings > 2 g m⁻² 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 ∼ 2 %, 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 which consist of a large core particle (> 90 µm) covered by a thick shell of smaller particles (Brown et al., 2012; Bagheri et al., 2016) may have lower densities than single grains of ash of equivalent size (Bagheri et al., 2016; 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⁻³ (Taddeucci et al., 2011; Gabellini et al., 2020; Rossi et al., 2021). As the aggregation scheme does not track explicitly the mass fraction represented by aggregates versus single grains for a given size bin, we must also make an assumption about how much of the mass released is represented by aggregates with lower density. Here we consider the case where 25 %, 50 %, and 75 % of the mass on each size bin, for ash with diameters ≤125 µm, is represented by aggregates. Assigning a lower density to the aggregates reduces 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, then, when they are assigned a density of 1000 kg m⁻³, the simulated ash cloud with mass loadings >2 g m⁻² is 152 308 km². This increases to 159 966 km² when they are assigned a density of 500 kg m⁻³. Figure 10 shows the relative increase in the area of the ash cloud with concentrations >2 g m⁻² 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 is not considered (area 143 459 km²). The largest modelled ash cloud is ∼ 1.1 times bigger. This is achieved when we use the aggregated GSD, assign the aggregates a density of 500 kg m⁻³, and assume that aggregates constitute 75 % of the total mass released in NAME (ash ≤125 µm).

Figure 7 The GSD of the Eyjafjallajökull (2010) non-aggregated tephra (dark grey bars; from Bonadonna et al., 2011), used to initialize the aggregation scheme, and the modelled aggregated GSD at the top of the plume (light grey bars), at 19:00 UTC on 4 May 2010. The density distribution of the non-aggregated particles, taken from Bonadonna et al. (2011), is also shown.

Table 5 Input parameters for the NAME runs.

Figure 8 Modelled 1 h averaged total column mass loadings of the Eyjafjallajökull ash cloud at 00:00 UTC on 5 May 2010, using (a) the measured GSD of the non-aggregated ash and (b) the time-varying aggregated GSD. The measured density distribution of the non-aggregated ash grains is applied in both cases. The area of the ash cloud with mass loadings >2 g m⁻², which is significant for aircraft operations, is shown.

Figure 9 Modelled 1 h averaged total column mass loadings of the Eyjafjallajökull ash cloud at 00:00 UTC on 5 May 2010 when 25 %, 50 %, and 75 % of the mass is on aggregates with densities of 1000 and 500 kg m⁻³. The area of the ash cloud with mass loadings >2 g m⁻², which is significant for aircraft operations, is shown.

Figure 10 Relative areas of the Eyjafjallajökull ash cloud with concentrations >2 g m⁻² at 00:00 UTC on 5 May 2010. The area of the ash cloud, when aggregation is not considered, has a relative radius of 1. The modelled areas, using aggregated GSDs when 25 %, 50 %, and 75 % of the mass released is assumed to be on aggregates with a density of 1000 and 500 kg m⁻³, are compared.

5 Discussion

We have integrated an aggregation scheme into the atmospheric dispersion model NAME. The scheme is coupled to a one-dimensional buoyant plume model and uses the fixed pivot technique to solve the Smoluchowski coagulation equations to simulate aggregation processes in an eruption column. The time-evolving aggregated GSD 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 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., 2002, 2003). This suggests that 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 (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 aggregated GSD is also sensitive to the GSD, and the density, of the non-aggregated particles 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. When particle densities are low, for the same mass flux, there are higher number concentrations and hence more aggregation.

Dispersion model simulations are influenced by the interplay between the size and density distributions assigned to the particles. 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 smaller grain sizes, the mode of the aggregated GSD did not differ from the source GSD of the erupted non-aggregated particles; for example, the output aggregated GSD at 19:00 UTC on 4 May 2010 has lost mass from ash ≤16 µm but the mode remains at 64–125 µm (Fig. 7).

Our dispersion model set-up in this study reflects the choices used by the London VAAC; as such we examine the transport and dispersion of ash with diameters ≤125 µm, and we consider the implications for the modelled extent of the ash cloud with mass loadings of significance to the aviation industry. We found that using the time-varying aggregated GSD to initialize 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⁻³ and make up 25 %–75 % of the total mass of the simulated aggregated GSD, we found that the area of the ash cloud with concentrations significant for aircraft operations (>2 g m⁻²) varied by a factor of just ∼ 1.1. 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 particles at the source than due to any change to the GSD or density due to aggregation. Osman et al. (2020) compared NAME simulations initialized 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 ash cloud in 2010 with concentrations >2 g m⁻² 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 the Eyjafjallajökull ash cloud with concentrations $>\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ g m⁻³ varied by a factor of 5. When generating operational forecasts, uncertainty on the plume height, vertical distribution, MER, and GSD of the non-aggregated particles at the source could, therefore, outweigh any error associated with not representing aggregation processes.

The grain size distribution of the non-aggregated Eyjafjallajökull tephra, determined using ground sampling and satellite retrievals, indicated that ∼ 20 wt % of the total mass erupted was ash with diameters ≤16 µm (Fig. 7; Bonadonna et al., 2011). Given their relatively low fall velocities, ash of this size can travel significant distances, e.g. >3000 km, given the plume heights and meteorological conditions during the Eyjafjallajökull eruption (Beckett et al., 2015). However, Bonadonna et al. (2011) observed that ∼ 50 % of this very fine ash was deposited on land in Iceland (within 60 km from the vent) and, therefore, must have fallen out quicker than their settling velocity would allow, due to either particle aggregation or gravitational instabilities or both. Our simulated aggregated GSDs at the top of the eruption column have formed very few larger grains; as such, 18 % of the total mass erupted is still represented by grains with diameters ≤16 µm. The small mass reduction in the fraction of grains with diameters ≤16 µm predicted by our model could be due to various limitations in our scheme and approach, for example, a poor description of dry aggregation, of regions of water saturation, and of particle collision in turbulent flows. We now discuss these limitations in detail.

5.1 Limitations

To be considerate of the computational costs for operational systems, we have limited aggregation processes to the eruption column 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 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, i.e. the point at which viscosity dominates, and the turbulent kinetic energy is dissipated into heat. In the treatment of the collision kernels in our scheme, we have assumed that the Saffman–Turner limit is satisfied and 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).

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 Lister, 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 (Costa et al., 2010). Furthermore, this approach could neglect the presence of particle clusters (Brown et al., 2012), which usually require less water to form, and so our approach might also be underestimating the formation of these aggregates.

Using scaling analysis (Sect. 3.1), we show that the modelled aggregated GSDs are particularly sensitive to the parameters used in the aggregation scheme to control the sticking efficiency of two colliding particles, i.e. 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 our aggregation scheme 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. (35), we have neglected any effect of the particle inertia induced by the background turbulent flow which represents a further source of uncertainty.

For the eruption considered in this study, liquid water is only present in top ∼ 1 km of the plume, and in some instances, no liquid water was formed (e.g. 13:00 UTC on 6 May 2010 and 12:00 UTC on 7 May 2010; Fig. 1). Folch et al. (2010) found, using their one-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 Helens in 1980, which generated a plume which rose 32 km, and only a ∼ 45 s window during the less vigorous eruption of Crater Peak in 1992. Atmospheric conditions in the tropics can generate taller eruption plumes, which entrain more water, than these eruptions in drier environments and, as such, may promote conditions more ideal for aggregation (Tupper et al., 2009).

Our one-dimensional treatment of the plume does not fully represent the three-dimensional turbulent flow and may be missing local pockets of liquid water. Initial experimental studies also suggest that aggregation can occur at relatively low humidity (Telling and Dufek, 2012; Telling et al., 2013; Mueller et al., 2016). As such, in our approach, we allow sticking to occur in regions where the relative humidity is <100 % and liquid water is not yet present. Experimental data which better constrain the influence of the ambient conditions, such as the relative humidity, on liquid bonding of ash aggregates could improve our simulations of aggregate formation in volcanic ash clouds.

When liquid water and ice did form, mass mixing ratios suggest that our modelled plumes are liquid water/ice rich; the maximum mass mixing ratio of liquid water (at the top of the plume) was $\mathrm{8.3}\times {\mathrm{10}}^{-\mathrm{4}}$ kg kg⁻¹ at 19:00 UTC on 4 May 2010, and the maximum mass mixing ratio of ice was $\mathrm{3.3}\times {\mathrm{10}}^{-\mathrm{4}}$ kg kg⁻¹ at 13:00 UTC on 6 May 2010. In comparison, mid-level mixed-phase clouds typically have liquid water mixing ratios of $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{4}}$–$\mathrm{4}\times {\mathrm{10}}^{-\mathrm{4}}$ kg kg⁻¹ and ice mixing ratios of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}}$–$\mathrm{4}\times {\mathrm{10}}^{-\mathrm{5}}$ kg kg⁻¹ (Smith et al., 2009). Atmospheric conditions in the tropics would likely ensure even higher quantities of ice in volcanic plumes (Tupper et al., 2009). Our scheme does not account for interactions between the hydrometeors formed and the ash particles; as such, we can neither represent the role of ash as an effective ice-nucleating particle (Durant et al., 2008; Gibbs et al., 2015), nor can we account for the process of ash-laden hailstones acting to preferentially remove fine ash from the atmosphere (Van Eaton et al., 2015).

Fine ash could also be preferentially removed from both the plume and dispersing ash cloud due to other size-selective processes currently not described in NAME, such as gravitational instabilities, which represent a dominant process for this eruption (Durant, 2015; Manzella et al., 2015). Ash aggregation might be also enhanced by the formation of fingers as a result of gravitational instabilities due to an increase in both ash concentration and turbulence (e.g. Carazzo and Jellinek, 2012; Scollo et al., 2017).

Finally, our one-dimensional treatment of the Smoluchowski coagulation equations 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 bin size. 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. To implement this change effectively would also require a better understanding of the structure (porosity) of aggregates.

6 Conclusions

We have integrated an aggregation scheme into the atmospheric dispersion model NAME. The scheme uses a 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 aggregated GSD in discrete bins. The modelled aggregated GSD 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 the viscous dissipation of surface liquid layers on the ash, and scale analysis indicates that our output-aggregated GSD is strongly controlled by under-constrained parameters which attempt to represent these liquid layers. The modelled aggregated GSD is also sensitive to the physical characteristics assigned to the particles in the scheme, namely 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 aggregated GSD in dispersion model simulations of the Eyjafjallajökull (2010) eruption had very little impact on the modelled extent of the distal ash cloud with mass loadings significant for aviation. However, our scheme neither represents 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.

Appendix A: Scaling analysis

In order to gain more insight into the dependence of the collision kernel K_j,k on q, St_cr, and ρ_s, we carry out a scale analysis of α_k,j and β_k,j in turn. Starting with the collision rate, we can write Eq. (33) as follows:

$$\begin{array}{}\text{(A1)}& {\displaystyle }\begin{array}{rl}{\mathit{\beta }}_{k,j}=& \mathcal{B}{\displaystyle \frac{(d_k+d_j{)}^{\mathrm{2}}}{d_k{d}_j}}\\ +& \mathcal{S}(d_k+d_j{)}^{\mathrm{3}}+\mathcal{I}(d_k+d_j{)}^{\mathrm{2}}|d_k^{\mathrm{2}}-d_j^{\mathrm{2}}|\\ +& \mathcal{D}(d_k+d_j{)}^{\mathrm{2}}|d_k^{\mathrm{2}}-d_j^{\mathrm{2}}|,\end{array}\end{array}$$

where, in the following:

$$\begin{array}{}\text{(A2)}& {\displaystyle }\begin{array}{rl}& \mathcal{B}={\displaystyle \frac{\mathrm{2}{k}_{\text{B}}T}{\mathrm{3}{\mathit{\mu }}_{\text{a}}}}\\ & \mathcal{S}={\displaystyle \frac{\mathrm{1}}{\mathrm{8}}}{\left({\displaystyle \frac{\mathrm{1.7}\mathit{ϵ}}{{\mathit{\nu }}_{\text{a}}}}\right)}^{\mathrm{1}/\mathrm{2}}\\ & \mathcal{I}={\displaystyle \frac{\mathrm{1}}{\mathrm{72}}}{\displaystyle \frac{\mathit{\pi }{\mathit{ϵ}}^{\mathrm{3}/\mathrm{4}}}{{\mathit{\nu }}_{\text{a}}^{\mathrm{1}/\mathrm{4}}}}{\displaystyle \frac{{\mathit{\rho }}_{\text{s}}}{{\mathit{\mu }}_{\text{a}}}}\\ & \mathcal{D}={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}{\displaystyle \frac{g{\mathit{\rho }}_{\text{s}}}{\mathrm{18}{\mathit{\mu }}_{\text{a}}}}\end{array}\end{array}$$

are taken to be constant (including ρ_s). Here we have assumed that the particles settle with Stokes' terminal velocity (i.e. we neglect the second term on the right-hand side of Eq. (30); this will lead to quantitative discrepancies with the collision kernel calculated in Sect. 3 for larger diameters, but the qualitative behaviour will be correct). We have also assumed that ${\mathit{\beta }}_{i,j}^{\text{TS}}>{\mathit{\beta }}_{i,j}^{\text{LS}}$; this assumption does not affect our conclusions below. Since aggregation is associated with the presence of liquid water or ice, and α_k,j only depends on q and St_cr in the presence of liquid water, we choose values of the constituent parameters in Eq. (A2) that are appropriate for this case. Thus, with T=300 K, ϵ=0.01 m² s⁻³, ρ_a=1.297 kg m⁻³, and ρ_s=2000 kg m⁻³, the constants in Eq. (A2) have the following orders of magnitude:

$$\begin{array}{}\text{(A3)}& {\displaystyle }\begin{array}{rl}& \mathcal{B}\sim {\mathrm{10}}^{-\mathrm{16}}{\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}\\ & \mathcal{S}\sim \mathrm{1}{\mathrm{s}}^{-\mathrm{1}}\\ & \mathcal{I}\sim {\mathrm{10}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}\\ & \mathcal{D}\sim {\mathrm{10}}^{\mathrm{7}}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}.\end{array}\end{array}$$

As in Sect. 3, we restrict attention to diameters in the range [1,10 000] µm. Figure 4 shows the variation in β_k,j, given by Eq. (A1), 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. (29)–(31) becomes clear for large diameters. Note that β_k,j is symmetric in the indices j and k.

In the special case that d_k=d_j, Eq. (A1) becomes the following:

$$\begin{array}{}\text{(A4)}& {\mathit{\beta }}_{j,j}=\mathrm{4}\mathcal{B}+\mathrm{8}\mathcal{S}{d}_j^{\mathrm{3}}.\end{array}$$

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 obtain ${\mathit{\beta }}_{j,j}\sim {\mathrm{10}}^{-\mathrm{14}}$ m³ s⁻¹, for d_j∼100 µm we obtain ${\mathit{\beta }}_{j,j}\sim {\mathrm{10}}^{-\mathrm{11}}$ m³ s⁻¹, and for d_j∼1000 µm we obtain ${\mathit{\beta }}_{j,j}\sim {\mathrm{10}}^{-\mathrm{8}}$ m³ s⁻¹.

In the case that d_k≪d_j, Eq. (A1) becomes the following:

$$\begin{array}{}\text{(A5)}& {\mathit{\beta }}_{k,j}\approx \mathcal{B}{\displaystyle \frac{d_j}{d_k}}+\mathcal{S}{d}_j^{\mathrm{3}}+\mathcal{I}{d}_j^{\mathrm{4}}+\mathcal{D}{d}_j^{\mathrm{4}}.\end{array}$$

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\in [\mathrm{1},\mathrm{10}\mathrm{000}]$ µ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 obtain ${\mathit{\beta }}_{k,j}\sim {\mathrm{10}}^{-\mathrm{13}}$ m³ s⁻¹, for d_j∼100 µm we obtain ${\mathit{\beta }}_{k,j}\sim {\mathrm{10}}^{-\mathrm{9}}$ m³ s⁻¹, and for d_j∼1000 µm we obtain ${\mathit{\beta }}_{k,j}\sim {\mathrm{10}}^{-\mathrm{5}}$ m³ s⁻¹. These are consistent with what is observed in Fig. 4 for d_k≪d_j. Furthermore, we note that these values are all larger than the values of β_k,j in the special case d_k=d_j, and that this difference increases as d_j increases in magnitude. This explains the kink in Fig. 4, when d_k=d_j, and why it becomes sharper as d_j increases.

For d_k≫d_j, the scale analysis shows that, in the following:

$$\begin{array}{}\text{(A6)}& {\mathit{\beta }}_{k,j}\approx \mathcal{B}{\displaystyle \frac{d_k}{d_j}}+\mathcal{S}{d}_k^{\mathrm{4}}.\end{array}$$

For $d_j\in [\mathrm{1},\mathrm{10}\mathrm{000}]$ µm and d_k≫d_j, the second term dominates. Thus, this leads to the order of ${\mathit{\beta }}_{k,j}\propto d_k^{\mathrm{4}}$ for d_k≫d_j, which is consistent with what is observed in Fig. 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. (29)–(31).

We now turn to the sticking efficiency. On making use of Eqs. (34) and (35), we can write St_v as follows:

$$\begin{array}{}\text{(A7)}& S{t}_{\text{v}}=V|d_j-d_k|d_j{d}_k+{\displaystyle \frac{B}{d_j+d_k}}+S{d}_j{d}_k,\end{array}$$

where, in the following:

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}& V={\displaystyle \frac{\mathrm{4}g{\mathit{\rho }}_{\text{s}}^{\mathrm{2}}}{\mathrm{81}{\mathit{\mu }}_{\text{l}}{\mathit{\mu }}_{\text{a}}}}\\ & B={\displaystyle \frac{\mathrm{64}{\mathit{\rho }}_{\text{s}}{k}_{\text{B}}T}{\mathrm{27}\mathit{\pi }{\mathit{\mu }}_{\text{l}}{\mathit{\mu }}_{\text{a}}}}\\ & S={\displaystyle \frac{\mathrm{4}{\mathit{\rho }}_{\text{s}}}{\mathrm{9}\mathit{\pi }{\mathit{\mu }}_{\text{l}}}}{\left({\displaystyle \frac{\mathrm{1.7}\mathit{ϵ}}{{\mathit{\nu }}_{\text{a}}}}\right)}^{\mathrm{1}/\mathrm{2}}\end{array}\end{array}$$

are assumed to be constant, and we have also assumed that the two colliding particles have the same density (as in Sect. 3). Using the same values of the parameters as above (with constant ρ_s), the constants have the following orders of magnitude:

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}& V\sim {\mathrm{10}}^{\mathrm{14}}{\mathrm{m}}^{-\mathrm{3}}\\ & B\sim {\mathrm{10}}^{-\mathrm{10}}{\mathrm{m}}^{-\mathrm{1}}\\ & S\sim {\mathrm{10}}^{\mathrm{7}}{\mathrm{m}}^{-\mathrm{2}}.\end{array}\end{array}$$

Note that St_v and, hence, α_k,j are symmetric in the indices k and j. The ranges of q and St_cr we consider are the same as those in Sect. 3, i.e. $q\in [\mathrm{0.4},\mathrm{1.6}]$ and $S{t}_{\text{cr}}\in [\mathrm{0.65},\mathrm{2.6}]$.

Consider first the special case d_k=d_j. Then, Eq. (A7) becomes the following:

$$\begin{array}{}\text{(A10)}& S{t}_{\text{v}}={\displaystyle \frac{B}{\mathrm{2}{d}_j}}+S{d}_j^{\mathrm{2}}.\end{array}$$

The first term dominates for d_j≲1 µm; the second term dominates for d_j≳10 µm. For d_j⩽100 µm, St_v≪St_cr and so, in the following:

$$\begin{array}{}\text{(A11)}& {\mathit{\alpha }}_{k,j}\approx \mathrm{1}-{\left({\displaystyle \frac{S{t}_{\text{v}}}{S{t}_{\text{cr}}}}\right)}^q\approx \mathrm{1},\end{array}$$

for $q\cancel{\ll }\mathrm{1}$, whereas for d_j⩾1000 µm, St_v≫St_cr and so, in the following:

$$\begin{array}{}\text{(A12)}& {\mathit{\alpha }}_{k,j}\approx {\left({\displaystyle \frac{S{t}_{\text{cr}}}{S{t}_{\text{v}}}}\right)}^q\left(\mathrm{1}-{\left({\displaystyle \frac{S{t}_{\text{cr}}}{S{t}_{\text{v}}}}\right)}^q\right)\ll \mathrm{1},\end{array}$$

for $q\cancel{\ll }\mathrm{1}$.

In the case d_k≪d_j, Eq. (A7) becomes the following:

$$\begin{array}{}\text{(A13)}& S{t}_{\text{v}}=V{d}_j^{\mathrm{2}}{d}_k+{\displaystyle \frac{B}{d_j}}+S{d}_j{d}_k.\end{array}$$

Again, we fix d_j and allow d_k to vary. The smallest admissible values of d_j and d_k that are in the range [1,10 000] µm are d_j∼10 µm and d_k∼1 µm; for these values, a scale analysis shows that the first term on the right-hand side is dominant. For d_j∼10⁴ µm, the largest admissible value, and $d_k\in [\mathrm{1},\mathrm{1000}]$ µm, a scale analysis also shows that the first term on the right-hand side is dominant. Thus, Eq. (38) becomes the following:

$$\begin{array}{}\text{(A14)}& {\mathit{\alpha }}_{k,j}\approx {\left({\displaystyle \frac{S{t}_{\text{cr}}}{V{d}_j^{\mathrm{2}}{d}_k}}\right)}^q\left(\mathrm{1}-{\left({\displaystyle \frac{S{t}_{\text{cr}}}{V{d}_j^{\mathrm{2}}{d}_k}}\right)}^q\right),\end{array}$$

if $V{d}_j^{\mathrm{2}}{d}_k/S{t}_{\text{cr}}>\mathrm{1}$, whereas, in the following:

$$\begin{array}{}\text{(A15)}& {\mathit{\alpha }}_{k,j}\approx \mathrm{1}-{\left({\displaystyle \frac{V{d}_j^{\mathrm{2}}{d}_k}{S{t}_{\text{cr}}}}\right)}^q,\end{array}$$

if $V{d}_j^{\mathrm{2}}{d}_k/S{t}_{\text{cr}}<\mathrm{1}$. Since St_cr is always of the order of unity, if $V{d}_j^{\mathrm{2}}{d}_k\ll \mathrm{1}$, then α≈1. As d_j increases, the range of d_k values for which $V{d}_j^{\mathrm{2}}{d}_k\ll \mathrm{1}$ decreases. For d_j≳100 µm, we see that $V{d}_j^{\mathrm{2}}{d}_k\mathit{\gtrsim }\mathrm{1}$ for all admissible values of d_k.

In the case d_k≫d_j, a similar scale analysis to that above shows that the first term on the right-hand side of Eq. (A7) is again the dominant term; it now takes the form $V{d}_k^{\mathrm{2}}{d}_j$. For d_j∼1 µm, $V{d}_k^{\mathrm{2}}{d}_j\mathit{\gtrsim }\mathrm{1}$, if d_k≳100 µm. For d_j∼1000 µm, it follows from the condition d_k≫d_j that d_k∼10 000 µm, and so $V{d}_k^{\mathrm{2}}{d}_j\mathit{\gtrsim }\mathrm{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}_k^{\mathrm{2}}{d}_j\mathit{\gtrsim }\mathrm{1}$ and α<1 increases (e.g. if d_j∼100 µm, then d_k≳10 µm for $V{d}_k^{\mathrm{2}}{d}_j\mathit{\gtrsim }\mathrm{1}$ to hold).

Consider first the variation in α_k,j with St_cr (Fig. 5a–c): when d_j=10 µm (Fig. 5a), then, for both d_k<d_j and $d_k\cancel{\gg }{d}_j$, it can be seen that α≈1, as expected, since here $V{d}_j^{\mathrm{2}}{d}_k\ll \mathrm{1}$ and $V{d}_k^{\mathrm{2}}{d}_j\ll \mathrm{1}$, for $d_k\cancel{\gg }\mathrm{10}$ µm. For increasing d_k≫d_j, we see α_k,j decreasing for all values of St_cr, since here $V{d}_k^{\mathrm{2}}{d}_j\gg \mathrm{1}$. For a given value of d_k≫d_j, Fig. 5a shows that α_k,j increases by approximately 4^q over the range of St_cr shown. A similar pattern can be seen for d_k≫d_j in Fig. 5b and c. Compared with Fig. 5a, Fig. 5 b and c show more variation with St_cr for d_k≪d_j and decreasing values of α_k,j; this occurs because $V{d}_j^{\mathrm{2}}{d}_k$ increases with increasing d_j, and, of course, the range of d_k values satisfying d_k≪d_j also increases with increasing d_j.

Turning now to the variation of α_k,j with q shown in Fig. 5d–f. Raising $V{d}_k^{\mathrm{2}}{d}_j>\mathrm{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}_j^{\mathrm{2}}{d}_k$). Thus, for example, as shown in Fig. 5d when d_j=10 µm and d_k≫d_j, $V{d}_k^{\mathrm{2}}{d}_j>\mathrm{1}$ for d_k≳100 µm, and so, for q>1, α_k,j is smaller than it would be for q=1, whereas, for q<1, it is larger. These patterns hold true in Fig. 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 Fig. 5f, for example, that, for d_k≪d_j, $V{d}_j^{\mathrm{2}}{d}_k>\mathrm{1}$ for all values of $d_k\in [\mathrm{1},\mathrm{100}]$ µm, and so α_k,j is closer to unity for q<1 and vice versa for q>1.

It should be noted that α_k,j, computed with Stokes' terminal velocity, shows a larger variation than that shown in Fig. 5. Since Stokes' terminal velocity is larger than that calculated from Eqs. (29)–(31) for large diameters, and α_k,j is dominated by differential settling, then, for large diameters, α_k,j becomes smaller than the values shown in Fig. 5. Conversely, Fig. 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 Fig. 6.