Magma accumulation at depths of prior rhyolite storage beneath Yellowstone Caldera

The Yellowstone volcanic system has fueled some of the largest explosive caldera-forming eruptions in the geologic record, including three catastrophic eruptions in the past 2.1 million years (1, 2). Explosive silicic eruptions on this scale can have widespread environmental impacts, including continent-wide ash falls, global climate disruption, and extinction events (3, 4). At Yellowstone, the most recent Lava Creek eruption (0.64 million years ago) emplaced >1000 km³ of rhyolitic material and blanketed much of the western United States and Great Plains in ash (1, 5). The subsequent collapse of the magma reservoir shaped the current Yellowstone Caldera in northwestern Wyoming, which has since been filled with rhyolite flows as young as 70,000 years old (1). Although it is clear from geophysical observations that the modern Yellowstone magmatic system remains active (6, 7), questions persist about the volume and distribution of melt and how it compares with conditions that preceded prior eruptions.

An emerging view of continental magma reservoirs is that a crystal mush zone (a crystal-dominated, partially molten body) can persist in the crust over long time scales (100,000 years or longer) but that eruptible melt-rich zones are likely to be short-lived (<5000 years) (8–10). From this perspective, layers of eruptible silicic melt rapidly accumulate near the top of the crystal mush zone before eruption. Thus, the presence or absence of a melt-rich zone at or near the top of the silicic magma reservoir could be an important indicator of where Yellowstone currently sits in its eruptive life cycle.

Seismic tomography provides one of the best tools for inferring the presence of melt in crustal magmatic systems, and numerous previous studies have produced images of the subsurface below Yellowstone that reveal a contemporary magma reservoir in the mid- to upper crust (11–15). The magma reservoir in these studies is typically imaged as a slow shear wave speed (V_S) anomaly of up to 10%, suggesting a relatively melt-poor system (<10% partial melt fraction). However, spatially isolated observations of scattered teleseismic body waves recorded at seismometers in or near Yellowstone Caldera indicate that the degree of partial melt could be much higher (16).

Seismic tomography has yielded important clues into Yellowstone’s magmatic system, but imaging melt-rich zones remains challenging because small-scale magma bodies with severely reduced seismic wave speed are unlikely to produce substantial travel time delays of first arrivals due to wavefront healing (17). Additionally, strong low–wave speed anomalies may be further diminished in seismic images because of assumptions such as locally one-dimensional (1D) or ray-based seismic propagation and inversion regularization (18). Advances in tomographic inversions based on 3D numerical modeling of seismic waveforms, sometimes referred to as “full waveform inversion” (FWI), can overcome some limitations of conventional methods that rely on asymptotic ray-based approximations (19–21). The 3D sensitivity kernels used in FWI are able to account for complex wave propagation effects and thus more accurately map the location and amplitude of seismic wave speed anomalies.

We present new images of the V_S below Yellowstone based on FWI of ambient noise correlations. To take full advantage of the rich and diverse seismic datasets available in the Yellowstone region, our images combine data from numerous broadband deployments over the past 20+ years, including the EarthScope Transportable Array, several dense temporary deployments, and a recently updated seismic network within Yellowstone National Park (Fig. 1). Our tomographic inversion approach uses vertical component noise correlation functions (NCFs) from 4991 interstation pairs and minimizes frequency-dependent travel-time differences between NCFs and 3D synthetic waveforms in six overlapping period bands between 5 and 30 s (22). The final model m₁₀ was achieved after 10 adjoint iterations, and it reduces the total misfit by ~50% compared with the starting model m₀₀, which is based on conventional inversion techniques (22).

The tomographic images (Fig. 2) illuminate a strong V_S anomaly corresponding to the magma reservoir in the mid- to upper crust centered below the Yellowstone Caldera, with peak V_S reductions of >30%, which is substantially stronger than previously recognized. The slowest seismic wave speeds (V_S < 2.3 km/s) are present at depths between 3 and 8 km, with a minimum of 2.15 km/s at 5 km depth. Previous ray-based seismic tomography studies have imaged a mid- to upper crustal reservoir at ~5 to 15 km depth [for example, (13, 14)]; the peak velocity anomaly typically lies between 7 and 10 km depth, which is deeper than most petrological estimates of the storage depth of previously erupted rhyolitic magmas. For example, a recent petrological study of the Lava Creek Tuff suggested melt storage pressures of 80 to 150 MPa, corresponding to a depth range of ~3 to 6 km (23). Similarly, the storage pressure of CO₂-rich magmas from the Central Plateau Member Rhyolites (eruption ages 175,000 to 70,000 years) is estimated to be between 90 and 150 MPa (24). Thus, our tomographic images suggest contemporary magma storage in a depth range overlapping with the storage zone of magmas that have supplied both explosive and effusive silicic eruptions at Yellowstone.

In map view, the maximum V_S reduction is offset from the center of the caldera toward the east (Fig. 2A) and overlaps with a cluster of seismicity below Yellowstone Lake. The maximum depth extent of the low V_S region below the caldera is ~30 km, although the anomaly is more subdued below 10 km, which suggests that melt is most concentrated at shallower reservoir depths. In addition to the low-velocity anomaly below Yellowstone Caldera, two other low-V_S regions are notable. First is a region of low V_S in the lower crust (~35 to 40 km depth) that extends to the southwest of Yellowstone along the Snake River Plain and appears to connect to the anomaly below the caldera (Fig. 2B). There, V_S reaches 3.5 to 3.6 km/s (approximately –8 to –9% slower than the regional average), which is consistent with previous studies [for example, (12)]. This region lies above exceptionally slow mantle (12) and could represent a deep crustal reservoir of basaltic melt, although how it connects with the shallow silicic reservoir is unclear. Second, a region of low V_S to the southeast of Yellowstone Caldera is imaged in the mid-crust with a minimum V_S of ~3.4 km/s (Fig. 2C), which is enhanced compared with that in previous studies [for example, (15)].

Improved waveform fitting compared with our tomographic starting model demonstrates that the extreme V_S anomaly below the caldera is required by the data (Fig. 3 and figs. S8 to S10). Short-period waveform fits along paths that directly traverse the caldera are most notably improved. Shown in Fig. 3, A to C, is a comparison between observed and synthetic NCFs filtered between 6 and 9 s for a virtual source at station H18A. Although the starting model m₀₀ produces good waveform fits for most of the paths, a large travel-time delay (>5 s) was observed at station IPID, which is located directly opposite from Yellowstone Caldera so that most of the path samples the magma reservoir. The waveform observed along this path is well explained by the final FWI model m₁₀ (Fig. 3C). As further validation, we show waveform fits from a local moment magnitude (M_w) 4.2 earthquake that occurred to the northeast of the caldera on 25 March 2008 (Fig. 3, D to F). Data from this event was not used in the tomographic inversion; however, Rayleigh wave travel-time misfits were noticeably improved with model m₁₀. The improved fit is most appreciable at stations LKWY and YFT, where paths most directly sample the magma reservoir (Fig. 3F).

The presence of partial melt in the crust has the effect of reducing seismic wave speed, although the relationship between the melt fraction and the spatially averaged V_S structure as seen with seismic tomography is difficult to constrain and depends on temperature, composition, and the geometrical organization of melt in the crust. We estimated the melt fraction of Yellowstone’s upper crustal reservoir using a theoretical model of a solid-liquid composite with ellipsoidal melt inclusions defined by their aspect ratio (22, 25). We show the modeled relationship between V_S and melt fraction for various aspect ratios, calibrated for a rhyolitic composition at 5 km depth (Fig. 4A). Silicic partial melt in textural equilibrium exhibits a dihedral angle of 20° to 40°, corresponding to an aspect ratio of 0.1 to 0.15 (26, 27), which implies a crystal mush with a melt fraction of 16 to 20% near the top of the upper-crustal reservoir (Fig. 4B). Under these assumptions, we estimated the total volume of silicic melt in the upper-crustal reservoir to be >1600 km³ (22). On the basis of this melt fraction scaling, previous shear wave tomography models of Yellowstone’s magmatic system (12, 15) would map to ~10% melt or less. If melt is organized in networks of thin crystal-poor sills (15), the aspect ratio could be less than 0.1. However, sills may not contain 100% melt, and silicic partial melt is likely to exist at grain boundaries in crystal-rich portions of the magma reservoir. An aspect ratio of the magmatic system of <0.1 would imply a lower melt fraction but a stronger organization of melt into layered structures, which could decrease the stability of the system because an eruptible body could rapidly assemble from interconnected sills (28).

Mobilization and eruption of a crystal mush is possible when the melt fraction exceeds the critical threshold that marks the transition from a crystal-supported framework to a fluid suspension of crystals, which is accompanied by a dramatic viscosity decrease. Estimates of the critical melt fraction range from ~35% (16) to ~50% melt (29); thus, the melt fraction we estimated is substantially lower than what would be expected if a large fraction of the Yellowstone reservoir were in the eruptible stage of its life cycle. However, the presence of small subset volumes of concentrated silicic melt cannot be ruled out. For example, features smaller than the minimum seismic wavelength (in this case, ~15 km) may not be well resolved, suggesting that high–melt fraction bodies of several hundred cubic kilometers or more could be present in Yellowstone’s magma reservoir. Such subset volumes of the magma reservoir would be capable of supplying eruptions comparable in size with those of the 170,000- to 70,000-year-old Central Plateau Member Rhyolites, which included at least 17 flows with volumes of ~5 to 50 km ³ (1).

Rapid migration of silicic magma is difficult because of its high viscosity and the low permeability of crystal mushes; however, compaction-driven viscosity reductions from renewed melt injection can facilitate melt extraction on 1000-year time scales (30), and viscous “unlocking” of crystal mushes can be triggered by strain events that lubricate grain boundaries (31). Although our results indicate that Yellowstone’s magma reservoir contains substantial melt at depths that fueled prior eruptions, our study does not confirm the presence of an eruptible body or imply a future eruption. Strain events such as new magma intrusions or tectonic deformation that could begin to mobilize and concentrate magma would likely be accompanied by a host of dynamic processes evident to ongoing geophysical and geochemical monitoring (1, 6).