Seismo-stratigraphic framework

The seismo-stratigraphy of the Nyegga area is characterized by the presence of six pronounced reflections, which are well-defined in echosounder data throughout the study area43. We traced these reflection in a dense network of echosounder profiles collected with the TOPAS echosounder system on-board RV G.O. Sars by University in Bergen in 2004 and a Parasound echosounder system on board of RV Meteor by Geomar in 2012 (Fig. 6). The stratigraphic section defined by the deepest and shallowest interpreted seismic reflections consists of spatially homogenous, laminated sediments and shows no unconformities indicating phases of erosion. The thickness of specific intervals within this stratigraphic section varies depending on the proximity to glacigenic depositional centres to the West7. The seismo-stratigraphy represents a continuous record of the sedimentation history of the study area and allows constraining a stratigraphic framework for the Nyegga pockmark field.

Fig. 6 Reconstruction of the seismo-stratigraphic framework of the Nyegga area. a Map showing the location of the echosounder profiles used for the seismo-stratigraphic framework, pockmarks, sediment core locations, model area, and area with mapped BSR (Bottom-simulating reflector).The map was created with IHS Kingdom Suite and data sourced from bathymetry grids created by Norsk Hydro AS55. b–g Sediment thickness maps for different intervals of the seismo-stratigraphic framework Full size image

Radiocarbon age calibration

The calibration of radiocarbon dates from marine sources need to be corrected for the global and local reservoir effects. The global reservoir effect is set to 405 years within the Marine13 calibration curve44. The local reservoir effect depends on time-varying oceanographic parameters, which are not well-constrained for the study area during the analysed period. Numerical modelling of the local reservoir effect for the North Atlantic suggests values between 200 and 600 years45, while the comparison of terrestrial and marine 14C reservoir ages revealed variations of the local reservoir effects between 100 and 400 years for the Younger Dryas46. To acknowledge this uncertainty, we used a local reservoir effect of 400 years with an uncertainty of ± 200 years for the calibration tool Calib7.1 using the Marine13 calibration curve for all radiocarbon based dates in this study43,47 to make them self-consistent.

Age-depth-model

The age-depth model builds on integrating radiocarbon dates from the five sediment cores MD99-2289, MD99-2291, HM119-02GC, HM119-03GC and GS08-155-47GC (Fig. 7; Table 1). Sediment core MD99-2291 is located close to the depositional centre in the West of the study area and has a high-resolution record of the shallow to intermediate interval of the seismo-stratigraphic framework, while sediment core MD99-2289 is located more than 60 km away close to the Storegga Slide side-wall and provides age information for the intermediate to deep interval of the seismo-stratigraphic framework. Sediments cores HM119-02GC, HM119-03GC and GS08-155-47GC targeted the shallowest seismic reflections of the seismo-stratigraphic framework.

Fig. 7 Sedimentation rate reconstruction based on combined sediment cores from the Nyegga area. a–d Echosounder profiles showing the locations of radiocarbon dates from the different sediment cores. The intervals of the seismo-stratigraphic framework are coloured. e Locations of the radiocarbon dates in the combined core. f Age-depth model from the combined core (calibrated ages) with the 2 sigma standard deviation in grey Full size image

Table 1 Radiocarbon dated samples used for the age-depth model and sedimentation history reconstruction. Radiocarbon ages are not corrected for reservoir effect Full size table

The analysis of the seismo-stratigraphy of the modelled area suggests that the shape of the sedimentation curve is approximately constant throughout the study area and only the magnitude of sedimentation varies depending on the proximity to the depositional centres. This implies that the relative depth of a specific date within an interval of the stratigraphic framework at its coring location should be approximately constant within the stratigraphic framework. Therefore, we transferred the depth from cores MD99-2289, MD99-2291, HM119-03GC and GS08-155-47GC, to the coring location of HM119-03GC with the following equation:

$${\rm{depth}}_k^0 = {\rm{top}}_i^0 + \frac{{{\rm{top}}_i^j - {\rm{depth}}_k^j}}{{{\rm{base}}_i^j - {\rm{top}}_i^j}} \times \left( {{\rm{base}}_i^0 - {\rm{top}}_i^0} \right),$$

whereas \({\rm{depth}}_k^0\) is the transferred depth of sample k from core site j, which lies in coring depth \({\rm{depth}}_k^j\) within the stratigraphic interval i. The radiocarbon ages, including their measurement errors and their transferred depth within the virtual core were the input for the Bayesian statistics based age-depth modelling script Bacon48, which provided a high-resolution age-depth model for the Nyegga pockmark field (Fig. 3d; Fig. 7). The cores MD99-2289, MD99-2291 and GS08-155-47GC provide overlapping age-depth information for the section between 8 and 33 m of the combined core (7 f). The age-depth information from MD99-2289 and MD99-2291 reveal very similar sedimentation rates. The age-depth information provided by GS08-155-47GC is in good agreement with the trend from MD99-2289 and MD99-2291. The shallow-most radiocarbon age from MD99-2291 is about 1500 years younger than the shallow-most radiocarbon age from GS08-155-47GC. However, the Bayesian age-depth modelling prefers the age-depth information from GS08-155-47GC leading to a smoother sedimentation curve than following the shallow-most radiocarbon age from MD99-2291. The “too young” radiocarbon age of the shallow-most radiocarbon age from MD99-2291 is most likely the result of reworking as the sample was taken directly below an erosional surface (between purple and red intervals in 7b).

4D sedimentation reconstruction

The local age-depth model and the seismo-stratigraphic framework were then integrated into 4D sedimentation history reconstruction. The age-depth-model provided absolute age information for each seismic reflections and the local sedimentation rate evolution. The sedimentation rate constrained by the age-depth model was then adjusted for every location within the seismo-stratigraphic framework by scaling the sedimentation rate for each interval and each location of the seismo-stratigraphic framework with the thickness of that interval at a specific position divided by the thickness of the interval in the combined core. This approach was not possible for the stratigraphic record younger than 16.2 ka BP and for the western-most part of the Nyegga pockmark field, where buried glacigenic debris lobes were deposited in between the laminated sediments.

Analyses of carbonate crusts sampled from Nyegga pockmark area

A carbonate crust sample was collected during a dive on submersible MIR (on board RV Akademik Mstislav Keldysh during GEOMAR cruise 40 in 1998). The sample was analysed for C and O isotope composition by dissolving three finely powdered sub-samples with 100% H 3 PO 4 at 72 °C in the Carbo Kiel IV preparation line and thereby releasing CO 2 , which was isotopically analysed using a Finnigan MAT 253 mass spectrometer. 18O/16O and 13C/12C ratios were determined against a working calcite standard and are expressed in the common δ-notation relative to Vienna Pee Dee Belemnite (VPDB). The uncertainty of this sample set accompanying standard measurements (n = 8) was 0.16‰ (2σ) for δ18O and 0.04‰ (2σ) for δ13C. Table 2 shows the results with their individual measurement uncertainties on 2 SD level. No further corrections for different acid fractionation factors of calcite and dolomite were applied. Mineral composition of fine-grained carbonaceous material was determined by using X-ray diffraction analysis (Philips PW 1820, Co-Kα, 0.02°/s). X-ray diffraction patterns were interpreted by using the XPowder® software. Identified carbonate minerals (i.e. calcite and dolomite) were quantified by Rietveld refinement. In addition, we applied U–Th geochronology methods for methane-derived authigenic carbonates to constrain the precipitation age of the carbonates49,50. However, the analysis indicated different formation ages for the calcitic and dolomitic dominated subsamples of the same carbonate sample. A too high contribution of detrital Th in both carbonate phases prevented deducing reasonable isochron-based age information. The finding of different formation ages for the different carbonate phases is in agreement with the light stable isotope findings. The 14C data from the carbonate crust sample Mie 1-3729-1 represent a rather robust maximum formation age of 18,000 years BP (15,710 14 C). This maximum age interpretation is based on the assumption that any carbon contributed from the sediment system below would be relatively depleted in 14C by ongoing decay during burial, resulting in apparent higher ages.

Table 2 Mineralogical and geochemical characterisation of carbonate concretion Full size table

Finite difference temperature modelling

The temperature profile modelling builds on Fourier’s law for heat flow and the assuming homogenous thermal properties of the deposited glacigenic sediments, which allows using the thermal diffusivity κ (4.2 × 10−7 m2/s; ref. 15) for calculating the heat flow: \(q = \kappa \,

abla T\). We transferred this simple differential equation in a 1D finite difference algorithm:

$$T_i^{n + 1} = \frac{{\kappa {\mathrm{\Delta }}t}}{{{\mathrm{\Delta }}x^2}}\left( {T_{i + 1}^n - 2T_i^n + T_{i - 1}^n} \right) + T_i^n,$$

whereas Δt is the temporal discretization of 1 s and Δs the spatial discretization of 0.1 meter. The simulation started with equilibrium temperature profile with a length of 300 m following the local temperature gradient of 55 °C/km (ref. 15). This profile was then extended at the top by the amount of sediments deposited within a time step of 200 years according to the 4D sedimentation history reconstruction. The sediments had the bottom water temperature of −1 °C, which was also used as boundary condition, the lower boundary was defined as a gradient and the heat flow adjustment was simulated over 200 years using an implicit approximation. The resulting temperature profile was then the starting model for the following time step.

Gas hydrate stability calculation

The gas hydrate stability calculations followed the approach by Tishchenko22, which build on thermodynamic stability and solubility calculation for methane hydrates. The dissociation temperature for specific depths below sea floor z is calculated following the empirical function: T diss (z) = a + b × log(P(z)) + c × (d − log(P(z)))3, whereas T diss is the dissociation temperature, P is the hydrostatic pressure and a, b, c, d are empiric parameters calculated by fitting the stability curve of methane hydrates for a given salinity of 3.5% and pure methane (a = −25.1391, b = 7.9832, c = −0.0959 and d = 6.0687). The hydrostatic pressure at the sea floor and within the top 300 m within the sediments for a given point at a given time was calculated analytically by implementing the sea level fluctuation, subsidence and the sedimentation reconstruction. The T diss profile was then compared with the temperature profile from the finite difference heat flow simulations to determine the BGHSZ, which is defined by the shallowest depth in the modelled temperature is higher than T diss .

Sensitivity analysis

In order to test the robustness of our simulations, we performed sensitivity analyses for different sedimentation rates, bottom water temperature profiles, the subsidence rates affecting the local sea level curve, bulk thermal conductivity values and the geothermal gradients (Fig. 8; the input parameter and simulations results used in the main part of the manuscript are marked with dashed black lines). We tested the impact of bottom water temperature on the gas hydrate dynamics by applying different temperature evolution profiles. These simulations included scenarios with constant temperatures of −1 °C (dashed black line in Figs. 8a) and −2 °C (yellow line) and repeatedly fluctuating temperatures between −1.75 and −2.25 °C (orange line). These input parameters result in very similar curves for the dissociation rate of gas hydrates. In addition, we performed simulations with four potential warming scenarios with an absolute warming of ~0.5 °C (green line) and ~1 °C (red, light blue and dark blue). All simulations result in similar shaped gas hydrate dissociation curves with a peak of dissociation after 18,000 years before present. The absolute thickness of gas hydrate dissociation at the end of the simulations lies between 35 and 47 m, which shows that bottom water temperatures had only a secondary impact on gas hydrate dynamic during the LGM at Nyegga. The solution space of the Bayesian age-depth reconstruction allows different sedimentation rate reconstructions (Fig. 8b). In addition to the statistically most likely reconstruction (dashed black line), we tested the sedimentation rate reconstructions, which bound the solution space (blue and yellow lines; marked in grey in Fig. 7f). These sedimentation rate reconstructions result in very similar gas hydrate dissociation rates and absolute dissociation curves. The analysis of different local subsidence rate (Fig. 8c), which influence the local sea level curve reveals that subsidence had only a minor impact on gas hydrate dynamics, even when neglecting subsidence (yellow line) or assuming a twice as high subsidence rate (blue line). The analyses of the thermal diffusivity of the sediments and the geothermal gradients reveal that the effect of sedimentation on gas high dynamics becomes more important for sediments with a high thermal diffusivity and in areas with high geothermal gradients (Fig. 8d, e). However, the timing of enhanced gas hydrate dissociation remains stable for various thermal property values. The sensitivity analysis revealed that the gas hydrate dissociation rate and the absolute value of dissociated gas hydrates are robust and the observed gas hydrate dynamics can be observed using a wide range of parameters.

Fig. 8 Results of the sensitivity analyses for different simulation input parameters. a Bottom water temperature profiles and the resulting gas hydrate dissociation rates and absolute gas hydrate dissociation curves. b Sedimentation rates: minimum, maximum and best-fit sedimentation rate reconstructions and the resulting gas hydrate dissociation rates and absolute gas hydrate dissociation curves. c Subsidence rate as well as their impact on the local sea level anomaly and the resulting gas hydrate dissociation rates and absolute gas hydrate dissociation curves. d Thermal diffusivity and the resulting gas hydrate dissociation rates and absolute gas hydrate dissociation curves. e Geothermal gradient and the resulting gas hydrate dissociation rates and absolute gas hydrate dissociation curves56 Full size image

Seepage potential calculation

The seepage potential can be calculated by multiplying the average thickness of dissociated gas hydrates (40 m from our simulations) with the area of gas hydrate occurrence A (model area: 540 km2), the gas hydrate concentration (3–12%, ref. 35), the sediment porosity Φ (0.38–0.76, ref. 51), the gas hydrate expansion factor under atmospheric pressure conditions (164, ref. 34) and the density of methane under atmospheric conditions ρ CH4_atmospheric (~ 0.656 kg/m3):

$${\mathrm{SP}} = H_{\rm{hydrate}} \times A \times s_{{\rm{hydrate}}} \times \Phi \times 164 \times \rho _{{{\rm{CH4}}}_{\rm {atmospheric}}}$$

The resulting seepage potential for the modelled area lies between 26 and 212 Mt.

Data availability

The data sets analysed during the current study are available from the corresponding author on reasonable request.