Multi-proxy approach to unravel methane emission history of an Arctic cold seep

Arctic Ocean sediments contain large amounts of methane in the form of free gas and gas hydrate. This highly dynamic methane reservoir is susceptible to be modified by bottomwater warming. The warming may lead to gas hydrate destabilization releasing elevated methane fluxes to the seafloor and seawater. Reconstructing past methane dynamics can be achieved by using specific proxies left in the geological record. In this study, we apply a multi-proxy approach for paleo seepage reconstruction from sediment records at gas hydrate mounds (GHMs) in Storfjordrenna (south of the Svalbard archipelago). These shallow water (~380 m water depth) systems are potentially vulnerable to global warming related temperature changes. 14C dating of foraminifera shells indicated an onset of deglaciation in the Storfjordrenna region at ~20 kyr BP and allowed us to establish a stratigraphic context based on sediment Zr/ Rb and Fe/Ca ratios. Several major (between 15 and 17 kyr BP) and minor methane venting phases were identified and interpreted to be related to gas hydrate instability triggered by isostatic adjustment right after the onset of the deglaciation. The detection of all major methane releases was only possible by combining data sets of stable carbon isotope compositions of foraminifera, mineralogy and d13C values of authigenic carbonates, and abundance and stable carbon isotope signatures of lipid biomarkers. The most robust single proxy in this study was provided by the d13C values of archaeal biomarkers. In contrast, the sediment Ba/Ti ratios recorded only the major events. Our results highlight the complexity and heterogeneity of methane dynamics in a small area of some hundred meters across. © 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). 1. General introduction Current global warming raises concerns about greenhouse gas emissions. One of the most potent greenhouse gases is methane. In addition to anthropogenic methane production, methane is also released from natural sources on land (e.g. from wetlands and permafrost) and in the ocean. In the ocean, hot spots of natural methane release are cold seeps, which are ubiquitous in shallow and deep-water settings on convergent and passive continental margins (Judd and Hovland, 2007). At seep sites, methane-rich fluids are expulsed from sediments into the water column at ambient temperature. Methane at cold seeps can be of microbial, thermogenic or abiogenic origin (Whiticar, 1999). Microbial ier Ltd. This is an open access artic methane is (usually) the end product of anaerobic organic matter degradation, thermogenic methane is produced during early diagenesis of buried sedimentary organic matter, and abiotic methane is produced through magmatic and gas-water-rock reactions. The latter formation is very rarely observed and poorly understood (Etiope and Sherwood Lollar, 2013). In surface sediments, methane can occur as free gas, dissolved gas and gas hydrates (GH) (MacLeod, 1982; Bohrmann and Torres, 2013). GH are ice-like crystals consisting of methane molecules that are caged by water molecules; they are stable under high pressure and low temperatures (i.e. conditions that are typically met in Arctic sediments at ~400 m water depth), and where methane availability exceeds its solubility in pore fluids. Bottom water temperatures may rise as a result of climate warming (Westbrook, 2009). Warmer seafloor conditions will ultimately lead to shift of the GH stability limit to greater depth (i.e. higher pressure), and GH that are stable at present will become le under the CC BY license (http://creativecommons.org/licenses/by/4.0/). H. Yao et al. / Quaternary Science Reviews 244 (2020) 106490 2 unstable and dissociate, and methanewill be mobilized (Bohrmann et al.,1998; Suess et al., 2001). Indeed, already a seasonal increase of temperatures in the Arctic Ocean affects sedimentary shallow GH systems (Berndt et al., 2014; Ferr e et al., 2020) and it is expected that climate warming in the near future may accelerate methane release from GH as has been shown for the geological past (Andreassen et al., 2017; Serov et al., 2017). To reconstruct the history of methane emission in the past, we cannot observe the critical variables in a system that no longer exists; but instead, we can use proxies as measurable descriptors for a desired but unobservable environmental variable (Wefer et al., 1999). In this article, we applied a set of the most common (bio) geochemical proxies to investigate benthic methane dynamics in the shallow-water GH area of Storfjordrenna (South of Svalbard), in order to reconstruct the post glacial methane emission history of this Arctic cold seep. We investigate (i) sedimentological parameters; (ii) geochemical imprints of gas/fluid migration; (iii) the oxygen and carbon isotope compositions of carbonate precipitates and foraminifera test (shells), and (iv) past methanotrophic communities and their stable carbon isotope composition. 1.1. Anaerobic oxidation of methane (AOM) In ocean sediment, most methane is consumed under anoxic conditions through the anaerobic oxidation of methane (AOM; Eq. (1)). CH4 þ SO4 / HCO3 þ HS þ H2O Eq. 1 This process is mediated by a microbial consortium of anaerobic methanotrophic archaea (ANME-1, -2, -3) and sulfate-reducing bacteria (SRB, SeepSRB1 cluster and Desulfobulbus sp.) (Boetius et al., 2000; Knittel and Boetius, 2009). Uprising methane and sulfate diffusing from the ocean water are consumed in a distinct sediment horizon, the sulfate-methane transition zone (SMTZ) (Barnes and Goldberg, 1976; Iversen and Jørgensen, 1985), where the highest AOM reaction rates usually take place. 1.2. Sedimentary proxies (barite front and magnetic susceptibility) indicate SMTZ position In sulfate-depleted sediments below the SMTZ, barite is destabilized and dissolved as barium and sulfate ions. As the fluids migrate upwards, the dissolved barium is transported across the SMTZ and once it comes into contact with sulfate above the SMTZ, barite precipitates and builds a ‘barite front’ (Torres et al., 1996, 2003a; Dickens, 2001; Solomon and Kastner, 2012). Indeed, barite is an established proxy for reconstructing the vertical positioning of the SMTZ and can easily be detected by X-ray Fluorescence (XRF) core scanning as Ba/Ti ratios, or Ba counts. Once the flux of methane increases, the SMTZ shifts upward, as does the barite front. The magnitude of the old barite front reduces, but when the Badissolution is slower than the upward flux of methane, a double peak can occur (Kasten et al., 2012). Alternatively, the old barite can disappear if the dissolution is faster than the upward flux of methane. When the methane flux decreases, the SMTZ shifts downward: if there are still enough barium ions in the ambient water, a new barite front forms at the new SMTZ below the old barite front, but if the barium ions are completely consumed, no barite front forms at the new SMTZ. The old barite front is not expected to disappear under decreasing methane flux conditions. Furthermore, the appearance of the barite front also depends on the equilibrium between SMTZ stabilization and barite dissolution/ reprecipitation. The reducing conditions at the SMTZ generate Fe-sulfides (Canfield and Berner, 1987; Peckmann et al., 2001; Riedinger et al., 2006; Dewangan et al., 2013) and cause an alteration of the initial sediment composition and magnetic properties due to the replacement of magnetic Fe-oxides by paramagnetic authigenic Fesulfides, thus reducing the magnetic susceptibility. 1.3. Methane-derived authigenic carbonates as AOM archive During AOM, the production of 13C-depleted bicarbonate increases pore water alkalinity, that results in the precipitation of methane-derived authigenic carbonates (MDACs). Authigenic carbonates are characterized by negative d13C values, often below 30‰ VPDB (Peckmann and Thiel, 2004). At cold seeps, AOM is typically fueled by already 13C-depleted methane, and experiences large isotope fractionation (ε 1⁄4 12e39‰; Holler et al., 2009). However, the d13C values of seep carbonates are isotopically heavier than the parent methane due to admixture of seawater dissolved inorganic carbon (DIC) and AOM-derived DIC during precipitation. Authigenic carbonates represent a good geochemical archive to study the history of methane seeps (Bohrmann et al., 1998; Aloisi et al., 2000; Peckmann et al., 2001; Mazzini et al., 2004; Birgel and Peckmann, 2008; Berndt et al., 2014). The mineralogy of the carbonates also provides some insights into the conditions of the precipitating environment. For example, aragonite is believed to precipitate close to the seafloor because there, the relatively high sulfate concentrations would inhibit the formation of calcite; whereas high magnesium calcite and dolomite would form in deeper sediments where sulfate is typically depleted (Bohrmann et al., 1998; Aloisi et al., 2000, 2002; Han and Aizenberg, 2003). 1.4. Foraminifera shells: templates for authigenic carbonate coating The 13C-depleted bicarbonate produced during AOM can precipitate not only as authigenic carbonate concretions (Aloisi et al., 2002; Blumenberg, 2010) but also as a coating on foraminiferal shells (Panieri et al., 2009). Once dead, both benthic and planktonic species contained in the sediments can record the 13C signature from the AOM process in carbonate layers coating the original tests (Panieri et al., 2016, 2017; Schneider et al., 2017). The coating layers usually exhibit depleted 13C values down to 36‰ with different states of shell preservation (Torres et al., 2003b; Hill et al., 2004; Panieri et al., 2009, 2017; Martin et al., 2010; Schneider et al., 2017). In regular present-day conditions in the Barents Sea, d13C values of the planktonic foraminifera Neogloboquadrina pachyderma ranges between 0.5 and 0.5‰ (Knies and Stein, 1998), while the stable carbon isotope composition of the benthic Cassidulina neoteretis ranges between 0 and 1‰ (Wollenburg et al., 2001). However, d13C values as low as 20‰ can only be explained by diagenetic alteration of foraminifera shells after death. Both C. neoteretis and N. pachyderma are good templates for MDAC precipitation at the SMTZ that cumulatively add a second or third layer of 13C-depleted carbonate around the foraminiferal tests (Schneider et al., 2017; Panieri et al., 2017). Living benthic foraminifera can also record the 13C signature of AOM, probably by using the 13C-depleted bicarbonate to build their tests, utilizing the AOM-microbes as a food source, resulting in slightly negative d13C values to the extent of 5.6‰ (Panieri, 2006; Panieri and Sen Gupta, 2008). Recently, livingMelonis barleeanus in gas-hydrate sediments were found living in close association with putative methanotrophic bacteria (Bernhard and Panieri, 2018). Therefore, foraminifera represent another good geochemical proxy to trace the past methane history in marine sediments. No endemic foraminiferal species have been found at methane seeps so far. Particular species distribution suggests that a high abundance of H. Yao et al. / Quaternary Science Reviews 244 (2020) 106490 3 endobenthic foraminifera may indicate methane seeps (Panieri and Sen Gupta, 2008) and agglutinated foraminifera are less abundant in sediments influenced by methane seepage, suggesting that this group of foraminifera does not tolerate the geochemical conditions at seeps (Panieri and Sen Gupta, 2008; Dessandier et al., 2019). 1.5. Lipid biomarkers of AOM communities Lipid biomarkers are constituents of the cellular membranes originating, for the case of AOM, mostly from ANMEs, and SRB (Niemann and Elvert, 2008, and references therein). Major components of archaeal cell membranes are isoprenoidal glycerol dior tetraether lipids, which are distinct from the acyl monoor di-ester lipids found in bacteria. The different sulfate-reducing partner bacteria in the AOM consortia in modern and paleoenvironments have been identified mostly through the analysis of fatty acids fingerprints as well as monoalcyl glycerol ethers (Hinnrichs et al., 1999; Niemann and Elvert, 2008; Lee et al., 2018; Yao et al., 2019). Furthermore, methane of microbial and thermogenic origin, the most common methane sources, is usually13C-depleted (typically 50 to 100‰; Whiticar, 1999). AOM further discriminates against the heavy CeCH4 (Holler et al., 2009) and the isotope enrichment factor ε during AOM (methane lipids) is ~39‰ Fig. 1. Map of Svalbard margin showing the bathymetry of the study area and the locatio Storfjordrenna. (Niemann and Elvert, 2008) Consequently, AOM associated lipids often display extremely depleted 13C signatures, which allows to distinct them from background sources (e.g., organic matter (OM) of photosynthetic ormethanogenic origin). The partner SRB in AOM also displays 13C-depleted lipid signatures, most probably because they incorporate AOM-derived 13C-depleted CO2 (Wegener et al., 2008). The lipid biomarkers of ANME and SRB can identify past AOM communities, making lipids a valuable tool for understanding methane cycling in anaerobic sedimentary environments (Hinrichs et al., 1999; Niemann et al., 2006a,b; Lee et al., 2018; Yao et al., 2019). However, depending on the environmental setting, lipids are also subjected to degradation (microbial and diagenetic), leading, for example, to shifts in the saturation patterns of alkyl chains (Volkman and Smittenberg, 2017) or to an overall reduction in concentration over time. Also, the diagnostic d13C signature of AOM lipids may be overprinted by lipids related to processes other than AOM, for example, methanogenesis (Lim et al., 2012; Kaneko et al., 2013). 1.6. Living biofilm associated with AOM The most promising proxy for AOM in sediments would be n of the investigated cores. The black square indicates the location of the study area, H. Yao et al. / Quaternary Science Reviews 244 (2020) 106490 4 living AOM-biomass, because the microbial community can be identified by genetic markers such as DNA. Such biofilms were described in subseafloor fractures in methane seeps from Cascadia margin and Indian Ocean margin (Briggs et al., 2011), and recently in the Storfjordrenna GHM5 (Gründger et al., 2019). However, finding of biofilms in pockets, cracks or fractures within the sediment matrix is by-chance and thus very rare.


General introduction
Current global warming raises concerns about greenhouse gas emissions. One of the most potent greenhouse gases is methane. In addition to anthropogenic methane production, methane is also released from natural sources on land (e.g. from wetlands and permafrost) and in the ocean. In the ocean, hot spots of natural methane release are cold seeps, which are ubiquitous in shallow and deep-water settings on convergent and passive continental margins (Judd and Hovland, 2007). At seep sites, methane-rich fluids are expulsed from sediments into the water column at ambient temperature. Methane at cold seeps can be of microbial, thermogenic or abiogenic origin (Whiticar, 1999). Microbial methane is (usually) the end product of anaerobic organic matter degradation, thermogenic methane is produced during early diagenesis of buried sedimentary organic matter, and abiotic methane is produced through magmatic and gas-water-rock reactions. The latter formation is very rarely observed and poorly understood (Etiope and Sherwood Lollar, 2013). In surface sediments, methane can occur as free gas, dissolved gas and gas hydrates (GH) (MacLeod, 1982;Bohrmann and Torres, 2013). GH are ice-like crystals consisting of methane molecules that are caged by water molecules; they are stable under high pressure and low temperatures (i.e. conditions that are typically met in Arctic sediments at~400 m water depth), and where methane availability exceeds its solubility in pore fluids.
Bottom water temperatures may rise as a result of climate warming (Westbrook, 2009). Warmer seafloor conditions will ultimately lead to shift of the GH stability limit to greater depth (i.e. higher pressure), and GH that are stable at present will become unstable and dissociate, and methane will be mobilized (Bohrmann et al., 1998;Suess et al., 2001). Indeed, already a seasonal increase of temperatures in the Arctic Ocean affects sedimentary shallow GH systems (Berndt et al., 2014;Ferr e et al., 2020) and it is expected that climate warming in the near future may accelerate methane release from GH as has been shown for the geological past Serov et al., 2017).
To reconstruct the history of methane emission in the past, we cannot observe the critical variables in a system that no longer exists; but instead, we can use proxies as measurable descriptors for a desired but unobservable environmental variable (Wefer et al., 1999). In this article, we applied a set of the most common (bio) geochemical proxies to investigate benthic methane dynamics in the shallow-water GH area of Storfjordrenna (South of Svalbard), in order to reconstruct the post glacial methane emission history of this Arctic cold seep. We investigate (i) sedimentological parameters; (ii) geochemical imprints of gas/fluid migration; (iii) the oxygen and carbon isotope compositions of carbonate precipitates and foraminifera test (shells), and (iv) past methanotrophic communities and their stable carbon isotope composition.

Anaerobic oxidation of methane (AOM)
In ocean sediment, most methane is consumed under anoxic conditions through the anaerobic oxidation of methane (AOM; Eq. (1)).

Sedimentary proxies (barite front and magnetic susceptibility) indicate SMTZ position
In sulfate-depleted sediments below the SMTZ, barite is destabilized and dissolved as barium and sulfate ions. As the fluids migrate upwards, the dissolved barium is transported across the SMTZ and once it comes into contact with sulfate above the SMTZ, barite precipitates and builds a 'barite front' (Torres et al., 1996(Torres et al., , 2003aDickens, 2001;Solomon and Kastner, 2012). Indeed, barite is an established proxy for reconstructing the vertical positioning of the SMTZ and can easily be detected by X-ray Fluorescence (XRF) core scanning as Ba/Ti ratios, or Ba counts. Once the flux of methane increases, the SMTZ shifts upward, as does the barite front. The magnitude of the old barite front reduces, but when the Badissolution is slower than the upward flux of methane, a double peak can occur (Kasten et al., 2012). Alternatively, the old barite can disappear if the dissolution is faster than the upward flux of methane. When the methane flux decreases, the SMTZ shifts downward: if there are still enough barium ions in the ambient water, a new barite front forms at the new SMTZ below the old barite front, but if the barium ions are completely consumed, no barite front forms at the new SMTZ. The old barite front is not expected to disappear under decreasing methane flux conditions. Furthermore, the appearance of the barite front also depends on the equilibrium between SMTZ stabilization and barite dissolution/ reprecipitation.
The reducing conditions at the SMTZ generate Fe-sulfides (Canfield and Berner, 1987;Peckmann et al., 2001;Riedinger et al., 2006;Dewangan et al., 2013) and cause an alteration of the initial sediment composition and magnetic properties due to the replacement of magnetic Fe-oxides by paramagnetic authigenic Fesulfides, thus reducing the magnetic susceptibility.

Methane-derived authigenic carbonates as AOM archive
During AOM, the production of 13 C-depleted bicarbonate increases pore water alkalinity, that results in the precipitation of methane-derived authigenic carbonates (MDACs). Authigenic carbonates are characterized by negative d 13 C values, often below À30‰ VPDB (Peckmann and Thiel, 2004). At cold seeps, AOM is typically fueled by already 13 C-depleted methane, and experiences large isotope fractionation (ε ¼ 12e39‰; Holler et al., 2009). However, the d 13 C values of seep carbonates are isotopically heavier than the parent methane due to admixture of seawater dissolved inorganic carbon (DIC) and AOM-derived DIC during precipitation. Authigenic carbonates represent a good geochemical archive to study the history of methane seeps (Bohrmann et al., 1998;Aloisi et al., 2000;Peckmann et al., 2001;Mazzini et al., 2004;Birgel and Peckmann, 2008;Berndt et al., 2014). The mineralogy of the carbonates also provides some insights into the conditions of the precipitating environment. For example, aragonite is believed to precipitate close to the seafloor because there, the relatively high sulfate concentrations would inhibit the formation of calcite; whereas high magnesium calcite and dolomite would form in deeper sediments where sulfate is typically depleted (Bohrmann et al., 1998;Aloisi et al., 2000Aloisi et al., , 2002Han and Aizenberg, 2003).

Foraminifera shells: templates for authigenic carbonate coating
The 13 C-depleted bicarbonate produced during AOM can precipitate not only as authigenic carbonate concretions (Aloisi et al., 2002;Blumenberg, 2010) but also as a coating on foraminiferal shells (Panieri et al., 2009). Once dead, both benthic and planktonic species contained in the sediments can record the 13 C signature from the AOM process in carbonate layers coating the original tests (Panieri et al., 2016Schneider et al., 2017). The coating layers usually exhibit depleted 13 C values down to À36‰ with different states of shell preservation (Torres et al., 2003b;Hill et al., 2004;Panieri et al., 2009Panieri et al., , 2017Martin et al., 2010;Schneider et al., 2017).
In regular present-day conditions in the Barents Sea, d 13 C values of the planktonic foraminifera Neogloboquadrina pachyderma ranges between À0.5 and 0.5‰ (Knies and Stein, 1998), while the stable carbon isotope composition of the benthic Cassidulina neoteretis ranges between 0 and 1‰ (Wollenburg et al., 2001). However, d 13 C values as low as À20‰ can only be explained by diagenetic alteration of foraminifera shells after death. Both C. neoteretis and N. pachyderma are good templates for MDAC precipitation at the SMTZ that cumulatively add a second or third layer of 13 C-depleted carbonate around the foraminiferal tests Panieri et al., 2017).
Living benthic foraminifera can also record the 13 C signature of AOM, probably by using the 13 C-depleted bicarbonate to build their tests, utilizing the AOM-microbes as a food source, resulting in slightly negative d 13 C values to the extent of À5.6‰ (Panieri, 2006;Panieri and Sen Gupta, 2008). Recently, living Melonis barleeanus in gas-hydrate sediments were found living in close association with putative methanotrophic bacteria (Bernhard and Panieri, 2018). Therefore, foraminifera represent another good geochemical proxy to trace the past methane history in marine sediments. No endemic foraminiferal species have been found at methane seeps so far. Particular species distribution suggests that a high abundance of endobenthic foraminifera may indicate methane seeps (Panieri and Sen Gupta, 2008) and agglutinated foraminifera are less abundant in sediments influenced by methane seepage, suggesting that this group of foraminifera does not tolerate the geochemical conditions at seeps (Panieri and Sen Gupta, 2008;Dessandier et al., 2019).

Lipid biomarkers of AOM communities
Lipid biomarkers are constituents of the cellular membranes originating, for the case of AOM, mostly from ANMEs, and SRB , and references therein). Major components of archaeal cell membranes are isoprenoidal glycerol di-or tetraether lipids, which are distinct from the acyl mono-or di-ester lipids found in bacteria. The different sulfate-reducing partner bacteria in the AOM consortia in modern and paleoenvironments have been identified mostly through the analysis of fatty acids fingerprints as well as monoalcyl glycerol ethers (Hinnrichs et al., 1999;Niemann and Elvert, 2008;Lee et al., 2018;Yao et al., 2019). Furthermore, methane of microbial and thermogenic origin, the most common methane sources, is usually 13 C-depleted (typically À50 to À100‰; Whiticar, 1999). AOM further discriminates against the heavy 13 CeCH 4 (Holler et al., 2009) and the isotope enrichment factor ε during AOM (methane -lipids) is~39‰  Consequently, AOM associated lipids often display extremely depleted 13 C signatures, which allows to distinct them from background sources (e.g., organic matter (OM) of photosynthetic or methanogenic origin). The partner SRB in AOM also displays 13 C-depleted lipid signatures, most probably because they incorporate AOM-derived 13 C-depleted CO 2 (Wegener et al., 2008).
The lipid biomarkers of ANME and SRB can identify past AOM communities, making lipids a valuable tool for understanding methane cycling in anaerobic sedimentary environments (Hinrichs et al., 1999;Niemann et al., 2006a,b;Lee et al., 2018;Yao et al., 2019). However, depending on the environmental setting, lipids are also subjected to degradation (microbial and diagenetic), leading, for example, to shifts in the saturation patterns of alkyl chains (Volkman and Smittenberg, 2017) or to an overall reduction in concentration over time. Also, the diagnostic d 13 C signature of AOM lipids may be overprinted by lipids related to processes other than AOM, for example, methanogenesis (Lim et al., 2012;Kaneko et al., 2013).

Living biofilm associated with AOM
The most promising proxy for AOM in sediments would be living AOM-biomass, because the microbial community can be identified by genetic markers such as DNA. Such biofilms were described in subseafloor fractures in methane seeps from Cascadia margin and Indian Ocean margin (Briggs et al., 2011), and recently in the Storfjordrenna GHM5 . However, finding of biofilms in pockets, cracks or fractures within the sediment matrix is by-chance and thus very rare.

Study area
Storfjordrenna is located between the islands of Spitsbergen to the north and Bjørnøya to the south. The Storfjorden trough is a broad structure that extends from the mouth of the fjord to the shelf edge. The trough was generated by glacial erosion of former ice streams of the Svalbard-Barents Sea ice sheet (Vorren et al., 1998;Dowdeswell and Elverhoi, 2002).
During a research cruise in May 2015 with R/V Helmer Hansen (CAGE 15e2 cruise), a group of mounds actively seeping methane gas from the seafloor to the water column was discovered in the Storfjordrenna area (~380 m water depth, Serov et al., 2017). These mounds are about 500 m in diameter and about 10 m in height above the seafloor. GHs were recovered from these mounds; therefore, they are referred to as gas hydrate pingos Waage et al., 2019) or gas hydrate-bearing mounds (GHM) (Hong et al., 2017. Real-time visually guided observations with the TowCam-Multicore system (TC-MC) (CAGE 15e2 cruise) and ROV dives (CAGE 16e5 cruise) revealed methane bubble streams rising from the mounds GHM 3 (west) and other discovered mounds, and only GHM 5 (east) (Fig. 1) showed no sign of bubble streams. GHM 3 and GHM 5 are the objectives of this study. GHs were recovered at depth as shallow as 40 cm below seafloor (bsf) from GHM 3 during various cruises. Because of our observations of bubble streams, we classify GHM 3 as active, whereas GHM 5 was acoustically inactive. No flares in the water column could be found during various surveys (in several years). We thus consider GHM 5 as inactive, and the motivation for this definition is explained in more detail in the discussion.
At Storfjordrenna, repeated growth and retreat of grounding glaciers shaped the seafloor bathymetry (Patton et al., 2015). Glacial-isostatic events may have also reactivated fractures and faults, imposing a temporal variability on spatially heterogeneous fluid flows (Andreassen et al., 2007;Wallmann et al., 2018)). In both the active (GHM 3) and inactive (GHM 5) mound, Waage et al. (2019) observed chimneys connected to fractures and faults within the underlying sedimentary rocks, suggesting that faultcontrolled Paleocene hydrocarbon pool was probably responsible for charging the GHMs with free gas.

Core collection and non-destructive analyses (MSCL and XRF)
We selected six sediment gravity cores collected during three CAGE cruises 15e2, 15e6, and 16e5 (Table 1) for the multi-proxy approach. Three cores (CAGE 15e2 940GC, CAGE 15e6 1520GC, and 1522GC) were retrieved from the most acoustically active gas hydrate mound (GHM3), and two cores (CAGE 15e2 920GC, CAGE 16e5 1070GC) were recovered from mound GHM 5 (here we did not record flares with our multibeam and thus considered it as acoustically inactive). Finally, one reference core (CAGE-15-2 921GC) was retrieved from an area with no indication of seep activity or evidence of gas hydrate occurrence.
Upon recovery, the cores were cut into 1 m sections, split longitudinally into working and archive halves. The working half was subsampled onboard for stable isotope analyses of foraminifera and authigenic carbonates, lipid biomarkers, and radiocarbon dating, and the archive half was stored at 4 C for the nondestructive analyses of magnetic susceptibility and elementgeochemical data.
Magnetic susceptibility was measured at 1 cm intervals using GEOTEK Multi Sensor Core Logger (MSCL-S) system. The X-radiography was obtained using a GEOTEK X-ray core imaging system (MSCL-XCT 3.0) using an X-ray intensity of 120 kV (Core 1520GC was not measured due to abundant carbonate nodules causing low surface contact with the sensor, therefore impeding the measurements). Element-geochemical data were acquired at 1 cm interval with an Avaatech XRF Core Scanner. Zr and Rb were quantified with 30 kV, 2000 mA, at 10 s using the Pd filter. The raw data were processed with the software WinAxil. For this study, we show the zirconium (Zr), rubidium (Rb), calcium (Ca), titanium (Ti), iron (Fe) and barium (Ba) counts ratios.

Oxygen and carbon isotope analyses of foraminifera shells and authigenic carbonates
Sediment samples were collected with variable spacing (5 cm, 10 cm and 20 cm) throughout the cores for carbon and oxygen stable isotope analyses of foraminifera shells. Samples were freezedried and wet sieved through 63 and 125 mm mesh size sieves, and the sieved residue was oven-dried at 40 C. The dry residue was examined under a light microscope to pick out planktonic and benthic foraminifera for stable isotope analyses. Planktic foraminifera Neogloboquadrina pachyderma (Ehrenberg, 1861) and benthic foraminifera Cassidulina neoteretis (Seidenkrantz, 1995), and Melonis barleeanum (Williamson, 1858) were picked from the >125 mm size fraction and examined using light microscopy to determine the state of preservation of foraminifera carbonate tests. Special attention was paid to evidence of dissolution, diagenesis, or postdepositional alteration due to the presence of authigenic carbonates (secondary overgrowth) on the foraminifera tests to be analyzed for stable isotope analyses (d 13 C and d 18 O). Approximately 10e30 individuals were selected for each analyses.
The carbon and oxygen stable isotope compositions of the N. pachyderma from core 921GC were measured on a Thermo Finnigan MAT252 Isotope Ratio Mass Spectrometer (IRMS) coupled to a CarboKiel-II carbonate preparation device at the stable isotope laboratory at Oregon State University. For the other cores, d 13 C-and d 18 O-measurements of foraminifera were conducted at the stable isotope laboratory at The Arctic University of Norway in Tromsø UiT using a Thermo Scientific MAT253 IRMS coupled to a Gasbench II system. Foraminiferal shells were placed in 4.5 mL vials and flushed with He gas. Five drops of water-free H 3 PO 4 were added manually. After equilibration (>3 h at 50 C), the samples were analyzed. Normalization to the VPDB for carbon and oxygen isotopes was done using in-house standards (1.96‰, À10.21‰, and À48.95‰ for d 13 C, and À2.15‰ and À18.59‰ for d 18 O). Analytical precision was estimated to be better than 0.07‰ for d 13 C and 0.08‰ for d 18 O against the certified standard NBS-19 (Dessandier et al., 2019). Some core intervals (Figs. 4h, 8g and 9d) contain carbonate nodules (0.5e2 cm in diameter) that were selected for elemental composition and stable oxygen and carbon analyses. The d 13 C and d 18 O were measured using the same method as the foraminiferal shell described above at UiT.

Radiocarbon dating
Radiocarbon dating was carried out at the Beta Analytic Radiocarbon Dating facilities in Miami, US. Conventional radiocarbon ages were obtained from N. pachyderma sampled from 921GC (depth interval 325e326 cm; laboratory code Beta-492,581) and bivalve shell from 1070GC (depth interval 316e317 cm; laboratory code Beta-492,580). A regional reservoir age correction DR of 67 ±34 was applied (Mangerud et al., 2006) for the calibration of marine carbonate samples, Marine 13 is the database used for the calibration of carbon 14 dating results. The age model is based on the calibrated ages obtained from the peaks of the probability curves within the 2s range.

Lipid biomarker analyses
Samples for biomarker analyses were collected immediately after the recovery of the sediment cores on-board with a methanolcleaned spatula at a spatial resolution of 2e10 cm, wrapped in aluminium foil, and subsequently stored frozen at À20 C until further analyses.
Lipid biomarkers were extracted and analyzed according to procedures previously described in Niemann et al. (2005) and Blees et al. (2014). Briefly, a total lipid extract (TLE) was obtained by ultrasonication of~20 g wet sediment in four extraction steps with solvent mixtures with decreasing polarity: dichloromethane (DCM)/methanol (MeOH) 1: 2; DCM/MeOH 2: 1; and DCM for the last two extraction steps. The TLE was then saponified with KOH and a neutral lipid fraction was extracted. The remaining polar fraction, which contained fatty acids was methylated yielding fatty acid methyl esters (FAME) for gas chromatographic (GC) analyses. Double bond positions of FAMES were determined by analyzing dimethyl-disulfide (DMDS) adducts. The neutral fraction was further separated into hydrocarbons, ketones, and alcohols, of which the latter was further derivatized to form trimethylsilyl (TMS) adducts for gas chromatographic analyses.
Individual lipid molecules were analyzed with a Gas Chromatograph (GC; Thermo Scientific TRACE™ Ultra), equipped with a capillary column (Rxi-5ms, 50 m, 0.2 mm id, 0.33 mm df) with He as a carrier gas at a constant flow of 1 mL min À1 . The initial oven temperature was 50 C, held for 2 min and then increased to 140 C at a rate of 10 C/min, held at 140 C for 1 min, then further increased to 300 C at 4 C/min. Final hold time was 63 min to analyze FAMEs or 160 min to analyze larger/high boiling point lipids (in the hydrocarbon and alcohol fractions), respectively. Concentrations were determined by flame-ionization detection (FID) against internal standards. Unknown compounds were identified with a quadrupole mass spectrometry unit at the chromatography periphery (Thermo Scientific DSQ II). Similarly, compound-specific stable carbon isotope ratios were determined using a IRMS unit (Thermo Scientific Delta V Advantage) coupled to a gas chromatography setup as outlined above. d 13 C values reported here have an analytical error of ±1‰.

Carbonate mineralogy
Mineralogy of carbonate nodules found in cores 1520GC, 920GC, and 940GC was determined by X-ray diffraction (XRD) of homogenized powders, prepared from samples grounded in isopropanol cleaned McCrone mill. Unoriented specimens of the dried powders were prepared by side loading. XRD were realized on a Bruker D8 Advance X-ray diffractometer (Cu Ka radiation in 3e75 2q range) (Sauer et al., 2017). Quantitative data were obtained with the Rietveld algorithm-based code, TOPAS 5, provided by Bruker. Following a displacement correction of the spectrum made on the main quartz peak, the displacement of calcite d 104 was used to estimate the amount of MgCO 3 mol% (Goldsmith et al., 1958).

Chronology of the cores
To establish the chronology of the cores of the studied site, we use the 365 cm long reference core 921GC as a template. This core was collected in-between the gas hydrate mounds (GHMs) where seismic and echo sounder surveys during cruises CAGE 15e2, 15e6 and 16-5 did not find indications of seepage. We used magnetic susceptibility, Zr/Rb and Fe/Ca ratio and foraminiferal d 18 O from core 921GC to correlate our sedimentary record with literature data (Lucchi et al., 2013;Jessen et al., 2010;Knies et al., 2018). After that, the Zr/Rb and Fe/Ca ratios of the other investigated cores were matched against core 921GC. Table 2 shows the dating results from different cores.
The Zr/Rb ratio is generally used to indicate the particle size, and it increases as sediment grain-size increases (Dypvik and Harris, 2001). Here we use it to identify ice rafted debris (IRD) content in the sediment cores (Lucchi et al., 2013). Importantly, the Zr/Rb ratio is not affected by methane-derived diagenesis (Dypvik and Harris, 2001). The Fe/Ca ratio can differentiate terrestrial source input vs. marine origin as higher Fe/Ca indicate terrestrial source with higher organic input (Rothwell and Croudace, 2015). For our study, all cores were recovered from a relatively small area. The primary organic matter source is thus expected to be similar, so that we use Fe/Ca ratios similar to Zr/Rb ratios for establishing stratigraphic correlation.
The magnetic susceptibility values in core 921GC agreed well  (Lucchi et al., 2013: Jessen et al., 2010. We used the stratigraphic sequence from these known cores, and the dating results from 921GC to establish three stratigraphic marker horizons. The three stratigraphic marker horizons for the Storfjordrenna area cover a period from post-LGM to Early Holocene. The stratigraphic marker horizon (i) is highlighted in grey at 325 cm bsf in Figs. 2 and 3, and coincides with a foraminifera rich event. One direct dating point of the planktonic foraminifera N. pachyderma at its base (325e326 cm) revealed an age of 20,609 ± 205 cal BP (Fig. 3). The marker horizon (ii) corresponds to Heinrich event H1, it is a foraminifera-rich and reddish oxidized event at 220e240 cm, characterized by increasing magnetic susceptibility, and Zr/Rb and Fe/Ca ratios (highlighted with a green shaded bar in Figs. 2 and 3). The marker horizon (iii) corresponds to the Younger Dryas to early Holocene transition and is characterized by a decreasing Zr/Rb ratio together with increasing Fe/Ca ratio (highlighted with a yellow shaded bar in Figs. 2 and 3). Lucchi et al. (2013), suggested that the onset of deglaciation on the Storfjorden TMF occurred around 20,000 cal. BP and no later than 18,000 cal. BP based on the observation of an oxidized event (assigned as OX-2), located at the base of the coarse-massive-IRD Fig. 2. Reference core 921GC, a. X-ray image, b. magnetic susceptibility, c. XRF sediment ratios of Zr/Rb, d. Fe/Ca, e. stable carbon and f. oxygen isotope records from N. pachyderma. g. In green, partial stable oxygen isotope from core 921GC; in black, stable oxygen isotope records from the reference core obtained by Knies et al. (2018). Red arrows show radiocarbon dating from N. pachyderma, and dash lines show correlation from Jessen et al. (2010) in the area based on magnetic susceptibility. Yellow, green and grey shaded bars are stratigraphic horizon markers. All isotope data are available in the supplementary file. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) facies. The oxidized event (OX-2) is believed to be derived from near seabed oxidation of detrital Fe-minerals under interglacial well-ventilated bottom water conditions marking the inception of deglaciation with the release of fresh oxygenated waters from ice melt mixing with bottom ocean waters in the Storfjorden TMF. We found the same oxidized event in 921GC below the marker horizon (i) at 330e340 cm bsf (Fig. 2). Together with the signals of a foraminifera bloom and high Zr/Rb-ratios, this observation supports the earlier interpretation of the deglaciation onset in the study area at~20,000 cal. BP, or even earlier. A second bloom of planktonic foraminifera follows the starting deglaciation in the area, suggesting a short, renewed surface productivity under warmer, oxygen/ nutrient-rich conditions (Lucchi et al., 2013). We observed decreasing Zr/Rb together with increasing Fe/Ca ratios in the XRF profiles of the reference core 921GC (Fig. 3). The changes in Zr/Rb ratios are independent of methane dynamics but rather associated to grain size (Dypvik and Harrison, 2001). Fe/Ca ratios are influenced by the input of terrestrial organic matter (Rothwell and Croudace, 2015), which is not to be expected to vary on the spatial scale of the area investigated here. Hence, we use these two ratios to wiggle-match the Zr/Rb and Fe/Ca ratios of the cores against core 912GC for a more robust chronostratigraphic correlation of the other cores. We did not use magnetic susceptibility (Mag. Sus.) for correlation purposes because iron oxides react with H 2 S (which is abundant at seeps because of AOM, see Eq. (1)) forming paramagnetic pyrite (FeS 2 ) (Canfield and Berner, 1987). The pyrite can cause a reduction in the Mag. Sus. signal in cores from sites with past/present AOM. As Mag. Sus. is overprinted by AOM, we used Mag. Sus.to infer the vertical migration of the SMTZ (note that the depth of the SMTZ may be independent from the stratigraphic chronosequence). We also exclude the use of d 18 O of N. pachyderma for correlation purpose because in our cores we recovered GHs and it has been shown that oxygen isotopes of foraminiferal calcareous shells from GH settings may be affected by GH dissociation (Dessandier et al., 2019).
In the five cores recovered from active and inactive GHMs in Storfjordrenna, at least one of the stratigraphic marker horizons can be recognized in all sites allowing a correlation among the cores (Fig. 3). In addition to Core 912GC, the first foraminifera rich event horizon marker (i) was also found in core 1522GC from GHM3 at the bottom of the core (305e310 cm). The reddish oxidized event (OX-1) horizon marker (ii) was recognized in core 1522 GC at 220 cm (Fig. 7), 200 cm in 920GC (Fig. 4) and 300 cm in 940GC ( Fig. 9) (marker indicated in green color). The radiocarbon dating of 16,035 ± 90 cal BP obtained from N. pachyderma in core 1522GC (Table 2) confirms the upper boundary of the oxidized event. The radiocarbon dating of bivalve shell (Table 2) at the bottom of core 1070GC corresponds to the beginning of Heinrich event 1.
The correlation of all cores indicates that the sites might have experienced variable erosion; for instance, cores 940GC and 1070GC retained more early Holocene sediment than the other cores. The proximity of the cores excludes the possibility of erosion caused by a strong bottom water current that would have affected the entire area. Because in Storfjordrenna TMF early Holocene sediments appeared to be well-preserved (Lucchi et al., 2013), we suggested that the local erosion observed in our cores might be due to trawling activity in the area, as we could also see trawl marks on the seafloor during our video camera and ROV surveys. Nonetheless, the relatively similar depositional characteristics in the studied area exclude the influence of major sedimentation events, such as mass transport deposits (MTDs).

Multi-proxy approach applied to inactive GHM 5
Previous investigations on gas hydrate mound GHM5 revealed a low methane flux Sen et al., 2018), confirmed by seismic surveys showing no prominent acoustic blanking or gas chimneys beneath the mound (Waage et al., 2019) and absence of gas flares in the water column . Waage et al. (2019) suggested that GHM5 was more active in the past when the gas reservoir beneath the mound was probably better connected to Paleocene sedimentary rock as shown from interpretation of seismic profiles. Porewater geochemistry indicates that today, the GHM5 is in a post active stage where the gas reservoir is exhausted, and only dissolved methane is delivered to the subsurface . The absence of chemosynthetic megafauna (e.g., frenulate aggregations) at the seafloor on top of the mound suggests a potential low H 2 S flux, which is indicative of a low AOM activity and thus low methane flux within the sediment (Sen et al., 2018). Therefore, we classify this mound as 'inactive' relative to other GHMs in the area.

Core 920GC
Applying the multi-proxy approach to core 920GC, we were able to identify four major seeping events (A to D) during the last 18 kyr (Fig. 4), and to estimate duration and relative seeping intensity of these events. Event A refers to the most recent event with a proxy record in shallow sediments, while B to D are successively older events with proxy records stored at increasing sediment depth.
Sulfate and methane profiles in core 920GC show that the present-day SMTZ is located at~150 cm (event A, orange, Fig. 4). At this depth, there is a small Ba/Ti peak, M. barleeanum tests show a low d 13 C value (À5.1‰), and there are increasing concentrations of archaeal and bacterial lipid biomarkers associated with decreasing d 13 C values of these biomarkers. The straight-line shape of the sulfate profile indicated that the pore water geochemistry reached steady-state conditions, and the small Ba/Ti peak implies a recent establishment of the steady-state conditions.
Low d 13 C values of foraminiferal tests occur also in the other three intervals (B, C, and D). In events B and C, respectively the benthic foraminifera carbonate shells (-10.72‰ and À13.7‰ for C. neoteretis, À8.63‰ and À8.5‰ for M. barleeanum) have lower values than the planktonic foraminifera (À5.49‰ and À10.85‰ for N. pachyderma). These negative values are due to secondary authigenic carbonate overgrowth on foraminiferal shells (Panieri et al., 2016;Schneider et al., 2018). The 13 C signature of foraminifera is a sensitive and stable proxy to infer a past seepage activity; however, due to the secondary overgrowth, it is difficult to determine the exact precipitation timing of the diagenetic alternation. No Ba/Ti peak was observed in event B. A possible scenario for the absence of a barite front could be that event B was a rapid and short-lasting methane seepage episode. It may not have been long enough to influence the barite dissolution and precipitation equilibrium substantially because the time needed to produce a barium peak at the SMTZ should be in the order of thousands of years (Riedinger et al., 2006). Event C is characterized by two Ba/Ti peaks of different intensities. The upper larger Ba/Ti peak suggests a period of enhanced indicates that the precipitation of MDACs occurred within deeper sediment, because high sulfate concentration inhibit the formation of Mg-calcite (Burton, 1993;Bohrmann et al., 1998). The exclusive presence of MDACs in event C, and their absence in other events, indicate pervasive and high methane flux during event C. In addition to the depleted foraminiferal d 13 C values, event C was characterized by high abundance of the lipid biomarkers archaeol, which also displayed low d 13 C values (down to À100‰) in this horizon.
The bacterial FA C16:1u5c shows the highest concentration in event C in core 920GC (Fig. 4g), and the corresponding d 13 C value is moderately depleted to À50‰ (Fig. 4h). This phenomenon might be due to the overprint from other non-AOM processes, such as organic matter degradation (Elvert et al., 2003). Yet both archaeal and bacterial biomarkers indicate the presence of AOM communities in this event in the past. Event D is located just below the current SMTZ and it is characterized by highest archaeol and relatively high FA C16:1u5c abundances coupled to low d 13 C-values. Similar to event B, a barium front is also absent in event D.

Core 1070GC
Core 1070GC (Fig. 5), similar to 920 GC, was recovered from the inactive GHM 5 and showed only two major seeping events. The current SMTZ (event A at ca 60 cm bsf) was accompanied by a large Ba/Ti peak above it. We also observed a biofilm at 68 cm bsf in a sediment cavity,~10 cm below the current SMTZ, which is composed of the AOM-related microbes . This indicates a steady methane flux that was able to sustain high biomass of AOM microbes, but it was too low to escape from the sediments. This is confirmed by the absence of active gas bubbling at the inactive GHM5 and the bright spot found in seismic profiles, which indicates that gas supply migrating from deeper sediments is trapped in the subsurface (Waage et al., 2019).
The d 13 C values measured on foraminifera at the current SMTZ are not very low (À0.35‰ for N. pachyderma and À1.18‰ for C. neoteretis, Fig. 5d). The pore water sulfate and methane profiles imply steady-state conditions and the formation of an AOM-biofilm may need decades to hundreds of years to form with a steady methane flux because of the extremely slow growth of the AOM communities (Nauhaus et al., 2007;Puglini et al., 2019). Hence, the steady-state was stable for an extended period of at least decades. The formation of secondary overgrowth of MDAC on the foraminiferal tests may also take decades to hundreds of years, depending on the different test-types. Smooth and imperforate tests of C. neoteretis are less likely to accommodate crystalline overgrowth than the N. pachyderma tests with a higher porosity and surface area, as already observed by Cook et al. (2011) and Consolaro et al. (2015). The presence of the AOM-biofilm 10 cm below the current SMTZ at the 'inactive' GHM5 implies that this kind of biofilm does not necessarily occur at high methane flux sites, but rather at places where methane flux is low but persistent , as here at the GHM 5.
The event B is characterized by negative d 13 C excursions of foraminifera in the interval 185e280 cm. The d 13 C values of foraminifera are À2.8‰ (planktonic) and À7.15‰ (benthic) on average. These values would agree with both, a diffuse steady-state or a long-lived focused fluid upflow. Two events were identified in core 1070GC, the current event A, and an older event B where the past SMTZ located at 185 cm, indicates that an increase in the methane flux shifted the SMTZ upward.
All proxy records from the two cores from the inactive GHM 5 show a complex methane seepage history, characterized by significant variability of methane discharge in space and time. Fig. 6. An illustrative summary of core 920GC methane emission history. Core 920GC was used for the illustration as this core had a complete profile of all the proxies used in this paper. Current event A marked in orange, other events in light brown. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) Site 920GC experienced a high and intense flux in the past, and then lower and fluctuating fluxes (Fig. 6) as a system dominated by diffusion. Whereas site 1070GC became more active than 920GC, as indicated by an upward shift of the SMTZ. Assessing the timing of seepage events is challenging, the host sediment age at the past SMTZ depth does not necessarily correspond to the age of one seeping event. The age of host sediments in the event D in core 920GC and event B in core 1070GC ranges from 17 to 12 Ka, which is consistent with the suggested high fluid activity in the Barents Sea (Cr emi ere et al., 2016), that followed the local deglaciation that started at 18 Ka based on the GHSZ and glacial isostatic modeling .
From our interpretation, the methane flux increased after event D at site 920GC, and abundant carbonate crusts were formed. We hypothesize that the oldest events (event D in 920GC and event B in 1070GC) in both cores could be associated with high fluid activity triggered by the Barents Sea deglaciation and ice stream retreat. The deglaciation, and the following isostatic adjustment, caused changes in the GH stability zone and therefore the methane expulsion occurred . The methane flux increased and decreased as the ice stream retreated and re-advanced, specifically at the GHM 5 area.

Multi-proxy approach applied to active mound GHM 3
GHM 3 is an active methane seep site as shown by acoustic flares (cruise report CAGE 15e2, CAGE 16e5, Serov et al., 2017) and GH layers interbedded in sediments recovered from core 1520GC (at 150 cm bsf). The multi-proxy approach on the three cores collected at GHM 3 (1522, 1520 and 940GC) confirms the activity of the methane system, characterized by highly complex heterogeneous fluid and gas dynamics.

Core 1522GC
The core 1522GC, 320 cm in length, was retrieved from the western rim of GHM 3 (Fig. 1) and records a period of 20 kyr (Fig. 7). The trend of the pore water sulfate (Fig. 7a) suggested that the SMTZ was located at the very bottom of the core (~320 cm) or slightly below, indicating a very low methane flux at this site (Hong et al., 2017). There is no Ba/Ti anomaly and background concentrations and d 13 C values of archaeal and bacterial lipid biomarkers in the upper part of this core (Fig. 7) indicate a low methane flux settings. Only at the bottom of the core (280e320 cm) the low d 13 C value of À75‰ of biomarkers (Fig. 7), provides indications for AOMrelated archaea at the SMTZ. Similar to lipids, d 13 C-values of N. pachyderma were in the range typical for normal, non-seep marine environments (À0.5 to 0.5‰) in the Barents Sea (Knies and Stein, 1998). In the lower section at 300e320 cm, the d 13 C values of N. pachyderma became low (down to À3.5‰; Fig. 7d), however, we could not exclude that the more negative values of N. pachyderma are due to changes in the ambient water.

Core 1520GC
The geochemical profiles of core 1520GC (Fig. 8a) collected at the southern rim of GHM 3 (Fig. 1) showed pore water sulfate depletion located at a sediment depth of~100 cm (event A). Authigenic carbonate nodules were present throughout the entire core, indicating that this site was affected by continuous, pervasive, diffuse methane emissions for a long time. This is in line with the lack of prominent Ba/Ti anomalies (Fig. 8b). Peaks in Ba/Ti did not develop because of the highly dynamic methane flux, which did not leave enough time for the generation of a stable barite front (although unlikely, the missing peak could also be related to a technical artefact, i.e. that the XRF sensor lacked proper contact to the sediment core during the XRF acquirement because of the abundant carbonate nodules). Even if this core has been subject to a continuous and pervasive methane seepage since the LGM, three distinct events (B, C, and D) can be identified in addition to the present-day SMTZ (event A), that might indicate more advective methane flow. The event B, located above the current SMTZ shows high concentrations of both archaeal (up to 1 mg/g dry sediment) and SRB lipid biomarkers (up to 1.2 mg/g dry sediment), negative excursions of d 13 C values of planktonic and benthic foraminifera down to À18‰, and visual observation of abundant carbonate nodules (Fig. 8). The events C and D below the SMTZ are characterized by the presence of numerous carbonate nodules, d 13 C Fig. 7. Core 1522GC, a. pore water sulfate from Hong et al. (2017), b. sediment elemental ratio Ba/Ti, c. magnetic susceptibility, d. δ 13 C of planktonic foraminifera (N. pachyderma), e. archaeal lipid biomarker concentration, f. archaeal lipid biomarker δ 13 C, g. bacterial lipid biomarker concentration and h. bacterial lipid biomarker δ 13 C. A tentative event of past methane emission in brown. Yellow and green shaded bars highlight the horizon markers in this core. All isotope and lipid biomarker data are available in the supplementary file.
(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) negative excursions of benthic and planktonic foraminifera (Fig. 8c), and 13 C-depleted archaeal lipid biomarkers (Fig. 8e). The archaeal and SRB lipid biomarkers concentration are much lower in events C and D (Fig. 8d), and the d 13 C values of SRB lipid biomarkers are not as depleted as in events A and B (Fig. 8e,g). These differences are likely due to degradation of the bacterial lipid biomarkers (Niemann et al., 2013).
The sequence of events C/D, events B, and A suggests that the methane flux was high in the past and has decreased recently as the SMTZ shifted from B to A. Hong et al. (2017) modeled the methane flux based on the pore water profile and suggested that it took 290 years of an increasing methane flux to achieve the current geochemical profiles at site 1520GC. This is not contradictive to our interpretation of long term decreasing methane flux, as the time scale is very likely different. The time scale in the modeled result from Hong et al. (2017) is hundreds of years, whereas the time scale for generating events A and B could be thousands of years, based on the carbonate nodule formation, lipid biomarker concentrations etc. The methane flux fluctuated during these thousands of years to produce even smaller events that cannot be resolved by our approach due to sampling resolution (every 5 cm) compared to the pore water modeling was based on the sampling resolution of 2 cm.

Core 940GC
In core 940GC (Fig. 9), from the eastern rim of GHM 3, we could only apply the Ba/Ti sediment ratio, magnetic susceptibility, and d 13 C values of carbonate nodules and planktonic foraminifera due to sampling availability. The sulfate obtained from pore water measurement is depleted at 210 cm bsf, indicating the current position of the present day SMTZ (assigned as event A). The presence of Ba/Ti peak in event A suggests a continuous decreasing methane flux at this site, as only during a decreasing methane flux, the previous Ba/Ti peaks above can be preserved as we outlined in the introduction (Kasten et al., 2012).
Two major Ba/Ti peaks indicate two other events, B and C. In both events B and C, authigenic carbonate nodules co-occur with negative excursions of d 13 C values of N. pachyderma down to À19‰ in event B and À16‰ in event C. However, the d 13 C values of N. pachyderma was À0.3‰ at the current SMTZ (event A), and it cooccurs with carbonate nodules with d 13 C values of À33.5‰. One of the possible reasons for the discrepancy could be related to the selection of foraminifera picked for isotopic investigation, because the glassy specimens are less affected by diagenetic alteration compare to the frosty ones . Nevertheless, a glassy appearance of foraminifera does not guarantee any absence of diagenetic influence Schneider et al., 2017).
From both mounds, our analyses revealed a major methane flux increasing phase, e.g. event D in core 920GC, event B in core 1070GC, and event D in core1520GC. We speculate that the major seeping events were in phase with the deglaciation that triggered changes in the GHSZ  as the host sediment in these events co-occur with the inception of the deglaciation. The increased or decreased fluid flow within the mounds was most probably controlled by the variations in the depth of the GHSZ and the associated GH dissociation Waage et al., 2019), which is commonly seen in other areas in the Arctic (Cr emi ere et al., 2016;Wallmann et al., 2018). The complexity of spatial heterogeneity in a small area as investigated here is likely caused by carbonate precipitation and/or GHs that can block upward transport routes of gasses and fluids (Treude et al., 2020). Thus, each coring site might have experienced a different fluid/gas discharge history.

Summary, conclusion and perspective
A profound understanding of methane emissions from Arctic Ocean sediments in the past and at present, and the associated drivers and trigger mechanisms can be the key to predict methane emission scenarios in a future Arctic Ocean impacted by climate change. This can be achieved by tracing evidence for past methane dynamics preserved in fossil records e i.e., by analysing proxy indicators allowing to reconstruct past methane dynamics. In this study, we provide an overview of a set of (bio) geochemical proxies (pore water sulfate and methane profiles, sedimentary ratios of Ba/ Ti, magnetic susceptibility, d 13 C records from benthic and planktic foraminifera, mineralogy and d 13 C of authigenic carbonates, abundance and d 13 C of lipid biomarkers). We applied this multi-proxy approach to the case of a shallow water (~380 m water depth) GH area, Storfjordrenna, to reconstruct the emission history of two GHMs in an Arctic cold seep system. This combined approach allowed to establish stratigraphic correlations between the cores and to identify several major and minor methane venting phases at each coring site, from both the active and inactive GHMs. Our data revealed a major methane venting phase at both mounds in sediments that are~21 kyr old. i.e. deposited at the time of the glacial ice sheet retreat from the area of Storfjordrenna. Also previous studies from other Arctic GH systems found indications for elevated methane seepage related to the retreat of the glacial ice shield (Wallmann et al., 2018). For example, in Prins Karls Forland (PKF) which has a similar water depth than Storfjordrenna, seismic and modeling results suggest the seepage association with the deglaciation and changes in GHSZ due to isostatic rebound (Wallmann et al., 2018). In another GH system, Vestnesa Ridge, where the water depth is 1200 m, it was postulated that the methane seepage was controlled by the last deglaciation that triggered changes in the Table 3 Summary of conventional proxies used for methane seepage reconstruction and applied in this study. the discrepancy at the current SMTZ Panieri et al., 2014Consolaro et al. (2015); this paper etc.

Proxy
Abundances and species distribution High abundance of endobenthic foram and high absolute abundance of benthic foram were interpreted as indicating methane seepage. Panieri and Sen Gupta, 2008;Heinz et al., 2005;Bernhard and Sen Gupta 1999;Bernhard et al., 2001 Shell preservation Diagenetically altered shells indicate precipitation of authigenic phases at a depth of the SMTZ at sites of methane seepage fault system and favoured gas migration (Schneider et al., 2018;Himmler et al., 2019). The multi-proxy approach that was followed in the present study provides a base and directions for future investigations in cold seep systems, giving insights and guidance for choosing the appropriate proxy to use depending on the environmental conditions. Table 3 summarizes the principle of each proxy applied in this study and evaluates the advantages and disadvantages of each proxy for reconstructing methane seepage. This guideline will be useful for future studies to reconstruct methane seepage, to sample and apply proper available proxies during the investigation. Applying our multiproxy approach revealed a differential seepage history across small spatial scales. This further highlights the need to investigate seepage history from multiple cores.

Author statement
HY collected the samples and did the experiments and analyses under supervision with GP and HN. GP provided support with sample collection, foraminiferal analyses and HN provided support for lipid biomarker experiments and data analyses. All authors contribute to the writing of the manuscript at all stages.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.