Digital Commons @ Michigan Tech Digital Commons @ Michigan Tech Adjusting particle-size distributions to account for aggregation in Adjusting particle-size distributions to account for aggregation in tephra-deposit model forecasts tephra-deposit model forecasts

. Volcanic ash transport and dispersion (VATD) models are used to forecast tephra deposition during volcanic eruptions. Model accuracy is limited by the fact that ﬁne-ash aggregates (clumps into clusters), thus altering pat-terns of deposition. In most models this is accounted for by ad hoc changes to model input, representing ﬁne ash as aggregates with density ρ agg , and a log-normal size distribution with median µ agg and standard deviation σ agg . Optimal values may vary between eruptions. To test the vari-ance, we used the Ash3d tephra model to simulate four deposits: 18 May 1980 Mount St. Helens; 16–17 September 1992 Crater Peak (Mount Spurr); 17 June 1996 Ruapehu; and 23 March 2009 Mount Redoubt. In 192 simulations, we systematically varied µ agg and σ agg , holding ρ agg constant at 600 kg m − 3 . We evaluated the ﬁt using three indices that compare modeled versus measured (1) mass load at sample locations; (2) mass load versus distance along the dispersal axis; and (3) isomass area. For all deposits, under these inputs, the best-ﬁt value of µ agg ranged narrowly between ∼ 2.3 and 2.7 ϕ (0.20–0.15 mm), despite large variations in erupted mass (0.25–50 Tg), plume height (8.5– 25 km), mass fraction of ﬁne ( < 0.063 mm) ash (3–59 %), atmospheric temperature, and water content between these eruptions. This close agreement suggests that aggregation may be treated as a discrete process that is insensitive to eruptive style or magnitude. This result offers the potential for a simple, computationally efﬁcient parameterization scheme for use in operational model forecasts. Further research may this physical on in the evolving cloud.

Abstract. Volcanic ash transport and dispersion (VATD) models are used to forecast tephra deposition during volcanic eruptions. Model accuracy is limited by the fact that fine-ash aggregates (clumps into clusters), thus altering patterns of deposition. In most models this is accounted for by ad hoc changes to model input, representing fine ash as aggregates with density ρ agg , and a log-normal size distribution with median µ agg and standard deviation σ agg . Optimal values may vary between eruptions. To test the variance, we used the Ash3d tephra model to simulate four deposits: 18 May 1980 Mount St. Helens;16-17 September 1992 Crater Peak (Mount Spurr); 17 June 1996 Ruapehu; and 23 March 2009 Mount Redoubt. In 192 simulations, we systematically varied µ agg and σ agg , holding ρ agg constant at 600 kg m −3 . We evaluated the fit using three indices that compare modeled versus measured (1) mass load at sample locations; (2) mass load versus distance along the dispersal axis; and (3) isomass area. For all deposits, under these inputs, the best-fit value of µ agg ranged narrowly between ∼ 2.3 and 2.7ϕ (0.20-0.15 mm), despite large variations in erupted mass (0.25-50 Tg), plume height (8.5-25 km), mass fraction of fine (< 0.063 mm) ash (3-59 %), atmospheric temperature, and water content between these eruptions. This close agreement suggests that aggregation may be treated as a discrete process that is insensitive to eruptive style or magnitude. This result offers the potential for a simple, computationally efficient parameterization scheme for use in operational model forecasts. Further research may indicate whether this narrow range also reflects physical constraints on processes in the evolving cloud.

Introduction
Airborne tephra is the most wide reaching of volcanic hazards. It can extend hundreds to thousands of kilometers from a volcano and impact air quality, transportation, crops, electrical infrastructure, buildings, water supplies, and sewerage. During eruptions, communities want to know whether they may receive tephra and how much might fall. Volcano observatories typically forecast areas at risk by running volcanic ash transport and dispersion (VATD) models. As input, these models require information including eruption start time, plume height, duration, the wind field, and the size distribution of the falling particles. Of these inputs, the particle-size distribution is perhaps the hardest to constrain.
Particle size (along with shape and density) determines settling velocity, which controls where particles land in a given wind field. For different eruptions, the total particlesize distribution (TPSD) can vary. Large eruptions produce more fine ash than small ones for example; and silicic eruptions produce more than mafic . The TPSD is difficult to estimate (e.g., , hence, estimates exist for only a handful of deposits. Even in cases where the TPSD is known, the raw TPSD, entered into a dispersion model, will not accurately calculate the pattern of deposition (Carey, 1996). This inaccuracy results from the fact that complex processes, not considered in models, cause particles to fall out faster than theoretical settling velocities would predict. These processes include scavenging by hydrometeors (Rose et al., 1995a), gravitational instabilities that cause dense clouds to collapse en masse (Carazzo and Jellinek, 2012;Schultz et al., 2006;Durant, 2015;Manzella et al., 2015), and aggregation, in which ash particles smaller than a few hundred microns clump into clusters. The rate of aggregation, as well as the type and size of resulting aggregates, depends on atmospheric processes such as ice accretion, electrostatic attraction, or liquid-water binding, whose importance varies from place to place. Although one VATD model, Fall3d, calculates aggregation during transport for research studies Costa et al., 2010), no operational models consider it. Instead, aggregation is accounted for by either setting a minimum settling velocity in the code (Carey and Sigurdsson, 1982;Hurst and Turner, 1999;Armienti et al., 1988;Macedonio et al., 1988), or, in the model input, adjusting particle-size distribution by replacing some of the fine ash with aggregates of a specified density, shape, and size range (Bonadonna et al., 2002;Cornell et al., 1983;Mastin et al., 2013b). These strategies will probably prevail for at least the next few years, until microphysical algorithms replace them.
These adjustments are mostly derived from a posteriori studies, where model inputs have been adjusted until results match a particular deposit. It is unclear how well the optimal adjustments might vary from case to case. For model forecasts during an eruption, we need some understanding of this variability. This paper addresses this question, using deposits from four well-documented eruptions. We derive a scheme for adjusting TPSD to account for aggregation, optimize parameter values to match each deposit, and then see how much these optimal values vary from one deposit to the next.

Background on the deposits
The IAVCEI Commission on Tephra Hazard Modeling has posted data from eight well-mapped eruption deposits, available for use by modeling groups to validate VATD simulations (http://dbstr.ct.ingv.it/iavcei/). Of these, we focus on eruptions that lasted for hours (not days), where the TPSD included at least a few percent of ash finer than 0.063 mm in diameter, and where data were available from distal (> 35 km) sample locations. Four eruptions met these criteria: the 18 May 1980 eruption of Mount St. Helens, the 16-17 June 1996 eruption of Ruapehu, and the 16-17 September and 18 August 1992 eruptions of Crater Peak (Mount Spurr), Alaska. The August Crater Peak eruption was already studied using Ash3d (Schwaiger et al., 2012)   . (Mastin et al., 2013b), aggregation has been unusually well characterized in recent years (Wallace et al., 2013;Van Eaton et al., 2015). Below are key observations of these events. Deposit maps are shown in Fig. 1, digitized from published sources.  Sarna-Wojcicki et al., 1981;Waitt and Dzurisin, 1981;Rice, 1981). This 9 h eruption expelled magma that was dacitic in bulk composition but contained about 40 % crystals and 60 % rhyolitic glass (Rutherford et al., 1985). The eruption start time (15:32 UTC) and duration are well documented (Foxworthy and Hill, 1982); the time-changing plume height was tracked by Doppler radar (Harris et al., 1981) and satellite (Holasek and Self, 1995) ( Table 2). The deposit was mapped within days, before modification by wind or rainfall, to a distance of ∼ 800 km and to mass load values as low as a few hundredths of a kilogram per square meter (Sarna-Wojcicki et al., 1981). Estimated volume of the fall deposit in dense-rock equivalent (DRE) is 0.2 km 3 (Sarna-Wojcicki et al., 1981) based on what fell in the mapped area. A TPSD was estimated by Carey and Sigurdsson (1982) and later by  to contain about 59 % ash < 63 µm in diameter (Table S1 in the Supplement), with a modal peak in particle size that coincided with the median bubble size of tephra fragments (Genareau et al., 2012). Some fine ash may have been milled in pyroclastic density currents on the afternoon of 18 May and in the lateral blast that morning. A secondary maximum in deposit thickness in Ritzville, Washington (∼ 290 km downwind) was in-ferred by Carey and Sigurdsson (1982) to have resulted from fine-ash aggregating and falling en masse, perhaps as the cloud descended and warmed to above-freezing temperatures . Wind directions that were more southerly at low elevations combined with elutriation off pyroclastic flows in the afternoon to feed low clouds, producing a deposit that was richer in fine ash along its northern boundary than in the south Eychenne et al., 2015). Aggregates sampled by Sorem (1982) in eastern Washington consisted mainly of dry clusters 0.250 to 0.500 mm in diameter, containing particles < 0.001 mm to more than 0.040 mm in diameter, though no aggregates were visible in the fall deposit except at proximal locations (e.g., Sisson, 1995). The eruption began under clear weather conditions. Clouds increased throughout the day. Some precipitation in the form of mud rain was noted within tens of kilometers of the vent (Rosenbaum and Waitt, 1981), probably due to entrainment and condensation of atmospheric moisture in the rising plume. But no precipitation was recorded at more distal locations during the event.
2. The 16-17 September 1991 eruption from Crater Peak, Mount Spurr, Alaska, was the third that summer from this vent. The eruption start time (08:03 UTC, 17 September) and duration (3.6 h; Eichelberger et al., 1995) were seismically constrained. The maximum  (Harris et al., 1981). We calculated a preliminary eruptive volume for each eruptive pulse using the duration and the empirical relationship between plume height and eruption rate (Mastin et al., 2009). This method underestimated the eruptive volume, as noted in previous studies (Carey et al., 1990). Hence, we adjusted the volume of each pulse proportionately so that their total equals the 0.2 km 3 DRE estimated by Sarna-Wojcicki et al. (1981).  (Rose et al., 1995b) increased for the first 2.3 h and then fluctuated between about 11 and 14 km above mean sea level (a.m.s.l.) until the plume height abruptly decreased at 11:10 UTC. The andesitic tephra consisted of two main types -tan and gray, which were both noteworthy for their low vesicularity (∼ 20-45 %) and high crystallinity (40-100 %) (Gardner et al., 1998). The deposit was mapped rapidly after the eruption (Neal et al., 1995;McGimsey et al., 2001) to a distance of 380 km and mass loads as low as 0.050 kg m −2 . This deposit displays a weak secondary thickness maximum 260-330 km downwind.  derived a TPSD for this deposit, estimating about 40 % smaller than 0.063 mm. Milling in proximal pyroclastic flows that accompanied this eruption (Eichelberger et al., 1995) could have contributed fine ash. The eruption occurred at night under clear skies (Neal et al., 1995).
3. The 17 June 1996 eruption of Ruapehu produced a classic weak plume that was modeled by , Hurst and Turner (1999), Scollo et al. (2008), Liu et al. (2015), and Klawonn et al. (2014), among oth-ers. The main phase involved two pulses, one beginning 16 June at 19:10 UTC and lasting 2.5 h and the second at 23:00 UTC and lasting approximately 1.5 to 2 h. Ash-laden plumes reached to about 8.5 km a.m.s.l. based on satellite infrared images (Prata and Grant, 2001). The deposit was mapped out to the Bay of Plenty (190 km), sampled at 118 locations to mass loads less than 0.01 kg m −2 , and yielded a total mass of about 0.001 km 3 DRE . Ejecta consisted mainly of scoria containing 75 % glass and 25 % crystals, with glass containing about 54 wt % SiO 2 (Nakagawa et al., 1999). A TPSD estimate based on the Voronoi tessellation method  suggested that ash < 0.063 mm composed only about 3 % of the deposit. A minor secondary thickness maximum was constrained by mapping at about 160 km downwind   (Fig. 1c). Although some witnesses at distal locations observed loose, millimeter-sized clusters falling, no aggregates or accretionary lapilli were present in the deposit (Klawonn et al., 2014). The eruption was not accompanied by significant pyroclastic density currents and occurred during clear weather.
4. Event 5 of the 23 March 2009 eruption of Redoubt Volcano, Alaska, erupted through a glacier and entrained a variable amount of water into a high-latitude earlyspring atmosphere. It began at 12:30 UTC, lasted about 20 min on the seismic record (Buurman et al., 2013), and sent a plume briefly to about 18 km as seen in both National Weather Service NEXRAD Doppler radar from Anchorage, and a USGS mobile C-band radar system in Kenai, Alaska (Schneider and Hoblitt, 2013). Within a few days after the eruption, the deposit was mapped by its contrast with underlying snow in satellite images (NASA MODIS), and sampled for mass load and particle-size distribution at 38 locations, at distances up to ∼ 250 km and mass loads as low as 0.01 kg m −2 (Wallace et al., 2013). During Ash3d modeling of this eruption, Mastin et al. (2013b) found that wind vectors varied rapidly with both altitude and time, making the dispersal direction highly sensitive to both the plume height (which varied from ∼ 12 to 18 km during the 20 min eruption) and the vertical distribution of mass in the plume. In the deposit, Wallace et al. (2013) described abundant frozen aggregates with size decreasing with distance from the vent, from about 10 mm at 12 km distance. Schneider et al. (2013) attributed the high (> 50 dBZ) reflectivity of the proximal plume in radar images, and a rapid decrease in maximum plume height over a period of minutes, to formation and fallout of ashy hail hydrometeors in the rising column. Van Eaton et al. (2015) combined analysis of the aggregate microstructures with a three-dimensional (3-D) largeeddy simulation to show that the ash aggregates grew directly within the volcanic plume from a combination of wet growth and freezing, in a process similar to hail formation.

The Ash3d model
We model these eruptions using Ash3d (Schwaiger et al., 2012;Mastin et al., 2013a), an Eulerian model that calculates tephra transport and deposition through a 3-D, time-changing wind field. Ash3d calculates transport by setting up a 3-D grid of cells, adding tephra into the column of source cells above the volcano, and distributing the mass in the column following the probability density function of Suzuki (Suzuki, 1983), modified by Armienti et al. (1988): where Q m is the mass eruption rate, H v is plume height above the vent, z is elevation (above the vent) within the plume, and k is a constant that adjusts the mass distribution. Suzuki (1983) defines this function as a "probability density of diffusion" of mass from the column as particles fall out.
Here we regard it as a simplified parameterization of mass distribution with no implication for physical process. At each time step, tephra transport is calculated through advection by wind, through turbulent diffusion, and through particle settling. For wind advection, simulations of Mount St. Helens, Crater Peak, and Redoubt use a wind field obtained from the National Oceanic and Atmospheric Administration (NOAA) NCEP/NCAR Reanalysis 1 model (RE1) (Kalnay et al., 1996). For the Ruapehu simulations we used a local 1-D wind sounding, which gave more accurate results as detailed below. The RE1 model provides wind vectors on a global 3-D grid spaced at 2.5 • latitude and 2.5 • longitude, and 17 pressure levels in the atmosphere (1000-10 hPa), updated at 6 h intervals. Ash3d calculates turbulent diffusion using a specified diffusivity D (Schwaiger et al., 2012, Eq. 4). D is set to zero for simplicity, though later we show the effect of different values of D.
Settling rates are calculated using relations of Wilson and Huang (1979) for ellipsoidal particles. Wilson and Huang define a particle shape factor ≡ F (b + c)/2a, where a, b, and c are the maximum, intermediate, and minimum diameters of the ellipsoid, respectively. Wilson and Huang measured a, b, and c for 155 natural pyroclasts. From data published in Wilson and Huang, we calculate an average F of 0.44, which we use in our model. For aggregates we use F = 1.0 (round aggregates).
Other model inputs include the extent and nodal spacing of the model domain; vent location and elevation; the eruption start time, duration, plume height, erupted volume, diffusion coefficient D, and a series of particle-size classes and associated densities. The size classes may represent either individual particles or aggregates. These input values are given in Tables 1 and 2. 3.2 Adjusting particle-size distributions to account for aggregation In deriving a particle-size adjustment scheme we found it necessary to prioritize the type(s) of processes and products we wish to replicate. The rate and type of ash aggregation are known to vary with both eruptive conditions and meteorology. Large aggregates, including frozen accretionary lapilli, form near the source and are abundant in phreatomagmatic deposits Brown et al., 2012;Houghton et al., 2015). They are associated with particles colliding in moist, turbulent updrafts within a rising plume ( Fig. 2) or an elutriating ash cloud. These near-source aggregates commonly exceed 1 cm diameter (Wallace et al., 2013;Swanson et al., 2014;Van Eaton and Wilson, 2013). In contrast, the low-density aggregates that produced the Ritzville Bulge, 230 km downwind from Mount St. Helens, are thought to have been triggered by mammatus cloud instabilities ) as the cloud descended, warmed, and ice melted into liquid water (red line, Fig. 2). These aggregates tend to be smaller than a millimeter, and form in the cloud hundreds of kilometers downwind from the source (Sorem, 1982;Dartayat, 1932). At Mount St. Helens and perhaps other places, investigators found evidence for both large, wet, proximal accretionary lapilli (Sisson, 1995) and distal, dry aggregates (Sorem, 1982). The latter type deposited over a larger area, involved a greater fraction of the total erupted mass, and affected a greater population. Thus, it is the latter process whose deposits we wish to reproduce. Aggregation is also a highly size-selective process. The threshold size below which most particles aggregate and above which they do not varies with moisture and electrical charge, ranging from several tens of microns under dry conditions, to hundreds of microns when liquid water is present (Gilbert and Lane, 1994;Schumacher and Schmincke, 1995;Van Eaton et al., 2012). Our aggregation scheme is too crude to distinguish the threshold size as a function of atmospheric conditions; hence, we use a broad range such that for φ > = 4, all ash aggregates; for φ < = 2, no ash aggregates. For 4> ϕ > 2, the mass fraction that aggregates varies linearly with ϕ from 1 (when φ = 4) to 0 (when ϕ = 2).
Path taken by aggregates that form in the cloud    Figure 2. Illustration of the path taken by coarse aggregates that fallout in proximal sections, less than a few plume heights from the source (left), and fine aggregates that fall out in distal sections (right). Among distal fine aggregates, we show the path taken by those that might have formed within or below the downwind cloud as hypothesized by  (red dashed line), and those that were transported downwind without changing size, as calculated by Ash3d (blue dashed line). Also illustrated are some key processes that might influence the distribution of fine, distal ash, including development of gravitational instability and overturn within the downwind cloud (Carazzo and Jellinek, 2012), and the development of hydrometeors as descending ash approaches the freezing elevation ). The TPSDs used to model these four eruptions are listed in Table S1 and illustrated as gray bars in Fig. 3. Particle sizes that do not aggregate according to this scheme are illustrated as black bars. We assume that the aggregates collect into clusters having a Gaussian size distribution of mean µ agg , and standard deviation σ agg (insets, Fig. 3). For deposit modeling, we ignore the small fraction of the erupted mass that goes into the distal cloud, typically a few percent (Dacre et al., 2011;Devenish et al., 2012).
In our study, the aggregated ash mostly deposits as a secondary thickness maximum. Different choices of a threshold size for particle aggregation would influence the mass building the secondary maximum. For Mount St. Helens, about 10 % of the erupted mass lies between ϕ = 2 and ϕ = 4. For Spurr, Ruapehu, and Redoubt, the percentages are 28, 6, and 11 %. These values reflect the variability in mass of the secondary maximum that could result from different choices of the aggregation-size threshold.

Aggregate density: different processes influence aggregate density
Wet ash (> 10-15 wt % liquid water) rapidly produces subspherical pellets with density > 1000 kg m −3 (Schumacher and Schmincke, 1991;Van Eaton et al., 2012); drier conditions lead to electrostatically bound clusters (Schumacher Table 3. Statistical measures of fit used in this paper.

Name Formula Explanation
Point-by-point method The mass load m o,i observed at each sample location i is compared with modeled mass load m m,i at the same location. Squared differences are summed to the total number of sample points N, and normalized to the sum of squares of the observed mass loads.
Downwind thinning method The log of modeled mass load m m,j at a point j on the dispersal axis is compared with the observation-based value m o,j expected at that location based on a trend line drawn between field measurements along the axis (Fig. 4). Differences between m m,j and m o,j are calculated on a log scale, squared, and summed.

Isomass area method
This method calculates the area A m,i of the modeled deposit that exceeds a given mass load i by summing the area of all model nodes that meet this criterion. It then takes the difference between A m,i and the area A o,i within same isomass line mapped from field observations. The sum of the squares of these differences, normalized to the sum of the squared mapped isopach areas, gives the index 2 area .
and Schmincke, 1995;Van Eaton et al., 2012) with density in the hundreds of kilograms per cubic meter range (James et al., 2002;Taddeucci et al., 2011). Taddeucci et al. (2011) estimated densities ranging from < 100 to > 1000 kg m −3 in dry aggregates photographed falling 7 km from the Eyjafjallajökull vent. For simplicity, we hold ρ agg constant at 600 kg m −3 , toward the middle of the observed range but higher than that of some dry aggregates. Optimal aggregate sizes that we derive later in this paper are determined by this assumed density, and may be larger or smaller than actual aggregate sizes.

Statistical measures of fit
For each eruption, we have done a series of model simulations, first using the TPSD without considering aggregation, and then systematically varying σ agg and µ agg to include the effects of aggregation. We compare the resulting modeled deposit with the mapped deposit using three methods presented in Table 3. Each has advantages and disadvantages.
1. The point-by-point index 2 compares model results with sample data collected at specific locations (dots, Fig. 1). It offers the advantage that the comparison is made directly with measured values, not with interpreted or extrapolated contours of data. But 2 can be influenced by errors in the wind field, which cannot be adjusted in the model. More importantly, 2 can be dominated by differences in proximal locations where mass per unit area is greatest, and where near-vent processes, such as fallout from the vertical column, are not accurately simulated. For these reasons, we exclude proximal data, within a few column heights distance from the vent, from the calculation of 2 .
2. The downwind thinning index 2 downwind compares modeled mass per unit area along the downwind dispersal axis with values expected at that distance based on a trend line drawn from field measurements (Fig. 4). The comparison is not made directly with measured values (a disadvantage). However, the method does not suffer the limitation of over-weighting proximal data, and, more importantly, it still provides a useful comparison when wind errors cause the modeled dispersal axis to diverge from the mapped one.
3. The isomass area index 2 area compares the area within modeled and mapped isomass lines. It is based on traditional plots of the log of isopach thickness versus square root of area (Pyle, 1989;Fierstein and Nathenson, 1992;Bonadonna and Costa, 2012), which are assumed to accurately depict the areal distribution of tephra while minimizing the effects of 3-D wind on the distribution (Pyle, 1989). Figure 5 shows plots for our four eruptions, using the log of isomass rather than isopach thick-ness to avoid problems introduced by varying deposit density.
The index 2 area is assumed to be insensitive to effects of wind (an advantage). However, model results are compared with isopach lines that are interpretive and may not be well constrained, depending on the distribution and number density of sample locations.

Sensitivity to various input values
We ignore complex, proximal fallout and concentrate on medial to distal areas, about 100 to ∼ 500 km downwind at Mount St. Helens, for example. There, under the average wind speed (15.1 m s −1 ) that existed below about 15 km, tephra falling from 15 km at average settling velocities of 0.4-1.5 m s −1 would deposit within this range (Fig. 6a). Tephra falling at 0.66-0.78 m s −1 would land 290-340 km downwind, the distance of the secondary maximum at Ritzville. A wide range of aggregate diameters d could fall at this rate depending on density ρ agg (Fig. 6b).
Other factors listed below can also affect the results.
-Aggregate shape: aggregate shape can strongly affect the settling velocity and thus where deposits fall, as illustrated in Fig. 7. For simplicity, we use round aggregates (F = 1.0).
-Suzuki k: simulations of Mount St. Helens (Fig. 8) show that increasing the Suzuki factor from 4 to 8 increases the prominence of a secondary thickness maximum. But at k > ∼ 8, the proximal deposit becomes unrealistically thin. Our simulations use k = 8 to replicate the known prominent secondary thickening while minimizing unrealistic thinning of proximal deposits.
-Aggregate size: the transport distance is highly sensitive to aggregate size. Reducing aggregate diameter d from 0.250 to 0.217 to 0.189 mm increases transport distance at Mount St. Helens from 300 to 366 to 448 km, respectively (Fig. 6a). In simulations that use a single, dominant aggregate size, these variations produce conspicuous changes in the location of a secondary maximum (Fig. 9). Decreasing size also decreases the percent of erupted mass that lands in the area shown in Fig. 9: from 63 to 35 to 15 % for d = 0.165, 0.143, and 0.125 mm, respectively (ϕ = 2.6, 2.8.3.0). At d = 0.1 mm (ϕ = 3.3), only 4 % of the erupted mass lands in the mapped area.
This constrains the range of aggregate sizes we may use in our simulations. Sparse observations suggest that > 90 % of erupted mass falls as an observable deposit while less than several percent is transported downwind as a distal cloud (Wen and Rose, 1994;Devenish et al., 2012). To ensure a similar relationship in our simulations, nearly all of the aggregate-size distribution must be coarser than about 0.1 mm. At the proximal end, for Mount St. Helens,  found that most fine ash fell at distances were then adjusted proportionally so that their sum added to 1.
-Fall-velocity model: different fall-velocity models are used in different tephra dispersion models. These models give slightly different results, and it should be noted that our results are specific to our choice of the Wilson and Huang fall model.
Finally, we note that key parameters, such as particle density, shape, Suzuki k, etc., are held constant for all four eruptions even though they may vary from one eruption to another. Such parameters cannot easily be scrutinized when setting up simulations during an eruption. An objective is to see how well "standard" values, even if locally unrealistic, can reproduce observations.

Results
We ran simulations at µ agg = 1.9, 2.0, 2.1. . .3.1ϕ, and σ agg 0.0, 0.1, 0.2, and 0.3ϕ. The latter used 1, 5, 7, and 11 aggregate size classes, respectively, in each simulation, with the percentage of fine ash assigned to each bin given in Table 4. Our calculations of 2 and 2 downwind only included sample points, whose downwind distance lies within the range indicated by the trend lines in Fig. 4. Figure 10 shows contours of 2 , 2 downwind , and 2 area as a function of σ agg and µ agg for each of these four deposits.  Values are given in Tables S3-S6. Although the three indices compare different features of the deposit, they provide roughly similar optimal values of µ agg . For Mount St. Helens, for example, the best-fit value of µ agg is about 2.4ϕ using 2 (Fig. 10a), 2.5ϕ using 2 downwind (Fig. 10b), and 2.7ϕ using 2 area (Fig. 10c). Optimal values of σ agg are 0.1, 0.1, and 0.2, respectively. For Crater Peak, optimal µ agg values are 2.6ϕ, 2.5ϕ, and 2.0ϕ, respectively, while for Ruapehu they are 2.3ϕ, 2.5ϕ, and 2.5ϕ. For both Crater Peak and Ruapehu, optimal values of σ agg range from 0.0 to 0.2. For Redoubt, optimal values are disparate: µ agg = 2.5ϕ, 2.5ϕ, and <2ϕ, respectively. The Redoubt deposit is the least constrained by field data and the most difficult to match due to the complex wind conditions. Figures 11-14 show results for each of these eruptions using µ agg = 2.4ϕ (0.19 mm) and σ agg = 0.1ϕ. The sizes of particles and aggregates used to generate these figures is given in Table S2. For all deposits these values are close to optimal, depending on which criterion is used. Similar figures for other values of µ agg and σ agg are provided as Figs. S005-S212. Figures S001-S004 show simulations using the original particle-size distribution, with no aggregation. Tephra fall beyond a few tens of kilometers is strongly underestimated in all these runs, especially for the three eruptions that contain more than a few percent fine ash. Values of 2 , 2 downwind , and 2 area are also higher than most simulations that use aggregates (Tables S3-S6). For Mount St. Helens, Crater Peak, Ruapehu, and Redoubt, the percentages of the erupted mass landing in the mapped area are very low: 29, 42, 88, and 59 %, respectively.
Optimal aggregates obtained from our study are similar in size but denser than those found optimal by Cornell et al. (1983) for the Campanian Y-5 (µ agg = 2.3ϕ,

Mount St. Helens
For the Mount St. Helens case, the modeled deposit follows a dispersal axis (solid black line, Fig. 11a) that matches almost exactly with the mapped one (dashed line). The agreement reflects both the faithfulness of the numerical wind field to the true one and the appropriateness of other inputs, such as k, that influence dispersal direction. The measured mass loads in Fig. 11a, indicated by the color of markers, agree reasonably well with modeled mass loads indicated by col-     Figure 11b shows that modeled and measured mass loads generally agree within a factor of 3 or so, except for those same distal, low-mass-load measurements, to the lower left of the legend label (those where modeled values are truly zero do not show up on this plot). Figure 11c shows that the modeled mass load (black line with dots) contains a secondary thickening at about the same location mapped (dashed line). It also has roughly the same downwind shape, in contrast to results using σ agg = 0.2 and 0.3 (Figs. S027-S028), in which the secondary thickening is broader and thinner than observed. However, the modeled mass load is consistently less than measured, especially at the most distal sites. In Fig. 11d, the log of modeled mass load versus square root of area shows reasonable agreement with mapped values until mass loads are less than about 1 kg m −2 , where they diverge.
Notably, modeled mass loads somewhat underestimate the measured values along the dispersal axis in Fig. 11c. The underestimate reflects the fact that the input erupted volume of 0.2 km 3 DRE (Table 1) was based on estimates by Sarna-Wojcicki et al. (1981), which lies within the mapped area in Fig. 11a, yet only about 78 % of the modeled mass landed within this area. Reducing the mean aggregate size to 2.6ϕ (0.164 mm, Fig. S036) improves the fit somewhat along distal parts of the transect but degrades it near Ritzville; the finer size moves the secondary maximum too far east and reduces the percentage deposited to ∼ 65 %.
In Fig. 11a, the modeled deposit is also slightly narrower than the mapped one. Adding turbulent diffusion, with a diffusivity D of about 3 × 10 2 m 2 s −1 (Fig. 15) visually improves the fit, and was likely important during this eruption due to high crosswind speeds that increased entrainment (Degruyter and Bonadonna, 2012;Mastin, 2014). But adding diffusion slightly increases 2 , improving fit on deposit margins at the expense of the axis. Ignoring turbulent diffusion also decreases run time by ∼ 3×, from ∼ 30 to 10 min, yielding faster results under operational conditions. Results with other models may vary depending on model setup and configuration.

Crater Peak (Mount Spurr)
At Crater Peak (Mount Spurr), results in Fig. 12a also show good agreement between the modeled dispersal axis and the mapped one (which is constrained by fewer sample locations than the Mount St. Helens case). The isomass lines in this plot are jagged and irregular due to effects of topography in this mountainous region. The modeled location of secondary thickening in Fig. 12c agrees with the mapped location, about 250-300 km downwind. Although Fig. 12c shows a tendency for the model to underestimate the mass load along the dispersal axis, there is less tendency to underestimate the mass load in the most distal locations as has occurred at Mount St. Helens. In Fig. 12d, the areas covered by modeled isomass lines are comparable to the mapped values, down to mass loads approaching 0.1 kg m −2 .

Ruapehu
For Ruapehu (Fig. 13), simulations using the NCEP Reanalysis 1 numerical winds produced an odd double dispersal axis, whose average did not correspond well with the mapped direction of dispersal (Fig. 1c). To improve the fit we used the 1-D wind sounding provided for this eruption at the IAVCEI Tephra Hazard Modeling Commission web page (http://dbstr.ct.ingv.it/iavcei/). Use of a 1-D wind sounding seems justified in this case because this deposit covers a smaller area than the others, making a 3-D wind field less important in calculating transport. The resulting dispersal axis (Fig. 13a) agrees with the mapped one out to about 140 km distance, beyond which it strays eastward, reaching the coast, 180 km downwind, about 10 km east of the mapped axis. This slight difference is enough to cause misfits in point-topoint comparisons at measured mass loads of ∼ 10 −1 kg m −2 (Fig. 13b).
The modeled mass load along the dispersal axis (Fig. 13c) agrees with measurements to about 60-90 km distance. At A large discrepancy is also apparent at distances of less than 60 km, where mass load along the dispersal axis (Fig. 13c) and the area covered by thick isomass lines (Fig. 13d) are greater than that for the mapped deposit. The implication is that too much mass is dropping out proximally in the model. Underestimates of isomass area at greater than or equal to 10 −1 kg m −2 (Fig. 13d) also show that too little is falling distally. Simulations (not shown) that raise the plume height or increase k to concentrate more mass high in the plume do not improve the fit. The discrepancy may reflect the coarse TPSD -50 % of which is coarser than 1 mm (compared with 2, 12, and 8 % for the other three deposits in Table S1). An additional simulation used the TPSD derived from technique B of   (Table S1), which divides the deposit into arbitrary sectors, and calculates a weighted sum of the size distributions in each sector following Carey and Sigurdsson (1982). Tech-nique B yields a finer average particle size than technique C, which uses Voronoi tessellation to sectorize the deposit. But the finer particle size of the technique B TPSD does not improve the fit. Further exploration of this discrepancy is beyond the scope of this paper, but other possible causes could include release of different particle sizes at different elevations, or complex transport in the bending of the weak plume that cannot be accommodated in this model.
A second, smaller discrepancy is that the modeled deposit is narrower than the mapped one (Fig. 1c). As at Mount St. Helens, deposit widening due to cross-flow entrainment is likely. Increases in entrainment resulting from cross flow is widely known to both increase plume width and decrease its height for a given eruption rate (Briggs, 1984;Hoult and Weil, 1972;Hewett et al., 1971;Woodhouse et al., 2013). Adding turbulent diffusion, we get a visually improved fit when D =∼ 3 × 10 2 m 2 s −1 (Fig. 16), consistent with findings by  based on the rate of downwind widening of isomass lines. This diffusivity is also similar to the visual best-fit value for Mount St. Helens (Fig. 15).   Figure 13. Results of the Ruapehu simulation that provide a good best fit to mapped data (µ agg = 2.4ϕ and σ agg = 0.1ϕ). The features in the sub-figures are as described in Fig. 11.
Despite the uncertainty in TPSD, simulations that systematically vary µ agg and σ agg fit best in Fig. 10g, h, and i when µ agg is about 2.3 to 2.5. Results similar to those presented in Fig. 13c use other values of µ agg (Figs. S109-S160) and show a secondary maximum migrating downwind as µ agg increases, coming into agreement with the mapped distance at µ agg = 2.3 to 2.5ϕ (0.20-0.18 mm), where errors in Fig. 10g, h, and i are the lowest.

Redoubt
This deposit is the second smallest in our group, the least well constrained by sampling, and the only one in our group not known to include a secondary thickness maximum. Mastin et al. (2013b) modeled this deposit using numerical winds from the North American Regional Reanalysis model (Mesinger et al., 2006). During that eruption, the winds at 0-4 km, 6-10, and > 10 km elevation were directed toward the northwest, north, and northeast, respectively, with the highest speeds at 6-10 km. Mastin et al. (2013b) found that the modeled cloud developed a north-oriented, northward-migrating wishbone shape with the west prong at low elevation and the east prong at high elevation. Mastin et al. (2013b) also found that the modeled dispersal axis and the mass load distribution roughly agreed with mapped values for a plume height of 15 km, k = 8, and a particle-size adjustment that involved   Figure 14. Results of the Redoubt simulation that provide a reasonable fit to mapped data (µ agg = 2.4ϕ and σ agg = 0.1ϕ). The features in the sub-figures are as described in Fig. 11.
taking 95 % of the fine ash (< 0.063 mm) and distributing it evenly among the coarser bins. In this study we use the same plume height and k value, a different wind field (RE1), and explore a different parameterization for particle aggregation. In Fig. 14a, the modeled dispersal axis diverges about 20 • westward from the mapped axis. We do not correct this divergence by adjusting mass height distribution, since the optimal values of µ agg and σ agg can still be obtained from 2 downwind and 2 area . As with the Crater Peak (Spurr) simulations, the isomass lines are jagged and patchy, an artifact of high relief. (The most distal sample location lies at 4.3 km elevation on the west shoulder of Mount Denali.) Although the value of µ agg (2.4ϕ, 0.19 mm) portrayed in Fig. 14 is close to optimal in Fig. 10j, many sample points do not plot in Fig. 14b because the modeled mass load is zero, and most values of 2 are high (0.99) largely because of the disparity in axis dispersal directions and the consequent fact that sample points lie outside the modeled deposit. The reason that 2 shows a clear minimum, around µ agg = 2.4ϕ (0.19 mm) in Fig. 10j, is apparent from Figs. S161-S212, which show that, as µ agg decreases in size, the modeled deposit extends farther north and takes a clear turn to the northeast, overlapping more with the mapped deposit. These figures also illuminate why 2 downwind is optimal at µ agg = 2.3; i.e., modeled and mapped loads come into best agreement along the dispersal axis for aggregates of this size. 2 area is optimized at µ agg < 2 because the area of the 1 kg m −2 isomass diverges below the mapped value, and the area of the 0.01 kg m −2 iso-  Helens, the freezing elevation was also checked using data from the North American Regional Reanalysis model (Mesinger et al., 2006), available at the same NOAA site, and found to be 3.3 km, similar to that given below by the RE1 model.

Mount St. Helens Crater Peak (Spurr) Ruapehu
Redoubt  Figure 16. Modeled mass load of the Ruapehu eruption for four cases using µ agg = 2.4ϕ, σ agg = 0.1ϕ, and different diffusion coefficients: (a) D = 0 m 2 s −1 , (b) 1× 10 2 m 2 s −1 , (c) 3 × 10 2 m 2 s −1 , and (d) 1 × 10 3 m 2 s −1 . Other inputs are as given in Table 1. Lines are isomass contours of modeled mass load and colored dots are sample locations. Colors of the dots and lines give the mass load corresponding to the color table.
on other inputs, such as particle density, shape factor, and Suzuki factor. Values assigned here may not always be representative. Aggregate density for example is frequently less than 600 kg m −3 ; different assumptions on particle or aggregate shape could significantly change our results. Moreover, our result is partly an artifact of our choice to optimize fit to deposits at medial distances of several tens to hundreds of kilometers. Including more proximal sample points may have given optimal aggregate sizes that spanned a wider range, as used for example in aggregation schemes for Vesuvius (Bar-sotti et al., 2015) or Iceland (Biass et al., 2014). Despite these considerations, the similarity in optimal values of µ agg between these four eruptions is noteworthy.
The overall agreement in modeled mean aggregate size (µ agg ) suggests that accelerated fine-ash deposition may be treated as a discrete process, insensitive to eruptive style or magnitude. It seems unlikely that these varied eruptions would produce aggregates of the same size, density, and morphology. A combination of processes removed ash, and our approach captures these processes implicitly, ignoring the microphysics.
What sort of processes could evolve in the cloud? Some possibilities are illustrated in Fig. 2. The evolution starts with ejection of particles from the vent, with size ranging from microns to meters. For an eruption having the TPSD of Mount St. Helens, the rising plume would have contained 10 6 -10 8 particles per cubic meter with diameter between 10 and 30 µm that collided with larger particles many times per second. High collision rates and the availability of liquid water in the plume would have led to rapid aggregation. Freezing of liquid water and riming would have shifted the maximum possible size of aggregates towards millimeter to centimeter sizes. Mud rain, observed falling at Mount St. Helens (Waitt, 1981), and ice aggregates collected near the vent at Redoubt , are evidence of these processes.
In the downwind cloud particle concentrations were lower, turbulence was less intense, a smaller range of particle sizes existed, and, for all four eruptions, atmospheric temperatures near the plume top were well below freezing (Table 5), leading to presumably slow aggregation rates. However, at least two other processes may help settle ash from downwind clouds. One is gravitational overturn. Experiments (Carazzo and Jellinek, 2012) have observed that fine ash settles toward the bottom of ash clouds as they expand and move downwind, accumulating gravitationally unstable particle boundary layers that eventually overturn and cause the entire air mass to settle rapidly. At Eyjafjallajökull in 2010, gravitational convective instabilities formed within 10 km of the vent, presumably as a result of accumulation of coarse ash over a period of minutes (Manzella et al., 2015). The development of fine-ash particle boundary layers presumably takes longer, perhaps hours, although the underlying processes remain a subject of active research.
A second process is hydrometeor growth. In some cases, magmatic and (or) externally derived water in the eruption cloud may condense on ash particles and initiate hydrometeor growth. Both hydrometeor growth and gravitational overturn have been suggested to produce the mammatus clouds that developed in mid-day over central Washington on 18 May 1980 and signaled mass settling (Durant, 2015;Carazzo and Jellinek, 2012). Mammatus descent rates are typically meters per second (Schultz et al., 2006), much faster than the settling rate of individual ash particles (< 0.1 m s −1 ) or even of ash aggregates (< ∼ 1 m s −1 , Fig. 6).
The extent to which these processes operated at Crater Peak, Ruapehu, and Redoubt is unknown. Cloud structures were not observed during the nighttime eruptions of Redoubt and Crater Peak (Spurr). Although virga-like structures can be seen in some near-vent photos of Ruapehu (Bonadonna et al., 2005, Fig . 9a), we have seen no documentation of such instabilities farther downwind.
For operational forecasting, these mechanisms cannot be considered in any case, because no operational model has the capability to resolve these processes. The fact that these eruptions can all be reasonably modeled using similar inputs for aggregate size is convenient, even if the processes involved are not specified in the model. The agreement suggests that model forecasts can still be useful during the coming years. Future work will focus on the development of more sophisticated algorithms that account for cloud microphysics.
The Supplement related to this article is available online at doi:10.5194/acp-16-9399-2016-supplement.
Author contributions. L. Mastin conceived the study, did the model simulations, and wrote most of the paper. A. Van Eaton provided advice on aggregation processes. A. Durant provided the data for Mount St. Helens and Crater Peak, and advice on aggregation processes that occurred during those two eruptions.