The Cascadia Subduction Zone (CSZ) is an active plate boundary where the Juan de Fuca, Gorda, and Explorer plates subduct beneath the North American plate. This convergent boundary extends over 1,300 km from the Nootka fracture zone offshore British Columbia (Canada) in the north (Audet et al., 2008) to the Mendocino triple junction offshore northern California (United States) in the south (Figure 1). A large accretionary wedge composed of both terrigenous and marine sediments has been accreted to the North America plate along the CSZ. This wedge contains a series of folded and thrusted ridges striking subparallel to the continental slope that form large anticlines (Cochrane and Lewis, 1984; Goldfinger et al., 1991). There is also significant along-strike geophysical and geologic heterogeneity that has been linked to different modes and rates of stress accumulation and release along the margin (see summary in Walton et al., 2021). Such variability could reflect differences in sediment and pore-fluid properties that ultimately influence fluid migration due to rapid compaction and dewatering rates (Tréhu et al., 2004; Han et al., 2017). For example, high glacial sediment flux during the Pleistocene significantly expanded the size of the outer wedge in northern Cascadia, between Juan de Fuca and Astoria Canyons (Silver, 1972; Barnard, 1973).

Broadly, the offshore accretionary wedge can be divided into two margin-parallel zones separated by the outer arc high (OAH), which is defined as a broad structural high that bounds the seaward edge of the shelf forearc basins (e.g., von Huene and Scholl, 1991; Clift and Vannucchi, 2004; McNeill and Henstock, 2014). While the OAH typically follows the shelf break along the 200-m isobath in Cascadia, deviations from the shelf break, both landward and seaward are observed, and may be linked to regional variations in the width of the outer wedge. Watt and Brothers (2020) identified four distinct morphotectonic regions based on careful examination of geophysical examination of seismic profiles to describe variable accretionary wedge width, shape, and structural vergence. These four regions include: (1) Vancouver Island, British Columbia, Canada (average width, linear wedge, seaward and mixed vergence); (2) Washington, United States (higher width, concave wedge, landward and mixed vergence); (3) northern and central Oregon, United States (average width, linear and convex wedge, mixed and seaward vergence); and (4) southern Oregon and northern California, United States (lower width, convex wedge, seaward and mixed vergence). Regional patterns and differences in the tectonic evolution of the outer accretionary wedge are also linked to Cascadia margin-wide heterogeneity in the distribution and character of the BSR (Tréhu and Phrampus, 2022). For example, BSRs were observed beneath anticlines where internal deformation of the wedge replenishes the supply of methane to the GHSZ at the crest of the anticlines by compacting sediments and providing pathways to the shallow subsurface (Tréhu and Phrampus, 2022). Therefore, tectonic factors influencing the occurrence and disruption of the BSR and associated fluid expulsion are also important in constraining the spatial distribution of cold seeps along the Cascadia margin.

The seep database includes previously published data documenting over 2,850 active seep flares on the seafloor using ship-based multibeam echo-sounding (MBES) sonar (Johnson et al., 2015; Riedel et al., 2018; Merle et al., 2021), and 902 new seep flares identified during surveys by the NOAA Ships Rainier in 2018 (surveys H13117, H13118, H13119, H13137, and H13206) and Fairweather in 2019 and 2021 (surveys W00474, W00475, and H13549) (Conrad and Rudebusch, 2023; Supplementary Table S1). Bubble plumes associated with seeps were mapped in the MBES water column backscatter data using QPS FMMidwater software, using techniques and identification criteria similar to those described by Merle et al. (2021). The geologic setting of the bubble plume locations was determined from interpretation of the morphologic and bottom characteristics of the seafloor from multibeam bathymetry and backscatter data (https://www.ncei.noaa.gov/maps/bathymetry/) and from subbottom data derived from seismic reflection profiles (Balster-Gee et al., 2023a). In addition, evidence for both active and relict (or dormant) seep indicators including authigenic carbonates, vesicomyid clams, tubeworms, and microbial mats, were included from seafloor observations and collections made using remotely operated vehicles (ROVs) (Supplementary Table S1). The inclusion of seep evidence from ROVs proved useful for revealing additional seep locations that were not detected by traditional ship hull-mounted multibeam echosounders. Therefore, seep-distribution reported here is a presence-pseudoabsence dataset, not presence-absence. Diffuse seepage, where gas concentrations are not sufficient for gas phase emissions, was not included given the difficulty in identifying diffuse seepage acoustically. Following previously published methods of Johnson et al. (2015) for removing possible double-counted seeps due to overlapping surveys, we aggregated seeps found within 300 m of each other into single seep “emission site,” yielding a dataset of 880 seep emission sites for a new composite database (Supplementary Table S1). To account for variance in surface area of the depth contour polygons and the effect of unequal survey effort across depth bins, we included only seeps detected via MBES so that we could normalize the number of seeps with the corresponding area mapped and determine seep density.

Survey footprints derived from the bathymetry of the corresponding MBES surveys used to acquire seep locations were loaded into a GIS and merged into a single polygon layer. This layer was then used as a mask to extract a continuous 100-m resolution depth surface from a global multiresolution bathymetric elevation model (Ryan et al., 2009), representing the area mapped corresponding with the seep location data. Depth contour polygons were created from this layer at 100-m intervals and seep emission sites were spatially joined to their respective depth contour polygon. Seep density (per sq. kilometer) was calculated for each 100-m depth contour. A smoothed kernel density estimation curve was fitted to the data for the margin-wide (Figure 2), as well as at the regional scale (Figure 3), in order to investigate seep distribution per morphotectonic region. To test for normality in the seep depth distribution we used normal Q-Q plots and a Shapiro-Wilk test (Shapiro and Wilk, 1965; Wilk and Gnanadesikan, 1968). Patterns in seep distribution were further explored in order to identify significant seep clusters along the margin (Optimized Hot Spot Analysis tool in ArcGIS Pro 3.0, Esri). Optimized Hot Spot Analysis was used to spatially interrogate the dataset to optimally aggregate coincident features into hexbins based on the Global Moran’s I statistic of spatial autocorrelation, and then calculate a Getis-Ord Gi* statistic for each hexbin to compare the aggregated counts of seeps within each hexbin against those of its nearest neighbors. Cluster predictions are determined by the aggregation of seeps within neighboring hexbins as compared to the null hypothesis of a random distribution. We used an analysis cell size of 10 km and a fixed distance band of 20 km, as determined by the initial steps run by the Optimized Hot Spot Analysis tool to identify the optimal scale of analysis. All aggregation hexbins had a minimum of 8 neighbors, and on average 16 nearest neighbors were used with a moving window analysis to determine statistical clusters. The resultant z-scores for each aggregation hexbin identified the intensity of the clustering, with larger positive z-scores (p-value<0.05) indicating more significant seep clusters. We included only the clustering for features significant at the 95%–99% confidence levels.

In order to determine potential controls on focused fluid flow in southern Cascadia, an area with previously sparse coverage, we systematically assigned a primary (and in some cases, secondary) geologic/geomorphic attribute to each seep emission site in region 4 of Watt and Brothers (2020) Seafloor expressions of fluid expulsion or other indicators of fluid flow, such as pockmarks or authigenic carbonates but otherwise lacking evidence of active seepage, were not included. Instead, we integrated new multibeam bathymetry, backscatter data, and seismic reflection profiles to identify attributes previously inferred to be associated with substrate fluid flow: anticline, fault, failure scarp, shelfbreak, channel, bedding outcrop, and canyon (e.g., Skarke et al., 2014; Crutchley et al., 2018; Riedel et al., 2018; Prouty et al., 2020). Locations of these attributes were determined from their geomorphic expression on the seafloor (Supplementary Figure S1) and supplemented with seismic reflection profiles where available. For example, bubble plume locations were combined with the multibeam bathymetry and backscatter data (https://www.ncei.noaa.gov/maps/bathymetry/) and seismic reflection profiles (Triezenberg et al., 2016) for visualization and interpretation of the local geologic/geomorphic setting within a 50-m buffer of the seep occurrences and logged using a geographic information system (GIS) (ArcMap 10, Esri). Information on anticline locations was derived from previously published neotectonic maps in southern Cascadia (Clarke and Carver, 1992; Goldfinger et al., 1992; Goldfinger et al., 1997) but is currently unavailable for northern Cascadia.

In addition to examining seep distribution in the context of regional controls in southern Cascadia, margin-wide seep distribution was evaluated with respect to the landward limit of the GHSZ and the OAH. The OAH was designated as the location of broad structural high with a ±10-km buffer that runs along the seaward edge of the shelf forearc basins (Watt and Brothers, 2020). Joung et al. (2022) showed the landward limit of gas hydrate determined for the northern Cascadia margin based on bottom water temperatures recorded by CTD data. This type of analysis is hard to adapt for the whole margin due to lack of complete data coverage (C.D. Ruppel, USGS, written commun., 07/04/23) and we instead adopted the 500-m isobath as the nominal the landward limit of the GHSZ (Tréhu and Phrampus, 2022).

Multichannel seismic reflection datasets were utilized to reveal subsurface structure and stratigraphy under mapped seeps where available (e.g., Triezenberg et al., 2016). Additional seismic surveys were collected in 2018 aboard the R/V Coral Sea using an 88-channel Geometrics GeoEel digital hydrophone streamer and Applied Acoustic Delta Sparker sound source (Balster-Gee et al., 2023a) and in 2019 onboard the R/V Rachel Carson and utilized a 72-channel Geometrics GeoEel streamer and the Delta Sparker sound source (Balster-Gee et al., 2023b). Both sparker datasets were processed following processing procedures laid out in Kluesner et al. (2019) and included SEGD to SEGY conversion, geometry correction, UTM conversion, common mid-point (cmp) binning, FK filter, spike removal, bandpass filtering, velocity analysis, normal move-out correction, trim static corrections, cmp stacking, post-stack migration, gap deconvolution, water column mute, and automatic gain control. In addition to the workflow above, the 2019 sparker dataset included prestack deterministic deconvolution using shot-to-shot source signature information recorded during acquisition.

To optimize the detection of gas and fluid migration pathways in the vicinity of active seep sites with focused fluid flow (Heggland, 2005; Ligtenberg, 2005; Kluesner et al., 2013; Brothers et al., 2014; Connolly, 2015), a fully connected multi-layer perceptron neural-network approach was applied to 2-D high-resolution multichannel seismic (MCS) profiles at seep sites in region 4. Within the OpendTect software package, twenty-one 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 Kluesner and Brothers (2016). The chimney meta-attribute results were projected onto seismic cross-sections using a gradational color scale with transparency that reveals only the highest (∼80% and above) chimney probabilities, yielding a measurement of the probability between 0 (lowest) to 1 (highest) for presence of a chimney structure. Two types of features were identified, 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 in areas of focused fluid flow (Kluesner and Brothers, 2016; Prouty et al., 2020).

Powdered carbonate samples for stable carbon (δ¹³C) and oxygen (δ¹⁸O) isotopic analysis were collected using a hand-held pneumatic drill to sample various components of authigenic carbonates collected along the Cascadia margin as a proxy for potential fluid sources. Isotopic composition was determined via ThemoScientific Kiel IV carbonate device interfaced to ThermoScientific MAT-253 dual-inlet isotope ratio mass spectrometer (IRMS) at the University of California, Santa Cruz, Stable Isotope Laboratory. Stable isotope values are reported in per mil (‰) relative to the international reference Pee Dee Belemnite (PDB). Analytical uncertainties (1σ) are 0.05‰ for δ¹³C and 0.10‰ for δ¹⁸O. 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 High Score search-and-match function to identify minerals. The XRD 100 intensity peaks at 20 for calcite: 29.4, dolomite: 30.8, quartz: 26.6, and aragonite: 26.2 in order to determine major (>25%) carbonate phase. Results from new analysis were combined with previously published results from Cascadia (Torres et al., 2009).

Margin-wide seep depth distribution is shown in Figure 2, based on results from the smoothed density estimation curve of seep sites per 100 m depth intervals overlaid on the relative percent-area surveyed of each depth contour (Table 1). The MBES survey coverages varies with depth and is lowest in the very shallow water on the continental shelf from 0 to 200 m, where the narrowing of the swath width makes collecting large surface areas exceedingly difficult for the typical MBES systems. In most other depth bins, however, coverage is between 60%–95%. The margin-wide seep depth distribution is left-skewed with bimodal peaks at 200–300 m and 500–600 m and long-tailed showing numerous observations of deep seeps up to 3,100 m. Non-normality of the distribution was confirmed using normal Q-Q plots, which show significant deviation of the sample points from the linear model especially in the lower and upper quantiles, and a Shapiro-Wilk test (W=0.5457, p-value<0.05; Supplementary Figure S2).

Given that margin-wide plots of seep depth distribution can mask along-margin variation in seep depth trends, we evaluated the depth distribution of the seep emission sites within four morphotectonic regions of the Cascadia margin (Watt and Brothers, 2020). The depth-distribution of each of the four morphotectonic regions shows considerable variation (Figure 3), especially compared to the margin-wide depth distribution (Figure 2). In addition to displaying the median and quartiles (i.e., inset box-and-whisker plots), the violin plots use a kernel density estimation to show the distribution shape of the data (i.e., position and relative amplitude) and probability of seep distribution along the depth range (Table 1). In regions 1 and 2, 39%–42% of seeps are found on shallow continental shelf (0–200 m), with 23%–26% on the upper slope and landward limit of the GHSZ (200–500 m), and 35% found in deeper waters and within the GHSZ (>500 m) (Figure 3). In region 3, 42% are located on the upper slope up to the landward limit of GHSZ, characterized by an interquartile range of 300–600 m, but equally 42% are found within the GHSZ, and a smaller portion on the shelf (15%). In contrast, seeps in region 4 were found in deeper water, such that 68% of the seeps are within the GHSZ (interquartile range = 600–1,000 m), 15% on the upper slope and landward limit GHSZ, and the remaining 17% on the shelf. In summary, only the median seep depth in region 3 coincides with the depth range of the landward limit GHSZ. The probability of seeps occurring within the depth range of the landward limit GHSZ is only valid in region 3. In contrast, the median seep depth in regions 1 and 2 are shallower with respect to the landward limit of the GHSZ and deeper in region 4, suggesting a margin-wide seep distribution trend toward deeper water southward along the margin (Figure 3). Similar to the margin-wide analysis, seep depth-distribution displayed non-normality in each of the respective regions, as highlighted by the shape of the violin plots (i.e., multimodal data distribution), illustrating scatter rather than concentrated around the median (Figure 3). This nonuniform distribution is consistent with previous observations (Johnson et al., 2015; Riedel et al., 2018).

Testing for patterns in spatial distribution margin-wide revealed that seeps exhibit spatial clustering, as opposed to random or dispersed distributions. Along the margin seven distinct clusters (i.e., hot spots) of high-density seeps were identified (Figure 3). All of the major clusters, except the cluster at Heceta Bank, are in areas where the OAH and landward limit of the GHZS coincide (Figure 3), as discussed in more detail below. In region 1, a single cluster was identified from 200 to 900 m depth at the head of Nitinat Canyon. At this site the OAH and the landward limit of the GHSZ converge (Figure 3 and Supplementary Figure S4), in contrast to the majority of region 1 where the landward limit of the GHSZ has a seaward offset relative to the OAH (Supplementary Figure S3). Region 2 exhibited the highest incidence of clusters, where each of the three clusters occurred at the heads of major submarine canyons, Quinault and Quillayute, Grays, and Astoria Canyons, that incise the shelf at 200 m. However, seeps in Astoria Canyon extend into the canyon and occur at a range of depths along the canyon thalweg. In contrast to regions 1 and 2, region 3 is devoid of canyons and clusters in this region occur on local structural highs, at Cascade Bench and Heceta Bank. In addition, the highest probability (i.e., peak in kernel density estimation curve) of seeps occurring within the depth range of the landward limit GHSZ is only valid in region 3, as discussed in detail below. A cluster at Coquille Bank was identified in region 4 where the OAH and landward limit of the GHSZ converge. Notably, this cluster is located where an inferred backstop boundary fault (Watt and Brothers, 2020) intersects a unique seafloor seep characterized by gas bubbles containing mantle-derived helium suspected to be sourced from the subducting Gorda plate (Baumberger et al., 2018). In comparison, no cluster was detected directly offshore Coos Bay, despite convergence of the OAH and the landward limit of the GHSZ. However, pockmarks were observed in the bathymetry in this area, suggesting potential fluid seepage in the past. South of Rogue Canyon, the OAH deviates seaward to a depth range of 800 to 2,000 m, where seeps appear to be strongly associated with the OAH rather than the landward limit of the GHSZ (Figure 3). A second cluster was identified around the head of Eel Canyon at the southern terminus of the Cascadia margin and is also where the OAH rejoins with the landward limit of GHSZ.

The seep depth distribution in southern Cascadia is unique, as captured in the violin plot of region 4 (Figure 3), where most seeps (62%) occur within the GHSZ (>500 m deep), and in an area where there is a distinct seaward divergence of the OAH into deeper water (Supplementary Figure S3). Given the lack of overlap between the OAH and landward limit of the GHSZ in southern Cascadia, the deeper distribution of seeps was investigated in the context of site specific geologic-geomorphic controls. Seeps could have associations with multiple geologic-geomorphic attributes and were therefore assigned up to two primary attributes (Supplementary Table S2). Of the seven types of geologic-geomorphic attributes evaluated for facilitating fluid flow, more than half (52%) of the 278 seep emission sites in southern Cascadia occur on anticlines, with the majority of anticlines located within or seaward of the OAH (Figure 3). In comparison, 16% of seeps were associated with canyons, 14% with the shelfbreak, 11% with seafloor failure scarps, 7% in channels or gullies, and 7% at sites of bedding outcrop. While assigning site specific geologic attributes to all seeps on the Cascadia margin is outside the current scope of the study, we did find that 44% of the seep emission sites in region 3 were located on anticlines, suggesting the strong association of seeps on anticlines carries to morpho-tectonic regions other than region 4.

To further investigate the link between the occurrence of seafloor seepage and underlying geologic structure in southern Cascadia, we combined high-resolution bathymetry, water-column, and two-dimensional MCS data at two seep flare sites in region 4 as representative examples of detection of gas reservoirs and fluid migration pathways. The first example is based on water-column data from 2018 NOAA Ship Rainier MBES survey, where three seep flares were detected at ∼638 m depth, with fluid emissions reaching a height of 260 m above the seafloor (Figure 4A). Bathymetry from the associated MBES survey revealed an uplifted and rough seafloor character underneath the seep cluster and the high-intensity backscatter values accompanying these patches suggests the rough seafloor consists of authigenic carbonate precipitation. The seafloor morphology is mostly likely controlled by a broad anticlinal structure, as depicted in the MCS data (line 39 from the R/V Coral Sea). Within the anticline and beneath the active seep there is a zone of chaotic reflections and acoustic wipe-out located below the hinge line (Figure 4B). Chimney analysis reveals a broad high-probability zone below the seep field, indicating probable fluid-pathways (Figure 4B). Additional pipe-like zones of disturbance on both sides of the seep field are also highlighted as high-probability chimney zones. Fluids appear to be focused toward the hinge line where vertical pipes or faults provide a pathway to the seafloor, resulting in seafloor emissions.

The second example is from a seep site identified from water-column data derived from a 2017 E/V Nautilus MBES survey at a water depth of ∼995 m depth and flare rising 430 m into the water column. By rotating the 3D perspective from north-south orientation to east-west, a long, thin “bubble curtain” with a horizontal extent of 500 m across is revealed (Figure 5A). Bathymetric shaded relief reveals uplifted and rough seafloor character underneath seep cluster. As depicted in the MCS data (line 38 from the R/V Rachel Carson), this seep site is above a normal fault on the flank of a broad anticlinal structure (Figure 5B). As described above, chimney analysis indicates the presence of probable fluid-pathways beneath a hydrate cap (BSR) and likely updip fluid migration through a chimney zone to the seafloor (Figure 5B). The seismic profile shows that at depth, disruption of the BSR coincides with the location of a deeper fault, as well as the presence of a relic (“double”) BSR.

Authigenic carbonate δ¹³C values range from −56.73‰ to 10.03‰ with an average value of −25.21‰ ± 15.39‰ (n=324), whereas δ¹⁸O values range from −1.85‰–9.15‰ with an average value of −3.70‰ ± 1.67‰ (n=324) (Supplementary Table S3). Aragonite is the dominant carbonate phase of the authigenic carbonates and is the major component in 57% of the samples, with low and high-Mg calcite accounting for 29% and dolomite present at 13%. There is a δ¹³C range of over 60‰ from samples collected along the margin, and δ¹³C values correlate with both depth and latitude, with statistically significant (p<0.01) Pearson-Product Correlation Coefficient (p<0.01) of −0.49 and −0.48, respectively. In contrast, no correlation was found between authigenic carbonate δ¹⁸O values with either depth or latitude. Based on the aragonite-temperature equation of Grossman and Ku (1986) and seawater δ¹⁸O value of 0.08‰ SMOW (Bohrmann et al., 1998), the average δ¹⁸O-derived temperature is 5.0°C ± 3.7°C, overlapping with CTD bottom water temperatures at a depth range between 150 and 1,600 m (Baker and Prouty, 2022; Prouty and Baker, 2022).

The identification of new seep sites in central and southern Cascadia in this study extends our understanding of seep depth distribution along a subducting margin. The margin-wide depth distribution presented here is consistent with other findings that seep density peaks around 500-m depth (Merle et al., 2021), but the non-normality contradicts the assumption of symmetrical probability distribution around ∼500 m. Instead of exhibiting symmetrical probability distribution around ∼500 m, the margin-wide depth distribution is left-skewed with a long-tail highlighting numerous observations of deep seeps up to 3,100 m. On a regional basis, the probability of seeps occurring within the depth range of the landward limit GHSZ is only valid in region 3 (Figure 3). While the seep distribution is non-random, the clustering described above is not restricted to the 500-m isobath but rather at distinct clusters with unique geologic and geomorphic characteristics. Whereas Merle et al. (2021) normalized seeps to percent area mapped for only a subset of eight multibeam surveys, we normalized seep occurrence to area mapped using all the co-located footprints from 30 multibeam surveys. Between 200–1,500 m depth, MBES mapping coverage was consistently high (80%–95%), such that our results and interpretation were not significantly hindered by data gaps in survey coverage at different depth bins. Increasing multibeam data acquisition in the 0–200 m water depth range would further refine our estimation of seep density in shallow waters along the margin. The opportunistic discovery of seeps and seep-indicators (i.e., chemosynthetic biota) using ROVs at sites where MBES did not detect seeps is a reminder that, in addition to the temporal ephemerality of seeps releasing detectable gas bubbles, an absence of seeps recorded in MBES surveyed area does not indicate a true absence of seeps.

In regions 1 and 2, the majority of seeps are on the shallow continental shelf (<200 m), with the seep clusters at the heads of canyons in northern Cascadia (Figure 3). Given that the canyons incise the shelf at (or less than) 200 m, it is not surprising that seeps in regions 1 and 2 have the highest probability of occurring at depths less than 200, consistent with previous observations (Johnson et al., 2019). Seep clusters are located at the heads of major submarine canyons except in Astoria Canyon, where seeps span a depth range from 109 to 300 m at the canyon head to 480–1,026 m in the canyon. As suggested by Merle et al. (2021), the large seep distribution in Astoria is likely linked to a combination of factors, with Astoria Canyon deeply incising into the shelf, resulting in steep canyon walls, facilitating massive sidewall and channel levee failures and exposing bedding planes on the mid to lower slope (Hill et al., 2022), thus promoting seepage pathways within the permeable sediment layers. In addition to a high occurrence of seep clusters in canyons in region 2, this region is also characterized by landward migration of the OAH, most likely in response to active extensional faulting. As a result, the OAH overlaps with the landward limit of the GHSZ (Figure 3), which most likely further facilitates focused fluid in these clusters.

Seep depth distribution in the south is progressively deeper but also signifies the only region (region 3) devoid of canyons and where the probability of seeps occurring within the depth range of the landward limit GHSZ is valid (Figure 3). The predominant cluster in region 3 is at Cascade Bench, a topographic high situated between Cascade Head and Tillamook Bay at depth of 400–600 m and inferred to be a paleo–shelf edge formed during the Pleistocene lowstand (McNeill et al., 2000). A second cluster was identified at Heceta Bank, a structural high along the outer shelf caused by subducting seamounts (Tréhu et al., 2012). Within this region, the OAH is co-located with the landward limit of the GHSZ except at Heceta Bank, suggesting some discrepancy of cluster behavior on structural highs with regards to convergence of the OAH and landward limit of the GHSZ.

The overlap between the OAH and the landward limit of the GHSZ may help explain the occurrence of the two clusters in region 4, one at Rogue Canyon and Coquille Bank, a north-south trending double plunging, asymmetrical anticline (Kulm and Bales, 1969), and the other around the head of Eel Canyon. Both sites are places where the OAH and landward limit of the GHSZ converge, along with the presence of submarine canyon heads (the importance of which was discussed above). Seafloor fluid emission at Coquille Bank may also be further facilitated by a backstop fault formed from the subducting Gorda Plate (Watt and Brothers, 2020) with access to mantle-derived fluids (Baumberger et al., 2018), and associated with a mafic ridge buried beneath the sediments of the accretionary complex from about 43° to 45°N (Fleming and Tréhu, 1999). In the rest of region 4, the OAH and the landward limit of the GHSZ are separated from each other by an average of 23 ± 6.1 km as the OAH deviates seaward around the Eel River Basin. As noted in this study, many seeps occurred in the deep water in association with the OAH, and prolific seep activity was also observed proximal to the channel heads of Trinidad Canyon where a north-south trending fault cuts across the channel heads and the shelfbreak (Balster-Gee et al., 2023a), but there was a dearth of seeps in the basin between. We posit that the lack of seeps here may be due to high sedimentation rates in the Eel River Basin that cause deposition of mud and silt to drape over the Basin and Trinidad Canyon (Nittrouer, 1999; Sommerfield and Nittrouer, 1999), which may limit sufficient porefluid methane saturation. Whereas gas hydrates tend to concentrate in coarser-grained sandy/silty turbidite sediments (Torres et al., 2008), high sedimentation rate of fine-grained sediment may preclude sufficient porefluid methane saturation. This scenario is consistent with observations by Tréhu and Phrampus (2022) to explain the absence of BSR and extension of the seaward edge of the Eel River Basin well to the west of the theoretical landward edge of the GHSZ. Therefore, we suggest that the lack of overlap between the OAH and the landward limit of the GHSZ and sedimentation history precludes development of seep cluster anywhere else in region 4.

The discussion above provides an interpretive framework for understanding and potentially predicting seep clusters along a subducting margin and helps address the unresolved origin of seeps at 500 m by identifying processes and characteristics that facilitate focused fluid independently of gas hydrate dissociation in response to contemporary ocean warming. While canyons and structural highs are critical geomorphic features facilitating seep occurrence, overlap between the OAH and the landward limit of the GHSZ may be required for statistically significant clustering of seep emissions. In Cascadia, there is co-occurrence of the OAH (and shelf break) with the downdip end of seismic ruptures and interseismic coupling (McNeill et al., 2000; Booth-Rea et al., 2008; Malatesta et al., 2021). Therefore, there may be links between seep distribution and seismogenic patterns given the release and drainage of fluids from clay dehydration in the transition between locked and seismogenic zone (Saffer and Tobin, 2011). Low permeability and porosity of marine sediment could lead to significant overpressure in these environments and trapping of fluids (Sahling et al., 2008; Saffer and Tobin, 2011). In other words, there is upward/updip fluid migration from zones of peak dehydration and/or elevated pore pressure toward that OAH along permeable stratigraphic layers and fault conduits. The OAH acts as a first-order anticlinal trap that focuses fluid flow along individual faults and bedding plane (Figure 6), as supported by sub-surface data showing where widespread folds channel fluids and gases toward the anticlinal crests, as seen in the structurally controlled seeps offshore Costa Rica (Kluesner et al., 2013). This may not apply to seeps away from the deformation front since dehydration may occur at a shallow depth and near the deformation front (Hyndman and Wang, 1993; Hyndman, 2007), given the high temperatures of the very young hot subducting plate (e.g., Hyndman, 2007). In contrast, near Hydrate Ridge, Torres et al. (2004) and Teichert et al. (2005) found an increase in deep freshening eastward (i.e., relative to the toe), whereas accretion at the westernmost sites is too recent for the sediments to have undergone significant illitization. While differentiating fluids from dehydration relative to hydrate dissociation requires analysis of porefluids, the co-location of the seep clusters with the OAH and inferred downdip end of interseismic coupling in Cascadia is intriguing and warrants further investigation.

As a proxy for the composition of fluids, the authigenic carbonate isotope values suggest a range of sources, including thermogenic decomposition of organic matter and biogenic mediation from microbial activity (Claypool et al., 1985). The average authigenic carbonate δ¹³C values (−53.53‰ ± 1.95‰) at the Astoria deep sites (>800 m) are depleted relative to carbonate δ¹³C values at the shallow site (500 m; −28.51‰ ± 2.26‰), yielding a difference of 25‰, similar in magnitude between gas samples at the shallow versus deeper Astoria sites (Baumberger et al., 2018). As described in Baumberger et al. (2018), anerobic oxidation of methane (AOM), as well as incorporation of dissolved inorganic carbon (DIC) and co-existing CO₂ and its associated shift in isotope composition, can produce heavy δ¹³C values that mimic a thermogenic origin (Pohlman et al., 2009). Given gas isotopic and hydrocarbon composition indicative microbial methane at Astoria (Baumberger et al., 2018), the heavier carbonate δ¹³C values at the shallow site are in line with ¹³C-enrichment by AOM.

In contrast, at Coquille the presence of only thermogenic methane (Collier and Lilley, 2005; Baumberger et al., 2018) provides an opportunity to calibrate the range in carbonate δ¹³C values (−26.50 to −17.76‰) in the absence of mixing with microbial sources. Taken together, data from Coquille and Astoria can be used as endmember values for estimating contribution from microbial and thermogenic methane sources to carbonate δ¹³C values at other sites where both sources are identified, such as at Heceta Bank (Baumberger et al. (2018). For example, carbonate δ¹³C values greater than the range reported at Coquille can be explained by oxidation of methane and DIC incorporation, whereas carbonate δ¹³C values less than this range reflect mixing with microbial methane. Assuming endmember microbial and thermogenic δ¹³C values from Baumberger et al. (2018), mixing with a microbial methane at Heceta Bank could contribute to a depletion of 15‰, equivalent to a contribution of less than 10% from microbial methane.

The inverse correlation between δ¹³C values and depth supports earlier views that seeps from within the gas hydrate stability zone are of microbial origin and seeps from the upper slope and the shelf have a thermogenic source (e.g., Collier and Lilley, 2005; Heeschen et al., 2005; Baumberger et al., 2018). However, the discussion above also demonstrates how processes impacting fluid chemistry, such as mixing and AOM, can impact the fluid isotopic composition (Pohlman et al., 2009) and ultimately alter the archived fluid-source signature in the authigenic carbonates. In addition, gas hydrate dissociation can provide a potential source of ¹⁸O-enriched fluids, which releases hydrate water and contributes to enrichment by up to 3.5‰ (Maekawa, 2004). Samples with ¹⁸O-enriched carbonate values yield unrealistically cold temperatures (<2°C) relative to CTD bottom water data (Baker and Prouty, 2022; Prouty and Baker, 2022), indicating influence from ¹⁸O-enriched fluid source. Margin-wide carbonate δ¹⁸O values and mineralogy, however, support shallow precipitation of aragonite driven by AOM and in isotopic equilibrium with seawater.

While the occurrence of seep clusters and the justification for their distribution may help explain the spatial distribution of seeps along the Cascadia margin, it does not explain the range or probability of seep depth distribution across the four regions. Instead, the site-specific investigation of region 4 reveals quantitative links between seafloor geomorphology/geology and seep occurrence. As reported above, over half of the seeps in region 4 occur on anticlines. This characteristic of region 4 helps to explain the probability depth distribution of seeps within the GHSZ. As shown in Figure 3 the OAH diverges seaward, representing a unique feature of southern Cascadia not observed elsewhere along the margin. According to Tréhu and Phrampus (2022), this region is also marked by the greatest distance between the landward-most BSR pick and 500 m contour (Supplementary Figure S4). Given that the OAH is an actively deforming portion of the margin, younger anticlines are formed from west to east (i.e., seaward or within the OAH) by uplift and horizonal compression as the deformation front steps seaward (Carson et al., 1974; McNeill et al., 1999; McCrory et al., 2002). In this region, the BSR is also tens of km seaward of the landward limit of the GHSZ where accretion of trench sediments forms a thick accretionary wedge (Tréhu and Phrampus, 2022). This spatial pattern highlights the potential role of the OAH on not only seep distribution but hydrate formation, where accretionary wedge tectonics and seafloor morphology serve as structural controls on fluid reservoirs, trapping, and transport. The 2-D conceptual model shown in Figure 6 demonstrates fluid flow paths and seafloor seepage features that are dependent on folds and fault patterns. Whereas convergence rates are similar along the Cascadia margin, the structural style and width of the actively deforming outer wedge varies along-strike, with a narrower wedge in the south relative to wider wedge in the north (Watt and Brothers, 2020). As a result, the narrower wedge forms a steeper accretionary prism with closely spaced faults/anticlines that stack on top of one another (i.e., imbricate faulting). In contrast, the wide wedge in northern Cascadia forms widely spaced anticlines. Region 3 is also characterized by a relatively high occurrence of seeps on anticlines characterized by mixed or dual vergence (Watt and Brothers, 2020), resulting in more-symmetrical fold structures and steeper seafloor gradients along the fold limbs evident in seismic data (Booth-Rea et al., 2008; Phrampus et al., 2017).

The occurrence of seeps on anticlines follows the conceptual models of focused fluid flow and gas migration along the crest of an anticline described by Paull et al. (2008), such as observed offshore north Panama (Reed et al., 1990), Costa Rica (Kluesner et al., 2013), and along the Queen Charlotte Fault (Prouty et al., 2020), where folding, thrust faulting, and sedimentation along a deformation front may facilitate fluid flow and seafloor venting. This relationship has been previously noted on the Cascadia margin, where the development of anticlinal ridges controls the location of seepage and are strongly correlated with vertical fluid migration paths (Carson et al., 1991; Carson et al., 1994; Johnson et al., 2003). Similarly, Tréhu et al. (1999) speculated positive feedback between gas hydrate formation and anticlines. Merle et al. (2021) also noted that deeper seeps (>525 m) were coincident with major compressional anticlinal (and diapiric ridges) within the accretionary prism. The efficiency of anticlines to trap and focus fluid flow is clearly seen in the MSC data and chimney analysis (Figures 4, 5), where broad high-probability zones of fluid migration and areas of chaotic reflections and acoustic wipe-out are located below the anticline hinge line. At depth, disruption of the BSR coincides with the location of a deeper fault, which is likely channeling fluids and causing breakdown of the hydrate cap. Upon intersecting the first fault up-dip, the gas/fluid presumably migrates upward along a path of least resistance to the seafloor. The seismic profile also shows the presence of a relic (“double”) BSR (Figure 5), possibly caused by episodes of uplift leading to shoaling/reestablishment of the BSR due to changing thermal conditions. This analysis also illustrates how faults disrupt the BSR (and presumably hydrate stability) and provide a pathway for fluids and gas to reach the seafloor, resulting in seafloor emissions. Hydraulic fracturing can also serve as a gas conduit. For example, along the southern Hikurangi margin in Aotearoa/New Zealand, anticlines act as a focusing mechanism for the fluids, where there is sufficient upward pressure to overcome the lithostatic load of the overburden and blow-out, or fracture, leading to a vertical pathway to the seafloor (Crutchley et al., 2021).

Taken together, this quantitative approach has identified characteristics/features as well as processes along the margin that facilitate seep emissions that operate independently of gas hydrate dissociation in response to contemporary ocean warming. In doing so, this study provides an expanded interpretive framework for understanding the controls and influencing factors on the depth and spatial distribution of seeps along the Cascadia margin and is consistent with stable gas hydrate occurrences even under future scenarios of warming bottom waters (Ruppel, 2011; Ruppel and Kessler, 2017). This interpretive framework may also be applied to other subduction margins that are dominated by accretion, such as the southern Hikurangi margin offshore Aotearoa/New Zealand (Barnes et al., 2010), where high fluid pressures play an important role in maintaining thrust wedges (Saffer and Bekins, 2002). For example, the majority of seeps on the southern Hikurangi margin were detected on crests of thrust-faulted anticlinal ridges at 700–1,200 m and where a network of near-surface faults and chimneys facilitates fluid expulsion and a breakdown in BSR (Mountjoy et al., 2009a; Barnes et al., 2010; Law et al., 2010; Crutchley et al., 2018; Watson et al., 2020). While seeps along the Hikurangi margin are consistently located above the deforming backstop, seep distribution based on MBES water column backscatter data is not uniform, and similar to the Cascadia margin, no evidence was found to suggest clustering of seeps near the landward limit of the GHSZ (650 m; Supplementary Figure S5). Instead, seeps in the north are concentrated near the shelf break within the Tuaheni Landslide Complex, where a combination of factors have contributed to slope instability, including earthquake ground shaking and gas release (Mountjoy et al., 2009b; Mountjoy et al., 2014; Gross et al., 2018; Carey et al., 2022; Crutchley et al., 2022). In comparison, seeps in the south occupy deeper water depths and are found to cluster along the crests of parallel ridges (Supplementary Figure S6), where gas is trapped beneath large sedimentary folds/anticlines formed by deep-seated thrust faults, and focused fluid flow is facilitated by both faults and hydraulic fractures (Crutchley et al., 2021).

This comparison highlights the role of regional geologic structures influencing clusters of focused fluid flow, as well as detailing margin-wide variability in seep distribution along both margins. For example, the largest depth range of seeps is located in the relatively narrow southern end of both margins where there is a transition from subduction to strike-slip deformation (Barnes et al., 2010; Watt and Brothers, 2020). At this transition, strike-slip faults that trend across the slope (Gulick et al., 2002) may provide additional near-vertical pathways for fluid flow across a range of water depths. The importance of the geometry of geologic structures to influence spatial location and character of seeps has been previously documented (Kulm et al., 1986; Le Pichon et al., 1992). However, results from our study, which demonstrate the importance of both geometry of geologic structures and the role of the margin-wide morphotectonics, may expand previous models of the distribution of fluid expulsion along subduction zones.