Natural marine gas emissions have been reported at numerous regions along ocean continental margins across the globe (Fleischer et al., 2001; Mazurenko and Soloviev, 2003; Phrampus et al., 2020), representing areas where hydrocarbons from the sediments enter the hydrosphere. A better understanding of gas bubble fluxes – consisting mainly of methane – is important to determine the global inventory of carbon in marine sediments (Klauda and Sandler, 2005; Wallmann et al., 2012; Ruppel and Kessler, 2017) and to discern processes related to ocean chemistry and biology (Wallmann et al., 2006). An increasing number of studies focus on the quantification of gas released from the seafloor to the hydrosphere and to the atmosphere. These studies often use short-term observations acquired during research expeditions, which are limited in documenting spatial and temporal variability. In fact, various studies (e.g. Tryon et al., 1999; Boles et al., 2001; Varadharajan and Hemond, 2012; Kannberg et al., 2013) have shown the highly transient nature of gas emissions in a wide range of time scales. Repeated measurements during research expeditions are limited in detecting changes over time, as the speed and dynamic range of observed changes remain unknown. Continuous long-term observations of seep systems are needed to better understand these dynamic environments (Suess et al., 2001; Heeschen et al., 2003). To date, the most advanced and effective way to continuously monitor the seafloor is through permanent seafloor observatories that provide high power and data bandwidth to the deep sea, allowing for the operation of various experiments in addition to optimizing recording parameters in reaction to seafloor events (Barnes et al., 2011). Ocean Networks Canada has been operating cabled observatories since 2006 and has been providing continuous data from the northern Cascadia continental margin since 2009. Using real-time data from NEPTUNE observatory’s widely distributed sensor networks (Barnes et al., 2011), the scientific community is able to investigate the dynamic changes of seep environments along this margin (Scherwath et al., 2019).

The northern Cascadia continental margin offshore Vancouver Island lies along the subduction zone of the Juan de Fuca Plate. The sediments deposited on the incoming oceanic crust are accreted and progressively folded and faulted, forming elongated anticlinal ridges up to 700 m in height (Davis and Hyndman, 1989). The occurrence of seafloor seepage and mud volcanism is typical at accretionary margins, as the degradation of organic matter often leads to high amounts of light hydrocarbons in the sediments within the prism, and fluid flux, driven by buoyancy, may be facilitated by tectonic forces during the accretion process (Kopf, 2002; Zühlsdorff and Spieß, 2004; Judd and Hovland, 2007; Suess, 2010). The Cascadia accretionary prism, however, is seismically locked (Hyndman and Wang, 1993; Obana et al., 2015; McGuire et al., 2018) and almost no seafloor motion is observed (Stone et al., 2018). There are several known areas of seepage along the Cascadia continental margin (Riedel et al., 2018; Scherwath et al., 2019) including the seep area at Clayoquot Slope, investigated in this study. Several studies focused on quantification and variability of fluid fluxes released along this margin (e.g. Heeschen et al., 2003; Kannberg et al., 2013; Hautala et al., 2014; Johnson et al., 2015; Philip et al., 2016; Römer et al., 2016; Riedel et al., 2018; Marcon et al., 2021). However, most of these studies are based on small-scale or short-term observations. Long-term variations of fluid fluxes over seasonal timespans or over decades are still widely unknown.

The Clayoquot Slope is part of the accretionary prism off Vancouver Island, BC, and located at a mid-slope basin in about 1,260 m water depth, and approximately 20 km landward of the deformation front. It hosts several gas emission sites characterized by various concentrations of gas in the sediments. In this work, we study a zone called Gastown Alley, a SW-NE elongated seep zone defined by several aligned flares (Römer et al., 2016). Gastown Alley extends for about 400–500 m from the Bullseye Vent area (a seafloor depression with shallow gas accumulations and hydrate occurrences surrounded by thin carbonate crusts) towards the Bubbly Gulch vents (a slope failure with strong gas ebullition from exposed sedimentary layers) (Riedel et al., 2006a, 2006b; Lapham et al., 2013; Paull et al., 2015). Gas and/or hydrate accumulations at Gastown Alley are located deeper than at Bullseye Vent and occur about 20 mbsf (Riedel et al., 2006a, 2006b; Willoughby et al., 2008; Römer et al., 2016). However, gas venting was detected at Gastown Alley but not at Bullseye Vent, where carbonate and/or hydrate capping of the sediments likely prevents ebullition at the seabed. Methane hydrates are occurring at the Bullseye Vent depression close to the seafloor down to at least 40 mbsf (Riedel et al., 2006a) but were not encountered in the shallow sub-seafloor outside of the depression. The presence of gas hydrates below Gastown Alley has not been confirmed. Gastown Alley lies well within the gas hydrate stability zone and gas hydrates should be expected to be stable. However, the sulfate-methane interface (SMI), above which gas hydrates cannot form (Paull et al., 2005; Bhatnagar et al., 2011; Riedel and Collett, 2017), is located about six to eight mbsf at Gastown Alley (Solem et al., 2002; Riedel et al., 2006b). At Bubbly Gulch, hydrates were encountered within 30 cm of the seafloor right at the seepage zone but were absent in the sediments around it (Paull et al., 2015).

A year-long acoustic monitoring of gas ebullition in an area located close to Bubbly Gulch (Figure 1) between July 2012 and July 2013 showed that decreasing bottom pressures during falling tides facilitated the gas migration in the sediments and the onset of gas ebullition at the seabed (Römer et al., 2016). Three gas ebullition phases, each lasting several months, were observed: phases 1 and 3 were characterized by transient gas release, whereas phase two was characterized by intense ongoing gas release. These three phases were postulated to be caused by variable gas supply rates in the subsurface. The study also found no link between the gas emissions variability (both short-term variations and activity phases) and the seismicity, the wave height variations or the seasonal oceanographic variability.

Using long-term acoustic monitoring we expand on previous work and now offer a 4 years long timeseries of data to address important questions on the variability of gas venting. In particular, we look into the role of tides on the timing and intensity of gas emissions, investigate if there are seasonal variations in bubble release over the 4 years of observations, or if other drivers can also be identified that explain the observed variability.

Flares were detected in 28,430 out of 33,650 scans, that is, at least 84% of the monitoring period (22 June 2017 to 30 May 2021). Times of inactivity of gas emissions (hereafter referred to as “inactive times”) were non-uniformly distributed throughout the monitoring period. Some months contained numerous inactive times – up to 43 stop-and-start events in March and June 2018, accounting for almost 50% of inactivity – whereas other months had no inactive time at all (Figure 2; Table 1).

Frequency analyses of the mean backscatter of detected flares show that the timeseries contains diurnal and semi-diurnal constituents corresponding to those of the local mixed tide regime (Figure 2), suggesting a tidal influence on the gas emissions. The local mixed tide is defined by two tidal cycles of unequal amplitude per day, which result from the sum of both diurnal and semi-diurnal constituents. The power of the tidal frequencies in the sonar data varies throughout the timeseries and appears to be strongest during periods when the venting stops and starts intermittently (Figure 2). Frequency analyses of selected timeseries segments, one with discontinuous venting (June 2018) and two with continuous venting (August 2017 and April 2020), showed that the semi-diurnal and diurnal tidal harmonic constituents were clearly present in the June 2018 timeseries but weak or even absent in the August 2017 and April 2020 timeseries (Figure 3). This implies that tidal cycles control the alternation of active/inactive times but to a lesser degree the variation in intensity of the gas release during active venting phases. This is confirmed by the power spectral density of the binary timeseries of active and inactive times (0: inactive, 1: active), which also reveals distinct peaks at the dominant tidal frequencies (Figure 4). This result clearly demonstrates that tides influence the timing of the venting activity in the monitoring area. However, fluctuations of backscatter magnitude during active venting phases correlate neither with the bottom pressure data [Pearson correlation coefficient r = -0.01, p-value = 0.1, n = 28,402, and 95% confidence interval (-0.02, 0.00)] nor with the rate of bottom pressure change [Pearson correlation coefficient r = 0.01, p-value = 0.04, n = 28,402, and 95% confidence interval (0.00, 0.02)], another indication that the vigor of active gas emissions may not be influenced by bottom pressure variations.

The bottom pressure is largely controlled by the mixed local tide and averaged around 1,284.16 dbar (median = 1,284.21 dbar) during the monitoring period (Supplementary Figure S1). Each venting activity phase (start, stop, active, inactive) was observed at all tidal phases and was not uniquely restricted to any particular tidal phase (Table 2). Venting start times were observed predominantly during falling tides (64.7%) (Figure 5A) and at pressures higher than the mean pressure (61.7%) (Figure 5B). By contrast, venting stop times were almost equally distributed around the mean pressure but significantly more common (72.6% vs. 27.4%) during rising tides compared to falling tides. Active venting times appeared to be independent from both absolute bottom pressure and sign of pressure change, whereas inactive venting was slightly more common at pressures above average (61.8%) and during rising tides (55.8%). The amplitude of the falling tides (ebb amplitude: min = 0.24 m, max = 4.17 m, mean = 2.03 m, standard deviation = 0.84 m) and of the rising tides (flow amplitude: min = 0.62 m, max = 3.47 m, mean = 2.03 m, standard deviation = 0.58 m) did not appear to influence the activation and cessation of gas emissions (Figure 5C).

The distributions of bottom current velocities in the lowest ADCP depth-bin, corresponding to water depths from 1,234 m to 1,242 m, are similar during active and inactive times, and during ebb and flow tides (Figure 6). This suggests that bottom current velocities do not affect the venting activity.

Bottom current directions in the lowest ADCP bin show a different distribution pattern depending on whether gas emissions are active or inactive (Figure 7, Figure 8 and Supplementary Figure S2). The bottom current directions during active times are strongly correlated to the bottom current directions of ebb tides [Pearson correlation coefficient r = 0.96, p-value = 2.4 × 10⁻²⁰, n = 36, and 95% confidence interval (0.92, 0.98)] and the bottom current directions during inactive times are correlated to the bottom current directions of flow tides [Pearson correlation coefficient r = 0.91, p-value = 1.4 × 10⁻¹⁴, n = 36, and 95% confidence interval (0.83, 0.95)]. By contrast, the current directions during active times are less correlated to the bottom current directions of flow tides [Pearson correlation coefficient r = 0.69, p-value = 3.1 × 10⁻⁶, n = 36, and 95% confidence interval (0.47, 0.83)] and the current directions during inactive times show no correlation to the current directions during ebb tides [Pearson correlation coefficient r = 0.22, p-value = 0.2, n = 36, and 95% confidence interval (−0.12, 0.51)]. Furthermore, the distribution of current directions during inactive times shows a higher proportion of westward (flowing towards west) bottom currents than the distribution of current directions during flow tides (Figure 7): almost 66% of inactive gas emissions co-occurred with bottom current flow directions ranging from 210° to 320° relative to North. By contrast, dominantly westward bottom currents within the same range of directions occurred during only half (51%) of the flow tidal phases. In other words, the venting activity pauses preferentially when bottom current flow towards southwest and west, not only during flow tides but also during ebb tides. This suggests that there is a connection, with or without causality, between the current direction and the cessation of venting.

Bottom water temperatures between June 2017 and May 2021 (Figure 9) ranged from 2.64 to 3.15°C (mean: 2.90°C, standard deviation: 0.067°C). A clear semi-diurnal periodicity shows that the intraday temperature variations are related to the tidally controlled bottom currents (Figure 8 and Supplementary Figure S2). Amplitudes of daily variations ranged from 0.03 to 0.31°C (mean: 0.14°C, standard deviation: 0.039°C). Longer multiday temperature variations are also present (Figure 9, Supplementary Figure S1, Supplementary Figure S2), which could not be related to the tidal cycles and are likely linked to longer-term changes of the bottom water masses possibly caused by seasonal variations or variability in the regional deep-ocean circulation. Fluctuations of backscatter magnitude during active venting phases do not correlate with the temperature variations of bottom waters [Pearson correlation coefficient r = 0.07, p-value = 3.5 × 10⁻³⁴, n = 28,418, and 95% confidence interval (0.06, 0.08)], an indication that the vigor of active gas emissions may not be influenced by temperature variations. Venting stop times and inactive times were associated with bottom temperatures 0.01 and 0.03°C above average (Figure 5D). Median temperatures were 2.90°C during start times and active times, 2.92°C during stop times and 2.91°C during inactive times. However, high temperatures were not consistently associated with inactive times indicating that temperature alone is unlikely to be responsible for pausing gas emissions.

Storms occurred every year mainly between November and April (Figure 10). The data indicate that the November-to-April months do not show a pattern of co-occurrence either with periods of elevated/reduced venting or with intervals of continuous/discontinuous venting.

The number of detected flares varied during the monitoring period and ranged from 0 to 8 with a median value of 2 (Figure 11). The outlets of the bubble release on the seabed cannot be located precisely with this sonar but the spatial distribution of their estimated positions reveal two main areas of recurrent seepage around the sonar (Figure 12). The southern cluster is active most often, but it is widespread and likely coalesces several smaller clusters that could not be resolved due to the uncertainties in the flare basepoint estimations. The northeastern cluster is smaller and less continuously active than the southern cluster.

Venting within our monitoring area is active more than 84% of the time but affected by recurrent interruptions related to local tides. It is evident that the onset and the cessation of bubble emissions are largely controlled by the diurnal and semi-diurnal frequencies of the local mixed tidal regime. While vent activation is more common during falling tides, cessation is observed predominantly during rising tides. This result corroborates previous results from similar monitoring efforts in deep sea settings near Bubbly Gulch in 2012–2013 and at Southern Hydrate Ridge in 2018, where bubble emission activation was found to peak during falling tides (Römer et al., 2016; Marcon et al., 2021).

Interestingly, the no-venting phases are not regularly distributed through the monitoring period and there are numerous week- and months-long intervals of continuous venting during which no venting interruptions were recorded. It is unclear why the bubble release is at times continuously active and does not stop and start with the tidal cycles. For instance, our results showed that the distribution of active and inactive times is unaffected by the amplitude of the pressure variations related to the spring and neap tidal cycles. Furthermore, seismic tremors or seasonal storms were not found to be related to the activity of gas emissions on the Clayoquot Slope (Lapham et al., 2013; Römer et al., 2016) and at Southern Hydrate Ridge (Marcon et al., 2021). The accretionary prism of the Cascadia Subduction Zone is fully locked (Hyndman and Wang, 1993; Obana et al., 2015), especially in our study area (McGuire et al., 2018), and local shaking levels are practically nonexistent (Scherwath et al., 2011; Stone et al., 2018). Hence, seismicity is unlikely to be responsible for the frequent alternation of transient and continuous ebullition phases that we observed. In this study, we also did not observe any clear link between the continuity or intermittence of the gas emissions and the seasonal winter storms (Figure 10). Accumulations of gas and possibly hydrates are known to occur at Gastown Alley at about 20 mbsf (Römer et al., 2016). One hypothesis is that the fluid supply at depth is transient and not always sufficiently high to ensure continuous bubble release at the sediment surface independently to the bottom pressure variations. Such a mechanism was postulated by Römer et al. (2016) to explain the observed activity phases of the gas ebullition. Unfortunately, we cannot test this hypothesis with our current results, as information about subsurface fluid supply over time is not available. Although the timeseries may show small decreases in venting intensity prior to some intervals of intermittent venting, this is not consistent throughout the timeseries (Figure 2). Another possibility could be that the bubble release in the Gastown Alley area is in fact continuous but so spatially variable that it does not always occur within the range of the sonar. This is a reasonable assumption given that the Gastown Alley is an elongated seep area about 400–500 m long (Römer et al., 2016), large parts of which lie outside of the monitoring area of the rotary sonar (Figure 1). Continuous but spatially variable venting was observed at the Southern Hydrate Ridge summit using an overview sonar that monitored the whole active area of the summit (Marcon et al., 2019). There, single vents were found to pause frequently but rarely, if ever, all simultaneously (Marcon et al., 2021). Similarly, the number of active flares within our monitoring area varies considerably over time (Figure 11) and is non-zero more than 84% of the time. Hence, it is likely that the venting over the entire active area of Gastown Alley rarely pauses completely. At Southern Hydrate Ridge, hydrate formation in the shallow subsurface is considered to be responsible for the short-term spatial and temporal variability of the bubble release (Marcon et al., 2021). Although Gastown Alley lies deep within the gas hydrate stability zone, the presence or absence of gas hydrates in the shallow sediments at Gastown Alley was never confirmed. However, the sulfate-methane interface (SMI) at Gastown Alley is located about 6–8 mbsf (Solem et al., 2002; Riedel et al., 2006b), and gas hydrates can only occur below the SMI and possibly much deeper, where methane exists above solubility (Paull et al., 2005; Bhatnagar et al., 2011; Riedel and Collett, 2017). Riedel and Collett (2017) compared data from 58 drill sites from ten different geological regions, including in the northern Cascadia Margin, and found that in the absence of advection, the top of the gas hydrate occurrence zone was always deeper than 30 mbsf. Potentially, hydrate formation in the pore and fracture spaces below this depth could block pathways for gas migration and cause individual gas emissions to temporarily pause.

The current understanding is that low bottom pressures associated with falling tides facilitate the activation of bubble emissions by reducing the total stress on the sediments (Tryon et al., 1999, 2002; Boles et al., 2001; Leifer and Boles, 2005; Liu and Flemings, 2009; Scandella et al., 2011). However, the vent activation times at Gastown Alley are more common at higher than average bottom pressures. Hence, it seems to be the pressure decrease rather than the absolute pressure that helps trigger ebullition. Because of fracture resistance, a pressure drop can cause a temporary dynamic disequilibrium between the pore pressure within gas-rich sediments and the hydrostatic pressure, favoring the advection of gas (Leifer and Boles, 2005; Sultan et al., 2020). Römer et al. (2016) found that bubble emission activation is more common at low pressures close to the low tide turning point. This difference with our results indicates that activation of gas emissions at Gastown Alley is sensitive to smaller pressure decreases than in the 2012-2013 study area. This is interesting since the study area from Römer et al. (2016) is located only 530 m from our monitoring area in the direction of Bubbly Gulch.

It is likely that gas fluxes are higher at Gastown Alley than in the 2012-2013 study area. This is supported by the findings from sub-bottom profiles, which show that the sedimentary layers below Gastown Alley contain large gas accumulations at about 20 mbsf (Römer et al., 2016). By contrast, gas occurrence in the sediments below the study area of the 2012-2013 monitoring is more diffuse and likely related to lower and more transient upward fluxes. Furthermore, flux estimations derived from acoustic flares detected with a calibrated vessel-mounted splitbeam echosounder indicated that instantaneous fluxes at Gastown Alley (min = 0.1 L/min, mean = 4.4 L/min, max = 17.9 L/min) were generally stronger than near the 2012-2013 study area (min = 0.3 L/min, mean = 2.7 L/min, max = 4.7 L/min) (Riedel et al., 2018). Although acoustic quantification of gas emissions is subject to large uncertainties, the flux estimations from Gastown Alley were the largest methane emission fluxes calculated along the Cascadia Margin (Riedel et al., 2018).

Inactive times at Gastown Alley are also strongly associated with westward bottom currents. This association is partly explained by the fact that these are the most represented current directions during rising tides (Figure 7), when inactive times are more common. However, westward bottom currents are in higher proportion during inactive times than they are during flow tides, which raises the question if and how bottom current directions could also play a role in the activation or cessation of gas emissions. Seabed morphology can generate high form drag at the sediment surface (Nash and Moum, 2001) leading to pressure and velocity fluctuations within the sediment pore system (Higashino et al., 2009). Furthermore, bottom currents are known to facilitate ebullition onsets in shallow lakes by causing shear stress at the sediment surface (Joyce and Jewell, 2003). Our results do not show any indication that the activity of gas emissions is influenced by current velocities. However, given the particular mid-slope location of Gastown Alley we cannot exclude that certain current directions might cause increased or reduced form drag that would affect the pore pressure and, thus, the release of gas. Another possibility is that current-induced temperature variations might affect the methane solubility. For instance, a temperature decrease could be expected to pause ebullition by increasing the solubility of methane in the bottom water (Riedel et al., 2021). Our results do not support this possibility since positive and negative temperature variations did not consistently coincide with activation and cessation of venting (Supplementary Figure S1). Also, we observed that inactive times were not more frequent at lower temperatures, which seems to rule out that temperature-driven solubility variations cause the onset or cessation of gas emissions. At Gastown Alley, a temperature decrease of 0.14°C (mean amplitude of daily variations) would result in a methane solubility increase of 2 × 10⁻⁴ mol of methane per kilogram of seawater (Duan and Mao, 2006; Kossel et al., 2013). This corresponds to 3.6 µl of methane (at the seabed conditions of pressure and temperature). Assuming that the volume of pore water that is affected by the temperature increase is 1 m³ (a large overestimation used for the sake of simplification), the volume of extra dissolved methane would reach 3.6 ml, which corresponds to about 15 bubbles of a typical bubble-size distribution (Römer et al., 2016). Hence, we consider it unlikely that the daily temperature variations alone would cause the gas ebullition to pause. Furthermore, pore water in the shallow sediments above the SMI is unlikely to be close to saturation.

The flare source distribution showed at least two main clusters within the monitoring area, possibly reflecting distinct pathways for gas migration within the sediments (Figure 1, Figure 12). The AUV-based bathymetry (Figure 1 and Paull et al., 2015) does not reveal any obvious local topographical features around those vents that might be associated with the bubble release, as is known from other seep areas both in the immediate vicinity of Gastown Alley and farther along the Cascadia margin. The Bullseye Vent area located 250 m SW of the sonar location is characterized by large 5–6 m deep depressions surrounded by thin carbonate crusts, whereas Bubbly Gulch (650 m NE of the sonar location) coincides with a slope failure and hosts hydrate mounds (Paull et al., 2015). Hydrate Ridge is associated with topographic highs, hydrate mounds and depressions shaped by focused methane venting (Suess et al., 2001; Tréhu et al., 2004; Bangs et al., 2011; Philip et al., 2016; Marcon et al., 2021). By contrast, the seabed within the sonar monitoring area at Gastown Alley is gently sloping towards the southeast (Figure 1) and it is covered with soft sediments from which bubble release occurs (Link to seabed footage: https://data.oceannetworks.ca/SeaTube?resourceTypeId=1000&resourceId=1001&diveId=410&time=2014-05-24T07:20:10.000Z). White microbial mats occur within the monitoring area but not where the bubble outlets were observed. The bubble release does not visibly alter the local shape of the seabed around the outlets. This might be due to the low-permeability fine silty and clayey sediments (Riedel et al., 2006a) and the absence of gas hydrates in the shallow sediments, which may favor a more direct transfer of gas through fractures. Alternatively, this seepage area may have formed as a result of the sealing of Bullseye Vent with carbonates or hydrates (Riedel et al., 2002, 2006b; Römer et al., 2016) and may represent a shift of the gas migration from Bullseye Vent towards the northeast. This would be consistent with previous work, which showed that new fractures had formed at Bullseye Vent between 2000 and 2005 (Riedel, 2007), and that the presence of extensive microbial mats and the absence of carbonate crusts at the sediment surface in the area northeast of Bullseye Vent indicate venting in this area may be more recent (Furlong, 2013). With time, the venting activity might also shape the seabed at Gastown Alley. However, more data about the sub-seabed plumbing network at Gastown Alley are needed to investigate these hypotheses further.

The 4 years monitoring of gas emissions at Gastown Alley showed clearly that tides influence the timing of gas emissions start and stop times. However, it could not confirm that the vigor of active gas emissions is modulated by tides. Furthermore, the occurrence of month-long intervals of uninterrupted venting within the monitoring area suggests that tidally-induced bottom pressure variations alone were not sufficient to pause the venting. We cannot exclude that external factors (e.g. bottom currents, temperature, transient gas supply) might influence the onset and cessation of bubble plumes, but it is likely that the venting at Gastown Alley during the 4 years monitoring period has been sufficiently intense to remain mostly active independently of the tidal cycles, but not always occurring within range of the sonar.