Deglacial bottom water warming intensified Arctic methane seepage in the NW Barents Sea

Changes in the Arctic climate-ocean system can rapidly impact carbon cycling and cryosphere. Methane release from the seafloor has been widespread in the Barents Sea since the last deglaciation, being closely linked to changes in pressure and bottom water temperature. Here, we present a post-glacial bottom water temperature record (18,000–0 years before present) based on Mg/Ca in benthic foraminifera from an area where methane seepage occurs and proximal to a former Arctic ice-sheet grounding zone. Coupled ice sheet-hydrate stability modeling shows that phases of extreme bottom water temperature up to 6 °C and associated with inflow of Atlantic Water repeatedly destabilized subsurface hydrates facilitating the release of greenhouse gasses from the seabed. Furthermore, these warming events played an important role in triggering multiple collapses of the marine-based Svalbard-Barents Sea Ice Sheet. Future warming of the Atlantic Water could lead to widespread disappearance of gas hydrates and melting of the remaining marine-terminating glaciers. Phases of high bottom water temperature in the northwestern Barents Sea caused repeated destabilization of methane gas hydrates since the last glacial, according to a foraminifera Mg/Ca bottom water temperature record and hydrate stability modelling

T he Arctic has experienced dramatic changes during recent decades in response to climate warming. In particular, the Barents Sea has shown signs of an ongoing "Atlantification", with an increased inflow of warm Atlantic Water (AW) and a significant loss of sea ice 1 . The incursion of warm AW onto the Barents Shelf can also destabilize buried gas hydrates-frozen compounds of water and methane-releasing greenhouse gases into the water column [2][3][4] . Although methane released from the seafloor is microbially consumed in the sediment by anaerobic oxidation of methane 5,6 or dissolved/consumed by aerobic methane oxidation in the water column 7,8 , a major methane venting event induced by dissociation of gas hydrates could amplify the effects of the current ocean acidification 9 and, if escaping to the atmosphere could lead to further warming.
Past periods with the subsurface inflow of warm AW have been recorded in sediment cores from the Nordic Seas during Heinrich Stadials (HS) [10][11][12] . These are millennial-scale events related to the extensive discharge of icebergs and meltwater into the North Atlantic that occurred in the last glacial period 13 . During the Last Glacial Maximum (ca. 26,000-19,000 calendar (cal) years before present (BP)), the Svalbard-Barents Sea Ice Sheet (SBIS) reached its maximum extent, covering the entire Barents Sea shelf 14 . Paleoceanographic studies 15,16 and ice-ocean modeling 17,18 have suggested that subsurface warming was a major factor in driving the collapse of the marine-based SBIS. The western Barents Sea margin during the last deglaciation is thus a highly relevant analog for the ongoing subsurface warming-induced retreat of the Western Antarctic Ice Sheet (WAIS) 19 and marine-terminating glaciers in Greenland 20 .
We have reconstructed bottom water temperatures using samples of core HH18-1059GC (hereafter core 1059), recovered from Storfjordrenna (Storfjorden Trough) in the northwestern Barents Sea at a water depth of 382 m (Fig. 1). Warm Atlantic water derived from the North Atlantic Current flows from the southern Nordic Seas to the Fram Strait along the western Svalbard margin as the West Spitsbergen Current 21 (Fig. 1). A branch of the Atlantic Water enters Storfjordrenna and Storfjorden in the eastern part. The East Spitsbergen Current, a cold polar current, enters Storfjordrenna from the Arctic Ocean and the northeast Barents Sea and flows in the western part of Storfjorden and along the inner Spitsbergen shelf (Fig. 1). During the late winterspring season, the AW generally flows above 200 m in the upper part of the water column with a temperature ca. 3°C and a salinity above 34.95 psu, below a polar, low-salinity, surface water layer 21 . This is because (depending on the prevailing wind directions and sea-ice conditions), a polynya may form in late winter and spring in Storfjorden resulting in the formation of brine-enriched shelf water (−1.9°C, >34.8 psu) that flows along the bottom of Storfjordrenna on its way to the shelf edge 21 . In summer-autumn, the AW occupies the entire water column. In July 2018, when core 1059 was retrieved, the AW reached the seafloor with a temperature of 2.5°C (Fig. 1C).
Previous attempts to estimate paleo-BWT in Storfjordrenna have involved converting δ 18 O of benthic foraminifera 15 into BWT by assuming constant seawater δ 18 O similar to modern values 22 or calculating BWT using benthic foraminiferal transfer functions 16 (Fig. S2). The latter method quantifies BWT as an average of the whole foraminiferal assemblage and might therefore be biased towards species-specific ecological preferences other than temperature (e.g., food and oxygen availability, salinity, and water depth). Here, we quantify BWT based on Mg/Ca in the benthic foraminiferal species Cassidulina neoteretis, in  order to more accurately constrain past temperature changes in the western Barents Sea. In addition, the reconstructed BWT is used to investigate the past dynamic of the gas hydrate stability zone (GHSZ) by modeling variations in its thickness to provide an oceanographically constrained history of methane-venting from the northwestern Barents Sea area. Our record spans the early deglaciation (since the Late Glacial 18,000 cal years BP) to the late Holocene, including the cold atmospheric periods Heinrich Stadial 1 (HS1; 17,500-14,600 cal years BP) and Younger Dryas (YD; 12,800-11,700 cal years BP) and the warm interstadial phases Bølling-Allerød (14,600-12,800 cal years BP) and the Holocene (since 11,700 cal years BP) 23 (Fig. 2). The core site is located close to an area of methane release from gas hydrate mounds ("Pingo Area") in the northwestern Barents Sea 22 , allowing correlation and comparison between methane-influenced records 24 and records from unaffected areas (e.g., core 1059; see "Methods" section). The gas hydrate stability zone outcrops at this water depth in the outer Storfjordrenna 25 in an area also strongly affected by the inflow of the warm Atlantic Water (Fig. 1).

Results and discussion
Bottom water temperatures during the last 17,500 years. The BWT in Storfjordrenna varied between 1.5 to 5.5°C since the start of HS1 at 17,500 years BP (Fig. 3A). The warmest BWT occurred during HS1 and the following deglaciation (Bølling-Allerød interstadials, 14,600-12,800 years BP), with an average of 4.4 ± 1°C (Fig. 3A). Our results generally agree with previous BWT estimates by transfer functions from this area 16 (Fig. S2). Quantifications of BWT based on benthic foraminiferal Mg/Ca from intermediate water depth (down to at least 1273 m water depth) in the southern Nordic Seas 11 and the western Svalbard margin 12 show similar trends with a warming of the bottom water of up to 5.5°C during HS1. During HS1 and older HS events during the last glacial period, the AW could flow northwards beneath a polar meltwater layer with an extended sea-ice cover and a strong halocline preserving most of its heat 10,11 . The meltwater came from the melting of icebergs released from the Northern hemisphere ice sheets and the formation of cold deep water stopped or became very weak due to reduced oceanic convection in the Nordic Seas.
The end of HS1 is marked by a progressive increase in BWT up to 5.3 ± 1°C, shortly before the beginning of the Bølling-Allerød interstadials (Fig. 3A). The transition to a warmer climate is characterized by laminated fine clays deposited from meltwater plumes 26 (Fig. 2). The benthic foraminiferal fauna is dominated by Elphidium excavatum 15 , a benthic foraminifera species that tolerates turbid meltwater, highly variable environmental conditions, and low salinities 27,28 indicating rapid ice retreat and meltback of the SBIS 16,29 .
During the YD, the BWT was relatively low (2.5 ± 1°C) and was closer to modern values (Fig. 2). The YD was linked to a slowdown of the Atlantic Meridional Overturning (AMOC) 30 but differs from other HS events as it occurred during warm conditions 10 . The YD cold spell coincided with the deglaciation of the Arctic Ocean sending sea-ice-loaded cold water into the Nordic Seas and North Atlantic 31 . Seasonal sea-ice cover and an increased seasonal brine formation throughout most of the YD 32 , probably kept the seafloor relatively cool similar to today (Fig. 3A).
Bottom water temperature progressively increased up to ca. 5°C by the early Holocene between 11,700 and 9000 years (Fig. 3A). During this period, benthic foraminiferal assemblages show an increase in relative abundances of Cibicides lobatulus and Melonis barleeanus supporting the presence of stronger advection of AW and a reduced cover of sea ice 15 . During the Mid-Late Holocene, the BWT stabilized with an average of 3.3 ± 1°C, except for a short warming event of around 6000 years. Bottom water temperature and gas-hydrate dissociation. Bottom water temperature is an effective regulator of the stability of sub-seafloor gas hydrates, with minor fluctuations on a seasonal basis directly impacting the flow of greenhouse gases from the seabed to the water column 22 . Using our BWT time series, we modeled changes in the thickness of the GHSZ for the past 18,000 years (see "Methods" section). The measured benthic foraminiferal species, C. neoteretis, thrives under the influence of AW 27,28 , and therefore, the thickness of the GHSZ presented in this work should be considered as a narrowest limit (see "Methods" section).
During the late glacial from 18,000-17,500 years BP, the rapid removal of the glacial overburden followed by the intrusion of the warm AW caused the ca. 200 m thick GHSZ to disappear entirely within ca. 1000 years (Figs. 3B, 4, 5A, B and Fig. S4). The combined effects of glacial uplift (that decreases the pressure) and sea-level rise (which exceeds the uplift and increases the pressure) heighten the gas hydrate stability (Fig. 4). The most likely factor in our model to drive the rapid changes in the postglacial hydrate stability is thus the recurrent increases in BWT in the Pingo Area (Fig. 4). Indeed, the GHSZ continued to outcrop at the sediment surface in the area during the HS1 and warm Bølling and Allerød interstadials with high BWT. The GHSZ reappeared and thickened at the Allerød-YD transition following a decrease in BWT (Fig. 3B). The climatically unstable periods of the YD and earliest Holocene with highly variable BWT resulted in rapid fluctuations in the thickness of the GHSZ (0 to between 50 and 80 m) (Fig. 3B). Thereafter, it disappeared again during the Holocene Thermal Maximum, before increasing to a thickness of 60-75 m beginning from 4500 years BP in the late Holocene, when BWT stabilized to modern values (Figs. 1C and 3A). Seismic data from the area shows that the base of the modern GHSZ (i.e., the bottomsimulating reflector) occurs today between 85 and 150 m below the seafloor 25 , validating the results of our model.
We have correlated nearby core CAGE15-2-920GC (hereafter core 920) from a hydrate pingo to our core 1059 ( Fig. 1 and Fig.  S3). Core 920, analyzed for its content of archaeobacteria and their 13 C signals, indicates three major episodes of methane venting 24 . Based on our correlation and placing core 920 into our age model (see "Methods" section), increased release of gas occurs in the late HS1, during the Bølling interstadial and Allerød interstadial.
In the Pingo Area methane is released through several gas hydrate mounds connected to the Hornsund fault system in Storfjordrenna 25 . The record from core 920 shows when the methane flux in the Pingo Area in Storfjordrenna increased during the deglaciation. There is no signal in core 920 for most of HS1, which may have been expected given the rapid disappearance of the GHSZ and the following high BWT (Fig. 3B, C). Data from the SW Barents Sea indicated a maximum in seepage immediately after the retreat of the SBIS 18,000-16,000 years BP 33 . The venting right after depressurizing from the ice burden may not have been recorded because (1) it was a very rapid and violent event, (2) the gas escaped through cracks in the gas phase, or (3) it escaped through cracks and faults elsewhere. Indeed, studies from the SW Barents Sea 33,34 also indicated an increase in methane venting during the Bølling and Allerød intervals supporting the hypothesis of the protracted nature of methane release from gas hydrate dissociation during the deglaciation over millennial timescales, where a high flux of methane follows the Significance of bottom water temperature changes for methane release and ice retreat. Our BWT record for the western Barents Sea provides quantifiable evidence for bottom warming events in the northwestern Barents Sea since the Last Glacial Maximum, providing a direct means of comparison of its impact on the ice sheet and hydrate dissociation. The core site in Storfjordrenna is located in front of a former major ice stream of the SBIS 35 , which had retreated to the central part of the shelf at least before 18,150 cal years BP (16,750 ± 110 14 C years 15 calibrated with Normarine18 29 ) (see "Methods" section). The collapse and retreat of this ice stream from the shelf edge have been linked to oceanographic forcing 36 , enhanced by the seafloor geometry and substrate 37,38 . The three major periods of high BWT (with temperatures ca. 4-5°C) identified after 18,000 years BP correlate closely with major postglacial retreat phases of the SBIS in Storfjordrenna and Storfjorden 16,38 (Fig. S4).
The early retreat of the SBIS from the outer part of Storfjordrenna was probably accelerated by the presence of warm (4.4 ± 1°C) AW flowing eastwards during the late glacial and HS1 15 (Fig. 3A and Fig. S4) followed by an increased release of methane from the seafloor 24 (Fig. 3C). The contemporaneous and rapid reduction in thickness of the GHSZ in Storfjordrenna has previously been attributed to decompression associated with deglaciation of the SBIS 39 . Our data suggest that the high BWT most likely also played an important role, driving the retreat of the SBIS, and preconditioning and accelerating the postglacial thinning of the GHSZ much faster than previously suggested 22 .
We observe a decrease in seawater δ 18 O following the peak in BWT during HS1, potentially indicating a freshening of the bottom waters (Fig. 3A, E). This indicates that the advection of warm AW into the western Barents Sea probably resulted in basal melting of the SBIS, forcing a rapid retreat of the ice stream in Storfjordrenna and a significant supply of freshwater (Fig. 5B) 15,16 . Similar chronologies of the retreat patterns of the Storfjordrenna Ice Stream and north Norwegian Ice Streams likely suggest that the retreat was controlled by a common forcing (i.e., atmosphere and/or oceanic control) rather than by local factors only 16 (Fig. 5A-D). Subsurface melting has been suggested by model studies to be the major control on the retreat of the marine-based SBIS 17 , in contrast to the minor role played by atmospheric forcing and sea-level rise 18 . Our BWT, with three main phases of AW warm pulses, agrees with the stepwise retreat pattern of the SBIS previously documented in Storfjordrenna 37 and the SW Barents Sea 40 .
The BWT increase at the beginning of the Holocene between 11,300 and 9000 years BP, may have driven the final retreat of the SBIS in Storfjorden 16 and the disappearance of the GHSZ (Fig. 3B). By the Mid-Late Holocene, modern BWT values were reached (Fig. 3A). The lower summer BWT and the increased pressure related to sea-level rise during the Holocene probably thickened the GHSZ again (Figs. 3B and 5F). Seasonal variations in bottom water temperature today and during the late Holocene 32,41 were pronounced due to increased brine formation in Storfjorden creating dense, cold outflow water via Storfjordrenna (winter with strong brine formation: T = −1.3°C, S = 35.3 psu 15 ; summer T = 2.5-3°C, S = 34.9-35 psu; Fig. 1C). This strong seasonal variability is likely to accelerate the destabilization of the GHSZ and cause variable seasonal patterns of methane seepage on Arctic continental shelves [2][3][4] (Fig. 5F).
A modeling exercise on the fate of the GHSZ in the southwestern Barents Sea, under warming ocean conditions since 1960, indicates that the inflow of Atlantic water played a major role in the thinning of the GHSZ between 1985 and 2010 42 . Furthermore, a linear increase of 1°C of BWT from 2010 to 2060 would cause a further thinning of the GHSZ, allowing shallow (<80 mbsf) methane hydrates to dissociate and release between 1 and 8 Gt of carbon into the ocean 42 . Under a very high baseline scenario of greenhouse gas emissions (RCP8.5) 43 , CMIP5 climate predictions 44 of BWT show that the temperature could increase up to 5°C by the end of the 21st century in the western Barents Sea at water depths of 390 m 45 , similar to the estimated BWT in Storfjordrenna during HS1 and the Bølling-Allerød interstadials (Fig. 3).
The mechanisms underlying the intrusion of warm AW into the western Barents Sea were different in the past (i.e., when AW could flow with minimal heat loss beneath a thick halocline 10,11 ). However, our study shows that the rate of thinning of the GHSZ in the Pingo Area was at least six times faster than previously estimated, occurring over a period of 1000 years compared to the suggested 6000 years 22 . The high BWT during HS1 (Fig. 5B) drove the rapid thinning of the GHSZ, and probably the enhanced gas hydrate dissociation and methane seepage reconstructed for the Bølling-Allerød interstadials (Figs. 3A-C, and 5C, D). These results highlight the important role of BWT changes in  46 . Records of past BWT at different time scales underline the risk that the current "Altantification" process in the Arctic bears in triggering an increase in seepage of methane into the ocean, causing an increase in ocean acidification and potentially amplifying the effects of current climate change if reaching the atmosphere.

Methods
Core handling. Core HH18-1059GC (core 1059) (76°06.117′N; 15°58.077′E, 382 m water depth) was retrieved from the southwestern Barents Sea during a cruise in July 2018 with R/V Helmer Hanssen (Fig. 1). The 4.15 m long core was split into 1-m sections, capped and taped at both ends, and stored at 4°C right after retrieval. Prior to opening, the core sections were X-rayed and logged with a GEOTEK 7.9 Multi-Sensor Logger. The core was split longitudinally into two halves. The work half was color imaged with a Jai L-107CC 3 CCD RGB line scan camera installed on an Avaatech XRF. The archive half was scanned at 10 and 30 kV on an Avaatech XRF for bulk element ratios.
Thereafter, the work half was sampled in 1-cm thick slices. Samples were weighed, freeze-dried, and weighed again. Samples were selected following the main focus of our study: (a) from 3.9 to 2.9 m samples were selected every 2 cm, except at 3.66-3.64 m and 3.83-3.80 m where every cm was selected; (b) from 2.9 to 0.1 m samples were selected every 10 cm; (c) from 0.1 to 0 m, every 2 cm was selected. The samples were wet-sieved over 63 μm, 100 μm, and 500 μm mesh-sizes. The residues were dried at 40°C, weighed and the weight percent of each grain size was calculated.
Lithology, radiocarbon dating, and construction of the age model. The lithological log is based on visual examination, together with the records of X-ray scanning, grain-size distribution, and magnetic susceptibility (Fig. 2).
Nine AMS-14 C dates were acquired on samples of the planktic foraminiferal species Neogloboquadrina pachyderma and bivalve samples at the 14 Chrono Centre, Queens University, Northern Ireland (Table 1 and Fig. 2). Two additional dates from the upper and lower boundaries of the laminated layer with a low concentration of foraminifera in core 1059 were added from correlation to dates of the same laminated layer from the well-dated nearby located cores JM02-460 15 and HH15-1282GC 16 . Dates older than 12,436 ± 66 14 C years were calibrated using the Normarine18 calibration curve which uses variable reservoir ages ranging from 420 14 C years to 1620 14 C years prior to the Bølling-Allerød warming 29 . Dates younger than 10,827 14 C years were calibrated using the Marine20 calibration curve 47 . , which corresponds to the lower limit of the GHSZ. The uppermost limit (top of gas hydrate occurrence zone) is at the seafloor 22 . Inflow and relative temperature of the Atlantic water is shown in red, where stronger red indicates a relatively higher temperature. A thick white arrow marks the position of core 1059. Age for each interval is shown as calibrated ages before present (BP). Note that the water depth is not drawn to scale. See Fig. 4 for details.
Individual samples were calibrated using the CLAM 2.3.5 package 48 ( Stable isotope analyses. Oxygen and carbon isotopes were analyzed on pristine tests of the planktic foraminiferal species Neogloboquadrina pachyderma and on the benthic foraminiferal species Cassidulina neoteretis, in the size fraction 150 to 250 µm. The measurements were performed on 10 to 20 specimens (only 8 samples out of 117 samples contained less than 10 specimens) of each species using a Thermo Scientific MAT253 IRMS and Gasbench II at the Department of Geosciences, UiT The Arctic University of Norway, Tromsø. The precision of the instrument is 0.1‰ for oxygen and carbon isotopes and the results are reported versus the in-house Vienna Pee Dee Belemnite standard. The measured δ 13 C values of C. neoteretis range from −1.77‰ to −0.51‰ (Fig. 3F). These values are within the expected range of δ 13 C values for C. neoteretis 50 indicating that core 1059 has not been affected by methane seepage and any associated diagenetic coatings. We, therefore, concluded that core 1059 was suitable for benthic foraminiferal Mg/Ca-based reconstructions of BWT.
Mg/Ca analysis and bottom water temperature calculation. A total of 15 to 30 pristine tests of the dominant benthic foraminiferal species in the core Cassidulina neoteretis were picked for Mg/Ca analyses from the 150 to 250 µm size fraction. The oxidative-reductive approach was used to clean the samples prior to the analysis [51][52][53] . This cleaning approach includes clay removal, reductive cleaning with hydrous hydrazine, oxidative cleaning with oxygen peroxide, and weak acid leach. After cleaning the samples were dissolved in HNO 3  Elemental ratios Mn/Ca, Fe/Ca, and Al/Ca were used in combination to evaluate contamination, and two samples indicating potential contamination were excluded. Thereafter, we applied Grubb's test to identify any outlier in the Mg/Ca data and two samples were excluded. The remaining samples showed a low correlation between Mn/Ca, Fe/Ca, and Mg/Ca (r 2 = 0.14 for Mn/Ca-Mg/Ca; r 2 = 0.06 for Fe/Ca-Mg/Ca; aluminum concentrations were below the detection limit indicating no clay contamination; Fig. S5 and Fig. S6).
Mg/Ca values were converted into temperature using the calibration equation for C. neoteretis from Kristjánsdóttir et al. 55 : Mg=Ca ¼ 0:864 ± 0:07 expð0:082 ± 0:020 BWTÞ The calibration error (±0.62°C 55 ), the analytical error (±0.017 mmol/mol equivalent to 0.18°C using Kristjánsdóttir et al.'s calibration), and the standard deviation of the replicates of four samples (±0.078 mmol/mol equivalent to 0.78°C) were used to calculate the error, that is the result of the squared root of the sum of the squared individual errors. The analytical error was calculated as two times the mean standard deviation of the repeated measurements on the in-house standard solution. This gave an estimated propagation error of ±1.01°C for C. neoteretis. When an average BWT was presented in the main text, the error was calculated taking into account the propagation error of the measurement (±1.01°C) and the standard deviation of the mean, and the higher "error" value was chosen.
Considering the affinity of C. neoteretis for Atlantic Water (AW) 27,28 and the hydrography of Storfjordrenna, our BWT record probably represents the temperature range of this water mass. Today, the AW has a temperature of ca. 3°C in Storfjordrenna and is representative of years with low brine production (in years with strong brine formation temperature may reach down to −1.3°C and salinity increase to 35.3) 21 . We, therefore, speculate that in the past, during periods with strong stratification, such as the early deglaciation, with low seasonality and no or low brine flow 32 , the AW reached the deepest part of Storfjordrenna almost all year and therefore that our BWT represents an annual mean. During periods of strong seasonality (e.g., Late Holocene 15 ), our BWT would represent seasons of AW inflow (i.e., modern summer-autumn conditions).
Ice volume correction for stable isotopes. We used Cassidulina neoteretis (δ 18 58 . The diffusive heat flow model was designed in 1D with 2000 cells at a resolution of 1 m with an upper boundary at the seafloor and basal boundary 2000 meters below the seafloor. Initial boundary conditions at 35,000 years were set assuming average present-day BWT of 2°C ( 22 , this study), sediment thermal diffusivity of 3.71 × 10 −7 m 2 s −1 (derived using a thermal conductivity of 1.41 Wm −1 K −1 59 , bulk sediment density of 1900 kgm −3 , and specific heat capacity of 2000 Jkg −1 K −1 3 , and a linear thermal gradient of 0.035°Cm −1 25 . The ice sheet was assumed to be present in the study area from 35,000 to 18,500 years BP 15,22 . During this period, a constant ice sheet thickness of 900 m and a constant basal ice temperature of −2°C 22 were used to determine pressure and temperature conditions within the sediments. After the ice sheet retreat at 18,500 years BP, the reconstructed BWT presented in this study was used as the upper boundary condition until the present day. The transient diffusive heat transport in sediments was then estimated using an explicit finite-difference numerical solution of the Fourier heat equation 60 . The subsurface thermal profile over the past 35,000 years generated by the heat flow model was then integrated with pressure changes resulting from variations in sea level 61 and topographic changes due to glacial isostatic effects. The maximum subsidence generated by the ice loading during the period 35,000 to 18,500 years BP was ca. 85 m 22 . The subsidence and the subsequent uplift were linearly distributed over the modeled period with peak subsidence reaching 18,500 years BP and uplift reaching zero by the present day. The pressure and temperature conditions were then compared with the theoretical hydrate stability phase diagrams generated with the CSMHYD program 58 to estimate the thickness of the gas hydrate stability zone at any time step. The gas hydrate phase boundary was generated for a feed gas composition containing 99.54% methane, 0.41% ethane, and 0.05% propane 22 as well as an assumed pore-water salinity of 35 psu. The sensitivity of GHSZ to the input parameters of the model was analyzed assuming a plausible range of values (see Fig. 4 and Table S1 for details).
Correlation of cores. Sediment core 1059 was taken from a non-methane affected area close to core CAGE 15-2 920GC (core 920) taken from a methane hydrate mound (termed "pingo") 24 . The two cores were correlated using minima and maxima in Zr/Rb and Fe/Ca ratios (Fig. S3). In Yao et al. 24 the age-depth model of core 920 is established by correlation to the reference core CAGE 15-2 921GC, with one radiocarbon date and three tie-points to other cores in the region. In order to allow a direct comparison between core 920 and our core 1059, a new and improved age-depth model of core 920 was established using our much more detailed chronology.
The concentration in µg g −1 of Archaea and bacterial lipids indicating anaerobic oxidation of methane is presented in Yao et al. 24 . However, with the new age model that takes the highly variable sedimentation rates in the area into consideration 15,16 , we calculated the flux of bacterial lipids to indicate the productivity (cm µg g −1 ky −1 ) using the archaeol concentrations 24 and the sedimentation rate. Only flux and δ 13 C of archaeols are shown in Fig. 3 (see Yao et al. 24 for more data).

Data availability
The data is stored at the UiT Open Research Data Repository: https://doi.org/10.18710/ XFYDFL.