A well-recorded Mw 7.8 megathrust earthquake occurred on 27 October 2012 under the Queen Charlotte terrace off the west coast of Haida Gwaii, western Canada. In this study, we supplement limited earlier seismic refraction work on the offshore velocity structure off Haida Gwaii with data from ocean bottom seismometers (OBS) operating between 6 December 2012 and 5 January 2013. The OBS recorded a portion of the aftershock sequence, and an active-source seismic survey was conducted in January 2013 to acquire seismic refraction data in the region of the Haida Gwaii earthquake across the Queen Charlotte terrace. P-wave velocity analyses using first-arrival tomography showed relatively shallow (2.0–3.0 km below seafloor) high-velocity material with values up to 4.0 km/s beneath the terrace. At the one OBS station seaward of the deformation front on the abyssal plain, refraction velocities of ∼4.5 km/s indicated the top of the oceanic plate at ∼1–2 km below the seafloor. At several OBS stations, converted shear-wave velocities were determined within the sediment cover using reflected arrivals. The S-wave velocities ranged from 0.5 to 1.5 km/s, and the corresponding P/S velocity ratio was between 3.0 and 4.2. The new refraction data confirm earlier interpretations of high-velocity material in the shallow terrace that may indicate fractured oceanic crustal material lies significantly above the location where a subducted slab is thought to occur under the terrace. Transpressive deformation of the Pacific plate may explain these observations.

This seismic structure study was conducted offshore of the Haida Gwaii archipelago on the west coast of Canada, mostly seaward of the Queen Charlotte fault. The Queen Charlotte fault is a primarily right-lateral transform boundary between the continental North America plate and oceanic Pacific plate. The Queen Charlotte fault is one of the most seismically active strike-slip faults globally, with a slip rate of ∼55 mm/yr (DeMets and Merkouriev, 2016; Brothers et al., 2020). The Queen Charlotte fault and its northern extension, the Fairweather fault, have had repeated large strike-slip earthquakes of magnitude Mw 7 and above, notably the M = 8.1 strike-slip event in 1949 (e.g., Rogers, 1982; Ristau et al., 2007). Along the southern portion of this transform margin, in the region of Haida Gwaii (formerly Queen Charlotte Islands), the relative motion of the North America and Pacific plates has a convergent component up to 15° at 15–20 mm/yr over the past 6–10 m.y. (Rohr and Dietrich, 1992; Hyndman and Hamilton, 1993; Atwater and Stock, 1998; Rohr et al., 2000; DeMets and Merkouriev, 2016; Brothers et al., 2020).

Shallow-penetration seismic reflection data have imaged compressional structures in the Queen Charlotte terrace, which have led to comparisons to sedimentary structures above subduction zones (e.g., Davis and Seemann, 1981; Riedel et al., 2020). Larger-size seismic sources were last used in 1983 off Haida Gwaii for analyzing the velocity structure of the region during a crustal-scale onshore-offshore refraction experiment (see location on Fig. 1; Horn et al., 1984; Dehler and Clowes, 1988). The seismic sources were a 32 L (∼1950 in3) air gun across the Queen Charlotte terrace, and 32 explosive shots in the offshore area west of the Queen Charlotte terrace. Using these large sources together with two ocean bottom seismometers (OBS) and eight land-based seismographs, Dehler and Clowes (1988) were able to construct a detailed crustal-scale velocity model extending to 20 km depth, an extended version of the initial model published by Horn et al. (1984). Yet, these two seismic refraction studies in the southern Queen Charlotte terrace, the area of greatest convergence, did not image a simple slab beneath the terrace but found a zone of complexly faulted blocks with velocities well above those expected for a purely sedimentary section (Horn et al., 1984; Dehler and Clowes, 1988).

Substantial studies have been conducted on the Queen Charlotte margin since the 2012 Mw 7.8 thrust earthquake (James et al., 2013, 2015), which may have been related to oblique underthrusting of the Pacific plate underneath North America (e.g., Hyndman, 2015), and the 5 January 2013, Mw 7.5, strike-slip Craig earthquake to the north (e.g., Aderhold and Abercrombie, 2015; Holtkamp and Ruppert, 2015; Tréhu et al., 2015; Walton et al., 2015). The oblique convergence at the margin appears to be partitioned into nearly pure strike slip on the Queen Charlotte transcurrent fault just off the coast and nearly orthogonal convergence on a thrust fault beneath the Queen Charlotte terrace (Hyndman, 2015), called the Haida Gwaii thrust. The Mw 7.8 event has been interpreted to have occurred on this thrust, with rupture dipping landward and extending to approximately beneath the coastline (Kao et al., 2015), although its depth is not well constrained. This event was Canada’s second-largest instrumentally recorded earthquake. It occurred at 8:05 p.m. (local time), 27 October 2012, off the west coast of Moresby Island, Haida Gwaii (Figs. 1 and 2). In the first month after the earthquake, more than 20,000 aftershocks were recorded (Mulder et al., 2013; Farahbod and Kao, 2015), the majority of which were less than magnitude 2. Most of the aftershocks were located in an area ∼120 km in length and ∼80 km in width along the margin (Fig. 1). They extended seaward from beneath the west coast of Moresby Island to >50 km offshore just seaward of the trench-like bathymetric feature called the Queen Charlotte trough (Figs. 1 and 2).

The region of the 2012 aftershocks was the target of a marine seismic expedition using 14 OBS. The deployment was in December 2012, with recovery in January 2013 (Riedel et al., 2014a, 2014b). The OBS station distribution was centered on the 2012 epicenter and covered most of the area of the aftershocks beneath the Queen Charlotte terrace and trough region (Fig. 2). The objectives of the OBS experiment were twofold: (1) to use OBS to obtain better aftershock locations, and (2) to use active-source seismic refraction data to obtain a detailed seismic velocity structure of the Queen Charlotte terrace. One of the critical questions off Haida Gwaii is the potential processes of underthrusting (or subduction) of the Pacific plate underneath North America, or if convergence is accommodated by other processes, for example, internal shortening or basin inversion within the two plates (e.g., Brothers et al., 2020; Schoettle-Greene et al., 2020). Constraining the velocity structure and style of deformation of the Queen Charlotte terrace could help in answering this critical question. The more recent active-source seismic experiment in 2013 consisted of two surveys conducted with a single smaller 8.5 L (520 in.3) G-gun and a single-channel streamer. The first survey was not recorded on the regional grid of OBS stations deployed to record the aftershock sequence, because of their limited battery capacity. A second survey was then carried out with six refurbished OBS (Fig. 2). Four crossing lines were acquired across the six OBS stations with the same seismic source and streamer. Four of the six OBS were located along the main SW-NE–oriented transect across the Queen Charlotte terrace (Fig. 2).


Seismic structure studies along the Queen Charlotte and Fairweather margin have been conducted since the early 1960s due to the recognized earthquake hazard potential in the region (e.g., Chandra, 1974; von Huene et al., 1979; Rogers, 1982; Carlson et al., 1988; Nishenko and Jacob, 1990). Seismic reflection surveys have imaged the striking morphology of the slope (Chase and Tiffin, 1972; Davis and Seemann, 1981; Rohr, 2015), including a margin-parallel trough at the base of the slope, the Queen Charlotte trough, and the 25-km-wide bench-like feature of the Queen Charlotte terrace, where folded sediments were imaged within the upper 500–800 m below the seafloor (Fig. 3). There are no reflection seismic data available that image the top of the oceanic crust beneath the Queen Charlotte terrace, including data from the multichannel seismic (MCS) EW9412 survey (Scheidhauer et al., 1999, 2015; Tréhu et al., 2015).

Landward, the Queen Charlotte terrace is bounded by a steep slope rising to the ∼1-km-wide continental shelf (Barrie et al., 2013); seaward, the Queen Charlotte terrace descends into a sediment-filled trough (Fig. 3). The seaward slope is steep and incompletely imaged by the legacy seismic data (e.g., as compiled by Davis and Seemann, 1981) and our newly acquired single-channel seismic data from the 2013 survey (Fig. 3). This morphology is similar to that of subduction zone trench and accretionary sedimentary prisms (Hyndman and Hamilton, 1993). An alternative explanation is that transpressive deformation has led to folding, uplift, and deformation (similar to an accretionary wedge) that loads the Pacific plate under the Queen Charlotte terrace and generates a similar morphology (Harris and Chapman, 1994; Prims et al., 1997; Rohr et al., 2000; Tréhu et al., 2015; Walton et al., 2015). Refraction models across the Queen Charlotte terrace have not detected a simple dipping slab; instead, a complex suite of arrivals necessitates a model of complex structure in the Queen Charlotte terrace. The landward extension of the oceanic plate dipping beneath the islands has been interpreted by seismic receiver functions using distant earthquake sources (Bustin et al., 2007). Several structural models explaining the complex tectonics have been proposed (Horn et al., 1984; Dehler and Clowes, 1988; Mackie et al., 1989; Rohr et al., 2000; Tréhu et al., 2015; Hyndman, 2015; Walton et al., 2015; Brothers et al., 2020). We show the original model by Dehler and Clowes (1988) in Figure 4 for comparison and discuss the interpretation of the units imaged.

The modeling by Dehler and Clowes (1988) and Mackie et al. (1989) showed an ∼1-km-thick sedimentary unit on the oceanic plate with velocities of 1.8–2.2 km/s. Oceanic crust on the Pacific plate (west of the trough, profile distance 0–105 km; Fig. 4) has two layers, with an upper moderate-velocity unit (3.8 km/s) overlying high-velocity material (6.7 km/s). The upper-crustal layer was interpreted to consist of “basaltic pillow lavas and sheet flows grading to dike swarms” (Dehler and Clowes, 1988, p. 1868, their table 2). The lower-crustal unit was interpreted to consist of normal oceanic interlayered gabbro. The terrace region (105–140 km profile distance; Fig. 4) is marked by steeply dipping faults, which vertically displace individual units of moderate velocity similar to the upper-crustal layer seen further to the west. These moderate-velocity units are either displaced crust or consist of highly compressed sediments. They are overlain by a sedimentary unit of variable thickness that increases to 3 km across the terrace, where the sediment layers are folded and faulted (Fig. 4). At the Queen Charlotte fault near the coast, two units of moderate-velocity are juxtaposed, and, based on the higher-resolution air-gun data, they are separated by a narrow, low-velocity zone (3.4 km/s), possibly indicating fault shear deformation (Fig. 4). The deeper portion of the model beneath the terrace is composed of a high-velocity unit (5.3 km/s) with a steep velocity gradient. The thickness of this unit increases to the east, with velocity approaching 7.3 km/s at its base. This unit may consist of deformed (sheared) gabbros (Dehler and Clowes, 1988). However, the internal structure was not resolved (Fig. 4). The continental portion of the transect (150–200 km distance of the profile shown in Fig. 4) consists of an upper, ∼5-km-thick unit with an average velocity of 5.0 km/s, interpreted as the rocks outcropping on Haida Gwaii (Sutherland Brown, 1968; Horn et al., 1984). Underneath this, there is the main uniform continental (silicic) crust with an average velocity of 6.1 km/s (Fig. 4). More recent work in the adjacent Queen Charlotte Basin (Spence and Asudeh, 1993) has led to revised interpretations of velocities underlying Haida Gwaii to higher values (6.4–6.8 km/s), which are more suggestive of mafic compositions.

Seismic Reflection Data from 2013 Survey and Legacy Seismic Data

The 2013 OBS survey employed a single 8.5 L (520 in.3) G-gun. The compressor capacity allowed the air gun to be fired at an interval of 30 s, with a pressure of 1800 psi (12.4 MPa). The shot spacing was defined by the time interval and an average vessel speed of ∼4 knots, giving a shot spacing of ∼70 m.

The streamer consisted of 48 individual Teledyne B-1 acceleration canceling hydrophones, bundled to six channels each with eight hydrophones. These six channels were then combined at the amplifier to a single channel output. In order to maintain a uniform tow depth, three passive guide birds were attached to the streamer. The tow depth was maintained to between 1 and 3 m. The streamer had an initial ∼9-m-long dead section, followed by the ∼45-m-long active section, and a second ∼9-m-long dead section. Seismic processing included geometry definition, band-pass filtering (6–120 Hz), and amplitude recovery to account for spherical divergence. The single-channel seismic data were used to constrain the structural model and the approximate location of the main refraction events, as described below. Several single-channel seismic reflection profiles acquired across the study region by the Geological Survey of Canada from 1967 to 1978 (Davis and Seemann, 1981; Hyndman et al., 1982; Rohr, 2015; Riedel et al., 2020) also were incorporated into the structural model. Original legacy data that were recorded only on printed paper were scanned and converted to trace format for interpretation. No further processing was applied to the digitized records. The locations of these older lines were poorly known due to low navigational accuracy, having been obtained prior to the global positioning system (GPS). However, these data served as a guideline for picking basement reflections and structure across the Queen Charlotte terrace in the region of the new OBS experiment.

OBS Relocation, Preprocessing, and Traveltime Picking

The 14 OBS instruments (each equipped with a three-component short-period geophone and one hydrophone) were provided from the instrument pool maintained by Dalhousie University, Nova Scotia, and the Geological Survey of Canada. The sample rate was 4 ms, giving a Nyquist frequency of 125 Hz. For the active-source OBS experiment, we used six of these 14 instruments. We had no time-server with which to autocorrect the times of the shot-logging computer and trigger, so a manual fix of the shot-time logger was applied. Manual synchronization of the shot-time logger was done every 10 min (600 s). The OBS instruments accumulated a relative clock-drift, which was corrected after instrument recovery through synchronization to the GPS clocks (Riedel et al., 2014b), assuming a linear drift. After an initial shot time definition, direct arrival times were picked from all OBS stations to carry out an OBS relocation process using the technique introduced by Zykov (2006). This relocation method allows inclusion of piece-wise linear unknown time-drift corrections not captured in the (onboard manual) shot-time synchronization, as well as air-gun trigger time delays, unknown bias in the water-column sound velocity profile, and uncertainties in the shot (air gun) position. Details on the relocation analysis and inversion have been described in Riedel et al. (2014b). The average OBS instrument drift was ∼200 m to the SW relative to the initial drop location from the vessel (water depths ranging from 1850 to 2880 m), likely due to local currents. After relocating the OBS stations, new offset (distance) information was calculated and used in the seismic refraction velocity analyses presented in this study.

The OBS hydrophone and geophone receivers all showed abundant but mostly low-frequency (<6 Hz) seismic arrivals from aftershocks, so the data were band-pass filtered to increase the signal-to-noise ratio (S/N) and enhance the air-gun–shot arrival frequencies. No deconvolution operator was applied to suppress the air-gun bubble. First arrivals for the P-wave analyses were determined from the hydrophone and/or from the vertical geophone component. Often, identification (and thus picking) of arrivals was visually easier on sections processed with a reduced velocity (2–4 km/s). Additional attributes (e.g., phase only) helped in identifying lower-amplitude but laterally consistent arrivals.

S-wave arrivals were determined from the horizontal and vertical geophone components in comparison to the hydrophone. The converted S-wave arrivals (typically with a lower-frequency content than P-wave arrivals) were identified based on the characteristics of phase reversals of reflections across the vertical incident point beneath each station. Prior to the final S-wave arrival picking and velocity analyses, the data were rotated into radial, transverse, and vertical orientations using Seismic Unix, with the average geophone orientation determined using the technique by Rosenberger (2010), described in detail by Riedel et al. (2014b).

Refraction First-Arrival Tomography

We built an initial structural model using constraints from the models by Dehler and Clowes (1988) and Mackie et al. (1989) and all available seismic reflection data (Fig. 2). These features included top of crust on the abyssal plain, location of main basin-and-ridge structures of the terrace, and associated sediment layering. Initial P-wave velocities for these main structural features were determined by applying a suite of reduction velocity values (Figs. 5 and 6) until the arrivals were flattened. The velocity associated with a flattened arrival was assigned as a starting value to the corresponding layer. Velocities (and the velocity gradient) between the seafloor and this first layer with a refracted arrival were constrained using reflections seen on the OBS data as described in Riedel et al. (2020).

Forward modeling was then conducted for first arrivals at the OBS stations, and a best-fit model was iteratively determined using a layer-stripping approach. This best-fit model from forward modeling was then used in an inversion to improve the misfit between observed and calculated first arrivals. The method of damped least squares from the RAYINVR package (Zelt and Smith, 1992) was used for this inversion. Picking uncertainties for individual layers were dependent on OBS sampling rate and overall S/N ratios and varied between ±2 ms for the direct arrival (through the water column) and ±25 ms for the deepest arrivals. The RAYINVR package was also used to verify the P-wave velocity field by predicting through forward modeling the first seafloor receiver multiple (using OBS 2 and OBS 6) as well as one peg-leg (using OBS 6). However, those raypaths were not incorporated into the inversion.

The OBS stations also recorded arrivals from shots during turns and parallel profiles (Fig. 2) not along the main two-dimensional (2-D) profile. We picked the events, calculated a projected offset, and superimposed those arrivals on the main 2-D profile as additional constraints on the velocity and depth structure, assuming that the underlying structures do not vary much over lateral distances of a few kilometers.

Velocity Analyses

The OBS experiment from the 2013 survey provided important new structure data but was limited by the small air-gun size (limited depth penetration), and station and shot spacing (offset limitation), as well as high-amplitude noise from the many aftershocks of the 2012 Haida Gwaii Mw 7.8 event. An additional challenge, given the low ray coverage, was the steep rise in topography from the trench to the terrace (∼900 m over a distance of ∼6.5 km) and the complex X-ray paths. All parameters and fit of the model to picked arrivals are listed in Table 1. A graphical representation of the goodness of fit of the model and rays to the picked events is illustrated in Figure 7. In the figure, selected rays (every 10th) are shown for each of the types of rays modeled to avoid cluttering. The data suggest a refracted arrival from within the upper sediment layer (>500 m sub-seafloor depth) with a velocity around 2 km/s. A second consistent refracted arrival with a velocity of ∼3.7–4.7 km/s was seen at all stations at a depth of ∼1.0 km sub-seafloor. At the westernmost station (OBS 6; see location in Fig. 2), a far-offset arrival toward the western portion of the profile across the stations gave a velocity of ∼4 km/s, indicating arrivals from the top of oceanic crust, and a velocity of ∼6 km/s, indicating arrivals from upper oceanic crust (Fig. 6).

After initial forward modeling and a simple visual fit of the model to the picked arrivals, we updated the velocities by applying the damped least-square inversion. However, this improved the fit only slightly. A detailed depiction of the modeled and picked arrivals for each of the rays listed in Table 1 is given in the Supplemental Material1 (Figs. S1 to S6). Due to the limited number of rays identified overall, we did not perform a resolution (checkerboard) test, and we limited the interpretation to the laterally consistent events described above.

The final velocity model is shown in Figure 8 as a smooth color grid. The structural interpretation and velocities along the major boundaries are shown in Figure 9. We resolved two sedimentary layers with consistent velocities of 1.6–3.0 km/s. Refraction from the uppermost sediment layer 1a was best identified at OBS station 2. At all other stations, the arrivals for these layers were masked and close to that of the direct arrival, so they could not be picked with certainty. However, reflections from the bottom of sediment layer 1a were identified at OBS 1, 2, and 6 (Fig. S3) and incorporated into the model. A velocity reduction near the Queen Charlotte fault was identified, as previously observed by Dehler and Clowes (1988).

In general, vertical incidence reflections provide good control on the layer depths, whereas the refractions constrain the layer velocities. Average picking uncertainty for reflection events in this study was 12 ms, which corresponds to an uncertainty of ∼15 m in layer boundary depth for an average velocity of 2.3 km/s for the sediment layers. Similarly, the maximum depth uncertainty in the crustal layer was ∼100–125 m for an average velocity of 4–5 km/s. This is a reasonable estimate considering that the model extends to 6–8 km beneath the seafloor.

Validating Model with Peg-Leg and Multiple Arrivals

In the velocity inversion described above, we only incorporated the regular first arrivals. However, multiple arrivals were abundant in the data, especially at OBS 2 and OBS 6 (Figs. 5 and 6). Using the velocity model from first arrivals, we predicted arrivals for the receiver-multiple (extra reflection in the water column near the OBS) and compared them to the picked arrivals of this multiple (Fig. 10). These raypaths sampled a slightly different portion of the subsurface and thus could be used to validate the appropriateness of the velocity model for those portions of the model not equally covered by the regular first arrivals. The multiple was predicted reasonably well for both OBS stations, indicating a robust velocity model. Additionally, we predicted a peg-leg where the rays reflected back into the subsurface near the OBS station and traveled an extra distance through sediment layers 1a and 1b (Fig. 11). Here, we only used the eastern portion of OBS 6, where a clear peg-leg was identified in the data. Again, the prediction of the peg-leg was within the picking uncertainty, indicating a robust velocity estimation. We also predicted a peg-leg for only upper sediment layer 1a at this OBS, but no arrival could be seen in the data, suggesting an overall weak impedance contrast between the two sediment layers.

Converted-Wave S-Wave Velocity Analyses

The OBS data also allowed us to determine sediment S-wave velocities. The S-wave velocity analyses used the ray-tracing code RAYINVR (Zelt and Smith, 1992) by defining boundaries for P-to-S conversion from the initial P-wave model and assigning an initial Vp/Vs ratio to the layers. For the uppermost sedimentary layer 1a, we started with a value of 4 for the Vp/Vs ratio. We iteratively modified this value (by changing the S-wave velocity) until arrivals matched those of picked S-wave reflection events. The horizontal geophone components were often ringy; i.e., they suffered from reverberations of incoming S-wave arrivals (see, e.g., Fig. 12). Thus, picking of S-wave reflection events was done sparsely, and a new converted-wave reflection was identified only if the move-out of that arrival was different from the reflection arrival immediately identified above. An example of S-wave arrivals on the three geophone components (after rotation) and hydrophone of OBS 4 is shown in Figure 12. Note, the hydrophone data do not show any S-wave arrivals, which serves as verification of our interpreted S-wave arrivals on the other geophone components. We only defined S-wave arrivals above the first (P-wave) multiple reflection. Across the terrace, four consistent S-wave reflection events were identified on three OBS stations (OBS 4, 2, 1) and used for a 2-D velocity model with some crossing raypaths (Fig. 13). However, station spacing was too large to define any meaningful lateral variations in S-wave velocity. At OBS stations 4 and 1, additional events were identified and included in one-dimensional velocity models at these locations (Fig. 14). Overall, we divided the first sediment layer (0–0.8 km depth below seafloor) across the terrace into an upper sublayer with S-wave velocity values of 0.38–0.4 km/s, and a lower sublayer with values of 0.46–0.6 km/s. The second sedimentary layer (>1 km depth below seafloor) showed velocity values from 0.8 to 1.0 km/s.

Combining the S-wave and corresponding P-wave velocities, we calculated Vp/Vs ratios for the upper sedimentary layers (Figs. 13 and 14). The shallowest sediments showed a Vp/Vs ratio of up to 4.2, and the ratio gradually decreased with depth to values of ∼3.0 for the layer corresponding to S-arrival S5 at ∼1 km depth below seafloor. At OBS station 6 on the abyssal plain, a very similar S-wave velocity structure was seen overall, but a jump to slightly higher S-wave velocities occurred within the shallowest sediment layer already at 0.3 km depth below seafloor. Otherwise, no significant difference existed for the S-wave velocity-depth structure from the abyssal plain to the terrace within these upper sedimentary layers.

The 2013 experiment confirmed and extended findings from previous studies conducted more than three decades earlier (Dehler and Clowes, 1988; Mackie et al., 1989). The first-arrival data that we were able to extract from the OBS records allowed rays to be traced up to 6 km below the seafloor, although ray density was sparse. We augmented the 2-D transect by incorporating arrivals from crossing lines and turns; overall, those arrivals supported a relatively robust P-wave velocity structure. Using waterborne multiple and seafloor peg-leg arrivals, we showed that the P-wave velocity structure is robust, because those arrivals are predicted within picking uncertainty, although they have quite different raypaths compared to the first arrivals. By combining observations from the studies in the 1980s (Dehler and Clowes, 1988; Mackie et al., 1989) with our new survey, two consistent features of the Queen Charlotte terrace are evident (Figs. 8 and 9): (1) a sharp sub-seafloor velocity discontinuity exists across the trough and into the terrace, and it has shallow (<2 km sub-seafloor depth) high-velocity material with P-wave velocities of ∼4 km/s, and (2) there is a narrow zone of reduced velocities around the Queen Charlotte fault.

The low ray coverage and steep change in topography from the offshore air-gun shots to the onshore stations created challenges in the construction of a broad-scale deeper velocity model. Our results were limited by the air-gun size used, OBS station spacing, and constraints on the relocation of OBS stations and time drift (outlined in Riedel et al., 2014b), as well as the low S/N ratio resulting from high-amplitude, low-frequency earthquake noise and previous shot noise.

Sediment compaction, porosity reduction, and associated dewatering over the short cross section of the terrace cannot easily produce seismic velocities as high as 4 km/s at <2 km sub-seafloor depth. A recent review of oceanic crustal velocities has shown that velocities of 3–4.5 km/s are typical of upper oceanic crust aged 0–10 Ma (Christeson et al., 2019). In the trough, the interpretation of these velocities as the top of igneous oceanic crust layer 2 is straightforward (Figs. 8 and 9). In the terrace, however, these high velocities do not fit well with the oft-reproduced cartoon of thick accreted sediments overlying an intact oceanic plate (e.g., Smith et al., 2003; Hyndman, 2015). The high velocities suggest a different kind of structural framework, one in which igneous crustal material may have been incorporated in the shallow terrace by compressional to transpressional deformation.

The uppermost sediments (<1 km depth) of the terrace (∼15–20 km landward of the trough) are lower in velocity than sediments of the Cascadia accretionary prism at a distance of 15–20 km landward of the trench (Fig. 14; Riedel et al., 2020). We interpret this as evidence that sediments of the Queen Charlotte terrace have experienced less compression and therefore less dewatering than those in the accretionary wedge. Gas hydrates, which increase P-wave velocity of the sediments, only have an effect within the upper 200–300 m of sediment below seafloor, as seen in studies off northern Cascadia (Yuan et al., 1994; Dash and Spence, 2011; Yelisetti et al., 2014). A bottom-simulating reflector has been identified along the entire margin off Haida Gwaii (Riedel et al., 2020), but gas hydrate saturations must be lower than those off central northern Cascadia, because no significant increase in P-wave (or S-wave) velocity was detectable off Haida Gwaii with the currently available data.

We compared velocities deeper than the sediment cover (∼2 km thickness) at the central Queen Charlotte terrace to other known velocity-depth functions (Fig. 15), and our new study is similar to results from Dehler and Clowes (1988). Both these functions are similar to velocities at the Cascadia accretionary prism but higher than those in the Winona basin (∼350 km farther south; Davis and Clowes, 1986) and the reference profile for normally compacted sediments proposed by Han et al. (2017). Our data from the Queen Charlotte terrace between 2 and 4 km depth, and the earlier results between 2 and 6 km depth by Dehler and Clowes (1988), respectively, represent the critical depth range for which an alternate interpretation to a compacted sediment lithology may be required.

It seems unlikely that a process that would preferentially dewater only sediments deeper than 2 km below seafloor would shift velocities to higher values than those proposed as normally compacted sediments. We, therefore, think it may be possible that these high velocities represent fractured and uplifted oceanic crust. We interpret these results to indicate that convergence between the Pacific and North American plates has faulted Pacific oceanic crust and uplifted it into the terrace, as previously hypothesized by Dehler and Clowes (1988). The interpretation of a simple subducted slab has been driven in part by large amounts of convergence predicted by plate circuit calculations of relative plate motion (DeMets and Merkouriev, 2016). Evidence for underthrusting and the presence of a (low-velocity) oceanic plate beneath Haida Gwaii from receiver function analyses (Smith et al., 2003; Bustin et al., 2007) does not contradict our findings.

A different approach to estimates of relative plate motions using geomorphology of the Queen Charlotte fault concluded that much less convergence has occurred in this area: 6 mm/yr (Brothers et al., 2020) compared to 16 mm/yr (DeMets and Merkouriev, 2016) for regional relative plate motion between the Pacific and North American plates. This reduced amount of convergence may be accommodated by deformation of the oceanic plate in the terrace and may not require underthrusting of the Pacific plate (Brothers et al., 2020). Transpression provides a mechanism for fracturing, buckling, uplifting, and accreting crustal material underneath the Queen Charlotte terrace, as also suggested by Schoettle-Greene et al. (2020). Our refraction results can be interpreted to substantiate this interpretation, although we note that along with deformation of the oceanic plate in the terrace, some part of the plate may have underthrust North America as interpreted from receiver function analyses (Smith et al., 2003; Bustin et al., 2007; Gosselin et al., 2015).

Around the Queen Charlotte fault, reduced velocities were found in the layer that had velocities >4 km/s in surrounding material. This was also found by Dehler and Clowes (1988) and interpreted to arise from sheared sedimentary and basaltic rocks. Further north and at a deeper level, reduced velocities immediately west of the Queen Charlotte fault were attributed to faulting and/or alteration of the Pacific plate (Walton et al., 2019).

New data in this region are needed to address the question of anomalous higher velocities within the Queen Charlotte terrace and lower velocities around the Queen Charlotte fault. A combination of densely spaced OBS refraction and MCS reflection data with large enough sources to image down to depths of >4 km below seafloor would also have the lateral resolution necessary to resolve these structures. The steep topography and associated changes in the velocity structure at the seaward limit of the Queen Charlotte terrace likely also require a dense shot spacing for the MCS and OBS data acquisition and special data processing (e.g., prestack depth migration).

Active-source wide-angle refraction data acquired in 2013, 4 months after the Haida Gwaii Mw 7.8 thrust earthquake, confirm the presence of high-velocity material (values ∼4.0 km/s) at shallow depth (2.0–3.0 km) underneath the Queen Charlotte terrace. The data confirm results from earlier velocity analyses based on the seismic data acquired across two parallel transects in 1979 at southern Moresby Island and 1983 at northern Moresby Island. Porosity reduction due to sediment compaction and associated dewatering over the cross section of the terrace cannot easily produce seismic velocities high enough to match our model results at such shallow depth. Upthrusting or imbrication of upper oceanic crustal material could explain the observations. Together with the observation of a heavily faulted incoming oceanic plate on a series of single-channel profiles on the abyssal plain west of the deformation front, we suggest a complex process of compressional to transpressional deformation resulting in disruption of the Pacific plate beneath the terrace.

The work conducted as part of the Haida Gwaii earthquake response would not have been possible without the support of many individuals. We would like to express special thanks to the Haida Nation, for permission to conduct the survey off Gwaii Hanaas, and to George Wesley, who joined the second expedition in January 2013, during which the active-source seismic refraction experiment was conducted. Additional thanks go to the Canadian Coast Guard (CCG) captains and crew of the two cruises onboard the CCGS John P. Tully, as well as staff from the CCG Regional Operation Centre in Victoria for their support in scheduling a vessel to conduct the research outside of the usual science patrols. Researchers critical to the success of the work were the two marine mammal observers, Rhonda Reidy and Jacklyn Barrs, who ensured that the survey was compliant with the Fisheries and Oceans Canada (DFO) seismic operation permit. We also want to thank those members of the Geological Survey of Canada who were instrumental in securing the substantial funding required for the offshore component of the earthquake response, including Adrienne Jones, Carmel Lowe, and Director General Daniel Lebel. We are also grateful to S. Han, who provided us with a digital version of the normal compaction curve published in 2017. This is NRCan contribution number (Numéro de contribution de RNCan) 20200295.

Global multiresolution topography (GMRT; Ryan et al., 2009) was used to generate topographic and bathymetric maps using ArcGIS V.10.2. Data are available for map-generation at: https://www.gmrt.org/GMRTMapTool/. The Canadian seismic reflection data used in this analysis are curated by Natural Resources Canada and are available online at the following site: http://ftp.maps.canada.ca/pub/nrcan_rncan/raster/marine_geoscience/Seismic_Reflection_Scanned/NRCan_Seismic_Reflection_Scanned.kmz. Additional data can be requested from scientists at the Geological Survey of Canada (curator Michelle Côté, michelle.cote@canada.ca) or through the Open Government Portal, https://open.canada.ca/data/en/dataset?organization=nrcan-rncan&q=geoscience.

1Supplemental Materials. Modeled ray paths and comparisons to picked arrivals are shown here for the refractions in sediment layers 1 (Fig. S1) and 2 (Fig. S2), the reflections from the bottoms of sediment layers 1 (Fig. S3) and 2 (Fig. S4), the refraction from the bottom of the crustal velocity layer (Fig. S5), and the head-wave from the bottom of the crustal-velocity layer (Fig. S6). Please visit https://doi.org/10.1130/GEOS.S.13286780 to access the supplemental material, and contact editing@geosociety.org with any questions.
Science Editor: Shanaka de Silva
Guest Associate Editor: Daniel S. Brothers
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.