Subduction geometry of the Yakutat terrane, southeastern Alaska

The active tectonic processes responsible for the spectacular topography of Alaska are a product of the relative motion between the Pacific plate and North America. The Pacific plate moves in a north-northwest direction at a rate of ∼50 mm/yr relative to North America (Plattner et al., 2007; Elliott et al., 2010) in southeastern Alaska. That convergent motion creates the Aleutian subduction zone in western Alaska (Page et al., 1989). West of Prince William Sound the subduction system is unambiguous, as evidenced by a well-defined trench, an active Wadati-Benioff zone of earthquake generation, and the Aleutian Islands volcanic arc (Murray, 1945; Van Wormer et al., 1974; Kay, 1978). East of Prince William Sound where the Yakutat Block microplate is located, the situation is more ambiguous (Fig. 1). At this point the Aleutian Trench terminates at the edge of the Yakutat Block and Benioff zone seismicity rates drop dramatically (Page et al., 1989). This geometry has led to questions regarding the nature of the subduction processes operating in southeast Alaska. Competing views of the tectonics in this corner of Alaska include a model of gradual transition from subduction to continental collision (Doser and Lomas, 2000) as well as tearing and underplating slab material (Fuis et al., 2008). A major complication in the past has been that southeast Alaska represents one of the most remote regions of North America. Until the St. Elias Erosion/Tectonics Project (STEEP; http://www.ig.utexas.edu/steep/), the seismic array in this area consisted only of a sparse collection of relatively poor quality short period sensors. We incorporate new data from STEEP to provide some new seismic imaging results to constrain the geometry of this system.

A key component of this system is the Yakutat Block microplate (Fig. 1). The Yakutat Block is thought to be composed of crust that was rifted from the western margin of North America at a latitude between that of present-day Washington or British Columbia (Davis and Plafker, 1986) and northern California (Bruns, 1983) ∼55 m.y. ago. The rifting process that formed the Yakutat Block produced a terrane of mixed crustal composition; the western portion of the terrane appears to be oceanic in composition (Bruns, 1983; Plafker, 1987). The 24–27 km thickness of the Yakutat Block is indicative of an oceanic plateau structure (Christeson et al., 2010). Previously, some suggested that the eastern portion of the Yakutat Block had the composition and thickness typical of continental crust (Bruns, 1983; Plafker, 1987). Refraction data from STEEP (Christeson et al., 2010; Worthington et al., 2012), however, showed that the crust of the Yakutat Block south of the Pamplona zone (Fig. 1) is nearly constant, but the upper crust and lower crust are very different. The upper crust of the Yakutat Block is now known to be a wedge of sediments with a thickness of ∼15 km at the western edge thinning to near zero at its southeastern edge. The lower crust is a complementary wedge that has velocities consistent with an oceanic plateau (Christeson et al., 2010). The wedge is complementary because the Yakutat Block crust seems to be a nearly constant thickness south of the deformation front called the Pamplona zone (Worthington et al., 2012).

If we assume the Yakutat Block originated along western North America, it must have been transported north along the Queen Charlotte–Fairweather fault system (Bruns 1983). The Queen Charlotte–Fairweather fault system forms the boundary between western North America and the eastern portion of the Pacific plate. The collision of the Yakutat Block with North America initiated ∼10 m.y. ago (Davis and Plafker, 1986; Bruns, 1983) and Yakutat collision continues to this day. Evidence of continued deformation is documented by many indicators, including numerous earthquakes (Doser et al., 1997; Estabrook et al., 1992; Boatwright, 1980; Stephens et al., 1980), deformation as measured by global positioning system (Sauber et al., 1997; Fletcher and Freymueller, 1999, 2003, Elliott et al., 2010), thermochronology (Berger et al., 2008; O’Sullivan and Currie, 1996), volcanism (Richter et al., 1990), and the recent uplift history (Plafker and Thatcher, 2008; Chapman et al., 2007).

Although evidence of deformation surrounding the Yakutat Block is prodigious, the subsurface structure of the Yakutat Block and its relation to the mantle was relatively unknown before this study. Three active source seismic profiles have been collected in the region of southeastern Alaska that are important to this discussion. First, the Trans-Alaska Crustal Transect (TACT) project followed the Richardson Highway north out of Valdez with cross lines along the Chitna River valley, as shown in Figure 2 (e.g., Page et al., 1986). The STEEP project also included an extensive offshore reflection and refraction active seismic component. We include important constraints from STEEP02 by Christeson et al. (2010) and STEEP01 by Worthington et al. (2010, 2012) to develop the model for the geometry of the subduction zone described herein.

In addition to these active source experiments, two passive seismic experiments have been conducted in southern Alaska. The Broadband Experiment Across the Alaska Range (BEAAR; http://research.flyrok.org/beaar.html) employed a passive array similar to that used by STEEP and provided additional constraints for the geometry of this subduction system (Ferris et al., 2003; Ai et al., 2005; Veenstra et al., 2006). The BEAAR experiment imaged subducted lithosphere beneath central Alaska (Ferris et al., 2003). We include a downdip segment of subducted lithosphere imaged by the BEAAR experiment to validate the quality of our data in the region north of Anchorage. Our results present the first passive, wavefield imaging methods applied to lithosphere associated with the subducted Yakutat Block.

A second passive imaging result that is an important foundation for this paper is the P-wave and S-wave tomography models developed from local earthquake data described in Eberhart-Phillips et al. (2006); their tomography models and results from BEAAR provide important cross-validation of our results. It is important to recognize, however, that neither of these previous efforts provides significant constraints on upper mantle structure in the STEEP area. BEAAR was deployed hundreds of kilometers from our study area. The Eberhart-Phillips et al. (2006) models are derived from local earthquake data that are composed of ray paths that mostly pass through the crust. Their model provides limited information regarding mantle structure in our study area due to the lack of significant mantle seismicity and sparse station coverage in the STEEP study area. This paper is important because we provide the best images yet obtained of the lower crust and upper mantle structure in southeast Alaska.

This paper presents P-wave and S-wave receiver functions computed from data collected via the STEEP passive seismic network and the Alaska Earthquake Information Center. We calculated both P- and S-wave receiver functions using Wiener filters stabilized by prewhitening. The P-wave migration methods we employed are a variant of those found in Pavlis (2011b). Here we present only a summary of key points about this particular methodology.

The seismic waveform data for this study came from stations installed as part of the passive seismic component of STEEP and additional stations affiliated with the Alaska Earthquake Information Center (AEIC; http://www.aeic.alaska.edu/). The STEEP network includes 22 broadband seismometers in the Chugach–Saint Elias region of south-central Alaska. An additional 11 AEIC broadband stations were included in the P-wave data as well as 33 AEIC stations in the regional S-wave data. The stations that contributed to the P-wave data formed an array with dimensions ∼525 km by 250 km. Station topography varied between sea level and 2.1 km.

The teleseismic waveform data for events with body wave magnitude >5.5 occurring between 2005 and 2007 were extracted from the original data and resampled at 20 samples s^–1. P-wave results used events extracted between 30° and 85°. S-wave results utilized S- and SKS-wave distance ranges of 55°–125°, and SKSdf (deep focus) events between 160° and 180°. The epicenter locations of events used for the P-wave data are shown in Figure 3A and epicenter locations of events used for S-wave results are shown in Figure 3B. These criteria resulted in 1453 P-phase waveforms and 3396 S-phase waveforms. All extracted waveforms were reviewed interactively in order to determine quality before further processing. Events with an estimated signal to noise ratio of 4:1 or greater were passed on for further receiver function calculation processing.

All stations included in this study incorporated broadband seismometers with a minimum response period of 30 s and matched response on all three components. Converting data from velocity to displacement required integration of the gain corrected data for each channel and band pass filtering. The band pass filter parameters included corners of 15 and 0.5 s for P-wave phases and 30 s and 3 s corners for S-wave phases.

The P-wave and S-wave receiver functions were calculated in a similar manner with some minor differences. All waveforms were rotated from the ZNE geographical coordinate system into the LQT ray coordinate system. The LQT coordinate system is oriented such that L is in the direction of the upward moving wave, the Q component is perpendicular to the L component in the vertical direction, and the T component is perpendicular to both the L and Q directions. This rotation was accomplished by interactively selecting the angle of incidence that minimized the amplitude of the deconvolved T component at the time of the P-wave arrival. For the S waves, the deconvolved L component was minimized at the time of the S-wave arrival (Kumar et al., 2005). This procedure was performed in order to more closely match the measured LQT coordinates rather than the calculated LQT coordinates. This difference in orientation was significant for a small number of waveforms. Receiver functions were calculated using an inverse wavelet computed via a Wiener (spiking) filter from a fixed time window of 10 s before the direct arrival and 90 s after the direct arrival (Kind et al., 1995; Yuan et al., 1997). The inverse wavelet was stabilized by employing a prewhitened autocorrelation matrix (Yuan et al., 1997; Treitel and Wang, 1976). The Wiener filter was calculated from the L component for the P-wave data and from the Q component for the S-wave data. In order to minimize noise, the same 15–0.5 s P-wave filter and 3–30 s S-wave band pass filters that were applied before deconvolution were applied again after the deconvolution calculation.

Several additional steps were required to make the S-wave receiver functions comparable to the P-wave receiver functions. Since the polarity of P-wave conversions is the opposite of S-wave conversion for the same velocity contrast interface (Sodoudi et al., 2006), the polarity of the S-wave receiver functions was inverted. In addition the converted P-waves generated from an S-wave conversion travel faster than the S-wave that generated them. As a result, converted waves are recorded before the direct S-wave arrival with the lead time scaled to depth of the conversion. S-wave data are therefore reversed in time to be projected in a manner similar to P-wave data.

The S-wave receiver function data were imaged using a variant of the common conversion point (CCP) stack method that has been developed over a number of years (Dueker and Sheehan, 1997; Yuan et al., 1997; Jones and Phinney, 1998; Gilbert et al., 2003; Wilson et al., 2003; Kumar et al., 2005). This direct imaging method divides the study area into bins. The method we used was built with a focus on imaging of a feature commonly interpreted as the lithosphere-asthenosphere boundary (LAB) often seen with this type of data at a depth of ∼100 km (e.g., Kumar et al., 2005). In this method each receiver function is placed in spatial bin at a position computed from the ray piercing point for an S-wave to P-wave conversion at a depth of 100 km. A time-to-depth mapping operator is computed to map each sample of each S-wave receiver function into a vertical line at the piercing point. We then apply a Gaussian smoother to simultaneously stack and interpolate an irregular distribution of waveform data into a regular grid in map view. The end result is a three-dimensional (3-D) volume comparable to a common mid-point stack of 3-D seismic reflection data. All the S-receiver function volumes presented here used the Preliminary Reference Earth model (Dziewonski and Anderson, 1981) to define the time-to-depth mapping operator.

There is a shortcoming of this method that arises in situations where station coverage is sparse and irregular such as in southeast Alaska. Wide station spacing requires large smoother widths to retain stability. Large smoothing areas mix ray paths that have passed through differing portions of the subsurface that can degrade the stack.

The most important example in this paper is the dipping horizon problem. The CCP method cannot image even low-angle horizons as a single, continuous feature when the stations are widely spaced compared to the vertical resolution length. Under these circumstances dipping horizons imaged by such large bins will be broken into a series of discrete horizontal segments. This segmentation turns a dipping horizon into a stair-step shape in the images we construct.

The P-wave receiver functions were imaged using a unique form of wavefield imaging employing plane wave decomposition known as plane wave migration. We did this via a suite of three computer programs (described by Pavlis, 2011a) collectively known as plane wave migration. The plane wave method employs a regularization scheme to create uniform station geometry by creating a regular Earth surface grid of virtual stations. The node locations of this grid system are referred to as pseudostations. The surface grid employed in this study was centered at 60.5N, 142.8W, with 460 × 200 km dimensions, and 10 km nominal node spacing. The long axis of the surface grid is oriented at 110° azimuth. The nodes of the surface grid are populated with receiver functions interpolated from the receiver functions of real stations via the pseudostation method (Neal and Pavlis, 2001; Poppeliers and Pavlis, 2003a) by the pwstack program.

An important feature of the pwstack program is the ability to group stations based on spatial relationships to geological features such as major faults and terrane boundaries. Processing stations geographically preserves geologic elements near boundaries that would be smoothed if all stations were lumped together. We refer to this form of processing as the split-station method. Our approach was to segregate stations located on the Yakutat Block from stations located on other terranes before the pwstack program. This processing flow produced better image quality at the northern boundary of the Yakutat Block than other processing schemes (Bauer, 2014).

The ability of the plane wave migration method to image dipping layers is achieved by adjusting the lag of each receiver function by the pwstack program. Pwstack calculates a moveout correction for each pseudostation trace relative to the zero lag trace. This lag calculation adjusts each trace in a manner that allows for constructive summation of traces produced by dipping layers. Therefore, pwstack produces a suite of stacked receiver functions that approximate plane waves for each slowness increment of each event.

The pwmig program combines plane wave components to construct 3-D back-projected ray paths. These are used to form a raygrid (Pavlis, 2011a) to locate each sample of each plane wave component at the theoretically correct image point for each sample. The pwmig program takes the data projected into the raygrid geometry to yield an image as a weighted sum of plane wave components. The weights in the summation are defined by the inverse generalized radon transform (Beylkin and Burridge, 1990), which uses the ray path geometry to define this analytic inverse using the geometry of isochrones that are tangential to the ray path at the location of each image point.

Plane wave components in pwmig are propagated to depth using ray-based propagators, i.e., we compute line integrals along ray paths to compute lag times as a function of position along each backward-projected ray path. We modified the 1-D IASPI-91 (International Association of Seismology and the Physics of the Earth’s Interior) Earth velocity model (Kennett and Engdahl, 1991) using seismic velocities in the upper 30 km of portions of the model over the Yakutat region based on active seismic studies (Christeson et al., 2010). This velocity model modification has the effect of correcting the thickness and velocity of the crust portion of the IASPI-91 to conform to the thickness and velocity within the Yakutat crust as constrained by the active source data of Christeson et al. (2010). The velocity model for the region covered by the Yakutat Block was modified by dividing it into two distinct velocity zones, the sediment cover and the oceanic plateau basement directly below the sediment. All velocity model cells that were defined as sediment or basement were given a single velocity value that reflected average velocities obtained from the active seismic portion of the STEEP project (Christeson et al., 2010; Worthington et al., 2012). The maximum sediment thickness found on the Yakutat Block was 15 km, a value found near Kayak Island. Sediment cover thinned linearly to both the east and south from Kayak Island. Cells below the bottom sediment boundary were designated as basement to a depth of 30 km. Below depths of 30 km, standard IASPI-91 velocities were used.

The pwmig program is a prestack migration method. Migrated events are then stacked to produce a single image volume of regular dimensions using the gridstacker program. We assigned equal weight to each event in the final stack.

We were able to use the S-wave data to successfully image a larger portion of southeastern Alaska than the P-wave data, for two main reasons. First, the highest frequency we can obtain with the deconvolved S-wave data is approximately a factor of 3 lower than the P-wave data. When mapping to depth this translates to a spatial scale three times larger than the P-wave data. This makes the very coarse data sampling away from the STEEP region less extreme and allows us to draw some inferences from sparse measurements. The second factor that helps the situation is a consequence of the sign difference in refraction of P-wave to S-wave conversions compared to S-wave to P-wave conversions. Snell’s Law dictates that for a P-wave to S-wave conversion from a flat horizon the converted S-wave ray path is steeper than the incident P-wave. This is reversed for an S-wave to P-wave conversion, leading to S-wave conversion points that are at a much greater lateral distance from a given seismic station than the corresponding P-wave conversion points. As a result, common conversion points for S-wave receiver functions sample Earth structure a greater distance from a seismic station than a corresponding P-wave receiver function. For example, the sparse station spacing between the Cook Inlet region and the STEEP area produced unresolvable holes in the P-wave data. We were able to partially overcome this limitation with S-wave data due to the wider lateral conversion points in conjunction with a relatively large smoother width. Nonetheless, this sparse sampling still yields artifacts that limit our interpretations over much of the study area, and the station geometry must be considered when evaluating the results.

Because of these issues the P-wave imaging volume was limited in both scope and depth. The P-wave image volume extends only between the eastern margin of Prince William Sound and Yakutat Bay, as shown in Figure 2. The depth of the P-wave data migration was limited to 200 km due to the limited array aperture of the densely covered STEEP study area. The S-wave data were processed in two different modes. The first is a regional scale date that encompasses all of southeastern Alaska (Fig. 1). A subset of the regional data was extracted along the boundaries of the P-wave data using a slightly smaller aperture Gaussian smoother. The Gaussian smoother applied to the regional S-wave data employed a 1σ distance of 30 km and a cutoff distance of 60 km. This smaller scale image volume is referred to as the local STEEP S-wave image (Fig. 2). The denser station coverage in the STEEP region allowed us to reduce the size smoother width values to 25 km (1σ) and a cutoff of 50 km.

The P-wave images used the following parameters discussed in Pavlis (2011a); 7 × 7 slowness grid with 0.01 s/km increments, 50 km Gaussian surface aperture, and stations split between the Yakutat Block and terranes to the north. (See Bauer, 2014, for details on the data processing.) A key point is that the slowness parameters in the pwstack program allow for the imaging of features with a maximum dip of ∼20°. Splitting stations at the pwstack stage and then later combining the data after migration allowed the preservation of features across the terrane boundary between the Yakutat Block and adjacent terranes to the north, reducing the impact of sedimentary free-surface multiples.

The thick sediment cover of the Yakutat Block introduced multiples that complicate the identification of subsurface features in the P-wave data. An important feature of S-wave to P-wave conversions is the results are not contaminated by multiples (Kumar, 2005). However, the lower frequency content of the original S-wave data limits the vertical resolution of the S-wave data compared to the P-wave data. Due to these variations in technique, our interpretations rely heavily on the identification of features common to our P-wave image and our local and regional S-wave images, and where they overlap the tomography models of Eberhart-Phillips et al. (2006).

Our goal is to determine the geometry of the convergent plate boundary between the subducting Yakutat Block and the overriding North American Plate and the Yakutat Block and Pacific plate in southeastern Alaska. We work from the premise that subduction is taking place along this boundary and constructed a 3-D model of the subducting slab top to help understand the geometry of the region.

We used our data and the U.S. Geological Survey slab model (http://earthquake.usgs.gov/research/data/slab/) to develop a kinematic model of the top surface of the subducting plate. The purpose of this model is to provide a geometric framework to understand the subduction geometry. The method used to construct this model is similar to that described in Pavlis et al. (2012a) for the Farallon slab under the western United States. The objective is to produce a 3-D surface parameterized by two generalized coordinates: the x₁ axis is a curve following flow lines defined by a prescribed model for rigid plate motion, and the x₂ axis defines a position marker for the plate at some constant time (the present). Curves defined by constant x₁ are thus approximate time contours and curves of constant x₂ are flow lines in the limit of rigid plate motion. Here, these generalized coordinates are related to Earth coordinates through three choices.

1. We use the plate motion model of Doubrovine and Tarduno (2008) for motion of the Pacific plate relative to North America. We experimented with other plate models including current motion models defined by a single pole of rotation. Variations in plate motion are small between different models for the past 5 m.y. Animation 1 shows that the 5 m.y. contour includes the entire STEEP study area. Therefore, this choice affects flow lines only in central Alaska.

2. The time 0 marker (x₁ = 0 curve) is was defined in three sections. (1) A western section is defined by the trench axis. This section is defined as a 3-D curve extending westward from an origin southeast of Middleton Island. The depth of this curve at digitized positions was defined from bathymetry as the water bottom. The eastern limit of this segment (59.058°N, 145.73°W) is defined by the mapped convergence of faults of the Pamplona zone into the updip limit of the trench. (2) An eastern section extends from a point under the Samovar Hills (60.2°N, 140.7°W) to a point south of the Bering Glacier (59.38°N, 143.63°W). The eastern limit of this segment is defined by the mapped eastern limit of the Malaspina fault. The west endpoint is defined by where the southern limit of the Pamplona zone meets the edge of the Yakutat block. The actual 3-D curve is define by a downdip projection from the mapped surface trace of the southern edge of active faults in the Pamplona zone identified from STEEP data (Worthington et al., 2012; Bruhn et al., 2012). We took that mapped surface curve and projected it downdip at 45° in the northward direction perpendicular to local strike. We chose to project this position to a 15 km depth for two reasons. First, to the west data from STEEP01 (Worthington et al., 2012) show this as the approximate depth of sediments. We know from those data, however, that within the Yakutat Block the sediment package thins to the east to near zero at the southern end of STEEP01. The eastern side of the model is thus based on a conjecture that the faults terminate near the brittle-ductile transition, which can be expected around this general depth. The point is that the model geometry is solid at the western edge of this segment where it is constrained by data from STEEP01, but is more questionable to the east in the vicinity of Mount Saint Elias. (3) The final segment was produced like the eastern segment (anchored by Pamplona zone mapped faults), but the projection depth was tapered from 15 km to 0 at the eastern edge of the trench (section 1) segment.

3. Control points used to define the top surface of this model within the subsurface were derived from three sources. (1) We defined a set of control points using the digitized contours in the U.S. Geological Survey slab model, which extends only to about Prince William Sound where the mantle seismicity falls off dramatically. (2) We used our data to extend the model eastward. For this purpose we used a custom picking procedure to define control points consistent with the P-wave and S-wave receiver function data and the Eberhart-Phillips et al. (2006) tomography model. The clearest constraints come from the S-wave receiver function data (Animation 2). The P-wave receiver function data are corrupted by multiples in Yakutat sediments that create an ambiguity in what peak was the Moho. In addition, the migration method we used can create artifacts with the irregular coverage of this network (Pavlis, 2011a). These are clearly present in Animation 2 as so-called “smile” artifacts. However, because of a difference in frequency content of the P-wave and S-wave data the P-wave data improve vertical resolution by about a factor of 3 compared to the S-wave data. We thus adopted the approach to use P-wave data picks where the two were consistent and not prone to coverage artifacts (smiles in the data). In other areas the S-wave data were used as the main constraint. The tomography data were used in somewhat the opposite way. That is, the vertical resolution of the tomography data is much poorer than either receiver function images. Consequently we used the tomography model image largely as an auxiliary constraint to guide each pick. The reliability of these results can be appraised by examining Animations 2–4. These control points were assumed to mark the depth of the Moho at the base of the downgoing lithosphere linked to the Yakutat block. We then defined top of slab as an upward projection of each control point by 15 km. This was used based on results from STEEP01 (Worthington et al., 2012). (3) The points used to define the time 0 curve (see preceding) were treated as control points. The complete ensemble of control points were then interpolated with a mesh defined by a masked Delaunay triangularization. The mask was necessary to avoid artifacts from the convex body requirement of Delaunay triangularization and is why the model has features that are not outlined by a convex polygon.

The model is illustrated in Animation 1 and is used as a framework for other animations and most of the figures. Flow lines in the model were computed at 10 km intervals along the time 0 curve and traced to follow the surface defined by the Delaunay triangularization through all the control points. The flow lines are computed at 0.1 m.y. intervals. The animations and figures show only 1-m.y.-interval, constant-time contours to avoid clutter. The result can be thought of as a regular mesh in the x₁-x₂ generalized coordinates sampled at 10 km along the x₂ axis and 1 m.y. intervals along the x₁ axis.

The geometry of this system is unambiguously three-dimensional. We urge readers strongly to begin by carefully examining Animations 1–5. Figures 4–5 were constructed by annotating frames from these animations. Figure 4 presents five cross-sectional slices parallel to the direction of Yakutat plate motion extracted from the 3-D data volumes at key locations. The cross-section A-A′ shown in Figure 4 is located clearly on the Pacific plate. The B-B′ location represents an area of crust at the transition between the Pacific to Yakutat plates. Cross-sections C-C′ through E-E′ are located on the Yakutat Block progressively closer to the eastern boundary of the block. Figure 5 presents a cross-sectional slice oriented roughly perpendicular to the direction of Pacific plate motion and is located far enough from the trench to image the subducting slab and crust as separate entities. These locations were chosen for the following reasons. Kayak Island is on the westernmost portion of the Yakutat Block that we can image well with our P-wave method. Icy Bay is located on Yakutat Block crust midway from the end of the Aleutian Trench and the Fairweather fault. The Yakutat Bay cross section represents the easternmost portion of the Yakutat Block that we can image well and is located near the eastern edge of the Yakutat block.

We suggest the following observations are apparent from our imaging results.

1. The Pacific plate and Yakutat Block are imaged well in the area east of Prince William Sound (Animations 2–5; Figs. 4 and 5). This area contains relatively dense station coverage and other data constrain the geometry of the slab. In particular, these data show that the U.S. Geological Survey Slab 1.0 plate model correlates positively with the dipping horizons imaged by the regional S-, local S-, and P-wave data. However, this cross-validation shows that the positive peaks in the P-wave and S-wave images that are our anchor points are linked to the Moho of the subducting lithosphere.

2. Additional cross-validations of our results are illustrated in Figures 6–8 and Animation 5. The subducting Yakutat Block imaged in the tomography model of Eberhart-Phillips et al. (2006) correlates well with the geometry of the slab imaged in our regional S-wave data (Animation 5; Fig. 6). Note that the result illustrated in Figure 7 is an example of the impact of widely spaced stations. The piercing points of contours from the U.S. Geological Survey slab model intersect a nearly flat positive conversion at a depth of ∼70 km. The stair step is almost certainly an artifact of low station density in the vicinity of that section.

3. The BEAAR experiment produced a seismic profile (Ferris et al., 2003) that is also within the area incorporated in our regional S-wave image volume. The location of a low-velocity zone identified in the BEAAR experiment has been interpreted as the top of the subducting Pacific slab and alternately as the top of the subducting Yakutat block. Figure 9 compares the regional S-wave image cut along their profile. Our data image the slab poorly in this region due to sparse coverage. Updip from the BEAAR profile, however, there is a dipping feature consistent with the assertion that this feature is the Moho of the Yakutat Block and not the top of the slab.

4. The TACT experiment collected active source seismic data from stations to the east and north of Prince William Sound (Fuis et al., 2008). This is near the edge of our local P-wave and S-wave images, but well within our regional S-wave image. The correlation between the subducted Yakutat slab in our regional S-wave data and the feature identified as the megathrust in Fuis et al. (2008) indicates the high quality of our data and imaging methods as shown in Figure 8.

5. Another problem induced by poor station coverage is present in the northeastern portion of the local STEEP area image volume (Fig. 4). In this region the only stations east of Icy Bay are concentrated along the coast. This geometry creates a hole in our images that is especially bad in the P-wave data. The lack of coverage in this region makes the exact location and geometry of the eastern edge of the Yakutat Block impossible to identify with these data; however, we resolve a systematic change in structure from west to east (Fig. 4). In the west a single dipping surface is well imaged in the S-wave data with irregular sampling of the same feature seen in the P-wave data. To the east this transitions cleanly to an abrupt crustal thickening where the section intersects the Pamplona zone.

6. At the northern edge the STEEP image volume of the subducting Yakutat Block appears to extend beyond the limits of the volume. The Wrangell volcanic field is found at the northern boundary of this volume. Projection of the top of slab surface under the Wrangell volcanic field suggests that the depth to the top of slab there is ∼80 km.

7. The crust is a continuous surface across this entire study area with no apparent offsets (Animations 1, 4, 5; Fig. 5).

A primary implication of our results is the ability to define the geometry of the subducting Pacific plate and Yakutat Block and the lateral relationship between them. A primary issue is the continuity of the slab. In particular, we can at least partially address the tectonic underplating conjecture of Fuis et al. (2008), who defined that region by a translation of a map figure of the model for the original geometry of the Yakutat Block suggested by Plafker and Thatcher (2008). Fuis et al. (2008) suggested the presence of a slab tear along a line northwestward (parallel to plate motion) from the eastern side of Prince William Sound, as seen in Figure 2; they defined this boundary based on the rapid drop in seismicity east off this line. Fuis et al. (2008, p. 269) stated “…that the Yakutat terrane has since been removed from this position atop the Pacific oceanic lithosphere in the region of the Wrangell volcanoes by tectonic underplating.” Our results support some parts of this model and reject others.

We first address areas of agreement. Examine the outline of the Yakutat Block shown in Figure 2 of this paper or Figure 1 of Fuis et al. (2008) and play Animations 1 and 5, and pay careful attention to the blue band commonly used as a measure of the LAB (Kumar et al., 2001) that is seen everywhere in the S-wave receiver function volume at depths between 50 and 150 km. The conjecture that this feature is measuring LAB is reinforced by the surface projected 70 km below the top of slab surface. That is, the animation shows that this feature (LAB) tracks downdip, loosely following the top of surface projection. The surface is smooth where coverage is good, but gets shingled in areas of sparse coverage and steeper dip (e.g., under Cook Inlet). The key observation, however, is that seen in Figure 5, where there is an intriguing overlap of LAB surfaces in the region predicted by Fuis et al. (2008) under Prince William Sound. This observation of overlapping LAB is tempered by the lack of direct coverage in that area. However, our results indicate the LAB surface is relatively flat in this region and piercing point coverage for S-wave to P-wave conversions is good under Prince William Sound due to the dominance of sources from the western Pacific. The large number of piercing points taken in conjunction with the inferred change in velocity seen in this region in the tomography model of Eberhart-Phillips et al. (2006) (seen in parallel sections in Animations 2 and 3) suggest that this region outlines the western limit of the Yakutat Block in the subsurface.

The negative side of our inferences for Fuis et al. (2008) is that we can definitively reject the conjecture of a discontinuity in the slab under the Chugach Mountains and the hypothesized tear northwest of Prince William Sound. The best evidence for this is seen in Animation 5 and the interpretation made from it in Figure 5. We see a continuous Yakutat Moho and related top of slab surface from the coast to the south side of the Wrangells. The east-west boundary marking a transition from Yakutat to Pacific affinity lithosphere Fuis et al. (2008) hypothesized is simply not there. All our results show that the Moho and linked top of slab are continuous across this transition with no resolvable change in the geometry of the subducting slab. Animation 5 also illustrates that seismic receiver function (SRF) images show the Moho splitting into the downgoing slab Moho and base of the Alaskan crust throughout most of the STEEP study area. We note that this feature is perhaps better understood as a wedge of low-velocity material (thickening to the north and east) embedded between the base of the Chugach and Peninsular Ranges terrane crust north of the Yakutat Block surface suture and the top of the subducting lower crust of the Yakutat block. There is no lateral tear at the western limit of this wedge, but a clean transition into the pure flat slab region of central Alaska where this low-velocity wedge becomes invisible (if present) in our data.

We conclude from this that the western limit of the Yakutat Block inferred from other studies (Zweck et al., 2002; Eberhart-Phillips et al., 2006; Fuis et al., 2008) is well supported by our data (Fig. 5). The LAB inferred from the SRF data provides some support of the east-west lithospheric convergence suggested by the Fuis et al. (2008) model of lithospheric underplating; our results, however, indicate that their suggestion the slab is torn and that Yakutat Block lithosphere extends only to the vicinity of the Chugach Range is not correct. More uniform coverage of the region north of the STEEP area is needed to extend this conclusion northward. That is, it must be true that somewhere downdip along the subduction zone we should find lithosphere of Yakutat Block affinity. We suggest that the northern boundary of the subducted Yakutat Block is to the north of the area of our study.

The dip of the subducting slab increases from west to east (Animations 1–3). Near Prince William Sound and the Copper River Delta we measure a dip of 11° that increases to 16° near Icy Bay. The increase in dip is gradual and results in a smooth surface with no clear discontinuities at the scale that we can resolve. This shallow dip is consistent with observations of surface deformation noted to the north and east of the Yakutat block. This shallow slab dip is consistent with deformation features noted in regions also associated with flat slab subduction (Haschke et al., 2002). The 11°–16° dip of the subducting Yakutat Block that we image is similar to the 9° dip reported in the region under Prince William Sound (Brocher et al., 1994). It is, nonetheless, much less than the nominal 60° dip along the slab in the western and central Aleutians (England et al., 2004).

A second aspect of the subducting slab geometry that can potentially be addressed by these data is the location of the eastern edge of the subducting plate. All the animations and Figures 5 and 10 illustrate that we are imaging a continuous, dipping slab feature between Prince William Sound and Icy Bay. However, our station coverage ends at Mount Saint Elias with a huge hole in coverage in southwestern Yukon Territory. This limits strongly our ability to image the lithosphere farther east with these methods. As a result, the slab geometry is not well resolved in the region northeast of Mount Saint Elias where the plate model flow lines (Animations 2–5) would predict that the eastern edge of the Yakutat lithosphere should be present. The data suggest that a dipping horizon is still present near Icy Bay (Fig. 4, line D-D′), consistent with our subduction model for the Moho of the Yakutat block. Animations 1–5 show that the S-wave data indicate that this feature may project east as far as Yakutat Bay. There are hints of changes in the geometry of the LAB east of Yakutat Bay (Animation 5), but this should be viewed with caution as this could be explained as a smoothing artifact of the CCP stacking method we used interacting with very sparse coverage to the east. The P-wave data algorithm is more conservative in this regard and defines a null solution in this region. Therefore, we conclude that the slab can clearly be traced at least as far as Mount Saint Elias, but the precise location of the eastern edge is somewhere east of Mount Saint Elias and cannot be determined using these data. This issue, however, should be well defined by the planned future deployment of the USArray (http://www.usarray.org/) in this region.

Animations 2–4 show clearly that there is a major change from west to east in the geometry of dipping slab defined by the Moho of the downgoing Yakutat block. Our results indicate that in the west the Moho defines a relatively simple surface dipping along flow lines that steepen slightly from west to east. In Animations 2 and 3 the geometry changes strongly as the section lines moves from a coastline intersection near the Bering Glacier (Fig. 4, C-C′) to Mount Saint Elias (Fig. 4, E-E′). An intermediate state is highlighted in Figure 4 (cross-section D-D′). Figure 5 overlays our interpretations on these sections based on three factors.

1. Christeson et al. (2010) and Worthington et al. (2010, 2012) demonstrated that the Moho of the Yakutat Block is essentially flat south of the deformation front defined by the Pamplona zone. Furthermore, the reflection data demonstrate that active deformation is largely concentrated north of the Pamplona zone and with little to no active deformation on the Transition fault (Gulick et al., 2013).

2. Chapman et al. (2012) described details of the structural geology of active faults in the Saint Elias region. Simplified versions of their results are integrated in Figure 10.

3. Our Moho and top of slab models are included with minor modification in these interpreted sections.

4. A primary picture of the overall structural setting of the Yakutat Block is that the sedimentary package of Yakutat Block sediments north of the Pamplona zone have all been intensely deformed by a multitude of active thrust faults (Pavlis et al., 2012b). At the scale of our imaging this can be reduced to a package of sedimentary rocks being stripped from the surface at the top of the downgoing Yakutat lithosphere. The northern limit of that package or rocks is bound by the mapped suture in the vicinity of the Bagley Ice Field.

5. Worthington et al. (2012) demonstrated that although the Yakutat Block crust has a nearly constant thickness, the geometry has the shape illustrated in Animations 2 and 3. That is, the sedimentary package is ∼15 km thick near the coastline intersection with STEEP01, but that package thins dramatically at the southern end of that refraction line. At the same time the lower crust is thin (∼15 km) to the west and thickens to ∼30 km of oceanic plateau material to the east.

6. Shennan et al. (2009) presented evidence that prehistoric megathrust earthquakes comparable to the 1964 great Alaska earthquake (Page et al., 1989) ruptured as far east as Icy Bay, indicating a continuity of a coherent top of slab to that vicinity.

7. A primary feature of this orogen is easily missed: the distance from the leading edge of the Pamplona zone to the northern limit of Yakutat sediments (measured along the flow line direction) to the west (Copper River Delta) is ∼150 km. The comparable measure approaches zero at Mount Saint Elias. This implies a much larger strain rate in the region of Mount Saint Elias than the western margin of the Yakutat block.

Our interpretations build on these concepts and can be summarized by two end members. At the western edge of the Yakutat Block the onshore geology and offshore reflection data show that we are dealing with a fold-thrust belt formed by sediments stripped from a décollement formed by a megathrust, as seen in Figure 5 (C-C′). The Shennan et al. (2009) results indicate that this megathrust surface is continuous with that farther west, which no one would dispute is linked to the Aleutian Trench. To the east the sediment package thins and STEEP01 shows that the package of crustal rocks the system is trying to subduct consists of 30 km of oceanic plateau style crust (Fig. 5, line E-E′). Our results suggest that the geometry in this region is better viewed as lower crustal thickening in response to this thickened subducting crust geometry. The style of deformation is approaching, but not exactly, continent-continent collision. We suspect that the deformation of the active structure on the south side of Mount Saint Elias modeled by Elliot et al. (2010) penetrates a region of ductile deformation in the lower crust, where the crust is markedly thickened (Fig. 10C). Figure 10B shows simply a representative transitional state. Whether this occurs abruptly or as a transition, as implied here, is actually not resolvable by these data, but we suggest that it is a rational way to expect the system to behave, given that it is likely driven by the wedge shape of the lower crust of the Yakutat block.

An issue for this area is the relationship of the plate convergence to the Wrangell volcanic field. Given that the deepest portions of the slab project out of our image volume due to lack of station coverage north of the Chitna River, our data provide no constraints on the geometry of the subducted slab north of the Wrangell volcanic field; however, we can state that the slab is continuous to a point below the field. Immediately beneath the Wrangell volcanic field the slab is at a depth of ∼80 km, as shown in Figure 4 (lines C-C′, D-D′). This is consistent with the standard model that the subducting slab is associated with the magma that feeds the volcanoes. The depth we estimate to the slab top under the Wrangell volcanic field is at the upper limit of the 80–160 km range of normal depths to slab under volcanic arcs (England and Katz, 2010).

England and Katz (2010) developed a relationship between the dip and convergence rate of subducting slabs and the depth of the slab under a given volcanic arc. Because adjacent portions of volcanic arcs have similar geometry and convergence rates, a natural result is that consistent depths to the slab are observed along individual segments of arcs (England and Katz, 2010). The depths to the slab under the easternmost Aleutian volcanoes, Spurr and Redoubt, are ∼120 km and 150 km, respectively. Given such a large discrepancy in depth to the slab, it is clear that the Wrangell volcanic field must be considered a distinct segment of the volcanic arc in spite of similar slab dips and convergence rates. It is not surprising that the processes of magma generation differ along portions of the subducting Pacific plate and the adjacent Yakutat block.

Another feature of both the S-wave and P-wave data is the difficulty imaging the deeper portions of the subducting slab. The lack of a distinct slab signature at depths between 50 and 100 km is observed in the region to the north and west of Prince William Sound (Fig. 4, lines A-A′ and B-B′). The slab is not imaged well in our data at depths below 100 km. The plate motion models show that convergence beyond that depth cannot be denied unless something happened before the Yakutat Block arrived that has not been recognized in the geologic record. At the depth where the slab is poorly imaged, a planar, subhorizontal, high-velocity feature also appears. The horizontal feature is located ∼25 km to the north of the point where the slab image becomes indistinct. Several potential explanations for this phenomenon must be explored.

One possible explanation for this observation is that the northern edge of the slab coincides with the boundary of the image volume and with an area of poor station coverage. If so, we could simply dismiss this as an imaging artifact. However, the slab also disappears in the regional S-wave data near Anchorage where coverage is reasonable. Furthermore, Ferris et al. (2003) analyzed the BEAAR data and saw a similar loss in P-wave receiver function amplitude that they interpreted as the result of low seismic velocity contrast between the subducted slab and adjacent mantle. These observations lead us to suggest that the inability to image the slab well is in fact due to a geological process and not merely a data collection and processing artifact. As a result, it is important to give examples of possible physical explanations for this feature.

The deep portions of subducted slabs have been difficult to detect with receiver function methods in subduction zones outside of Alaska. One notable example of a location where the phenomenon of poorly imaged slabs has been observed is in the Cascadia subduction zone of western North America. A horizontal, planar feature has also been noted in conjunction with the point where the slab exhibits little seismic impedance contrast in Cascadia as well as Alaska (Bostock et al., 2002).

The physical cause of the partial transparency of the deep slab in the Cascadia subduction zone has been attributed to physical and chemical changes within the crust of the subducting slab. It has been hypothesized that the crust of basaltic composition has been eclogitized due to the interaction of fluids released by the slab with the heat and pressure at depth (Peacock et al., 2011). The pressure and temperature conditions associated with eclogite facies metamorphism occur at relatively shallow depths of 45 km. The eclogitization of the slab in conjunction with raised temperatures could reduce the seismic impedance contrast between the slab and surrounding mantle and crust. This reduced seismic impedance contrast would weaken seismic conversions at that boundary.

An important difference between Cascadian and Yakutat subduction geometry is the depth at which the slab becomes transparent. The transition between a distinct slab and transparent slab occurs at a depth of 45 km in Cascadia versus a depth of ∼75 km in southeast Alaska. The depth difference between Cascadia and Yakutat could indicate different thermal structures or fluid contents.

An alternative explanation for the transition to a transparent slab in the receiver function images could be linked to the subduction of a small portion of sediment with the slab. The slab becomes transparent at a depth that coincides with a horizontal feature that could be the signature of subducted sediments that have detached from the slab. Numerical models indicate that buoyant forces acting on sediments are capable of detaching deposits as thin as 350 m from the slab (Currie et al., 2007). Under normal subduction zone thermal structures, the detachment takes place at depths near 100 km. Although the Currie et al. (2007) model predicts that increased temperatures would result in sediment detachment at shallower depths, the temperature required to detach at 75 km, as found in our data, is not reported. This process could be involved in the formation of this horizontal feature that serves to obscure the slab in our data. Although the majority of Yakutat sediment cover does not appear to be subducted with the slab, given the resolution of our data it is not possible to state conclusively that <350 m of sediment is subducted.

This study incorporates migrated receiver functions to characterize the geometry of subducting Yakutat Block displaced terrane under southeast Alaska. These new data allow us to trace the top of a continuous slab from its established location under Prince William Sound to the edge of our coverage near Yakutat Bay. This extends the North Pacific subduction zone from the end of the Aleutian Trench an additional 300 km to the vicinity of the Fairweather–Queen Charlotte fault system. The slab dips 11° near Kayak Island and steepens toward the east with a maximum dip of 16°. This dip produces ∼100 km of separation between the Yakutat Block and the North American Plate under the Alaska Range. This amount of separation implies that the deformation found along the Denali fault is not caused by coupling between the subducted Yakutat Block and North America. Due to poor seismic station coverage we are unable to locate the eastern edge of the subducting plate with certainty. The surface expression of the boundary between the Yakutat Block and North American Plate extends from the eastern terminus of the Aleutian Trench, along the front of the Pamplona deformation zone, through Icy Bay, and to a point under the Seward Glacier on the east side of Mount Saint Elias. The downdip portion of the subducting slab system extends to the north at least to the bounding edge of the STEEP image volume and includes the Wrangell volcanic field. Our kinematic plate model clarifies that the subducted Yakutat slab is present under the Wrangell volcanic field and is at the proper depth to be the possible source of magma that feeds the Wrangell volcanic field. We image the crust and upper mantle at lithospheric scales to argue that the geometry of the crust-mantle system varies systematically from west to east. To the west sediments on the upper part of the Yakutat Block are being stripped to form active structures seen on land along a megathrust. We argue from the geometry of the system that the megathrust under this region is continuous with that associated with Pacific plate subduction to the west. Moving to the east, the subduction system transitions into a zone of intense shortening in the vicinity of Mount Saint Elias. There the distance over which shortening occurs drops dramatically, and the system appears to be deforming through ductile deformation and thickening of the lower crust. We suggest that this change in the style of deformation is driven by crustal structure of the Yakutat Block imaged by Worthington et al. (2012) from STEEP01.

Professor Terry Pavlis of the University of Texas-El Paso had the vision to bring together the wide array of scientists that contributed to the STEEP project. We also thank Sean Gulick, Gail Christeson, and Ben Hooks for their constructive suggestions. The generous access to the Alaska Earthquake Information Center (AEIC) database was invaluable.