Focused fluid flow and methane venting along the Queen Charlotte fault, offshore Alaska (USA) and British Columbia (Canada)

Fluids play an important role along plate boundaries because they influence almost every chemical, physical, mechanical, and thermal process that occurs in the lithosphere (Judd and Hovland, 2007; Kastner et al., 1991; Saffer and Tobin, 2011; Tryon et al., 2001). Such processes are all directly impacted by the transport of mass, heat, and chemical species in a hydrogeologic framework (Tryon et al., 2001). In turn, these processes can influence fluid and gas budgets and geochemical cycling (Kopf, 2002) as well as support chemosynthetic communities on the seafloor (e.g., Levin et al., 2016). Fluid seepage by cold seeps, where hydrogen sulfide, methane, and other hydrocarbon-rich fluids are released into the sediment and water column, plays an important role in transferring methane carbon from storage in ocean-floor sediments into the ocean and atmosphere. For example, Judd (2004) estimated that global methane fluxes from the seafloor to the ocean and atmosphere vary between 0.4 and 48 Tg/yr. In contrast to hydrothermal vents, temperatures of cold seeps (referred to herein as seeps) are similar to surrounding seawater temperatures. Some mechanisms of fluid seepage at seeps include fluid flow driven by excess pore pressure in areas of high sedimentation (Berndt, 2005), seismic activity (e.g., Fischer et al., 2013), links to sea-level lowstands (Feng et al., 2010; Liebetrau et al., 2010; Teichert et al., 2003), hydrological and tidal pumping, and warming of bottom water (Suess, 2014). These temporally recent processes are superimposed on geologic controls of fluid migration, whereby the geologic setting defines seep characteristics and the driving mechanisms of fluid expulsion (Judd and Hovland, 2007; Suess, 2014). For example, rapid sediment loading along passive margins promotes seep development due to pore fluid overpressure (Dugan and Flemings, 2002), and at active convergent margins, plate motions facilitate faulting and deformation, leading to sediment compaction, the formation of fluid pathways, and traps for deeply sourced fluids that control seepage patterns (Kluesner et al., 2013; Sahling et al., 2008; Watson et al., 2019). Seeps documented along strike-slip fault systems are most commonly observed along zones of transpression and folding (Dupré et al., 2015; Fossen et al., 2010; Kluesner and Brothers, 2016; Paull et al., 2008). For example, recent investigations along the strike-slip margin in southern and central California (USA) suggest that localized areas of transpression along transform faults may generate tectonic features such as fault bends and stepovers, which appear to act as loci of substrate fluid flow (Kluesner and Brothers, 2016) and in turn support chemosynthetic habitats and authigenic carbonate formation (Conrad et al., 2018; Maloney et al., 2015; Stakes et al., 1999). Similar patterns have been observed in the Sea of Marmara (Turkey), where seeps are preferentially clustered along zones of localized uplift, folding, and transpressional bending along the North Anatolian fault (Dupré et al., 2015), as well as along the Mendocino transform fault (offshore California; Stakes et al., 2002) and Guaymas Basin (offshore Mexico; Hensen et al., 2019; Paull et al., 2007), where active seafloor fluid seepage along transform margins may represent important conduits for facilitating vertical fluid flow.

Authigenic carbonates commonly form at seeps as a result of anaerobic oxidation of methane (AOM) via sulfate reduction (CH₄ + SO₄²⁻ → HCO₃⁻ + HS⁻ + H₂O). The production of bicarbonate (HCO₃⁻) increases the pore-fluid alkalinity and favors carbonate precipitation. These carbonates represent a proxy record of the local and regional controls on the source and flux of methane carbon, the conditions and environment under which carbonates formed, and information regarding fluid-sediment and fluid-rock interactions (Campbell, 2006; Formolo et al., 2004; Magalhães et al., 2012; Naehr et al., 2007). For example, differences in the isotopic composition of carbonate cements along the Cascadia margin (northeastern Pacific) highlight the contrasting hydrogeologic settings between shallow-sourced fluids in a thrust-fault setting versus a deep origin of fluids from a strike-slip fault associated with an accretionary wedge (Sample and Reid, 1998). Integrating geochemistry, mineralogy, and petrology of authigenic carbonates with multichannel seismic reflection (MCS) profiles as well as multibeam bathymetry, backscatter, and water-column imaging data offers a valuable opportunity to investigate the relationship between fluid flow, sedimentation patterns, and tectonics along a proglacial transform margin. To evaluate potential flow pathways and sources of fluid flow along the Queen Charlotte fault (QCF; offshore southeastern Alaska [USA] and western British Columbia [Canada]), a series of authigenic carbonates was collected and analyzed in a collaborative research program between the U.S. Geological Survey (USGS) and the Geological Survey of Canada. Data collected between 2009 and 2018 included geophysical surveys (multibeam bathymetric imaging, backscatter, multichannel seismic reflection profiling, and 3.5 kHz CHIRP subbottom profiling), piston coring, seafloor photography, grab sampling, and conductivity-temperature-depth (CTD) casts along the continental margin of western British Columbia and southeastern Alaska (Barrie et al., 2013, 2018; Brothers et al., 2019; Greene et al., 2019). This paper presents the first comprehensive examination of seabed fluid-flow indicators along the QCF margin, including gas bubbles observed in water-column backscatter data, the presence of chemosynthetic communities such as bacterial mats and bivalves, precipitates of authigenic carbonate rocks and pavement, and mud volcanoes (Barrie et al., 2020). We also evaluate structural and stratigraphic controls on fluid-flow pathways using seismic reflection profiles. Results from this study therefore expand our understanding of the possible origin and occurrence of fluid seepage along an ocean-continent transform boundary and test the assumption that seeps along strike-slip fault systems are focused along faults and zones of transpressional folding in a poorly studied geographic region.

The offshore QCF and its northern onshore extension, the Fairweather fault, form an ∼1150-km-long transform boundary that separates the Pacific and North American tectonic plates (Fig. 1). Early studies noted evidence for progressive right-lateral offset on the QCF at the southern edge of the Yakobi Sea Valley (Fig. 1) and mapped its general location along the ocean-continent boundary for another ∼800 km southward to the Queen Charlotte triple junction (Bruns and Carlson, 1987; Hyndman and Hamilton, 1993; Rohr and Furlong, 1996; Von Huene et al., 1979). While the QCF carries ∼50 mm/yr of right-lateral transform motion, portions of the plate boundary are also transpressive (Brothers et al., 2020; Hyndman and Hamilton, 1993; Rohr et al., 2000; Tréhu et al., 2004). Convergence rates are expected to increase gradually southward from the Yakobi Sea Valley, where the QCF is closely aligned with plate motion, to a maximum near the southern end of the fault and west of Haida Gwaii (British Columbia; Fig. 1; Brothers et al., 2019, 2020; Hyndman, 2015; Rohr et al., 2000; Tréhu et al., 2015). The oblique convergence expected to the west of Haida Gwaii has led to the buildup of the Queen Charlotte terrace, a ∼25-km-wide, ∼150-km-long topographic platform located between the QCF on the east and the abyssal ocean basin to the west (Fig. 1). Within this region, local transpressional bends along the fault may generate even higher rates of convergence (as much as 15–20 mm/yr), which is evident in the tectonic geomorphology observed along the southern QCF (Brothers et al., 2020). The formation of the terrace has been described by competing end-member models: one attributes the terrace to an accreted sediment wedge from oblique subduction of the Pacific plate beneath Haida Gwaii (e.g., Hyndman, 2015; Hyndman and Hamilton, 1993), the other to transpressional deformation and crustal thickening of the eastern edge of the Pacific plate (Brothers et al., 2020; Rohr et al., 2000; Tréhu et al., 2015). Although these competing models are expected to generate distinct differences in crustal fluid flow and deformation patterns, the terrace contains relatively sparse geophysical observations. South of Haida Gwaii, plate-boundary deformation becomes distributed (e.g., Rohr, 2015; Rohr and Furlong, 1996), but the transform boundary ultimately links with Explorer Ridge, where nascent oceanic crust is translated to the northwest along the QCF (Rohr et al., 2000).

A series of recent earthquakes generated substantial interest in the QCF margin beginning with a M7.8 thrust event near Haida Gwaii in 2012 and followed by a M7.5 strike-slip event west of Craig (Alaska) in 2013 (Brothers et al., 2017; Lay et al., 2013; Yue et al., 2013). However, at the time of the 2012 event, the morphology of the QCF and surrounding continental slope within U.S. waters had not been comprehensively surveyed using high-resolution multibeam sonar or modern seismic reflection profiling techniques, and only a narrow strip of seafloor along the QCF had been mapped this way in Canadian waters (Barrie et al., 2013). Therefore, detailed information about the physical properties of the offshore margin was sparse, and little evidence for fluid seepage had been directly documented. Furthermore, recently documented submarine landslides along the margin attest to the potential for earthquake triggering of slope failures (Brothers et al., 2019; Greene et al., 2019). Taken together, high rates of sediment accumulation, fluid overpressure from sediment compaction, and tectonic deformation can generate a variety of sources for focused fluid flow and seafloor seeps. An active seafloor seep was discovered in 2011 offshore Haida Gwaii (Barrie et al., 2013) and another in 2015 along the crest of a ridge located near the international border (Greene et al., 2019). Greene et al. (2019) proposed that fluid flow is prominent along most of the central and southern QCF and that fluid flow indicators occur close to, or within, the fault zone, which in turn may be acting as a conduit for the venting of gases and fluids. In addition, the gas could be accumulating in antiforms that form many of the conjugate ridges described by Tréhu et al. (2015) along the continental slope. Gas plumes could originate from multiple sources, such as from the compaction and degradation of organic-rich sediment, continuously forming hydrocarbon reservoirs (Tinivella and Giustiniani, 2012), and dissolution of clathrates or frozen gas-hydrate systems (Depreiter et al., 2005; Dickens and Quinby-Hunt, 1994; Dimitrov, 2002; Sauter et al., 2006; Tinivella and Giustiniani, 2012).

Seafloor mapping and morphological analyses presented in this study are based on several multibeam echosounder (MBES) data sets. The QCF and surrounding slope were mapped between 2009 and 2018 using high-resolution MBES systems, including a Reson Seabat 7160 (44 kHz) aboard the R/V Medeia, a Reson Seabat 7111 (100 kHz) aboard the R/V Solstice, and a Kongsberg Simrad EM710 (70 kHz) aboard the U.S. National Oceanic and Atmospheric Administration (NOAA) Ship Fairweather and the Canadian Coast Guard Ship (CCGS) Vector (Balster-Gee et al., 2017a; Barrie et al., 2018; Brothers et al., 2019; Greene et al., 2019). These data were merged with lower-resolution multibeam bathymetry data sets and coastal relief models spanning the Gulf of Alaska region (available at https://www.ngdc.noaa.gov/mgg/bathymetry/relief.html). Bathymetry data were processed using the Teledyne CARIS HIPS software package and included offset corrections, tidal corrections, sound velocity correction, ping editing, and gridding. Raster surface derivatives (hillshade, gradient, and aspect) for the high-resolution data were created from 20 m bathymetric elevation models and used to interpret seabed morphology along the entire length of the QCF (Brothers et al., 2020). The MBES surveys along the continental slope in U.S. waters recorded water-column backscatter data (blue polygon on Fig. 1A), and bubble plumes were detected using the FM Midwater module within the QPS Fledermaus software package. Mapped bubble plumes were then combined with multibeam bathymetry and seismic reflection profiles for visualization and interpretation. MBES water-column data were not recorded during surveys in Canadian waters. However, single-beam 18 kHz profiles were collected using a Simrad EK60 onboard the R/V Ocean Starr in 2017 and the CCGS John P. Tully in 2011, 2015, and 2017, leading to discoveries of several bubble plumes in the water column. Although some seafloor ridges were suspected to be sites of seepage and so were specifically targeted for survey by EK60, we consider the southern half of the margin to be less thoroughly surveyed for seeps than areas to the north.

High-resolution multichannel seismic data were also collected in conjunction with bathymetry data onboard the R/V Medeia in the spring of 2016 (Balster-Gee et al., 2017a). The seismic shots were generated using a SIG 2Mille 50-tip mini-sparker operated at 1200 J. Reflected arrivals were recorded using a 32-channel solid-core, digital GeoEel streamer system. Additional MCS data were collected in 2016 onboard the R/V Norseman using an Applied Acoustic Delta Sparker Sled operated at 2400 J and a 64-channel GeoEel streamer system (Balster-Gee et al., 2017b). Recorded SEGD shots (SEGD file format is one of several standards developed by the Society of Exploration Geophysicists for storing geophysical data) were processed using the academic software packages SIOSEIS and Seismic Unix, as well the commercial Paradigm Echos software package. Seismic processing largely followed the approach described by Kluesner et al. (2018). After processing, seismic imagery was merged with multibeam bathymetry and mapped water-column gas plumes using the QPS Fledermaus software package.

Seismic attribute calculations can provide quantitative justification for the interpretation of a particular geologic feature. A powerful approach to detect geologic targets and isolate their acoustic response involves the combination of multiple attributes with structurally orientated calculations and discontinuity-attribute calculations (e.g., Gersztenkorn and Marfurt, 1999; Kluesner and Brothers, 2016; Marfurt et al., 1999; Tingdahl and de Groot, 2003; Tingdahl et al., 2001). This study applied a fully connected multilayer perceptron neural-network approach to two-dimensional high-resolution MCS profiles to optimize the detection of gas and fluid migration pathways in the vicinity of active seep sites (Brothers et al., 2014; Connolly, 2015; Heggland, 2005; Kluesner et al., 2013; Ligtenberg, 2005). Within the OpendTect software package, 36 attributes were used as input nodes into the supervised neural-network chimney calculation, and each node was weighted during the neural-network training (for details, see Table S1 in the Supplemental Material¹; Kluesner and Brothers, 2016; Brothers et al., 2014). The single hidden layer of the neural network consists of 15 nodes, and training was stopped when the normalized root-mean-square (RMS) error of the training set, as well as the misclassification percentage, reached a minimum value. Training was stopped before the RMS error trend shifted upward to prevent overfitting of the data and overtraining of the neural net. The generated meta-attribute represents a measurement of probability for the presence of a chimney structure, with 0 representing the lowest probability and 1 representing the highest. The chimney meta-attribute results were projected onto seismic cross-sections using a gradational color scale that reveals only the highest (∼80% and above) chimney probabilities. The chimney analysis was used to identify two features: zones of probable fluid migration (e.g., vertical pipes, faults, and chimneys) and gassy sediments. The frequency-dependent reflectivity of gassy sediments due to scattering and absorption has been previously documented (e.g., Wood et al., 2008). The chimney meta-attribute analysis described above therefore represents a powerful method to discern subtle patterns of gas-related attenuation and their spatial relationships with the surrounding structure and stratigraphy.

Seep sites (seeps 1–18; Fig. 1A) identified in acoustic data were targeted for carbonate and sediment sampling during three cruises on the CCGS John P. Tully in 2011, 2015, and 2017 and are assigned station numbers in Table 1 (e.g., 2017-STN15). Given weather and time restraints, 11 authigenic carbonate samples were collected at seven of the 18 seep sites. Authigenic carbonate samples were collected at water depths ranging from 507 to 1003 m, either by a grab sampler developed by Institutt for Kontinentalsokkelundersøkelser (Trondheim, Norway), capable of collecting a 1 m³ sample of relatively undisturbed sediment or rock, or from a piston corer (Table 1). In the northern QCF region, a large carbonate sample (2017-STN15) was collected in 2017 on a seep site (seep 3; Table 1) situated on a bathymetric ridge where a gas plume was present, extending 250 m from the seafloor in height (Figs. S3A, S3B [^{footnote 1}]). Within the central QCF region, two authigenic carbonate samples were collected in 2017 (2017-STN40 and 2017-STN41) at an active seep site (seep 8; Table 1; Figs. S8A, S8B) near the mouth of the Chatham Sea Valley (Alaska), and one was collected at the active Chatham Fan gas seep site (seep 7; 2017-STN04; Table 1; Figs. S7A, S7B). Based on in situ photographs, bacterial mats and bivalve shells were also observed in this region. Farther south in the central QCF region, three carbonate samples (2015-STN34, 2015-STN35, and 2017-STN03) were collected in 2015 and 2017 in the vicinity of the mud volcano identified by Barrie et al. (2020) west of the Dixon Sea Valley (seep 11; Figs. 1, 2B), where a plume was detected in the 18 kHz data and previously detected extending 700 m into the water column in 2015 and 350 m in 2017 (Barrie et al., 2020; Figs. S11A, S11B, S18). Two additional carbonate samples (2017-STN65 and 2011-STN43) were collected from the southern QCF region located near Haida Gwaii, south of the Queen Charlotte terrace and between the Dellwood fault and the QCF (Fig. 1) at seeps 17 and 16, respectively (Fig. S16). Based on in situ photographs, these sites are characterized by carbonates, lithified mudstone, gravel, and clam beds. Sediment recovered from a piston core suggests a mixture of grain sizes from clay to gravel and igneous material, including angular granitic or granodiorite clasts.

Thin sections were investigated microscopically for textural, compositional, and morphological features of carbonate precipitation. Mineralogy was determined by X-ray diffraction (XRD) using a Philips XRD with graphite monochromator at 40 kV and 45 mA as described in Prouty et al. (2016). Step scans were run from 5° to 65° 2θ with 0.02° steps, using CuKα radiation and a count time of 2 s per step following Hein et al. (2013). XRD digital scan data were analyzed using the Philips X’Pert HighScore search-and-match function to identify minerals. The XRD 100 intensity peak at 2θ for calcite was 29.4 degrees; dolomite, 30.8 degrees; quartz, 26.6 degrees; and aragonite, 26.2 degrees. Relative contributions (relative percentage of carbonate minerals within each sample) are reported as major (>25%), moderate (5%–25%), and minor (<5%). Calcium carbonate content, reported as weight percent (wt%), was determined using a coulometer at the USGS Pacific Coastal and Marine Science Center (Santa Cruz, California).

A handheld pneumatic drill (Dremel) was used to sample various components of the authigenic carbonates. The resulting powdered carbonate was collected and analyzed for stable carbon (δ¹³C) and oxygen (δ¹⁸O) isotopic composition using a ThermoScientific Kiel IV carbonate device interfaced to a ThermoScientific MAT-253 dual-inlet isotope ratio mass spectrometer at the University of California, Santa Cruz (UCSC) Stable Isotope Laboratory. Stable isotope values are reported in per mil (‰) relative to the international reference Vienna Peedee belemnite. Analytical uncertainties (1σ) are 0.05‰ for δ¹³C and 0.10‰ for δ¹⁸O.

Subsamples from the authigenic carbonate samples were analyzed for radiogenic strontium isotope (⁸⁷Sr/⁸⁶Sr) composition. In brief, using an agate mortar and pestle, the carbonate was homogenized, and dissolved in 2 N HNO₃ and filtered to minimize contamination from associated insoluble fractions (e.g., clays and other silicates) entrained in the carbonate. Samples were digested in sealed Teflon vessels following the protocol from Bullen et al. (1996), whereby Sr was separated from other ions using a Bio-Rad AG-502-X8 cation exchange resin with HCl as the eluent. The purified Sr was then converted to nitrate form, taken up in 30 μL of 0.15 M H₃PO₄, and loaded onto a Ta ribbon. The isotopic composition was measured on a Finnigan MAT 261 multicollector mass spectrometer using a static collection mode at the USGS facility at Menlo Park, California (Bullen et al., 1996). The ⁸⁷Sr/⁸⁶Sr values have been corrected for analytical fractionation to the standard ⁸⁸Sr/⁸⁶Sr ratio of 8.37521, and measurements are precise to ±0.00002 at the 95% confidence level.

Dating of select authigenic carbonates was conducted at the British Geological Survey Natural Environment Research Council (NERC) Geochronology and Tracers Facility (Keyworth, UK) using an analytical protocol outlined in Prouty et al. (2016) and modified from Edwards et al. (1987) and Shen et al. (2002). Powdered carbonate samples were processed via total dissolution techniques, with isotope ratios measured on a Thermo Neptune Plus multicollector inductively coupled plasma mass spectrometer (ICP-MS) relative to a mixed ²²⁹Th-²³⁶U tracer calibrated against gravimetric solutions of CRM 112a uranium reference material and Ames Laboratory high-purity Th. The authigenic carbonates can incorporate detrital material that carries ²³²Th and an associated amount of initial ²³⁰Th that is not related to the in situ decay of ²³⁴U. Therefore, a correction is required to calculate a reliable carbonate precipitation age. In lieu of site-specific detrital isotopic composition (e.g., Bayon et al., 2015; Teichert et al., 2003), the detrital isotopic composition was based on measured sediment values from Prouty et al. (2016) from 385 and ∼1600 m below sea level along the U.S. Atlantic margin, with the detrital ²³⁰Th/²³⁸U value for each sample or site interpolated based on water depth. This correction accounts for both the ²³⁰Th contained within the detrital grains themselves and the so-called hydrogenous ²³⁰Th produced by the decay of U dissolved in the water column, which is adsorbed onto the surfaces of sediment particles. Using the decay constants of Cheng et al. (2013), U-Th age calculations were performed using an in-house Microsoft Excel spreadsheet.

The void gas from sample 2017-STN03 at seep 11 was sampled directly from gas trapped in the core liner at 160 cm downcore using a gas-tight syringe, and injected into an evacuated 30 ml serum vial equipped with septa for geochemical analysis, thereby keeping it sealed from mixing with ambient air. Sediment collected for headspace gas analyses was sampled between 205 and 215 cm downcore and was collected by placing a 10 cm section of core from the core cutter into a septa-equipped 1 L sample can. Seawater was added to the brim of the can, and then 200 ml of this water was removed to create a 200 ml headspace. The can was sealed and shaken for 5 min. In preparation for geochemical analysis, gas was extracted using a syringe. Both void gas and headspace gases were analyzed by gas chromatography following the methods of Kvenvolden (1993) at the USGS Pacific Coastal and Marine Science Center. Gas samples collected directly as gas (void gas) and as sediment gas are reported in parts per million by volume (ppm); methane (C₁), ethane (C₂), propane (C₃), isobutane (iC₄), normal butane (nC₄), isopentane (iC₅), and normal pentane (nC₅). Methane carbon isotopic composition was performed on a laser-based Los Gatos Research methane isotope analyzer (DLT-100). Samples were injected into a known-volume cell where the length of the laser path is determined as a function of the carbon isotope. The results are reported as per mil (‰) relative to Vienna Peedee belemnite. The accuracy of the measurements is ±0.5‰.

A total of 18 active seeps were discovered in this study based on detection of bubble plumes in the water column (Fig. 1; Table 1). Backscatter anomalies were not observed at active seep sites or at sites where authigenic carbonates were collected. Instead, in areas of rugged topography along the continental slope, backscatter anomalies were observed presumably due to changes in sediment lithology such as grain size and sand content in and around channels, galleys, and submarine canyons. In this section, seep distribution is described in terms of the seafloor geomorphology and substrate characteristics along distinct kinematic sections of the QCF (following Brothers et al., 2020). The northern section extends from the Yakobi Sea Valley to the Chatham Sea Valley, the central section from the Chatham Sea Valley to the Dixon Sea Valley (to the USA-Canada border), and the southern section from the Dixon Sea Valley to the southern end of Haida Gwaii (Fig. 1).

Six seeps were identified along the northern section of the QCF and to the west of Baranof and Chichagof Islands (Figs. 1A, 2A). The QCF in this region crosses the upper slope and outer shelf, and deformation appears to be concentrated within a narrow zone along the seafloor fault trace (Brothers et al., 2019, 2020). The continental slope to the west of the QCF is characterized by networks of submarine canyons and submarine landslide scars but displays little or no evidence of active tectonic deformation away from the mapped trace of the QCF. Seeps 1, 2, and 6 are located along the shelf edge at distances of 1.5–3.7 km away from the QCF, seeps 3 and 4 are located near the crests of ridges that separate submarine canyons or gullies, and seep 5 appears to be on the QCF (Fig. 2A). MCS profiles crossing the canyonized morphology (e.g., Fig. 3) reveal stratigraphic layering beneath the ridgelines that has been disrupted and show bottom-simulating reduction in acoustic frequency and amplitude, or acoustic wipeout, that we infer to be associated with sediments containing free gas trapped below the base of the gas hydrate stability zone (GHSZ) (Chand and Minshull, 2003). The bottom-simulating change in character, which is most evident in the chimney meta-attribute calculation, cross-cuts reflectors under the ridgelines to either side of seep 4 but is not evident directly below the seep site. A narrow zone of vertically disrupted strata, or pipe (Fig. 3B inset), extends vertically above a cone of acoustic wipeout beneath the ridgeline.

Six seeps were identified along the central section of the QCF (Fig. 2B). The continental slope along this section is more rugged and morphologically diverse and displays a greater degree of tectonic influence than the northern section (Brothers et al., 2019; Greene et al., 2019). Although the seafloor trace of the QCF remains confined to a narrow zone of deformation, a series of elongate ridgelines to the west of the QCF are bounded by steep escarpments and are either incised by or have deflected the processes of submarine canyon erosion (Fig. 2B). The elongate ridgelines are of predominantly tectonic origin (e.g., Tréhu et al., 2015), and all seeps identified along the central section are located on ridgelines. Seeps 7 and 8 are along a ridge crest that defines the edge of the westward-facing escarpment, but the ridge appears to be buried by sediment accumulation along its eastern flank (Figs. 4–6). MCS profiles crossing both seep sites reveal complex patterns of acoustic wipeout, and noticeable shoaling of the wipeout zone occurs beneath both seep sites. The chimney analysis shows evidence for disrupted reflectors, or pipes, extending vertically from the wipeout regions to the seep sites (Figs. 4 and 5). Although these high-resolution MCS profiles do not show definitive evidence for bottom-simulating changes in character, cross-cutting reflections are observed along the flanks of the anticlinal ridge and may be related to the base of the GHSZ (Figs. 4–6; Fig. S17 [^{footnote 1}]). Seeps 9, 10, and 11 are located along a fault-controlled anticlinal ridgeline that extends across the international border (see Tréhu et al., 2015, their figure 6). The ridgeline is 10–20 km west of the QCF and does not appear to be directly related to the QCF in its present-day configuration. Sites 10 and 11 appear to be clustered sets of seeps along linear alignments ∼500 m long (Fig. 7). The sparker MCS data did not provide adequate penetration and imagery of the internal structure of this steeply flanked ridgeline. Seep 12 is located along an anticlinal structure similar to that of seeps 9–11, and a nearby airgun MCS profile (EW9412 line 1263; Fig. 7) contains a bottom-simulating reflector beneath the projected seep site that was previously inferred to represent a zone of free gas trapped beneath the base of the GHSZ (Riedel et al., 2020).

Based on this study, a total of six seeps were identified along the southern section of the QCF, with all seeps except 13 and 18 located along structurally controlled morphology within the Queen Charlotte terrace (Fig. 2C). Along the southern section of the QCF, the slope morphology is complex and rugged. This most likely reflects an increase in the degree of oblique convergence accommodated by the QCF in the south. Seep 13 lies ∼4 km southwest of the QCF and is associated with a slide mass located on a bathymetric high formed on the west side of the fault, possibly emanating adjacent to a slide block within the slide zone (Greene et al., 2019). Seep 14 occurs along the trace of the QCF, seep 18 is located along the shelf edge of southernmost Haida Gwaii, and seeps 16–17 occur at the top of two broad, rounded bathymetric promontories (Fig. 2C) that were constructed through the transport of both sediment (mud) and fragments of mafic rock to the surface. Previous observations from Barrie et al. (2020) documented mud rivulets on the upper flanks of the cone-like structures and the presence of basaltic fragments with manganese crusts as well as authigenic carbonate from the summits of these edifices. These observations are consistent with evidence used by Paull et al. (2015) to identify mud volcanoes in the Beaufort Sea (Arctic Alaska and Canada). Seep 15 occurs near the top of the western flank of a similar-sized round, but flat-topped, guyot-like mass (e.g., wave-base eroded; Barrie et al., 2013) ∼10 km northwest of these cone-like features. The Pacific crust is only ca. 5 Ma at this latitude and gets progressively younger southward to the Explorer Ridge. Seismic reflection data along the southern section of the QCF are limited, but previous studies have documented evidence for transpressional deformation, anticlinal folding, elevated geothermal anomalies, and gas hydrate (Hyndman and Hamilton, 1993; Riedel et al., 2020; Rohr et al., 2000). For example, the occurrence of gas hydrates along the Haida Gwaii margin was recently described in Riedel et al. (2020) using seismic reflection data to identify gas hydrate–related bottom-simulating reflectors.

The carbonates sampled from along the QCF represent a suite of authigenic carbonates consisting of fossils, peloids, and intraclasts embedded in a microcrystalline matrix. Aragonite dominates the authigenic carbonates and is the major component in all the samples except for 2017-STN03 (seep 11) where dolomite is the major carbonate phase (Table 2). Dolomite is present in two other samples (2015-STN35-2 and 2017-STN41-3) as moderate and minor (2017-STN40-3SH) contributions. The contribution from calcium carbonate (CaCO₃) ranged between 7% and 97%, suggesting a range of contributions from sedimentary detritus. The carbonates are a combination of highly compacted, low-porosity mudstones and breccias with clastic grains and carbonate pore filling. Moderate contributions were from quartz, low-Mg calcite, and plagioclase, and minor contribution from kaolinite. In addition to quartz, framework grains include K-feldspar and plagioclase, as well as amphibole, mica, and pyroxene (Table 2). Both foraminifera and shell fragments were identified in thin section (samples 2017-STN40 and 2017-STN41; Figs. 6D, 6F). Samples identified as containing fragments of bivalve shell material typically were enriched in CaCO₃ (e.g., sample 2017-STN40-2SH), whereas anomalously low CaCO₃ characterized quartz-rich subsamples (e.g., sample 2017-STN40-4-1). As described below, several samples contained fractures, veins, and detrital components, with carbonate cement occupying these pore-filled zones. The veins are typically filled with acicular or microcrystalline aragonite, representing two phases of aragonite where the latter signifies diffuse flow of fluids (Peckmann et al., 2009).

The textural, compositional, and morphological features of carbonate precipitation based on thin sections of six individual authigenic carbonate hand specimens include the following: (1) macroporosity with voids and fractures, with cross-cutting veins and acicular pore filling of highly cemented, matrix-supported mudstone with a clotted texture in groundmass (sample 2017-STN15; Figs. 8B, 8C), indicative of a high-energy environment and multiple stages of precipitation; (2) highly compacted, detrital-rich breccia with intact foraminiferal shells (sample 2017-STN40; Fig. 6D) with the lack of veins or cross-cutting representing a less-developed precipitation stage; (3) poorly sorted, fine-grained aragonite breccia with secondary infilling by microcrystalline and acicular pore filling, and inclusion of bivalves and diatoms (sample 2017-STN41; Fig. 6F); (4) lithified, well-sorted dolomite mudstone with scattered clastic grains and highly cemented “rhombs” in a fine-grained crystalline matrix (sample 2017-STN03; Figs. 7B, 7C); (5) poorly sorted and highly compacted, microcrystalline aragonite with low porosity, and angular quartz detrital components (sample 2011-STN43; Fig. 9D); and (6) breccia with cross-cutting veins, acicular pore filling, and micrite and/or microcrystalline pore-filling cement, suggestive of multiple generations of carbonate precipitation, and clotted microfabric consistent with peloid formation from microbial metabolic processes (Kiel et al., 2013) (sample 2017-STN65; Fig. 9F).

Authigenic carbonate δ¹³C values ranged from −46.82‰ to −2.96‰ (Fig. 10). Samples that were identified as containing shell material were consistently enriched in ¹³C (Table 2). The presence of shell material and seawater-like ⁸⁷Sr/⁸⁶Sr values, as discussed below, suggests that the authigenic carbonates most likely formed at or near the sediment interface. As such, incorporation of seawater dissolved inorganic carbon during precipitation could help explain the high δ¹³C values in these samples relative to methane-derived carbon. However, additional sources of ¹³C enrichment could include AOM, which can generate ¹³C enrichment of the residual CH₄ due to preferential microbial oxidation of lighter CH₄ (Claypool et al., 1985; Whiticar and Faber, 1986), carbonate reduction methanogenesis from highly enriched ¹³CO₂ sources, as well as incorporation of marine carbonate (0‰) and sedimentary organic carbon (−20‰). The carbonate δ¹⁸O values ranged from −6.89‰ to 6.44‰, with samples 2017-STN03 and 2017-STN41-3 accounting for the ¹⁸O-depleted values (Fig. 10). Bottom-water δ¹⁸O values are not available from these sites. However, the δ¹⁸O range yields unrealistically cold bottom-water temperatures (<0 °C) based on the aragonite-temperature equation of Grossman and Ku (1986). In situ temperatures from the CTD rosette casts range between 3.03 and 4.10 °C, suggesting an aragonite δ¹⁸O value between 2.6‰ and 2.85‰. Only 12% of the aragonite δ¹⁸O values fall within this range, with the majority of aragonite δ¹⁸O values suggesting input from an ¹⁸O-enriched fluid source. In contrast, the ¹⁸O-depleted samples 2017-STN03 (seep 11) and 2017-STN41-3 (seep 8) may reflect warmer temperatures and/or the influence of ¹⁸O-depleted fluids.

The range in ⁸⁷Sr/⁸⁶Sr of the carbonates was from 0.708922 to 0.709158 (Table 2). Compared to modern seawater, with ⁸⁷Sr/⁸⁶Sr of 0.709175 (Paytan et al., 1993), two of the authigenic carbonate samples, 2017-STN03 (seep 11) and 2011-STN43-1 (seep 16), yielded less-radiogenic values. These ⁸⁷Sr/⁸⁶Sr values were lower than modern seawater by as much as 0.000253, whereas the majority of samples yielded seawater-like ⁸⁷Sr/⁸⁶Sr values with an offset of <0.00008 compared to modern seawater.

Discrete subsamples of cavity-filling authigenic carbonate weighing between 0.81 and 75.14 mg were hand-drilled from five of the QCF specimens. These samples were chosen for U-Th dating given evidence of relatively pure late-stage cavity-filling carbonate cements. The specimens contain 1.99–5.62 ppm U and 0.09–0.90 ppm Th. (²³⁰Th/²³²Th) activity ratios were between 2.6 and 6.2, at the lower end of the range reported for other methane-related authigenic carbonates in the literature (Bayon et al., 2009; Crémière et al., 2013; Feng et al., 2015; Teichert et al., 2003). Therefore, a detrital correction is required, as discussed above. Following this correction, several of the specimens yielded initial (²³⁴U/²³⁸U)_i values offset from the modern seawater ²³⁴U/²³⁸U of 1.1466 (Robinson et al., 2004), an indication that the carbonates have exchanged U and/or Th with their environment since precipitation, which results in inaccurate calculated ages. Therefore, only two specimens, for which the (²³⁴U/²³⁸U)_i from late-stage carbonate precipitates was within uncertainty of the modern seawater value, were deemed to have reliable dates: samples 2017-STN65 (ca. 19 ka) and 2017-STN15 (ca. 3.6 ka) (Table S2 [^{footnote 1}]).

The composition of the void gas from this site was dominated by methane (99%), with ethane and propane contributing at 690 ppm and 6 ppm (Table 3). The δ¹³C value of the void gas was −47‰ and yielded a C₁ / (C₂ + C₃) ratio of 1230. Sediment gas extracted from the headspace of a canned core-cutter sample was also dominated by methane and had a similar methane isotopic composition (−44‰) to the void gas and a lower hydrocarbon C₁ / (C₂ + C₃) ratio of 240, a result consistent with the preferential loss of methane due to sediment degassing during core recovery.

The seismic data and associated attribute calculations along the QCF suggests that both active and dormant venting along the margin are stratigraphically and structurally controlled by anticlinal features located away from the QCF instead of narrowly focused along the QCF. The fact that seeps along the QCF were not observed in abundance may suggest that the QCF strike-slip fault is generally not a zone of gas accumulation; instead, methane traps develop away from the fault in structural and sedimentary anticlines. The occurrence of seeps away from the QCF is compatible with conceptual models of focused fluid flow and gas migration along the crest of an anticline (Paull et al., 2008), such as observed offshore northern Panama (Reed et al., 1990), where folding, thrust faulting, and sedimentation along a deformation front may facilitate fluid flow and seafloor venting (Fig. 11). The absence of abundant fluid flow along the QCF could also be partially due to the lack of transpressional and transtensional deformation along the majority of the fault’s length (Brothers et al., 2020). This more pure translational motion along the fault differs from other more obliquely striking offshore strike-slip faults and/or plate boundaries (e.g., Dupré et al., 2015), limiting the vertical deformation of surrounding strata (extensional and compressional). This lack of deformation prevents the development of structures that promote updip fluid flow and trapping. However, the active seeps and subsurface fluid-flow indicators presented in this study are extracted from data sets with limited spatial coverage and detection abilities and thus likely represent just a fraction of the total seep distribution along the QCF margin. For example, two sites located directly on the active trace of the QCF (seeps 5 and 14) were not sampled because of weather constraints, but both likely represent fluid migration along the fault zone.

Although the continental slope along the northern section of the QCF appears to be dominated by marine sedimentary processes (incised networks of submarine canyons and emplacement of fan aprons) rather than by tectonic deformation (Brothers et al., 2019), the seeps are in several different geomorphic environments. This includes localized areas of transpression along a transform fault and sedimentary loading due to rapid sedimentation from Quaternary glacial processes. In some cases, these processes manifest as seeps at topographic crests where gas is trapped, possibly beneath the low-permeability base of the GHSZ (e.g., Bünz et al., 2012; Petersen et al., 2010). Therefore, seeps sampled in this study appear to preferentially occur along the crests of ridgelines formed by (1) incision and/or aggradation surrounding submarine canyons and (2) folding and uplift along contractional structures.

We propose that sedimentary aggradation along canyon interfluves coupled with incision in the thalwegs generates lenticular beds that resemble structural anticlines and promote updip fluid migration and trapping. We propose updip migration of fluids from underlying hydrocarbon reservoirs. This interpretation is further warranted given that the authigenic carbonate precipitates and areas of active venting are located on the spines or crests of ridges (Figs. 1, 2, 7A, 8A). In cross-section, these ridges are antiformal structures that are flanked by growth faults and growth strata depicting recent deformation (Walton et al., 2018). The development of anticlinal ridges is a key feature along the QCF that facilitates focused fluid flow where the plate boundary is obliquely convergent. This scenario is similar to that observed in Southern California where fluid migration is controlled by fault-segment boundaries where folding and convergent deformation occur along strike-slip faults (Kluesner and Brothers, 2016; Maloney et al., 2015) and along other active margins characterized by folding and shortening within the accretionary wedge (Carson et al., 1991, 1994; Johnson et al., 2003; Reed et al., 1990). Upward fluid migration is further enhanced where pore-fluid pressure gradients are increased as a result of sediment loading and compaction, which facilitates transportation to the seafloor.

The high-resolution MCS data in the study region show bottom-simulating changes in character and acoustic wipeout (Figs. 3–5), which are commonly seismic indicators for the hydrate–free gas interface (Chand and Minshull, 2003) and sometimes referred to as gas wipeout. We also observe gas wipeout along the upper limbs and axes of lenticular bedding planes, suggesting updip fluid migration within a chimney structure (Løseth et al., 2009). Similar observations are found along the structurally controlled ridgelines of the central and southern sections of the margin. The neural-network chimney attribute calculations on seismic profiles help delineate broad chimney features and more focused and/or narrow zones of fluid migration referred to as “pipes” (e.g., Fig. 3B inset). The attribute calculation is focused on highlighting vertical zones of discontinuity and regions with significant vertical change in the frequency content of the data (Kluesner and Brothers, 2016). Such zones highlighted by the attribute calculation in the data commonly simulate the seafloor bottom reflection, likely delineating the isothermally controlled transition from frozen hydrate to free methane gas. Due to the high-frequency content of the sparker data, the presence of free gas can result in washout of higher frequencies due to scattering attenuation of the P waves. In addition, narrow vertical zones of discontinuity (pipes) highlighted by the chimney attribute also extend vertically toward the seafloor near the apex of ridges and anticlinal folds (Figs. 3, 4, 5). The pipe features commonly root into the broad chimney zones, likely acting as focused pathways for fluid and/or gas escape (e.g., Fig. 5B inset). Such pipe-like features have been observed in seismic data from other regions and have been interpreted as fluid-migration pathways and episodic gas blowout structures (Løseth et al., 2009).

The observation of chimneys and disappearance or shoaling of the inferred free gas–hydrate interface below seep sites may be due to heat transport by ascending fluids that contribute to zones of hydrate destabilization. The warmer ascending fluids can trigger destabilization by heating the adjacent sediments, promoting the release of hydrate-derived methane (Davis et al., 1995; Mann and Kukowski, 1999; Wood et al., 2002). This scenario provides a potential source of ¹⁸O-enriched fluids through hydrate dissociation, which releases hydrate water and contributes to enrichment by as much as 3.5‰ (Maekawa, 2004). For example, along the Cascadia margin in the northeastern Pacific (Johnson et al., 2015, and references therein), active areas of methane emission occur at the upper limit of the GHSZ near a water depth of 500 m, suggesting that gas hydrate dissociation is a major driving force of methane release along this margin (Johnson et al., 2015). Despite the observation that the upper limit of the GHSZ is near a water depth of 500 m and hydrate should be stable at our sampling sites, it is possible that heating from deeply sourced fluids is exerting a thermal control on hydrate stability and shifting the hydrate phase boundary (Chand and Minshull, 2003). It should also be noted that a direct comparison to the Cascadia margin is confounded by a lack of thermal data along the QCF, which varies in age and therefore thermal structure, which also exerts a primary control on the position of the GHSZ. The QCF carbonate δ¹⁸O values are ¹⁸O enriched by between 0.13‰ and 3.59‰ relative to predicted values based on the aragonite-temperature equation of Grossman and Ku (1986) using CTD temperatures. This offset suggests isotopic disequilibrium and the influence of an ¹⁸O-enriched fluid source. For example, given the age of the authigenic carbonate sample 2017-STN65 (seep 17; 17.55–20.00 ka), a shift in the GHSZ during the Last Glacial Maximum (LGM) can explain the ¹⁸O-enriched carbonate values. Present-day conditions indicate that seep 17 is well within the GHSZ (800 m water depth and 3.9 °C), requiring either a warming of ∼6 °C or a sea-level drop of ∼400 m, or possibly a combination. Given that an isostatic forebulge under Haida Gwaii resulted in an early post-glacial relative sea level of ∼−150 m below present (Clague, 1983; Josenhans et al., 1997; Shugar et al., 2014), a warming of only 4 °C would have been required to shift the GHSZ. There is evidence of surface warming in the area, with Globigerina bulloides recording a range of temperatures from 7.3 to 8.8 °C from 21 to 14.7 ka (Taylor et al., 2014). Therefore, it is possible that warming acted in combination with an isostatic crustal forebulge during the LGM to drive hydrate destabilization along the QCF near Haida Gwaii. However, additional pore-fluid chemistry is needed to tease apart the different processes impacting the fluid chemistry, including hydrate dissociation (Davidson et al., 1983) as well as the potential role of clay-mineral dehydration during the smectite-illite transformation that releases interlayer water (e.g., Hensen et al. (2004).

Compaction and degradation of organic-rich sediment, continuously forming hydrocarbon reservoirs (Tinivella and Giustiniani, 2012), and dissolution of clathrates or frozen gas-hydrate systems suggest that fluids could originate from multiple sources as well as evolution pathways along the QCF. The range in δ¹³C values from the QCF authigenic carbonates captures mixing between carbon sources (Campbell, 2006). There are three clusters of δ¹³C values: δ¹³C > −20‰, −30‰ > δ¹³C > −40‰, and δ¹³C < −40‰ (Fig. 10). Broadly, there are three potential carbon pools with distinct δ¹³C values: a heavy marine inorganic carbon pool, enriched thermogenic methane, and depleted microbial methane (Whiticar, 1999). Depleted microbial methane (<−50‰) results from degradation of microbial organic matter within the methanogenic zone (Tissot and Welte, 1984), whereas thermogenic methane results from thermal decarboxylation during burial of organic matter, producing relatively enriched methane (>−40‰) (Bernard et al., 1978). Along the Cascadia margin to the south, both microbial methane and thermogenic hydrocarbon sources have been identified (Sample and Reid, 1998; Teichert et al., 2005), helping to explain the large spatial heterogeneity and complex hydrology of the margin (e.g., Torres et al., 2009). While limited in scope, the single void-gas δ¹³C measurement from sample 2017-STN03 suggests methane of thermogenic origin, whereas the hydrocarbon composition, defined as the C₁ / (C₂ + C₃) ratio, is more consistent with microbially produced methane (Milkov and Etiope, 2018). However, the paired hydrocarbon composition and δ¹³C-methane value from sample 2017-STN03 are outside the anticipated mixing line between thermogenic and microbial end members (Bernard et al. 1978), indicating that additional processes should be considered, including anaerobic methane oxidation, which can generate ¹³C-enriched methane (Claypool et al., 1985), or, for example, carbonate reduction methanogenesis from highly enriched ¹³CO₂ sources like those measured at Integrated Ocean Drilling Program Sites U1327 and U1329 along northern Cascadia (Pohlman et al., 2009). Additional gas samples are needed to further determine the role of microbial alteration and to identify distinct carbon sources. However, based on the geophysical data, we propose that along the QCF, a deep-rooted methane signature can occur via the presence of active venting and formation of authigenic carbonates on the spines or crests of ridges (Fig. 11). For example, vertical pipes evident in the chimney meta-attribute calculations overlie regions of gas wipeout beneath the ridgelines and suggest that fluid and/or gas migrates from a deeper reservoir zone along a focused pathway to the seafloor (e.g., Fig. 4).

The majority of authigenic carbonate samples yield seawater-like ⁸⁷Sr/⁸⁶Sr values, consistent with previously recovered authigenic carbonates from the Cascadia slope (Joseph et al., 2013), where shallow precipitation of aragonite triggered by AOM forms authigenic carbonates in anoxic but seawater-ventilated environments, such as observed at Hydrate Ridge (offshore Oregon, USA; Greinert et al., 2001). This scenario is consistent with aragonite formation in open systems at or close to the sediment-water interface and in situ brecciation of weakly consolidated sediment. The incorporation of shells also indicates in situ cementation of chemosynthetic communities that live at the sediment-water interface. In contrast, the depleted radiogenic Sr-isotope composition from samples 2011-STN43 (seep 16) and 2017-STN03 (seep 11) suggests mixing with fluids other than modern seawater, such as deep-sourced mud breccia flow (e.g., Klasek et al., 2019) and meteoric fluids from groundwater circulating at depth due to changes in hydraulic head linked to ice-sheet dynamics (e.g., Person et al., 2007).

Based on an integrated approach using geology, geophysics, and geochemistry, we demonstrate that anticlinal features located away from the QCF stratigraphically and structurally control active seepage and are indicators of past fluid expulsion. This is in contrast to the assumption that seeps along strike-slip fault systems are primarily focused along faults. Upward fluid migration of methane is facilitated by sediment loading, compaction, and tectonic and sedimentary processes that result in folded strata and anticlinal traps. Fluids are subsequently transported upward through chimneys, eventually oxidized by archaea carrying out AOM and favoring precipitation of authigenic carbonates. Carbonate microstructures in the thin sections including brecciation and pore-filling cement suggest rapid release of fluids and/or gases into weakly compacted sediment at or near the sediment interface. The carbonate chemistry of the authigenic carbonate samples highlights both microbial and thermal degradation of organic matter in continental-margin sediments along the QCF and subsequent fluid seepage. The majority of aragonite carbonates formed in the shallow, anoxic subsurface environment with elevated bicarbonate concentration as a result of AOM and fed by the rise of methane-rich fluids. Future studies focused on localized folding and faulting will further test our assumption that anticlinal structures along the QCF create high-permeability regions where conduits form and facilitate upward migration of fluids. Ideally, future work would include submersible dives to collect discrete gas samples and targeted coring to further investigate fluid sources and the presence of hydrates. In doing so, we can greatly expand our predictive capabilities of the geographic extent and geologic origin of focused fluid flow and methane venting along transform boundaries.

The authors thank officers and crew of the CCGS Vector and CCGS John P. Tully; M. Baker (USGS), R. Garrison (UCSC), J. Fitzpatrick, (USGS), N. Vokhshoori (UCSC), and C. Maupin (Texas A&M University, College Station, Texas, USA) for laboratory and sample assistance; and J. Pohlman (USGS) for helpful comments. Input from two anonymous reviewers substantially improved the manuscript. The USGS Coastal and Marine Hazards and Resource Program funded this study. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government. Additional geochemical and geophysical data to support this project can be found in Prouty et al. (2019) and Balster-Gee et al. (2017a, 2017b), respectively.