Autonomous methane seep site monitoring offshore western Svalbard: hourly to seasonal variability and associated oceanographic parameters

Abstract

Improved quantification techniques of natural sources are needed to explain variations in atmospheric methane. In polar regions, high uncertainties in current estimates of methane release from the seabed remain. We present unique 10- and 3-month time series of bottom water measurements of physical and chemical parameters from two autonomous ocean observatories deployed at separate intense seabed methane seep sites (91 and 246 m depth) offshore western Svalbard from 2015 to 2016. Results show high short-term (100–1000 nmol L${}^{-\mathrm{1}}$ within hours) and seasonal variation, as well as higher (2–7 times) methane concentrations compared to previous measurements. Rapid variability is explained by uneven distribution of seepage and changing ocean current directions. No overt influence of tidal hydrostatic pressure or water temperature variations on methane concentration was observed, but an observed negative correlation with temperature at the 246 m site fits with hypothesized seasonal blocking of lateral methane pathways in the sediments. Negative correlation between bottom water methane concentration (and variability) and wind forcing, concomitant with signs of weaker water column stratification, indicates increased potential for methane release to the atmosphere in fall and winter. We present new information about short- and long-term methane variability and provide a preliminary constraint on the uncertainties that arise in methane inventory estimates from this variability.

1 Introduction

Unexplained changes in atmospheric methane (CH${}_{\mathrm{4}}$) mole fraction motivate research in understanding and quantifying non-anthropogenic sources . The atmospheric forcing of CH${}_{\mathrm{4}}$ is particularly sensitive to changes in emission rates due to its high warming potential and short lifetime. Improved knowledge about atmospheric CH${}_{\mathrm{4}}$ fluxes is therefore crucial to constrain future climate projections . These properties of atmospheric CH${}_{\mathrm{4}}$ also make reducing anthropogenic CH${}_{\mathrm{4}}$ emissions a potential solution for rapid climate change mitigation . A global effort to cut greenhouse gas emissions through international agreements is, however, dependent on precise estimates of sources and sinks to verify contributions from different nations.

Seabed seepage is considered a minor source of atmospheric CH${}_{\mathrm{4}}$, but there is high uncertainty in current and predicted emission estimates . Current estimates suggest a total contribution of 7 (5–10) Tg yr${}^{-\mathrm{1}}$ , which is $\sim$ 1 % of the total CH${}_{\mathrm{4}}$ emissions to the atmosphere. Methane is released from the seabed as free gas (bubbles) and dissolved gas in sediment pore water. Bubbles rise quickly towards the sea surface, but most CH${}_{\mathrm{4}}$ dissolves near the seafloor because of gas exchange across the bubble rims and bubble dissolution . Dissolved CH${}_{\mathrm{4}}$ is dispersed and advected by ocean currents and is continuously transformed to carbon dioxide (CO${}_{\mathrm{2}}$) by bacterial aerobic oxidation . These processes significantly limit the lifetime of CH${}_{\mathrm{4}}$ in the water column, and the amount of CH${}_{\mathrm{4}}$ that can reach the atmosphere is highly dependent on the depth where the seepage occurs . Intense CH${}_{\mathrm{4}}$ seepage at shallow depths in coastal areas and on continental shelves is therefore the main potential source of seabed CH${}_{\mathrm{4}}$ to the atmosphere.

The shallow continental margins of the Arctic Ocean store large amounts of CH${}_{\mathrm{4}}$ as free gas, gas dissolved in pore water fluid, and gas hydrates , i.e., clathrate structures composed of water trapped by hydrocarbon molecules formed and kept stable at low temperature and high pressure . Increasing bottom water temperature has the potential to liberate methane from these reservoirs via various mechanisms, potentially resulting in a positive climate feedback loop .

Studies on CH${}_{\mathrm{4}}$ inventory, distribution and release in the Arctic Ocean are mainly based on research cruise data from late spring to early fall, when ice and weather conditions allow fieldwork in the region , whereas winter data are sparse. Bottom water temperature , water mass origins , micro-seismicity and hydrostatic pressure have all been proposed to be linked with sources and sinks of CH${}_{\mathrm{4}}$ in the water column. These processes act on a wide range of timescales, from hours (e.g., hydrostatic pressure) to decades (bottom water temperature). Without a better understanding of the spatial and temporal variability of CH${}_{\mathrm{4}}$ in Arctic seep sites, it is challenging to untangle these processes. Unconstrained local variability in CH${}_{\mathrm{4}}$ seepage and concentration also imposes a high degree of uncertainty on CH${}_{\mathrm{4}}$ inventory estimates . The combination of climate-sensitive CH${}_{\mathrm{4}}$ storages, vast shallow ocean regions and limited data availability highlights the need for more understanding of seabed CH${}_{\mathrm{4}}$ seepage on Arctic shelves.

To assess the aforementioned challenges, we have obtained, analyzed and compared two unique long-term underwater multi-parameter time series from seafloor observatories deployed at two distinct intense CH${}_{\mathrm{4}}$ seep sites on the western Svalbard continental shelf (Fig. ) where no CH${}_{\mathrm{4}}$ measurements have previously been done in winter season. We combine high-frequency physical (ocean currents, temperature, salinity, pressure) and chemical (O${}_{\mathrm{2}}$, CO${}_{\mathrm{2}}$, CH${}_{\mathrm{4}}$) data to perform hypothesis testing and provide new insights on CH${}_{\mathrm{4}}$ distribution, content, and variability on short (minutes) and long (seasonal) timescales and potential implications.

1.1 Regional settings

Two observatories (O${}_{\mathrm{91}}$ and O${}_{\mathrm{246}}$) were deployed from June 2015 (CAGE 15-3 cruise) to May 2016 (CAGE 16-4 cruise) from R/V Helmer Hanssen at the inter-trough shelf region between Isfjorden and Kongsfjorden, west of Prins Karls Forland. The O${}_{\mathrm{91}}$ observatory was deployed at 91 m water depth on the continental shelf (78.561${}^{\circ }$ N, 10.142${}^{\circ }$ E), and the O${}_{\mathrm{246}}$ observatory was deployed at 246 m water depth further offshore close to the shelf break (78.655${}^{\circ }$ N, 9.433${}^{\circ }$ E, Fig. ).

Figure 1

Bathymetry of the study area with the locations of the observatories O${}_{\mathrm{91}}$ and O${}_{\mathrm{246}}$ offshore western Svalbard. Flares detected by single-beam echo sounder survey prior to recovering the observatories (May 2016, cruise CAGE 16-4) are indicated with red dots, and ship tracks are shown as brown lines. The inset map shows the working area (red square) offshore Svalbard. WSC and CC refer to the warm West Spitsbergen Current and cold Coastal Current, respectively.

[Figure omitted. See PDF]

Both sites were located in areas with thousands of previously mapped CH${}_{\mathrm{4}}$ gas seeps (e.g., ; this work; see Fig. ), often referred to as “flares” due to the appearance of bubble streams in echo sounder data. Nonetheless, atmospheric sampling in this region suggests that any emissions to the atmosphere are small . Gas accumulation at the O${}_{\mathrm{246}}$ seep site has been suggested to be a result of gas migration in permeable layers within the seabed from deeper free gas or hydrate reservoirs , while seepage at site O${}_{\mathrm{91}}$ has been attributed to thawing sub-sea permafrost due to ice sheet retreat at the end of the last glaciation . Water sampling has indicated high temporal variability, with bottom water concentrations (average) changing from 200 nmol L${}^{-\mathrm{1}}$ within 1 week in July 2014 at O${}_{\mathrm{91}}$ and $\sim$ 80 nmol L${}^{-\mathrm{1}}$ within 20 h (two single-point measurements) at O${}_{\mathrm{246}}$ in August 2010 . A consistent pattern of decreasing concentrations from the seafloor to the sea surface at both sites (400 to $<$ 8 nmol L${}^{-\mathrm{1}}$ at O${}_{\mathrm{91}}$, , and from $>$ 500 to $<$ 20 nmol L${}^{-\mathrm{1}}$ at O${}_{\mathrm{246}}$, ) has also been observed. Further offshore, continuous measurements from a towed fast-response underwater laser spectrometer also revealed very high spatial CH${}_{\mathrm{4}}$ variability .

The local water masses are characterized by exchange and convergence of warm, saline Atlantic water (e.g., temperature $T>\mathrm{3}$ ${}^{\circ }$C and salinity $S_{\mathrm{A}}>\mathrm{34.9}$; ) in the West Spitsbergen Current and colder, fresher Arctic water (e.g., $T<\mathrm{0}$ ${}^{\circ }$C, $\mathrm{34.3}<S_{\mathrm{A}}<\mathrm{34.8}$; ) in the Coastal Current combined with seasonal cooling, ice formation and freshwater input from land (Fig. ). Local mixing rates can be strongly affected by synoptic-scale weather systems, causing upwelling and disruption of the front between the two ocean currents . Freshwater input in summer stratifies the water column, while cooling, storm activity and sea ice formation can facilitate vertical mixing in winter .

2 Methods

The “K-Lander” ocean observatories were designed to monitor CH${}_{\mathrm{4}}$ release and associated physical and chemical parameters in challenging environments (see Appendix ). A launcher equipped with camera and telemetry allowed for safe deployment at a site selected by visual control. Observatory O${}_{\mathrm{91}}$ recorded data from 2 July 2015 to 6 May 2016, while O${}_{\mathrm{246}}$ recorded data from 1 July until 3 October 2015, when data recording ceased due to an electrical malfunction.

Both observatories were equipped with an acoustic Doppler current profiler (ADCP), a conductivity–temperature–depth (CTD) profiler with an oxygen optode, and Contros HydroC CO${}_{\mathrm{2}}$ II and HydroC Plus CH${}_{\mathrm{4}}$ sensors (Fig. a; details about the instrumentation are provided in Appendix ). The deployed HydroC CH${}_{\mathrm{4}}$, being a younger iteration of the sensor, relies on a tunable diode laser absorption spectrometry (TDLAS) detector (rather than non-dispersive infrared spectrometry (NDIR)), while the CO${}_{\mathrm{2}}$ sensors use NDIR detectors. Both sensors were equipped with polydimethylsiloxane (PDMS) membranes and Seabird SBE 5M pumps (see Appendix ).

Table 1

Correlation coefficients between variables at O${}_{\mathrm{91}}$. The terms RTC CH${}_{\mathrm{4}}$ and raw CH${}_{\mathrm{4}}$ refer to response-time-corrected and untreated CH${}_{\mathrm{4}}$ data, respectively (Sect.  and Appendix ).

	RTC CH${}_{\mathrm{4}}$	Raw CH${}_{\mathrm{4}}$	Temperature	Salinity	Oxygen	Pressure	Solubility	Wind speed	CO${}_{\mathrm{2}}$
	mol L${}^{-\mathrm{1}}$	mol L${}^{-\mathrm{1}}$	${}^{\circ }$C		mol L${}^{-\mathrm{1}}$	dbar	mol L${}^{-\mathrm{1}}$	m s${}^{-\mathrm{1}}$	$\mathrm{\mu }$atm
RTC CH${}_{\mathrm{4}}$	1	0.91	$-$0.06	0.23	0.03	0.08	0.06	$-$0.33	$-$0.25
Raw CH${}_{\mathrm{4}}$	0.91	1	$-$0.07	0.27	0.03	0.10	0.06	$-$0.37	$-$0.31
Temperature	$-$0.06	$-$0.07	1	0.69	$-$0.94	$-$0.01	$-$0.99	0.37	0.29
Salinity	0.23	0.27	0.69	1	$-$0.78	$-$0.06	$-$0.58	0.06	0.46
Oxygen	0.03	0.03	$-$0.94	$-$0.78	1	0.02	0.85	$-$0.33	$-$0.67
Pressure	0.08	0.10	$-$0.01	$-$0.06	0.02	1	0.16	0.00	$-$0.10
Solubility (CH${}_{\mathrm{4}}$)	0.06	0.06	$-$0.99	$-$0.58	0.85	0.16	1	$-$0.35	$-$0.30
Wind speed	-0.33	$-$0.37	0.37	0.06	$-$0.33	0.00	$-$0.35	1	0.52
CO${}_{\mathrm{2}}$	$-$0.25	$-$0.31	0.29	0.46	$-$0.67	$-$0.10	$-$0.30	0.52	1

The high power consumption of the Contros HydroC CH${}_{\mathrm{4}}$ and CO${}_{\mathrm{2}}$ sensors required a power cycling mode to allow for long-term monitoring while simultaneously capturing rapid short-term variability. Partial pressure of CH${}_{\mathrm{4}}$ and CO${}_{\mathrm{2}}$ was therefore measured continuously for 24 h every 21 d and for 1 h every day (see Table ). Methane concentration data were corrected for slow response time onto a 3 min interval grid and converted to absolute concentration following , which is the default “CH${}_{\mathrm{4}}$ concentration” discussed and described in this text (see Appendix ). Faulty pumps in the CO${}_{\mathrm{2}}$ sensors ambiguously increased the response time, which prevented response time correction, making CO${}_{\mathrm{2}}$ data suitable only for long-term qualitative analysis.

Uncertainty ranges for the CH${}_{\mathrm{4}}$ sensor data are reported as 95 % confidence intervals and typically vary between 5 % and 20 % (Fig. b). We did not perform any post-validation and/or intermittent validation. Although always an advantage for all sensors in long-term deployments, such validation is not a requirement for the TDLAS-based sensor (as opposed to NDIR), due to its high long-term stability. Standard post-processing (e.g., inspection of metadata such as internal pressure and temperature) and evaluation of fit residuals in the response time correction procedure (see Appendix B and Dølven et al., 2021) also indicated consistent sensor behavior throughout the deployments. It is also worth noting that the current paper is concerned with large changes and high concentrations, and we are confident that the quality of the response-time-corrected Contros HydroC CH${}_{\mathrm{4}}$ data is sufficient to support the inferences described herein.

We calculated correlation coefficient ($R$) matrices to give a first-order overview of the linear relationships between the measured parameters. We mapped the flares in the area using single-beam echo sounder data collected during the observatory recovery cruise in 2016 (CAGE 16-4, Fig. ) and estimated gas flow rates using the FlareHunter software . Additionally, we obtained 10 m wind reanalysis data from the ERA-Interim database.

We calculated seawater density and CH${}_{\mathrm{4}}$ solubility using the CTD data. A CTD cast (SBE plus 24 Hz) prior to the O${}_{\mathrm{91}}$ recovery (6 May 2016) showed a salinity drift in the conductivity sensor of around $-$0.4 (here and elsewhere in the paper, salinity values are practical salinity). Post-calibration inspection of the conductivity signal and potential water mass mixing end-members indicates that this might have been caused by mud pollution occurring in late 2015 or early 2016.

3 Results

3.1

Time series at site O${}_{\mathrm{91}}$

Dissolved CH${}_{\mathrm{4}}$ concentration at site O${}_{\mathrm{91}}$ ranged from 5 $\pm$ 3 nmol L${}^{-\mathrm{1}}$ (6 December in 2015) to 1748 $\pm$ 142 nmol L${}^{-\mathrm{1}}$ (20 August in 2015) (Fig. a and Appendix ), with 2.5 and 97.5 percentiles of 16 and 785 nmol L${}^{-\mathrm{1}}$. The data follow a nearly lognormal distribution, with a mean and median of 227 and 165 nmol L${}^{-\mathrm{1}}$, respectively, and interquartile range of 88–334 nmol L${}^{-\mathrm{1}}$. Large variations ($>$ 100 up to almost 1000 nmol L${}^{-\mathrm{1}}$) in CH${}_{\mathrm{4}}$ concentration occurred on short timescales ($<$ 1 h) throughout the measurement period (see Fig. a, d and all 24 h periods in Appendix ), with an average range for all 24 h periods of 840 nmol L${}^{-\mathrm{1}}$ and a median rate of change (ROC) of 3.2 nmol L${}^{-\mathrm{1}}$ min${}^{-\mathrm{1}}$. We also observe a long-term trend of decreasing running median (2-week window) concentrations towards winter, from 495 nmol L${}^{-\mathrm{1}}$ in July–August 2015 to 53 nmol L${}^{-\mathrm{1}}$ in January 2016 (Fig. ). There was a relatively weak but significant negative correlation between the wind speed and CH${}_{\mathrm{4}}$ concentration ($R_{\mathrm{RTC}}=-\mathrm{0.33}$) but otherwise weak to non-existent linear relationships between CH${}_{\mathrm{4}}$ concentration and the measured ocean parameters (Table ).

CO${}_{\mathrm{2}}$ averaged 403 $\mathrm{\mu }$atm with an increase towards mid-November 2015 ($\sim$ 410 $\mathrm{\mu }$atm) and then a decrease until 6 May ($\sim$ 391 $\mathrm{\mu }$atm) in 2016 (Fig. a). CO${}_{\mathrm{2}}$ dropped to $\sim$ 305 $\mathrm{\mu }$atm on 24 August, concurrent with a rapid decrease in salinity ($-$0.5), increase in temperature and oxygen, and high CH${}_{\mathrm{4}}$ concentration. The increase in oxygen rules out methanogenesis. Instead, there might be at least two explanations for the reduction of CO${}_{\mathrm{2}}$ and enrichment of CH${}_{\mathrm{4}}$: (i) water column mixing, which brings oxygen-rich, warm and fresh surface water to deeper depths (and with it CO${}_{\mathrm{2}}$-depleted water), or (ii) methane enrichment by zooplankton following the summer bloom.

Figure 2

Time series from (a) O${}_{\mathrm{91}}$ and (b) O${}_{\mathrm{246}}$ showing response-time-corrected (see Appendix ) CH${}_{\mathrm{4}}$, CO${}_{\mathrm{2}}$, temperature, salinity, pressure, CH${}_{\mathrm{4}}$ solubility, oxygen and wind speed (10 m) data. The O${}_{\mathrm{246}}$ data are truncated due to an electrical malfunction in the system on 3 October. (c) A 2 d histogram and 1 SD variance ellipse of bottom current velocity (81 m depth at O${}_{\mathrm{91}}$ and 236 m depth at O${}_{\mathrm{246}}$) and (d) an example of 24 and 1 h (2 and 3 August) CH${}_{\mathrm{4}}$ concentration measurement period from O${}_{\mathrm{91}}$ (green box). All 24 h measurement periods are shown in Appendix . Note the different scales between O${}_{\mathrm{91}}$ and O${}_{\mathrm{246}}$.

[Figure omitted. See PDF]

Bottom water temperature increased steadily from $\sim$ 3 in July 2015 to $\sim$ 5.5 ${}^{\circ }$C in October–November 2015, with occasional sharp shifts ($T\pm \mathrm{1}$ ${}^{\circ }$C) occurring within hours to days (Fig. b). Temperature then decreased from the beginning of December 2015 to $\sim$ 1.8 ${}^{\circ }$C at the end of the deployment in May 2016, showing more frequent and stronger episodes of rapid temperature shifts ($T\pm \mathrm{2}$ ${}^{\circ }$C also occurring over hours or days). Despite uncertainty in salinity data, it is worth noting that these rapid shifts in temperature and salinity were reproduced by the Svalbard 800 model in the same area by eddy activity.

Hydrostatic pressure was mostly governed by tides (94.5 % of variance) with a dominant semi-diurnal M2 tide (M2 refers to a tidal constituent with period 12.42 h; see, e.g., ). Amplitudes varied from $\sim$ 1.2 to 1.5 m during neap and spring cycles (Fig. c).

The calculated CH${}_{\mathrm{4}}$ solubility decreased from 0.016 mol L${}^{-\mathrm{1}}$ in July 2015 to 0.015 mol L${}^{-\mathrm{1}}$ at the end of November 2015 and increased to almost 0.017 mol L${}^{-\mathrm{1}}$ in May 2016 (Fig. c). This long-term trend was mainly caused by temperature variability ($R=-\mathrm{0.99}$), while tidal pressure changes caused a semi-diurnal variation of $\pm \sim$0.005 mol L${}^{-\mathrm{1}}$.

Dissolved O${}_{\mathrm{2}}$ decreased from $\sim$ 385 $\mathrm{\mu }$mol L${}^{-\mathrm{1}}$ in July 2015 to $\sim$ 350 $\mathrm{\mu }$mol L${}^{-\mathrm{1}}$ at the beginning of December 2015 and increased to $\sim$ 400 $\mathrm{\mu }$mol L${}^{-\mathrm{1}}$ towards 6 May 2016 (Fig. d) and followed temperature inversely ($R=-\mathrm{0.94}$), with similar long- and short-term variability.

The averaged bottom water current (81 m above the seafloor) was 4 cm s${}^{-\mathrm{1}}$ in a northwestward direction (321${}^{\circ }$ N) (Fig. c). The current usually had one counterclockwise rotation every 23.93 h period, corresponding to the diurnal K1 tidal constituent (tide with period 23.93 h; see ) with a secondary semi-diurnal (M2) modulation.

3.2

Time series at site O${}_{\mathrm{246}}$

CH${}_{\mathrm{4}}$ concentration at site O${}_{\mathrm{246}}$ ranged from 10 $\pm$ 3 nmol L${}^{-\mathrm{1}}$ on 21 September 2015 to 2727 $\pm$ 182 nmol L${}^{-\mathrm{1}}$ on 18 August 2015, with 2.5 and 97.5 percentiles of 107 and 1374 nmol L${}^{-\mathrm{1}}$. The data approximately follow lognormal distribution with average and median of 577 and 600 nmol L${}^{-\mathrm{1}}$, respectively, and an interquartile range of 293–721 nmol L${}^{-\mathrm{1}}$. The median ROC of CH${}_{\mathrm{4}}$ (31 nmol L${}^{-\mathrm{1}}$ min${}^{-\mathrm{1}}$) was almost 20 times higher than site O${}_{\mathrm{91}}$ (Fig. b and Appendix ). There was also clear diurnal periodicity in CH${}_{\mathrm{4}}$ concentration at O${}_{\mathrm{246}}$. The long-term trend (2-week running mean) shows decreasing concentrations until 3 October 2015 (end of the measuring period, Fig. b). Dissolved O${}_{\mathrm{2}}$ decreased from $\sim$ 380 to $\sim$ 300 $\mathrm{\mu }$mol L${}^{-\mathrm{1}}$ and was negatively correlated with water temperature ($R=-\mathrm{0.61}$; see Table  for the complete correlation matrix).

Table 2

Correlation coefficients between variables at O${}_{\mathrm{91}}$. The terms RTC CH${}_{\mathrm{4}}$ and raw CH${}_{\mathrm{4}}$ refer to response-time-corrected and untreated CH${}_{\mathrm{4}}$ (see Sect.  and Appendix ).

	RTC CH${}_{\mathrm{4}}$	Raw CH${}_{\mathrm{4}}$	Temperature	Salinity	Oxygen	Pressure	Solubility	Wind speed	CO${}_{\mathrm{2}}$
	mol L${}^{-\mathrm{1}}$	mol L${}^{-\mathrm{1}}$	${}^{\circ }$C		mol L${}^{-\mathrm{1}}$	dbar	mol L${}^{-\mathrm{1}}$	m s${}^{-\mathrm{1}}$	$\mathrm{\mu }$atm
RTC CH${}_{\mathrm{4}}$	1	0.78	$-$0.31	$-$0.24	0.30	0.15	0.33	$-$0.29	$-$0.13
Raw CH${}_{\mathrm{4}}$	0.78	1	$-$0.45	0.26	0.48	0.10	0.45	$-$0.44	$-$0.09
Temperature	$-$0.31	$-$0.45	1	0.87	$-$0.61	$-$0.02	$-$0.99	0.38	0.22
Salinity	$-$0.24	$-$0.26	0.87	1	$-$0.22	$-$0.03	$-$0.87	0.07	0.13
Oxygen	0.30	0.48	$-$0.61	$-$0.22	1	0.06	0.59	$-$0.65	$-$0.41
Pressure	0.15	0.01	$-$0.02	$-$0.03	0.06	1	0.16	$-$0.05	0.14
Solubility (CH${}_{\mathrm{4}}$)	0.33	0.45	$-$0.99	$-$0.87	0.59	0.16	1	0.38	$-$0.20
Wind speed	$-$0.29	$-$0.44	0.38	0.07	$-$0.65	$-$0.05	0.38	1	0.18
CO${}_{\mathrm{2}}$	$-$0.13	$-$0.09	0.22	0.13	$-$0.41	0.14	$-$0.20	0.41	1

Temperature and salinity increased from $\sim$ 2.5 to $\sim$ 4.0 ${}^{\circ }$C and $\sim$ 34.85 up to $\sim$ 35.0, respectively, from the deployment until October 2015 (Fig. b), with Atlantic water being dominant throughout the measuring period. Rapid shifts of around $\pm$ 1 ${}^{\circ }$C and 0.05 salinity occurred occasionally over a period of hours to days.

Variance in hydrostatic pressure was mainly explained by the tides (95.2 %), which were mainly governed by the semi-diurnal M2 tide, with weaker diurnal and fortnightly modulation (Fig. b). Changes in pressure varied from $\sim$ 1.2 to $\sim$ 1.5 m during periods of neap and spring tide.

Being governed mainly by temperature ($R=-\mathrm{0.99}$), CH${}_{\mathrm{4}}$ solubility dropped from 0.042 to 0.040 mol L${}^{-\mathrm{1}}$ from the deployment in July until October 2015, with a semi-diurnal variation of $\sim$ 0.005 mol L${}^{-\mathrm{1}}$ due to tidal changes in hydrostatic pressure.

The averaged current was $\sim$ 10 cm s${}^{-\mathrm{1}}$ northward (7${}^{\circ }$ N) (Fig. c). Variability in the along-slope current ($-$10${}^{\circ }$ N direction) was strongly related to the semi-diurnal M2 tidal component, while the cross-slope currents were governed by the diurnal K1 frequency. The bottom water current rotated counterclockwise with a period of 23.93 h (K1 tidal constituent) and semi-diurnal modulation in the along-slope component. Dissolved CH${}_{\mathrm{4}}$ concentration was weakly anti-correlated with wind speed ($R=-\mathrm{0.29}$), temperature ($R=-\mathrm{0.31}$), and salinity ($R=-\mathrm{0.24}$) and positively correlated with CH${}_{\mathrm{4}}$ solubility ($R=\mathrm{0.33}$) and oxygen ($R=\mathrm{0.3}$).

4 Discussion

4.1

CH${}_{\mathrm{4}}$ variability

Combining mapped flares and flow rates from the recovery cruise (May 2016) with bottom water current velocity (9 m above the seafloor) reveals that CH${}_{\mathrm{4}}$ concentration was strongly affected by whether water was advected from areas where we mapped strong or weak seepage in May 2016 (Fig. ). Strong seeps (flow rate $>$ 200 mL${}^{-\mathrm{1}}$ min${}^{-\mathrm{1}}$) were mainly located between $\sim$ 30 and 80 m to the north or northeast of site O${}_{\mathrm{91}}$, and only weak and more distant seepage was observed southwest of the observatory (Fig. a). Consequently, the averaged CH${}_{\mathrm{4}}$ concentration from water coming from northeast was $\sim$ 440 nmol L${}^{-\mathrm{1}}$, while water from the southwest averaged $\sim$ 100 nmol L${}^{-\mathrm{1}}$. Similarly, a strong CH${}_{\mathrm{4}}$ seep (flow rate $\sim$ 1200 mL min${}^{-\mathrm{1}}$) was mapped $\sim$ 40 m north of site O${}_{\mathrm{246}}$, making water advected from this direction highly elevated in CH${}_{\mathrm{4}}$, with an average of $\sim$ 1400 nmol L${}^{-\mathrm{1}}$ compared to the overall average of 577 nmol L${}^{-\mathrm{1}}$ (Fig. b). The rapid changes in dissolved CH${}_{\mathrm{4}}$ can to a high degree be explained by this relationship, due to the high variability in ocean current velocity. That this relationship holds for most of the measuring period also shows that even though observed average concentrations are lower in winter months, the seep configuration did not change significantly from July 2015 to May 2016, and dissolved CH${}_{\mathrm{4}}$ was efficiently dispersed in relatively high concentrations in the whole seepage area.

Furthermore, daily CH${}_{\mathrm{4}}$ concentrations at site O${}_{\mathrm{91}}$ were higher on average than the 24 h measurements (313 vs. 200 nmol L${}^{-\mathrm{1}}$). This can be explained by the comparable measurement periodicity (24 h) and tidal periodicity (23.93 h) in the ocean currents, resulting in predominantly eastward advection during daily measurements, thus systematically transferring water from a weak seepage area (Fig. ). We did not observe this effect at site O${}_{\mathrm{246}}$, most likely due to less tidal variance in the current direction (Fig. b). Nonetheless, this systematic tide-induced bias on the daily measurements at site O${}_{\mathrm{91}}$ highlights the importance of taking the oceanographic conditions into account to avoid misinterpretation of variability.

Figure 3

O${}_{\mathrm{91}}$ (a) and O${}_{\mathrm{246}}$ (b) location (yellow dot) and flow rates from flares mapped in their vicinity during CAGE 16-4 (color scale). Background color (green and blue) illustrates seafloor bathymetry. The compass diagram shows the relationship between ocean current direction (angle) and CH${}_{\mathrm{4}}$ concentration (distance from center; black is response-time-corrected (RTC) data, and raw data are in blue).

[Figure omitted. See PDF]

Since currents are mostly northward and seepage is mostly located to the north of both observatories, averaged measured CH${}_{\mathrm{4}}$ concentrations are likely lower than the average over the immediate surrounding area (Fig. ). Despite this, the observatory data show higher average CH${}_{\mathrm{4}}$ concentrations than previously reported. In the area surrounding site O${}_{\mathrm{91}}$, reported average concentrations of 92, 70 and 61 nmol L${}^{-\mathrm{1}}$ in June 2014, July 2015 and May 2016, respectively, based on discrete water sampling. Averaged CH${}_{\mathrm{4}}$ concentrations measured at site O${}_{\mathrm{91}}$ in July 2015 and May 2016 were 566 and 110 nmol L${}^{-\mathrm{1}}$, respectively, i.e., around 8 and 2 times higher than values reported by Silyakova et al. (2020), respectively. The maximum CH${}_{\mathrm{4}}$ concentration at O${}_{\mathrm{91}}$ of 1748 $\pm$ 142 nmol L${}^{-\mathrm{1}}$ on 20 August 2015 also significantly exceeds the previously maximum recorded concentration in the area of 480 nmol L${}^{-\mathrm{1}}$ (July 2014, Silyakova et al., 2020). At site O${}_{\mathrm{246}}$ the August 2016 average (564 nmol L${}^{-\mathrm{1}}$) was 8 times higher than what found in August 2010 (70 nmol L${}^{-\mathrm{1}}$) using an altimeter-controlled CTD towed at 2 m above the seafloor. The maximum concentration in August 2016 also significantly exceeded previous observations, with 2661 $\pm$ 163 nmol L${}^{-\mathrm{1}}$ compared to 524 nmol L${}^{-\mathrm{1}}$ measured by Gentz et al. (2014).

These differences could be a result of temporal, local or regional differences in CH${}_{\mathrm{4}}$ concentration. However, strong vertical gradients in dissolved CH${}_{\mathrm{4}}$ are well documented at both seep sites , and our sensors measured closer to the seafloor (1.2 m above seafloor) compared to Gentz et al. (2014) (2 m above seafloor) and Silyakova et al. (2020) (5 to 15 m above seafloor). Additionally, the observatories were deployed close to seeps using a launcher, as opposed to “blind” water sampling from a shipborne rosette. Methane was also measured in situ, thereby avoiding potential CH${}_{\mathrm{4}}$ outgassing after retrieval of water samples .

Dissolved CH${}_{\mathrm{4}}$ within shallow seep sites where gas can bypass the oceanic sinks often present heterogeneous distribution and rapid temporal variability . Our results show that the temporal variability at the two seep sites is higher than previously reported and that changing ocean currents and configuration of nearby seeps are major contributors. This high short-term variability introduces a conceptual error in studies relying on discrete water sampling (e.g., to calculate inventories) because the time required to conduct the survey (measured in days) is much longer than large temporal variations in concentration (up to an order of 10${}^{\mathrm{3}}$ nmol L${}^{-\mathrm{1}}$ within hours).

Figure 4

Relative standard error of the mean for different numbers of samples $N$ for O${}_{\mathrm{91}}$ 24 h data, data presented in Silyakova et al. (2020) (“June-14”, “July-15”, “May-16”), and O${}_{\mathrm{246}}$ 24 h data (in black). Relative standard deviation (corresponding to the standard error with $N=\mathrm{1}$) is given in the legend (${\mathit{\sigma }}_{\mathrm{rel}}$). * Data from Silyakova et al. (2020) calculated assuming that the sample distribution resembles the underlying distribution (see Appendix ).

[Figure omitted. See PDF]

We can obtain a first-order constraint on errors caused by short-term variability in a hypothetical water sampling survey using the 24 h time series from the observatories. We assume the hypothetical survey seeks to find the average concentration in the bottom layer of the seep site. The expected error can then be found by calculating the standard error of the mean (SEM) for a given number of samples $N$ using the 24 h time series as an underlying distribution representing the sub-daily variability of the seep site (Fig. ; Appendix  contains a detailed outline of the methodology). Even though surveys often require more than 24 h to complete (2–3 d in Silyakova et al., 2020), a majority of processes causing short-term variability have periods below or at $\sim$ 24 h (for instance tides and many turbulent eddies; see, e.g., Sect. 3.2 and 3.1 and ), likely making the daily distribution relevant also for surveys with longer duration. We compared SEM calculations based on the observatory 24 h time series with SEM calculations for the bottom water ($\sim$ 5 m above the seafloor) discrete water sample data used for average and inventory estimates of the O${}_{\mathrm{91}}$ seep site in Silyakova et al. (2020) (also included in Fig. ).

The absolute SEM (in nmol L${}^{-\mathrm{1}}$) is generally higher for time series with higher averaged concentrations, making the relative SEM cluster well, with gradually diminishing range for increasing $N$ (an inherent property of the SEM, e.g., 12 %–45 % for $N=\mathrm{10}$, 9 %–30 % for $N=\mathrm{30}$, Fig. ). The SEM of the data from Silyakova et al. (2020) is similar to the SEM of the 24 h time series, with a common range of 5 %–15 % expected error for surveys with $N\sim \mathrm{60}$ samples ($N=\mathrm{64}$, 62, and 63 in Silyakova et al. 2020). It should be noted that the comparison with data from Silyakova et al. (2020) has caveats, e.g., that the observatory data do not contain errors due to spatial variability and an assumption of representative short-term temporal variability at the observatory sites (see also Appendix ).

Evidently, detailed surveys of individual seep sites, such as the study by Silyakova et al. (2020), can provide reasonable estimates of local inventories ($<$ 15 % uncertainty) despite high short-term temporal variability. However, it is important to note that the area investigated in Silyakova et al. (2020) was densely mapped and homogeneous in the sense that it is an area where seepage is well documented . Interpolation or averaging across larger regions where the amount of seepage is mostly unknown can result in considerable errors due to false interpolation assumptions and amplification of individual measurement errors, which can be large (expected errors up to $\sim$ 140 % for single measurements; see listed standard deviations in Fig. ). These effects can potentially explain some of the discrepancies in estimates of oceanic CH${}_{\mathrm{4}}$ inventories and fluxes.

Our findings stress the importance of sufficiently dense mapping and knowledge about the underlying seep condition when collecting water samples for inventory estimates. They also highlight the advantage of towed or autonomous instrumentation capable of providing continuous CH${}_{\mathrm{4}}$ data, giving considerably better coverage and representation of the CH${}_{\mathrm{4}}$ distribution in less time e.g.,. Assuming a distribution that better reflects the uneven spread of CH${}_{\mathrm{4}}$ when applying interpolation and extrapolation techniques could also limit estimation errors. Future studies should investigate how initial errors due to short-term and small-scale variability propagate via different upscaling techniques and how these errors can be mitigated.

4.2 Hydrostatic pressure

Tidal changes in hydrostatic pressure can trigger CH${}_{\mathrm{4}}$ release by build-up of CH${}_{\mathrm{4}}$ in sediment pore water at rising tide and subsequent release when pore pressure decreases at falling tide, as observed at the Hikurangi Margin and Clayoquot Slope . Our study sites differ from these sites in depth (they are $>$600 m) and in tidal amplitude (4 m at Clayoquot Slope compared to 1.5 offshore Prins Karls Forland). and also observed bubbles hydro-acoustically, while we measure dissolved CH${}_{\mathrm{4}}$, which is strongly affected by the (also tidally dependent) current direction (Fig. ).

To evaluate the effect that hydrostatic pressure changes have on the in situ concentration, we need to constrain the variance caused by changing current directions (since they operate in the same frequency domain). To do this, we first binned the CH${}_{\mathrm{4}}$ concentration data into overlapping bins defined by the current direction at the time when the measurement was obtained and calculated standard scores (the number of standard deviations each value deviates from the sample mean; see, e.g., ) for the data in each bin. We used larger current direction intervals for O${}_{\mathrm{246}}$ due to the shorter data set, with a 12${}^{\circ }$ window for O${}_{\mathrm{91}}$ and a 30${}^{\circ }$ window for O${}_{\mathrm{246}}$. This resulted in a data set (i.e., the standard scores from all bins) effectively unrelated to the current direction. We then binned all the standard-scored CH${}_{\mathrm{4}}$ data according to when the data were collected in relation to the M2 governed tidal cycle peak using overlapping 30 min bins (the M2 tide explains 79.2 % and 80.3 % of the pressure variance at O${}_{\mathrm{91}}$ and O${}_{\mathrm{246}}$, respectively). Average and median values were calculated for each bin, giving the averaged or median normalized dissolved CH${}_{\mathrm{4}}$ value (standard score) for each current velocity defined data bin as a function of the M2 tidal cycle (Fig. ). This partial decoupling of variability in hydrostatic pressure and current direction was possible since the bottom water current and hydrostatic pressure changes had different dominant tidal constituents; i.e., the current was mainly dominated by the diurnal K1 constituent ($\sim$ 23.91 h period), while the M2 tide is semi-diurnal (12.42 h period).

A strong effect of the hydrostatic pressure on local seepage should elevate the standard scores at decreasing pressure (from 0 to 6.2 h, i.e., in the right half of Fig. ), which we observe at both observatories. However, we observe stronger peaks at increasing hydrostatic pressure ($-$3 h) at site O${}_{\mathrm{91}}$ and at the M2 peak (0 h) at site O${}_{\mathrm{246}}$, which contradicts this hypothesis. This does not mean that there is no effect of hydrostatic pressure changes but rather that the seepage in the area is widespread at both falling and rising tide conditions. The high variability caused by the strong effect of current direction also makes it particularly challenging to detect moderate changes in seepage intensity.

Figure 5

Median and averaged standard scores of CH${}_{\mathrm{4}}$ binned according to bottom water current direction and according to where the data were sampled in relation to the phase of the M2 pressure tide.

[Figure omitted. See PDF]

4.3 Bottom water temperature

Bottom water temperature can affect CH${}_{\mathrm{4}}$ release by altering hydrate stability and CH${}_{\mathrm{4}}$ solubility in pore water and the water column . Seasonal CH${}_{\mathrm{4}}$ release variability resulting from temperature variations in the bottom water has been linked to migration of the gas hydrate stability zone (GHSZ) and hydrate dissociation further offshore at $\sim$ 390 m water depth . Our observatories were deployed in areas too shallow for gas hydrate to form. However, inversely varying seepage intensity between seepage at the GHSZ depth (390 m) and site O${}_{\mathrm{246}}$ can suggest that these areas are fed by the same hydrocarbon source and that hydrates seasonally block the lateral pathways between these seep sites . This is in agreement with the observed long-term ($\sim$ 3 months) negative correlation between bottom water temperature and dissolved CH${}_{\mathrm{4}}$ at site O${}_{\mathrm{246}}$ ($R=-\mathrm{0.31}$). It should be noted that the same relationship is observed at O${}_{\mathrm{91}}$; however, no geophysical data are available from this area due to the shallow depth.

Tidal pressure variations can affect CH${}_{\mathrm{4}}$ release via pore water solubility (Sect. ), but on longer timescales CH${}_{\mathrm{4}}$ solubility is almost exclusively a function of water temperature. Higher CH${}_{\mathrm{4}}$ solubility implies more CH${}_{\mathrm{4}}$ dissolved in pore water and within bubble streams, potentially increasing the amount of CH${}_{\mathrm{4}}$ dissolved in bottom water. A small but significant ($R=\mathrm{0.33}$) positive correlation between CH${}_{\mathrm{4}}$ solubility and concentration at site O${}_{\mathrm{246}}$ and site O${}_{\mathrm{91}}$ (considering the same time period, i.e., until 3 October in 2015) could indicate such an effect. This is also an alternative explanation for the negative correlation between temperature and CH${}_{\mathrm{4}}$ concentration at site O${}_{\mathrm{246}}$.

4.4 Pore water seepage

Short-term temperature increase further offshore (390 m depth) has been linked with release of warm, CH${}_{\mathrm{4}}$-rich fluids from the sediments triggered by short-duration seismic events . This means that increased CH${}_{\mathrm{4}}$ concentration should be accompanied by increased water temperature and reduced salinity due to admixture of warmer, less saline pore water. We compared short-term anomalies (i.e., deviations from daily means) in these three variables in the 24 h data sets at both seep sites but found no corroborating evidence for this hypothesis. Instead, the covariance between current velocity and temperature and salinity anomalies indicates that short-term variability is mainly caused by cross-shelf exchange of Atlantic water in the West Spitsbergen Current and the colder, fresher Arctic water in the Coastal Current due to eddies . It also indicates that CH${}_{\mathrm{4}}$ release comes mainly from bubble dissolution and not from pore water seepage.

Figure 6

(a) Bulk velocity shear ($\mathrm{\Delta }H=\mathrm{8}$ m) and two-dimensional correlation with wind stress (contours). Relationships between 5 d accumulated wind stress and median (b) and median absolute deviation (c) of CH${}_{\mathrm{4}}$ concentration for 24 h data periods. Persistent wind events with speeds of more than 10 m s${}^{-\mathrm{1}}$ in periods over 6 h are indicated with blue stars along the $x$ axis of panel (a). Blue highlights fall and winter water column conditions as described in the text.

[Figure omitted. See PDF]

4.5

Seasonal variation of CH${}_{\mathrm{4}}$ distribution at site O${}_{\mathrm{91}}$

Low release of CH${}_{\mathrm{4}}$ to the atmosphere from the O${}_{\mathrm{91}}$ seep area during summer despite high seabed influx has been explained by suppression of vertical mixing by strong stratification or absence of mechanical forcing such as wind stress . However, in fall and winter, the water column offshore Prins Karls Forland is expected to have more horizontal and vertical mixing due to weaker stratification from cooling or sea ice formation , baroclinic instability in the frontal structures of the West Spitsbergen Current , and more frequent storms .

We expect lower CH${}_{\mathrm{4}}$ variability and lower CH${}_{\mathrm{4}}$ concentration during periods of high mixing and dispersion due to weaker horizontal and vertical gradients and more efficient dispersion of CH${}_{\mathrm{4}}$ away from sources. We use three sets of parameters to evaluate long-term changes in the amount of mixing in the water column (see Appendix ): (i) the 4-week averaged bulk velocity shear ($S_{\mathrm{b}}$), (ii) the two-dimensional correlation between wind stress and current velocity ($R_{\mathrm{WC}}$), and (iii) the number of stormy days as defined by persistent winds $>$ 11 ms${}^{-\mathrm{1}}$ lasting longer than 6 h (Fig. ). Calm weather and low $S_{\mathrm{b}}$ and $R_{\mathrm{WC}}$ until mid-September 2015 indicate a stable water column with limited mixing in the bottom waters. From mid-September, $S_{\mathrm{b}}$ increased and stayed high until mid-November, together with a gradual increase in $R_{\mathrm{WC}}$ that can be attributed to a gradual breakdown of stratification and increasing number of storm events (Fig. a). $R_{\mathrm{WC}}$ remained high ($R_{\mathrm{WC}}>\mathrm{0.5}$ at 60 m depth) until March 2016, indicating a significant effect of wind forcing in the water column. From March until observatory retrieval, $R_{\mathrm{WC}}$ decreased to $<\mathrm{0.2}$ below 50 m depth while $S_{\mathrm{b}}$ increased below 60 m depth, indicating available energy for mixing in the bottom waters.

We quantified CH${}_{\mathrm{4}}$ variability during the 24 h measurements using the median absolute deviation (MAD) and used the median as a measure of the amount of dissolved CH${}_{\mathrm{4}}$. The three 24 h periods collected during the calmer period prior to mid-September had high median concentration ($>$ 300 nmol L${}^{-\mathrm{1}}$) and the overall highest variability (MAD $>$ 160 nmol L${}^{-\mathrm{1}}$), as expected for low mixing conditions (Fig. b and c). From mid-September until the end of March (i.e., fall and winter seasons), the 24 h CH${}_{\mathrm{4}}$ concentration time series had generally lower MAD and median concentration. In this period, CH${}_{\mathrm{4}}$ variability and median also showed a good statistical relationship with the 5 d accumulated wind stress ($R=-\mathrm{0.82}$ for MAD and $R=-\mathrm{0.61}$ for median concentration), indicating that wind forcing has a deep impact on mixing and redistribution of CH${}_{\mathrm{4}}$ in the water column (which also fits well with a high $R_{\mathrm{WC}}$). The two last 24 h CH${}_{\mathrm{4}}$ time series (10 April and 1 May) had low median concentration, which could be explained by the absence of stratification and generation of mixing from the observed increase in $S_{\mathrm{b}}$.

Accumulated wind stress, $S_{\mathrm{b}}$ and $R_{\mathrm{WC}}$ are only limited indicators of water column dispersion and mixing. Nonetheless, the relationship between these parameters and the MAD and medians of the 24 h period CH${}_{\mathrm{4}}$ time series gives a good indication of the seasonal cycle of distribution and vertical transport of CH${}_{\mathrm{4}}$: strong stratification, less wind forcing and eddy activity in summer limit mixing and prevent CH${}_{\mathrm{4}}$ from reaching the atmosphere. However, in fall and winter, reduced stratification makes the water column more prone to mixing, and distribution of CH${}_{\mathrm{4}}$ seems to be strongly linked with wind forcing from September to April.

5 Conclusions

Time series of dissolved CH${}_{\mathrm{4}}$ at both lander locations show considerably higher CH${}_{\mathrm{4}}$ concentrations (up to 1748 $\pm$ 142 nmol L${}^{-\mathrm{1}}$ at O${}_{\mathrm{91}}$ and 2727 $\pm$ 182 nmol L${}^{-\mathrm{1}}$ at O${}_{\mathrm{246}}$) than previously found in ship-based water sampling surveys (maximum of 482 near O${}_{\mathrm{91}}$ and of 564 near O${}_{\mathrm{246}}$). The time series also uncover high CH${}_{\mathrm{4}}$ variability (up to $\sim$ 1000 nmol L${}^{-\mathrm{1}}$) within short timescales ($<\mathrm{24}$ h), highlighting the potential uncertainty of flux and inventory estimates based on interpolation and extrapolation techniques relying on, e.g., an assumption of linearity. We calculated the standard error of a mean estimate based on a hypothetical discrete water sampling survey based on a range of samples by using the 24 h time series as the underlying distribution. The results aligned well with previous discrete water sampling surveys in the area, giving a standard error of the mean of 5 %–15 % for $\sim$ 60 samples.

Variability can be linked to directional ocean current variations occurring at tidal timescales, which shows the importance of taking the current direction and seep locations into account when interpreting intense seep site observations. The persistent relationship between current direction and location of seeps during recovery shows that there was seepage throughout the year and that the seep configuration was relatively constant.

We did not observe a direct effect of tidal pressure variations on CH${}_{\mathrm{4}}$ release, but this could be hidden by the strong effect of variations in current direction. A negative (long-term) correlation between temperature and dissolved CH${}_{\mathrm{4}}$ at O${}_{\mathrm{246}}$ is in agreement with the hypothesized seasonal blocking of lateral CH${}_{\mathrm{4}}$ pathways in the sediments but could also be explained by increased CH${}_{\mathrm{4}}$ solubility in the water column.

Short-term, small-scale variations in temperature and salinity were not linked with increased amounts of dissolved CH${}_{\mathrm{4}}$, but rather with cross-frontal exchange of water masses due to eddies.

We observed a seasonal cycle in the characteristics of the 24 h time series that fits with seasonal changes in dispersion and mixing characteristics of the water column. Higher CH${}_{\mathrm{4}}$ concentration and variability in early fall, when stratification was strong, was followed by lower median concentrations and variability in late fall and winter when the water column was more affected by mixing. In late fall and winter, wind forcing was statistically coupled to the concentration and variability of CH${}_{\mathrm{4}}$, probably due to weaker water column stratification.

When estimating the atmospheric impact of a particular CH${}_{\mathrm{4}}$ source based on sparse measurements, it is crucial to have some constraints on the temporal and spatial variability. These constraints can either be direct knowledge about variability itself or how inventory and fluxes are affected by related physical and/or chemical parameters. We observed considerable temporal and spatial variability at the two seep sites that need to be taken into account to obtain meaningful estimates of CH${}_{\mathrm{4}}$ fluxes or inventories. That no strong direct link was found with other oceanographic parameters illustrates the nonlinearity of the system, making careful interpretation of measurements important. Future studies should aim to identify the errors that arise via different upscaling and interpolation techniques and how these errors can be mitigated. Based on our observations, we suggest that uncertainties in CH${}_{\mathrm{4}}$ inventory and seep estimates can be mitigated by taking the local seep configuration, ocean currents and mixing rates into account and employing autonomous instrumentation capable of resolving the steep horizontal gradients in dissolved CH${}_{\mathrm{4}}$. This, alongside direct measurements of seepage by, e.g., acoustic instrumentation, can help constrain future estimates of CH${}_{\mathrm{4}}$ flux to the atmosphere from seabed seepage.

Appendix A The K-Lander

Figure A1

(a) The K-Lander is a 1.6 m high and 3.6 m wide trawl-proof stainless steel frame with multiple instrument mounts and batteries. The side panels are perforated to allow unobstructed water flow to the instruments inside the structure. See Appendix  for details on instrumentation. (b) One of the K-Landers during deployment with a launcher mounted on top and camera system mounted on a boom for visual control of landing area. (c) The two K-landers before deployment.

[Figure omitted. See PDF]

Appendix B Measurement intervals, general post-processing and data

The CTD (oxygen sensor) and ADCP conducted measurements every 4 and 9 min, respectively, during the continuous monitoring of CH${}_{\mathrm{4}}$ and CO${}_{\mathrm{2}}$ measurements and every 21 and 29 min, respectively, during the rest of the deployment period (see Table  for acronyms, description and measurement accuracy). Salinity was measured on the practical salinity scale.

The upward-mounted ADCP measured ocean currents in 1 m bins with a bottom 7 m blank distance, where the topmost 20 % of the water column was disregarded due to side lobe interference. The high resolution, relatively short ensemble time (1 min) and potential presence of CH${}_{\mathrm{4}}$ bubbles in the water resulted in noisy data. We dampened the noise by first removing any data points with error velocities exceeding one short-term (1 week) standard deviation, smoothed the data using a second-order Butterworth low-pass filter with a 3 h cutoff period and a spatial (i.e., vertical) moving average filter with a 5 m Hann window (increasing the blank distance to 10 m). The accuracy of the ADCP data is therefore not explicitly constrained and is based on comparing current velocity frequency spectra before and after filtering, combined with averaged error velocity of the raw data (Table ).

Since sensors were recording at different frequencies, chronological alignment of the data was carried out by identifying nearest neighbor data points or by resampling. For correlation coefficients, histograms and Fourier analysis, the data sets were resampled to a uniform 15 min or 1 h measuring interval depending on the sample frequency of the raw data, using a poly-phase anti-aliasing filter. Due to the power-cycling mode of the CH${}_{\mathrm{4}}$ and CO${}_{\mathrm{2}}$ sensors and differing sampling frequencies, some statistics were based on more data points than others (outlined in Table ). Daily measurements of CH${}_{\mathrm{4}}$ were excluded from these statistics due to the high probability of systematic errors induced by periodic diurnal effects.

Harmonic analysis of hydrostatic pressure and ocean currents was done using t_tide see and the fast Fourier transform.

We calculated the rate of change (ROC) in CH${}_{\mathrm{4}}$ concentration using the response-time-corrected CH${}_{\mathrm{4}}$ data and the absolute value of the three-point (9 min) finite differences to limit the effect of noise on the calculation.

The absolute concentration of CH${}_{\mathrm{4}}$ in the water (nmol L${}^{-\mathrm{1}}$) was estimated from the partial pressure of CH${}_{\mathrm{4}}$, pressure, temperature and salinity using Henry's law and Henry constants obtained from and the practical molar volume and gamma term from .

The CH${}_{\mathrm{4}}$ sensors were calibrated to relevant water temperatures prior to deployment. The TDLAS detectors provide measurements with good selectivity (fit for purpose) and high long-term stability (intermittent calibration not necessary) and are unaffected by dissolved oxygen content (unless via complete depletion). Biofouling was also minimal at retrieval (due to the cold water and local setting), and the PDMS membranes are almost unaffected by cold water. Generally, we did no observations indicating issues with any of the sensors except for what has already been mentioned regarding the conductivity probe and electrical malfunction of O${}_{\mathrm{246}}$. Furthermore, we discarded all data recorded during instrument warm-up (i.e., when internal temperature was below correct operating temperature) before the individual measurement periods (the instruments were turned on $\sim$ 35 min prior to recording the data used in the analysis).

In the Contros HydroC CH${}_{\mathrm{4}}$ and CO${}_{\mathrm{2}}$ sensors, dissolved gases diffuse through a hydrophobic membrane into a gas chamber and equilibrate with the ambient environment. This results in the slow response time (e.g., ${\mathit{\tau }}_{\mathrm{63}}\sim$50 min under certain conditions for our membrane and pump setup for the CH${}_{\mathrm{4}}$ sensor) and poor representation of the rapid changes in CH${}_{\mathrm{4}}$ we expected in our study area . We therefore performed a response time correction for the dissolved CH${}_{\mathrm{4}}$ data following the methodology presented in , modulating the response time using the temperature data (effects of salinity on membrane permeability were not taken into account since these are negligible for the local ranges; see ). The CO${}_{\mathrm{2}}$ sensors had a faulty pump, which ambiguously increased the response time of the sensors, making response time correction impossible.

The response time correction was performed for each period individually (1 and 24 h, i.e., 377 periods), using the stated measurement accuracy of the instrument (2 $\mathrm{\mu }$atm or 3 % of measured value, whichever is higher) as input uncertainty. We first identified the ideal $\mathrm{\Delta }t$ according to the maximum curvature point in the L curves of the 24 h measurement periods. These varied slightly between each measurement period but averaged close to 180 s (176.4 s). To keep the same measuring interval for all the CH${}_{\mathrm{4}}$ data, we therefore corrected all the data with a specified $\mathrm{\Delta }t$ of 180 s, which falls well within the bend of the L curve and should therefore safeguard a good balance between noise and model error (Fig. a). Inspection of model fit residuals showed a slight modulation following the variance in the signal, which is explained by our choice to use the same 3 min measurement grid across a relatively wide variance range, but the residuals were otherwise Gaussian. Although this is expected, it indicates that errors might be slightly overestimated for low-variance sections of the time series and vice versa for high-variance sections.

The uncertainty estimate varies depending on the amount of CH${}_{\mathrm{4}}$ measured by the TDLAS unit in the measurement chamber of the instrument. The distribution of the uncertainty estimates is shown as percentages in Fig. b. Estimated uncertainty ranged from 3 to 205 nmol L${}^{-\mathrm{1}}$ (95 % confidence, high for high concentrations in the measurement chamber and vice versa) or usually between 5 % and 20 %, although there were some outliers when the concentration was low and uncertainty estimate was high (Fig. b).

Table B1

Instruments mounted on O${}_{\mathrm{91}}$ and O${}_{\mathrm{246}}$ (see Fig. ), measured parameters, height in meters above seafloor (masf), and stated accuracy. ADCP stands for acoustic Doppler current profiler. $N$ shows the number of data points used for later multi-variable analysis of O${}_{\mathrm{91}}$ and O${}_{\mathrm{246}}$.

Instrument	Parameter(s)	masf	$N$	Accuracy
Teledyne RDI ADCP WHLM75-3000	Current velocity profile	1.6	$\mathrm{17}\mathrm{438}/\mathrm{4731}$	$\sim$ 3 cm s${}^{-\mathrm{1}}$
Contros HydroC CH${}_{\mathrm{4}}$	$p$CH${}_{\mathrm{4}}$ (instrument output)${}^{*}$	1.2	$\mathrm{1491}/\mathrm{281}$	$\sim$ 5 %–20 % (RTC)${}^{***}$
	xCH${}_{\mathrm{4}}$ (reported)${}^{**}$
Contros HydroC CO${}_{\mathrm{2}}$	$p$CO${}_{\mathrm{2}}$	1.2	$\mathrm{1491}/\mathrm{281}$	n/a (no pump)
Seabird SBE16 plus V2	Conductivity–Temperature–Depth	1.2	$\mathrm{29}\mathrm{660}/\mathrm{9065}$	0.0005 Sm${}^{-\mathrm{1}}/\mathrm{0.005}$ ${}^{\circ }$C,
				/0.02 % of range
Seabird SBE63 oxygen optode	Dissolved oxygen	1.2	$\mathrm{29}\mathrm{660}/\mathrm{9065}$	3 $\mathrm{\mu }$mol kg${}^{-\mathrm{1}}$ or $\pm$ 2 %

${}^{*}$ The Contros HydroC CH${}_{\mathrm{4}}$ outputs partial pressure from the internal gas chamber. ${}^{**}$ We report absolute concentration in seawater (nmol L${}^{-\mathrm{1}}$) using Henry's law. ${}^{***}$ We report accuracy only for response-time-corrected (RTC) concentrations (see Fig. ) since the accuracy for untreated CH${}_{\mathrm{4}}$ concentration data is ambiguous due to the slow response time. n/a: not applicable.

Figure B1

(a) L curve for response time correction of CH${}_{\mathrm{4}}$ data showing the location of the chosen $\mathrm{\Delta }t$ (180 s) for 6 May at O$\mathrm{91}$. (b) Estimated relative (%) uncertainty for response-time-corrected CH${}_{\mathrm{4}}$ data (both observatories).

[Figure omitted. See PDF]

Appendix C

The 24 h measurements of CH${}_{\mathrm{4}}$

Figure C1

All 24 h periods of CH${}_{\mathrm{4}}$ concentration at O${}_{\mathrm{91}}$ using response-time-corrected data (black) with the uncertainty estimate (grey shading, 95 % confidence) and raw data (blue) from O${}_{\mathrm{91}}$.

[Figure omitted. See PDF]

Figure C2

All 24 h periods of CH${}_{\mathrm{4}}$ concentration at O${}_{\mathrm{246}}$ using response-time-corrected data (black) with the uncertainty estimate (grey shading, 95 % confidence) and raw data (blue) from O${}_{\mathrm{246}}$.

[Figure omitted. See PDF]

Appendix D Standard error of mean estimate due to temporal variability

To obtain the (theoretical) true dissolved CH${}_{\mathrm{4}}$ average or inventory for an area requires a known concentration everywhere at a single point $T_{\mathrm{0}}$ in time. Considering a hypothetical ship-based discrete water sampling survey, any small-scale spatial variability not resolved by the sampling grid or localized (not seep-wide) short-term temporal variability occurring during the survey time can be considered a measurement errors for the purpose of the survey. Assuming that the water samples are sufficiently spaced out to be considered independent samples, the estimated average concentration from $N$ samples in a particular depth layer at a seep site can be expressed as follows:

D1$${\displaystyle }E\left(m,{\mathit{ϵ}}_{\mathrm{t}},{\mathit{ϵ}}_{\mathrm{s}}\right)={\displaystyle \frac{{\sum }_{n=\mathrm{1}}^N\left(m+{{\mathit{ϵ}}_{\mathrm{t}}}_n+{{\mathit{ϵ}}_{\mathrm{s}}}_n\right)}{N}},$$ where $m$ is the average of the seep site at $T_{\mathrm{0}}$, ${\mathit{ϵ}}_{\mathrm{t}}$ is error due to temporal short-term deviation from $m$ at sampling time $T_{\mathrm{0}}+\mathrm{\Delta }t$ and ${\mathit{ϵ}}_{\mathrm{s}}$ is spatial deviation in concentration from $m$. The expected standard error of $E(m,{\mathit{ϵ}}_{\mathrm{t}},{\mathit{ϵ}}_{\mathrm{s}})$ from the short-term temporal and spatial variability is then given by D2$${\mathit{\sigma }}_{E\left(m,{\mathit{ϵ}}_{\mathrm{t}},{\mathit{ϵ}}_{\mathrm{s}}\right)}={\displaystyle \frac{\mathit{\sigma }}{\sqrt{N}}},$$ where $\mathit{\sigma }$ is the standard deviation of the distribution we sample from . From Eqs. () and () we obtain D3$${\mathit{\sigma }}_{E\left(m,{\mathit{ϵ}}_{\mathrm{t}},{\mathit{ϵ}}_{\mathrm{s}}\right)}={\mathit{\sigma }}_{E\left(m,{\mathit{ϵ}}_{\mathrm{t}}\right)}+{\mathit{\sigma }}_{E\left(m,{\mathit{ϵ}}_{\mathrm{s}}\right)}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{t}}}{\sqrt{N}}}+{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{s}}}{\sqrt{N}}},$$ where ${\mathit{\sigma }}_{\mathrm{t}}$ and ${\mathit{\sigma }}_{\mathrm{s}}$ are the distribution standard deviation related to temporal (${\mathit{ϵ}}_{\mathrm{t}}$) and spatial (${\mathit{ϵ}}_{\mathrm{s}}$) variability and ${\mathit{\sigma }}_{E(m,{\mathit{ϵ}}_{\mathrm{t}})}$ and ${\mathit{\sigma }}_{E(m,{\mathit{ϵ}}_{\mathrm{s}})}$ are the corresponding contributions to the standard error of the mean. Assuming the daily variance at the observatory is representative for the seep site, we can describe the expected error caused by sub-daily variability (all ${\mathit{ϵ}}_{\mathrm{t}}$) in a scenario where a seep site is being sampled $N$ times using the 24 h time series as the underlying distribution. In essence, we treat every measurement as having an associated probability distribution that is represented by the 24 h time series (which gives the sub-daily variability).

In the discrete water sample data presented in Silyakova et al. (2020), the underlying distribution is unknown, and we can only assume that the sample distribution resembles the underlying distribution, i.e., that D4$${\mathit{\sigma }}_{E\left(m,{\mathit{ϵ}}_{\mathrm{t}},{\mathit{ϵ}}_{\mathrm{s}}\right)}\approx {\widehat{\mathit{\sigma }}}_{E(m,{\mathit{ϵ}}_{\mathrm{t}},{\mathit{ϵ}}_{\mathrm{s}})}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{sampled}}}{\sqrt{N}}},$$ where ${\widehat{\mathit{\sigma }}}_{E(m,{\mathit{ϵ}}_{\mathrm{t}},{\mathit{ϵ}}_{\mathrm{s}})}$ is the standard error estimate of the mean based on the sample distribution and ${\mathit{\sigma }}_{\mathrm{sampled}}$ is the standard deviation of the measurements. All three data sets, i.e., “June-14” ($N=\mathrm{64}$), “July-15” ($N=\mathrm{62}$) and “May-16” ($N=\mathrm{63}$), have a similarly skewed distribution compared to what is found in the observatory data (see Fig. ), which supports this assumption. The survey in Silyakova et al. (2020) required 2–3 d to complete, while the observatory data only concern sub-daily variability (24 h time series). Nonetheless, we believe the comparison is valid, since the known major contributors to short-term (timescales shorter than weeks) variability act on sub-daily (or at least $\le$ daily) scales, such as the dominant frequencies in the ocean currents and pressure changes.

There is a clear relationship of increasing ${\mathit{\sigma }}_{E(m,{\mathit{ϵ}}_{\mathrm{t}},{\mathit{ϵ}}_{\mathrm{s}})}$ with increasing daily average, making relative ${\mathit{\sigma }}_{E(m,{\mathit{ϵ}}_{\mathrm{t}},{\mathit{ϵ}}_{\mathrm{s}})}$ a meaningful quantity to use, as opposed to absolute ${\mathit{\sigma }}_{E(m,{\mathit{ϵ}}_{\mathrm{t}},{\mathit{ϵ}}_{\mathrm{s}})}$. Additionally, for simplicity, we have not differentiated in the notation of the standard error of the mean (SEM) in the main text of the paper, referring to it as simply SEM in all situations.

It is also enlightening to consider the distribution of average estimates and how the skewed underlying distribution affects the distribution of average estimate errors for smaller $N$. We did this by simulating hypothetical surveys by random sampling from the 24 h data sets (Fig. ), which shows the elevated probability of underestimating the average for estimates based on few samples ($N\mathit{\lesssim }\mathrm{30}$), i.e., that the median error is smaller than the average error. This is caused by an inheritance of the skewed underlying distribution in the CH${}_{\mathrm{4}}$ concentration data (see Fig. a). This also allows for severe overestimates due to the long right-hand-side tail of the distribution. For larger values of $N$ ($N\mathit{\gtrsim }\mathrm{30}$), average estimates tend towards being normally distributed, thus avoiding these effects (see Fig. b).

Error estimates of more complicated properties, such as the total CH${}_{\mathrm{4}}$ content in a volume of water based on interpolation techniques, require an assessment of the individual uncertainties of each measurement and how these errors propagate via, e.g., linear interpolation in the spatial domain. While not being explicitly applicable to inventory estimates, the ${\mathit{\sigma }}_E$ still describes how random errors cancel out for larger values of $N$ in evenly sampled grids, assuming this variability is representative for the seep site.

Figure D1

Distribution of CH${}_{\mathrm{4}}$ concentrations from the (a) June 2014, (b) July 2015 and (c) May 2016 data in Silyakova et al. (2020) and (d) from the 24 h data (all periods) at O${}_{\mathrm{91}}$. Note the different scale for the $y$ axis between (a)–(c) and (d).

[Figure omitted. See PDF]

Figure D2

Histograms of simulated average estimates based on $N=\mathrm{10}$ (a) and $N=\mathrm{30}$ (b) samples from the 24 h data set from 23 August at O${}_{\mathrm{91}}$ showing the median and mean as vertical lines.

[Figure omitted. See PDF]

Appendix E Bulk velocity shear and wind stress correlation

We calculated bulk wind stress using 10 m above sea level ERA-interim reanalysis wind data . Water column bulk velocity shear $S_{\mathrm{b}}$ (see, e.g., ) was calculated as follows: E1$$S_{\mathrm{b}}^{\mathrm{2}}={\left({\displaystyle \frac{u_{\mathrm{u}}-u_{\mathrm{l}}}{h_{\mathrm{diff}}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{v_{\mathrm{u}}-v_{\mathrm{l}}}{h_{\mathrm{diff}}}}\right)}^{\mathrm{2}},$$ where $u_{\mathrm{u}}$, $u_{\mathrm{l}}$, $v_{\mathrm{u}}$ and $v_{\mathrm{l}}$ refer to the eastward and northward ADCP velocity components in the upper (subscript u) and lower (subscript l) layer and $h_{\mathrm{diff}}$ is the vertical distance between layers. The direct effect of wind stress is usually confined to surface water, although indirect effects such as Ekman transport or overturning and the formation of eddies can facilitate currents and mixing at deeper depths . The two-dimensional correlation coefficient $R_{\mathrm{WC}}$ between the wind and ocean currents was calculated using and the complex representations of the wind stress and de-tided current velocity vectors (${\mathit{\tau }}_{\mathrm{c}}$ and ${\mathit{u}}_{\mathrm{c}}$, respectively) as follows: E2$${\mathrm{R}}_{\mathrm{WC}}={\displaystyle \frac{⟨{\mathit{\tau }}_{\mathrm{c}}^{*}{\mathit{u}}_{\mathrm{c}}⟩}{⟨{\mathit{\tau }}_{\mathrm{c}}^{*}{\mathit{\tau }}_{\mathrm{c}}{⟩}^{\frac{\mathrm{1}}{\mathrm{2}}}⟨{\mathit{u}}_{\mathrm{c}}^{*}{\mathit{u}}_{\mathrm{c}}{⟩}^{\frac{\mathrm{1}}{\mathrm{2}}}}},$$ where $⟨\mathrm{\dots }⟩$ gives the normalized inner product of the vectors and ${}^{*}$ annotates the complex conjugate. We allow time lags of up to 15 h to account for the gradual and indirect effects of wind stress on the ocean currents. Both properties were estimated throughout the valid current velocity profile, but this was only done only down to 80 m depth due to the 8 m vertical distance between the defined layers used in the bulk velocity shear calculation.

Code and data availability

All data presented in this paper can be obtained upon request to the authors and are also available in the platform Open research Data at the University of Tromsø – The Arctic University of Norway (https://dataverse.no/dataverse/uit, ). All computer code used can be obtained upon request to the corresponding author.

Author contributions

The study was conceptualized by KOD, BF, AS, PL, and PJ. The data were curated by KOD, BF, and MM. Formal analysis was done by KOD, BF, and MM. Funding was acquired by BF. Investigation was done by KOD, BF, and AS. Methodology was implemented and developed by KOD and BF. The project was administrated by BF and AS. Resources was obtained by BF. Any software used was developed by KOD. The project was supervised by BF. Visualizations were made by KOD and MM. KOD wrote the manuscript (original draft), and KOD, BF, AS, PL, PJ, and MM contributed to reviewing and editing the manuscript.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

We thank the crew of R/V Helmer Hanssen for their assistance during the deployment (CAGE 15-3) and recovery (CAGE 16-4) cruises. This study is a part of CAGE (Centre for Arctic Gas Hydrate, Environment and Climate), Norwegian Research Council (grant no. 223259). We thank Nicholas Warner for proofreading the article. We the thank two anonymous reviewers and the editor for their constructive feedback on the manuscript.

Financial support

This research has been supported by the Norwegian Research Council (grant no. 223259).

Review statement

This paper was edited by Mario Hoppema and reviewed by two anonymous referees.