Estimating the 3D shape of volcanic ash to better understand sedimentation processes and improve atmospheric dispersion modelling

Abstract The sedimentation rate of volcanic ash through the atmosphere influences its travel distance, with important implications for aviation and health. The fall velocity of a particle depends on its size and density, but also shape, and volcanic ash is not spherical. To capture the sedimentation of ash, atmospheric dispersion models use empirical drag equations calibrated using geometric shape descriptors. However, particle shape data are scarce and there is no standard method of shape measurement. In addition, shape measurements are not always available during an eruption, when dispersion models are used operationally to forecast ash hazard. We assess the variability in the shape of volcanic ash from Icelandic eruptions using X-ray computed tomography. To consider how good different drag equations and shape descriptors are at representing the sedimentation of volcanic ash we compare calculated fall velocities to measured fall velocities of volcanic ash in air in a settling column. We then suggest the best drag equations and shape descriptors for use in atmospheric dispersion models. We find that shape-dependent drag equations produce more accurate results than a spherical approximation. However, accurate drag calculations based on the shape descriptor sphericity, which is a function of surface area, require the imaging resolution to be within the range of 102 - 105 voxels per particle (where a voxel is a volumetric pixel) as surface area is sensitive to imaging resolution. We suggest that the large-scale form of the particle impacts sedimentation more than small-scale surface roughness. Shape descriptors based on ratios between principal axis lengths are more practical as they are less variable among particle size classes and much less sensitive to imaging resolution. Finally, we use particle shape data from this study and literature sources to make recommendations on default values for use with atmospheric dispersion models where no shape data are available.


Introduction 9
Volcanic ash (tephra particles with diameters < 2 mm) can remain in the atmosphere 10 from minutes to days or longer after a large explosive eruption (Durant et al., 2010); the 11 rate of removal from the atmosphere depends on meteorological processes and particle ter-12 minal fall velocity. Accurate calculation of terminal velocity requires a parameterisation of 13 particle shape (e.g., Riley et al., 2003;Bagheri and Bonadonna, 2016a,b). This means that 14 forecasts of the movement of volcanic ash clouds generated using atmospheric dispersion  The question of how to measure, and model, particle shape is relevant in many other 24 disciplines. Many atmospheric particulates are non-spherical, meaning that shape is an behaviour of particles settling in water is also a function of shape; this is important, for 29 example, in analysing assemblages of marine species such as foraminifera, which are used 30 to estimate palaeodepth (e.g., Speijer et al., 2008;Jorry et al., 2006). Particle shape also 31 affects the drag of non-spherical particles in streams (Komar and Reimers, 1978) as well 32 as particle sorting during flow (Oakey et al., 2005). Finally, shape affects crystal settling 33 velocities in magma, with implications for magma evolution and rheology (Higgins and 34 Roberge, 2003). 35 Another consideration relates to the effect of particle shape on grain size measurement. 36 For example, particles passing through a sieve mesh may have one dimension larger than 37 the mesh aperture, and so the results of sieving are dependent on the shape of the particle 38 (Arason et al., 2011); shape also affects the results of particle size distributions calculated 39 from 2D images (Higgins and Roberge, 2003). Particle volume is a shape-independent 40 measure of size, but volume is not straightforward to measure for irregular particles.
ters. We use the resulting insights to present a database of shape descriptors for use with 73 semi-empirical drag equations (Ganser, 1993

Modelling the terminal velocity of non-spherical ash
Particle terminal velocity is defined as the velocity reached by a falling object when the 78 drag force is equal to the gravitational force and acceleration is zero. Volcanic ash particles 79 reach terminal velocity in the atmosphere over distances which are small compared to the 80 distance required to sediment from a plume, and so it is reasonable to neglect acceleration 81 in modelling. Terminal velocity w t can be calculated as a function of drag: (1) where d is the particle size (for a sphere, diameter), g is gravitational acceleration, ρ is fluid 83 density, ρ P is particle density, and C D is the drag coefficient, a dimensionless coefficient 84 which is a function of particle shape and flow regime. Volcanic ash falling in air can be 85 subject to several flow regimes, defined by the dimensionless Reynolds number Re; the 86 flow around a particle is classed as laminar when Re < 0.1 and turbulent at Re > 1000. 87 The drag of spheres can be calculated analytically with high accuracy for all flow regimes 88 (e.g., White, 1974). Solutions for non-spherical particles, which are characterised by higher 89 C D than spheres of equivalent size and density, are generally empirical or semi-empirical 90 correlations which relate C D to one or more geometric shape descriptors. Therefore, such 91 correlations are valid for finite Re and limited to particle shape ranges which are covered 92 by the experimental conditions and the formulation used. The complexity in modelling 93 volcanic ash shape means that many operational atmospheric dispersion model setups use 94 a spherical approximation by default (Hort, 2016

118
Prior to analysis, all samples were manually dry sieved at 1 or 0.5 φ intervals (where φ 119 is a measure of grain size defined as −log 2 D D 0 , where D is the diameter of the particle in 120 millimetres and D 0 is a reference diameter of 1mm). We used sieve mesh diameters of 4 to 121 -1 φ (62.5 to 2000 µm), apart from sample KVE which was sieved at half-φ intervals using 122 sieve mesh diameters of 3.5 to -0.5 φ (88 to 1414 µm). 123 We measured the shape of 19557 particles from KVE, KSU, and HEK using X-ray CT 124 to provide a large volcanic ash shape database, including shape distributions for each sieve 125 fraction, as shape can vary with size (e.g., Mele and Dioguardi, 2018). For this analysis the 126 samples were analysed in bulk, including crystals and lithics. We chose only three samples 127 for this analysis due to the significant scan time needed; the samples we chose span a range 128 of qualitative shape characteristics (see Table 1).

129
To determine the effectiveness of our measured shape descriptors when used to calcu-130 late terminal fall velocity, we also selected 46 individual particles of juvenile glass from the 131 0 φ and -0.5 φ sieve fractions of the KSM, KVE, HEK and EYJ samples, which were indi-132 vidually scanned using X-ray CT. Particles chosen ranged from 1.0 -2.6 mm, sufficiently 133 large to image and track; again grains were selected to include a wide range of shapes. As 134 density is a crucial input in drag equations, we also measured the density of these particles. 135 We then considered how well we could calculate the fall velocity of our particles using   particle's trajectory. The error arising from the relative positioning of the camera, particle, 158 and ruler was corrected assuming that each particle was falling in the centre of the 15 cm 159 main settling column. Each particle velocity measurement presented here represents the

166
We calculated particle density using volumes from X-ray CT scans and particle mass 167 measured on a balance with a precision of 0.0001 g. Particle dimensions calculated using 168 Avizo CT software were checked using digital calipers, to ensure that CT data gave accurate

198
It is important to note that benchtop X-ray CT cannot be used to accurately quantify 199 the shape of the finest volcanic ash fractions relevant to aircraft hazard (< 30 to 60 µm; to constraints on imaging resolution. Although synchrotron X-ray CT systems can achieve 202 resolutions of 1 µm or less, in our system the minimum voxel edge length is ∼ 3 µm.

203
For this reason, particles smaller than 62.5 µm (4 φ) were not used, and our X-ray CT  to calculate drag by these equations, we measured the shape descriptors used in their design.

219
The drag equations are all applicable for the range of flow regimes expected for volcanic 220 ash falling in air and are all calibrated using 3D geometric shape descriptors. We calculate Material. Since the drag of spheres can be determined analytically to high accuracy, a 225 popular approach in defining shape descriptors is to use a ratio of a particle parameter to 226 that of a volume-equivalent sphere. The Ganser (1993) drag equation uses sphericity ψ G , 227 the ratio of surface area of a volume-equivalent sphere to surface area of the particle being 228 described: where A surf is a measure of 3D surface area and therefore an effective descriptor of rough- with equal principal axes to the particle: where z = 1.6075. The Dioguardi et al. (2018) drag equation is calibrated using an ellipsoid 235 approximation; their shape descriptor Ψ D , which we term the Dioguardi shape factor, is the 236 ratio of 3D sphericity to 2D circularity: and 239 X = P proj P c ; (6) P proj = maximum projected perimeter and P c = the perimeter of a circle with equal projected 240 area to the particle being described. In this study we focus on 3D shape measurement and 241 do not measure X, as this 2D parameter is a function of particle perimeter, one of the 2D The shape descriptors we measure (sphericity, the Dioguardi shape factor, elongation Particle surfaces are reconstructed from raw CT data by separating 3D regions rep-257 resenting particles from the surrounding epoxy; this requires the selection of a threshold 258 greyscale value. As the particle edges are characterised by a gradient (over ∼ 3 -5 voxels) 259 rather than a sharp boundary, the choice of threshold is subjective and so we determined 260 the sensitivity of particle volume and shape to this choice. We did this by increasing and 261 decreasing our best estimate greyscale threshold by 10, which covered the particle bound-262 ary gradients; the results and example greyscale images are given in Supplementary Figure   263 A1. We calculate a maximum 6% error on mean d v and 4% error on mean sphericity arising 264 from the selection of the particle boundaries. As shape can be sensitive to imaging resolution, we investigated the impact of voxel 267 size on measured shape factors using X-ray CT data for 16 individual particles. Scans were 268 conducted at the scanner's maximum resolution (a voxel edge length of ∼ 3 µm, giving 269 between 6.5 × 10 6 and 1.3 × 10 8 voxels per particle). We progressively resampled the scan 270 data from 2 to 32 times the original voxel edge lengths, giving a maximum voxel edge 271 length of ∼ 96 µm; after each resampling we recalculated surface area, volume, sphericity, 272 the Dioguardi shape factor, elongation, and flatness.

273
The results are shown in Figure 3 and highlight the sensitivity of surface area, and the 274 surface area-based shape factor sphericity, to resolution. For a large, rough particle (KSM- Higher scan resolutions result in higher surface area measurements, which give lower sphericity values and therefore calculated velocity which is too low.

Effectiveness of published drag equations 318
The results of our comparison between drag equations are shown in Figure 5. The an-   (Ganser, 1993) where sphericity is calculated using CT scan resolutions of 96 µm voxel edge length. Note that the Ganser (1993) scheme performs worse than even a spherical particle approximation when sphericity is calculated using much higher image resolutions. minal velocity. In practice, when forecasting ash dispersion operationally, information on 336 particle morphology is unlikely to be available. Therefore, the use of these schemes re-337 quires a database from which to choose a default ash shape. It is important to assess the relationship between particle shape and size, as well as 340 to obtain shape data for particles smaller than the mm-sized particles we used in settling 341 column experiments. Operationally, Volcanic Ash Advisory Centres use volcanic ash dis-342 persion models with a particle size distribution (PSD); for most operational systems the 343 bulk of the modelled erupted mass is restricted to particles with diameter < 100 µm (Hort,

A database of the shape of volcanic ash
As illustrated above, another important consideration is the resolution used to obtain 346 sphericity ( Figure 3). Accurate calculation of terminal velocity using the Ganser (1993) 347 scheme requires sphericity to be calculated from CT data with a resolution which gives 348 between 10 2 and 10 5 voxels / particle, and so our database can only include sphericity data 349 in this resolution range. Because the Dioguardi shape factor, elongation and flatness are not 350 sensitive to imaging resolution, it is possible to directly compare our data to other studies 351 measuring the same shape factors, regardless of experimental conditions. We calculated 352 these shape factors using equations 7, 8, and 9, where only axis lengths or approximate 353 ellipsoid sphericity ψ e were reported. We also exclude studies which use different shape 354 equations; for example, many studies calculate approximate sphericity using 2D images. 355 We prefer to limit our analysis to the exact shape descriptors by which the drag equations 356 were calibrated. 357 We obtain shape descriptors for every sieve fraction of samples KVE, KSU, and HEK, 358 adjusting voxel size, with a resolution of 4728 -20,000 voxels / particle. The resulting par-359 ticle shape distributions are shown in Figure 6; we compare them in Figure 7 to published 360 data from other eruptions that match our criteria. All particle shape data are available in  Figure 6: Shape distributions for bulk sieved samples HEK, KVE, and KSU. We chose three samples due to the scan time needed; the samples span a range of qualitative particle shape characteristics (Table 1). We shade each plot to highlight shape factors < 0.5 to aid visual comparison of the proportion of highly non-spherical particles in each sample. Red lines indicate the median; boxes show the interquartile range; crosses indicate outliers. 386 We compare our shape data to data from previous studies, to expand our morphology 387 database to include different eruptions and particle size fractions, and to determine whether 388 the eruptions studied here show a 'typical' range of volcanic ash shape. All ash shape data 389 are given in Figure 7. As data vary between studies, we plot only mean and standard 390 deviation of shape for each sample. Some studies report the mean shape factor for each 391 size fraction of a sample; others give the mean shape of a bulk sample; where the particles 392 vary in size, we indicate the size range using the X-axis error bars. We note that methods 393 of measuring grain size differ between studies, and so the size ranges shown here should 394 be considered approximate. For the shape descriptor sphericity, we include only data from 395 studies that use CT data with a resolution of between 10 2 and 10 5 voxels/particle, which 396 is our recommendation. Despite this limitation we still find a weak correlation between 397 image resolution and sphericity (Figure 7b).  from the other ash shape data in any mean shape factor. 411 6. Discussion 6.1. Measuring shape 413 X-ray CT is an accurate and efficient method of assessing particle size and shape in 3D, 414 and allows imaging of hundreds of particles, and multiple shape factors, relatively rapidly.

415
The analytic error resulting from manual selection of greyscale values is low (< 4 % on 416 sphericity, which is insignificant compared to its sensitivity to image resolution, and < 6 417 % on diameter d v ). The high resolution makes it an invaluable tool for examining detailed 418 structures in volcanic rocks. For shape quantification, however, high resolution surface area 419 measurements result in very low sphericity and therefore underestimate terminal velocity 420 by the drag equation of Ganser (1993). We recommend using a resolution between 10 2 and 421 10 5 voxels/particle to calculate sphericity. The best agreement between measured veloc-422 ity w t, rec and calculated velocity w t, Ganser is reached at a range of image resolutions; we 423 suggest this is due to the range of particle shapes as well as uncertainty on other param- resolution cannot be kept constant. We conclude that above our lower resolution limit of ∼ 428 10 2 voxels/particle, imaging resolution is not a concern for calculation of these shape de-429 scriptors, meaning that for calculation of elongation, flatness or the Dioguardi shape factor 430 it is practical to sacrifice higher resolution in favour of speed.

431
As shape parameters based on axis lengths are less sensitive to resolution, we assess Shape factors based on principal axis lengths (elongation and flatness) do not change significantly with particle size or between eruptions (Figures 6 and 7). The exceptions are  For volcanic ash we found that an imaging resolution of 10 2 -10 5 voxels per parti-502 cle is required for determining surface-area-dependent shape parameters for accurate drag 503 calculation. This range may extend to higher resolutions if particles are smoother. 504 We suggest that shape-dependent drag equations should also be evaluated for mod-

Recommendations for including ash shape in dispersion models
We present a table of default shapes to be used with shape-dependent drag equations for 510 modelling atmospheric ash concentrations and travel distance ( Table 2). The shape ranges 511 given are based on the mean and standard deviation of ash shape (Figure 7). Bagheri and forecasting the dispersion of ash in particle size ranges typically modelled by VAACs; and 524 our shape data are calculated from axis lengths, which are insensitive to image resolution, 525 or surface area at resolutions we find to be effective for drag calculation. We recommend 526 the use of these shape values as defaults in place of a spherical approximation in volcanic 527 ash dispersion models.
A surf , or sphericity of an approximate ellipsoid: