Adjusting particle-size distributions to account for 1 aggregation in tephra-deposit model forecasts 2 3

3 Larry G. Mastin1, Alexa R. Van Eaton1, and Adam J. Durant2,3 4 [1] [U.S. Geological Survey, Cascades Volcano Observatory, 1300 SE Cardinal Court, Bldg. 5 10, Suite 100, Vancouver, Washington, USA] 6 [2] [Section for Meteorology and Oceanography, Department of Geosciences, University of 7 Oslo, Blindern, 0316 Oslo, Norway] 8 [3] [Geological and Mining Engineering and Sciences, Michigan Technological University, 9 1400 Townsend Drive, Houghton, MI 49931, USA] 10 11 Abstract 12 Volcanic ash transport and dispersion models (VATDs) are used to forecast tephra deposition 13 during volcanic eruptions. Model accuracy is limited by the fact that fine ash aggregates, 14 altering patterns of deposition. In most models this is accounted for by ad hoc changes to model 15 input, representing fine ash as aggregates with density agg ρ , and a log-normal size distribution 16


177
where Qm is the mass eruption rate, Hv is plume height above the vent, z is elevation (above the 178 vent) within the plume, and k is a constant that adjusts the mass distribution. 179 At each time step, tephra transport is calculated through advection by wind, through turbulent 180 diffusion, and through particle settling. For wind advection, simulations of Mount St. Helens,181 Crater Peak, and Redoubt use a wind field obtained from the National Oceanic and Atmospheric 182 Administration's (NOAA's) NCEP/NCAR Reanalysis 1 model ("RE1") (Kalnay et al., 1996).

Adjusting particle size distributions to account for aggregation 199
The TPSD used to model these four eruptions are listed in Table S1 and illustrated in Fig. 2. 200 We aim to adjust the TPSD in our model to better match the mapped deposits. In doing so, we 201 assume that some fraction (magg) of ash smaller than some size

Statistical measures of fit 212
For each eruption, we have done a series of model simulations, first using the TPSD without 213 considering aggregation, and then systematically varying agg σ and agg µ to include the effects of 214 aggregation. We compare the resulting deposit with the mapped deposit using three methods 215 presented in Table 3. Each has advantages and disadvantages. 216 1) The point-by-point index ∆ 2 compares model results with sample data collected at specific 217 locations (dots, Fig. 1). It offers the advantage that the comparison is made directly with 218 measured values, not with interpreted or extrapolated contours of data. But ∆ 2 values are 219 dominated by differences in proximal locations where mass per unit area is greatest; and values 220 of ∆ 2 can be influenced by errors in the wind field, which cannot be adjusted in the model. 221 2) The downwind thinning index 2 downwind ∆ , compares modeled mass per unit area along the 222 downwind dispersal axis with values expected at that distance based on a trend line drawn from 223 field measurements (Fig. 3). The comparison is not made directly with measured values (a 224 disadvantage). However the method does not suffer the limitation of over-weighting proximal 225 data. And, more importantly, it still provides a useful comparison when wind errors cause the 226 modeled dispersal axis to diverge from the mapped one. 227 3) The isomass area index 2 area ∆ compares the area within modeled and mapped isomass 228 lines. It is based on traditional plots of the log of isopach thickness versus square root of area 229 (Pyle, 1989;Fierstein and Nathenson, 1992; Bonadonna and Costa, 2012), which are assumed 230 to accurately depict the areal distribution of tephra while minimizing the effects of 3-D wind 231 on the distribution (Pyle, 1989). Fig. 4 shows plots for our four eruptions, using the log of 232 isomass rather than isopach thickness to avoid problems introduced by varying deposit density. 233 The index 2 area ∆ is assumed to be insensitive to effects of wind (an advantage). However, 234 model results are compared with isopach lines that are interpretive and may not be well 235 constrained, depending on the distribution and number density of sample locations. 236

Sensitivity to various input values 237
We ignore complex, proximal fallout and concentrate on medial to distal areas, about 100 to 238 ~500 km downwind for example at Mount St (Fig. 5a). In simulations that use a single, 256 dominant aggregate size, these variations produce conspicuous changes in the location of a 257 secondary maximum (Fig. 8). Decreasing size also decreases the percent of erupted mass lands 258 in the mapped area: from 70% to 53% to 39% for d=0.165, 0.143, and 0.125mm respectively. 259 Our simulations limit µagg to values of 1.8-3.1φ (0.287-0.117mm), and σagg to 0.1-0.3φ, to 260 ensure that most deposits fall in the region of interest. 261 262

Results 263
We ran simulations at agg µ =1.8, 1.9, 2.0 . . . 3.1φ, and agg σ 0.1, 0.2, and 0.3φ. The latter used 264 1, 3, and 5 aggregate size classes respectively, in each simulation, with the percentage of fine 265 ash assigned to each bin given in Table 4 Figure 10b shows that modeled 305 and measured mass loads generally agree within a factor of three or so, except for those same 306 distal, low-mass-load measurements, to the lower left of the legend label (those where modeled 307 values are truly zero do not show up on this plot). Figure 10c shows that the modeled mass 308 load (black line with dots) contains a secondary thickening at about the same location mapped 309 (dashed line). However, the modeled mass load is consistently less than measured, especially 310 at the most distal sites. In Figure 10d, the log of modeled mass load versus square root of area 311 shows reasonable agreement with mapped values until mass loads are less than about 1 kg m -2 , 312 where they diverge. 313 Notably, modeled mass loads somewhat underestimate the measured values along the dispersal 314 axis in Fig. 10c. The underestimate reflects the fact that the input erupted volume of 0.2 km 3 315 DRE (Table 1) Fig. 10a; yet only about 79% of the modeled mass landed within this area. 317 Reducing the mean aggregate size to 2.7φ (0.153mm, Figs. S032-S034) improves the fit 318 somewhat along distal transect near the dispersal axis but not along the entire transect length. 319 And the finer size moves the secondary maximum too far east and reduces the percentage 320 deposited to 50-60%. 321 In Fig. 10a, the modeled deposit is also narrower than the mapped one. Adding turbulent 322 diffusion, with a diffusivity D of about 3×10 2 m 2 s -1 (Fig. 14) visually improves the fit, and was 323 likely important during this eruption due to high crosswind speeds that increased entrainment 324 Helens. In Fig. 11d, the areas covered by modeled isomass lines are comparable to the mapped 336 values, down to mass loads approaching 0.1 kg m -2 . 337

Ruapehu 338
For Ruapehu (Fig. 12), simulations using the NCEP Reanalysis 1 numerical winds produced an 339 odd double dispersal axis whose average did not correspond well with the mapped direction of 340 dispersal (Fig. 1c). To improve the fit we used the 1-D wind sounding provided for this eruption 341 at the IAVCEI Tephra Hazard Modeling Commission web page (http://dbstr.ct.ingv.it/iavcei/). 342 Use of a 1-D wind sounding seems justified in this case because this deposit covers a smaller 343 area than the others, making a 3-D wind field less important in calculating transport. The 344 resulting dispersal axis (Fig. 12a) agrees with the mapped one out to about 140 km distance, 345 beyond which it strays eastward, reaching the coast, 180 km downwind, about 10 km east of 346 the mapped axis. This slight difference is enough to cause misfits in point-to-point comparisons 347 at measured mass loads of ~10 -1 kg m -2 (Fig. 12b). 348 The modeled mass load along the dispersal axis ( Fig. 12c) agrees with measurements to about 349 60-90 km distance. At 100-200 km, modeled values level off and show a hint of secondary 350 thickening at ~180 km, in agreement with the mapped deposit ( Fig. 1c and 11c), although the 351 mapped secondary thickening is more prominent. 352 A large discrepancy is also apparent at distances of less than 60 km, where mass load along the 353 dispersal axis (Fig. 12c) and the area covered by thick isomass lines (Fig. 12d)  model. Underestimates of isomass area at <=10 -1 kg m -2 (Fig. 12d) also show that too little is 356 falling distally. Simulations (not shown) that raise the plume height or increase k to concentrate 357 more mass high in the plume do not improve the fit. The discrepancy may reflect the coarse 358 TPSD-50% of which is coarser than 1mm (compared with 2%, 12%, and 8% for the other 359 three deposits in Table S1). An additional simulation used the TPSD derived from technique 360 B of Bonadonna and Houghton (2005) (Table S1), which divides the deposit into arbitrary 361 sectors, and calculates a weighted sum of the size distributions in each sector following Carey 362 and Sigurdsson (1982). Technique B yields a finer average particle size than technique C, 363 which uses Voronoi tessellation to sectorize the deposit. But the finer particle size of the 364 technique B TPSD does not improve the fit (Fig. S173). Further exploration of this discrepancy 365 is beyond the scope of this paper; but other possible causes could include release of different 366 particle sizes at different elevations, or complex transport in the bending of the weak plume that 367 can't be accommodated in this model. 368 A second, smaller discrepancy is that the modeled deposit is narrower than the mapped one 369

Redoubt 382
This deposit is the second smallest in our group, the least well-constrained by sampling, and 383 the only one in our group not known to include a secondary thickness maximum. Mastin et al. 384 (2013b) modeled this deposit using numerical winds from the North American Regional 385 load distribution roughly agreed with mapped values for a plume height of 15km, k=8, and a 391 particle size adjustment that involved taking 95% of the fine ash (<0.063mm) and distributing 392 it evenly among the coarser bins. In this study we use the same plume height and k value, a 393 different wind field (RE1), and explore a different parameterization for particle aggregation. of fine (<0.063mm) ash (3-59%), atmospheric temperature, and water content between these 415 eruptions. The value of this narrow range depends strongly on other inputs, such as particle 416 density, shape factor, and Suzuki factor. But, holding those factors constant, the similarity in 417 this range between these four eruptions is noteworthy. 418 The overall agreement in modeled mean aggregate size ( agg µ ) suggests that accelerated fine-419 ash deposition may be treated as a discrete process, insensitive to eruptive style or magnitude. 420 It seems unlikely that these varied eruptions would produce aggregates of the same size, density, 421 and morphology. A combination of processes removed ash. Our approach captures these 422 processes implicitly, ignoring the microphysics. 423 What sort of processes could evolve in the cloud? Some possibilities are illustrated in Fig. A1. 424 The evolution starts with ejection of particles from that vent whose size ranges from microns 425 In the downwind cloud particle concentrations were lower, turbulence was less intense, a 433 smaller range of particle sizes existed, and, for all four eruptions, atmospheric temperatures 434 near the plume top were well below freezing (Table 5) Fig. A1). These types of distal aggregates tend to be smaller than a 478 millimeter, forming in the downwind cloud up to hundreds of kilometers from source (Sorem, 479 1982;Dartayat, 1932).         the freezing elevation was also checked using data from the North American Regional 775 Reanalysis (NARR) model (Mesinger et al., 2006), available at the same NOAA site, and 776 found to be 3.3 km, similar to that given below by the RE1 model. 777 778