Shear Velocity Evidence of Upper Crustal Magma Storage Beneath Valles Caldera

Valles Caldera was formed by large rhyolitic eruptions at ∼1.6 and 1.23 Ma and it hosts post‐caldera rhyolitic deposits as young as ∼69 ka, but the contemporary state of the magmatic system is unclear. Local seismicity beneath Valles Caldera is rare and shear‐velocity (Vs) structure has not been previously imaged. Here, we present the first local Vs tomography beneath Valles Caldera using ambient noise Rayleigh dispersion from a ∼71 km transect of nodal seismographs with mean spacing of ∼750 m. An ∼6 km wide low‐Vs anomaly (Vs < 2.1 km/s) is located at ∼3–10 km depth within the 1.23 Ma caldera's ring fracture. Assuming magma in textural equilibrium, the new tomography suggests that melt fractions up to ∼17%–22% may be present within the upper crustal depth range where previously erupted rhyolites were stored.

Seismic imaging provides insight into the contemporary abundance of magma in the upper crust beneath Valles Caldera. Early P-wave studies indicated low P-velocity (Vp) and elevated attenuation beneath the caldera (Ankeny et al., 1986;Roberts et al., 1991Roberts et al., , 1995. Teleseismic P-wave data from the 1993-1994 Jemez Tomography Experiment (JTEX) enhanced resolution within the caldera and revealed a low-Vp anomaly between ∼5 and 20 km depth in the shape of a vertically elongated ellipsoid with a Vp reduction of −23% (Steck et al., 1998). Following JTEX, there was a long hiatus in data collection for local imaging. Denser arrays and application of newer seismic methods such as ambient noise surface wave tomography could add valuable S-velocity (Vs) constraints with complementary sensitivity to melt, improve depth-resolution, and better facilitate comparison of Valles Caldera to other systems that may be in a similar life-cycle stage (e.g., Schmandt et al., 2019).
Here, we present the first local Vs tomography beneath Valles Caldera by applying ambient noise Rayleigh wave tomography to data from a new dense seismic transect (Figure 1). Short-period Rayleigh wave dispersion constrains absolute Vs in the middle to upper crust and the mean seismograph spacing of ∼0.75 km provides the ability to recover local variations in structure along the ∼71-km transect. The new Vs tomography results are used to estimate the potential depth interval and concentration of magma beneath Valles Caldera.  (Morgan et al., 1996) and squares are borehole locations. Borehole B-12 is the deepest at ∼3.2 km (Nielson & Hulen, 1984). A yellow hexagram shows the virtual source used in panel (c). Line, A-A′, delineates the tomographic cross section. The left inset shows regional physiographic provinces: Colorado Plateau, Basin and Range, Rocky Mountains, and Rio Grande Rift (RGR). Semi-transparent tan fill shows the RGR. Solid black fill shows Cenozoic volcanic fields of the Jemez lineament. Solid red fill shows the Jemez volcanic field which includes Valles Caldera and the black square outlines the extent of the main figure. (b) Noise cross-correlations are shown filtered from 3 to 9 s and stacked at a distance interval of 1 km. (c) A virtual source gather is shown using a northern station and 3-5 s bandpass filter. Black dotted lines denote the 1.23 Ma caldera's ring fracture.

Data
Continuous seismic data were recorded with 97 three-component Magseis-Fairfield nodal seismometers. The stations were deployed along a ∼71 km NNE striking linear transect across Valles Caldera ( Figure 1) with a spacing of ∼750 m and operated between 29 September 2019, and 9 November 2019. The nodal seismographs were coupled to the ground with stakes but not buried to minimize environmental impact. Ten seismographs tipped during deployment, presumably due to wildlife interactions based on frequent observations of elk and cattle. Three of the tipped stations had >10 days of data that were recovered by identifying the day that inter-station noise correlations abruptly changed. Tipped nodes are noted in Supporting Information S1 (Table S1).

Ambient Noise Correlations
Ambient noise cross-correlation functions were calculated from ambient seismic noise following Bensen et al. (2007). The data were down-sampled to 10 Hz, bandpass filtered from 0.02 to 2 Hz, normalized in the time domain using a running absolute mean, and then whitened before cross-correlation. Correlations were computed for 4-hr half-overlapping time windows throughout the continuous data (e.g., Seats et al., 2012) and then all correlations for each station pair were linearly stacked to preserve phase (Yang et al., 2023). We focus on the vertical-vertical (ZZ) component cross-correlations which show clear fundamental mode Rayleigh waves on the positive and negative lag portions of the symmetric cross-correlations ( Figure 1b).

Phase Velocity Dispersion
Phase velocity dispersion curves were calculated using the automated frequency-time analysis method (Bensen et al., 2007;Levshin et al., 1972). We focus on 3-9 s periods that are sensitive to the upper and middle crust ( Figure  S1 in Supporting Information S1). Below 3 s period clear fundamental mode Rayleigh waves were not observed. Beyond 9 s period the ∼71-km long transect provides few measurements with inter-station spacing greater than our minimum of 1.5 wavelengths. To boost the signal-to-noise ratio as well as to reduce the potential effects of inhomogeneous noise source distribution, we averaged the positive and negative portions of the cross-correlation functions. We retained measurements with phase velocities between 1 and 4.6 km/s, signal-to-noise ratio > 5, and wavelengths > 1.5, resulting in 2,749 total dispersion curves available for phase velocity tomography ( Figure S2 in Supporting Information S1).

Phase Velocity Tomography
A damped least-squares method was used to invert inter-station phase velocity dispersion measurement for phase velocity along the transect for periods between 3 and 9 s (e.g., Wilgus et al., 2020). At 4 s period there are ∼1,300 inter-station phase velocity measurements used for tomography. This number decreases with increasing period as the 1.5 wavelength requirement becomes a larger fraction of array length ( Figure S3 in Supporting Information S1). The inversion was conducted with straight rays on a 0.5 km grid along linear transect A-A′ shown in Figure 1. To reduce the influence of potential outliers among the dispersion measurements, a two-stage inversion approach was used in which measurements with travel time residuals beyond two standard deviations after the first inversion were removed and then the inversion was repeated (e.g., Wang et al., 2017).

Shear Velocity Inversion and Modeling
Phase velocity dispersion curves from 3 to 9 s period were extracted for 144 locations along the transect to invert for Vs as a function of depth from the surface. We closely followed the workflow of Wilgus et al. (2020) using a Bayesian Markov chain Monte Carlo (BMMC) approach to obtain an ensemble of Vs models capable of fitting the dispersion measurements (Shen et al., 2013). The subsurface Vs structure is represented by a total of 9 parameters, consisting of 7 b-splines for Vs in the middle to upper crust, a depth transition to an underlying half-space at 20-23 km, and Vs in the half-space. To accommodate potentially strong heterogeneity in the upper crust and diminishing resolution with depth, the prior distribution is wider in the upper crust and narrows with depth ( Figure S1 and Table S2 in Supporting Information S1). Phase velocity sensitivity kernels show that 90% of the sensitivity for the longest period Rayleigh wave, 9 s, is located above 21 km depth. Consequently, the parameterization transitions from 7 b-splines to a half-space within a depth range of 20-23 km ( Figure S1 in Supporting Information S1). The BMMC inversion explores the model space, iteratively predicting dispersion (p i ), and evaluating the fit to the observed dispersion (d i ) with a reduced Chi-square misfit equation where is uncertainty and n is the number of discrete periods. A total of 2 million iterations were used at each point along the transect. The best 2,000 models at each location are used as the posterior distribution and the mean of the ensemble is used for constructing the final Vs profile. All forward calculations of phase velocity were conducted with software from Computer Programs in Seismology (Herrmann, 2013) using empirical crustal rock scaling relationships between Vs, Vp, and density (Brocher, 2005).
Multiple estimates of phase velocity uncertainty were tested in the Vs inversions. We chose to use a fixed value of 25 m/s, which results in a mean χ 2 of 1.7 ( Figure S4 in Supporting Information S1). An alternate approach using bootstrap resampling and repeated phase velocity tomography (e.g., Jiang et al., 2018) resulted in smaller uncertainties and a greater mean χ 2 of 14.8. Results based on different uncertainty choices show that the geometry of major features of the model is stable but small variations in the amplitude of velocity anomalies are present ( Figure S5 in Supporting Information S1). For instance, using the smaller uncertainties from the bootstrap approach results in a slightly lower minimum Vs of ∼1.95 km/s rather than ∼2.0 km/s.

Teleseismic P-Wave Relative Delay Times
The dense spacing of the nodal seismograph array provides an opportunity to observe teleseismic P-wave residual times as a simple metric of consistency with prior P-wave studies and consistency between any major Vp and Vs anomalies. However, the brief deployment did not provide many high-quality events. One of the clearest teleseismic P-waves observed when most nodes were operating is shown in Figure 2. The event occurred in Japan and the P-wave approaches Valles Caldera from the northwest. The seismograms were bandpass filtered from 0.25 to 0.75 Hz and aligned based on travel time predictions for the AK135 velocity model (Kennett et al., 1995). Seismograms with signal-to-noise amplitude ratios <3 were removed. Relative alignment of the seismograms with multi-channel cross correlation was used to identify Vp anomalies sampled by steeply incident P-waves  (Steck et al., 1998). Missing nodes were either not recording during the event or had signal-to-noise amplitude ratios <3. (d) Period averaged (3-5 s) phase velocity is plotted along A-A′. (VanDecar & Crosson, 1990). Since there is ∼1.8 km topographic relief along the transect, we applied elevation corrections assuming an upper crustal Vp of 5.5 km/s based on estimates from controlled source and earthquake travel time tomography (Ankeny et al., 1986).

Phase Velocity Pseudo Cross-Section
The depth of peak sensitivity for Rayleigh wave phase velocity increases with period such that plotting the phase velocities beneath each point along the transect provides a pseudo cross-section perspective on local crustal structure (Figure 3b). The most prominent feature of the pseudo cross-section is a low-velocity anomaly for periods between ∼3 and 5 s beneath the Redondo Peak resurgent dome, with the lowest velocity found for ∼4 s period (Figures 2d and 3c). Beneath the central caldera across a width of ∼6-8 km within ring fracture, 3-5 s phase velocities are ∼12%-15% lower than the array mean, ∼2.3-2.4 km/s (Figures 2d and 3b). The magnitude of the low-velocity anomaly makes it visible in virtual source gathers of noise correlations as a deflection in the Rayleigh wave arrival as it crosses the central caldera within the ring fracture (Figure 1c). Phase velocity tomography resolution tests demonstrate that a prominent low-velocity anomaly within the caldera's ring fracture is resolvable with the available data coverage, but the magnitude of the velocity anomaly would be slightly underestimated ( Figures S6 and S7 in Supporting Information S1). A test with an input anomaly of −20% in phase velocity across a 6 km width resulted in a recovered minimum velocity of −17.5% ( Figure S6 in Supporting Information S1). Other phase velocity features include low velocities (∼2.5 km/s) at 3-4 s period on the northern flank of the caldera (∼45-55 km transect distance) and high velocities (∼3.1 km/s) at 3-4.5 s on the southern flank of the caldera (∼15-23 km transect distance; Figure 3b).

Shear Velocity Cross-Section
Inversion for Vs provides constraints on absolute Vs as a function of depth. Extremely low Vs (<2.1 km/s) is found from ∼3 to 10 km depth across a width of ∼6 km within the 1.23 Ma caldera's ring fracture (Figure 3c). The low-velocity anomaly within the Vs < 2.1 km/s contour corresponds to an ∼32% Vs reduction compared to the mean across the array. The anomaly's location under the Redondo Peak resurgent dome is slightly offset from the area of highest geothermal gradients, however most boreholes were drilled on the west side of Redondo Peak, whereas the nodal array crossed closer to the center of the dome (Figure 1). The highest χ 2 misfit values (>4) are situated on either side of the low-Vs anomaly ( Figure S4 in Supporting Information S1), suggesting that phase velocities near the edges of the anomaly are difficult to fit with locally 1D velocity structure. Aside from the low-velocity anomaly beneath Redondo Peak, there is a low-velocity anomaly with Vs of ∼2.2-2.5 km/s on the northern flank of the caldera (transect distance of ∼45-55 km), but unlike the central caldera anomaly its depth extent is restricted to the uppermost ∼2 km (Figure 3c, Figure S8 in Supporting Information S1).

Teleseismic P-Wave Lag Times
P-wave travel time lags of up to ±0.65 s were measured for a clearly recorded teleseismic earthquake (Figure 2c). A deflection in the P-wave arrival can be seen in the waveforms recorded within the caldera's ring fracture (Figure 2b). Delayed arrivals, indicative of low Vp at depth, are located on the Redondo Peak resurgent dome within the Valles Caldera ring fracture (Figure 2b). The location of the most delayed arrivals is consistent with previous P-wave studies that used seismographs with larger inter-station spacing but distributed over the area of the caldera rather than in one transect (Roberts et al., 1991;Steck et al., 1998). The along-transect distance of the most delayed arrivals, ∼0.5-0.65 s, coincides with the area of highest geothermal gradients (Morgan et al., 1996). More moderate lag times of ∼0.3 s extend south from Redondo Peak to ∼10 km beyond the southern edge of the caldera (Figure 2c). The broad width of the ∼0.3 s delays is consistent with a deeper origin located west of the nodal transect based on the ∼315° back-azimuth of this event (Figure 2c).

Discussion
The new Rayleigh wave tomography advances insights into local Vs structure within and surrounding Valles Caldera. The primary result is strongly reduced Vs beneath the resurgent dome, Redondo Peak. A secondary low-Vs anomaly, ∼2.2-2.5 km/s, located in the uppermost 2 km on the northern flank of the caldera is more likely related to the history of volcanic deposition in the area (Figure 3c). Beginning in the mid-Miocene there was intermediate-to-mafic volcanic activity on the north side of the caldera and the resulting extrusive rock deposits may cause lower Vs in the uppermost crust that contrasts with the southern flank of the caldera (Gardner et al., 1986;Goff et al., 2011). The central caldera low-Vs anomaly and its potential implications for the contemporary state of the magmatic system are the focus for the remainder of the discussion.
Within the caldera the lowest Vs volume is concentrated between ∼3 and 10 km depth over a width of ∼6 km within the 1.23 Ma ring fracture, where Vs is ∼2.0-2.1 km/s (Figure 3c). The relative Vs anomaly within that volume is −32%. Prior teleseismic P-wave tomography estimated a relative velocity anomaly of −23% within a more vertically elongated ellipsoidal anomaly (Steck et al., 1998). Reflectors previously identified by P-wave coda migration at ∼4 km and ∼9-14 km below the surface may represent the vertical boundaries of the magma reservoir, but the upper reflector may alternatively be related to the contact between tuff deposits and the underlying basement rock (Aprea et al., 2002). The availability of absolute Vs in the upper to middle crust from this study provides valuable new constraints for estimating the origin of the low-velocity anomaly. Given the Vs anomaly's location under the resurgent dome and that its minimum velocities are located between ∼3 and 10 km depth, it cannot be explained by unconsolidated caldera fill. The depth range of the anomaly overlaps petrologically estimated storage depths of erupted rhyolites, ∼2.5-9 km (Boro et al., 2020;Spell & Kyle, 1989;Wilcock et al., 2013). So, we proceed to interpret silicate partial melt and magmatic volatiles as probable contributors to the low-velocity anomaly.
We first consider a base scenario in which melt fraction is estimated assuming a composition like that of the ∼69 ka Banco Bonito rhyolite flow (Table S3 in Supporting Information S1) and that partial melt in the subsurface today is in textural equilibrium. Then we proceed to discuss uncertainties that could lead to over-or under-estimation of the melt fraction. To predict Vs as a function of the melt fraction, we used the theoretical model of Berryman (1980) for an elastic medium with ellipsoidal fluid inclusions (e.g., Paulatto et al., 2019). Elastic properties of the solid were calculated with Perple_X assuming bulk composition of the Banco Bonito rhyolite, pressure of 170 MPa (∼5 km depth), and a temperature of 700°C (Connolly, 2009). The velocity reduction due to partial melt depends on the assumed aspect ratio of intergranular melt pockets and aspect ratios of ∼0.1-0.15 are expected for textural equilibrium (Takei, 2002). In this scenario, Vs of 2-2.1 km/s would correspond to melt fractions of 17%-22% ( Figure S9 in Supporting Information S1). We consider textural equilibrium a reasonable assumption because the system has not erupted since ∼69 ka and it hosts little seismicity, so any deformation and magmatic recharge are expected to be slowly evolving processes.
Uncertainties in seismic imaging and the multi-phase structure of the magma reservoir could bias the estimated melt fractions. This study benefits from a dense local array, but simplifying assumptions include a 2D phase velocity inversion and straight ray paths. Resolution tests using these assumptions indicate that ∼85% of the input velocity anomaly amplitude could be recovered for a low-velocity anomaly like that imaged beneath Redondo Peak ( Figures S6 and S7 in Supporting Information S1). A surface wave tomography resolution study using 3D synthetic waveforms (e.g.,  would provide a more realistic assessment but it is not considered feasible within the scope of this study. Insights from 3-D full wave synthetic tests conducted by  suggest that conventional surface wave travel time tomography is likely to underestimate the true magnitude of Vs reduction in crustal magma reservoirs, but the problem is more subdued for magma reservoirs that are large with respect to the inter-station spacing. In this study, the velocity anomaly of interest is ∼6 km wide in the upper crust and the mean inter-station spacing is <1 km so we do not expect severe underestimation. If the seismic properties are influenced by magmatic volatiles that could bias our interpretation of the melt-fraction toward over-estimation. Valles magmas may have several percent dissolved volatiles based on past eruptions (Boro et al., 2020;Waelkens et al., 2022). As magma cools in the upper crust buoyant volatiles may accumulate in a thin low-velocity layer atop the magma reservoir (e.g., Seccia et al., 2011). Such a scenario is plausible at Valles Caldera given that boreholes encountered high geothermal gradients up to 3.2 km deep ( Figure 1) but did not reach magmatic fluids or recently cooled intrusions (Goff & Gardner, 1994;Nielson & Hulen, 1984), so some volatiles could be trapped beneath a low-permeability boundary (De Siena et al., 2017;Vanorio et al., 2005). Resolution of a potential layer of magmatic volatiles atop a silicate melt reservoir may be possible with higher frequency Vp/Vs imaging, such as that conducted at Campi Flegrei, because lower Vp/Vs is expected for a trapped volume of exsolved volatiles compared to silicate melt (Calò & Tramelli, 2018;Chiarabba & Moretti, 2006). However, the paucity of seismicity beneath Valles Caldera hinders the ability to conduct similar imaging. Comparison of teleseismic P-wave tomography and Rayleigh tomography is complicated by differing resolution, but the existing results do not suggest a low Vp/Vs anomaly since the Vs reduction (−32%) is greater than the Vp reduction (−23%) reported by Steck et al. (1998). A dominant role for exsolved volatiles in creating the low-velocity anomaly further seems unlikely because the high geothermal gradients in Valles Caldera are consistent with continued presence of melt (Morgan et al., 1996;Nielson & Hulen, 1984).
Comparison of the Vs structure beneath Valles Caldera to that of other active silicic volcanic fields suggests an active magmatic system. The Laguna del Maule volcanic complex, which has hosted many Holocene rhyolite eruptions, exhibits similar Vs characteristics with a minimum Vs of ∼2.0 km/s at ∼4 km depth based on ambient noise Rayleigh wave tomography (Wespestad et al., 2019). Beneath Yellowstone Caldera the minimum Vs of ∼2.15 km/s at ∼3-8 km depth is comparable to Valles Caldera , but the width of Yellowstone's upper crustal anomaly is up to ∼60 km in comparison to the ∼6 km width imaged here. Long Valley Caldera's seismically imaged reservoir appears deeper with a top at ∼8 km and an underlying Vs anomaly extending to ∼20 km depth with a minimum Vs of ∼2.5 km/s (Flinders et al., 2018;Nakata & Shelly, 2018). Valles Caldera's Vs structure is more like the examples of Laguna del Maule and Yellowstone where the depth interval of the low-Vs anomaly largely overlaps typical pre-eruptive rhyolite storage depths of ∼4-10 km (Bachmann & Bergantz, 2008). Given the 2D geometry of our study we refrain from detailed estimates of the 3D volume of magma beneath Valles Caldera, but a simple approximation of a cylindrical volume with radius of 3 km and depth interval of 7 km would enclose most of the low-Vs anomaly. An ∼20% melt fraction within the cylinder would correspond to ∼40 km 3 of magma. Only a fraction of this volume (∼10%-25%, e.g., Annen et al., 2008) would need to be mobilized to fuel eruptions analogous to those that produced the ∼4 km 3 El Cajete pyroclastic deposits from ∼74 ka or the ∼6.8 km 3 Banco Bonito rhyolite flow from ∼69 ka (Nasholds & Zimmerer, 2022;Wolff et al., 2011;Zimmerer et al., 2016).
Evidence for upper crustal magma storage beneath Valles Caldera highlights challenges for hazard monitoring and research. Regional scale geophysical studies show that seismicity is locally absent beneath Valles Caldera (Nakai et al., 2017) and the adjacent Rio Grande Rift deforms slowly with an extension rate of ∼0.1 mm/year (Berglund et al., 2012). Unlike similar settings in the United States such as Yellowstone Caldera or Long Valley Caldera, Valles Caldera does not have continuous open-access seismic or ground-based geodetic data. There is a local seismic network concentrated on the eastern flank of the caldera operated by Los Alamos National Laboratory (House & Roberts, 2020), but the data are not openly available and more spatially balanced coverage would be advantageous for detecting microseismicity beneath the caldera. Additionally, due to the lack of ground-based geodetic instruments, it is unclear whether the apparent lack of seismogenic deformation at Valles Caldera is accompanied by slower or ductile strain.

Conclusion
We have conducted the first local Vs tomography beneath Valles Caldera using ambient noise Rayleigh tomography with a dense linear array. A prominent low-Vs anomaly is focused within the 1.23 Ma caldera's ring fracture. It exhibits a Vs reduction of ∼32% and absolute Vs of ∼2-2.1 km/s at depths of ∼3-10 km coinciding with the depths of rhyolite storage for past eruptions. The upper crustal Vs reductions in the magma reservoir beneath Valles Caldera are similar or more severe than those at systems with more abundant evidence of seismicity or surface deformation. Our results indicate the potential importance of improved hazard monitoring capacity at Valles Caldera and, more generally, affirm that even seismically quiescent volcanic fields should be regarded as potential hosts of magma in the upper crust.

Data Availability Statement
The ambient noise cross-correlation functions and inter-station phase velocities used for tomography are openly archived (https://zenodo.org/record/7466411). The raw seismic data for network 4E in 2019 are openly available from the IRIS Data Management Center (Schmandt, 2019; https://doi.org/10.7914/SN/4E_2019). Coordinator Robert Parmenter, for facilitating land access for seismic data collection. Field workers are thanked for their efforts in collecting quality data. Fraser Goff, Peter Roberts, and Tobias Fischer engaged in helpful discussions throughout the study. This manuscript also benefited from helpful reviews by Luca De Siena and Arne Spang. This research was supported by NSF Grants EAR-2113315 (B. Schmandt) and EAR-2113367 (J. Chaput). Some of the seismic equipment used in the study was provided by the IRIS PASSCAL Instrument Center and all the raw seismic data are archived at the IRIS Data Management Center; both facilities are supported through the Seismological Facilities for the Advancement of Geoscience (SAGE) Award of the NSF under Cooperative Support Agreement EAR-1851048.