We used the volcanic ash transport and dispersion model Ash3d to estimate the distribution of ashfall that would result from a modern‐day Plinian supereruption at Yellowstone volcano. The simulations required modifying Ash3d to consider growth of a continent‐scale umbrella cloud and its interaction with ambient wind fields. We simulated eruptions lasting 3 days, 1 week, and 1 month, each producing 330 km 3 of volcanic ash, dense‐rock equivalent (DRE). Results demonstrate that radial expansion of the umbrella cloud is capable of driving ash upwind (westward) and crosswind (N‐S) in excess of 1500 km, producing more‐or‐less radially symmetric isopachs that are only secondarily modified by ambient wind. Deposit thicknesses are decimeters to meters in the northern Rocky Mountains, centimeters to decimeters in the northern Midwest, and millimeters to centimeters on the East, West, and Gulf Coasts. Umbrella cloud growth may explain the extremely widespread dispersal of the ∼640 ka and 2.1 Ma Yellowstone tephra deposits in the eastern Pacific, northeastern California, southern California, and South Texas.

1 Introduction The consequences of a future, caldera‐forming eruption from Yellowstone have been the subject of much speculation but little quantitative research in terms of regional ashfall impacts. Despite graphic and often fanciful media depictions of the devastation and the impact on human life that would result from a modern supereruption (producing >1000 km3 volcanic ash or >400 km3 DRE of magma), no historical examples exist from which to draw comparison [Self, 2006]. The largest eruptions of the past few centuries have produced a few to several tens of cubic kilometers of magma. Examples include Tambora volcano, Indonesia in 1815 [Oppenheimer, 2003], Krakatau in 1883 [Simkin and Fiske, 1983], the Katmai/Valley of Ten Thousand Smokes eruption, Alaska in 1912 [Hildreth, 1983], Quizapu volcano, Chile in 1932 [Hildreth and Drake, 1992], and most recently, Pinatubo, Philippines in 1991 [Newhall and Punongbayan, 1996]. These erupted volumes are much larger than the Mount St. Helens eruption in 1980 (0.2–0.4 km3) [Pallister et al., 1992], but at least an order of magnitude smaller than the largest Yellowstone events [Christiansen and Blank, 1972]. From geological evidence, we know that ash from the last large eruptions at Yellowstone (2.1 Ma, 1.3 Ma, and 640 ka) spread over many tens of thousands of square kilometers. These tephra deposits however have been remobilized in the millennia after emplacement, and estimates of primary ash thickness are challenging to obtain. Widespread Yellowstone‐derived deposits, known as the “Pearlette ash beds,” have been important stratigraphic markers throughout the central and western United States and Canada (Figure 1), even before their volcanic source location was recognized [Wilcox and Naeser, 1992]. Izett and Wilcox [1982] listed nearly 300 locations for Yellowstone ashes spread as widely as California, Texas, Iowa and Saskatchewan. Over 3 m of ash from the Lava Creek Tuff eruption (640 ka) are found in north‐central Texas. In the Gulf of Mexico, tens of meters of ash‐dominated deposits were emplaced immediately following the 2.1 Ma Huckleberry Ridge Tuff eruption [Dobson et al., 1991] and the later Lava Creek eruption. The great thickness of these and other deposits [Izett and Wilcox, 1982] partially reflects fluvial and (or) aeolian reworking, obscuring the primary distribution of ash thicknesses. In this study, we address the following questions: (1) given modern meteorological patterns, where and how much ash would be deposited by a Yellowstone supereruption? (2) how sensitive is the thickness distribution to changes in the season and duration of eruptive pulses? and (3) what is the long‐term (probabilistic) estimate of ashfall distribution? Figure 1 Open in figure viewer PowerPoint Location of sites where distal ashes from the Huckleberry Ridge Tuff (2.1 Ma) and Lava Creek Tuff (640 ka) eruptions have been identified [Izett and Wilcox, 1982], with some additional sites from A. Sarna‐Wojcicki, personal communication, March 2013). “DSDP 36” refers to Deep‐Sea Drill Hole 36, whose core contains Yellowstone ash deposits as described by Sarna‐Wojcicki et al. [1987].

2 Methodology We investigate these questions using the volcanic ash transport and dispersion model Ash3d [Schwaiger et al., 2012]. Ash3d is a finite‐volume Eulerian model that calculates tephra transport by dividing the atmosphere into a three‐dimensional grid of cells, placing tephra particles into cells above the volcano, and calculating their flux through cell walls as tephra is advected by wind and falls at a settling velocity determined by its size, density, and shape. Ash3d also calculates turbulent diffusion using a Crank‐Nicolson formulation with either a constant or variable diffusivity. Diffusivity has been adjusted in modeling studies to match the observed rate of downwind widening of a deposit or ash cloud, with different values required to match different observations [Bonadonna et al., 2005; Folch et al., 2009]. This procedure likely replicates the effects of more than one process, such as turbulent entrainment in weak plumes [Bonadonna et al., 2005], or umbrella cloud spreading in strong plumes [Costa et al., 2013]. Its effect is to increase the dispersal of ash; a process already accounted for to an exceptional degree by umbrella cloud spreading. Thus for simplicity we set diffusion to zero in these simulations. Suzuki, 1983 (1) The standard version of Ash3d places tephra into a column of source nodes above the volcano and distributes it vertically using a simple formula [] designed to concentrate mass near the plume top: Here is the mass flow rate into the plume, H T is the height of the plume top, z is the elevation at a particular point in the column, and k s is an adjustable constant that controls the degree to which mass is concentrated near the top of the column. Example curves using k s =4, 8, and 12 are shown in Figures 2a–2c (right side). Figure 2 Open in figure viewer PowerPoint Illustrations and model representations of (a) a weak plumes; (b) a strong plume that develops an anvil cloud (i.e., an umbrella cloud that is truncated on the upwind side); and (c) a major umbrella cloud. On the right‐hand side of each subfigure are curves showing the distribution of mass with elevation, calculated using equation 1, with values of k s as labeled. Gray nodes with magenta outlines are source nodes, to which ash is added at each time step in the model. Nodes with blue outlines are wind nodes, in which a radial wind field is added to simulate umbrella cloud expansion. Baines and Sparks, 2005 H T , Figure H u , Figure dV/dt is proportional to the volume rate at which the rising column feeds the cloud [Sparks et al., 1986 Sparks et al., 1997 (2) Approximating the plume as a column of source nodes is adequate for weak plumes (Figure 2 a), and those with moderate‐sized umbrella clouds (Figure 2 b), so long as they do not extend much farther upwind than a typical cell width. A supereruption however could drive an umbrella cloud thousands of kilometers upwind []. Clouds of this scale are fed by a column which rises buoyantly to a maximum height (, Figure 2 c), then collapses downward and outward as a density current that spreads radially at its neutral density elevation (, Figure 2 c). Their volume rate of growthis proportional to the volume rateat which the rising column feeds the cloud [.,.,, section 11.2]: has been inferred to follow the relationship [Morton et al., 1956 Sparks, 1986 (3) k e is the radial entrainment coefficient of the rising plume and N is the Brunt‐Väisälä frequency. C is a constant of proportionality, empirically determined through 3‐D modeling to be ∼0.5×104 m3 kg−3/4s−7/8 in tropical eruptions and 1×104 m3 kg−3/4s−7/8 in midlatitude and polar eruptions [Suzuki and Koyaguchi, 2009 k e =0.1, C=1×104 m3 kg−3/4s−7/8, and N=0.02 s−1 in our calculations. From dimensional considerations,has been inferred to follow the relationship [.,],whereis the radial entrainment coefficient of the rising plume andis the Brunt‐Väisälä frequency.is a constant of proportionality, empirically determined through 3‐D modeling to be ∼0.5×10kgin tropical eruptions and 1×10kgin midlatitude and polar eruptions []. We use=0.1,=1×10kg, and=0.02 sin our calculations. dP/dr driving outward motion results largely from the head drop between the cloud top (H T , Figure H U ) [Sparks et al., 1997 (4) R is cloud radius, is the mean cloud density, is ambient air density at the neutral buoyancy elevation, and g is gravitational acceleration. The pressure gradientdriving outward motion results largely from the head drop between the cloud top (, Figure 2 c) and the neutral buoyancy elevation () [.,, equation (11.5)],whereis cloud radius,is the mean cloud density,is ambient air density at the neutral buoyancy elevation, and g is gravitational acceleration. Sparks et al. [ 1997 t [Costa et al., 2013 Sparks et al., 1997 Suzuki and Koyaguchi, 2009 (5) Suzuki and Koyaguchi, 2009 Holasek et al., 1996 Sparks et al., 1997 Suzuki and Koyaguchi, 2009 Equations 2 and 4 account for conservation of mass and momentum in the cloud, respectively. They have been combined and integrated, with certain substitutions given in. [, chap. 11] to derive the following formula for cloud radius with time.,, equation 1 .,, equation (11.8);, equation 2 ],where λ is a factor reflecting cloud shape, whose value has been estimated at ∼0.2 from 3‐D modeling []. This equation has been shown to match the rate of growth observed in the 1991 Pinatubo eruption cloud [.,.,, Figure 11.3] and in 3‐D numerical simulations []. u R , can be estimated by taking the derivative of equation Costa et al., 2013 (6) The expansion speed of the cloud's outer margin,, can be estimated by taking the derivative of equation 5 with time [.,]: u r that approaches u R as r→R. Assuming a cylindrical cloud geometry and a nondivergent flow field, this leads to the relation [Costa et al., 2013 (7) Within the cloud, ash is assumed to spread at a radial velocitythat approachesas r→R. Assuming a cylindrical cloud geometry and a nondivergent flow field, this leads to the relation [.,]: We implemented these relationships in Ash3d using a modification of the method of Costa et al. [2013]. Specifically, an arrangement of source nodes consisting of a 3 × 3 matrix of cells (in plan view) extends over the upper 25% of the plume height, from the base to the top of the umbrella cloud (we do not include an overshooting top in our simulations) (Figures 2c, 3b, and 3c). From the vent elevation to the base of the umbrella cloud, the eruption plume source consists of a single column of nodes as in the standard version of Ash3d. At the beginning of each simulation, we calculate from equation 3. Then, at each time step, the existing tephra in each 3 × 3 layer of source nodes is averaged; new tephra is added to these cells according to the mass eruption rate, and distributed vertically using equation 1 with a top‐heavy k s of 12 (Figures 2c and 3d). Within each layer of cells, tephra mass is evenly distributed horizontally. Then, the cloud radius at that time step is calculated using equation 5 and, within the umbrella (purple region, Figure 3e), a radial wind field is calculated using equation 7 and added to the ambient wind field. Figure 3 Open in figure viewer PowerPoint (a) Illustration in Google Earth® of the model domain used for the Yellowstone model simulations. Our simulations used a grid spacing 0.5 degrees latitude and longitude, and 4 km cell height. (b) Top and (c) side view (looking north) of the source nodes above Yellowstone from which ash is dispersed. (d) Distribution of mass with height in the source nodes. (e) Radial wind field (white arrows) added to the ambient wind field in the umbrella region (purple). (f) Radial wind speed u r /u R in the cloud. Copyrighted images by Google (2011), Europa Technologies (2011), Tele Atlas, and Geocenter Consulting. Use of these images is consistent with usage allowed by Google (http://www.google.com/permissions/geoguidelines.html) and do not require explicit permission for publication. This method differs slightly from that of Costa et al. [2013] in that it uses equation 5 to calculate umbrella‐cloud radius at any given time rather than an equation that integrates dR/dt. This modification requires that the eruption be specified as a single pulse of constant plume height and eruption rate rather than a time series of heights and rates. The appendix Appendix A provides a validation of this approach by comparing Ash3d simulations to the well‐documented growth of the 1991 Pinatubo umbrella cloud. Ash3d does not resolve the dynamics of eruption plume ascent – instead, the height and vertical distribution of mass are held constant throughout the eruption. Despite a highly parameterized depiction of the plume, the model's key strength is its ability to examine the long‐distance umbrella expansion and dispersal of erupted material in a complex, time‐varying wind field.

3 Model inputs Inputs to the model include a 3‐D time‐varying meteorological wind field; volume of magma erupted; eruption start time; eruption duration; column height; and the size, density, and shape factor of erupted fragments. Our choice of inputs is described below. 3.1 Meteorological Inputs For meteorology, we use historical wind patterns represented in the NOAA NCEP/NCAR Reanalysis 1 model [Kalnay et al., 1996] (RE1). Model output gives the wind field every 6 h on a 3‐D grid spaced at 2.5° intervals in latitude and longitude, and at 17 pressure levels in the atmosphere, from 1000 mb (sea level) to 10 mb (about 34 km elevation). Wind vectors in the RE1 grid are linearly interpolated in space and time onto our grid at every time step. We have chosen representative time periods from the year 2001 for our simulations. This year was not strongly influenced by swings in either the El Niño Southern Oscillation (ENSO) [Di Lorenzo et al., 2010] or the Pacific Decadal Oscillation (PDO) [Mantua and Hare, 2002]. Wind speeds and directions above Yellowstone for this year (Figure 4, red dots) show reasonable agreement with longer‐term patterns between 2000 and 2010 (blue dots). Figure 4 Open in figure viewer PowerPoint Plots of wind direction and speed at Yellowstone at elevations of >24 km (a), 16–24 km (b), 11–16 km (c), 5–11 km (d), and <5 km (e). The radial position of each dot represents the wind speed, from zero (center) to 70 m s−1 (outermost circle); the azimuthal location represents the direction toward which the wind is blowing, with true north in the up direction. Blue dots represent all daily data from 1 January 2000 to 31 December 2010, whereas the red dots represent data for 2001. 3.2 Volcanic Inputs For erupted volume, we chose a fixed value of 330 km3 dense‐rock equivalent (DRE) of magma. The three major caldera‐forming eruptions from Yellowstone that produced the Huckleberry Ridge Tuff, Mesa Falls Tuff, and Lava Creek Tuff expelled about 2450, 280, and 1000 km3 DRE of magma, respectively [Christiansen, 2001]. But only a fraction of this volume rose in buoyant ash columns that could be carried by winds to form fall deposits; the remainder was emplaced either as ignimbrites that spread along the ground or as intracaldera fill that remained within the structural depression formed by collapse. The amount of these eruptions emplaced as tephra fall deposits is not well known, but is likely tens of percent based on better‐studied, recent caldera‐forming eruptions [e.g., Bacon, 1983; Wilson, 2001, 2008]. For our model simulations, the 330 km3 volume represents about one to two thirds of the total volume of a 500–1000 km3 eruption, which is above the threshold for a supereruption (1000 km3 bulk or ∼400 km3 DRE) [Sparks et al., 2005]. Expulsion of 330 km3 DRE of magma in a week would imply a mass eruption rate of 7×108 kg s−1, assuming a magma density of 2500 kg m−3. This is comparable to the ∼3–10×108 kg s−1 estimated for Pinatubo [Holasek et al., 1996; Koyaguchi and Ohno, 2001; Suzuki and Koyaguchi, 2009]. For eruption duration, we explore values from days to a month, reflecting durations that have been inferred and observed for moderate to large eruptions. The 1912 Novarupta eruption in Alaska, the largest volcanic eruption of the 20th century, ejected 15 km3 of magma over about 52 h [Hildreth, 1983]. The Tambora eruption of 1815, the largest in recent centuries, ejected about 50 km3 in 24 h [Self et al., 1984]. The climactic phase of the 1883 Krakatau eruption ejected 10 km3 of magma in about 36 h [Simkin and Fiske, 1983]. Field evidence from deposits of the 0.76 Ma Bishop Tuff from the Long Valley Caldera in California suggest that ∼700 km3 magma was emplaced in a matter of days [Wilson and Hildreth, 1997]. On the other hand, Wilson [2008] estimates that the ∼26 ka Oruanui eruption was a series of large‐scale outbreaks of increasing vigor, including multiple hiatuses of up to several weeks. Many smaller eruptions have persisted for weeks or even months. Eyjafjallajökull volcano in Iceland, whose eruption shut down air space over Europe in April and May of 2010, erupted about 0.2 km3 magma in 59 days [Gudmundsson et al., 2012]. For the height of ash injection, several factors are considered. Plume height is known to correlate with eruption rate, suggesting that a high‐flux Yellowstone eruption would produce a very high plume. But such correlations [e.g., Mastin et al., 2009, equation 1; Sparks et al., 1997, equation (5.1)] are based mainly on plumes that emanate from single, central vents, whereas a large, complex Yellowstone plume is more likely to rise from multiple vents, or as an elutriated ash cloud from pyroclastic flows. The 15 June, 1991 eruption of Pinatubo produce the highest historically observed plume of ∼40 km [Holasek et al., 1996], but its umbrella cloud spread at much lower heights of about 18–25 km [Fero et al., 2009; Holasek et al., 1996; Suzuki and Koyaguchi, 2009]. Based on these observations, most of our simulations use an umbrella cloud whose top is at 25 km. We also explore end‐member values of 15 km and 35 km. For a particle‐size distribution, we run most simulations using values listed in Table 1, referred to as GSD1. The size distribution primarily determines how the deposit is distributed with distance from the vent, with coarse particles settling near the vent and fine particles settling at greater distance. Ash3d calculates the settling rate of individual particles of ellipsoidal shape using relations of Wilson and Huang [1979], and a shape factor F (F=(b+c)/2a, where a, b, and c are the semimajor, intermediate, and semiminor axes of the ellipsoid, that is the average of tephra particles measured by Wilson and Huang [1979]. In real eruptions, particles smaller than about 0.063 to 0.125 mm aggregate into clusters and fall out faster than they would as individual particles [Van Eaton et al., 2012]. The physics governing the rate of aggregation is an area of active research and not yet explicitly resolved by Ash3d. Table 1. Grain‐Size Distribution GSD1, Used in Most Model Simulations Size (mm) Mass Fraction Density (kg/m3) 2 0.25 800 1 0.15 800 0.5 0.20 1000 0.25 0.15 1000 0.125 0.20 1800 0.0625 0.05 2000 To account for these effects, however, we consider a range of grain size distributions (GSD) that have been adjusted to match observed deposits. GSD1 has been adjusted to reproduce a heavily aggregated deposit—the 23 March 2009 (event 5) eruption of Redoubt Volcano, Alaska [Mastin et al., 2013b]. This eruption was small (∼1.6 M m3 DRE) [Wallace et al., 2013] and mixed with abundant external water from a crater‐filling glacier, resulting in extensive ash aggregation (A.R. Van Eaton et al., unpublished manuscript, 2014). Although a modern Yellowstone eruption may interact with external water (e.g., the 350 km2 Yellowstone Lake) magma‐water ratios would likely be higher than in the 2009 Redoubt eruption. Therefore, we consider two size distributions that represent less aggregation. The moderately aggregated GSD2 (Table 2) uses a modification of the total grain‐size distribution (TGSD) of the 18 May 1980 Mount St. Helens deposit estimated by Durant et al. [2009, supporting information file 2008jb005756‐txts03.txt], which integrated the deposit to a downwind distance of 671 km. To simplify calculations, we consolidated all tephra coarser than 1 mm into a single 1 mm size class. We then simulate aggregation by moving ash ≤0.125 mm diameter into four aggregate classes of size 0.177, 0.210, 0.250, and 0.297 mm with density 600 kg m−3, comparable to ash aggregates of low to intermediate density [Van Eaton et al., 2012]. This scheme can reproduce the 18 May 1980 Mount St. Helens tephra fall distribution in Ash3d simulations [Mastin et al., 2013a], and is reasonably consistent with field observations by Sorem [1982] that much of the very fine ash fell as clusters 0.25–0.5 mm diameter in the distal area. The weakly aggregated GSD3 (Table 2) starts with the Mount St. Helens TGSD of Durant et al. [2009], consolidates the coarser particles as above, and then moves 25% of the 0.031 mm particles, 75% of the 0.022 mm particles, and all particles smaller than 0.022 mm into a single aggregate class of 0.20 mm size and a density of 200 kg m−3. Cornell et al. [1983] used this scheme to reproduce 38 ka Campanian Y‐5 tephra at Campi Flegrei, Italy, a deposit of >60–100 km3 volume DRE [Cramp et al., 1989]. Table 2. Grain‐Size Distributions GSD2 and GSD3 Size (mm) GSD2 Moderate Aggregation Mass Fraction GSD3 Weak Aggregation Mass Fraction Density (kg/m3) Individual particles 1 0.043 0.043 1955 0.707 0.021 0.021 2209 0.5 0.046 0.046 2612 0.354 0.072 0.072 2634 0.25 0.053 0.053 2624 0.177 0.029 0.029 2690 0.125 0 0.033 2664 0.088 0 0.046 2697 0.063 0 0.060 2698 0.044 0 0.070 2640 0.031 0 0.060 2581 0.022 0 0.021 2570 Aggregate classes for GSD2 0.297 0.184 0 600 0.250 0.184 0 600 0.210 0.184 0 600 0.177 0.184 0 600 Aggregate class for GSD3 0.200 0 0.445 200

5 Discussion The simulations described above disperse tephra over distances that are qualitatively consistent with dispersal characteristics of known supereruption deposits. In most of our simulations, tephra thicknesses >1 cm cover a few to several million square kilometers. For comparison, the ∼170 km3 (DRE), ∼26 ka Oruanui fall deposit, the most recent and best preserved of any supereruption, covered more than a million km2 with ash thicknesses >1 cm [Wilson, 2001, Figure 9b]. The ∼800 km3 (DRE), 74 ka Toba fall deposit covered more than about 7 million km2 with >10 cm [Rose and Chesner, 1987]. T versus square root area, [Pyle, 1989 Fierstein and Nathenson, 1992 Pyle, 1989 (8) (9) (10) To examine how ash dispersal from the Yellowstone simulations compares with other deposits, Figure 12 a shows results on a plot of log thicknessversus square root area,] for the 21–27 April simulations using grain‐size distributions GSD1, GSD2, and GSD3. Using this method, the rate of thinning is expressed as either a single best‐fit line [e.g.,], according to:or as two line segments, using the following equations: Figure 12 Open in figure viewer PowerPoint (a) Log of deposit thickness (cm) versus square root of area covered (km) for simulations using grain‐size distributions GSD1 (blue triangles), GSD2 (red squares), and GSD3 (green diamonds) from 21 to 27 April. Dashed lines indicate best‐fit lines through these data, fit using equation 5. Dotted lines indicate a two‐line best fit in which the slopes and intercepts are automatically chosen to minimize squared residuals. The intersection of the two lines is indicated by a cross. Also listed are the best‐fit values of k and k 2 . (b) Plot of k or k 2 versus log bulk erupted volume for tephra deposits compiled in Fierstein and Nathenson [1992, Table 6], plus the value of k 2 estimated by Wilson [2001] for the ∼26 ka Oruanui tephra fall. Also shown are values of k and k 2 for the Yellowstone deposits illustrated in Figure 12a (green circles and crosses). Note that the x axis in Figure 12b represents bulk volume, consistent with values listed by Fierstein and Nathenson [1992] rather than DRE volume used elsewhere in this paper. In equation 8, T 0 is the y intercept and k is the negative slope of the line. In equations 9 and 10, T 1 and T 2 are the intercepts, and k 1 and k 2 are the negative slopes, of the proximal and distal lines respectively. The term represents the x coordinate at which the two lines intersect. Figure 12a plots a single best‐fit line through each data set as dashed lines, and a two‐line best fit as dotted lines, using an automatic procedure for picking the slopes and intercepts of these two lines that minimizes the squared residuals. Fierstein and Nathenson [1992] compiled values from the literature and found the slowest deposit thinning rates (lowest slopes) were associated with the largest, most widely dispersed eruptions. Figure 12b illustrates this trend by plotting the log of k and k 2 versus the log of bulk (not DRE) tephra volume, using data from Fierstein and Nathenson [1992, Table 6] along with one data point representing the distal slope of the Oruanui deposit [Wilson, 2001]. Regardless of whether one considers the total data set or separates out points for k or k 2 , a downward trend with increasing volume is visible at volumes >∼1 km3. At the right side of the plot, lying along the trend line, are values of k and k 2 from the three Yellowstone simulations in Figure 12a. Their location along this trend line suggests that the simulated tephra dispersal is reasonable for an eruption of this size. These results portray a dramatically different picture of the ashfall hazard from a caldera‐forming supereruption than one might expect based on traditional tephra transport models, which tend to neglect the physics of an expanding umbrella cloud. The distribution of tephra from such a large eruption is less sensitive to the wind field than that of smaller eruptions. It is also more widespread, and radially symmetric about the vent. The wide dispersal of tephra in our simulations results from rapid expansion of the umbrella rather from than a high plume or variability in the wind field at the umbrella spreading level during prolonged activity. 5.1 Immediate Ash Thickness Versus Long‐Term Impact North America's highest population density lies along its coastlines. Deposit thicknesses on the coasts from nearly all simulations is millimeters to a few centimeters. Thicknesses of this magnitude seem small but their effects are far from negligible. A few millimeters of ash can reduce traction on roads and runways [Guffanti et al., 2009], short out electrical transformers [Wilson et al., 2012] and cause respiratory problems [Horwell and Baxter, 2006]. Ash fall thicknesses of centimeters throughout the American Midwest would disrupt livestock and crop production, especially during critical times in the growing season. Thick deposits could threaten building integrity and obstruct sewer and water lines [Wilson et al., 2012]. Electronic communications and air transportation would likely be shut down throughout North America. There would also be major climate effects. Emission of sulfur aerosols during the 1991 Pinatubo eruption produced global cooling by an average of 1° C for a few years, while the 50 km3 Tambora eruption of 1815 cooled the planet enough to produce the famed “year without a summer” in 1816, during which snow fell in June in eastern North America and crop failures led to the worst famine of the 19th century [Oppenheimer, 2003]. Other indirect effects include wind reworking of tephra into migrating dunes that bury roads and structures; or increased sediment load to streams that exacerbates flooding and impedes river traffic. 5.2 Comparison With Reported Locations of Yellowstone Ash Deposits These simulations help us understand the extremely widespread distribution of the Huckleberry Ridge and Lava Creek eruption deposits from Yellowstone. Even under modern‐day prevailing winds (blowing dominantly to the east), the influence of a large umbrella cloud results in millimeters to centimeters of ash accumulation in southern California (Figures 6-8), which is consistent with field observations (Figure 1). Likewise, the presence of thick deposits of Huckleberry Ridge and Lava Creek Tuffs in Pleistocene Lake Tecopa in southeastern California [Izett and Wilcox, 1982] requires significant ashfall at least as far as southern Nevada to feed the ancestral Amargosa River. In the Los Angeles Basin, both eruptions delivered significant ash to the river systems in that area; whether windblown ash was postdepositionally remobilized to source these deposits is unknown, but we speculate that no deposits would be likely without at least several millimeters of ash fall in this region. Ash deposits in northern California and southern Oregon are easily explainable under the influence of a rapidly expanding umbrella cloud, even if prevailing winds were unfavorable to westerly transport (which is poorly known at 2.1 Ma and 640 ka). Without an umbrella cloud, westward transport may still occur within high‐level easterly winds (A. Sarna, written communication, 2014) (e.g., Figure 9b), but less frequently, at least under present‐day wind patterns (Figure 4). This again points to the important role of outward (umbrella) expansion in transport of ash from large‐scale eruptions. Yellowstone's Lava Creek B ash bed has also been found in cores from Tulelake on the California‐Oregon border [Rieck et al., 1992] (Figure 1), and Huckleberry Ridge deposits are located in Deep‐Sea Drill Site 36 in the eastern Pacific (Figure 1). The offshore deposits are not likely to have been transported fluvially from onshore [Sarna‐Wojcicki et al., 1987]. These occurrences are readily explained by an umbrella cloud extending at least 1500 km westward of Yellowstone, although transport to DSDP site 36 may require both an umbrella cloud and favorable winds (e.g., Figures 7a, 8a, and 8c).

6 Conclusions Geological activity at Yellowstone provides no signs that a supereruption will occur in the near future. Indeed, current seismicity, crustal deformation and thermal activity are consistent with the range and magnitude of signals observed historically over the past century [Lowenstern et al., 2006]. Over the past two million years, trends in the volume of eruptions and the magnitude of crustal melting may signal a decline of major volcanism from the Yellowstone region [Christiansen et al., 2007; Watts et al., 2012]. These factors, plus the 3‐in‐2.1‐million annual frequency of past events, suggest a confidence of at least 99.9% that 21st‐century society will not experience a Yellowstone supereruption. But over the span of geologic time, supereruptions have recurred somewhere on Earth every 100,000 years on average [Mason et al., 2004; Sparks et al., 2005]. As such, it is important to characterize the potential effects of such events. We hope this work stimulates further examination of ash transport during very large eruptions.

Acknowledgments We thank reviewers Andrei Sarna‐Wojcicki, Robert L. Christiansen, and Manuel Nathenson for comments that improved the paper. We would also like to express our gratitude to an anonymous reviewer who recommended incorporation of an expanding umbrella cloud into Ash3d, and pointed us to the approach of Costa et al. [2013]. Second author AVE acknowledges NSF Postdoctoral Fellowship EAR1250029 for funding at the time of her contribution to this work. The data and model results used in this paper have been posted as supporting information. The supporting information includes input files and ASCII or kmz output files that can be opened in ArcMap® or Google Earth® software (use of trade names does not imply endorsement of these products by the U.S. Government).

Appendix A: Simulating the Growth of the Pinatubo Umbrella Cloud We validate the umbrella cloud formulation by simulating the growth of the 15 June 1991 umbrella cloud at Pinatubo. This is the only large umbrella cloud observed with modern techniques and its growth has been examined in numerous studies. The model setup and atmospheric parameters are indicated in Table A1. For meteorological conditions we use the RE1 wind field. Typhoon Yunya prevented the launching of radiosondes during the eruption [Fero et al., 2009] and limited the accuracy of modeled wind fields. Thus we, like others [Costa et al., 2013; Fero et al., 2009], had to rotate wind directions 30 degrees counterclockwise to match the observed direction of cloud movement. Table A1. Model Parameters Used to Simulate Growth of the 15 June 1991 Umbrella Cloud at Pinatubo Parameter Value(s) Eruption start time (UTC) 15 Jun 1991, 0440UTC (1340 local time) Height of top of umbrella cloud 25 km Duration 9 h Erupted volume 6 km3 DRE (1.5×1013 kg) Model domain 89.558–132.454° longitude 1.227–28.846° latitude Model resolution 0.2° horizontally, 2 km vertically Grain sizes One grain size 0.01 mm diameter 1000 kg m−3 density Diffusion constant 0 m2 s−1 k e 0.1 N 0.02 s−1 C 0.5×104 m3 kg−3/4s−7/8 λ 0.2 Our umbrella‐cloud formulation requires a single eruptive pulse of constant rate and umbrella‐cloud height. We use a starting time at 1340 Pinatubo Daylight Time (0440 UTC) on 15 June 1991; a duration of 9 h; and an umbrella‐cloud height of 25 km, which corresponds to the cloud top observed in satellite images [Holasek et al., 1996; Koyaguchi and Tokuno, 1993], but is lower than the 35–40 km of the overshooting top [Holasek et al., 1996]. The erupted volume has been estimated 4.8–6.0 km3 DRE from deposits [Wiesner et al., 2004, 2005]. We use an erupted volume of 6 km3 DRE, and assume a magma density of 2500 kg m−3. Figure A1 (a) shows hourly outlines of the observed cloud perimeter starting at 1440 PDT, digitized from Fig. 5b of Holasek et al. [1996]. Figure A1 (b) shows outlines of the cloud at the same times, based on the Ash3d simulation. The result illustrates reasonably good agreement. Figure A2 shows cloud diameter versus time obtained from the observations of Holasek et al. [1996, Fig. 5b], from the Ash3d simulation, and from same theoretical method used in Fig. 5a, assuming an erupted volume of 6 (dashed line) or 10 (solid line) km3 DRE. For the Holasek et al. and Ash3d results, cloud diameter d was determined by measuring the cloud area A and using the formula d = 2sqrt(A/π). The Ash3d cloud diameters agree surprisingly well with those of Holasek et al., especially given that the simulation assumed a constant eruption rate with time.

Supporting Information Filename Description ggge20543-sup-0001-suppinfo01.txt2.6 KB Readme ggge20543-sup-0002-suppinfo01.zip561.4 KB Simulation results ggge20543-sup-0003-suppinfo02.zip783.1 KB Simulation results ggge20543-sup-0004-suppinfo03.zip315.3 KB Simulation results ggge20543-sup-0005-suppinfo04.zip511.1 KB Simulation results ggge20543-sup-0006-suppinfo05.zip451 KB Simulation results ggge20543-sup-0007-suppinfo06.zip579.5 KB Simulation results ggge20543-sup-0008-suppinfo07.zip773.3 KB Simulation results ggge20543-sup-0009-suppinfo08.zip522.4 KB Simulation results ggge20543-sup-0010-suppinfo09.zip554.2 KB Simulation results ggge20543-sup-0011-suppinfo10.zip558.9 KB Simulation results ggge20543-sup-0012-suppinfo11.zip549.4 KB Simulation results ggge20543-sup-0013-suppinfo12.zip587.6 KB Simulation results ggge20543-sup-0014-suppinfo13.zip562.1 KB Simulation results ggge20543-sup-0015-suppinfo14.zip762.1 KB Simulation results ggge20543-sup-0016-suppinfo15.zip539.8 KB Simulation results ggge20543-sup-0017-suppinfo16.zip567 KB Simulation results ggge20543-sup-0018-suppinfo17.zip768.2 KB Simulation results ggge20543-sup-0019-suppinfo18.zip518 KB Simulation results ggge20543-sup-0020-suppinfo19.zip561.2 KB Simulation results ggge20543-sup-0021-suppinfo20.zip3.6 MB Simulation results ggge20543-sup-0022-suppinfo21.zip168.5 KB Wind data ggge20543-sup-0023-suppinfoA1.pdf578.2 KB Supporting Information Figure A1 ggge20543-sup-0024-suppinfoA2.pdf397.3 KB Supporting Information Figure A2 ggge20543-sup-0025-suppinfo25.pdf43.6 KB Supporting Information Files ‐ Corresponding Text Figures Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.