Sensitivity of Volcanic Ash Dispersion Modelling to Input Grain Size Distribution Based on Hydromagmatic and Magmatic Deposits

: The size distribution of volcanic ash is rarely measured in real time and Volcanic Ash Advisory Centres (VAACs) often rely on a default particle size distribution (PSD) to initialise their dispersion models when forecasting the movement of ash clouds. We conducted a sensitivity study to investigate the impact of PSD on model output and consider how best to apply default PSDs in operational dispersion modelling. Compiled grain size data confirm that, when considering particles likely to be in the distal ash cloud (< 125 μ m diameter), magma composition and eruption size are the dominant controls on grain size distribution. Constraining the PSD is challenging but we find that the grain size of deposits from large hydromagmatic eruptions remains relatively constant with distance, suggesting that total (whole ‐ deposit) grain size distributions (TGSDs) for these eruptions could be estimated from a few samples. We investigated the sensitivity of modelled ash mass loadings (in the air and on the ground) to input PSDs based on coarse to fine TGSDs from our dataset. We found clear differences between modelled mass loadings and the extent of the plume. Comparing TGSDs based on ground ‐ only and ground ‐ plus ‐ satellite data for the Eyjafjallajökull 2010 eruption, we found that basing input PSDs on TGSDs from deposits alone (likely missing the finest particles) led to lower modelled peak ash concentrations and a smaller plume.


Introduction
Ash clouds from explosive volcanic eruptions can travel huge distances and cause severe disruption to air traffic (e.g., [1]). To mitigate against aircraft encounters with ash clouds, Volcanic Ash Advisory Centres (VAACs) issue Volcanic Ash Advisories (VAAs) and Volcanic Ash Graphics (VAGs) which indicate the expected location of the ash cloud up to 18 h ahead. In Europe, there is an additional requirement to provide forecasts of ash concentrations [2]. As part of their forecasting process, the VAACs use atmospheric dispersion models to simulate the transport, dispersion and removal of ash from the atmosphere [3][4][5]. Where ash falls as individual particles, residence time in the atmosphere is controlled by particle parameters (grain size, density and shape), vertical advection and turbulence. Grain size is typically the dominant particle parameter in controlling the ash settling velocity and whether the ash residence time is significantly affected by turbulence [6][7][8][9][10]. It is, therefore, important to select an appropriate particle size distribution (PSD) to initialise model runs, but data on the size distribution of ash are rarely available in near real time and so most VAACs use default PSDs (e.g., [3,11,12]). In addition, many assume that only the smallest particles survive nearsource removal processes, i.e., they only consider particles with diameters, e.g., ~< 100 μm (e.g., [3,4]). All of the ash in the air eventually falls to the ground and previous studies of the deposits show that grain size distribution varies greatly between eruptions, reflecting eruption intensity, magma composition and whether the eruption is wet or dry (hydromagmatic or magmatic) (e.g., [13,14]). It is still not clear how best to set a default PSD when measurements are not available. In this study, we test the sensitivity of modelled ash plumes to the input PSD and investigate how best to apply default PSDs in operational dispersion modelling.
The PSD of volcanic ash can be determined using a variety of methods but each approach has its limitations. The ash cloud can be directly sampled with research aircraft, but sampling is usually restricted to relatively dilute portions of the cloud over short periods (e.g., [15][16][17][18]). Remote sensing can provide data on plume evolution, but standard retrieval techniques are only sensitive to particles up to ~32 μm diameter (e.g., [19,20]). As the PSD of ash in the air is reflected in the PSDs of deposits at different distances from the source, ground sampling is also a useful technique. The total grain size distribution (TGSD) of the deposit can be determined from multiple samples (e.g., [21,22]), but compilation of a reliable TGSD requires well-preserved deposits and spatially extensive sampling, particularly when the range of grain sizes is large [23][24][25]. Sampling the whole extent of a deposit can be difficult because, in the short term, there can be safety concerns in very proximal areas, whereas over time, preservation will deteriorate in distal areas. In any case, deposit TGSDs do not include the most distal material that is inaccessible (e.g., in the sea) or too thin and fine-grained to be visible (e.g., [26]). It is traditional to use the term total grain size distribution (TGSD) when referring to the size distribution of grains of an entire volcanic deposit and to use particle size distribution (PSD) for dispersion model inputs. Of course an ash grain is a type of particle but we have kept the distinction in this paper to make clear whether we are referring to model inputs (PSD) or deposit data (TGSD); in particular, we use deposit TGSDs as the basis for model input PSDs.
We collate TGSDs from the literature and investigate dispersion model sensitivity to input grain size distribution, considering both hydromagmatic and magmatic eruptions. We use the fine fraction of TGSDs (< 125 μm diameter) to initialise our model, as these particles travel the furthest and are therefore the most important for modelling the long-range transport of volcanic ash. Deposit TGSDs are determined after an eruption has ended (unlike direct sampling in the air and remote sensing) and they can therefore be used to investigate ash emissions from past eruptions. Importantly, deposits provide information on grain size and transport of ash from eruptions from potentially active volcanoes that have not erupted historically and eruptions of higher intensity and magnitude than have occurred since satellite remote sensing of ash was possible.
Unlike magmatic eruptions, where the grain size becomes more fine-grained with increasing distance from the vent (reflecting sorting by particle terminal velocity until the ash is sufficiently fine that turbulence dominates), when water interacts with magma (hydromagmatic conditions), grain size in the deposit is often relatively constant with distance [14]. This raises the possibility that TGSDs for hydromagmatic eruptions could be constructed from a limited set of samples. The lack of a fining trend with distance from source may result from more efficient initial fragmentation conditions [27] and/or distinct sedimentation processes. For example, water can enhance aggregation, with small particles falling out in clusters close to source (e.g., [28,29]) or sedimentation of fine particles can be driven by gravitational instabilities in the plume [30][31][32].
We undertake a sensitivity study of dispersion model forecasts to the PSD used to initialise the simulation by using PSDs representing the range of coarse to fine grain size distributions from our hydromagmatic and magmatic eruption TGSD dataset. We compare results to simulations using the default PSD used by the London VAAC during operational response. Almost all published TGSDs, and so almost all our model PSDs, are determined from deposit samples only. We also investigate the impact of the ash beyond the mapped deposit for the Eyjafjallajökull 2010 eruption for which there is a TGSD based on a combination of deposit (ground) and satellite data. For this case, we compare model results for PSDs based on ground-only and ground-plus-satellite data to identify the impact of the finest particles, while acknowledging that satellite retrievals have their own limitations (discussed further in Section 5.2). Based on our results, we make recommendations for how to apply default PSDs in operational atmospheric dispersion models. Our choice of dispersion model and model source location in Iceland, and our emphasis on Icelandic eruptions in the TGSD compilation, make the results particularly useful for the London VAAC, but the conclusions are more broadly applicable.

Compilation of Published Grain Size Data and Selection of PSDs for Modelling
We compiled published grain size distributions from hydromagmatic and magmatic eruptions (Tables 1 and 2, respectively). The compilation includes eruptions for which at least one of the following datasets was available:  TGSD, compiled from samples collected at a range of distances from source and including grains in the size range used by the VAAC default PSD (≤ 100 μm).  Grain size distributions or median grain size at a range of distances from the source.
We have limited our investigations to vent-derived plumes and eruptions of juvenile magma and so we excluded TGSDs from phreatic eruptions. We have not interrogated the quality of the grain size data, but see, for example, [24] and [23] for recent discussions of how sampling strategy affects the reliability of TGSD.
The eruptions (detailed in Supplementary Materials Section S1) range in composition from basaltic (e.g., Hverfjall 2000 BP) to rhyolitic (e.g., Taupo 130 CE) and in eruption size (as measured by the Volcanic Explosivity Index (VEI)) from VEI 2 (e.g., Hverfjall 2000 BP) to VEI 8 (e.g., Taupo 25.4 ka BP). For the magmatic eruptions, we have focused on Iceland, as it is within the London VAAC's area of responsibility (e.g., basaltic Eldgja 10th Century, andesitic Hekla 1693 and rhyolitic Askja 1875 Layer D) but include some other eruptions. TGSDs have been compiled from ground samples except for Eyjafjallajökull 2010, where TGSDs from ground observations alone and from ground observations plus satellite data are both available [33].
Data on median grain size were included because they provide a simple measure of the stability of grain size with distance from source and are reported as a function of distance for a wider range of hydromagmatic eruptions than full grain size distributions.  Tables 1  and 2, respectively. The eruptions cover a wide range of magma compositions from mafic (basalt, ≤ 52% SiO2) and intermediate (basaltic andesite-andesite, 52%-63% SiO2) to silicic (dacite-rhyolite, > 63% SiO2). Hydromagmatic eruptions show a smaller range of grain sizes than magmatic eruptions, with medians from 0 to 4 Φ (62.5-1000 μm), with the larger eruptions (which are also all silicic), e.g., Taupo Oruanui, Rotongaio and Hatepe and Towada (all VEI ≥ 6) having the finest TGSDs. The magmatic eruptions have median values of −2-4 Φ (62.5-4000 μm). However, when we consider only particles likely to be in the distal ash cloud (> 3 Φ, < 125 μm), TGSD is not a significant indicator of the presence of external water, with a large overlap between hydromagmatic and magmatic eruptions (Figure 3a). Magma composition and VEI seem to be more significant, with lower VEI and mafic or intermediate eruptions having generally coarser TGSDs than silicic ones (Figure 3b). The link between eruption size and magma composition is expected, both because melt viscosity (which depends on composition) affects vesiculation kinetics (e.g., [13]), and the largest explosive eruptions tend to be silicic. A limitation of this analysis is that we are not aware of any published TGSDs for large (VEI ≥5) mafic or intermediate-composition hydromagmatic eruptions.    Table 1). Magma composition: ▽ silicic, o mafic/intermediate.  Table 2). AMS = Campi Flegrei, Agnano Monte Spina. Magma composition: ▽ silicic, o mafic/intermediate.

Change in Grain Size with Distances from Source
Many of the hydromagmatic eruptions show little variation in PSD with distance, particularly for larger eruptions ( Figure 4). Similarly, when considering locations beyond the influence of near-source processes (> ~20 km from source), median grain size is relatively constant (or even coarsens slightly with distance) at tens or even hundreds of km from source for larger hydromagmatic eruptions The only dataset which has sufficient samples to allow a robust comparison of the hydromagmatic and magmatic phases from the same eruption is from Askja 1875 (units C and D, respectively; [36]). Figure 6 shows a more constant median for the hydromagmatic phase of this eruption than for the magmatic phase, where there is a significant fining with increasing distance from the vent. Other median grain sizes for the hydromagmatic and magmatic phases from the same eruption are included in Supplementary Materials Section S2; for Eldgja 10th century and Ilopango ~1.5 ka BP, samples from > 20 km from source are available for only one phase and for Grímsvötn 2011, there are only three hydromagmatic data points, prohibiting robust conclusions. Hydromagmatic total grain size distributions (TGSDs) and particle size distributions (PSDs) at varying distance from source. Colours indicate distance of PSD from source: < 10 km red/yellow; 10-50 km greens; > 50 km blues. TGSDs are in black. Dotted line shows the median grain size (i.e., 50% mass fraction). PSDs for Towada and El Chichón are only for grain size ≤ 4 Φ. TGSD for Ilopango Unit A [43] was not included as it seems inconsistent (too coarse) with grain size data from the individual sites.

Methods
We perform a sensitivity study using the dispersion model Numerical Atmospheric Dispersion Modelling Environment (NAME) [67] to investigate how the PSD used to initialise a simulation affects the model output. NAME is used by the London VAAC and has been extensively validated against observations, including satellite retrievals, lidar data and measurements taken by research aircraft (see for example, [16,[68][69][70][71]). We use coarse to fine grain size distributions from across the range of our dataset to initialise simulations and compare the results to simulations initialised with the default PSD currently used at the London VAAC.

The London VAAC Default PSD
The London VAAC uses a default PSD, given in Table 3, which is based on airborne measurements of an ash cloud taken on January 8, 1990 following the eruption of Mount Redoubt, USA [4]. Figure 7 shows the size distribution of ash particles in samples taken 2.5 and 7 h after the eruption [15]. The simplified curve used for the default NAME modelling, based on these samples is also shown, and includes linear extrapolation for particle diameters > 30 μm. To obtain the default PSD, the mass distribution for each particle size bin was calculated assuming particles are spheres of constant density (i.e., volume is proportional to mass). In addition, a 'token' percentage of mass was added in the 30-100 μm bin to account for large particles in the distal ash plume [72].  eruption. Curve B: 4.0 km above sea level (asl), ~130 km downwind and ~2.5 h after eruption. Curve C: 2.6 km asl, ~170 km downwind and ~7 h after eruption [15]. Met Office simplified plot is used for NAME default particle size distribution [from data in 72].

The Phi (Φ) scale for PSDs
NAME defines particle size in microns (μm) but grain size measurements of ash samples taken from deposits often use the phi (Φ) scale [73]: where D is the particle size in mm, and D0 is a reference value of 1 mm. To make it easier to compare results from NAME simulations using TGSDs defined in the phi scale, we compiled micron equivalents of whole-Φ and half-Φ scales (Tables S1 and S2). We then calculated mass fractions equivalent to the VAAC default PSD and compared results from NAME simulations of the Eyjafjallajökull 2010 eruption using the new scales and the current VAAC default. Air mass loadings showed very little difference between the three scales for all cases where the plume is dispersed beyond the immediate source area. Further details of the calculations and the tests are shown in Supplementary Materials Section S3; statistical tests are defined in Supplementary Materials Section S4.

Normalising TGSDs to ≤ 125 μm and Selecting Eruptions for Modelling
We selected TGSDs from four hydromagmatic eruptions that cover the range of coarse to fine TGSDs (Eldgja, Grímsvötn 2004, Rotongaio and Oruanui). In addition, one fine (Eyjafjallajökull 2010) and one coarse (Hekla 1991) magmatic TGSD were included. The Eyjafjallajökull 2010 TGSD is compiled from both ground samples and satellite data [33] and so, to investigate the impact of including satellite retrievals in TGSDs, we also added the Eyjafjallajökull 2010 TGSD derived only from ground samples (Figure 8).
The TGSDs were then normalised to include only grains ≤ 125 μm diameter to assist comparison of model outputs with those using the VAAC default PSD (≤ 100 μm). On the Φ scale, we selected reported grain sizes ≥ 4 Φ  Modelling Environment (NAME) model runs.

NAME Simulations
NAME was run on the JASMIN scientific data analysis environment [74] with the input parameters shown in Table 4. The PSD used to initialise the run was varied and all other eruption source parameters and met data were kept the same. Particles were assumed to be spherical and to have a density of 2300 kg m −3 (see Appendix A, Table A1 for a list of the PSDs considered).

Met data
Global configuration of the Unified Model, ~25 km horizontal resolution (mid-latitudes) and 3 hourly temporal resolution [76] NAME represents the long-range transport of volcanic ash in the atmosphere [67]. It also includes schemes to represent near-source processes such as buoyant plumes [77] and umbrella clouds [78] but does not explicitly represent all near-source processes. Particles only fall out individually under gravity, with no collective settling as aggregates or from density-driven ('en masse' or 'streaking') instabilities [4].
As we are investigating the extent to which input PSD influences cloud size and ash concentration, we are interested in relative rather than absolute changes in the cloud. For each simulation, particles were released above the vent of Eyjafjallajökull volcano with a uniform distribution, from the vent to the plume top. The source strength was calculated from reported plume heights [79] using the method of Mastin et al. [6], and we took just 5% of the calculated MER to represent the total erupted material in the distal plume. The buoyant plume and umbrella cloud schemes were not invoked. This follows the London VAAC default approach; it accounts for the largest particles falling out close to source and the fact that NAME does not represent all near-source processes [4].
Modelled ash mass loadings in the air and on the ground were output onto a grid with a horizontal resolution of 33 km and a 3 h temporal resolution. These were chosen as a balance between time taken for the simulations to run and resolution required to identify the impact of changing PSD on modelled ash mass loadings.
Results were compared to runs using the VAAC default PSD for the period 4-8 May 2010, chosen as this was the period relating to the published Eyjafjallajökull 2010 TGSD [33]. However, we were not aiming to replicate the Eyjafjallajökull 2010 cloud and deposit, rather, we used it as a benchmark from which to consider the sensitivity of the simulated plume to the input PSD. Within this time frame, total column mass loadings in the air (g m −2 ) were selected for comparison on 2 days: one with a narrow plume (6 May 12:00 UTC) and one with a more extensive plume (8 May 00:00 UTC). Deposition for the period 4-12 May was considered, as ~95% of the mass was found to deposit within this timeframe.

Quantifying the Impact on Dispersion Model Forecasts
To investigate the sensitivity of the model output to the input PSD used, we compared NAME simulations initialised using the VAAC default to those using our other selected PSDs. We compared ash concentration within the model plumes, the different spatial extents of the model plumes, and differences in deposit concentrations. No one statistical test can comprehensively describe differences for dispersion model results and so we selected four tests (used in previous studies [8,80] and Statistical analyses are usually used to compare model results with a limited set of observed values, and hence to make the tests meaningful we applied thresholds before comparing model runs with different input parameters. We considered ash mass loadings > 0.2 g m −2 , which is the minimum threshold represented by ash concentration charts (assuming a 1 km deep plume) and ash deposits with mass loadings > 0.1 g m −2 , chosen as it covers a similar area to the air concentrations.

Results
This section presents our results on the sensitivity of dispersion model outputs to the initialising grain size distribution. We compare results from model simulations using the London VAAC default PSD to simulations initialised using varying PSDs which reflect different magma compositions and eruption styles.   Residual ash mass loading was calculated by subtracting the modelled ash concentrations simulated using the VAAC default PSD from the modelled concentrations when the input PSD is varied. Figure 11 shows the results for total column mass loadings on 6 May 2010 12:00 UTC and 8 May 2010 00:00 UTC, for the coarsest (Hekla 1991) and finest (Eyjafjallajökull 2010) PSDs. Values for the remaining PSDs are shown in Figures S5 and S6. Peak values of total column mass loading are up to 5-fold greater using the VAAC default PSD than for all other PSDs, with the most extreme differences observed when comparing the VAAC default to the coarsest PSDs (Eldgja and Hekla) ( Figure 10).  Figure 10 for equivalent non-residual maps); the bottom row is residual deposit mass loadings (on the ground) for the period 4-12 May 2010 (see Figure 12 for equivalent non-residual maps). Table 5 summarises statistical comparisons of the modelled total column mass loadings using the VAAC default PSD compared to simulations initialized with the other PSDs. All PSDs produced lower total column mass loadings than the VAAC default, as shown by the negative FB (−0.218-−1.522). Bias was more extreme when the plume was more dispersed on 8 May. PCC values were high for total column mass loadings, ranging from 0.846 to 0.992, with the highest values seen with the more concentrated plume. The pattern was similar for FoM, where values for total column mass loadings ranged from 38.68 to 90.25 and again the highest values were seen with the more concentrated plume. The modelled total column mass loadings using the finest PSD (Eyjafjallajökull 2010 including satellite observations) to initialise NAME were most similar to the modelled mass loadings using the VAAC default PSD across all measures. Using the coarsest PSD (Hekla 1991) resulted in the biggest differences in the modelled plume, except for KSP values for 6 May, where Eldgja had the highest value. It should be noted that the discrepancies in modelled plumes using different PSDs to initialise the simulation will grow with time and as the plume becomes more dispersed.  Figure 12 shows the modelled deposit mass loading for the period 4-12 May 2010. Using the VAAC default PSD extends the deposit further than when using any other input PSD, but proximal deposit mass loadings are lower (and so show positive residuals in Figure 11 and Figure S7), reflecting the low mass fraction of larger particles (≤ 5 Φ, ≥ 31.25 μm). Proximal mass loadings of these relatively coarse particles are high for all other input PSDs, covering the whole range of magma compositions from basaltic, e.g., Hekla 1991, to rhyolitic, e.g., Rotongaio ( Figure 13).

Deposit Mass Loading
Particles ≤ 5 Φ (≥ 31.25 μm) reached distances up to ~1500 km from source (Figure 13), highlighting a limitation of using remote sensing alone when compiling grain size distributions, as satellite retrievals are usually not configured to identify such large particles [19]. For the smaller particles (> 5 Φ, < 31.25 μm), using the VAAC default PSD results in a more extensive deposit with higher mass loading of ash, and the difference is largest when compared with using the coarsest (mafic/intermediate) PSD (Hekla 1991). In this case, the VAAC PSD is likely to substantially overestimate the extent of the plume and the distance travelled by the small particles before deposition, because there is too much mass on small particles.   (Table 5) show that using the Eyjafjallajökull PSD produced deposits with the lowest bias; the Oruanui PSD produced deposits that were most similar in the other measures, and the Hekla 1991 PSD produced a deposit that was most different. All PSDs resulted in higher deposit mass loadings compared to using the VAAC default PSD, although FB values were smaller (0.068-0.398) than when comparing column mass loadings. PCC values were 0.542 to 0.686, FoM values were between 35.97 and 82.04 and KSP ranged from 9.4 to 38.5.
Binned median values of the modelled deposit grain size, calculated at each grid point ( Figure  14), show that the median modelled grain size is ≤ 5 Φ (≥ 31.25 μm) to ~1400 km from source. The VAAC default deposit most closely matches deposits modelled with the finest input PSDs (Eyjafjallajökull 2010 and Oruanui ~25.4 ka BP).

Grain Size Distributions for Hydromagmatic Eruptions
Grain size distributions of deposits from hydromagmatic eruptions can be relatively constant over large distances from source (Figures 4-6). This is in contrast to magmatic eruptions, where the deposits (and therefore, median grain size) become finer as distance increases (e.g., Figure 6 and Figure S1). The lack of fining with distance in large hydromagmatic deposits suggests that it may be possible to obtain reliable estimates of TGSD from relatively few sampling locations, at least for the largest eruptions. This could be valuable in enlarging the dataset available to dispersion modellers and would aid studies of the fundamental physics of magma fragmentation and controls of eruption intensity. However, the consistency in grain size with distance indicates that the addition of water promotes aggregation and/or en masse settling of particles, which needs to be considered in forecasts of hydromagmatic ash clouds. Both of these processes are promoted by high ash concentrations and so tend to be more important close to source. Accretionary lapilli (compound particles much larger than their component ash particles) found in proximal samples (e.g., [46,51]) can help to explain the relative abundance of fine particles close to source and may account for slight increases in median grain size with distance as seen in Figure 5. Gravitational instabilities in the plume can cause en masse settling of a range of grain sizes including the finest particles, resulting in deposition where grain size is independent of distance [30]. Sedimentation by gravitational instabilities has been observed, for example during eruptions at Eyjafjallajökull in 2010 and Sakurajima in 2015 [17,32].
Questions remain about the extent of water-induced secondary fragmentation of particles in large explosive eruptions, especially 'phreatoplinean' (hydromagmatic Plinian) eruptions. It may be that in the largest hydromagmatic eruptions, water has relatively little effect on the ash size distribution (i.e., the ash is dominantly the result of primary fragmentation) [81]. If this is the case, the hydromagmatic TGSDs could also be useful for estimating TGSD of equivalent magmatic eruptions, particularly where the most distal (fine-grained and thin) portions of the deposits are not mapped and so are not readily included in TGSDs [82,83].

Sensitivity of Modelled Ash Mass Loadings to Input PSDs
Modelled total column and deposit mass loadings show significant differences depending on the input PSD used, with results using the VAAC default PSD being most similar to the finest PSDs tested, which are from silicic eruptions. However, the VAAC default PSD is based on measurements in the air [15,72], and the PSDs compiled in this study, with the exception Eyjafjallajökull 2010, were determined using only samples collected on the ground. Ground-based TGSDs will generally underestimate the mass fraction of the finest material because there is fine ash beyond the mapped deposit (e.g., [83]) and indeed satellite data for Eyjafjallajökull 2010 added significant mass on particles ≥ 7 Φ (≤ 7.8 μm) [33]. This is important because VAAC dispersion modelling aims to quantify distal airborne concentrations of the smallest particles.
Comparison of the results for Eyjafjallajökull 2010 using PSDs defined by (1) the TGSD based on ground-only, and (2) ground-plus-satellite data shows that using ground-based TGSDs produces lower peak ash concentrations and reduces the spatial extent of the plume. This accords with results of other studies highlighting the importance of compiling TGSDs from both ground and airborne measurements (e.g., [84,85]). However, satellite retrieval data can also be problematic as the particles are generally assumed to be dense spheres and the retrieval algorithm is set to interpret brightness temperature difference as coming from particles with effective radius < ~16 μm. This does not take into account that reflection off bubble walls of larger particles would also give the same brightness temperature difference [19]. Further work is required to understand how satellite data can best be used to improve TGSDs of modern historic eruptions. In addition, we have no satellite and aircraft sampling for some types and sizes of eruptions that have not occurred in recent decades (e.g., VEI 7 eruptions) and for these, deposit-only TGSDs are key for assessing suitable PSDs for simulations. Our study focused on particle size, but shape and density are also important factors when considering ash transport and deposition (e.g., [7,9,86]) and the mass fraction held within the finest particles also affects plume and deposit concentrations. For the eruptions selected for our sensitivity tests, we used the London VAAC assumption that particles in the distal ash cloud (which we have taken as ≤ 125 μm) represent 5% of the mass erupted. However the mapped deposits (TGSD) in our entire dataset have mass fractions of particles > 3 Φ (≤ 125 μm) ranging from 3% to 73%. The range for those that are the basis of the PSD for our simulations is 14% to 73%. So including only 5% of total mass erupted in the simulations implies that most of the fine ash is deposited sufficiently close to source that it is not of interest for the applications of the modelling results. There is debate about the mass fraction typically in the distal ash cloud, with comparisons of masses of distal ash based on satellite data to total masses based on ground data, suggesting that the distal ash fraction may be only 0.1% of erupted mass for the largest eruptions [87], but other analyses of the spatial distributions of tephra and fine ash contents of deposits suggest much higher distal fine ash fractions [26,83].

PSDs for Real-Time Forecasting During an Eruption
In this study, we ran simulations with PSDs based on the TGSDs of hydromagmatic and magmatic deposits of a range of compositions. The TGSDs for hydromagmatic eruptions (Figure 1) fall within the range of TGSDs for magmatic eruptions (Figure 2) but are less variable and tend to be finer (greater in Φ scale); however, when we consider only particles likely to be in the distal ash cloud (> 3Φ, < 125 μm), there is almost total overlap of hydromagamatic and magmatic grain size distribution data, with composition (and VEI) being the most important factor (Figure 3). Our NAME simulations using PSDs based on the < 125 μm portion of deposit TGSDs show clear differences between total column mass loadings for mafic/intermediate and silicic eruptions, with modelled concentrations using the current VAAC default PSD more closely matching results for silicic eruptions with the finest PSDs. This could be a problem particularly for forecasting ash plumes from coarse-grained, mafic/intermediate eruptions which occur frequently in Iceland (as shown in Tables 1 and  2).
Model PSDs derived only from ground samples are likely to underestimate the mass of the finest particles and this would need to be considered when choosing any new coarse-grained default PSD. One option would be to explore more sample eruptions where satellite data are also available, as done for Eyjafjallajökull 2010 [33]. Alternatively, a statistical distribution could be used, based on best fit with sampled data (e.g., [59]). Pioli et al. [23] suggest lognormal and Weibull (Rosin-Rammler) distributions are suitable for modelling deposits with fine tails.
The default PSD used by the London VAAC is from aircraft samples, taken from the margin of the ash cloud, during the 1990 eruption of Mount Redoubt and so it represents only a snapshot from that eruption. The ash cloud was generated when a (possibly hydromagmatic) vent explosion and dome collapse created a co-pyroclastic density current (co-PDC) andesitic-dacitic plume [15,72,88]. Because only the finest particles are entrained into co-PDC plumes, the resulting TGSD is likely to be finer than for a vent-derived plume [89].

Conclusions
The grain size distribution of an ash cloud is rarely measured in real time but operational dispersion model simulations used to forecast the location and concentration of ash in the atmosphere are highly sensitive to this parameter. VAACs often rely on using a default PSD to initialise model runs but it has not been clear how best to constrain these.
In this study, we have focused on constraints on PSD from the deposits resulting from ash in the air landing on the ground. We compiled published grain size data from the deposits of hydromagmatic and magmatic eruptions and found that, while grain size distributions for magmatic eruptions change with distance from source, for large hydromagmatic eruptions, grain size remains relatively constant with distance. This suggests that TGSDs for these eruptions could be compiled from fewer samples than are required for typical magmatic eruptions, potentially providing a larger dataset of TGSDs in the future. For the largest hydromagmatic eruptions, if water impacts transport and deposition of particles but has relatively little effect on the ash size distribution, the TGSDs could also provide useful estimates of PSDs for modelling ash dispersion from the equivalent magmatic eruptions.
Although TGSDs for hydromagmatic eruptions tend to be finer than for magmatic eruptions (c.f. Figures 1 and 2), they are similar for the fine portion (< 125 μm) that is likely to be in the distal ash cloud (Figure 3a). Previous studies, focusing on the whole range of grain sizes, have found that TGSD is dominantly controlled by magma composition (and eruption size), with silicic eruptions (> 63% SiO2) typically generating finer TGSDs than eruptions of more mafic magma (e.g., [13,24]). We have found that this relationship also holds true for this fine fraction.
We investigated the sensitivity of modelled ash mass loadings to input PSDs using end-members which represented a coarse PSD from a mafic/intermediate eruption and a fine PSD from a silicic eruption. There are clear differences between mass loadings for forecasts with PSDs based on mafic/intermediate and silicic eruption TGSDs. For the default PSD used by the London VAAC, simulations were most similar to those generated using the finest PSDs based on silicic eruptions and peak values of total column mass loading were up to 5-fold greater than for the coarsest PSD (based on Hekla 1991).
Our input PSDs were based on ground samples, which are likely to lack the most distal material and so adding mass on the finest particles could produce more representative results. Nevertheless, we have constrained end-members representing coarse-and fine-grained PSDs, which could be used by the VAACs to assess the sensitivity of their forecasts to the PSD applied (within the caveats of the limitations associated with the measurement methods).
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Section S1: details of eruptions considered; Section S2: Change in median grain size with distance from source for eruptions with samples for both the hydromagmatic and magmatic phases; Section S3: Comparison of current VAAC default particle size distribution (PSD) with equivalent values on whole-Φ and half-Φ scales; Section S4: Statistical tests and residual mass loadings; Section S5: References to supplementary material. Figure S1 Change of median grain size with distance from source for eruptions having both phreatomagmatic and magmatic phases, Figure S2 VAAC default particle size bins and Φ values plotted on a log scale Φ bins, Figure S3 Log values of Φ and micron scales a) general case. b) example for 8-9 Φ bin, Figure S4 Comparison of NAME output for total column mass loading using VAAC default PSD, whole-Φ scale PSD, half-Φ scale PSD, for times of high and low wind and high and low plume during the Eyjafjallajökull eruption 2010, Figure S5 Residuals for total column mass loading for 6 May 2010 12:00 UTC after removing VAAC default PSD values, Figure S6 Residuals for total column mass loading for 8 May 2010 00:00 UTC after removing VAAC default PSD values, Figure S7 Residuals for deposit mass loadings for 4-12 May 2010 after removing VAAC default PSD values. Table S1 Micron equivalents of whole-Φ particle size scale and corresponding mass fractions for the VAAC default particle size distribution, Table S2 Micron equivalents of half-Φ particle size scale and corresponding mass fractions for the VAAC default particle size distribution, Table S3 Fractional bias of whole-Φ and half-Φ results ( Figure S4) when compared with the VAAC default PSD. Funding: This research received no external funding. interesting conversations on NAME model physics and Jen Saxby for help with NAME and JASMIN. We also thank three reviewers, whose comments have assisted us in substantially improving the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.