Modeling ash fall distribution from a Yellowstone supereruption

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).

3.2 Volcanic Inputs

For erupted volume, we chose a fixed value of 330 km³ 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 km³ 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 km³ volume represents about one to two thirds of the total volume of a 500–1000 km³ eruption, which is above the threshold for a supereruption (1000 km³ bulk or ∼400 km³ DRE) [Sparks et al., 2005]. Expulsion of 330 km³ DRE of magma in a week would imply a mass eruption rate of 7×10⁸ kg s⁻¹, assuming a magma density of 2500 kg m⁻³. This is comparable to the ∼3–10×10⁸ kg s⁻¹ 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 km³ of magma over about 52 h [Hildreth, 1983]. The Tambora eruption of 1815, the largest in recent centuries, ejected about 50 km³ in 24 h [Self et al., 1984]. The climactic phase of the 1883 Krakatau eruption ejected 10 km³ 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 km³ 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 km³ 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.

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 m³ 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 km² 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⁻³, 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⁻³. Cornell et al. [1983] used this scheme to reproduce 38 ka Campanian Y-5 tephra at Campi Flegrei, Italy, a deposit of >60–100 km³ volume DRE [Cramp et al., 1989].

4.1 Calculated Growth of the Umbrella Cloud

Tephra distribution from a large-scale explosive eruption would be influenced both by growth of the umbrella cloud and by the ambient wind field. Umbrella clouds propagate upwind until the spreading velocity equals the ambient wind velocity, at which point the leading edge stagnates (Figure 2b) [Carey and Sparks, 1986]. Locations upwind of the stagnation point are unlikely to receive tephra fall deposits. Figure 5a plots the cloud radius with time using equation 5 for an eruption lasting 3 days (blue lines), 1 week (green), and 1 month (red). The expansion rate u_R at the cloud's leading edge, calculated from equation 6 is shown in Figure 5b and compared with the mean wind speed (black dashed line) at 16–24 km elevation above Yellowstone, obtained from RE1 data for 2000–2010 (Figure 4). The shaded gray region covers the range of wind speeds one standard deviation (σ) above and below the mean for this time period. Solid, dashed, and dotted line segments in Figure 5a indicate segments of the growth curve in which , , and , respectively.

These results show that, for the 3-day eruption, the expansion rate at the leading edge of the cloud could exceed the mean ambient wind speed for a period of days, producing a cloud that extends >2000 km upwind. For eruptions lasting a month or more, the expansion rate exceeds the mean wind speed for only several hours, producing a cloud whose maximum upwind extent is perhaps hundreds of kilometers.

4.2 Simulation Results

Figures 6-8 illustrate the simulated tephra thickness distribution in the conterminous United States and southern Canada for eruptions of 330 km³ DRE lasting a month (Figure 6), a week (Figure 7), and 3 days (Figure 8). Parts a, b, c, and d, in each figure illustrate seasonal effects, showing results in wind fields calculated for winter (a), spring (b), summer (c), and fall (d). For comparison, Figure 9 illustrates month-long (a), week-long (b), and 3-day (c) eruptions in winter using simulations that do not account for the presence of an umbrella cloud.

One feature is clear in Figures 6-8: proximal isopachs depicting >10 mm tephra thickness are mostly circular to elliptical and don't vary dramatically in size, shape, or coverage with season or with eruption duration. The shape of thinner isopachs (1 and 3mm tephra) vary more strongly in response to variations in the ambient wind field.

Comparison of simulations with an umbrella cloud (Figures 6-8) and without (Figure 9) reveal two important differences. First, neglecting the physics of umbrella propagation, patterns of tephra deposition are more sensitive to the spatiotemporal variations in the wind field. For example, the briefest eruption (3 days; Figure 9c) produces a relatively narrow tephra deposit that extends toward the southwest. Second, simulations with an umbrella cloud produce ash deposits that are much more widely dispersed. Without the umbrella cloud for example, the 10 mm isopach in the 3-day example (Figure 9c) covers 1.2 million km²; with the umbrella cloud (Figure 6a), it covers 5.4 million—4.5 times the area.

Table 3 lists the average, maximum, and minimum tephra-fall thicknesses for cities shown in Figures 6-8. Cities within 500 km such as Billings and Casper are covered by tens of centimeters to more than a meter of ash, upper Midwestern cities such as Minneapolis and Des Moines receive centimeters, and those in the East and Gulf Coasts receive millimeters or less. California cities receive millimeters to centimeters. And Pacific Northwest cities of Portland and Seattle receive up to a few centimeters.

The concentric pattern of isopachs suggests the possibility that tephra thicknesses might be expressed as a simple function of distance from the volcano. Figure 10 illustrates tephra-fall thickness as a function of distance from Yellowstone for cities listed in Table 3. A general decrease in thickness with distance is apparent, although cities west of Yellowstone (indicated by black arrows and labels) receive much less ash overall than cities farther east. There is no clear difference in thickness distributions between month-long, week-long, and 3-day eruptions.

4.2.1 Effects of Column Height and Grain-Size Distribution

Modern observations have supported the fundamental concept that eruption column height has a strong influence on tephra distribution [Carey and Sparks, 1986]. However, Figures 11a and 11b, which illustrate depositional patterns from an umbrella cloud at 15 and 35 km height, respectively, don't differ greatly from the 25 km column height (Figure 7b). The implication is that the dynamics of the powerfully spreading umbrella cloud dominates tephra dispersal in very large eruptions.

Grain-size distribution, however, does significantly affect tephra dispersal. This is illustrated in Figures 11c and 11d using grain-size distributions GSD2 and GSD3 for the 21–27 April simulation. The moderately aggregated GSD2 grain-size distribution, which is derived from the Mount St. Helens 18 May 1980 distal ash, shows similar dispersal to the heavily aggregated GSD1 (Figure 7b). In contrast, the weakly aggregated GSD3 shows much more widespread dispersal. Model results using GSD3 show that East Coast cities such as Washington, D.C., receive more than a centimeter, as opposed to the millimeters forecast under most GSD1 scenarios. The greater dispersal of the GSD3 size distribution likely results from the small size (0.2mm) and low density (200 kg m⁻³) of the single aggregate class, which contains 44% of the erupted mass. The true range of aggregate sizes and densities that might develop during such an eruption is difficult to anticipate, would depend strongly on atmospheric conditions, and is clearly an important factor to consider in long-distance ash dispersal.

The weakly aggregated GSD3 also places much more ash into distal clouds. The 21–27 April simulations using GSD1 and GSD2 deposited 97.5% and 99.7% of their erupted mass within the model domain. But the GSD3 simulation in Figure 11d sent 23.8% of the erupted mass beyond the model-domain boundary as an airborne cloud.

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 km³ 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).