Crustal structure and thickness along the Yellowstone hot spot track: Evidence for lower crustal outflow from beneath the eastern Snake River Plain

Receiver functions from seismic stations about the Yellowstone hot spot track are migrated to depth using a Vp/Vs map constructed from stacking of the direct and free surface Moho reverberations (i.e., H‐K analysis) and a shear velocity tomogram constructed from surface wave measurements. The thickest crust (48–54 km) resides in the Wyoming province beneath the sampled Laramide age blocks, and the thinnest crust (32–37 km) resides in the Montana Basin and Range province. The eastern Snake River Plain (ESRP) crust is thickest (47 km) at its NE end beneath the young calderas and thinnest (40 km) at its SW end beneath the older Twin Falls caldera. Two ESRP crustal thickness domains are found: (1) at the older Twin Falls and Picabo calderas, the mean ESRP crust is 4 km thicker with respect to its margins and (2) adjacent to the Heise caldera field, the mean ESRP crust is 4 km thicker with respect to its SE margin crust but no thicker with respect to its NW margin crust. This lobe of anomalously thick crust is explained as resulting from lower crustal outflow from beneath the Heise caldera field. Confirmation of these crustal thickness variations is provided by inspection of common conversion point (CCP) stacks that delineate several secondary features: the top of a thick high‐velocity (3.9 km/s) lower crust layer within the Wyoming province up to 17 km thick and a paired negative and positive amplitude arrival at 12 km depth and 18 km depth beneath the Yellowstone Caldera. This paired arrival would be consistent with a low‐velocity zone perhaps associated with magma staging beneath the caldera. Our most important finding is that the magmatic loads injected into the ESRP crust over the last 4–12 Myr, in tandem with the ESRP crustal viscosity structure, have been sufficient to drive significant outflow of the ESRP lower crust.


Introduction
[2] The Yellowstone hot spot track manifests a sequence of calderas that begin 16.9 Myr ago near the tristate region of Oregon-Idaho-Nevada with the calderas propagating to the NE (Figure 1). The sequence of calderas trend to the NE from the Bruneau-Jarbidge (12.7-10.5 Ma) to the Twin Falls (10.5-8.6 Ma) to the Picabo (10.2-9.2 Ma) to the Heise (6.6-4.4 Ma) caldera fields. All together, these calderas reside within the structural downwarp termed the eastern Snake River Plain (ESRP) [Perkins and Nash, 2002;Bonnichsen et al., 2007;Anders, 2009;Rodgers and McCurry, 2009;Leeman et al., 2008]. The most recent 2.1-0.6 Ma Huckleberry Ridge, Island Park, and Yellowstone calderas reside primarily upon the Yellowstone Plateau. Explanation of the ESRP downwarp requires knowledge of the crustal magma injection volumes and composition, time-integrated extension, and time-integrated lower crustal flow fluxes [McQuarrie and Rodgers, 1998;Stachnik et al., 2008;Rodgers and McCurry, 2009]. The largest uncertainty with respect to calculating an ESRP crustal mass balance is knowledge of the pre-hot spot magmatism and extension and the crustal thickness and density. The pre-hot spot conditions can only be roughly estimated from the tectonic history of this region [Hamilton, 1989;Dickinson, 2006;Foster et al., 2006]. However, the modern day crustal thickness and density structure can be estimated from the combination of the Earthscope Transportable array and previous PASSCAL seismic data.
[3] Previous petrologic and geochemical modeling suggests that 8-14 km of basaltic magma was injected primarily into the 8-18 km depth range beneath the individual calderas [Bonnichsen et al., 2007;Hanan et al., 2008;McCurry and Rodger, 2009;Leeman et al., 2008]. Heat budget modeling of the duration and volumes of the Rhyolite eruption volumes also requires about 10 km of basaltic magma injection into the mid to upper crust [Leeman et al., 2008]. The midcrust is the preferred Mixing-Assimilation-Hybridization (MASH) region where the Rhyolitic liquids are distilled via fractionation and modest levels of crustal assimilation [Hildreth and Moorbath, 1988]. Geophysical evidence for the fossil (crystallized and cooled) midcrustal MASH zone beneath the ESRP derives from gravity analysis [Sparlin et al., 1982;Lowry and Smith, 1995;Lowry et al., 2000;DeNosaquo et al., 2010], seismic refraction analysis Sparlin et al., 1982], and local earth-quake analysis [DeNosaquo et al., 2010]. In addition, our surface wave seismic tomogram finds an approximately 10 km thick high-velocity layer in the ESRP midcrust with the top of the highvelocity layer at 15-25 km depth [Stachnik et al., 2008].
[4] As the magma injections beneath the calderas cool on 1-2 Myr time scales [Anders and Sleep, 1992], the midcrustal sill complex (MCS) is predicted to become denser than the surrounding country rocks [McCurry and Rodger, 2009;DeNosaquo et al., 2010]. This high-density MCS would thus become a positive load on the lower crust which could force the lower crust to flow outward on million year time scales if the crustal viscosity is sufficiently low [Buck, 1991;McQuarrie and Rodgers, 1998;Royden et al., 2008]. The ESRP can thus be viewed as a magmatic time machine that loads the lower crust; this factor permits assessment of potential lower crustal flow if accurate crustal thickness maps are available.

Data and Methods
[5] The broadband seismic recording analyzed ( Figure 1) derive from Earthscope Transportable array data, six three-component short-period University of Utah Seismic Network seismometers in Yellowstone Park and five dominantly broadband PASSCAL seismic experiments: the 1993 eastern Snake River Plain line array [Saltzer and Humphreys, 1997], the N-S and NW-SE Deep Probe passive line arrays [Dueker and Yuan, 2004], the 2000-2001 Yellowstone array [Yuan and Dueker, 2005], the 1999-2000 Billings array  and the 2000-2001 Stanford Snake River Plain axis array [Walker et al., 2004]. Teleseismic P wave arrivals for receiver function analysis were windowed from the continuous data for all events with body wave magnitudes >5.3 and the seismic recording components rotated into the vertical, radial, tangential coordinate system. To source normalize the events, a multitaper spectral correlation method was used [Park and Levin, 2000;Helffrich, 2006]. The vertical component was used as an estimate of the source function and deconvolved with a dynamic water level [Clayton and Wiggins, 1976] derived from the preevent noise spectral amplitude. The resulting receiver functions were then culled of noisy traces by removing the 20% of the radial component receiver function with RMS amplitudes greater than three times the mean of the data set. In addition, a visual inspection was done to remove dead and/or harmonic traces. Finally, all the radial receiver functions were linearly stacked with moveout corrections to provide 134 station stack traces ( Figure 2).
[6] For each of the individual stations, a crustal thickness and bulk crustal V p /V s analysis (i.e., H-K analysis) is performed Zhu and Kanamori, 2000] to constrain the mean crustal thickness and V p /V s value in the seismic sampling cone beneath each station ( Figure 3). The weighting of S wave arrivals converted at the Moho in the H-K stacks is 0.6 (direct arrival), 0.3 (2p1s free surface reverberation), and 0.1 (2s1p free surface reverberation). The shear velocity at each station was specified to be the 1-D velocity profiles which were extracted at each station point (Figure 4) from the shear wave velocity tomogram [Stachnik et al., 2008]. The V p /V s and crustal thickness marginal probability density functions were estimated using bootstrapping with replacement [Efron and Tibshitani, 1986]. The probability functions derived from the bootstrapping are generally peaked unimodal functions (Figures 3b and 3c): indicating reasonable resolution of the trade-off between crustal thickness and V p /V s variations. The best estimate of the two model parameters was considered to be the mode of the distributions. The parameter errors were estimated from these probability functions by estimating the standard deviation about the mode of the probability function.
[7] To construct maps of the V p /V s ( Figure 5), crustal thickness and crustal thickness errors ( Figure 6), a least squares spline algorithm with a second derivative regularization term was used to interpolate the single-station results. The RMS difference between the spline predicted parameter values and the H-K measured crustal thickness and V p /V s values were small due to the spatial coherence of the single-station measurements. The crustal thickness errors (standard deviation) from bootstrapping the peak arrival time of the Moho arrival are <1.5 km ( Figure 6b). However, this crustal thickness error estimate does not include the migration velocity model uncertainties associated with the shear wave velocity model ( Figure 4) and the V p /V s map uncertainties used to migrate the receiver functions ( Figure 5). Sensitivity analysis finds that a 0.2 km/s variation in the bulk shear velocity would produce a 0.5 km variation in Moho depth [Zhu and Kanamori, 2000]. Our surface wave based shear velocity tomogram is a wellresolved image and has errors <0.2 km/s for the mean crustal velocity [Stachnik et al., 2008]; thus, the crustal thickness errors associated with our shear velocity tomogram migration velocities is estimated as <1 km. Sensitivity analysis also shows that a 0.1 change in V p /V s will produce about a 4 km change in Moho depth [Zhu and Kanamori, 2000]. Our V p /V s error map ( Figure A1) finds that the standard deviation of our bootstrapped V p /V s measurements is <0.03; thus, crustal thickness errors associated with our V p /V s map are estimated as <1.3 km. Assuming the three sources of error are independent, the maximum uncertainty in our crustal thickness maps is estimate as <3.7 km, more typically <2.5 km.
[8] The common conversion point (CCP) receiver function images [Dueker and Sheehan, 1997] were constructed using a three dimensional pixel parameterization with 2 km thick layers and 40 by 40 km wide CCP bins (Figures 8-11). The bin center points were spaced every 10 km so that adjacent bins had 75% data overlap, but no data overlap at three bin offsets. Time was mapped to depth using our shear velocity tomogram [Stachnik et al., 2008] and the V p /V s map ( Figure 5) for the crust. For the mantle below the Moho, the IASPEI91 velocity model mean upper mantle V p /V s value of 1.81 was used.

Single-Station Interpolated Maps
[9] As our prior seismic results reported [Stachnik et al., 2008], the mean crustal velocity map generally shows high (>3.6 km/s) mean velocity Wyoming province crust and low (<3.5 km/s) mean velocity Yellowstone caldera crust ( Figure 4). The high-velocity Wyoming crust is primarily due to a high-velocity lower crustal layer imaged by the shear wave tomogram ( Figure 11). This layer has been previously imaged as the so-called 7.x magmatic underplate layer (i.e., with a velocity >7.0 km/s) from active source studies [e.g., Gorman et al., 2002]. In addition, the CCP images find a positive arrival from the top of this high-velocity layer where the station density is highest within the Billings, Montana array (Figures 8, 9, and 11).
[10] The V p /V s map ( Figure 5) shows a reasonable range of 1.76-1.86 (ignoring the map edge values). The mean Vp/Vs is 1.81 compared to global estimate mean continental crustal value of 1.79 [Zandt and Ammon, 1995;Christensen, 1996]. A typical quality H-K analysis is shown in Figure 3 which finds a Vp/V s value of 1.79 ± 0.02. The maximum error in the V p /V s map is 0.03 as found by bootstrapping the H-K analysis. The principle anomaly observed in the V p /V s map is the relatively high values (>1.84) along the ESRP and normal values within the Yellowstone Park.
[11] The crustal thickness map (Figure 6a) [Zeiler et al., 2005], beneath the sampled Idaho Batholith [Kuntz et al., 2005], and to the south of the central ESRP. The Yellowstone Plateau region has crustal thick- ness of 47-52 km; this thick crust primarily manifests the Laramide age shortening associated with the Beartooth Mountains and the magmatic underplate that created the high-velocity lower crust beneath much of the Wyoming Province crust.

Common Conversion Point Images
[13] The CCP images (Figures 8-11) find that the direct Moho arrival is well imaged by migration of our radial receiver function data set. The ESRP/ Yellowstone Plateau parallel cross section (Figure 8, cross section A) shows the Moho thickening to the NE toward the Yellowstone Plateau and thick crust within the Wyoming province. At the eastern end of this cross section, the top of the Wyoming Province high-velocity lower crustal layer is imaged at 35 km depth with the Moho at 52 km depth. The NW-SE cross section through the Yellowstone caldera (Figure 8, cross section B) shows a sharp change in crustal thickness at the NE corner of Yellowstone Park from thin crust beneath the Montana Basin and Range to thick crust beneath the Yellowstone caldera and the Wyoming province. The ESRP perpendicular cross section across the Picabo caldera (Figure 8, cross section C) shows thick Wyoming province crust at the SE end of the image and thin crust beneath the Montana Basin and Range province. The ESRP crust beneath the Picabo caldera is found to be seismically transparent with upper crustal structure outside the ESRP being truncated at the ESRP margins. In general, the ESRP Moho is depressed by 2-4 km with respect to the adjacent NW and SE margins (see also Figure 7).
[14] The E-W Montana/northern Idaho cross section (Figure 8, cross section D) shows the greatest crustal thickness variation from 52 km beneath eastern Montana to 35 km near the Eocene age Bitterroot detachment and granitic batholith in northern Idaho [Foster et al., 2001]. The top of the high-velocity lower crustal layer is also found at the east end of this cross section. The N-S Wyoming/Montana cross section (Figure 8, cross section E) shows the thick (48-52 km) crust beneath the region shortening during the Laramide orogeny and not affected by late Cenozoic extension [Dickinson, 2004]. The top of the highvelocity lower crustal layer is imaged north of 44.5°latitude beneath the Billings array ( Figure 1). Two cross sections through the Billings array ( Figure 9) show the direct Moho arrivals from the top of the high-velocity lower crustal layer at 28-34 km depth. The thickness of this lower crustal layer is up to 17 km thick with the layer thinning to zero thickness at the NE and SW end of cross section B.
[15] A final notable feature in the CCP images is a paired positive and negative amplitude arrival at 12 and 18 km depth beneath the Yellowstone caldera (Figures 8 (cross sections A and B) 10). This paired arrival would be consistent with a low-velocity  zone: the 12 km negative polarity arrival manifesting a negative velocity gradient and the 18 km positive polarity arrival manifesting a positive velocity gradient. This paired arrival is directly under the two most volcanically active regions of the Yellowstone caldera between the Mallard Lake and Sour creek resurgent Rhyolitic domes [Lowenstern and Hurwiltz, 2008]. A similar finding of an upper crustal low-velocity zone is found by waveform modeling of teleseismic S wave data from three   4. Discussion

Lower Crustal ESRP Outflow
[16] Our most important new result is our new crustal thickness map which provides a synoptic scale sufficient to assess the mass balances associated with the magmatic inflation of the ESRP crust and potential lower crustal outflow forced by the midcrustal densification caused by this magmatic load. If the hypothesis that magmatic injection along the ESRP has stimulated crustal flow is correct, then the time history of the crustal flow should be correlated with the magmatic injection history and the postcaldera lower crustal thermal evolution which controls viscosity. Given the high (>90 mW/m 2 ) heat flow along the ESRP and its margins [Blackwell and Richards, 2004] and the finding of low flexural rigidity [Lowry and Smith, 1995;Lowry et al., 2000], it seems plausible that the lower crust is capable of flow driven by the midcrustal sill load. Estimates of lower crustal flow rates of 1-7 cm/yr (10-70 km/Myr) in regions of elevated heat flow have been proposed where sufficient pressure gradients in the lower crust exist [Buck, 1991;Royden et al., 1997Royden et al., , 2008. Dividing the 50-80 km lateral extent of the anomalously thick crust to the NW of the Heise caldera field by the mean age of the Heise caldera field (5 Ma) provides a maximal lower crustal flow rate of about 1 cm/yr (10 km/Myr).
[17] The simplest crustal thickness evolution scenario for the ESRP and its margins assumes that prior to formation of the hot spot track, the crustal thickness and density structure was uniform perpendicular to the ESRP between the Twin Falls and Heise calderas and its adjacent margins. In addition, post-hot spot track integrated crustal dilatation is also assumed to be relatively uniform perpendicular to the ESRP. Estimates of net dilatation at the NW and SE ESRP margins is estimated as 15% and 25% [Rodgers et al., 2002]. Thus, to first order, we assume that the ESRP initial conditions and integrated dilatation perpendicular to the ESRP have been relatively uniform. Given these assumptions, the only change in the ESRP crustal thickness with respect to its margins would result from mass additions due to magmatism and sedimentation or mass subtractions associated with caldera eruption atmospheric plumes. The sedimentation addition is estimated as <1 km and the atmospheric losses as <0.75 km [Rodgers and McCurry, 2009]. These two effects are opposite in sign and nearly cancel; thus, the dominant ESRP mass variable is the caldera forming crustal magmatic injections. In summary, the expected thickening of the ESRP crust with respect to its margins is dictated by the petrologic and caldera heat budget constraints that suggest 8-14 km of mantle derived basaltic magmas are required to fuel the calderas [Hanan et al., 2008;McCurry and Rodger, 2009;Leeman et al., 2008].
[18] The shear velocity tomogram finds a highvelocity (>3.6 km/s contour) midcrustal layer within the ESRP that starts at the NE end of the ESRP beneath the Island Park caldera and extends to the edge of our surface wave velocity sampling at the Picabo caldera (Figures 11a and 11b). Beneath the Heise caldera field, the 3.6 km/s shear velocity contour that outlines this layer resides between18-32 km depth (14 km thick). At the SE end of the ESRP, the top of this high-velocity midcrustal layer deepens by about 5 km. This layer is thought to manifest the midcrustal sill (MCS) complex that forms the mafic magmatic "tanks" from whence the Rhyolitic caldera magmas were derived via fractionation, assimilation, and hybridization. Beneath this MCS layer, a relatively low velocity lower crustal "wedge" is found that is thickest beneath the Island Park caldera and pinches out at the SW edge of the Heise caldera field.
[19] The ESRP perpendicular cross section (Figures 11c and 11d) shows two interesting features: the MCS layer dips to the NW and the lowvelocity (<3.6 km/s) lower crust extends up to 50 km to either side of the ESRP margins. With respect to the mantle velocity anomalies shown in the tomograms, the ESRP-YP physiographic province is underlain by very low (4 km/s) shear velocities at 75-125 km depth with a relatively high velocity mantle lid between the Moho and 75 km depth (Figures 11a and 11b). The low ESRP mantle velocity rapidly changes beneath the Beartooth Mountains to cratonic lithospheric velocity values of 4.9 km/s [Artemieva, 2009;Bedle and van der Lee, 2009]. In the ESRP perpendicular cross section (Figures 11c and 11d), the low-velocity ESRP mantle is found directly beneath the 90 km wide ESRP with higher-velocity mantle lithosphere beneath SW Wyoming (SE end) and western Montana (NW end).
[20] Based on the above observations, we believe a good circumstantial case is made for the flow of magmatically thickened lower crust from beneath the Heise caldera field into the adjacent NW margin crust. However, this statement begs the question of where the magmatically thickened ESRP crust has flowed from beneath the older Twin Falls and Picabo caldera fields. Two differences are noted between the Heise caldera field and the older Twin Falls and Picabo calderas. First, these two older calderas are spaced farther apart along the ESRP (Figure 1) indicating less mass flux into the ESRP per unit area. Second, these caldera fields have had more integrated dilatation with respect to the Heise caldera field because extension associated with caldera formation began earlier in this region [Anders et al., 1989;Pierce and Morgan, 1992]. In addition, the proximity of these older calderas to the concentrated extension of the western Snake River Plain graben [Cummings et al., 2000] and the Northern Nevada Rift [Glen and Ponce, 2002] is noted; these regions of concentrated extension could create lateral pressure gradients in the lower crust that would promote the flow of lower crustal mass from the magmatically thickened calderas areas.

Yellowstone Caldera Low-Velocity Zone
[21] Modern day gas fluxes near recent intracaldera basaltic eruptions [Lowenstern and Hurwiltz, 2008] and deformation monitoring Puskas et al., 2007] suggest that ongoing post 0.6 Ma Yellowstone caldera magmatic activity is occurring. The paired negative/positive amplitude arrivals at 12-18 km depth found beneath the Yellowstone caldera by the CCP images (Figures 8 (cross  sections A and B) and 10) are consistent with other geophysical data that suggest the Yellowstone caldera low-velocity and low-density anomalies extend to about 20 km depth: the low velocities beneath the caldera in our shear wave tomogram (Figures 11a  and 11b) and gravity modeling [DeNosaquo et al., 2010]. A tomogram constructed from measured local earthquake traveltimes finds a low (5.4 km/s) P wave velocity at 8 km depth, but cannot resolve structure below 12 km depth [Husen et al., 2004]. Waveform modeling of S-P precursors from three broadband stations near the Mallard Lake dome requires a substantial low-velocity zone at 5-8 km with a thickness of 3-5 km [Chu et al., 2009]. Thus, our finding of a 6 km thick low-velocity zone with its top at 12 km depth provides a depth range over which basaltic magma is being staged beneath the Yellowstone caldera. In the near future, this finding can be tested via higher-resolution ambient noise surface wave imaging with the deployment of new broadband seismometers within Yellowstone Park in 2010.