Seismic-reflection data along the Haida Gwaii margin collected from 1967 to 2013 were used to identify gas hydrate–related bottom-simulating reflectors (BSRs). The BSRs occur along the Queen Charlotte Terrace only, within more strongly folded and tectonically deformed sedimentary ridges. The BSRs are absent within well-bedded and sediment-filled minibasins. The BSR is modeled as the base of the phase boundary of the methane hydrate (structure I) stability zone and is used to estimate geothermal gradients. The P-wave velocity structure required to convert observed depths of the BSR in two-way time to meters below seafloor was constrained from ocean-bottom seismometers. The BSR-derived gradients are lower than data from heat-probe deployments in the region, as well as predicted values from previous modeling of the large-scale tectonic thermal regime. Lower values of the BSR-derived thermal gradients may be due to topographic effects across the ridges where BSRs were observed. The previously identified landward decrease in thermal gradients across the terrace was also identified to a lesser extent from the BSRs, in accordance with the effects of oblique convergence of the Pacific plate with the North American plate. Geothermal gradients decreased from south to north by a factor of two, which is likely an effect of plate cooling due to an increase in age of the underlying plate (ca. 8 Ma off southern Haida Gwaii to ca. 12 Ma at Dixon Entrance) as well as the fact that sediments triple in thickness over the same distance. This may be due to downward flexure of the underlying crust during transpression and/or a high flux of sediments through Dixon Entrance.
The Queen Charlotte margin is primarily a dextral transcurrent continent–ocean boundary between the North America and Pacific plates. Along the southern portion of this margin, in the region of Haida Gwaii (formerly named Queen Charlotte Islands), relative motion of the two plates has had a convergent component up to 20° at 15–20 mm/yr over the past ∼6–10 million years (Rohr and Dietrich, 1992; Hyndman and Hamilton, 1993; Atwater and Stock, 1998; Rohr et al., 2000; DeMets and Merkouriev, 2016). The Queen Charlotte fault (QCF) and its northern extension, the Fairweather fault, are regions of repeated large strike-slip earthquakes (e.g., Rogers, 1982; Ristau et al., 2007) of magnitude Mw 8 and above. The margin has seen a new wave of attention after the 27 October 2012, Mw 7.8 thrust earthquake (James et al., 2013), which may be related to oblique subduction 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.
Studies along the margin have been conducted since the early 1960s due to the recognized earthquake potential in the region (e.g., Chandra, 1974; Rogers, 1982; Nishenko and Jacob, 1990). Seismic-reflection surveys imaged a striking morphology of the slope (e.g., Chase and Tiffin, 1972; Davis and Seemann, 1981) as shown in Figure 1: a bench-like feature composed of folded sediments ∼25 km wide was named the Queen Charlotte Terrace (QCT). Landward, the QCT was bound by a steep slope rising to a continental shelf only ∼1 km wide and seaward by a steep slope that descended to a sediment-filled trough. This morphology has been linked to that of subduction zones (Hyndman and Hamilton, 1993); yet transpressive deformation alone can also load the Pacific plate in the QCT and generate a similar morphology (Rohr et al., 2000). Refraction models across the QCT to date have not imaged a simple slab; several structural models explaining the complex tectonics have been proposed (Horn et al., 1984; Dehler and Clowes, 1988; Rohr et al., 2000; Tréhu et al., 2015; Hyndman, 2015).
Thermal modeling to further constrain an oblique subduction process (Smith et al., 2003; Wang et al., 2015) was based on sparse heat-flow data from measurements with a heat probe off the southern portion off Haida Gwaii and in Hecate Strait (Hyndman et al., 1982) as well as onshore and offshore borehole temperature measurements (Lewis et al., 1991). To augment the available coverage of heat flow and thermal gradient data, we re-investigated seismic data across the margin (Fig. 2) to identify gas hydrate–related bottom-simulating reflectors (BSRs). The BSR is used as a seismic expression of the base of the structure I (sI) gas hydrate stability field. Here, we only consider a pure methane-seawater system. The gas hydrate stability or phase boundary is primarily governed by pressure (depth), temperature, gas composition, and pore-water salinity (Sloan, 1998). As such, the BSR is in principle an isotherm, and by converting the seismically observed reflection time of the BSRs to depth, thermal gradients can be determined. We adopt this technique from previous studies conducted off Vancouver Island (Ganguly et al., 2000; Riedel et al., 2010b) but modify required parameters (e.g., P-wave velocity function) to the conditions off Haida Gwaii. We utilize single-channel seismic (SCS) and multichannel seismic (MCS) reflection data to image BSRs (digital and vintage analog data), as well as ocean-bottom seismometer (OBS) data to constrain the velocity information. An OBS survey was conducted in 2013 to capture the aftershock sequence of the 27 October 2012, Mw 7.8 thrust earthquake (Riedel et al., 2013a, 2013b; Riedel et al., 2014); it allows us to improve the velocity function beyond that available from MCS-based stacking velocities (Scheidhauer et al., 1999).
DATA AND METHODS
For this study, we selected seismic-reflection data from the region off Haida Gwaii and farther north along Dixon Entrance (Fig. 2; Table 1). Vintage data off Haida Gwaii from surveys conducted by the Geological Survey of Canada (GSC) are available as original paper records and were digitized and converted to SEGY standard format. Absolute depth in two-way time (TWT) is determined from information available on the paper records; amplitude and phase are non-recoverable. Thus, these data records can only be used as images without any further processing applied (e.g., band-pass filtering). Digital data from the GSC database (Table 1) are SCS reflection data only; these data were processed with band-pass filters to increase the signal-to-noise ratio. Digital SEGY data for cruises F-7-89 (Bruns et al., 1992), L-5-77 (Snavely et al., 1981), L-3-78 (Bruns et al., 1987; Bruns and Carlson, 1987), and EW9412 (Hollister and Andronicos, 1997; Scheidhauer et al., 1999; Rohr et al., 2000; Tréhu et al., 2015; Walton et al., 2015, 2018) were downloaded “as is” from the United States Geological Survey (USGS) or University of Texas Institute of Geophysics (UTIG) databases. No additional processing was applied to these data sets for this study. Additional vintage SCS profiles were collected across the Fairweather fault (USGS survey S-5-78; Carlson et al., 1979, 1988) and QCF (surveys RC1109, 1967 and RC1407, 1971, available in image format only), but image quality is poor.
Where available, navigation for digital seismic data was used as provided through the various database sources. Vintage analog data have lower navigation accuracy, especially due to less accurate fixes to determine the position of the vessel at the time of acquisition. Line location was thus matched where possible with structural seafloor features (e.g., the eastern edge of the abyssal plain or prominent ridges) and otherwise linearly interpolated between start and end points of line locations.
From these various types of seismic data, we identified BSRs as indicators for the base of gas hydrate sI stability assuming a purely methane-dominated system. Examples of BSRs are given in Figure 3. The BSR must fulfil the following criteria: (1) mimicking seafloor topography; (2) crosscutting regular stratigraphy; and (3) indicating a polarity opposite to that of the seafloor (for original digital records).
The BSR may also be seen to be associated with amplitude brightening underneath the phase boundary (indication of free gas) or be associated with amplitude truncations at the phase boundary without forming a separate reflection.
Geothermal gradients (and associated heat flow) were determined from BSRs using the identified reflection times of seafloor and BSR following previous studies (e.g., Ganguly et al., 2000; Riedel et al., 2010b). Here, we adopt the following equations to define geothermal gradients:
- (a) Seafloor temperatures (Tsf) are assumed to follow a water depth (wd) dependent trend (adopted from Riedel et al., 2018):
- (c) Temperature at the BSR (Tbsr) is defined using the following polynomial equation, which is an approximation (R2 value of 0.99987) to the phase boundary of a pure seawater-methane system as defined in Sloan (1998):with P as pressure given by Equation 3, and the coefficients a1 to a5 defined as: a1 = −0.00036543, a2 = 0.02104761, a3 = −0.46787289, a4 = 5.37421234, a5 = −12.27751347. This equation was adopted from Riedel and Collett (2017).
- (d) Thermal conductivity (k) used to convert thermal gradients into heat flow values is estimated by Davis et al. (1990) to follow this depth-dependent trend:
P-wave velocity is an important constraint in obtaining thermal gradients from BSR depth variations. Along the Haida Gwaii margin, only sparse velocity information for the upper few hundred meters below seafloor (mbsf) exists. Previous MCS data acquisition is limited to a few lines from the survey L-5-77 (Rohr et al., 1993) and the survey EW9412 (e.g., Scheidhauer et al., 1999). A previous experiment using three OBS stations resulted in constant velocity values for the uppermost sediment section up to 1500 m below seafloor (mbsf) (Hyndman and Ellis, 1981). Velocities from OBS records using refracted arrivals also provide only rough estimates of the average velocity within the upper ∼1000 mbsf. The first detectable refracted arrivals from the 2013 expedition along OBS stations 1–5 show average values of 2000–2200 m/s for the upper 1000 mbsf (e.g., Riedel et al., 2014), similar to values reported by the refraction study of Dehler and Clowes (1988).
Therefore, we utilized reflections up to a depth of 800 mbsf recorded on the OBS stations from the 2013 survey (Riedel et al., 2013a) to define 1D velocity-depth profiles beneath each station. During the 2013 survey, a single 520 in3 (∼8 L) G-gun shot at a 30 second interval (equivalent to a ∼60 m shot spacing) was used along two main profiles across six OBS stations (Riedel et al., 2013a) near the epicenter of the 27 October 2012, Mw 7.8 thrust earthquake. The air gun was fired with a pressure of 1800 psi. The sample rate of the OBS recordings was 4 ms (equivalent to a Nyquist frequency of 125 Hz). The OBS instruments (equipped with a three component short-period geophone and one hydrophone) were provided from the instrument pool maintained by Dalhousie University and the GSC. Details on OBS data pre-processing (clock-drift definition, instrument re-localization) are described in Riedel et al. (2014).
First arrivals of up to eight reflections were determined on the OBS data records and forward modeled based on the ray-tracing code defined by Zelt and Smith (1992) to obtain best-fit models constraining both depth and velocity in the subsurface. Picking uncertainty of reflections on OBS records was set to 10 ms for OBS stations 2, 3, and 6. A picking uncertainty of 15 ms was set for stations OBS 1, 4, and 5 because these stations were recording data with a lower sampling rate and had an overall lower signal to noise ratio. P-wave velocity of each layer is constrained by the overall shape of the reflection hyperbola and can be determined within ∼50 m/s using this technique. With the velocity defined, reflection time of each layer beneath the OBS station can be converted to depth. The uncertainty in this depth calculation is a direct function of the velocity uncertainty and did not exceed 10 m in the analyses.
RESULTS OF OBS VELOCITY MODELING
Using the reflections observed on OBS stations, forward modeling was used to match picked arrivals with synthetic arrival times within the picking uncertainty of the arrivals. An example of this analysis is provided for OBS station 3 along the N-S crossing of the OBS stations 3–5 (Fig. 4). Station 3 allowed identification of eight reflections for which velocities were defined. They increase from near-seafloor velocities of 1500 m/s to 2280 m/s at 740 mbsf (Fig. 5). The structural model was defined using the SCS data as guidance. The acquisition geometry and subsurface velocity field allow construction of a model for ∼2 km to both sides of the OBS station only. With an OBS spacing of 5 km along the seismic line (Fig. 4), no crossing rays can be modeled. However, between OBS stations 3 and 4, the same upper four reflections were identified and modeled. Farther south along the seismic line, a BSR is observed beneath OBS station 5. The BSR is associated with a velocity reduction from ∼1700 m/s to 1500 m/s for the layer beneath the BSR. This is the only OBS station with a direct indication of a BSR and associated velocity reduction.
GEOTHERMAL REGIME FROM BSR DEPTH
BSRs were observed on a number of seismic lines and converted to thermal gradients (Fig. 7) with values lower than 100 °C/km. They range from as low as 28.9 °C/km to as high as 92.5 °C/km, and the median value is 44 °C/km. Combined uncertainties in the conversion of BSR depths to thermal gradients from P-wave velocity and seafloor temperature are estimated to be between 5 and 7 °C/km. Thermal gradients from heat-probe measurements (Hyndman et al., 1982) in comparison to the BSR-derived values are up to 30 °C/km higher for comparable geographic regions (see Fig. 7, region off Moresby Island). Reasons for this discrepancy are difficult to assess because the original data from the heat-probe measurements are no longer available to verify calculations. However, near-seafloor measurements of temperature can be affected by seasonal temperature variations, which can alter the determination of reliable thermal gradients. In contrast, the BSR, taken as the seismic expression of the lower limit of gas hydrates, represents the thermal state of the deeper sediment column where such near-surface effects are not an issue. Seismic data inherently result in a smooth representation of the thermal state when using the BSR as a proxy for heat flow, due to the low-frequency sources used for imaging (<120 Hz in the cases considered here). The BSR-derived values are comparable to each other, and thus, regional trends or local anomalies can be described and used for further interpretations, even if the total absolute value may be somewhat uncertain.
Overall, a northward decrease in BSR-derived thermal gradients is evident: values range from ∼70 to ∼90 °C/km off the southern portion of Haida Gwaii and from ∼30 to ∼50 °C/km off the northern portion. There is also a small decrease in the thermal gradient of ∼20 °C/km over a horizontal distance of 20 km across the terrace in an east-west trend (Figs. 8A and 8B) according to the expected reduction in thermal gradients as a result of deepening of the Pacific plate toward North America, as modeled by Wang et al. (2015). There may also be a change in thermal gradient at the QCF trace, where some data may suggest a local increase in gradients by up to 20 °C/km (on average 70–80 °C/km) relative to the expected reduction from oblique subduction (Fig. 8B). Uncertainties in the conversion of BSR depth to heat flow are roughly within the same order of magnitude, and thus, no conclusion can be drawn from this observation at this time. Additional seismic data acquisition to better image the BSR and determine the local velocity variation, as well as detailed heat-probe measurements, are required to confirm this local trend in the thermal structure.
The overall northward decrease (especially north of 52°53′N, or ∼115 km along the trough-parallel transect) in BSR-derived thermal gradients is depicted in Figure 8C. All data from BSR and heat-probe measurements are shown, irrespective of the individual distance to the trough, and data are therefore scattered. Furthermore, the data were not corrected for any local topographic effects that may reduce (ridges) or increase (troughs) thermal gradients and surface heat flow, which add to the scatter. The original heat-probe data have not been corrected for topographic effects.
Several assumptions were included in the conversion of BSR depth to thermal gradients and should be evaluated before any further interpretation. We first evaluate the impact of a varying P-wave velocity profile and changes in bottom-water temperatures. The same observed time differences in BSR and seafloor depths are converted using different P-wave velocity-depth functions defined by Equations (6A) and (7B). This results in a maximum of 5 °C/km difference in thermal gradients from using Equation (7B) after Scheidhauer et al. (1999) with a lower velocity gradient than Equation (6A) based on the OBS-determined velocities from the Queen Charlotte Terrace.
Additional uncertainties in heat-flow estimation from BSR depths arise from unknown thermal conductivity values. An equation describing thermal conductivity as a function of depth was proposed by Davis et al. (1990), although reducing this to a constant value of 1.0 W/(K·m) results in a reduction in heat-flow values by ∼8 mW/m2 but leaves the value of thermal gradient unchanged. Seafloor temperature may vary up to 2 °C and thus change the thermal gradient accordingly. These two effects result in a maximum uncertainty in the derived thermal gradients of 10 °C/km, much smaller than the observed regional thermal gradient trend.
The BSR has been taken as expression of the base of pure methane gas hydrate stability zone (BGHSZ) in this study. If higher hydrocarbons than methane are present (i.e., ethane and propane), the BSR could be substantially deeper because structure II (sII) gas hydrates are potentially present. Thus, if a BSR (and BGHSZ) from a mix of gases with higher hydrocarbons is falsely interpreted as being from pure methane only (sI hydrate), the calculations for thermal gradients result in lower values than what should be present. To explain the magnitude of the regional trend in thermal gradient by a difference in hydrocarbon gas mixtures (i.e., reducing the thermal gradient by an approximated factor of two over a distance of ∼200 km from south to north) requires unreasonably high concentrations of higher hydrocarbons (e.g., >20% ethane) off northern Haida Gwaii. There is currently no evidence for such a strong chemical gradient along the margin. None of the thermal gradient data obtained from seismic data or from heat probes have been corrected for any topographic effects. Thus, there may be local discrepancies and reduced gradients from a topographic effect that is de-focusing the (conductive) heat transport across the frontal ridges or focusing heat in troughs resulting in local apparently higher gradients. BSRs are generally only observed across ridges along the Haida Gwaii margin and not over flat topography and/or minibasins within the QCT. Seafloor topographic changes of ∼250 m associated with thrust ridges seen at the northern Cascadia accretionary prism result in a reduction of surface heat flow up to 20 mW/m2 (e.g., Riedel et al., 2010b). BSR-derived values reported in this study may therefore be underestimated, but because all data may suffer equally from this topographic effect, they remain comparable to each other.
P-wave velocities determined from OBS data across the QCT are significantly lower than those determined at the central northern Cascadia margin. As shown by Ocean Drilling Program (ODP) and Integrated Ocean Drilling Program (IODP) drilling and accompanying studies (see summary by Riedel et al., 2010a), the velocity increase seen in central northern Cascadia around drill sites U1327 and 889 (and included in Fig. 6) is associated with a lithological change from the upper ∼100 mbsf of stratified basin-fill sediments (with higher porosity and lower velocities) and deeper accreted sediments (lower porosity and higher velocity) that show almost no coherent seismic reflectivity. In addition to this general depth-trend in velocity are isolated spikes of high velocity attributed to gas hydrates. At the Haida Gwaii margin, no similar velocity increase in accreted sediments is observed. The OBS data across the QCT were acquired at a setting that is at a landward distance of 15–20 km from the trench (deformation front) and thus at a similar distance to the deformation front as at northern Cascadia margin, where the bulk of the velocity data have been acquired in the past.
The seismic-reflection data at the QCT show continuous seismic reflections, despite structural deformation and folding (see e.g., Figs. 3C and 4), in the area where OBS were deployed. Thus, the amount of structural deformation (and accompanying reduction in porosity and velocity increase) at the Haida Gwaii margin may not have reached the same amount as at the central northern Cascadia accretionary prism. It should be noted that convergence at Haida Gwaii is highly oblique, not margin-perpendicular as occurring at central Cascadia. Thus, a transect across the QCT perpendicular to the trench would not have experienced a similar amount of compression and deformation (with porosity reduction) over a distance of ∼20 km as off central Cascadia.
The lower amount of deformation from oblique convergence may be linked to lower volumes of fluid expulsion and methane advection. Thus, the generally lower velocities may also be a result of overall lower gas hydrate accumulations and the lack of widespread BSR observations (limited to folds and ridges) across the entire region off Haida Gwaii. They may also be a consequence of the lack of free gas from upward migrating fluid (as postulated by Hyndman and Davis, 1992, and further developed by Riedel et al., 2010a, for northern Cascadia).
It should be noted that this model is for ocean floor away from any exposed bathymetric highs and does not consider hydrothermal circulation. Crustal ages off the southern tip of Haida Gwaii are ca. 8 Ma (Fig. 7), which is equivalent to a surface heat flow of ∼170 mW/m2, using the simple Equation (8). At the northern end of the Haida Gwaii margin, crustal ages are ca. 12 Ma (Fig. 7), and heat flow is thus ∼146 mW/m2. Using this age-related cooling, a decrease in heat flow by ∼15% along a margin-parallel transect is predicted, which is insufficient to explain the observed trend. In the profile shown in Figure 8C, we adopted the age-related cooling trend assuming a starting condition of 80 °C/km at the south end of the profile to match our observations. The change in BSR-derived thermal gradients is especially noticeable north of 52°53′N, which is close to the proposed possible change in tectonic environment from oblique subduction to strike-slip, defined as 53.1°N by ten Brink et al. (2018) and Tréhu et al. (2015). Walton et al. (2015) also argued for no Pacific plate subduction north of 53.2°N based on seismic-reflection data.
One additional effect that may reduce the thermal gradients is related to the overall sediment thickness. Several previous studies have determined that sediment thickness on the incoming Pacific plate off the southern terrace ranges from 1.0 to 1.5 km (e.g., Hyndman and Ellis, 1981; Dehler and Clowes, 1988; Mackie et al., 1989; Spence and Long, 1995). We have used the available seismic data showing a clear reflection from the top of subducting crust in the trough (particularly time-migrated data from the EW9412 and L-5-77 surveys) to estimate sediment thickness at the foot of the terrace along the entire margin. The sediment thickness values (marked in white font in Fig. 7) were estimated using an average velocity of 1.9 km/s. Sediment thickness is about twice as large off the northern end of the Haida Gwaii margin (3.2–3.3 km) compared to off southern Moresby Island (1.1–2.0 km). This increase in sediment thickness has been explained by Walton et al. (2015) as the result of tectonic flexure of the Pacific plate when it was farther south. We also note the proximity of this crust to Dixon Entrance and Hecate Strait, which are significant sediment fairways of the Coast Mountains and Queen Charlotte margin. High sediment input to the margin from repeated glacial processes (Barrie and Conway, 1999) could have provided a source for higher sediment thickness off Hecate Strait and associated relatively thinner sediments farther south. Thus, the reduced thermal gradients farther north could be explained by a combination of age-related cooling of the underlying crust and higher sedimentation rates and inflated sediment thickness.
SUMMARY AND CONCLUSION
The extensive database of seismic-reflection surveys off Haida Gwaii allows for the estimate of thermal gradients from the observation of bottom-simulating reflectors, used as an indicator of the base of the sI gas hydrate phase boundary. Thermal gradients calculated from heat flow are useful constraints in modeling tectonic structure of the margin. We have determined thermal gradients along 26 seismic lines, which significantly augment the available data on the thermal state of the margin, particularly on the northern end of the margin. Surveys span from 1967 to 2013. Many reflection lines are only available as paper copies that were digitized and geo-referenced by the GSC.
Conversion of reflection times of seafloor and BSRs to depth and then computation of the thermal gradient required careful velocity control. We have utilized the available multi-channel seismic-reflection data and newly acquired digital ocean-bottom seismometer data to constrain the P-wave velocity across the Queen Charlotte Terrace for the uppermost few hundred meters. In comparison to sediments in the Cascadia accretionary prism, these sediments are lower in velocity and by inference have experienced less compression and dewatering.
Our data suggest a reduction in thermal gradients across the Queen Charlotte Terrace from west to east. This has been interpreted in previous studies as a result of the oblique subduction or deepening of the Pacific plate under the terrace. In addition, thermal gradients are reduced by a factor of two from south to north. We interpret this to be the result of increasing age of the underlying Pacific plate from ca. 7 Ma to ca. 12 Ma, as well as increased sediment thickness above the oceanic crust. Sediment thickness at the northern end of Haida Gwaii at the Dixon Entrance (∼3.2 km) is almost twice that off Moresby Island (1.1–2.0 km) and likely the result of sediments that have been carried through Dixon Entrance.
Data from offshore Haida Gwaii and Alaska are available through the U.S. Geological Survey seismic database and can be downloaded from https://walrus.wr.usgs.gov/namss/search/. Data from the UTIG database can be accessed through this portal: http://www-udc.ig.utexas.edu/sdc/map_asp.htm (Scheidhauer et al., 2015). 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. Heat-probe data were previously published, and the original data including navigation can be accessed in the paper by Hyndman et al. (1982).
The authors would like to acknowledge all people who participated in the collection of data around Haida Gwaii that are combined in this study. Additional thanks go to Roger McLeod and Robert Kung from the Geological Survey of Canada (Sidney) for their ArcGIS support.