Three-dimensional seismic velocity structure and earthquake relocations at Katmai, Alaska

Three-dimensional seismic velocity structure and earthquake relocations at Katmai, Alaska

Journal of Volcanology and Geothermal Research 276 (2014) 121–131 Contents lists available at ScienceDirect Journal of Volcanology and Geothermal Re...

3MB Sizes 1 Downloads 24 Views

Journal of Volcanology and Geothermal Research 276 (2014) 121–131

Contents lists available at ScienceDirect

Journal of Volcanology and Geothermal Research journal homepage: www.elsevier.com/locate/jvolgeores

Three-dimensional seismic velocity structure and earthquake relocations at Katmai, Alaska Rachel Murphy a,1, Clifford Thurber a,⁎, Stephanie Prejean b, Ninfa Bennington a a b

Department of Geoscience, University of Wisconsin-Madison, Madison, WI 53706, USA U.S. Geological Survey, Volcano Science Center, Anchorage, AK 99508, USA

a r t i c l e

i n f o

Article history: Received 29 September 2013 Accepted 26 February 2014 Available online 18 March 2014 Keywords: Katmai Tomography Volcanism Magma storage

a b s t r a c t We invert arrival time data from local earthquakes occurring between September 2004 and May 2009 to determine the three-dimensional (3D) upper crustal seismic structure in the Katmai volcanic region. Waveforms for the study come from the Alaska Volcano Observatory's permanent network of 20 seismic stations in the area (predominantly single-component, short period instruments) plus a densely spaced temporary array of 11 broadband, 3component stations. The absolute and relative arrival times are used in a double-difference seismic tomography inversion to solve for 3D P- and S-wave velocity models for an area encompassing the main volcanic centers. The relocated hypocenters provide insight into the geometry of seismogenic structures in the area, revealing clustering of events into four distinct zones associated with Martin, Mageik, Trident-Novarupta, and Mount Katmai. The seismic activity extends from about sea level to 2 km depth (all depths referenced to mean sea level) beneath Martin, is concentrated near 2 km depth beneath Mageik, and lies mainly between 2 and 4 km depth below Katmai and Trident-Novarupta. Many new features are apparent within these earthquake clusters. In particular, linear features are visible within all clusters, some associated with swarm activity, including an observation of earthquake migration near Trident in 2008. The final velocity model reveals a possible zone of magma storage beneath Mageik, but there is no clear evidence for magma beneath the Katmai-Novarupta area where the 1912 eruptive activity occurred, suggesting that the storage zone for that eruption may have largely been evacuated, or remnant magma has solidified. © 2014 Elsevier B.V. All rights reserved.

1. Introduction The Katmai volcanic cluster is a tightly spaced chain of eight andesite-to-dacite volcanoes located on the upper Alaska Peninsula (Fig. 1A), approximately 350 km northwest of the Aleutian megathrust and 100 km above the top of the subducting Pacific slab. From southwest to northeast, the volcanoes in the group are Alagogshak volcano, Mount Martin, Mount Mageik, Novarupta volcano, Trident volcano complex (henceforth Trident), Mount Katmai, Mount Griggs, and Snowy Mountain (Fierstein and Hildreth, 2001). All but two of the volcanoes form a relatively straight line trending at an azimuth of 66°; Novarupta and Mount Griggs are located to the north of the main line. While several of these volcanoes have been historically active, the vent-forming eruption of Novarupta in 1912 was the largest volcanic eruption, globally by volume, of the 20th and 21st centuries (Hildreth et al., 2003). Trident was the source of the most recent volcanic activity in the Katmai area, producing lava flows and ash plumes of Volcanic Explosivity Index 3 from 1953 to ⁎ Corresponding author. E-mail address: [email protected] (C. Thurber). 1 Present address: BP, Houston, Texas 77079, USA.

http://dx.doi.org/10.1016/j.jvolgeores.2014.02.022 0377-0273/© 2014 Elsevier B.V. All rights reserved.

1974 (Hildreth et al., 2003). Fumaroles on Martin, Mageik, Novarupta, Trident, Katmai, Griggs, and Snowy (Wood and Kienle, 1990; Hildreth and Fierstein, 2000) along with continued seismicity associated primarily with Martin, Mageik, Trident, and Katmai indicate that the volcanic system at large is still active and could produce eruptions (Ward et al., 1991; Jolly and McNutt, 1999; Jolly et al., 2001; Moran, 2003). The cataclysmic June 6th, 1912 eruption ejected 13 km3 of material (dense rock equivalent) over approximately 60 h (Hildreth et al., 2003) and sent an ash cloud around the world. Most of the material in this eruption was ejected from the Novarupta vent, but the summit of Mount Katmai, over 10 km away, collapsed to produce a caldera with a diameter of 2.5 km and a volume of ~ 5 km3 (Fierstein and Hildreth, 2001). Three magmas partially mingled and erupted together: quartzhypersthene rhyolite, pyroxene dacite, and black pyroxene andesite, with more than half being rhyolite (Wood and Kienle, 1990). There are two different hypotheses currently proposed to explain this observation. One interpretation is that because the three magmas vented concurrently and repeatedly, the compositional gaps between them imply the existence and interaction of separate reservoirs that catastrophically erupted when a large rhyolite dike intruded a shallow dacite–andesite magma chamber (e.g., Eichelberger and Izbekov,

122

R. Murphy et al. / Journal of Volcanology and Geothermal Research 276 (2014) 121–131

Fig. 1. Study area and seismic stations. A) Map of the region showing AVO's permanent seismic network stations and volcanoes in the Katmai area. B) Zoomed map of the study area showing both permanent (purple) and temporary (cyan) stations. Volcanoes are marked by red triangles: MN — Martin, Mk — Mageik, Tr — Trident, No — Novarupta, Ka — Katmai, Sn — Snowy, St — Steller. KP indicates Katmai Pass. The area of the model shown in subsequent figures is indicated by the box.

2000). The other interpretation is that a large compositionally zoned magma chamber beneath Mount Katmai erupted via a propagating sill and dike to Novarupta (e.g., Hildreth and Fierstein, 2000). Over the past century, considerable work has been done investigating the 1912 eruption of Novarupta in an effort to understand the processes that controlled both that eruption and large volcanic eruptions in general. The Katmai area is exceptional not only in the huge volume of the 1912 eruption, but also in the spatial concentration of volcanoes with compositionally different magmas. Studies aimed at understanding the 1912 eruption, as well as studies focused on the recent eruptions of Trident

volcano, raised many fundamental questions regarding the interconnectedness of the Mount Katmai, Novarupta, and Trident magmatic plumbing systems at depth (Hildreth and Fierstein, 2000; Coombs and Gardner, 2001; Hammer et al., 2002). The Katmai volcanic cluster is one of the most seismically active volcanic areas in the Alaska-Aleutian chain (Dixon et al., 2013). Approximately 850 earthquakes occur on average each year, with the vast majority of these being small, high frequency volcano-tectonic (brittlefailure) events (Moran, 2003; Dixon et al., 2010). Dixon et al. (2013) report that the average magnitude of completeness for the Alaska

R. Murphy et al. / Journal of Volcanology and Geothermal Research 276 (2014) 121–131

123

Volcano Observatory (AVO) catalog is 0.5 for the Katmai area for the years 2002–2012. Routine earthquake locations produced by AVO for the Katmai area at the time of our study were calculated using the HYPOELLIPSE algorithm of Lahr (1999) and a one-dimensional (1D) velocity model. The earthquakes primarily occur in three clusters associated with (1) Martin and Mageik, (2) Trident, and (3) Mount Katmai, with generally more diffuse seismicity to the northeast in the area around Snowy volcano (Ward et al., 1991; Jolly and McNutt, 1999; Jolly et al., 2001; Moran, 2003), but including a swarm of events near Snowy during 2002 (Dixon et al., 2003). Most of the events are at depths less than 6 km below sea level (bsl), but occasional deep long period events occur as well (Power et al., 2004). The seismicity is dominated by small magnitude events, the largest of which in the AVO catalog have a ML magnitude under 4.0. A previous seismic tomography study for the Katmai area was carried out by Jolly et al. (2007) using AVO network data. The dataset modeled in their study includes P-wave data from local earthquakes recorded by 18 AVO seismic stations distributed throughout the Katmai region, covering an area ~ 70 by 80 km. They produced a threedimensional (3D) P-wave velocity model of the area that images a prominent low velocity zone (~ 10 km across, reaching to ~ 6 km depth) centered on Katmai Pass, extending from Mageik to Trident volcanoes, along with a moderate low velocity zone stretching from Martin to Katmai. Their results support earlier studies that found anomalously high attenuation (Matumoto and Ward, 1967; Matumoto, 1971), delayed teleseismic arrivals (Ward et al., 1991), and weak high gravity and strong high magnetic anomalies (Goodliffe et al., 1991) in the same area. Furthermore, the depth extent of the velocity anomaly is also consistent with petrologic estimates of magma source depths for the 1974 Trident eruptions of 1.8 to 4.4 km below the summit (~0–2 km bsl) (Coombs et al., 2000; Coombs and Gardner, 2001). In addition to the anomalies found beneath Katmai Pass, several of the above studies identify Trident and a broader region between Martin and Katmai as potential magma storage areas. Moran (2003) used earthquake focal mechanisms to determine that the regional stress field was modified near these volcanoes, indicating that seismic events are occurring in response to local, not regional, processes, which could be interpreted as reflecting magmatic and hydrothermal activity. Here, we seek to define seismic structure and magma storage through the development of Vp and Vs models, as well as precise earthquake relocations. Our models achieve finer spatial resolution than Jolly et al. (2007) due to the inclusion of data from a dense, temporary array of 3-component instruments deployed by the Alaska Volcano Observatory (AVO) and the University of Wisconsin-Madison (UW).

arrival picks for both the permanent and temporary stations into a single database. For the velocity inversion, quality restrictions were imposed on the AVO and temporary array data separately; for time intervals where only the AVO data were available, the events were included that had P arrivals for at least 10 stations, S arrivals for at least 1 station, and at least 12 arrivals in total. For the time interval when the data from both networks are available, we only used the events with a magnitude of 1.0 or above as listed in the AVO catalog. All earthquakes were then screened for quality using the SIMUL2000 code (Thurber and EberhartPhillips, 1999) to eliminate events with calculated location uncertainties exceeding 2 km. This provided a grand total of 1584 events.

2. Dataset

3. Methods

AVO has a permanent network of 20 seismic stations in the Katmai area that have been installed incrementally since 1995 (Fig. 1A, Table 1); 15 of these stations have short-period, single-component sensors; 3 have short-period, 3-component sensors; and 2 have broadband, 3-component sensors (Dixon et al., 2010). A joint AVO and UW project deployed an array of 11 broadband, 3-component seismometers (2 Guralp CMG-3ESP instruments, 9 Guralp CMG-6TD instruments) in July 2008 to supplement the permanent network (Fig. 1B). Because the previous geophysical studies had identified Katmai Pass (located between Mageik and Trident) as a primary area of magma storage, the dense seismic array was focused in this area. For most of the temporary stations, the available data covers through April 2009. For this study, we use waveform time series and arrival times for both the P and S phases for the local earthquakes in the Katmai area that occurred from January 2004 to April 2009. For most of this time period, we only have data from AVO's permanent network. We manually picked the arrivals for the temporary array stations for the events documented in the AVO catalog using dbpick in the Antelope® software package (Boulder Real Time Technologies, 2002); we then integrated all

3.1. Waveform cross-correlation

Table 1 Name, location, year of installation, and sensor type (SP1 — short period 1-component; SP3 — short period 3-component; BB — broadband 3-component) for all stations used in this study. The KA stations are the temporary AVO–UW broadband stations. Station

Latitude (N)

Longitude (E)

Elevation (m)

Installation

Sensor

ACH ANCK CAHL CNTC KABR KABU KAHC KAHG KAIC KAKN KAPH KARR KAWH KBM KCE KCG KEL KJL KVT MGLS KA01 KA02 KA03 KA04 KA05 KA06 KA11 KA12 KA13 KA15 KA16

58.210667 58.198833 58.052500 58.264500 58.131167 58.270417 58.64900 58.494000 58.485000 58.296983 58.596833 58.497833 58.383667 58.275000 58.243333 58.307617 58.440017 58.054000 58.381667 58.134333 58.314300 58.251900 58.260200 58.222600 58.215700 58.211000 58.283400 58.231700 58.221000 58.192900 58.180000

−155.32600 −155.49400 −155.30150 −155.88367 −154.96917 −155.28223 −155.00600 −154.54633 −155.04583 −155.06113 −154.34683 −154.70333 −154.79917 −155.20167 −155.18333 −155.11140 −155.74070 −155.57317 −155.29500 −155.16083 −155.091500 −155.152000 −155.131300 −155.144200 −155.084700 −155.019100 −155.139300 −155.266700 −155.191800 −155.185700 −155.100000

960 869 807 1158 884 1065 1250 923 734 1049 907 610 777 732 777 762 975 792 457 472 810 999 1015 994 935 1003 1098 884 899 926 714

1996 1996 1996 1996 1998 2004 1998 1998 1998 2004 1998 1998 1998 1991 1991 1988 1988 1996 1988 1996 2008 2008 2008 2008 2008 2008 2008 2008 2008 2008 2008

SP3 SP1 SP1 SP1 SP1 BB SP1 SP1 SP1 BB SP3 SP1 SP1 SP1 SP1 SP3 SP1 SP1 SP1 SP1 BB BB BB BB BB BB BB BB BB BB BB

We employed a waveform cross-correlation (CC) technique in order to identify nearby events and determine differential times based on the CC results. To implement this approach, typically a threshold is determined for the CC coefficient. If this coefficient is set too low, the resulting

Table 2 Starting velocity model, modified from Jolly et al. (2007). Node depth (km)

Vp (km/s)

Vp/Vs ratio

0 2 4 6 8 10 12 15

3.0 4.2 5.2 5.8 6.2 6.4 6.5 6.8

1.85 1.83 1.81 1.78 1.75 1.75 1.74 1.74

124

R. Murphy et al. / Journal of Volcanology and Geothermal Research 276 (2014) 121–131

R. Murphy et al. / Journal of Volcanology and Geothermal Research 276 (2014) 121–131

125

data may be of poor quality because unrelated events will be linked and produce unreliable differential times. Conversely, if the coefficient is set too high, very little additional data will be gained (Du et al., 2004). To avoid such limitations of a “hard” threshold, we use the BCSEIS algorithm, which employs a flexible threshold approach. BCSEIS compares waveforms in the third-order spectral (bispectrum, BS) domain in order to reduce contamination by correlated Gaussian or low-skew noise (Du et al., 2004). The BS verification process is accomplished by using thirdorder spectral methods to calculate two additional time delays from both the raw and the band-pass filtered waveforms. These delay values are then used to either select or reject the time delay computed by a traditional CC technique using the filtered waveforms (Du et al., 2004). For this study, we apply the BCSEIS algorithm to the waveforms for events that occurred between 2004 and 2009. 3.1.1. Joint inversion for seismic velocity structure and earthquake locations Velocity models and seismic event locations are determined simultaneously using the double-difference (DD) tomography algorithm tomoDD (Zhang and Thurber, 2003), which is an extension of the DD location method (Waldhauser and Ellsworth, 2000) to solve for locations and velocity structure jointly. The DD tomography method has the ability to sharpen the velocity image near the source region because of the combination of the better accuracy of the differential time data and the concentration of the corresponding model derivatives in the source region (Zhang and Thurber, 2006), as well as the ability to produce precise relative earthquake locations. We apply this method to the AVO network and temporary array data for the quality-controlled events from 2004 to 2009. In total, about 27,000 absolute arrival times (73% P), 125,000 catalog differential times (76% P), and 54,000 CC differential times (68% P) were used in the tomographic inversion. The node spacing in the inversion grid is reasonably uniform across the study area, with nodes spaced 4 to 6 km apart in the ENE–WSW direction and 4 to 5 km in the NNW–SSE direction. In the vertical direction, nodes are positioned at 0, 2, 4, 6, 9, and 12 km depth (relative to mean sea level), although at 0 and 9, and 12 km depth the model is not adequately resolved. The initial P-wave velocity model used varies only in the vertical direction and is based on laterally averaging the model results from a preliminary series of tomographic inversions (Table 2). The initial Vs model is constructed by assuming a constant Vp/Vs ratio of 1.72, again based on an average value for Vp/Vs from preliminary inversions. A hierarchical weighting scheme is used for the inversion to define the velocity structure progressively from larger scale to smaller scale. During initial iterations of the inversion, catalog data are weighted more heavily than the differential times to image coarser-scale structure and obtain first-order improvements in locations. Next, differential catalog and finally CC differential times are increasingly weighted during the later iterations in order to make finer adjustments to the relative relocations and the source-region velocity structure. Iterations alternate between the joint inversions for velocity and hypocenter locations and inversions only for the hypocenter locations, due to the more nonlinear character of earthquake location compared to velocity inversion (Thurber, 1981). The damping value was set to 80 to produce a reasonable condition number, and the horizontal smoothing weights for Vp and Vs were set to 10 and 50, respectively, and were a factor of 2 lower in the vertical direction to allow for an expected greater velocity gradient with depth, on average. These values were selected by assessing the model norm and data residuals associated with different values for damping and smoothing parameters, and examining model restoration and checkerboard tests, for which the synthetic data were generated with Gaussian noise, using a standard deviation of 0.025 s to match the real data. The latter is shown in Fig. 2. The main features in the model

Fig. 3. Residual histograms. Comparison of the initial distribution of residuals in red to the final distribution of residuals in black.

are observed using a relatively broad range of values for these parameters, indicating that the model is robust. We compare the initial and final residual distributions in Fig. 3. The final model yields weighted catalog and CC RMS residual of 0.026 s and 0.025 s, respectively, a reduction of about 90% from the initial values. 4. Tomography and relocation results To provide a basis for our discussion of the tomography models, we briefly summarize some general factors that influence Vp, Vs, and Vp/Vs. Seismic velocities and their ratio depend primarily on lithology, porosity, pore geometry, saturation, and temperature. The presence of magma should be evidenced by low Vp and Vs anomalies and a relatively high Vp/Vs ratio, but gas-filled porosity would generally result in low Vp and Vs anomalies and a relatively low Vp/Vs ratio (Foulger et al., 2003; Lees, 2007). Horizontal slices through the Vp and Vs models are shown in Figs. 4 and 5, and in Fig. 6 we show WSW–ENE cross-sections through the Vp and Vs models along the volcanic axis, along with an interpretive diagram overlain on the Vp model. There are several velocity anomalies present in the 3D model resulting from our inversion (Figs. 4, 5, and 6). The model is not sufficiently resolved in the uppermost layer (Z = 0 km) to interpret, due to predominantly vertically traveling waves near the surface. At a depth of 2 km (Fig. 4A), there is a primary area of low P-wave velocity (Vp) extending approximately from Trident to Mageik. This anomaly matches reasonably well with the location and extent of the shallow part (0 to 4 km depth) of the low Vp anomaly identified by Jolly et al. (2007). At 4 km depth (Fig. 4B), the low Vp anomaly diminishes in amplitude and separates into two anomalies, a lower amplitude anomaly directly beneath Mageik and a larger anomaly below Trident. By 6 km depth (Fig. 4C), the anomaly weakens somewhat beneath Trident and vanishes near Mageik. This is in sharp contrast to the results of Jolly et al. (2007), who show the anomaly to be broadening and slightly intensifying at greater depth (4 to 8 km). We attribute this difference to the dense coverage of the temporary array surrounding Katmai Pass, allowing better localization of these anomalies. Our S-wave velocity results are somewhat more enigmatic. At 2 km depth (Fig. 5A), there are low S-wave velocity (Vs) anomalies between Mageik and Trident that correspond reasonably well to those seen in the Vp model. At 4 km depth (Fig. 5B), we find relatively low Vs anomalies near Mageik, Novarupta, and Katmai, and a smaller anomaly near

Fig. 2. Checkerboard test. Results of checkerboard test using the same data and parameters of our final inversion. Input structure is shown on the left and recovered structure is shown on the right. The left column shows the input model and the right column shows the recovered model. The top three panels are Vp, lower 3 panels are Vs, for the depths of 2, 4, and 6 km, respectively.

126

R. Murphy et al. / Journal of Volcanology and Geothermal Research 276 (2014) 121–131

Fig. 4. P-wave velocity model. Horizontal slices through the 3D P-wave velocity model at the depths of 2, 4, and 6 km. Velocity values in km/s. Inversion nodes are indicated by +. The white contour in the Z = 6 km layer denotes the contour of 5% of the maximum DWS value for Vp, a proxy for model resolution; DWS exceeds that value for almost the entire model region for the other two layers and thus is largely not visible. The DWS and checkerboard test results are compatible.

Trident. By 6 km depth (Fig. 5C), there is just a single resolved low Vs anomaly extending from Trident to Mageik. The low-Vs anomaly at 6 km depth approximately coincides with the location of the inflation pattern seen by Lu and Dzurisin (2014), for this area in the time period 1995 to 2010. They model it as a spheroid or sill centered near or east of Mageik's summit, with lateral dimensions of 10 to 20 km and a minimum depth of 5 km. Although direct determination of Vp/Vs is unwise due to the differing resolutions of the Vp and Vs models (Eberhart-Phillips, 1990), we can consider the general character of the anomalies present in the Vp and Vs models to provide some interpretations. At 2 km depth, the main, central anomaly has both low Vp and low Vs, suggesting a volume of high fracture density and/or magmatic or hydrothermal fluids, presuming that a significant volume of magma is unlikely at such shallow depths. At 4 km depth, Vp is low but Vs is not directly beneath Trident, which would be consistent with a volume of gas-filled porosity.

Both Vp and Vs are somewhat low beneath Mageik, which suggests the possible presence of magma or, equally or perhaps more likely, magmatic gases and/or geothermal fluids. Our earthquake relocations clearly show four distinct clusters of events associated with (1) Martin, (2) Mageik, (3) Trident-Novarupta, and (4) Mount Katmai (Fig. 7). Compared to the catalog, the depth scatter of the earthquake locations is reduced in all clusters. In particular, there appears to be a relatively sharp cut-off depth to the seismicity, which varies from cluster to cluster, that is not evident in the catalog locations. Moran (2003) observed variations in maximum seismicity depth and attributed it to variations in the depth to the transition from the brittle to ductile regime, and we concur with that interpretation. Within all of the clusters, many discrete linear features can also be seen. There are also many small swarms that we examined for potentially important characteristics such as evidence of earthquake migration with time.

R. Murphy et al. / Journal of Volcanology and Geothermal Research 276 (2014) 121–131

127

Fig. 5. S-wave velocity model. Horizontal slices through the 3D S-wave velocity model at depths of 2, 4, and 6 km. Velocity values in km/s. Inversion nodes are indicated by +. The white contour lines enclose the area exceeding 5% of the maximum DWS value for Vs, a proxy for model resolution. The DWS contours and the results from the checkerboard test (Fig. 2) are compatible.

In the Martin–Mageik region, the seismicity that was previously lumped into one cluster (Jolly and McNutt, 1999; Moran, 2003; Jolly et al., 2007) can be seen as certainly two and potentially three separate clusters (Fig. 8A), although we note that our study covers a different time period than these previous studies. The contrast between the two or three clusters that we identify and the five clusters identified by Dixon and Power (2009) is primarily due to differences in labeling rather than differences in epicentral relocations. Dixon and Power (2009) subdivide the Mageik cluster into two parts (which can be seen in Figs. 4B, 5B, 7, and 8) and they define a small group of small earthquakes near Martin as a separate cluster. In contrast to the Dixon and Power (2009) results, which show this series of clusters as dipping strongly to the southwest, our depths are centered at a constant value of about 2 km for all clusters. An interesting lineation is illuminated which links the Mageik cluster with the Martin cluster, but lying slightly southeast of the main concentration of events. An intense swarm in

January 2006 at Martin volcano is not spatially distinct from either the seismicity preceding or following that time period. O'Brien et al. (2012) attribute this swarm to shallow magma intrusion based on analysis of earthquake fault plane solutions; our results do not refute this hypothesis, but we do not see the evidence for tight clustering and/or migration of seismicity beneath Martin that would be indicative of a dike intrusion, as seen for example at Augustine volcano (Sumiejski et al., 2009; Syracuse et al., 2011). The lack of detectable migration, however, does not exclude the possibility of a dike intrusion (Roman and Cashman, 2006; Keir et al., 2009). The resolution of our Vp and Vs models is too coarse to shed light on this issue. The Trident-Novarupta epicenters are shifted north toward Novarupta from the catalog locations (Fig. 7). The earthquake depths cluster around 3 km compared to a diffuse spread between 0 and 5 km in the catalog. A swarm that began in late 2007 and peaked in October 2008 occurred in this area, which bounds the north side of

128

R. Murphy et al. / Journal of Volcanology and Geothermal Research 276 (2014) 121–131

Fig. 6. Velocity model cross-sections and interpretation. West-southwest–east-northeast cross-section through the (A) Vp model and (B) Vs model along with (C) an interpretive diagram superimposed on the Vp model.

the Trident-area velocity anomalies. The October 2008 activity appears to migrate from the northeast to southwest, generally away from Trident's summit and toward Katmai Pass (right-hand panels in Fig. 8B). In the absence of any reported evidence for volcanic tremor or

58.3° N

A’

Katmai Novarupta

5. Discussion

Trident

58.2° N Mageik Martin

A

Depth (km bsl)

Elevation (m asl)

155.4° W

155.2° W

A Martin

Mageik

155° W Trident

Katmai

A

2000 1000 -2 0 2 4 6 8

16

24

other signs of shallow magma intrusion, we assume that the 2008 swarm was driven by magmatic gases and/or geothermal fluids (e.g. Moran, 2003). The Mount Katmai cluster relocates into one primary linear feature striking 15°north-northeast with two smaller linear features nearby, one of which is a swarm of events that occurred in June 2006 (Fig. 8C). All of these features seem to concentrate around 2 km bsl except the 2006 swarm, which is slightly shallower. The largest of the linear features dips slightly to the northeast.

32

40

Distance (km) Fig. 7. Earthquake relocations. Comparison relocation map and cross-section. The light gray dots are the catalog locations and the red dots are our relocations. Note the tighter clustering, presence of linear features, and general deepening trend seen in the relocations. Artifacts in the catalog (e.g., the concentration of hypocenters at −1 km (i.e. 1 km above sea level)) are also removed in our relocations.

Important new results of this study are provided by the earthquake relocations. We find both diffuse zones of seismicity and more sharply defined earthquake clusters. The numerous discrete features that may be discerned in the precise relative relocations shown in Figs. 7 and 8 indicate that there are specific structures and zones along and within which seismicity is concentrated. A number of phenomena may cause seismicity to focus in these areas. Hydrothermal activity, gas release, and magma movement can all cause earthquakes in volcanic regions, and a variety of different pre-existing structures may control where these processes are likely to cause earthquakes (De Gori et al., 2005; Lees, 2007). Considering the relatively persistent seismicity in these clusters, magma movement seems an unlikely mechanism, as no volcanic tremor or other clear indicators of magma migration have been observed. The swarms of earthquakes span relatively short time periods, from several days up to a number of months, and are consistent with short periods of gas release or hydrothermal activity. “Halos” of seismicity at other Alaskan volcanoes, which are found at similar depths to the relocated events in the Katmai area, have been interpreted as overlying magma chambers (Sumiejski et al., 2009; Thurber, 2009). This makes determining the absolute depth of these earthquakes an important goal toward which significant progress has been made here. The results of the inversions for velocity structure do not reveal conclusive evidence identifying potential areas of magma storage in the upper 6 km bsl. We infer that there is no large zone of magma storage currently beneath the Katmai Pass area, in contrast to prior studies. Instead, there seem to be several localized anomalies that were smeared

R. Murphy et al. / Journal of Volcanology and Geothermal Research 276 (2014) 121–131

A

58.2° N

129

A

Mageik

58.18° N

Martin 58.16° N

A −155.4° W

Elevation (m asl)

A

−155.3° W

Martin

Mageik

A

2000 1000

Depth (km bsl)

0 2 4 6 4

8

12

Distance (km) 2008

2006

Time (years)

B

Novarupta

Novarupta

A

Trident

58.22 N

−155.1 W

A

A

−155.2 W

−155.1 W

A

A

Trident

2000

A

Trident

2000 1000

1000

0

Depth (km bsl)

Depth (km bsl)

Trident

58.22 N

A

−155.2 W

Elevation (m asl)

A

58.26 N

58.26 N

2 4

4

0

2

4

4

8

Distance (km)

8

Distance (km)

10/22

10/01

9/08

2009

2007

2005

Time (years)

Date in 2008

Fig. 8. Earthquake clusters. Zoomed-in map and cross-sections of earthquake clusters in selected sub-regions: (A) Martin–Mageik, (B) Novarupta-Trident, and (C) Katmai. Earthquakes (circles) are color-coded by time.

130

R. Murphy et al. / Journal of Volcanology and Geothermal Research 276 (2014) 121–131

C A 58.28° N

58.26° N

A −155.02° W

Katmai caldera

A

2000

(m asl)

Elevation

A

−154.94° W

1000

Acknowledgements

(km bsl)

0

Depth

catalog locations, the relative relocations display an abundance of finer scale concentrations of seismicity. Our results indicate that there are a number of seismic structures responsible for focusing seismicity in certain areas. Because earthquake depths affect the interpretation of subsurface features, the more reliable depth estimates for the events also provide important information regarding where magma is potentially stored. The 3D velocity models display noteworthy anomalies at several locations along the volcanic axis (Fig. 6). While the low velocity anomaly at shallow depths beneath Katmai Pass could be interpreted to support earlier hypotheses of the presence of a magma chamber in this region, we believe these anomalies are due to gas-filled cracks and/or a shallow geothermal system as there are no geological or geophysical indicators of magma storage at such shallow depths. We also do not see anomalies in the velocity structure that support or refute either of the mechanisms proposed to explain the 1912 eruption. Further imaging studies using combined body-wave and ambient noise data and other methods are needed to provide more insight into the deeper subsurface structure of the Katmai region. Earthquake source studies can also contribute to our improved understanding of the seismic and volcanic processes occurring in the Katmai region.

2

4

6 0

4

8

Distance (km)

2008

2006

Time (years) Fig. 8 (continued).

into one large anomaly in the previous studies due to limited resolution. We also conclude that there is no significant magma storage at seismogenic depths beneath Novarupta or Mount Katmai, suggesting that the 1912 eruption may have completely evacuated those zones of magma storage, or that any remaining magma has solidified. Although several low Vp and Vs anomalies are present, the P- and S-wave anomalies do not all occur in coincident areas and thus have several possible non-magmatic explanations, as discussed above. We suggest that the most likely zone of magma storage we image lies beneath Mageik. Ambient noise tomography, receiver function analysis, and potentially teleseismic or joint local-teleseismic tomography could potentially improve our ability to extend and further interpret our tomographic results. 6. Conclusions We perform a joint inversion for earthquake locations and seismic velocity structure to determine precise hypocenter relocations and 3D P-wave and S-wave velocity models for the Katmai region. Our precise earthquake relocations are an important contribution toward understanding the active areas in the Katmai-area volcanic systems. Whereas only relatively large-scale seismic structures may be discerned from the

This material is based upon work supported by the Alaska Volcano Observatory and the National Science Foundation under Grant No. EAR-0910674. Our thanks to John Paskievitch, Lee Powell, Cyrus Read, Summer Ohlendorf, and Helena Buurman for their assistance with field work and to Jim Dixon and Scott Stihler for their work on the AVO earthquake catalog. We thank Seth Moran and two anonymous reviewers for their constructive comments on the initial version of the manuscript. We also acknowledge Jeremy Pesicek and Ellen Syracuse for their assistance in carrying out our analyses. Several figures were generated using the GMT package (Wessel et al., 2013). References Boulder Real Time Technologies, 2002. The Antelope relational database system, Datascope: a tutorial (68 pp.). Coombs, M.L., Gardner, J.E., 2001. Shallow-storage conditions for the rhyolite of the 1912 eruption at Novarupta, Alaska. Geology 29, 775–778. Coombs, M.L., Eichelberger, J.C., Rutherford, M.J., 2000. Magma storage and mixing conditions for the 1953–1974 eruptions of Southwest Trident volcano, Katmai National Park, Alaska. Contrib. Mineral. Petrol. 140, 99–118. De Gori, P., Chiarabba, C., Patane, D., 2005. Qp structure of Mount Etna: constraints for the physics of the plumbing system. J. Geophys. Res. 110, B05303. Dixon, J.P., Power, J.A., 2009. The January 2006 volcanic–tectonic earthquake swarm at Mount Martin, Alaska. In: Haeussler, P.J., Galloway, J.P. (Eds.), Studies by the U.S. Geological Survey in Alaska, 2007. U.S. Geol. Surv. Professional Paper (1760-D). Dixon, J.P., Stihler, S.D., Power, J.A., Tytgat, G., Moran, S.C., Sanchez, J., Estes, S., McNutt, S.R., Paskievitch, J., 2003. Catalog of earthquake hypocenters at Alaskan volcanoes: January 1 through December 31, 2002. U.S. Geological Survey Open-File Report OF 03-0267 (58 pp.). Dixon, J.P., Stihler, S.D., Power, J.A., Searcy, C.K., 2010. Catalog of earthquake hypocenters at Alaskan volcanoes; January 1 through December 31, 2009. U.S. Geol. Surv. Data Series 531 (84 pp.). Dixon, J.P., Stihler, S.D., Power, J.A., Haney, M., Parker, T., Searcy, C.K., Prejean, S., 2013. Catalog of earthquake hypocenters at Alaskan volcanoes; January 1 through December 31, 2012. U.S. Geol. Surv. Data Series 789 (84 pp.). Du, W., Thurber, C.H., Eberhart-Phillips, D., 2004. Earthquake relocation using crosscorrelation time delay estimates verified with the bispectrum method. Bull. Seismol. Soc. Am. 94, 856–866. Eberhart-Phillips, D., 1990. Three-dimensional P and S velocity structure in the Coalinga Region, California. J. Geophys. Res. 95, 2156–2202. Eichelberger, J.E., Izbekov, P.E., 2000. Eruption of andesite triggered by dyke injection: contrasting cases at Karimsky volcano, Kamchatka and Mt. Katmai, Alaska. Philos. Trans. R. Soc. Lond. A 358, 1465–1485. Fierstein, J., Hildreth, W.J., 2001. Preliminary volcano-hazard assessment for the Katmai volcanic cluster, Alaska. U.S. Geol. Surv. Open File Rept. 00-489. Foulger, G.R., Julian, B.R., Pitt, A.M., Hill, D.P., Malin, P.E., Shalev, E., 2003. Threedimensional crustal structure of Long Valley caldera, California, and evidence for the migration of CO2 under Mammoth Mountain. J. Geophys. Res. 108, 2147. http://dx.doi.org/10.1029/2000JB000041. Goodliffe, A.M., Stone, D.B., Kienle, J., Kasameyer, P., 1991. The vent of the 1912 Katmai Eruption: gravity and magnetic measurements. Geophys. Res. Lett. 18, 1521–1524.

R. Murphy et al. / Journal of Volcanology and Geothermal Research 276 (2014) 121–131 Hammer, J.E., Rutherford, M.J., Hildreth, W., 2002. Magma storage prior to the 1912 eruption at Novarupta, Alaska, Contrib. Mineral. Petrol. 144, 144–162. Hildreth, W., Fierstein, J., 2000. Katmai volcanic cluster and the great eruption of 1912. Geol. Soc. Am. Bull. 112, 1549–1620. Hildreth, W., Lanphere, M.A., Fierstein, J., 2003. Geochronology and eruptive history of the Katmai volcanic cluster, Alaska Peninsula, Earth Planet. Sci. Lett. 214, 93–114. Jolly, A.D., McNutt, S.R., 1999. Seismicity at the volcanoes of Katmai National Park, Alaska, July 1995–December 1997. J. Volcanol. Geotherm. Res. 93, 173–190. Jolly, A.D., Stihler, S.D., Power, J.A., Lahr, J.C., Paskievitvh, J., Tytgat, G., Estes, S., Lockhart, A.B., Moran, S.C., McNutt, S.R., Hammond, W.R., 2001. Catalog of earthquake hypocenters at Alaskan Volcanoes: January 1, 1994 through December 31, 1999. U.S. Geological Survey Open-File Report 01-189 (24 pp.). Jolly, A.D., Moran, S.C., McNutt, S.R., Stone, D.B., 2007. Three-dimensional P-wave velocity structure derived from local earthquakes at the Katmai group of volcanoes, Alaska. J. Volcanol. Geotherm. Res. 159, 326–342. Keir, D., Hamling, I.J., Ayele, A., Calais, E., Ebinger, C., Wright, T.J., Jacques, E., Mohamed, K., Hammond, J.O.S., Belachew, M., Baker, E., Rowland, J.V., Lewi, E., Bennati, L., 2009. Evidence for focused magmatic accretion at segment centers from lateral dike injections captured beneath the Red Sea rift in Afar. Geology 37, 59–62. Lahr, J.C., 1999. HYPOELLIPSE: a computer program for determining local earthquake hypocenter parameters, magnitude, and first motion pattern. U.S. Geol. Surv. Open File Rept. 99-23. Lees, J.M., 2007. Seismic tomography of magmatic systems. J. Volcanol. Geotherm. Res. 167, 37–56. Lu, Z., Dzurisin, D., 2014. InSAR imaging of Aleutian volcanoes: monitoring a volcanic arc from space. Geophysical Sciences. Springer Praxis Books. ISBN: 978-3-642-00347-9, p. 388. Matumoto, T., 1971. Seismic body waves observed in the vicinity of Mount Katmai, Alaska, and evidence for the existence of molten chambers. Geol. Soc. Am. Bull. 82, 2905–2920. Matumoto, T., Ward, P.L., 1967. Microearthquake study of Mount Katmai and vicinity, Alaska. J. Geophys. Res. 72, 2557–2568. Moran, S.C., 2003. Multiple seismogenic processes for high-frequency earthquakes at Katmai National Park, Alaska: evidence from stress tensor inversions of fault-plane solutions. Bull. Seismol. Soc. Am. 93, 94–108.

131

O'Brien, J.F., Roman, D.C., Dixon, J.P., Power, J.A., Arnold, R., 2012. Multiple causes for noneruptive seismic swarms at Mt. Martin, Katmai Volcanic Cluster, Alaska (2004–2008). J. Volcanol. Geotherm. Res. 229–230, 13–22. Power, J.A., Stihler, S.D., White, R.A., Moran, S.C., 2004. Observations of deep long-period (DLP) seismic events beneath Aleutian arc volcanoes; 1989–2002. J. Volcanol. Geotherm. Res. 138, 243–266. Roman, D.C., Cashman, K.V., 2006. The origin of volcano-tectonic earthquake swarms. Geology 34, 457–460. Sumiejski, L., Thurber, C.H., DeShon, H.R., 2009. Relocation of earthquake families associated with the 2006 eruption of Augustine volcano using the Equal Differential Time method. Geophys. J. Int. 176, 1017–1022. http://dx.doi.org/10.1111/j.1365-246X. 2008.04037.x. Syracuse, E., Thurber, C., Power, J., 2011. The Augustine magmatic system as revealed by seismic tomography and relocated earthquake hypocenters from 1994 through 2006. J. Geophys. Res. 116, B09306. http://dx.doi.org/10.1029/2010JB008129. Thurber, C.H., 1981. Earth structure and earthquake locations in the Coyote Lake area, central California. (Ph.D. Thesis) Massachusetts Institute of Technology, Cambridge, MA (331 pp.). Thurber, C.H., 2009. Alaskan volcano seismicity and structure from sparse datasets. Seism. Soc. Am. Annual Meeting. Thurber, C., Eberhart-Phillips, D., 1999. Local earthquake tomography with flexible gridding. Comp. Geosci. 25, 809–818. Waldhauser, F., Ellsworth, W.L., 2000. A double-difference earthquake location algorithm: method and application to the Northern Hayward Fault, California. Bull. Seismol. Soc. Am. 90, 1353–1368. Ward, P.A., Pitt, A.M., Endo, E., 1991. Seismic evidence for magma in the vicinity of Mt. Katmai, Alaska. Geophys. Res. Lett. 18, 1537–1540. Wessel, P., Smith, W.H.F., Scharroo, R., Luis, J.F., Wobbe, F., 2013. Generic Mapping Tools: improved version released. EOS Trans. Am. Geophys. Union 94, 409–410. Wood, C.A., Kienle, J., 1990. Volcanoes of North America: United States and Canada. Cambridge University Press. Zhang, H., Thurber, C.H., 2003. Double-difference tomography: the method and its application to the Hayward fault, California. Bull. Seismol. Soc. Am. 93, 1875–1889. Zhang, H., Thurber, C., 2006. Development and applications of double-difference seismic tomography. Pure Appl. Geophys. 163, 373–403.