NEW CONSTRAINTS ON THE MAGMATIC SYSTEM BENEATH NEWBERRY
VOLCANO FROM THE ANALYSIS OF ACTIVE AND PASSIVE 
SOURCE SEISMIC DATA AND AMBIENT NOISE
by
BENJAMIN ALLEN HEATH
A THESIS
Presented to the Department of Geological Sciences
and the Graduate School of the University of Oregon
in partial fulfillment of the requirements
for the degree of
Master of Science 
December 2014
THESIS APPROVAL PAGE
Student: Benjamin Allen Heath
Title:  New Constraints on the Magmatic System beneath Newberry Volcano from the 
Analysis of Active and Passive Source Seismic Data and Ambient Noise
This thesis has been accepted and approved in partial fulfillment of the requirements for 
the Master of Science degree in the Department of Geological Sciences by:
Douglas Toomey Chairperson
Emilie Hooft Member
Gene Humphreys Member
Ilya Bindeman Member
and
J. Andrew Berglund Dean of the Graduate School
 
Original approval signatures are on file with the University of Oregon Graduate School.
Degree awarded December 2014
ii
© 2014 Benjamin Heath
iii
THESIS ABSTRACT
Benjamin Heath
Master of Science
Department of Geological Sciences
December 2014
Title: New Constraints on the Magmatic System beneath Newberry Volcano from the 
Analysis of Active and Passive Source Seismic Data and Ambient Noise
Using joint P-wave seismic tomography, receiver functions and ambient noise we 
image the magmatic structure beneath Newberry Volcano, located near Bend, Oregon. 
Use of active source and teleseismic events in a joint tomographic inversion provides the 
ray crossings necessary to resolve a low velocity body around 4 km depth. Receiver 
functions show large lateral heterogeneity and are consistent with the location of a low 
velocity body derived from the tomography but require a larger low velocity anomaly. 
Ambient noise autocorrelations are used to image a low velocity reflector, located at ~3 
km depth, shallower than the imaged low velocity body recovered using tomography and 
receiver functions. Ultimately, our results reveal a magma chamber at 3-4 km depth 
beneath Newberry caldera, with an overlying partially molten sill at ~3 km depth. These 
results show the usefulness of dense seismometer deployments over volcanoes.
iv
CURRICULUM VITAE
NAME OF AUTHOR:  Benjamin Heath
GRADUATE AND UNDERGRADUATE SCHOOLS ATTENDED:
University of Oregon, Eugene
Northwestern University, Evanston, IL
DEGREES AWARDED:
Master of Science, Geological Sciences, 2014, University of Oregon
Bachelor of Arts, Earth and Planetary Sciences, 2012, Northwestern University
AREAS OF SPECIAL INTEREST:
Seismology
GRANTS, AWARDS, AND HONORS:
Johnson Fellowship, University of Oregon, 2012
v
ACKNOWLEDGMENTS
I would like to thank Doug Toomey and Emilie Hooft for their guidance. I would 
also like to thank my friends and family for their support. 
vi
TABLE OF CONTENTS
Chapter Page
I. INTRODUCTION.................................................................................................... 1
II. TOMOGRAPHY..................................................................................................... 3
Introduction ........................................................................................................... 3
Geologic Setting .................................................................................................... 5
Previous Geophysical Studies ............................................................................... 8
Data ....................................................................................................................... 9
Methods ................................................................................................................. 11
Results ................................................................................................................... 15
Teleseismic Delay Times ................................................................................ 15
Joint Teleseismic and Active Source Tomographic Inversion ........................ 18
Discussion ............................................................................................................. 24
Conclusion ............................................................................................................ 34
III. RECEIVER FUNCTIONS .................................................................................... 36
Introduction ........................................................................................................... 36
Methods ................................................................................................................. 38
Results ................................................................................................................... 40
Discussion and Conclusion ................................................................................... 51
IV. AMBIENT NOISE ................................................................................................ 52
Introduction ........................................................................................................... 52
Methods ................................................................................................................. 53
Results ................................................................................................................... 56
vii
Chapter Page
Discussion ............................................................................................................. 62
Conclusion ............................................................................................................ 67
V. CONCLUSION ...................................................................................................... 68
REFERENCES CITED................................................................................................ 69
viii
LIST OF FIGURES
Figure Page
2.1. Map of Newberry showing seismic stations and explosive sources .................... 7
2.2. Plot stations from the 2008 experiment................................................................ 10
2.3. Plot of event locations used in our teleseismic tomography ................................ 12
2.4. Plot of the various components for the P arrival on Event 1102 ......................... 13
2.5. Plot of ray paths for an active source event and a teleseismic earthquake .......... 14
2.6. Plot of delay times by back azimuth..................................................................... 17
2.7. Plot of mean delay times by back azimuth .......................................................... 18
2.8. Relative delay times modeled for low and high frequency sources ..................... 19
2.9. Cross sections through our tomographic model at various depth ranges ............. 20
2.10. Plot of recovered tomographic velocities along line A to A' .............................. 23
2.11. Plot of various synthetic models......................................................................... 25
2.12. Solidification times for different sills ................................................................ 33
2.13. 3-D plot of tomographic model .......................................................................... 34
3.1. Plot of locations of events used for receiver functions ........................................ 39
3.2. Plot of radial receiver functions and transverse receiver functions for events    
coming from the North West; plot of vertical component receiver function........ 43
3.3. Plot of radial and transverse receiver functions for events from the South East;       
plot of vertical component receiver function....................................................... 44
3.4. Radial receiver function for event 1102 ............................................................... 45
3.5. Model with a thin magma chamber; predicted radial and vertical receiver       
functions .............................................................................................................. 48
ix
Figure                                                                                                                             Page
3.6. Model with a thick magma chamber; predicted radial and vertical receiver   
       functions ............................................................................................................... 49
3.7. Radial/vertical receiver functions for source frequency of 0.75 Hz ....................    50
4.1. Mean autocorrelation for station NWB35 ............................................................ 55
4.2. Plot of active source traces .................................................................................. 57
4.3. Stack of autocorrelations for all stations .............................................................. 58
4.4. Stack of autocorrelation for all stations using different filter frequency.............. 58
4.5. Amplitude relative to frequency for station NWB02 ........................................... 59
4.6. Autocorrelation for active source data ................................................................. 60
4.7. Grid searched fit to autocorrelations .................................................................... 62
4.8. Autocorrelation for event 1102............................................................................. 63
4.9. Synthetic amplitudes for the autocorrelation........................................................ 66
x
1  
CHAPTER I INTRODUCTION  The subsurface structure of Newberry Volcano, located south of Bend Oregon, has remained enigmatic despite intense study. Newberry, which lies at the Northwest edge of the Basin and Range, appears to be strongly influenced by several zones of extensional faulting [e.g. Gettings and Griscom, 1988], though the influence of the faults on the magmatic system, especially beneath the caldera, remains uncertain. Several fundamental questions, such as magma chamber geometry, location, and partial melt percentage also remain poorly constrained. Earlier studies imaged a low velocity body consistent with a shallow magma chamber at ~4 km depth [e.g. Achauer et al., 1988] and recently Beachly et al. [2012] were able to model a secondary arrival to constrain partial melt percentage and volume of the magma chamber ( > 20%, 2-8 km3) as well as its rough location. However, neither study was able to image the full 3-D structure of the magmatic system.  For volcanic systems in general there are a number of competing hypothesis for magma chamber structure ranging from spherical melt-filled chambers, to isolated, thin, melt lenses [e.g. Detrick et al., 1987]. These models imply different constraints on geochemical make-up and evolution of the melt.  Here we seek to illuminate the magmatic system beneath Newberry, specifically beneath the caldera, using several different seismic methods. First, we use a combination of active source travel time and teleseismic earthquake delay time measurements to create a 3-D tomographic model of the velocity structure beneath Newberry Volcano. Second, we use receiver functions to 
2  
image the subsurface and map abrupt velocity variations. Finally, we use ambient seismic noise to create reflection profiles of the subsurface that supplement the observations derived from receiver functions. By combining multiple methods, we are able to create a more realistic image of the subsurface magmatic system beneath this continental volcano.  
3  
CHAPTER II TOMOGRAPHY  Introduction:  Despite intense effort, the geometry and structure of magmatic systems remains poorly constrained. Magmatic models range from spherical melt filled chambers to melt sills [e.g. Detrick et al., 1987], to largely crystallized mush regions [e.g. Eppich et al., 2011]. Variability in models is partially due to variations in magmatic composition and its effect on magma viscosity, but also due to limited geophysical resolution capabilities. Large variation in geochemical models is due in part to the imperfect resolution of these geophysical experiments. Improved resolution of subsurface structure will lead to better estimates of melt within magmatic systems as well as better geochemical models of magmatic processes. These in turn will lead to a better understanding of the whole magmatic system and enhanced predictive capacity for eruption.    While exact geometry of magmatic systems remain poorly constrained, there is a well-known link between volcano placement and tectonic structures.  At Mount St. Helens a fault inferred from seismicity appears to run beneath the volcano [e.g. Waite and Moran, 2004; Weaver et al., 1987]. At Medicine Lake the shield volcano edifice seems to correlate well with the intersection of multiple fault zones [Donnelly-Nolan, 1987], with the faults providing the pathway necessary for basalt to reach the surface with minimal crustal contamination. Similarly, Newberry volcano, also a shield volcano, appears to lie at the intersection of multiple fault zones [Fitterman, 1988]. These studies suggest a tectonic influence on volcano location and potential influence on magma chamber 
4  
geometery.  Tomographic imaging of magma chambers helps place constraints on their size and location, information which is critical to understanding magmatic processes [e.g. Husen, 2004; Lees, 1992; Sinton and Detrick, 1992]. Despite these efforts, accurate resolution of magnitude of anomalies and hence precise measurement of bulk physical properties remains difficult, in part due to wavefront healing [e.g. Nolet and Dahlen, 2000], smoothing constraints and imperfect ray coverage. Large geophysically resolvable magma chambers remain relatively rare [e.g. Lees, 2007], suggesting that these chambers might either be short lived, or of smaller percentage melt than previously believed.   Both teleseismic and active source data sets have been used to tomographically image magmatic systems [e.g. Ritter and Evans, 1996; Beachly et al., 2012; etc], though traditionally not together in the same inversion. Teleseismic tomographic methods have commonly imaged mantle and lower crust structure and not shallower crustal structure because of the lack of ray crossings in the crust and relatively large station spacing of these experiments. These studies usually are limited to lower frequency (< 1 Hz) arrivals relative to active source data (~10 Hz). In contrast, active source tomography experiments often resolve upper crustal structures, with the depth of imaging being limited by the turning depth of the rays. The higher frequency content, higher ray and station density, and more numerous ray crossings, all contribute to the greater resolution power of active source tomography in the upper crust.   The utilization of a dense seismic array allows for recording of both active source and teleseismic passive source events on the same stations. Studies on dense arrays using teleseisms have been used to infer crustal structure through vertical-component receiver 
5  
functions [Schmandt and Clayton, 2013] and in coupled tomography experiments using local, regional and teleseismic seismicity [Biryol et al., 2013]. This shows that the decreased station spacing of small scale active and passive source experiments increases teleseismic imaging capacity due to the added ability to interpret lateral variations in these waveforms.  In this study we use a coupled tomographic method that combines active source and teleseismic datasets to investigate the magmatic structure beneath Newberry Volcano. By combining active source P-wave travel times with teleseismic P-wave delay times we can make use of the crossing of the rays allowing for better resolution and a better constraint on the size, location and partial melt volume of the inferred magma chamber beneath Newberry.  The relative strengths and weaknesses of each method complement each other, increasing seismic resolution with little added effort. While this tomographic method suffers from several of the short comings of other tomographic methods, the complications due to coupling the datasets are justified a posteriori.   Geologic Setting:  Newberry Volcano, located in Central Oregon ~60 km east of the High Cascades (Figure 2.1), lies at the end of a northwest trending chain of bimodal volcanism in the High Lava Plains which began ~16 Mya [Fitterman, 1988; Jordan et al., 2004]. Newberry itself rises 1.1 kilometers above the regional elevation. Its dimensions are roughly 60 km by 30 km, north-south and east-west, respectively and covers a rough area of 1300 km2 with a total eruptive volume of 500 km3 [Jensen et al., 2009; Macleod, 1995].  In contrast to the edifice of the volcano, which is elongated in the north-south 
6  
direction, the caldera is elongated in the east-west direction which is likely the result of regional extension along a north-south trending fault [e.g. Acocella et al., 2004].   Newberry Volcano lies at the intersection of the Brothers, Sisters and Walker Rim fault zones (Figure 2.1). The Brothers fault marks the northern boundary of Basin and Range tectonics [Lawrence, 1976]. The Sisters fault zone extends northward from Newberry [Fitterman, 1988], while the Walker Rim fault zone is located south of the volcano and appears to intersect the volcanic edifice [Fitterman, 1988]. Previous studies suggest the fault zones are connected underneath Newberry [MacLeod and Sherrod, 1988], thus explaining the general shape of the caldera.  At Newberry during the Holocene, there have been numerous basaltic andesite and rhyolite eruptions. Rhyolitic eruptions occurred in the caldera, primarily at locations on or in close proximity to the ring fractures [MacLeod and Sherrod, 1988]. In contrast, basaltic andesite eruptions have primarily occurred on the flanks of the volcano. Typical repose times for rhyolitic eruptions are on the order of 2,000-3,000 years [MacLeod and Sherrod, 1988]. The most recent eruption, the Big Obsidian Flow, was 1300 years ago, with a total eruptive volume of 0.16 km3. Combined Holocene rhyolitic eruptive output totals about 1 km3 [MacLeod and Sherrod, 1988]. Holocene erupted rhyolite is predominately aphyric with similar chemistry between eruptions [MacLeod and Sherrod, 1988]. These observations, combined with the observation of a consistent repose time has  
7  
Figure 2.1: Map of Newberry showing seismic stations and explosive sources. Elevation is contoured and shaded by the gradient. Triangles denote stations from the 2008 experiment (red), and stations from previous experiments (blue). The red stations recorded both active source as well as teleseismic events. Blue stations recorded only active source events. Yellow stars denote locations of all active source events (explosive sources) in this study. The Sisters, Brothers and Walker Rim fault zones [Fitterman, 1988] are shown by the Green shaded areas. The dashed box shows the outline for Figure 2.2. The orange lines mark caldera faults.     
8  
 led to the inference that there is a long-lived molten rhyolitic magma chamber beneath Newberry throughout the Holocene. The absence of Holocene basaltic volcanism inside the caldera, but frequent basaltic volcanism outside of it suggests a magma chamber shadow effect [MacLeod and Sherrod, 1988]. This is where basaltic volcanism can erupt around a felsic magma chamber, but not through it, leading to the surface expression of a zone of rhyolitic eruptive material encompassed by basalt. Basaltic volcanism may facilitate rhyolitic eruptions through magmatic underplating.  Using this interpretation, one predicts a felsic magma chamber with a maximum diameter of 5 km at shallow depths underneath Newberry caldera [Fitterman, 1988].   Previous Geophysical Studies:  There have been a number of geophysical studies conducted at Newberry. Using travel times from active sources, Achauer et al., [1988] tomographically imaged low velocities inside the caldera at ~3 km depth, which they interpreted as a magma chamber. Using the same data set as Achauer et al. [1988], Zucca and Evans [1992] inverted the data for seismic attenuation. They suggested that, due to the lack of a high attenuation zone in the location of the Achauer et al. [1988] low velocity magma chamber, the interpreted magma chamber either was part of a breccia pipe, or was a dry, solidified and cracked pluton. Beachly et al. [2012] combined the data set of Achauer et al. [1988] with a dense (~0.8 km station spacing) linear array of 81 seismic stations deployed to record a shot from the High Lava Plains experiment [e.g. Cox et al, 2013] and tomographically imaged the volcano (Figure 2.1 and Figure 2.2). Magma chambers ranging from 0 km3 melt to 60 km3 of melt were all consistent with their tomographic results. Utilizing a 
9  
secondary arrival from the shot, they were able to show that their observations were consistent with only 2 to 8 km3 of melt, a substantially smaller amount than what the tomography could uniquely resolve.   In addition to the active source studies, Stauber et al. [1988] inverted teleseismic delay times to image the magmatic system.  They did not observe a magma chamber, but noted that the dimensions of melt would just need to be smaller than the dominant seismic wavelength (~10 km) and a magma chamber was therefore not inconsistent with their results. Electric resistivity and gravity studies failed to conclusively detect a crustal magmatic system, though none of these studies exclude a magma chamber from existing beneath Newberry [Fitterman et al., 1988; Gettings and Griscom, 1988].  Data:  Our tomographic analysis uses a combination of active source travel time data and delay time data from teleseismic earthquakes. The combination of lower frequency (<1 Hz), subvertically incident, teleseismic data with higher frequency (~10 Hz), subhorizontally propagating active source data improves three-dimensional sampling of structure.  The active source data are from 3 previous seismic experiments [Cotton and Catchings, 1989; Dawson and Stauber, 1986; Beachly et al. 2012]. Figure 2.1 shows the distribution of sources and receivers.  The travel time data set is identical to that used by Beachly et al. [2012], and is comprised of 1006 P wave arrival times and their associated uncertainties, which average ~20 msec.   
10  
Figure 2.2: Plot stations from the 2008 experiment, shown as red triangles. Plotted in green and blue are earthquakes shallower than 1 km and deeper than 1 km respectively. Earthquakes were located by PNSN stations from Dec. 2012 to Sept. 2014. The line A to A', marks an area along which we take all of our cross sections.      The teleseismic data were recorded in 2008 on an array of 81 Mark Products L-22D short period (2 Hz) seismometers deployed in a roughly linear array across Newberry Volcano, trending from the south-west to the north-east (Figure 2.2).  Station spacing was 800 m on the flanks of the volcano and 300 m in the caldera, producing a line of approximately 40 km centered on the volcano. The primary objective of deploying these seismometers was to record waveform data from a one ton, 26 m deep explosion 
11  
from the High Lava Plains Experiment. Since the array was deployed for about 2-3 weeks after the shot, it recorded good quality data for 21 teleseismic events at distances of 38˚ to 92˚ (Figure 2.3); event magnitudes were 5 to 6.4 mb. Good quality three-component data were recorded for these events, with the best signal-to-noise ratio data being observed on the vertical channel (Figure 2.4). In some instances, the amplitudes observed on the radial or transverse channels were comparable to the vertical, possibly as a result of strongly heterogeneous or anisotropic structure [Langston, 1979]. In addition, due to the dense station spacing, continuous, short wavelength variations (~600 m) in the waveform could be observed in the teleseismic waveform. Here we measure delay times using only the vertical component data. Teleseismic delay times were measured relative to a standard Earth model [Kennett and Engdahl, 1991] using the cross correlation method of VanDecar and  Crosson [1990]. We measured delays for 3 different frequency bands, using Gaussian filters with center frequencies 0.3 Hz, 0.5 Hz, and 1 Hz and half-widths 0.15 Hz, 0.2 Hz, 0.4 Hz, respectively. Delay time uncertainties for all teleseismic measurements were set at 30 ms. Below we discuss potential problems of utilizing this method when the shallow structure (< 10 km depth) is strongly heterogeneous.  Methods:  The novelty of our analysis is the simultaneous inversion of active source travel times and teleseismic delay times to obtain a tomographic image of crustal structure. The combination of long period, low frequency (< 1 Hz), subvertically incident,   
12  
             Figure 2.3: Plot of event locations shown in gray used in our teleseismic tomography. Newberry is plotted as a black triangle in the center. Event 1102 is plotted as a black star. This event has high signal to noise ratio and is shown below.    teleseismic data with short period, high frequency (> 3 Hz), subhorizontally propagating active source data significantly improves the resolution of a shallow (< 10 km depth) 3-D magmatic system beneath Newberry (Figure 2.5).  For the calculation of the travel times and ray paths of the P-waves, we deepen the Beachly et al. [2012] starting model to a depth of 10 km below sea level (bsl). The starting velocity model is one dimensional with grid spacing of 200 m and 100 m in the horizontal and vertical dimensions respectively.  Further description is available in   
13  
 Figure 2.4: Plot of the various components for the P wave arrival on Event 1102 (Figure 2.3) as mapped to A-A' (Figure 2.2). Distance of zero is the caldera center with negative distances denoting closer proximity to A and positive closer to A'.  All traces are plotted on the same scale. You can see a clear P wave arrival on the vertical. Note the ringing in the P wave code on both the radial and transverse channels.     Beachly et al. [2012]. For depths greater than the Beachly et al. [2012] starting model, we utilize the fastest velocity from their 1-D model. We use a perturbational grid with spacing of 400 m and 200 m in the horizontal and vertical dimensions respectively.  The forward problem of calculating predicted arrival times uses a graph theory based 3-D ray tracer that makes a high frequency approximation. For teleseismic travel times, we calculate travel times and paths though a 3D model to 10 km bsl and then include the travel time from the event location to the bottom of our model assuming a radial earth velocity structure. Our 3-D ray tracing incorporates elevation into the forward problem by vertically shearing the velocity mesh [Toomey et al., 1994].   
14  
Figure 2.5: Plot of tomographic cross section from Beachly et al. [2012] with contours denoting fractional change in velocity. Cross section taken from A to A' (Figure 2.2). Plotted in red are the ray paths for an active source event. Plotted in blue are the ray paths for a teleseismic earthquake. Note that the active source ray paths are predominately horizontal and the teleseismic ray paths are predominately vertical. The crossing of the ray paths increases the resolution.         Because all active source event locations fall within our model, our inversion matches the travel time for these picks. We calculate event static corrections for nine of the furthest active source events [Beachly et al., 2012] and all of the teleseismic events. We do not employ station static corrections in the inversion since we are inverting for near station structure. Below we discuss how teleseismic delay time measurements influenced by reverberations can also map as unphysical anomalies.    We jointly invert the teleseismic delay times and the active source travel times for seismic velocity following the approach of Toomey et al. [1994]. The relative weighting of the active source and teleseismic data sets is achieved through normalization of residuals by their relative uncertainties. We use all of our teleseismic picks in the inversion. This means that for each station/event pairing we can have up to 3   
15  
measurements, corresponding to the different picking frequencies. This has a tendency to favor measurements that are consistent across the 3 frequency bands.  Because the active source tomography relies on travel times for sources without event terms, we are able to constrain absolute velocities. In contrast, teleseismic tomography uses relative delay times, which can only resolve relative variations in velocity.   We calculate synthetic seismograms derived from a 2-D model using E3D [Larsen and Harris, 1993] to test the ability of our tomographic model to reproduce the character of the coda observed in teleseismic arrivals as well as to quantify the accuracy of the ray theory approximation. We use a Ricker wavelet as a source (0.5 Hz) and initiate with a plane wave incident from left to right with ray parameter of 0.05 s/km.  Results: Teleseismic Delay Times  We observe systematic variations in the teleseismic delay times at a variety of length scales.  Figure 2.6 shows that teleseismic delays vary by ~0.4 s over the 40-km-long aperture of the seismic array.  At the longer spatial scales (~10 km), delays are consistently greater on the eastern flanks of the volcano.  At a smaller scale (~1 km) the observed delays are consistently negative (early) near the caldera rim, whereas arrivals within the caldera are more positive (late).  Negative delays near the caldera rim are consistent with the results of Beachly et al. [2012], who reported high velocity anomalies in this region. Figure 2.6 also shows that teleseismic delays systematically vary on the scale of the station spacing (300 m). The observed teleseismic delays also show 
16  
differences that are both frequency dependent and that vary with event backazimuth (Figure 2.6 and Figure 2.7). Variations by backazimuth are consistent with three-dimensional variations in seismic velocity structure.  Figure 2.7 shows that the variation in teleseismic delays is greater when measured at lower frequencies, a result that is not consistent with finite frequency effects [Nolet and Dahlen, 2000].  We attribute this observation to waveform interference between the direct arrival and a reverberation from a near surface interface. To support this inference, we calculated synthetic seismograms, assuming a one dimensional background model with a low velocity region at 3 km depth (2.3 km/s). For measurements made at higher frequencies, the first arrival and reverberation off the low velocity feature are isolated from one another, allowing for easy identification of the first full period of the arrival for cross correlation. In contrast, at lower frequencies, the first full period is not only a function of the first arrival, but also of the reverberation, which may vary spatially. An example of the clear frequency dependence of this effect for reverberations off of a shallow magma chamber is shown in Figure 2.8. Because the first order delay times that give our model its distinct characteristics are observed across all three frequency bands, we hold that these artifacts are negligible.  This allows us to investigate how structures at depth are mapped to delay times measured at the surface. Ideally, the variation in travel times observed by cross-correlation of waveforms should match the predictions of ray theory. Instead, we observe a distinct pattern, which ultimately biases how we map structures at depth. Ray theory predicts a more localized anomaly, whereas the delay times obtained from cross-correlation show a broader anomaly, consistent with theoretical results derived from  
17  
Figure 2.6: Plot of delay times by back azimuth. In blue are the delay times with respect to iasp91 [Kennett and Engdahl, 1991] mapped to the line A-A'. Red denotes the mean for the given back azimuth range.   wavefront healing [Nolet and Dahlen, 2000]. A broad anomaly in delay times generated by a narrow velocity perturbation at depth will tend to map to a broad anomaly in our tomographic model. While this modeling illustrates the frequency component bias in our forward problem, we cannot use E3D for further analysis of our observed delay times because the heterogeneous three dimensional structure of Newberry does not lend itself to two dimensional approximations. We consider our tomographic image to be a low resolution image that reconstructs anomalies over a broader region, in comparison with the actual structure.  
18  
Figure 2.7: Plot of mean delay times by back azimuth colored by frequency and mapped to the line A-A'. Note that the low frequency delay times (0.3 Hz) often have higher amplitudes than the other frequencies. This may be a result of interference from a reverberation off of a magma chamber.     Joint Teleseismic and Active Source Tomographic Inversion Our tomographic image reveals several interesting low and high velocity features.  First, we recover a low-velocity zone (LVZ) beneath the caldera at about 3-5 km depth. This LVZ represents a ~10%  reduction in velocity (~ 0.6 km/s), is ~1 km thick, and has  
19  
Figure 2.8: Above: Plot of model velocity shown with a magma chamber. Arrows denote approximate incidence angle of modeled teleseismic wave. Bottom: Relative delay times modeled for low frequency source (0.3 Hz, red), for the mid-range frequency source (0.5 Hz, blue) and for the high frequency source (0.7 Hz, green). In black are delay times modeled from the forward problem of our tomographic model. You can see the ray theory approximation in our forward problem predicts a smaller width for the anomaly     horizontal dimensions of approximately 4 km by 2 km in the north-south and east-west dimensions respectively (Figure 2.9). Compared to the Beachly et al. [2012] model, our LVZ is of larger magnitude (~7% greater)  and is more consistent with the low-velocity sill model proposed to explain a secondary arrival [Beachly et al., 2012], although our anomaly magnitudes are not large enough to produce the observed secondary arrival. Problems associated with resolving low velocity anomalies are addressed below. 
20  
                 Figure 2.9:  Cross sections through our tomographic model at various depth ranges. Images are contoured with respect to constant fractional change in velocity. White triangles denote the stations from the 2008 experiment and gray lines denote the caldera faults. Resolution is poor outside the caldera.    We also recover low velocities at shallow (< 1 km) depths beneath the caldera. This low velocity region is of slightly larger magnitude (~3%) than the Beachly et al. [2012] tomography, but also only extends to ~1 km depth, compared to the ~2 km depth of the Beachly et al. [2012] anomaly. Stretching from this surface anomaly to the LVZ is 
21  
a low velocity pipe of ~6% magnitude, which is observed in both tomographic images.   In both models, flanking the LVZ is a high velocity (~12%) anomaly extending from near the surface to ~5km depth. This anomaly is predominately observed in the north-east and south-west and compares favorably with the gravity anomalies of Getting and Griscom [1988] and seismic anomalies from previous studies [e.g. Achauer et al., 1988].  To check whether our model dimensions were influencing our inversion results, we tested the effect of varying the vertical extent of the tomographic model and concluded that a model from the surface to 10 km bsl is the smallest model that reasonably fits all the observations. This means, however, that delay times from deeper structures will be mapped into our model. These features will have longer wavelengths [e.g. Nolet and Dahlen, 2000]. Observed travel times (Figure 2.6 and Figure 2.7) do show long wavelength behavior, though this behavior likely has a purely surficial geologic explanation, as discussed later. The depth limitation of our model is primarily influenced by the poor illumination by teleseismic events.  We conducted several resolution tests, which show that velocity recovery is best in regions of high ray density as well as of complimentary ray coverage by both teleseisms and active source data. Although there are 81 densely-spaced stations, they form a roughly linear profile and most of the usable teleseismic events arrive perpendicular to the array. Due to the sparcity of events and imperfect event coverage, there are limited teleseismic ray crossings. The teleseisms themselves, with no active source data, contribute very little resolvable information to the inversion. It is only through the coupling of the active source and teleseismic data in the inversion that we can 
22  
resolve structures (Figure 2.10). The teleseisms, due to their roughly vertical incidence angle, provide horizontal constraints on structure. In contrast, the active source data provides primarily the vertical constraint. It is this combination of horizontal and vertical constraints, provided by the complimentary data sets, that allows for the imaging of the low velocity anomaly at 3 km depth. Neither the active source nor teleseismic tomography is able to accurately resolve this anomaly by itself.   The observed teleseismic waveforms exhibit complexity that modeled waveforms run through the tomographic model do not, indicating that our tomographic image does not fully reconstruct the actual structure. While our data suggest strong 3-D lateral heterogeneity, we make a 2-D simplifying assumption, though we consider the effect of this assumption to be largely negligible in the qualitative assessment of waveform character.  We run a synthetic teleseism through a cross-section of our model. The observed heterogeneity of the waveform on both the vertical and radial channels is noticeably lower than the observations. For the synthetics, the first arrival is always of greatest magnitude, while for our data this is not always true. Because much of the arriving energy is reverberations and reflections, we modify our model to include discrete interfaces to increase later arriving energy.  While models with interfaces produce more late arriving energy on the vertical and radial channels, it is not of the same character as the energy from our data (Figure 2.3). We then include a high Vp/Vs, low Vp anomaly in the location of our LVZ (Figure  
23  
Figure 2.10: Top: Plot of recovered tomographic velocities along line A to A'. Upper Middle: Plot of fractional change in velocity with respect to our starting model. Lower Middle: Plot of the fractional change in velocity for the Beachly et al. [2012] model, relative to our starting model. Bottom: Fractional change in velocity with respect to the recovered Beachly et al. [2012] model.        
24  
 2.11). This better reconstructs the character of the waveform. While entirely qualitative, these simple models suggest that our delay and travel-time measurements are insufficient to reconstruct the anomalies at depth beneath Newberry. This inability to reconstruct anomalies is in part limited by our station geometry and event locations. More fundamentally, our results suffer from waveform healing, where the first arrival time measurements do not fully sample the anomalously slow regions, as elaborated on by Beachly et al. [2012]. Despite the near vertical incidence angle, the minimum travel time still avoids the LVZ for stations in the caldera. Because our arrival times are not directly sampling the magma chamber we cannot accurately reconstruct this feature without using more complicated methods. We can, however, use this ray bending to infer large low velocities in the upper crust [Flecha et al., 2004].   Discussion:  Our inversion reveals large lateral variations in seismic velocity that may reflect changes in temperature, composition, porosity and partial melt. More importantly, our inversion reveals a low velocity body consistent with a magma chamber that compares well with both previous studies and the integrated tectonic history of the volcano and surrounding region. While we recover larger lateral velocity variations than previous tomographic images, the velocity model remains insufficient to predict the observed waveform heterogeneity. This ultimately implies that imaging fine scale structure in volcanic regions must also incorporate aspects of the waveform other than first arrival times. Low velocities at shallow depths (< 1 km) within the caldera are interpreted as 
25  
caldera fill deposits. In comparison with Beachly et al. [2012], we resolve a lower velocity anomaly at a more geologically reasonable, shallower, depth range (< 1 km as opposed to < 2 km) (Figure 2.10). Caldera drilling cores show that there is a predominance of ash flow and tuffs in the upper 500 m with interbedded basaltic flows   
Figure 2.11:  Plot of various synthetic models. First column records the synthetic seismograms as a result of a cross section through our tomographic model (A-A'). Middle row show the vertical traces, bottom row shows the radial (horizontal) traces. Second column shows the result for a more discretized (sharpened) version of our tomographic cross section. Last column denotes a model derived from our tomographic cross section with a 3 km/s (P wave) magma chamber superimposed. To recreate the heterogeneity observed in our data, especially the radial (Figure 2.4), we need to include a substantially lower velocity magma chamber. 
26  
and rhyolites at depth [Keith and Bargar, 1988]. The interbedded flows are largely brecciated, creating a lower than average seismic velocity. We expect the upper 500 m to contribute significantly to our observed late arrivals (low velocities) with partial contribution from deeper brecciated material.    Trending down from the low velocity in the upper kilometer of the caldera to ~3 km is a low velocity pipe-like feature (5% reduction in velocity)  (Figure 2.9 and Figure 2.10). This feature has been interpreted by Achauer et al. [1988] as a breccia dominated flow, which would presumably have high porosity and hence low velocities. In contrast, Macleod and Sherrod [1988] predict that, due to the relatively long-lived magma body at depth and the consistency of eruption, there might be a thermal conduit stretching from the magma body to the caldera surface. Assuming a granitic composition for the thermal conduit, utilizing a change in velocity with respect to temperature of 0.00039 (km/s)/1°C [Christensen and Mooney, 1995], and assuming a reference velocity of 4.25 km/s (our 1-D model at about 2km) we conclude that to achieve a 5% velocity reduction would require a ~500°C anomaly. This seems to be unreasonable, despite high temperatures at shallow depths (265°C at a depth of 930 m [Keith and Bargar, 1988]) and a relatively steep Cascade arc geotherm (~45°C/km) (Rothstein and Manning, [2003]).  We conclude that the dominant controls on seismic velocity in the anomalously slow pipe are porosity and composition and not temperature, in spite of large temperature variations.  We interpret the high velocities beneath the caldera, surrounding the magma chamber as cooled intrusives (Figure 2.9 and Figure 2.10). These features are probably mafic on the flanks, with increasing silicic content with decreasing horizontal distance to the caldera center. The absolute velocities are consistent with a combination of basaltic 
27  
and felsic components. The velocities are too slow to be entirely mafic without appealing to partial melt. Gettings and Griscom [1988] modeled the gravity field and came to a similar conclusion. They note that the structure is not dense enough to be purely solid mafic rock and consider the case of partial melt decreasing the density. However, because the density is different from solid basalt, the amount of melt (~30%) needed to explain this density anomaly would be seismically observable and would decrease the seismic velocities drastically lower than what we observe. We therefore appeal to felsic material to decrease density and velocity and note that small degrees of partial melt could still be present (< 5%).  In addition to the fast velocities beneath the caldera, our model recovers fast velocities on the southwest flank at shallow depths (Figure 2.10). This is important because the Beachly et al. [2012] model placed slow velocities there, but this low velocity also marks the area where ray paths from the 2008 High Lava Plains source go through the magma body and map to the surface. This potentially means that the low velocities to construct a more realistically sized magma chamber were evident in the Beachly et al. [2012] dataset, but the ability and resolution necessary to place the delay times in the correct spot were not.  The most interesting feature is the low velocity zone (LVZ), sitting beneath the center of the caldera and stretching northward ~2 km at about 3-5 km depth, flanked by high velocities in the NE and SW. Assuming a constant composition between the cooled intrusives adjacent to the LVZ and the LVZ itself, we infer a partial melt percentage of around 5%, assuming all variation in velocity is related to melt [Chu et al., 2009]. We interpret the LVZ to be a felsic magma chamber due to the scarcity of basaltic and 
28  
preponderance of rhyolitic volcanism historically inside the caldera during the Holocene. If this magma chamber is entirely molten, we would expect a drastically lower P-wave velocity than is recovered. Synthetic tests show that these low velocities required for melt cannot be tomographically recovered with our array geometry. While, our results do not necessitate a large melt fraction, we infer the presence of a partially molten magma chamber.   Notably, the magma chamber has similar shape to the caldera, except that it is rotated 90 degrees. Elongation of the caldera in the E-W direction and elongation of the magma chamber in the N-S direction are both consistent with an extensional fault running through the caldera with a general N-S orientation. The long axis of the caldera maps to the extension direction [Acocella et al., 2004]. This is related to the partial reactivation of pre-existing extensional faults, creating a larger than expected region of failure. In contrast, the magmatic system of such an environment could feasibly be oriented in the opposite direction. Magma tends to preferentially align in the direction of weakness, in this case perpendicular to the extension direction of the fault. Moreover, the inferred magma chamber's orientation roughly mimics the shape of the volcano itself. The ease of melt movement due to the fault would allow magma to erupt in one plane of direction (along the fault) over another (perpendicular to the fault). Seismically fast areas, interpreted as cooled intrusives are located predominately in the east and west and are roughly perpendicular to the inferred fault direction, consistent with the fault hypothesis.   Our results require only a single magma chamber, though secondary magmatic systems could still reasonably exist in the deep crust. Other geologic evidence is consistent with a single magma chamber.  First, Newberry rhyolitic eruptions during the 
29  
Holocene have erupted similar chemical and compositional aphyrric magmas [MacLeod and Sherrod, 1988]. This is suggestive of a uniform source and would not be expected as a result of many distinct magma plumbing systems erupting rhyolitic magma. We also note that the aphyrric magma, lacking significant phenocrysts, has not crystallized very much during the 2000-3000 year repose time between eruptions. This is suggestive of a long-lived magmatic system, and probably a large magma chamber. Finally, all best fitting models to the secondary arrival observed by Beachly et al. [2012] have a single magma chamber. All these results, while not excluding multiple magma chambers, require only a single magma chamber. Due to our poor teleseismic resolution and limited sampling depth of active source data, we cannot preclude large melt fractions at depths greater than ~7 km.   On the basis of a lack of attenuation at depth, Zucca and Evans [1992] suggested that the inferred magma chamber was solidified, but must also be hot, dry, and fractured as well to account for the observed thermal and seismic velocity anomalies. To test this hypothesis, we use porosity/velocity calculations to show that the porosities needed to achieve our velocities are unlikely at 3-5 km depth. Utilizing the Christensen and Wilkens [1982] basalt velocity/porosity relationships, we find that our velocities correspond to approximately 0% porosity and approximately 7% porosity for our fast velocities and slow velocities respectively. Using the Zamora et al. [1994] velocity/porosity/depth discrete measurements, which correspondingly line up well with the Christensen and Wilkens [1982] porosity/velocity curve, we see that at depths greater than 3 km we expect low porosities (~ 5%). Because the low velocity we recover is most likely under recovered, as evidenced by our synthetic calculations, we conclude that the porosities 
30  
needed to explain our velocity variations, assuming these porosity variations are not melt filled, are implausible though still possible. We also note the porosity velocity relationships derived from the Christensen and Wilkens [1982] study as well as the velocities from the Zamora et al. [1994] study include fluid filled cavities in the porosity calculation. This is opposed to the Zucca and Evans [1992] interpretation of a dry, cracked pluton.  In contrast to the Zucca and Evans [1992] interpretation, Beachly et al. [2012], prefer a pure melt sill as one of three possible models. Because this model seems to correspond favorably with our velocity structure (Figure 2.9 and Figure 2.10), we test its physical plausibility. We perform a simple 1-D calculation from Turcotte and Schubert [2002] to investigate whether or not it is feasible for a magmatic sill to have solidified since the Big Obsidian Flow eruption, 1300 years ago (Figure 2.12). Our values are 320 kJ kg-1 for latent heat of fusion, 1.2 kJ kg-1 K-1 for the specific heat, and thermal diffusivity of 0.5 mm2/s. We vary the temperature difference between the melt and the background medium as well as the thickness of the sill. Assuming the Beachly et al. [2012] waveform interpretation of 600 m sill currently underneath the Newberry Caldera and assuming the lack of new magmatic input, we can see that the inferred sill from the eruption of the Big Obsidian Flow would have been in excess of 800 m thick. For the cooling of the sill, the assumption of bottom up and top down symmetric cooling predicts an equal layer of cooled material above and below the magma chamber. Assuming a current thickness of 600 m and assuming no replenishment, this yields somewhere around 200-300 m of cooled magma on the outside of the magma chamber with about 100-150 m on the top and bottom individually, due to symmetric cooling. This small thickness 
31  
variation is not tomographically observable. We also note, given a current thickness of 600 m, it would take in excess of 1300 years to entirely freeze this sill. We consider the inferred thickness to be plausible and consistent with our and prior results. In spite of the fact that we do not require alternate heat sources to explain current inferred melt measurements, we hypothesize that, due to the consistent repose time and lack of crystallization, basaltic underplating might be keeping the magma chamber warm. Earthquakes recorded in and around the caldera also contribute to our understanding of the magmatic system. The thickness of the seismogenic layer, related to the depth of the brittle/ductile transition is largely temperature dependent [e.g. Ito, 1998 and citations therein], with a maximum temperature of brittle earthquakes of around 250-400°C. Because of this, we can use brittle earthquake depths and locations to infer crustal temperatures. We assume that all earthquakes recorded by the PNSN are caused by brittle ductile behavior and not directly by magmatic activity. We plot earthquakes inside the caldera and map them to our 3-D tomographic model (Figure 2.13). Earthquakes tend to be plotted in the direct vicinity of, or above, our inferred magma chamber. Located in the northern part of the caldera, earthquakes tend to occur on the western edge of Paulina Lake. Predominately in the 1-5 km depth range, the spatial distribution of the earthquakes compares well with our model and the inferred magma chamber. Our 3-D outline of the low velocity zone shows earthquakes occurring close to the inferred magma chamber (Figure 2.13). Large errors in event locations (~0.5 km horizontal, ~0.7 km vertical) make precise interpretations impossible. If earthquakes are not brittle but instead related to magmatic activity, the locations of the earthquakes are evidence for either substantial errors in our model or in the earthquake location. Locations of the earthquakes 
32  
correspond poorly with the location of inferred magma, but correspond well with regions of inferred high temperature gradient in the vicinity of the inferred magma chamber. Dzurisin [2008] noted that repeated levelings showed possible uplift inside the caldera sometime between 1931 and 1994. Using a simplified Mogi model, Dzurisin  [2008] calculated this magmatic input was approximately 0.06 km3 at 10 km depth. In contrast, the total volumetric output from the Big Obsidian flow 1300 years ago was around 0.16 km3 [MacLeod and Sherrod, 1988]. The input of magmatic activity was centered near the southern edge of Paulina Lake. Curiously, this is not where we would expect magmatic activity at depth. There are also a lack of earthquakes recorded in the vicinity of this area.  A spherical volume of 0.06 km3 is substantially smaller than our seismic wavelength, so it might also be that, if this is an isolated flow, we will not be able to resolve it teleseismically nor with our active source data because of the shallow turning depth of the rays. This magmatic input is also presumably basaltic, since the felsic magma chamber is located at ~3-5 km depth. Because one of the mechanisms for rhyolitic eruptions is basaltic underplating, and this magma is coming from under the volcano, it suggests that Newberry will likely continue to erupt in a similar manner to how it has historically throughout the Holocene, if this magmatic upwelling is common. However, we note that compared to the inferred dimensions of our magma chamber, 0.06 km3 is insignificant and unlikely by itself to alter the magmatic system beneath Newberry.  From geothermal borehole data which provides temperature, composition and porosity both inside and outside the caldera, it is apparent that the dominant effect on absolute seismic velocity and relative seismic velocity variations in the upper kilometers is a combination of compositional differences and porosity differences. The relatively 
33  
small derivative of velocity with respect to temperature means that very large variations in temperature are needed to explain small variations in velocity. We therefore consider the shallow thermal structure beneath Newberry largely unresolvable with seismic methods.   
Figure 2.12: Top: Plot of various solidification times for different thickness sills and different temperature differences between the magma chamber and background crust. Bottom: Thickness of sill remaining for after 1300 years (time since Big Obsidian Flow) for various starting thickness sills as a function of difference in temperature between magma chamber and background crust.      
34  
  
Figure 2.13: 3-D plot of our tomographic model, around the caldera. The blue denotes areas of 0.08 or higher fractional velocity incerase. The red denotes areas of -0.07 or lower fractional velocity decrease. Plotted in black are earthquakes in the caldera as recorded by PNSN stations. Topography is raised by 3 km, and contrast is increased by a factor of 3.   Conclusion:  Joint seismic tomography imaged the magmatic system beneath Newberry Volcano. The passive data, recorded on a dense, linear array of 81 stations that were deployed for ~2 weeks, showed variations in structure over short (~600 m) wavelengths. This system as imaged is indicative of low degrees  (< 7 % partial melt); however, 
35  
 synthetic tests reveal that larger degrees of partial melt could exist but are not well constrained by the tomography. Our results are consistent with constraints on magma chamber size from previous geophysical studies and match especially well with a magmatic sill model put forward by Beachly et al. [2012]. Our results indicate that the joint use of teleseisms and active source picks in tomographic inversions of shallow (<10 km depth) structure yield better, more robust results than either dataset individually. Synthetic tests reveal, however, that with the addition of teleseisms, increased care must be taken to ensure that the measured delay times are indicative of travel time structure and not reverberations off of distinct boundaries. Despite these shortcomings, this dataset shows that dense arrays will benefit from long instrument deployment times to allow for better recording of passive events, but also implies that in order to image fine scale structure in volcanic regions, other aspects of the waveform different from the first arrival time must be included.      
36  
CHAPTER III RECEIVER FUNCTIONS  Introduction:  Receiver functions are often used to study the subsurface by looking at P-s or S-p conversions. Higher order conversions and reverberations are often resolved using this technique. The technique is motivated by the idea that the seismic record is a convolution of the Green's function (i.e. response of the earth to an impulse) between two points and a source signature [e.g. Langston, 1979]. Deconvolution, in theory, removes the source component but not the Green's function. Because receiver function analyses work best on seismic wave conversions, they are most often used to map sharp boundaries at depth where discontinuities in the seismic velocity structure exist, in contrast to the smooth velocity images derived from seismic tomography. Moreover, because the receiver function is a relationship between the P wave and S wave, you can use them to estimate Vp/Vs, which provides a measure of partial melt.  In contrast to most studies, we use receiver functions to search for structures at shallow (< 5 km) depths; specifically on Newberry Volcano.   Receiver function studies on volcanos and volcanically active areas have yielded interesting results. In Yellowstone, Chu et al. [2009] use the P-s conversions of a teleseismic wave at a shallow, ~10 km deep magma chamber. By 1-D modeling the receiver function find evidence for a low velocity at depth consistent with a magma chamber.  They also use other characteristics of the seismic waveform to further constrain these low velocities. By looking at the polarity on the horizontals of teleseismic 
37  
waveforms, they are able to show that the wavefront wraps around a presumed magma chamber. This wrapping around only occurs if the magnitude of the anomalous velocity is large. In a different location, Schlue et al. [1995] observed a localized, anomalously large isolated pulse that was interpreted as a P-wave passing through a magmatic zone and then converting into an S-wave at the top boundary of the sill around ~20 km depth.   Leahy et al. [2013], in contrast, mapped structural topography at very shallow depths (~3 km), but in a non-volcanic environment. They were able to resolve such shallow structure, in part, due to high station density, as well as the higher frequency content of their arrivals. Increased station density increases the ability to observe station-to-station coherence in radial receiver function which gives a better idea of the spatial variability of the signal allowing for more robust subsurface characterization. The high frequency content (up to 15Hz) allowed for unique identification of arrivals, leading to a better constraint on subsurface structures and their topography.    Others have used variations in the arrival times on the radial, transverse and vertical channels to infer shallow, low Vs, high Vp/Vs structures at the surface [Owens and Crosson, 1988]. This means that instead of the initial P wave mapping energy to both the radial and vertical channels, the P-wave maps to mostly the vertical channel and a P-s conversion maps to the radial. This ultimately leads to a delayed pulse in the receiver function [Owens and Crosson, 1988].   Here we use a dense array to investigate shallow structure beneath Newberry Volcano. We utilize two methods of receiver functions: the radial/vertical receiver function primarily to map P-s conversions and the vertical component receiver function that can record P-s as well as various other wave conversions and reverberations [e.g. 
38  
Schmandt and Clayton, 2013]. We also investigate transverse/vertical receiver functions. The density of our experiment provides the necessary trace to trace coherence to map the magmatic plumbing at depth.    Methods:   Receiver functions derived from teleseisms recorded by the 2008 seismometer deployment [e.g. Beachly et al., 2012] illuminate the structure beneath Newberry Volcano. We are limited to 21 candidate teleseismic events (Figure 2.3) to use for receiver functions, due to the short, ~2 weeks, deployment of the array. We subset events based on predicted iasp91 [Kennett and Engdahl, 1991] arrival times of various phases, choosing to only use events with P-wave arrivals that appear 5 seconds or more before any other arrival. This ensures that any receiver function will have an isolated response, allowing for the resolution necessary to resolve near surface structure. Of the 21 possible events, we only use 9 for receiver functions (Figure 3.1). We are limited due to insufficient difference between arrival times of phases as well as poor signal to noise ratios on the seismic channels. To compute traditional receiver functions, we deconvolve the vertical from the radial and transverse components, respectively [e.g. Langston, 1979], looking for events that have a clear arrival on both the radial/transverse channel as well as the vertical. We use the frequency domain, water-level deconvolution method [Helmberger and Wiggins, 1971; Clayton and Wiggins, 1976] for calculation of our receiver functions. We use a water-level of 0.01, qualitatively chosen from analyses. We filter the receiver functions with a Gaussian filter of center frequency 0.75 Hz and 0.65 Hz, filtering both on the initial traces as well as post calculation of the receiver function.  
39  
 Figure 3.1: Plot of locations of events used for our receiver functions.            We also use vertical component receiver functions. This method assumes that an average trace recording of the event is representative of the source signature. The mean trace is derived from averaging of traces after alignment by cross correlation [e.g. VanDecar and Crosson, 1990], which can then be deconvolved from each individual trace [Langston and Hammer, 2001; Li and Nabelek, 1991]. We model our array technique after Schmandt and Clayton [2013], using spectral division to deconvolve the source estimate (mean trace) from each individual trace, and a water-level (0.01) to achieve stability. We calculate the mean trace by only averaging the traces that fall on the flanks of the volcano. By not including stations inside the caldera we seek to isolate reverberations off of an inferred magma chamber, because this signal will not be present in our mean trace. Moreover, we time limit the length of the mean trace to about 5 s. This allows for us to look for secondary reverberations that may or may not be consistent across the array [e.g. Yang et al., 2012]. If we do not time limit the traces, we will be  
40  
unable to resolve any laterally consistent structure across the array. For vertical component receiver functions we look for events that have a high signal to noise ratio on the vertical channel and that have a clear identifiable arrival in our frequency band of choice (center frequency ~ 0.75 Hz).   We perform synthetic tests using both vertical and radial receiver functions using E3D, a 2-D synthetic seismogram code [Larsen and Harris, 1993]. We test various simple models and compare the fit with respect to our results. Our E3D teleseism source is a Ricker wavelet.  Results:  We show that the teleseismic waveforms, specifically the coda, preserve information about the upper crustal subsurface and that variations in the waveforms observed on dense arrays can be used to infer shallow structure. Specifically, the tomographic velocity model which fits only the first arrival times is not sufficient to produce the observed heterogeneity in the waveform coda. Our results reveal that in spite of the fact that volcanoes are geologically complex, there are stationary results that preserve structural information.   We assess the character of the 3 seismic channels for all events. Qualitatively, the vertical component for all events has the largest signal to noise ratio. However, the P wave coda has large amplitudes on both the radial and transverse channels (Figure 2.4). Large radial channel amplitudes may be due to a P to s conversion resulting from the drastic velocity reduction near the surface [e.g. Owens and Crosson, 1988]. The extremely large transverse component observed in our raw seismic data, roughly the 
41  
same amplitude as the vertical in many cases, is indicative of complicated structures at depth and is most likely the result of scattering.  While our tomography model recovers the relative variation in arrival times due to seismic heterogeneity well, it does not recover well the heterogeneity nor amplitudes observed in the seismic coda (Figure 2.11). This inability to recover the coda is indicative of larger and sharper gradient velocity contrasts at depth that are not being mapped into our model. We recreate the coda observed in the seismograms by making synthetic seismograms for various 2-D models. Qualitatively, a low Vp and high Vp/Vs ratio body at ~3km depth beneath the caldera fits the characteristics of the observations. The body accurately reproduces the “ringy” character of the teleseismic waveform. This “ringy-ness” is the result of the low Vs velocities of the body which causes S waves to reflect off of the flat top of the body, producing a seismogram that has a dominant frequency related to the total S wave travel time to and from the interface. Different magma chamber geometries produce different ringing characteristics. For example, a magma chamber that does not have a flat top would not necessarily lead to a stable resonance signature. High Vp anomalies tend to have less energy reflected back, leading to a smaller signal post P wave arrival.   We stack receiver functions, derived from our dataset, by azimuth in bins to recover consistent structure (Figure 3.2 and Figure 3.3). We do not bin by ray parameter in large part due to poor statistics. Our event bins separate the azimuths into 4 quadrants, each representing 90° range of azimuth with the initial bin starting at 0°;  we only have two quadrants that produce results with more than a single measurement. We find that in both the radial receiver functions and the transverse functions the signal is not very 
42  
coherent. However, in the vertical component receiver function this signal is more consistent and we see a pulse at ~3 s located beneath the caldera that we cannot easily explain.  We note in our receiver functions that different frequencies record different structures. At low frequencies, we map out structures that are consistent with the Moho (Figure 3.4). However, the depth resolution using these frequencies (<0.5 Hz) is poor. Shallower, smaller structures, such as an inferred magma chamber, cannot be resolved with this frequency content. In contrast, at higher frequencies we increase our ability to resolve these structures, but find that large scale structures deep in the coda of the wave (e.g. Moho) are poorly resolved. This could be due largely to increased high frequency scattering which plays a large part in the dissipation of energy, rendering late arrivals largely insignificant. If this is true, it shows that scattering along the volcano is largely frequency dependent, as has been noted by previous studies at other volcanic and nonvolcanic localities [Chaput et al., submitted; Margerin et al., 2009]. Because there are multiple conversions/reflections from single layers, a common conversion point (CCP) [e.g. Dueker and Sheehan, 1997] stack is not useful. The low frequency content of the teleseisms combined with the shallow  (~3 km) depth of the magma chamber also combine to make CCP stacks a poor choice to resolve the magmatic system.   Due to the difficulty in interpreting our receiver functions using traditional methods, we turn to modeling. Because our synthetics are 2-D yet our data requires large 3-D heterogeneity, we are necessarily partially misconstructing anomalies. However, 
43  
                    Figure 3.2: Above: Plot of radial receiver functions (left) and transverse receiver functions (right) for ~6 events coming from the North West. Traces were filtered with a Gaussian filter of center frequency 0.75 Hz and 0.65 Hz halfwidth. Below: Plot of vertical component receiver function. Note the blue pulse at ~3 s ranging from a distance of ~-3 km to 0 km. Ringing (noted by vertical continuity of pulses) is seen following 6 s. This is a result of a single event exhibiting multiple arrivals in close proximity to the P wave.                   
44  
                       Figure 3.3: Above: Plot of radial (left) and transverse (right) receiver functions for ~3 events from the South East. Traces were filtered with a Gaussian filter of center frequency 0.75 Hz and 0.65 Hz halfwidth. Below: Plot of vertical component receiver function. Note the blue pulse at ~3 s consistent with Figure 3.2.          
45  
                Figure 3.4: Radial receiver function for event 1102 (Figure 2.1) filtered with a Gaussian filter of center frequency 0.3 Hz. The blue pulse across the array at ~5 s is consistent with the Moho from previous studies.     there are differences in the character of the transverse channel data with respect to the radial channel data (Figure 2.4) suggesting they reflect different structure, meaning our 2-D approximations will not be entirely erroneous.    We create 2-D models using various velocities as well as spatially varying Vp/Vs ratios, also including attenuation. By modeling various subsurface geometries we can test the effects of magma chamber velocities, shapes, sizes and depths. By testing a wide variety of magma chamber representations, we can constrain the magmatic system beneath Newberry Volcano.         Our synthetic tests show several interesting results. We compare the vertical receiver functions that were obtained from our observations to calculated receiver functions derived from the E3D synthetic models. We note that because of the 2-D 
46  
assumption, for all models we record stronger than observed conversions off of the magma chamber. The strongest conversions are in a plane parallel to the wave propagation direction, bisecting the magma chamber and for our 2-D models our stations lie along this plane. In reality, the incoming wave direction is never parallel to the plane defining our stations meaning we record smaller conversions. Therefore, when we compare the synthetics to reality, we find that the synthetics overpredict these conversions.   We choose a velocity structure taken from a cross section through our tomographic model and superimpose a low velocity magma chamber over it. The cross section offers several advantages over a 1-D velocity model. High velocity features surrounding the magma chamber focus and defocus energy, changing the overall structure of the coda. Low velocity near surface features also change the initial local arrivals. By accurately modeling these features, we can better resolve the effect of a low velocity body (magma chamber) on our results.   We run an incident plane wave through a cross section of our model, traveling from left to right with a ray parameter of 0.05 s/km. We aligned traces based on the arrival time on the vertical, and computed radial and vertical receiver functions (Figure 3.5). From the synthetics we see that while there is heterogeneity, we do not recover the ringing in the receiver functions derived from our dataset. We do, however, see a very slight delay in the initial arrival on the horizontal component with respect to the vertical, related to a P-s conversion near the surface as a result of a low velocity layer. The lack of ringing in our synthetic traces is reflected in small amplitudes for both the radial and vertical component receiver functions.  
47  
 With the knowledge that our tomographic model is not accurate, we test various magma chambers to study the effect of different heterogeneities. Here we try different models assuming that the depth to the top of the magma chamber is 3 km and that the magma chamber is consistent with large amounts of partial melt (Vp ~3 km/s), using a Ricker wavelet source of 0.5 Hz. We vary only the thickness of the magma chamber, with thicknesses from a thin 0.25 km, to 1.0 km thick (Figure 3.5 and Figure 3.6). Again, we plot the radial and vertical receiver functions. We notice several curious features among the results. First, we see a notch that occurs in the first arrival on the horizontal traces around ~3 km away from the caldera. This notch continues to grow in size with increased thickness of the magma chamber. This is due to a wraparound effect of the wave [e.g. Chu et al., 2009] where it becomes faster for the wave to ignore the low velocity and go around the other side. This creates a change in polarity on the radial but not on the vertical. In this case, due to the shallow depth of the magma chamber and relatively thin thickness, the wave is just starting to wrap around, coming in more or less around vertical. This makes the amplitudes on the radial much smaller than they otherwise would be. Interestingly, none of our observations show this phenomena (e.g. Figure 2.4). The lack of a notch on the horizontal traces in the data recorded at Newberry may be due to a thinner than assumed magma chamber, a deeper location of the magma chamber or an increased velocity of the magma chamber (ie. less partial melt). It may also be due partially to the 2-D approximation to a 3-D environment. For example, the wave could be wrapping around in a transverse direction instead.   
48  
 Figure 3.5: Left: Model with a thin 0.25 km thick magma chamber, about 4 km in width centered inside the caldera at 3 km depth with Vp of 2.3 km/s and Vs of around 0 km/s. Below: Predicted radial and vertical receiver functions from this model. Note the lack of a clear blue pulse at ~3 s on the vertical receiver function.           The synthetic radial receiver functions show several interesting results. The notch observed on the horizontal traces is readily observed in the radial receiver functions. We see energy propagating away from the magma chamber with an apparent velocity of ~4 km/s. and an actual velocity of ~3km/s after accounting for the alignment of the traces by first arrival. In general, the amplitudes of the respective arrivals are larger with a magma chamber than without. The vertical receiver function show a similar trend.  Introduction of a magma chamber produces larger amplitudes beneath the caldera. Curiously, the blue pulse seen in the vertical receiver functions at ~3 s for all of our observed events is absent 
49  
 Figure 3.6: Left: Model with a thick 1 km thick magma chamber, about 4 km in width centered inside the caldera at 3 km depth with Vp of 2.3 km/s and Vs of around 0 km/s. Below: Predicted radial and vertical receiver functions from this model. Note the lack of a clear blue pulse at ~3 s on the vertical receiver function. The anomaly amplitudes have noticeably increased with respect to Figure 3.3.           when using this source frequency (0.5Hz). Because this pulse is seen for events from different azimuths and deltas, the lack of a pulse is not the result our limited 2-D geometry. Testing of various source frequencies shows that a 0.75 Hz source better predicts the blue anomaly at ~3 s (Figure 3.7).   This blue pulse is isolated and centered 2.5 km away from the center of the caldera. Modeling shows that this pulse is a Ps conversion. The pulse originates from a conversion off of the magma chamber surface. Because the arrival is clear in our data   
50  
      Figure 3.7: Plot of radial and vertical receiver functions for a source frequency of 0.75 Hz. The model is the same as Figure 3.6 except the magma chamber thickness is 0.6 km instead of 1 km. Note the blue pulse at ~3 s in the vertical receiver functions at around -2 km.    from multiple different azimuths and deltas, we conclude that these stations mark the rough corner of the magma chamber. Qualitatively, a low velocity body of ~3 km/s with a high Vp/Vs ratio best fits the observations. Higher velocities do not produce as strong a pulse, however lower velocities (~2.3 km/s) can recreate the character and are still physically plausible.   We also grid searched over 1-D models to see if the low velocity body was resolvable. We found that even when using synthetic data, the approximation of fitting a 2-D model with a 1-D one rendered the interpretations largely useless. Because our actual data shows 3-D character, a 1-D model is unable to predict the heterogeneity of our observations.  Direct S-waves were poorly recorded on the array. This may be attributed to poor signal to noise ratio on the horizontal channel as well as poor frequency response of the seismometer at longer periods. The large amplitude post P-wave reverberations on both 
51  
the radial and transverse components imply that horizontal coupling is not the reason for poor S-wave resolution. It is unfortunate that the S-waves are poorly resolved because they represent the best opportunity for characterizing melt in the region, due to their sensitivity to partial melt percentage.  Discussion and Conclusion:  We place constraints on the velocity structure in the context of partial melt percentage, and investigate the implications of these findings. The depth and magnitude of the inferred magma chamber, 3 km and 3 km/s respectively, indicate that the velocity reduction is most likely partial melt of a high percentage. Porosity at these depths is too low to account for a velocity reduction of this magnitude [e.g. Zamora et al., 2004]. Assuming the whole velocity reduction is due to partial melt, we conclude that the velocity reduction is consistent with ~30% partial melt [Chu et al., 2009]. These results are consistent with results put forward by Beachly et al. [2012]. While these results do not place a constraint on the size or thickness of the magma chamber, they do show that the teleseismic waveform can be used to constrain upper crustal structure.  Our results show that the magmatic system beneath Newberry has a larger velocity contrast than is currently tomographically imaged. Because of the poor frequency distribution in our dataset, we are unable to uniquely resolve near surface structure, however, we can match the general character of the waveform by including a low velocity body. Longer deployment time of seismometers will lead to a better earthquake distribution and potentially higher frequency arrivals, increasing the resolution power of teleseismic waves in the upper crust. 
52  
CHAPTER IV AMBIENT NOISE  Introduction:  Ambient noise has become a highly popular tool for exploring the subsurface. The method relies on the assumption of a largely diffuse wavefield with noise from many different directions. By cross correlating noise from two stations (or auto correlating the noise from a single station) you can recover the impulse response between the two stations, which should be symmetric around zero lag due to reciprocity. Differences in the symmetry of the cross correlation function between stations is indicative of asymmetric noise distribution. In theory, this method can be used for both surface and body waves although in practice is has predominately been used for surface waves. This is in part because the dominant microseismic signal, the common signal used in cross correlation, propagates as Rayleigh waves, but also due to the relative differences in energy dissipation as a result of geometric spreading for surface waves and body waves (ie. 1/r vs. 1/r2).   Autocorrelation of the seismic wavefield can be used to resolve upper crustal features. It was shown by Claerbout [1968] that the autocorrelation of a normally incident planar wave in a planar media should yield the reflection response of a source at the station location. This theory was later expanded to 3-D and now includes more sophisticated geometries [Lobkis and Weaver, 2001; Wapenaar, 2004]. This method inherently relies on P-waves and S-waves as sources for seismic noise because of the necessity of noise propagation in the vertical direction. Other studies, such as Roux et al. 
53  
[2005], showed that the cross correlation of ambient noise does in fact produce P-waves as well as Rayleigh waves. Further studies showed recovery of what is believed to be reflections off of the Moho for both P and S waves at Transportable Array (TA) stations in Nevada [Tibuleac et al., 2012]. Pulses were observed by stacking autocorrelations of individual day segments for over a year’s worth of data. The use of many stations and consistency in the arrivals across stations was one of the main determining factors in identifying the arrivals. Draganov et al. [2007], used a dense array of stations to test the recovery of the Greens function. Using only 10 hours of data they were able to recover P-waves and Rayleigh waves at shallower depths. They employed common source gathers and common offset gathers to image the subsurface, which greatly increased the signal to noise ratios.   Here we utilize the properties of the autocorrelation to image the magmatic system beneath Newberry Volcano. We place our results in the context of previous studies and models. We utilize the vertical component to image reverberations in P-waves, indicative of reflections off of a magma chamber.   Methods:  We use data collected from the 2008 seismometer deployment on Newberry Volcano [e.g. Beachly et al., 2012], subsetting the data into day segments. We discard data segments for stations that did not record for the full day. Overall, we have a total of 16 possible days of data for each station, with an average of 14 days of useful data for any given station.  We compute Green's functions by autocorrelating noise in a manner similar to 
54  
Bensen et al. [2007]. We use one bit normalization for our traces, but note that because of the sparcity of earthquakes recorded on the array combined with the low frequency content of the arrivals, we see very little difference between one-bit normalization and no normalization (Figure 4.1). Because we are looking mostly at higher (~5Hz) frequencies than are observed from the earthquakes (~1Hz), one-bit normalization is probably excessive. We spectrally whiten by dividing the amplitude spectrum in the frequency domain by a smoothed version of the amplitude spectrum. This preserves local, but often not global maxima and minima. We iteratively smoothed the amplitude spectrum of the data in two hour segments. This made each individual day chunk more manageable and allowed us to see variations in the autocorrelation throughout the day. The smoothing window is artificially chosen to have length of 1/100 of the total number of samples, which was selected by visual inspection of various smoothing windows. We stacked our final results over all days for each station. We filter using a Gaussian filter focusing on two frequencies. We chose center frequencies of 6 Hz and 1.25 Hz with corresponding gaussian half-widths of 3 Hz and 0.75 Hz respectively, filtering both before one bit normalization as well as after the application of spectral whitening. The choice of center frequencies were based on visual comparison of the Green's function for a given station over various days, looking for a reliable and consistent signal. We also investigate the autocorrelation of the active source experiment coda. While we only have a single coda to autocorrelate, we gain amplitude information from this arrival. This is because the active source is impulsive, exciting a range of frequencies. Because of this broad frequency excitation we do not need to spectrally whiten, which alters amplitude information of arrivals. In the highly scattering  
55  
                     Figure 4.1: Plot of the mean autocorrelation for station NWB35 on Newberry volcano for all days, using two hour time increments for stacking. Plotted in black is the stack for traces that have had one bit normalization applied. In red is the stack for traces with no normalization techniques. In blue is the difference between the one bit normalization and no normalization traces. Amplitudes have been scaled by a factor of 10.    environment, the energy in the coda becomes equipartitioned, losing source signature and direction [e.g. Chaput et al., submitted; Chaput et al., 2013]. So, even though we only have one shot, because of the scattering it is as if noise is coming from all directions. We selected Z (vertical) component coda arriving between 15 s to 45 s post P wave arrival (Figure 4.2), although the results were not sensitive to the time window provided we did not include the initial P wave in the autocorrelation. We filtered the data with a Gaussian filter of center frequency 1.25 Hz and half-width of 0.75 Hz.  Despite evidence for high frequency arrivals, it is the low frequencies that seem to give the most robust and stable 
56  
results. This may be due to attenuation and apparent attenuation in the high frequencies. It is not due to the source because the frequencies used to pick data for the High Lava Plains experiment were filtered in the 2-15 Hz region [Cox et al., 2013]. We use this autocorrelation in a grid search over possible magma chamber models.   Results:  We plot record sections of the stack over all days of the autocorrelations for the two filter frequency ranges. The record sections shows various arrivals. First, using the low frequencies, we see a positive (blue) pulse somewhere in the range of 1-1.5 seconds. This blue pulse roughly extends from -7 km to about 3 km (Figure 4.3), although the bulk of the pulse is in a shorter window ranging from about -3 km to 0 km. This pulse represents what we believe is a reflection off of a magma chamber. In contrast, the higher frequency autocorrelation reveal more precise arrival times, however the spatial width of the blue pulse is smaller, stretching from -1.5 km to about 0.5 km and its polarity is harder to constrain. The isolation of the blue pulse in the high frequency plot is indicative of no other major structures being present (Figure 4.4). There is substantial ringing on the flanks of the volcano to the west in the high frequency autocorrelation, but not in the low frequency autocorrelation. This ringing is likely the result of a monochromatic noise source on the flanks of the volcano. This noise can be seen in frequency plots from various stations (Figure 4.5). Because of the narrow width of the monochromatic noise source and the relatively large smoothing window, this effect is not smoothed out by the spectral whitening, a result of the preservation of local maxima and minima when spectral smoothing. 
57  
 
Figure 4.2: Plot of the active source traces, filtered with Gaussian filter of center frequency 1.25 Hz and half-width 0.75 Hz. Traces are aligned based on arrival time and normalized by respective maximum amplitudes. Black lines at 15 s and 45 s correspond to the window used for the autocorrelation of the shot.      The active source coda autocorrelation compares favorably with the low frequency ambient noise autocorrelation. We see the same similar blue pulse at around 1.3 seconds (Figure 4.6). The exact pattern of arrivals is slightly more difficult to discern in the active source data, potentially because of a low signal to noise ratio due to the use of a single event. Similarity in structure resolved between the two data sets (ambient noise and   
58  
              Figure 4.3: Plotted on left is the stack of the autocorrelation at each station, plotted with respect to the position when mapped to line A-A' (Figure 1.1) on Newberry volcano. Plotted on the right is the stack over all autocorrelations for all stations. Data were filtered with a Gaussian filter of center frequency 1.25 Hz and half-width 0.75 Hz.                      Figure 4.4: Same as Figure 4.2 but using different filter frequencies. Data were filtered with a Gaussian filter of center frequency 6 Hz and half-width 3 Hz.  
59  
                      Figure 4.5: Plot of amplitude relative to frequency (Hz) for station NWB02. A monochromatic noise signal can be seen around 4-5 Hz. Plotted in black are the individual frequency day plots. Plotted in red is the mean of the day plots.     active source coda) suggests that the results are reliable. We focus in on the arrival beneath the caldera at about 1.3 seconds (Figure 4.3 and Figure 4.4). We stack the data from stations that show the arrival (-1.5 km around the caldera to 0.5 km) to get better constraints on the errors of our mean trace and hence the tolerable errors for synthetics modeling.  We derive the error bars from the measure of the standard error, assuming selected traces are in theory identical and that any differences are erroneous. Because the isolation of the arrival at around 1.3 sec as well as the lateral continuity of the arrival over a few kilometers, we can model the velocity field as a smoothly varying 1-D velocity  
60  
 
Figure 4.6: Plot of the autocorrelation for the active source data. Data were filtered with a Gaussian filter with center frequency of 1.25 Hz and half-width 0.75 Hz.     model, which does not produce large reflections, and iterate over various magma   chamber velocities, sizes and depths, which will produce a reflection. Our starting 1-D velocity model is derived from the Beachly et al. [2012] starting model. We then grid search over a range of possible magma chamber models. We constrain the depth of the magma chamber using ambient noise autocorrelations due to arrival time clarity, but 
61  
match the amplitude and phase of the blue pulse derived from the autocorrelation of the shot coda because we don’t need to spectrally whiten this data. To produce the Green's function, we autocorrelate the seismogram recorded at a station as a result of a near-vertically incident teleseism with a delta function source (ANIREC, [Levin and Park, 1997]). This gives us the reflection response of the seismometer to an impulsive source at the same location [Claerbout, 1968]. Because of the inherent frequency discrepancy between the theoretical delta-like Green's function and the observations, we convolve the Green's function with a Ricker wavelet source with dominant frequency equivalent to our center frequency. Our results are not very sensitive to the center frequency of the wavelet. That grid searching yields a best velocity of 4 km/s for the magmatic system, starting at around 3 km depth (Figure 4.7). This is about 1 km/s slower than the surrounding rock velocity in the area. Due to the high scattering potential of volcanoes, we also test whether teleseismic earthquakes satisfy the equipartitioning of the wavefield. We focus on Event 1102, which originated in South America (Figure 4.8) because this event has high signal to noise ratio We window around the coda of the wave, making sure not to include an obvious secondary phases. We filter using a Gaussian filter of center frequency 1.25 Hz and half width 0.75 Hz. The autocorrelation shows general agreement with prior results. However, the blue pulse arrives slightly earlier than for the prior noise and active source autocorrelation results. This could be an indication that the waveform is not entirely equipartitioned. Further study is needed to investigate this phenomenon and test when teleseisms satisfy the equipartitioning requirement for noise.  
62  
Figure 4.7: Results of grid searching to fit our model. Left: in blue is plotted the observed autocorrelation for the selected range of stations. In dashed blue are the respective standard errors associated with the pulse. In black is the modeled autocorrelation with no magma chamber. In red, the modeled autocorrelation with a magma chamber. Right: In black and red respectively are the 1-D model without and with a magma chamber. Note the bottom depth of the magma chamber is unresolvable considering the time ranges we are matching.      Discussion:  We investigated vertical variations in magma chamber location and variations in magnitude of partial melt. We place our results in the context of previous studies of the magma chamber, showing that the magmatic system beneath Newberry is more complicated than previously thought. Our results, combined with previous results suggest 
63  
that there may be a shallow (~3 km deep) magmatic system beneath Newberry in the upper crust.   The 4 km/s velocity at 3 km depth found by autocorrelating the active source is interpreted as being part of the magmatic system. We attribute this variation in velocity with respect to the ambient 1-D model to be a result of partial melt. This yields a partial melt percentage in the rock of about 20% at a depth of around 3km depth [Chu et al., 2009]. Interestingly, this does not correspond well with the magma chamber in our tomographic model, which has a magma chamber that sits at about 4km depth and a             Figure 4.8:  Autocorrelation of Event 1102 (see Figure 3.1). Data were filtered with Gaussian filter with center frequency 1.25 Hz and half-width 0.75 Hz. The blue pulse is slightly earlier than as seen using autocorrelations.   
64  
partial melt percentage of ~10%, although the partial melt percentage is poorly constrained as well as the exact depth of the magma chamber.  In contrast, the Beachly et al. [2012] waveform modeling came up with three reasonable interpretations for magma chamber location and geometry: a sill, a mush zone and a crystal suspension zone. All three models fit the tomography equally well. In contrast, the mush zone is the only model that is entirely consistent with the reflection profiling. The other models have the top of the magma chamber at 4 km depth, which is entirely different from the ambient noise results as well as the autocorrelation of the active source shot results. The mush zone model put forth by Beachly et al. [2012]  starts at about 3 km depth and in general has a low magnitude of velocity decrease, especially compared to the sill model. We consider the mush zone model most likely, because it is fit by both datasets. We note however the autocorrelation results could be the result of a thin sill residing ~1 km above a larger magmatic system, so we cannot entirely rule out the other models. The partial melt percentage needed to explain the inferred velocity contrast resulting from the blue autocorrelation pulse (~20%) is substantially less than 50% needed to prevent crystal lock up, rendering the system largely uneruptable in the short term unless magmatic underplating occurs [Bachmann and Bergantz, 2006]. With a deeper magma chamber, invisible to our noise studies, the partial melt percentage could be above 50% and the melt would therefore be considered eruptable.  We discuss implications of constructive and deconstructive interference in the context of reverberated arrivals, source azimuth, depth of interface and wavelength of the dominant noise frequency. We assume a halfspace model composed of two models with differing velocities. To compare with our results, we treat the case of a faster velocity 
65  
layer on top. To model this amplitude effect, we consider the case of plane waves form the lower layer incident from a variety of angles. We make a simplifying assumption that although the plane wave moves is incident at some angle it reaches all points on the bottom layer at the same time (i.e. we assume a ray parameter of 0 s/km outside of the modeled box). The equations we use are: 
 where λ is wavelength, d is depth of interface, c is some constant and vtop and vbottom are the velocities on the top and bottom layers respectively.  is angle of incidence of a planar wave on the bottom interface. Our equation for A0 is related to the angle of incidence on the surface because we are only concerned with the amplitude as measured on the vertical component. These equations were derived from analagous equations in Harmon et al. [2007]. We find that at 0 interface depth, all signals constructively interfere. As we go deeper, we find in general a decrease in the angles which contribute to constructive interference. The amplitude of recovery is a function of the ratio of the depth of interface to the wavelength of the seismic wave (Figure 4.9). Because we are not varying the depth of the interface, the amplitude is frequency dependent. When we go to lower frequencies, we tend to get larger amplitudes causing the prediction of larger anomalies. In contrast, for synthetics we are modeling these interfaces as frequency independent; the Green's function amplitude is only a function of structure and the amplitude of the signal that we are convolving with the delta-like Green's function. We  
66  
 model two cases for our synthetic amplitudes: the case where the top layer has a faster velocity than the bottom and the case where the top layer is slower than the bottom. These two cases exhibit different results primarily because in the first case, total internal reflection can occur causing a decrease in amplitude of coherent signal.   We note several differences and problems between autocorrelation and cross correlation. The autocorrelation function in the frequency domain is inherently zero phase. This is because in the time domain it is necessarily symmetric. If we spectral whiten too much and then stack autocorrelations, we end up with a white spectrum in the frequency domain and a phase of zero for all frequencies. This is effectively a delta function at the origin. In contrast, the cross correlation is only in theory zero phase. Therefore, when we stack the cross correlations, we preserve some complex component in the frequency domain. This means that even if we whiten everything, we preserve some aspect of amplitudes based on the stacking of the complex frequency components. When converting back to the time domain, the result looks less like a delta function.    Figure 4.9: Plotted are the synthetic amplitudes modeled for a range of plane waves. The blue (red) denotes the predicted amplitude of autocorrelation for a two layer model with top layer velocity of 4 km/s (6 km/s) and bottom layer velocity of 6 km/s (4 km/s). Results are normalized to the case where A0 is a constant and all energy constructively interferes. 
67  
 Conclusion:  We image the top of a magma chamber at ~3 km depth that has a velocity consistent with 30% partial melt. Overall, we consider the result to be a simplistic interpretation of reality. We apply these methods in the hope of recovering further information about the subsurface geology and to better constrain the magmatic system. We have shown that modeling of P wave reverberations is useful for understanding magmatic conditions on volcanoes. More importantly, we show that deployment of seismometers for extended periods of time over quiescent volcanoes is useful for constraining the magmatic system.    
68  
CHAPTER V   CONCLUSION     We used active source data, teleseismic earthquake data and ambient noise data to constrain the magmatic system beneath Newberry Volcano. Joint tomography using active source travel times and teleseismic delay times better recovers a low velocity anomaly beneath the caldera at 4 km consistent with magma chamber models from previous studies. Broad matching of an anomalous signal at ~3 s on the vertical component receiver functions is indicative of a magma chamber, with general reverberations on the radial channel suggestive of a high Vp/Vs anomaly. The ambient noise reveals a reflection consistent with a magma chamber depth of ~3 km. The apparent inconsistency between the ambient noise and tomography is probably related to poor tomographic resolution. We consider a magmatic system beginning at ~3 km depth as opposed to ~4 km depth as most plausible. Overall, we find that the dense deployment of stations across the volcano, even for short periods (~2 weeks) of time is useful for constraining subsurface structure. Our results show that through a variety of different methods, one can constrain the magmatic system beneath volcanoes. Further studies would benefit from longer deployment time of seismometers as well as denser, more three dimensional array deployment.       
69  
REFERENCES CITED    Achauer, U., J.Evans, and D. Stauber (1988), High-resolution seismic tomography of compressioinal wave velocity structure at Newberry Volcano, Oregon Cascade Range, J. Geophy. Res., 93(B9), 10,135-10,147, doi:10.1029/JB093iB09p10135.  Acocella, V., R. Funiciello, E. Marotta, G. Orsi, and S. deVita (2004), The role of extensional structures on experimental calderas and resurgence, J. Volcanol. Geotherm. Res., 129, 199-217, doi: 10.1016/S0377-0273(03)00240-3  Bachmann, O., and G. Bergantz (2006), Gas percolation in upper-crustal silicic mushes as mechanism for upward heat advection and rejuvenation of near-solidus magma bodies, J. Volcanol. Geotherm. Res., 149, 85-102, doi: 10.1016/j.jvolgeores.2005.06.002.  Beachly, M.W., E. E. E. Hooft, D. R. Toomey, and G. P. Waite (2012), Upper crustal structure of Newberry Volcano from P-wave tomography and finite difference waveform modeling, J. Geophys. Res., 117, B10311, doi:10.1029/2012JB009458.  Bensen, G. D., M. H. Ritzwoller, M. P. Barmin, A. L. Levshin, F. Lin, M. P. Moschetti, N.M. Shapiro, and Y.Yang (2007) Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, Geophys. J. Int., 169, 1239-1260, doi:10.1111/j.1365-246X.2007.03374.x.  Biryol, Cemal B., G. Leahy, G. M. Zandt, and S. L. Beck (2013) Imaging the shallow crust with local and regional earthquake tomography, J. Geophys. Res., 118, 1-18, doi:10.1002/jgrb.50115.  Catchings, R., and W. Mooney (1988), Crustal structure of east central Oregon: Relation between Newberry Volcano and regional crustal structure, J. Geophys. Res., 93(B9), 10,081-10,094, doi:10.1029/JB093iB09p10081.  Chaput, J., M. Campillo, R. C. Aster, P. Roux, P.R. Kyle, H. Knox, P Czoski Multiple scattering from icequakes at Erebus Volcano, Antarctica; Implications for imaging at glaciated volcanoes, submitted  Chaput, J. A., D. Zandomeneghi, R. C. Aster, H. Knox, and P. R. Kyle (2012)  Imaging of Erebus volcano using body wave seismic interferometry of Strombolian eruption coda, Geophys. Res. Lett., 39,  L07304, doi: 10.1029/2012GL050956.  Christensen, N., and W. Mooney (1995), Seismic velocity structure and composition of the continental crust: A global view, J. Geophy. Res., 100(B6), 9761-9788, doi:10.1029/95JB00259.  
70  
Christensen, Nikolas, I., R. H. Wilkens (1982) Seismic Properties, Density, and Composition of the Icelandic Crust Near Reydarfjördur, J. Geophy. Res., 87(B8), 6389-6395, doi: 10.1029/JB087iB08p06389.  Chu, R., D. Helmberger, D. Sun, J. Jackson, and L. Zhu (2010), Mushy magma beneath Yellowstone, Geophy. Res. Lett., 37, L01306, doi:10.1029/2009GL041656.  Claerbout, J. F. (1968), Synthesis of a layered medium from its acoustic transmission response, Geophysics, 33, 264–269 Clayton, R. W., and R. A. Wiggins (1976), Source shape estimation and deconvolution of teleseismic body waves, Geophys. J. R. Astron. Soc., 47, 151–177. Cotton, J., and R. Catchings (1989), Data report fro the 1983 US Geological Survey East-central Oregon Seismic Refraction Experiment, U.S. Geol. Surv. Open File Rep., 89-124, 60 pp.  Cox, C., G. R. Keller, and S. H. Harder (2013), A controlled-source seismic and gravity study of the high lava plains (HLP) of Eastern Oregon, Geochem. Geophys. Geosyst., 14, 5208–5226, doi:10.1002/2013GC004870. Dawson, P.B., and D. A. Stauber (1986), Data report for a three-dimensional high-resolution P-velocity structural investigation of Newberry Volcano, Oregon,  using seismic tomography, U.S. Geol. Surv. Open File Rep., 86-352, 90 pp.  Detrick, R. S., P. Buhl, E. Vera., J. Mutter, J. Orcutt, J. Madsen, and T. Brocher, Multi channel seismic imaging of a crustal magma chamber along the East Pacific Rise, Nature, 326, 35-41, 1987. Donnelly-Nolan, J.M., 1988. A magmatic model of Medicine Lake Volcano, California. J. Geophys. Res., 93 (5), 4412–4420. Draganov, Deyan, K. Wapenaar, W. Mulder, J. Singer, and A. Verdel (2007) Retrieval of reflections from seismic background-noise measurements, Geophy. Res. Lett., 34, Lo4306, doi: 10.1029/2006GL028735  Dueker, Kenneth G., and A. Sheehan (1997) Mantle discontinuity structure from midpoint stacks of converted P to S waves across the Yellowstone hotspot track, J. Geophy. Res., 102(B4), 8313-8327  Dzurisin, D. (1999), Results of repeated leveling surveys at Newberry Volcano, Oregon, and near Lassen Peak Volcano, California, Bull. Volcanol., 61(1-2), 93-91, doi:10.1007/s004450050264.  Eppich, Gary R., K. M. Cooper, A. J.R. Kent, and A. Koleszar (2012) Constraints on crystal storage timescales in mixed magmas: Uranium-series disequilibria in plagioclase from Holocene magmas at Mount Hood, Oregon, Earth  Planet. Sci. Lett., 317-318, 319-330, doi: 10.1016/j.epsl.2011.11.019 
71  
Fitterman, D. (1988), Overview of the structure and geothermal potential of Newberry Volcano, Oregon, J. Geophys. Res., 93(B9), 10,059–10,066, doi:10.1029/JB093iB09p10059. Fitterman, D., W. Stanley, and R. Bisdorf (1988), Electrical structure of Newberry volcano, Oregon, J. Geophys. Res., 93(B9), 10,119–10,134, doi:10.1029/JB093iB09p10119. Gettings, M., and A. Griscom (1988), Gravity model studies of Newberry Volcano, Oregon, J. Geophys. Res., 93(B9), 10,109–10,118, doi:10.1029/ JB093iB09p10109. Harmon, Nicholas, Donald Forsyth, and Spahr Webb. "Using ambient seismic noise to determine short-period phase velocities and shallow shear velocities in young oceanic lithosphere." Bulletin of the Seismological Society of America97.6 (2007): 2009-2023. Helmberger, D. V., and R. Wiggins (1971), Upper mantle structure of mid-western United States, J. Geophys. Res., 76(14), 3229–3245. Husen, S., Smith, R.B., (2004) Probabilistic earthquake resolution in three-dimensional velocity models for the Yellowstone National Park region, Wyoming, Bull. Seismol. Soc. Am. ,94 (3), 880–896. Kiyoshi Ito, (1999) Seismogenic layer, reflective lower crust, surface heat flow and large inland earthquakes, Tectonophysics,  306 (3-4), 423-433.  Jensen, R., J. Donnelly Nolan, and D. Mckay (2009), A field guide to Newberry Volcano, Oregon, in Volcanoes to Vineyards: Geologic Field Trips Through the Dynamic Landscape of the Pacific Northwest, Field Guide, vol. 15, pp. 53–79, Geol. Soc. of Am., Boulder, Colo., doi:10.1130/2009.fld015(03). Jordan, B. T., A. L. Grunder, R. A. Duncan, and A. L. Deino (2004), Geo- chronology of age-progressive volcanism of the Oregon High Lava Plains: Implications for the plume interpretation of Yellowstone, J. Geophys. Res., 109, B10202, doi:10.1029/2003JB002776. Keith, T., and K. Bargar (1988), Petrology and hydrothermal mineralogy of US Geological Survey Newberry 2 drill core from Newberry Caldera, Oregon, J. Geophys. Res., 93(B9), 10,174–10,190, doi:10.1029/ JB093iB09p10174. Kennett, B.L.N. & Engdahl, E.R., 1991. Traveltimes for global earthquake locations and phase identification, Geophys. J. Int., 105(2), 429–465. Larsen, S., and D. Harris (1993), Seismic wave propagation through a low-velocity nuclear rubble zone, Lawrence Livermore Natl. Lab., Livermore, Calif., doi:10.2172/10130414. Lawrence, R. D. (1976), Strike-slip faulting terminates the Basin and Range province in Oregon, Geol. Soc. Am. Bull., 87, 846–850.  
72  
Langston, C.A., 1979. Structure under Mount Rainier, Washington, inferred from teleseismic body waves, J. geophys. Res., 84(B9), 4749–4762. Langston, C. A., and J. K. Hammer (2001), The vertical component P-wave receiver function, Bull. Seismol. Soc. Am., 91, 1805–1819, doi:10.1785/0120000225. Leahy, G., R. Saltzer, and J. Schmedes (2012), Imaging the shallow crust with teleseismic receiver functions, Geophys. J. Int., 191(2), 627–636. Lees, J.M., (1992) The magma system of Mount St. Helens: non-linear high resolution P-wave tomography. J. Volcanol. Geotherm. Res., 53 (1–4), 103–116. Lees, J. M. (2007), Seismic tomography of magmatic systems, J. Volcanol. Geotherm. Res., 167(1–4), 37–56, doi:10.1016/j.jvolgeores.2007.06.008. Levin, V. & Park, J., 1997. P-SH conversions in a flat-layered medium with anisotropy of arbitrary orientation, Geophys. J. Int., 131(2), 253–266. Li, X.-Q., and J. L. Nabelek (1999), Deconvolution of teleseismic body waves for enhancing structure beneath a seismic array, Bull. Seismol. Soc. Am., 89, 190–201. Lobkis, O., and R. Weaver (2001), On the emergence of the Green’s function in the correlations of a diffuse field, J. Acoust. Soc. Am., 110, 3011–3017. MacLeod, N., D. Sherrod, L. Chitwood, and R. Jensen (1995), Geologic map of Newberry volcano, Deschutes, Klamath, and Lake Counties, Oregon, U.S. Geol. Surv. Misc. Invest. Ser., Map I-2455, 2 sheets, scale 1:62,500. MacLeod, N., and D. Sherrod (1988), Geologic evidence for a magma chamber beneath Newberry Volcano, Oregon, J. Geophys. Res., 93(B9), 10,067–10,079, doi:10.1029/JB093iB09p10067. Margerin, L., M. Campillo, B. van TIggelen, and R. Hennino (2009), Energy partition of seismic coda waves in layered media: theory and application to Pinyon Flats Observatory, Geophys. J. Int., 177, 571-585. Moser, T. J. (1991), Shortest path calculation of seismic rays, Geophysics, 56(1), 59–67. Nolet, G., and F. A. Dahlen (2000), Wave front healing and the evolution of seismic delay times, J. Geophys. Res., 105(B8), 19,043–19,054, doi:10.1029/2000JB900161. Owens, Thomas J., and R. S.Crosson (1988) Shallow structure effects on broadband teleseismic P waveforms,  Bull. Seismol. Soc. Am., 78(1), 96-108.  Schlue, J. W., R. C. Aster, and R. P. Meyer (1996), A lower crustal extension to a midcrustal magma body in the Rio Grande Rift, New Mexico, J. Geophys. Res., 101(B11), 25283–25291, doi:10.1029/96JB02464.   
73  
Schmandt, B., and R. W. Clayton (2013), Analysis of teleseismic P waves with a 5200-station array in Long Beach, California: Evidence for an abrupt boundary to Inner Borderland rifting, J. Geophys. Res. Solid Earth, 118, 5320–5338, doi:10.1002/jgrb.50370. Sinton, J. M., and R. S. Detrick (1992), Mid-ocean ridge magma chambers, J. Geophys. Res., 97(B1), 197–216, doi:10.1029/91JB02508. Stauber, D. A., S. M. Green, and H. M. Iyer (1988), Three-dimensional P velocity structure of the crust below Newberry Volcano, Oregon, J. Geophys. Res., 93(B9), 10,095–10,107, doi:10.1029/JB093iB09p10095. Ritter, Joachim R. R., and J. R. Evans (1997) Deep structure of Medicine Lake volcano, California, Tectonophysics, 275, 221-241  Rothstein, D.A., and Manning, C.E., 2003, Geothermal gradients in continental magmatic arcs: Constraints from the eastern Peninsular Ranges batholith, Baja California, México, in Johnson, S.E., Paterson, S.R., Fletcher, J.M., Girty, G.H., Kimbrough, D.L., and Martín-Barajas, A., eds., Tectonic evolution of north- western México and the southwestern USA: Boulder, Colorado, Geological Society of America Special Paper 374,  337–354. Roux, P., K. G. Sabra, P. Gerstoft, W. A. Kuperman, and M. C. Fehler (2005), P-waves from cross-correlation of seismic noise, Geophys. Res. Lett., 32, L19303, doi:10.1029/2005GL023803. Tibuleac, Ileana M., and D. von Seggern (2012) Crust-mantle boundary reflectors in Nevada from ambient seismic noise autocorrelations, Geophys. J. Int., 189, 493-500, doi: 10.1111/j.1365-246X.2011.05336.x  Toomey, D. R., S. C. Solomon, and G. Purdy (1994), Tomographic imaging of the shallow crustal structure of the East Pacific Rise at 9 30′N, J. Geophys. Res., 99, 24,135–24,157, doi:10.1029/94JB01942. Turcotte, D. L., and G. Schubert (2002), Geodynamics, 2nd ed., Cambridge Univ. Press, New York. VanDecar, J., and R. Crosson (1990), Determination of teleseismic relative phase arrival times using multi-channel cross-correlation and least squares, Bull. Seismol. Soc. Am., 80(1), 150–159. Waite, Gregory P. and S. C. Moran (2009) Vp Structure of Mount St. Helens, Washington, USA, imaged with local earthquake tomography, J. Volcanol. Geotherm. Res., 182, 113-122, doi: 10.1016/j.jvolgeores.2009.02.009  Wapenaar, K. (2004), Retrieving the elastodynamic Green’s function of an arbitrary inhomogeneous medium by cross correlation, Phys. Rev. Lett., 93(254301), 1–4.   
74  
Weaver, C. S., W. C. Grant, and J. E. Shemeta (1987), Local crustal extension at Mount St. Helens, Washington, J. Geophys. Res., 92(B10), 10170–10178, doi: 10.1029/JB092iB10p10170.  Yang, Z., A. F. Sheehan, W. L. Yeck, K. C. Miller, E. A. Erslev, L. L. Worthington, and S. H. Harder (2012), Imaging basin structure with teleseismic virtual source reflection profiles, Geophys. Res. Lett., 39, L02303, doi:10.1029/2011GL050035. Zamora, M., G. Sartoris, and W. Chelini (1994), Laboratory measurements of ultrasonic wave velocities in rocks from the Campi Flegrei volcanic system and their relation to other field data, J. Geophys. Res., 99(B7), 13553–13561, doi: 10.1029/94JB00121.  Zucca, J. J., and J. R. Evans (1992), Active high-resolution compressional wave attenuation tomography at Newberry Volcano, central Cascade Range, J. Geophys. Res., 97(B7), 11,047–11,055, doi:10.1029/ 92JB00492.