Methane fluxes in marine sediments quantified through core analyses and seismo-acoustic mapping (Bornholm Basin, Baltic Sea)

Accepted Manuscript Methane fluxes in marine sediments quantified through core analyses and seismo-acoustic mapping (Bornholm Basin, Baltic Sea) Karen Marie Hilligsøe, Jørn Bo Jensen, Timothy G. Ferdelman, Henrik Fossing, Laura Lapham, Hans Røy, Bo Barker Jørgensen PII: DOI: Reference:

S0016-7037(18)30426-5 https://doi.org/10.1016/j.gca.2018.07.040 GCA 10871

To appear in:

Geochimica et Cosmochimica Acta

Received Date: Revised Date: Accepted Date:

24 January 2018 16 July 2018 29 July 2018

Please cite this article as: Marie Hilligsøe, K., Bo Jensen, J., Ferdelman, T.G., Fossing, H., Lapham, L., Røy, H., Barker Jørgensen, B., Methane fluxes in marine sediments quantified through core analyses and seismo-acoustic mapping (Bornholm Basin, Baltic Sea), Geochimica et Cosmochimica Acta (2018), doi: https://doi.org/10.1016/ j.gca.2018.07.040

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Methane fluxes in marine sediments quantified through core analyses and seismo-acoustic mapping (Bornholm Basin, Baltic Sea)

Karen Marie Hilligsøe1, Jørn Bo Jensen2, Timothy G. Ferdelman3, Henrik Fossing4, Laura Lapham5, Hans Røy1 and Bo Barker Jørgensen1*

1

Center for Geomicrobiology, Department of Bioscience, Aarhus University, Ny Munkegade 116,

DK-8000 Aarhus C, Denmark. 2

Department of Marine Geology and Glaciology, Geological Survey of Denmark and Greenland,

C.F. Møllers Allé 8, DK-8000 Aarhus C, Denmark. 3

Department of Biogeochemistry, Max Planck Institute for Marine Microbiology, Celsiusstr. 1, D-

28359 Bremen, Germany. 4

Marine Ecology Section, Department of Bioscience, Aarhus University, Vejlsøvej 25, DK-8600

Silkeborg, Denmark. 5

Chesapeake Biological Laboratory, 146 Williams Street / 0038, Solomons, MD 20688, USA.

*Corresponding author E-mail address: [email protected] (B.B. Jørgensen)

Key words: Sediment echosounder, methane gas bubbles, diffusion, sulfate-methane transition, SMT, Holocene mud layer, free gas depth, HML model, FGD model, organic carbon degradation, methane budget

Abstract

Methane is a terminal product of anaerobic degradation of organic matter in subsurface marine sediments depleted of reactive oxidants. The depth and age of the sediment where sulfate is depleted determines the extent of methane production relative to the burial of organic carbon. We aimed to understand how this methane production is controlled and distributed in an apparently uniform sediment basin. We combined seismo-acoustic surveys and geochemical analyses of sediment cores to explore the distribution of methane fluxes in brackish-marine mud deposits of the Bornholm Basin, southern Baltic Sea. Geophysical mapping revealed the depth distribution of a) the thickness of the Holocene organic-rich mud layer overlaying organic-poor Postglacial clay, and b) the upper boundary of methane gas bubbles trapped in the Holocene sediment. By correlating these two parameters with the methane distributions in sediment cores from 44 stations we developed algorithms to estimate and map, at high spatial resolution, the diffusive methane fluxes in the 75-95m deep and >4,000 km2 large Bornholm Basin. The two approaches, termed the FGD (Free Gas Depth) model and the HML (Holocene Mud Layer) model, yielded similar budgets for the total upwards methane flux through the sediment column in the Bornholm Basin, about 17 ton methane C day-1. Complete sulfate depletion at depth, and thus onset of methane production, required a minimum threshold thickness of the Holocene mud layer of 4 m. Although the Bornholm Basin has an even depth contour and uniform surface sediments, methane production was strongly focused in hotspots where the HML has greatest thickness. This heterogeneity could not be predicted from bathymetry or from properties of the surface sediment but was related to the topography of the glacial landscape buried underneath the Holocene mud blanket. This demonstrates the importance of including seismo-acoustic mapping of subsurface stratigraphy for the interpretation and geographic extrapolation of sediment core data for biogeochemical processes such as methane production and methane flux. The two approaches, which were here combined for the first time, may thereby be applied to map methane production also in other coastal and shelf sediments with shallow gas and distinct Holocene deposits.

1. INTRODUCTION

Production of methane in the seabed accounts for an important fraction of global biological methanogenesis, yet most of the methane is retained and degraded to CO2 in the subsurface sediment (Reeburgh, 2007). Methane in the atmosphere is a strong greenhouse gas and abrupt warming events during the past 100,000 years have repeatedly been associated with strong increases in atmospheric methane (Brook et al., 1996). A potential source of this methane is the great quantities of gas hydrates bound in the seabed (Wallmann et al., 2012). Methane from gas hydrates or from shallow gas in the coastal seabed may thereby affect future climate development (Kvenvolden 1988; Best et al., 2006). The awareness of this methane source has been enhanced in recent years by observations of methane ebullition from Arctic shelf sediments, apparently induced by recent ocean warming (Shakhova et al., 2010). The geochemical and biological controls on methane production in the seabed therefore need to be understood and quantified, both on a regional and a global scale.

Organic matter, which is buried in marine sediments, is exposed to microbial degradation and mineralization at a rate which steadily decreases with age and depth beneath the sediment surface (Middelburg, 1989). Sulfate reduction is the prevailing terminal step in the anaerobic mineralization down to a depth where sulfate becomes depleted. Below that depth, methanogenesis is the main terminal pathway and leads to the accumulation of methane in the pore fluid. Due to the decreasing reactivity of the organic matter as it is buried and finally reaches the methane zone, methanogenesis is generally only <1-10% of the anaerobic organic matter mineralization in marine sediments (Jørgensen and Parkes, 2010; Flury et al., 2016; Egger et al., 2018). Yet, the opposed pore water gradients of downwards diffusing sulfate and upwards diffusing methane are very distinct and highlight the anaerobic oxidation of methane (AOM) with sulfate, which takes place at the sulfatemethane transition (SMT) (e.g. Niewöhner et al., 1998; Borowski et al., 1999). This SMT is generally defined as the zone where sulfate and methane gradients overlap, yet the transition from sulfate to methane can be quite variable and the zone is not always well constrained. For practical purposes, the specific depth of the SMT is in this study defined as the depth where sulfate and methane occur in equi-molar concentrations (cf. Dale et al., 2008; Flury et al., 2016). This corresponds roughly to the depth where the linear extrapolation of sulfate and methane concentration gradients above and below the SMT, respectively, hit zero.

A consequence of the steep and continuous decrease in reactivity of organic matter with depth in the seabed is that most methanogenesis takes place in the uppermost part of the methane zone (Beulig et

al., 2018). This is where the methane concentration gradient is the steepest and where the methane therefore diffuses upwards to the SMT. The methane flux to the SMT thereby comprises much or most of the methane produced in the entire sediment column below (Jørgensen and Parkes, 2010; Flury et al., 2016). Under conditions where gas bubbles or gas hydrates accumulate, these reservoirs store a variable fraction of the produced methane (Mogollón et al., 2012; Wallmann et al., 2012). Shallow free gas, i.e. small gas bubbles entrained in the unconsolidated upper sediment column, occurs widespread in coastal marine environments, in particular in areas with rapid accumulation of fine-grained mud rich in organic matter, such as estuaries, fjords, and bays (Fleischer et al., 2001).

Current estimates of the methane production in the global seabed are based primarily on data from multiple sediment cores in different geographic regions or depth zones of the ocean and from transport-reaction models based on such data (Reeburgh, 2007; Regnier et al., 2011; Egger et al., 2018). Similar to sulfate reduction, such estimates show a strong focusing of methanogenesis towards the ocean margins, especially towards near-coastal sediments with high deposition rates and high organic carbon content (Canfield et al., 2005). The Postglacial sediment accumulation on the inner continental shelf and in fjords and estuaries is temporally and spatially heterogeneous, however, which makes it difficult to perform quantitative estimates of the regional methane production rates (Egger et al., 2018).

Shallow-seismic transects through the Baltic Sea reveal its Postglacial history of sedimentation, which is particularly distinct in the stratigraphy of the deeper basins where deposition has been continuous (Moros et al., 2002; Andrén et al., 2015a). The combination of such transects with sediment coring has been used for mapping and quantification of methanogenesis in near-coastal sediments (e.g., Dale et al., 2009). The seismo-acoustic data may generate transects or maps showing sediment deposition rates, while sediment core analyses may couple such deposition rates to the occurrence and depth distribution of methane. Shallow gas has been recorded by acoustic methods in the Bornholm Basin and in many other areas of the Baltic Sea (Laier and Jensen, 2007; Schneider von Deimling et al., 2013; Tóth et al., 2014a). Such shallow gas in muddy sediments can be used to detect the geographic hot-spots of methane accumulation (Martínez-Carreño and Garcia-Gil, 2013), while the depth of the gas front can be used to model the rate of methane production (Mogollón et al., 2012, 2013).

We used such a combination of dense shallow-seismic surveys and data from multiple sediment cores to map subsurface methane fluxes in the Bornholm Basin. The goals were to establish a

methane budget for an entire basin and to understand the depositional controls on the distribution and intensity of methane production. The high geographic density of sediment coring and methane analyses enabled two different approaches for the quantification of the overall methane flux. One approach used an empirical relationship between the thickness of the Holocene mud layer and the calculated methane fluxes to make a geographic extrapolation using the shallow-seismic data. Another approach used a diffusion model that coupled the calculated methane fluxes in sediment cores to the depth of the gas bubble front. Both approaches yielded similar high-resolution maps of a very heterogeneous distribution of methane production, and both models resulted in similar estimates of the total methane flux in the basin. These approaches can be applied to a wide range of sediments of fjords, estuaries and inner continental shelf regions, where Holocene deposits overlaying glacial till and clay can be identified by shallow seismics and where the sediment may have accumulated shallow methane bubbles.

1.1.

The Bornholm Basin

The Baltic Sea with its 380,000 km2 is one of the largest brackish water bodies on Earth. The catchment area is 1,700,000 km2 with mostly forest and wooded land in the north and agricultural land to the south. The Baltic Sea is connected to the North Sea via a narrow passage through the Danish straits. The present Baltic Sea Basin was shaped during the last ice age by glaciers from the Fennoscandian Ice Sheet. As the ice started to melt and retreat about 14,000 yr BP (years before present), a large freshwater lake developed in front of the ice cap and deep deposits of Postglacial clay were laid down in the Baltic Ice Lake. A first marine transgression into the southern Baltic Sea occurred about 11,600 yr BP, resulting in the brackish Yoldia Sea. Due to the changing balance between eustatic sea level rise and isostatic uplift of the Scandinavian land mass as the ice cap melted (Jensen et al., 2017), the connection to the open sea was again cut off and a new limnic stage, the Ancylus Lake, gradually developed around 11,000 yr BP (Andrén et al., 2000). A second marine transgression finally established permanent brackish-marine conditions around 8,600 yr BP as the current gateways to the North Sea opened through the Danish Belts and gave rise to the Littorina Sea, i.e. the modern Baltic Sea (Björck, 1995).

An abrupt increase in the organic carbon content in the Bornholm Basin sediment shortly before 8,000 yr BP reflects a general increase in the primary productivity in the Baltic Sea in response to the stronger marine influence and to a warmer climate (Andrén et al., 2000). This lithological boundary marks the transition from limnic clays to the Holocene brackish-marine mud deposit, which is here

termed the Holocene Mud Layer (HML). High rates of sulfate reduction and methanogenesis occur in the HML whereas the underlying Postglacial clay has very little organic matter and low organic degradation rates (Holmkvist et al., 2014; Beulig et al., 2018).

We studied the central part of the Bornholm Basin, NE of the island Bornholm (Fig. 1A). The basin has a uniform bathymetry with water depths of mostly 75-95 m and an area of >4,000 km2. The organic-rich Holocene mud layer covers a rugged landscape of Quaternary deposits which reflect the earlier faulting and depositional patterns in this pull-apart basin formed during the late Cretaceous and early Tertiary (Jensen et al., 2017). The Basin is bounded by major faults and has been interpreted as a half-graben in which the Holocene mud deposits are thickest near the fault scarp to the southwest and are thinning towards the deepest, central part of the basin in the northeast (Jensen and Endler, 2012). During the last ice age, massive glaciers first eroded the bottom down to the Cretaceous bedrock. They then left glacial tills that covered the bedrock and reshaped the topography. In Lateglacial and Postglacial time, deep organic-poor sediments of glacio-lacustrine clay were deposited in the deep-water basin of the Baltic Ice Lake, which had only weak bottom currents. The subsequent depositional sequences have been studied through a combination of seismoacoustic mapping and analysis and dating of multiple sediment cores (e.g., Kögler and Larsen, 1979; Björck, 1995; Jensen et al., 1997; Moros et al. 2002; Andrén et al., 2011, 2015b). As a result, the chronology of major historical events during the Postglacial development of the Baltic Sea can now be identified from acoustic reflectors in the sediment. Current sedimentation rates in the Bornholm Basin have been determined by the 210Pb method (Christoffersen et al., 2007) and by the analysis of a 137Cs peak linked to the Chernobyl reactor catastrophe in 1986 (Povinec et al., 2003; Christoffersen et al., 2007). Sedimentation rates are highest, 500-1100 g m-2 yr-1 (2-4 mm yr-1), in the deep, central and western part of the basin and generally lower in the periphery of the basin. The present inflow of bottom water brings sediment material into the Bornholm Basin, in particular from the west from the shallower Arkona Basin (Fig. 1A) (Christoffersen et al., 2007). The inflow is focused along the southwestern margin of the Bornholm Basin where a trough cuts through the uppermost HML. The total organic carbon (TOC) content is 4-6% d.wt. (dry weight) in the deep basin with a slight increase over the last 20-40 years due to eutrophication (Christoffersen et al., 2007).

A halocline at 50-80 m depth persists between the shallower brackish water of the Baltic Sea and the more saline bottom water originating from the North Sea (Väli et al., 2013). Major inflows of saline

and oxygenated bottom water occur irregularly and the deep basin is generally hypoxic with seasonal anoxia (e.g., Carstensen et al., 2014a). Lack of lamination indicates that some bioturbation does occur, at least after periods of oxic seawater inflow (Andrén et al., 2015b). The annual mean temperature of the bottom water, and thus of the sediments, in the Bornholm Basin recorded during 1998 to 2010 was 7°C with little variation. The mean salinity was 17 (Rak and Wieczorek, 2012).

2. MATERIALS AND METHODS

Pore water data for sulfate, methane and other relevant chemical and physical parameters were obtained from sediment cores taken at 44 stations in the Bornholm Basin (Fig. 1C). The coring positions were chosen along acoustic transects by which the stratigraphy and the presence or absence of shallow gas had been recorded (Fig. 1B). Their locations were selected to provide a broad distribution of gas depths and methane gradients. The center of the basin was avoided during coring because of toxic munitions dumped here after World War II. Coring tools included gravity corer, Rumohr corer (Meischner and Rumohr, 1974), and multiple corer, often in combination in order to obtain samples and data from the deeper sediment (down to 6-12 meter below seafloor, mbsf) together with high-resolution data from undisturbed surface sediment (top 0-1 mbsf). The sulfate concentration profiles were aligned between the shallow and the deep cores in order to correct for the loss or disturbance of surface sediment by the core catcher in the gravity corer.

Cores and shallow-seismic transects were obtained during the following five research cruises:

1) RV Gunnar Thorson, August-September 2004 (GT04) 2) RV Poseidon, August-September 2009 (PO392) 3) RV Maria S. Merian, July-August 2010 (MSM16.1) 4) RV Aurora, June 2014 (AUBO14) 5) RV Aurora, June 2016 (AUBO16)

Geochemical data from expeditions PO392 and MSM16.1 are deposited in the public online database, PANGEA (https://www.pangaea.de/), and can be retrieved under the project names, METROL and BALTIC GAS. Data from expeditions AUBO14 and AUBO16 were published for a different purpose by Beulig et al. (2018).

The geographic positions and other relevant parameters for all stations are given in the Supplementary Information Table S1.

2.1.

Analyses

The procedures for sampling and analysis of the sediment cores changed in some of the details during the 12-year interval between the first and the last research cruise. Procedures used by the METROL project during the 2004 and 2009 cruises are described by Parkes et al. (2007) and Knab et al. (2008). Procedures used during the BALTIC GAS cruise in 2010 and the Aarhus University Center for Geomicrobiology cruises in 2014 and 2016 follow those described in Flury et al. (2016) and Beulig et al. (2018). In the following we only summarize the procedures.

Immediately after a gravity core was on deck, it was cut in 1-m sections, small holes were drilled in the PVC core liner at 25-30-cm depth intervals, and sediment samples for CH4 analysis were taken with 3 mL cut-off syringes. Methane was analyzed immediately as described below in order to direct the drilling of further holes for methane sampling at high depth resolution around and below the sulfate-methane transition for more accurate determination of the methane diffusion gradient. The first, fast sampling was also used to check for possible methane loss during the more detailed, but slower sampling. This comparison showed no systematic methane loss between the two rounds of sampling as long as the in situ methane concentration of the sediment was not above 1 bar partial pressure.

During the 2004 and 2009 cruises, the entire sediment core was gradually pushed out of the core liner sections with a piston and subsamples for all further analyses were taken from the core end, which became exposed step by step. Several 5 mL samples were taken with cut-off syringes for porosity and total organic carbon (TOC) while 3-cm wide and 20-cm long core liners were pushed in and the sediment used for pore water sampling. Pore water was obtained by squeezing through a 0.45 µm filter under N2 directly into a 1% w/w Zn-acetate solution. During the 2010, 2014 and 2016 cruises, pore water was obtained by Rhizon soil moisture samplers (Rhizosphere Research Products, Wageningen, The Netherlands) inserted through 4-mm wide holes drilled in the core liner. Windows of roughly 50 cm2 were then cut in the side of the core liner using a vibrating saw. A ca 1-cm deep layer of contaminated surface sediment was carefully scraped away and further subsampling of the exposed sediment was done for porosity, TOC and other parameters.

For methane analysis, the 3 mL of sediment was immediately transferred to 20 mL glass vials containing 6 mL 2.5% NaOH or saturated NaCl and crimp capped with butyl rubber stoppers. Samples were shaken thoroughly and stored upside down at 4°C until analysis within a few days. A 100 µL headspace gas sample was injected with a Hamilton glass syringe into a gas chromatograph equipped with a packed, stainless steel Poropak-Q column and a flame ionization detector (HewletPackard or SRI GC) (Knab et al., 2008; Flury et al., 2016). In cores, where sulfate had not been depleted, the methane concentration did not exceed 5-10 µM, i.e. about 1000-fold lower than in sulfate-depleted cores with methane accumulation. In the following we label this "low methane".

In pore water samples for sulfate analysis, the hydrogen sulfide (H2S) was immediately removed to prevent its oxidation to sulfate. Two different approaches were used. During the 2004 and 2009 cruises the sulfide was precipitated with Zn-acetate and the precipitated ZnS removed by filtration or centrifugation. During the 2010, 2014 and 2016 cruises the pore water was acidified with HCl and H2S was stripped out by bubbling with humidified air or by shaking on a Vortex mixer. Alternatively, humidified CO2 gas was used to simultaneously acidify the sample and strip out H 2S. The SO42- was analyzed by ion chromatography (Waters or Dionex IC). Further cores, which are described elsewhere, were taken for sedimentological analyses. Multisensor core logging of 6–12-m long gravity cores was used to determine the physical properties of sediments with and without gas bubbles. This information was used for the interpretation of the seismo-acoustic records (Jensen and Endler, 2012). Split cores were also used for general core description, sub-sampling and sedimentology. Temperature, porosity and total organic carbon (TOC) were measured as described in the Supplementary Information.

Whole sediment diffusion coefficients of sulfate and methane were calculated from tabulated molecular diffusion coefficients corrected for tortuosity (θ), in situ temperature and salinity (Iversen and Jørgensen, 1993; Schulz, 2006). The correction factor used for turtuosity was θ2 = 1-ln(φ2) (Boudreau, 1996), where φ is sediment porosity. Temperature and salinity of the sediment were measured during each expedition. However, the deep methane diffusion coefficient will depend also on long-term variations in temperature which were not recorded but can be estimated from published hydrographic records. The bottom water temperature in the Bornholm Basin varies throughout the year, depending on the season and, in particular, on the inflow of saline bottom water. Over the decade of our study the bottom temperature was in average 7°C with a standard deviation of ±2°C

(Rak and Wieczorek, 2012). This small temperature variation is further dampened at the depth of the SMT (Dale et al., 2008).

2.2.

Shallow seismics

A detailed mesh of shallow-seismic lines (Fig. 1B) was used to develop stratigraphic maps of the Bornholm Basin. The data were recorded during the five expeditions and supplemented with archive data from the Geological Survey of Denmark and Greenland (GEUS) and the Leibniz Institute for Baltic Sea Research Warnemünde (IOW). During the multiple cruises and by the archive data, a broad range of acoustic techniques were used: single-beam echo sounder (12-200 kHz), sediment echo sounder (5-40 kHz), single channel Boomer (2-4 kHz), single channel Sparker (1 kHz) and multichannel Airgun (0.2 kHz). The highest frequencies provide the highest resolution and enable the determination of the water depth with an accuracy of a few cm (Wunderlich and Müller, 2003). The lower frequencies penetrate deeper into the sediment and provide an acoustic image of the sediment strata. Methane gas bubbles effectively absorb and reflect the acoustic energy, and the attenuation of the return signal from below is seen as acoustic blanking of the underlying gascharged sediment (e.g., Thiessen et al., 2006; Toth et al., 2014b, 2015). Thus, the upper front of free methane gas can be observed as a conspicuous blanking on the seismic images. Low-frequency highresolution multichannel seismics, which could also penetrate the gas front, was used to provide information on the geological stratification of deeper structures (Jensen et al., 2017). Due to acoustic blanking it was generally not possible from the high-frequency acoustics to see how deep down in the sediment methane bubbles occurred.

2.2.1. Thickness of the Holocene mud (HML)

A detailed correlation of the acoustic data and the Holocene stratigraphy was performed utilizing numerous gravity cores, sediment analyses and radiocarbon dating (Jensen et al., 2017). The IODP Leg 347 Sites M0065 and M0066 (Andrén et al., 2015a,b) provided additional documentation of the Glacial, Lateglacial and Postglacial sedimentary units. A sound velocity of 1450 m/s was chosen for the conversion of two-way-travel-time (TWT) to depth, based on IODP average velocity measurements of the Holocene marine sediments (Andrén et al., 2015b). A map showing the thickness of the Holocene mud layer (HML) was constructed based on the digitized depth of the seabed and the bottom reflector of the Holocene marine mud. All the mentioned acoustic techniques were included, but in areas with high frequency acoustic blanking only the low frequency data were

used (Toth et al., 2014a, 2015). The vertical resolution in a map of the HML thickness generally depends on the thickness of the layer, the accuracy of the velocity determination, and the frequency of the acoustic system. In our case, with a maximum HML thickness of 20 m, with good IODP in situ velocity measurements, and with an equipment frequency of around 1 kHz, the HML thickness accuracy is estimated to be 0.5 m.

2.2.2. Free gas depth (FGD)

We define the upper front of free methane gas as the free gas depth (FGD) according to Dale et al. (2009). We digitized the FGD in seismic transects and developed a gas distribution map with indication of the depth from the seafloor to the top of the acoustic gas front. The field data were loaded onto a seismic work station and combined with seismic archive data from the same area. The total seismic dataset was interpreted and combined with physical characteristics of the sediments to compile a map of the gas distribution. The vertical resolution of the free gas depth determination based on sediment echosounder data (5-40 kHz) is 0.2-0.3 m (Wunderlich and Müller, 2003).

2.3.

Models

Our study combined multiple sediment cores, in which the depth distribution of methane was analyzed, with multiple acoustic transects, from which a) the thickness of the Holocene mud layer and b) the upper front of gas bubbles could be mapped. We used this information to calculate and map, with a 100×100 m grid, the methane fluxes in the Bornholm Basin using data from either a) or b). In the following, the two independent approaches are explained.

2.3.1. Holocene mud layer (HML model)

Based on the sediment stratigraphy, visible in cores and recognizable in acoustic transects, we could demonstrate a tight systematic relationship between the thickness of the HML and the degree of methane accumulation in the sediment. We used this statistical relationship, which we call the HML model, to calculate methane fluxes from the HML thickness for the entire study area.

2.3.2. Free gas depth (FGD model)

We analyzed methane distributions in such a large number of sediment cores that the calculated methane fluxes could be directly correlated to the upper boundary of gas bubbles (the FGD). Such a correlation is site specific, however, as it depends on water depth, salinity, temperature and other parameters. We therefore developed a mechanistic model for the calculation of methane fluxes in the sediment based on the measured position of the FGD that could be generalized also to other sediment locations (Fig. 2C). We thereby made the following simplifying assumptions about sulfate and methane and their diffusion above the gassy sediment: 1) Linear sulfate gradients from the sediment surface to the SMT. A survey of sulfate distributions in Bornholm Basin sediment cores showed that this is a realistic assumption in the gassy sediment as bioirrigation is low in the sediment areas with gas because of low oxygen concentration in the overlying bottom water (Carstensen et al., 2014b). Thus, a shallow SMT tended to correlate with linear gradients (Beulig et al., 2018; Egger et al., 2018). 2) Linear methane gradients from the SMT and down to the FGD. Due to high supersaturation and loss of methane after retrieval of the sediment cores on deck, it was only possible to analyze methane correctly in the upper part of the methane zone, just beneath the SMT. However, the methane gradient in this depth interval was well defined from our measurements and could be extrapolated linearly downwards to estimate the depth where it might exceed the in situ hydrostatic pressure and generate bubbles. 3) The upper limit of gas bubbles trapped in the mud should correspond to the depth where the total gas pressure equals the ambient hydrostatic pressure so that methane bubbles form in situ. The hydrostatic pressure was calculated from the bathymetry plus the free gas depth within the sediment, assuming 1 bar pressure increase per 10 m total depth. We thereby assume that the sediment did not have overpressure (e.g. Lafuerza et al., 2009). Other dissolved gases in the pore water, such as N2, will equilibrate with the gas phase of the bubbles. Out of the 1 bar atmospheric pressure at the sea surface, 0.8 bar is due to N2, which is also present in the methane bubbles. Of other relevant gasses, O2 was consumed at the sediment surface while the CO2 partial pressure, based on sediment pH and maximum DIC concentration at the FGD, was only 0.035 bar. The 0.8 bar of N2 corresponds to 8% of the total gas pressure in the methane bubbles at a water depth of 75-95 m. This was taken into account when calculating the methane partial pressure at the FGD and, thus, the pore water methane concentration in equilibrium with that partial pressure. If gas bubbles migrate upwards slowly from below, they expectedly dissolve again at the FGD and therefore do not affect the model calculations.

4) The downward flux of sulfate, JSO4, and upward flux of methane, JCH4, to the SMT (both expressed numerically) satisfy a 1:1 stoichiometry of anaerobic methane oxidation by sulfate:

(1)

(2)

Based on these four assumptions we developed a FGD model as follows. The upward flux of methane to the SMT was determined from Fick’s first law of diffusion:

(3) where φ is sediment porosity, DCH4 is the whole sediment diffusion coefficient of methane, [CH 4] is the methane concentration, z is the depth beneath the sediment surface, and d[CH4]/dz is the methane concentration gradient.

The downward flux of sulfate was similarly calculated as (numerically):

(4)

where DSO4 is the whole sediment diffusion coefficient for sulfate corrected for temperature, salinity and tortuosity (Schulz, 2006). Given the assumption about linear gradients, the sulfate gradient is the difference between the seawater sulfate concentration at the sediment surface, [SO42-]SW, and the near-zero sulfate concentration at the SMT depth divided by the depth of the SMT beneath the sediment surface, ZSO4 (Fig. 2C). Thus, the sulfate flux, which is here assumed equal to the methane flux, is:

(5)

The methane gradient is the difference between the near-zero concentration at the SMT and the in situ saturation concentration at the FGD, [CH4]FGD, taking the N2 partial pressure into account. This saturation concentration can be determined from the temperature, salinity and water depth plus the

depth in the sediment (Yamamoto et al., 1976; Duan et al, 1992). It is located at a distance, ZCH4, beneath the SMT (Fig. 2C). The methane flux is thus:

(6) The free gas depth, ΔZFGD, is composed of the depth interval of the sulfate zone, from the sediment surface and down to the SMT, ZSO4, plus the depth interval from the SMT and down to the FGD, ZCH4. We use ΔZFGD = ZSO4 + ZCH4 and add together Eq. (5) and (6). The methane flux, JCH4, can then be expressed as:

(7)

According to this simplified relationship, which we call the FGD model, the methane flux can be determined from the upper front of bubbles, i.e. the FGD as recorded by acoustic methods, provided that the following simple parameters are known: water depth, bottom water temperature and salinity, and sediment porosity.

According to the FGD model in Eq. (7), the methane flux in gassy sediments should be inversely proportional to the free gas depth (ΔZFGD), provided that other parameters are constant. These other parameters are: 1) Sediment porosity, which varied with depth in the sediment but did not show systematic variations between the different coring sites and water depths, as long as gas bubbles were not present (Fig. S1B). 2) Whole core diffusion coefficients, which relate to sediment porosity and are particularly sensitive to temperature. We recorded temperature and salinity of the sediment during each expedition and used the data to calculate diffusion coefficients for that period. 3) The sulfate concentration at the sediment surface, [SO42-]SW, which is a linear function of bottom water salinity, S: [SO42-]SW = 28.6 × S/35 mM. The mean bottom water salinity over the decade of the study was 17 with a standard deviation of ±1 (Rak and Wieczorek, 2012). 4) The methane concentration at the FGD, [CH4]FGD, which depends linearly on the ambient hydrostatic pressure, which again depends on the total depth beneath the sea surface, i.e. the sum of the water depth above the sediment and the free gas depth within the sediment. Thus, water depth, FGD and 0.8 bar N2 pressure (see above) were taken into account for each

calculation of methane concentration at the FGD. The [CH4]FGD also depends on the temperature (±2.6% solubility change by ±1°C temperature variation) and the salinity (±0.7% solubility change by ±1 salinity variation) (Yamamoto et al., 1976).

3. RESULTS

Figure 2 shows three examples of sulfate and methane profiles in different sediment cores. In the periphery of the Basin the shallow, 4.3-m thick Littorina Sea mud caused a sulfate minimum of ca 30 µM but no accumulation of methane (≤6 µM) at any depth in the 8-m long gravity core (Fig. 2A). In a core from the methane-rich sediment without gas bubbles (Fig. 2B), steep gradients of sulfate and methane defined an SMT at 0.35 m depth in the HML. Methane occurred throughout the Littorina Sea mud and penetrated down into the layers deposited >8,500 yr BP: the Ancylus Lake clay, the Yoldia Sea clay and the deep Baltic Ice Lake clay. The Baltic Ice Lake clay had low methane and still contained a small relic sulfate pool from which sulfate diffused upwards and generated a secondary sulfate-methane transition at 8-9 m depth. The third example is a core from the central basin with very shallow free gas (Fig. 2C). The linear gradients fitted to the sulfate and methane data define an SMT at 0.18 m depth. The linear extrapolation of methane concentrations indicates that the free gas depth (FGD) may be reached at 0.60 m depth, below which gas bubbles occurred. According to the acoustic data for this site, gas bubbles started at 0.51 m depth. The difference of 9 cm is within the resolution of the acoustic technique (see Section 2.2.2.).

The methane diffusion fluxes up towards the sulfate-methane transition (SMT) were calculated for each coring site. The magnitude of these fluxes are indicated on the map in Fig. 1D. The highest calculated flux was 3.7 mmol CH4 m-2 day-1.

3.1. Shallow-seismic transect

A representative shallow-seismic transect was selected to illustrate the relationship between sedimentary units, as determined by the acoustic survey, and methane distribution and flux, as determined by the analysis of sediment cores. Figure 3 shows data from the south-to-north oriented transect along the red line in Fig. 1B and 1C. The transect highlights the Postglacial and Holocene stratigraphy of the Bornholm Basin, which has been interpreted from seismo-acoustic transects and sediment core descriptions (cf. Moros et al., 2002; Jensen et al., 2017). The underlying Glacial to

Lateglacial deposits, studied by the IODP Leg 347 drilling expedition in the northern periphery of the Basin, were comprised of glacial diamictic sand and gravel overlaid by a 25-m thick deposit of light-brown, organic-poor lacustrine clay from the Baltic Ice Lake (Andrén et al., 2015b). The shallow-seismic transect in Fig. 3 identifies the uppermost part of this Ice Lake clay, which is overlain by thin clay deposits from the brackish Yoldia Sea and the limnic Ancylus Lake. Organicrich Littorina Sea mud from the modern, brackish-marine Baltic Sea covers these deposits with a very uniform surface depth contour throughout the basin. Thus, the bathymetry varies by only 3 m over the total distance of 33 km between the eight coring positions (Fig. 3A). The Holocene mud unit varies strongly in thickness, however, depending on the topography of the underlying formations and the spatial and temporal sedimentation patterns in the basin. In areas where methane bubbles caused acoustic blanking in the Holocene mud the stratigraphy of the underlying strata was recorded by lowfrequency seismic techniques (see Fig. 7 in Jensen et al., 2016).

The seismic transect in Fig. 3A shows a characteristic feature of the Bornholm Basin: the marine Holocene organic-rich mud sediments form a wedge shaped sediment body, thickest towards the southwest (Jensen and Endler, 2012). This is due to Lateglacial sedimentation, followed by contour current sedimentation dominated by inflowing currents during the marine Littorina Sea stage transgression. Figure 3A indicates an inflow channel to the south with a water depth down to 100 m and indications of sediment erosion on the flanks. On the southern flank of the channel there was no Littorina Sea sediment but to the north of the inflow there is net sedimentation and up to 10-m thick Littorina Sea mud. Where the thickness of the Littorina Sea mud exceeds 8 m, methane gas bubbles are trapped in the mud and cause blanking of the acoustic signal (labeled “Gas” in Fig. 3A). Core data from the two southern-most coring positions in Fig. 3A, shown to the left in Fig. 3B, are the same core data as presented in Fig. 2B and 2C.

Sediment coring was done at eight positions along the transect, as marked in Fig. 3A, and the depth distribution of methane was analyzed in each core (Fig. 3B). In the gassy sediment, the SMT occurred within a few dm from the sediment surface, gas bubbles reached up to about 0.5 mbsf, and methane fluxes peaked with 2.5 mmol CH4 m-2 d-1 (Fig. 2C). In the deep HML outside the gassy area methane also reached up to within several dm from the sediment surface (Fig. 2B and 3B). As the HML was thinning towards the north, the methane concentrations decreased and the SMT depth increased. At the northernmost station of the transect, where the HML thickness was only 4 m, methane did not exceed 10 µM concentration at any depth in the 9.5-m deep core and the methane flux was insignificant (Fig. 3B, far right panel).

Sampling and analysis of methane in the cores were performed fast, consistently and with great care. Yet, the data in Fig. 3B show distinct scatter at high methane concentrations for three of the cores. This is due to high super-saturation and loss of methane by the retrieval of the cores to 1 bar ambient pressure on the deck. At the mean in situ temperature of 7°C and salinity of 17 (Rak and Wieczorek, 2012), the solubility of methane at 1 bar pressure is 1.87 mM (Yamamoto et al, 1976). The data used to calculate diffusion gradients near the top of the methane zone show smooth depth trends, indicating no significant methane loss, even up to 5-6 mM concentration in some cases, corresponding to 3-4 bar of methane partial pressure or 40% of the methane concentration at the FGD. We have visually distinguished the flawed methane data in Fig. 2B and 2C by shaded symbols.

3.2. Models of methane fluxes

3.2.1. Holocene mud layer Figure 4A shows the calculated methane fluxes in individual sediment cores (“Core CH4 flux”) plotted against the HML thickness. Methane concentrations remained at a background level of <10 µM at HML thicknesses less than 4 m. In those cores, sulfate penetrated throughout the HML and into the underlying Postglacial clay. The example shown in Fig. 2A with 4.3 m HML thickness is a boundary case where sulfate dropped to a minimum of 30 µM, but there was still no methane accumulating zone. Cores with dissolved methane but no gas are clustered among HML thicknesses of 4.5-8 m (Fig. 4A). In sediments with HML thickness >8 m, the in situ methane concentrations exceeded the ambient hydrostatic pressure of 8-10 bar (75-95 m water depth) and gas bubbles were present, starting from different depths beneath the sediment surface.

We calculated a linear regression for sediments with HML thicker than the 4 m threshold (Fig. 4A) to derive the statistical HML model:

(8) where the methane flux is measured in mmol CH4 m-2 d-1 and the HML thickness in m. We used Eq. (8) for the mapping of methane fluxes by the HML model.

3.2.2. Free gas depth

According to the model in Eq. (7), the methane flux in gassy sediments should be inversely proportional to the free gas depth (FGD), provided that other parameters are constant. In Fig. 4B the modeled methane fluxes are plotted against 1/FGD of each site. The best fit to the data is:

(9)

Thus, there is a close linear correlation between the methane fluxes calculated from the FGD model for each site and the inverse FGD at that site. This linear correlation provides the best estimate of JCH4 from the acoustically determined free gas depths. We therefore used Eq. (9) for the mapping of methane fluxes based on the acoustically mapped FGD.

Figure 4C shows how the core CH4 fluxes calculated from the methane data measured in each core relate to the methane fluxes calculated from the FGD model for the same site. The core CH4 fluxes were calculated from Eq. 3 using a methane concentration gradient chosen as the best linear fit to the steepest slope of the CH4 concentrations just beneath the SMT. The linear regression through the origin is:

(10)

As seen in Fig. 4C, the two independent estimates of CH4 fluxes are close to a 1:1 ratio. Thus, considering the approximations made and other sources of uncertainty, which are discussed in section 4.2, the FGD model matches systematically the methane fluxes calculated directly from core data at individual coring sites.

The FGD model in Eq. (7) expresses the methane fluxes as a function of the acoustically determined FGD. The equation can be rewritten, however, to express instead for each site the modeled FGD as a function of the core methane fluxes (Supplementary Information, Eq. (S5)). In Fig. 4D this modeled FGD for individual sediment cores is related to the acoustically measured FGD at the same sites. The two independently determined sets of FGD are presented in Fig. 4D. The data are plotted as 1/FGD because a) this is the relation between the methane fluxes and b) this provides a more informative spread of the data. The linear regression through the origin is:

(11)

As seen from Fig. 4D and Eq. (11), the two independent estimates of the FGD based on modeling and on acoustic measurement show a 1:1 ratio. This supports that the FGD model and its underlying assumptions yield a realistic estimate of the FGD.

3.3

Mapping of methane fluxes

Figure 5A shows the thickness of the Holocene mud layer for the mapped area of the Bornholm Basin. In the periphery of the basin, at water depths of <80 m, the HML is generally thin, <2 m. At many places around the basin there has even been no net deposition of Holocene mud and the underlying Postglacial clay penetrates up to the sediment surface (Jensen et al., 2017). In the deepest part of the basin, including several sub-basins to the NW, the HML is deep due to valleys in the underlying glacial topography (Fig. 3A) and the HML reaches a maximum thickness of 20 m.

The distribution of methane fluxes, calculated from the HML model according to Eq. (8), is shown in Fig. 5B. The flux distribution obviously mirrors the HML thickness map in Fig. 5A. In the periphery of the basin, where the HML thickness is <4 m, methane does not accumulate (“low methane”) and the calculated flux is insignificant. Towards the central part of the basin, where the thickest HML deposits occur, the JCH4 increases steeply and reaches highest values of 4-7 mmol CH4 m-2 d-1. This even exceeds the highest JCH4 calculated in individual cores (Fig. 4) and results from an extrapolation of the model to the greatest HML thicknesses where coring was not done.

The geographic distribution of free gas, determined by acoustic mapping throughout the basin (Fig. 1B), is shown in Fig. 5C. Shallow gas occurred only in areas where the HML thickness exceeded 4 m (Fig. 4A). In the periphery of the gassy area the bubble front was deep, >5 mbsf, while in areas with the thickest HML the gas could be detected even within 0.5 mbsf. Fig. 2C shows an example of the sulfate-methane distribution at such a hotspot.

Figure 5D shows the distribution of methane fluxes, calculated from the FGD model using Eq. (9). The distribution is only partly similar to that calculated from the HML model (Fig. 5B). One main difference is that the FGD-modeled methane flux is zero in the periphery of the area where the HML model calculated low but positive fluxes. This is because the areas with HML thicknesses between

4.5 and 7-8 m had dissolved methane but no shallow gas and were therefore not included in the FGD model (cf. Fig. 4A).

This difference in methane fluxes calculated by the HML and FGD models is mapped in Fig. 5E (see Discussion section 4.2.). In most of the gassy area the two models calculate similar methane fluxes (purple signature in Fig. 5E). However, the FGD model calculates higher methane fluxes than the HML model at the central hotspots of methane flux. Only in a restricted area with up to 20 m HML thickness towards the west does the HML model calculate the higher fluxes.

3.3. Methane flux budget

The total area included in this analysis is shown by blue color contrast in Fig. 5 and comprises 4,261 km2. As seen from Fig. 5, the area covers most of the deep basin where methane is being produced. However, it misses the NW corner where the maps show a truncated area of thick HML and free gas but where seismo-acoustic data were not available. Out of the total gridded area, 2,857 km2 or 67% had HML thicknesses of 0-4 m and very low methane, <10 µM (Table 1). The other 1,404 km2 or 33% had methane with low to high areal fluxes. With the HML thickness increasing from 4 to >20 m, the area of each 2-m depth interval generally decreased. The main methane contribution came from the 8% of the area that had HML thicknesses >12 m and very high methane fluxes, thereby accounting for 56% of the methane flux in the Bornholm Basin. When integrated over the entire gridded area, the total methane flux calculated by the HML model was 14.4 × 105 mol CH4 d-1 or 17.3 metric ton C d-1.

By the alternative approach based on the free gas depth the highest areal fluxes were focused towards the shallowest FGD of <0.5 m and fluxes decreased with increasing FGD. When multiplied with the area of each FGD depth interval the main contributions to the methane flux, 67%, were concentrated to the shallowest FGDs of <1.0 m (Table 1). The total area with free gas was 943 km2 and corresponds to the area with HML thicknesses >8 m. This is two thirds of the area that had methane in the sediment. When integrated over the entire basin, the total methane flux calculated by the FGD model corresponds to 12.5 × 105 mol CH4 d-1 or 15.0 metric ton C d-1. Thus, the two estimates are of similar magnitude but the FGD model yields 13% lower total flux of methane than the HML model.

4. DISCUSSION

It was a striking observation that the methane distribution in the Bornholm Basin sediments was geographically extremely heterogeneous and was not related to bathymetry or sediment type in any simple manner. In fact, without detailed geophysical mapping of the sub-seafloor stratigraphy it would not have been possible to interpret the occurrences of methane observed by sediment coring. In spite of the very even and gently sloping bottom topography of the Bornholm Basin beneath the 70 m depth curve, some of the retrieved cores were highly super-saturated with methane while others had only low µM levels of methane (Fig. 3B). Acoustic mapping of the underlying Cretaceous bedrock, Glacial till and Postglacial limnic clay showed that the thickest Holocene mud layer occurred towards the southwest. A comparison of the methane fluxes and the thickness of the HML throughout the Bornholm Basin showed that an HML thickness of >4 m was a prerequisite for effective sulfate depletion and methane accumulation (Fig. 4A).

This depth threshold is related to the salinity and productivity history of the Baltic Sea basin over the past 12,000 years combined with the gradual deposition of Lateglacial and Postglacial sediment (Mogollon et al., 2012; Holmkvist et al., 2014). Thick deposits of silt and clay from glacial runoff were laid down in the Baltic Ice Lake as the glaciers retreated more than 12,000 years ago. After the first intrusion of seawater about 11,600 yr BP and over the following several thousand years, seawater ions such as sulfate diffused deep down into the organic-poor limnic clay in which no significant sulfate reduction took place due to its low organic carbon content. As the brackish Littorina Sea phase started, the mud deposits were initially thin and not yet effective barriers against continued deep penetration of sulfate in spite of their relatively high organic carbon content. During the mid-Holocene, about 7,000-5,000 yr BP, the salinity was also higher than today and the sulfate concentration may have reached 1.7-fold the present-day level in the southern Baltic (Gustafsson and Westman, 2002; Kortekaas, 2007).

Only as the HML grew to a sufficient thickness of >0.5 m did partial sulfate depletion turn around the gradient and sulfate started to diffuse back up from the underlying clay towards a sulfate minimum in the marine mud (Holmkvist et al., 2014). Such a transient situation is depicted in Fig. 2A where sulfate was nearly depleted but methane accumulation had not yet started in the HML. A further stage is shown in Fig. 2B where methane had a broad maximum in the HML but dropped to a low level in the glacio-lacustrine clay where sulfate was still present. Coring at positions with different HML thicknesses in the Bornholm Basin revealed different stages of this transition, from

complete sulfate penetration with only an intermediate sulfate minimum in the upper HML to complete sulfate depletion with high concentrations of methane and even free methane gas.

A similar trend of sulfate depletion and methane production that developed as Holocene mud accumulated was observed in the neighboring, 40-50-m deep Arkona Basin (Fig. 1A) (Thiessen et al., 2006; Mogollon et al., 2012). By a time-resolved reaction-transport model, Mogollon et al. (2012) simulated how the build-up of a >2-m deep HML over 2,500 years after the onset of the Littorina Sea was required for a sulfate minimum to reach depletion and for methane production to start. Another thousand years was required for depletion of the deep relic sulfate reservoir and yet another few thousand years for free methane gas bubbles to developed within the deeper part of the HML.

These Holocene time slices in the Mogollon et al. (2012) model show sulfate and methane zonations that are very similar to the diverse zonations observed at different locations in the Bornholm Basin. Yet, while the Arkona Basin model assumed a fixed sedimentation rate and used time as a variable, our observations in the Bornholm Basin have a fixed time (8,600 years since brackish-marine mud deposition started) and variable sedimentation rates. Since modern sedimentation rates in the Bornholm Basin, determined by multiple positions of 210Pb dating, are consistent with the HML thicknesses, sulfate and methane in the Bornholm Basin sediments will in the future expectedly develop along a similar trajectory as modeled for the Arkona Basin, yet at a slower pace towards the periphery of the basin where sedimentation rates are low.

4.1. The HML model of methane flux

A plot of calculated methane fluxes relative to the HML thickness showed a positive linear correlation for HML >4 m with an r2 of 0.69 (Fig. 4A; Eq. (8)). The level of methanogenesis was highly sensitive to the HML thickness. For example, at Station BB01 the HML was only 6% thicker than at BB02 (Fig. 3A), yet the methane flux was seven-fold higher (Table S1). The seismic reflectors pinched out somewhat in the upper part of the HML, which may explain that there was a lateral decrease in recent sedimentation rate at BB02 and therefore lower rates of organic matter mineralization. Yet, the observation also reflects a positive feed-back mechanism described recently by Flury et al. (2016). The feed-back is due to the steep decrease in organic matter degradability with depth and age in the sediment. Thus, Beulig et al. (2018) found in Bornholm Basin sediments a power law relationship between organic carbon oxidation rate and age (in years):

(12)

Importantly, this relationship was valid throughout the sulfate zone and into the methane zone. A small increase in methanogenesis at depth will thus push the SMT a bit upwards. Thereby, sediment layers with much more readily degradable organic matter become available for methanogenesis, which in turn leads to a higher overall methane production and flux and thereby pushes the SMT upwards even further. Although this is clearly not a run-away feed-back, it is indeed a mechanism that enhances the geographic heterogeneity in methane production and distribution.

While a minimum threshold thickness of the HML of 4 m was required for the sediment to have sulfate depletion and accumulate methane (Fig. 4A), a thicker HML of >6 m was required for sulfate depletion to also reach the clay deposits beneath the HML (Fig. 2B). The secondary SMT and the resulting deep sink for methane therefore constituted an additional variable in the relationship between methane flux and HML thickness. Furthermore, the deposition rate and distribution of Holocene mud has changed somewhat over the past 8,600 years depending on the flow patterns of bottom water currents carrying suspended material into the basin (Jensen et al., 2017). This material enters mainly from the northwest where flow velocities are higher and where the bottom flow has created a trough along the southwestern rim of the basin (Fig. 3). Along the flow path throughout the basin the settling material may gradually become finer and richer in organic matter (Christoffersen et al., 2007; Leipe et al., 2010). The same HML thickness will therefore lead to higher methane accumulation in the central and eastern part of the basin than in the western and the HML model may tend to overestimate methane fluxes in the western part.

We investigated whether this trend was significant by comparing the geographic distribution of methane fluxes calculated by the two models in Fig. 5B and D. A map of the difference in methane fluxes calculated from the HML model minus fluxes calculated from the FGD model shows that the two models compare well throughout most of the study area but that discrepancies exist in spots in the western part of the basin where the HML model calculates higher fluxes than the FGD model (orange color in Fig. 5E). This is pronounced in limited areas where deeply buried glacial troughs had high sedimentation during the early Littorina Sea period. The HML here even reaches down to 40 mbsf, but the organic matter is relatively old throughout the lower sediment column and therefore has lower rates of methanogenesis than predicted by the extrapolated HML model. In the large central basin further to the southeast, the HML model sometimes underestimates the methane fluxes

calculated by the FGD model, especially where the bubble front comes very close to the sediment surface and reveals hotspots of methane flux not captured by the HML model (Fig. 5D; light-blue color in Fig. 5E). This comparison shows the importance of using two independent models to obtain accurate maps the methane fluxes.

4.2. The FGD model of methane flux

The highly simplified FGD model expressed in Eq. (7) is based on a small number of parameters, each of which is straightforward to measure and map: a) The sulfate concentration of the bottom water, a simple function of salinity. b) The methane concentration at the FGD, a function of water plus sediment depth, salinity and temperature. c) The whole sediment diffusion coefficients of sulfate and methane, a function of temperature, salinity and tortuosity, as determined by the sediment porosity. d) The FGD which is determined from sediment echo sounding with appropriate acoustic equipment.

We tested the FGD model by a comparison of the model results with methane fluxes calculated directly from 17 cores taken at different positions within the area of gassy sediments. As seen from the linear regression in Fig. 4C (Eq. (10)), methane fluxes calculated from core data and from the FGD model yielded comparable results. The linear correlation through origin showed a slope of 1.02, i. e. close to unity, and an r2 of 0.67 (Fig. 4C). This supports the basic assumptions behind the FGD model in Eq. (7), but it also reveals significant scatter in the correlation.

One potential variable and source of scatter, which was not included in the FGD model, is the ratio of sulfate and methane fluxes to the SMT. We assumed a 1:1 ratio according to the stoichiometry of AOM. However, a global compilation of sulfate and methane flux data to the SMT shows that the flux ratio may vary and is in average: JSO4:JCH4 = 1:0.7 (Egger et al., 2018). The excess sulfate flux is due to organoclastic sulfate reduction within the SMT (Flury et al., 2016; Beulig et al., 2018). If we therefore assume that only 70% of the sulfate flux to the SMT is used for AOM to balance the methane flux, then the FGD model in Eq. (7) modifies to:

(13)

All other parameters in the equation are unaltered by this modification. Those parameters vary only little due to the small variations in bottom water temperature (which affects mostly DSO4 and DCH4), salinity (which affects mostly [SO42-]SW) and water depth (which affects mostly [CH4]FGD). If we take 30% organoclastic SRR in the SMT into account, then the methane fluxes calculated from the FGD model will in average be reduced by 8%.

The scatter in the correlation between gradient-calculated and FGD-modeled methane fluxes (Fig. 4C) is due to simplifications in the model that are fulfilled only to certain degree. Scatter may be due to temporal variations in bottom water salinity and temperature or limited accuracy in the acoustic determination of the free gas depth. Thus, the sediment echosounder resolution of the FGD determination near the sediment surface was at best 0.2 m. The model furthermore assumes only diffusional transport and no bioturbation or bioirrigation, i.e. no importance of advective transport of sediment or pore water due to the activity of burrowing macrofauna. While this is a reasonable assumption when the FGD is positioned deeply, >1 mbsf, it becomes increasingly questionable as the FGD moves close to the sediment surface. Bioturbation and bioirrigation will tend to compress the sulfate and methane gradients by a given FGD and thereby cause the model to underestimate the actual methane flux. The shallowest FGD was 0.5 m, yet there was no trend towards a systematic model underestimation of the methane flux for the shallow FGD (high methane flux in Fig. 4C). The reason may be that macrofauna was not abundant or active in the deeper Bornholm Basin during our sampling expeditions in summer due to suboxic or anoxic conditions (cf. Leppakoski 1969; Conley et al., 2009).

The FGD model assumes linear gradients of sulfate from the sediment surface to near the SMT and of methane from near the SMT down to the FGD. Sulfate gradients are often convex down due to sulfate reduction, but just as often they are essentially linear due to bioirrigation, to decreasing porosity with increasing depth, to non-steady state sedimentation, or to sulfide re-oxidation by buried iron(III) (Fig. 2B,C) (e.g., Schulz et al., 1994; Niewöhner et al., 1998; Flury et al., 2016; Riedinger et al., 2017). Similarly, the methane gradients are expectedly concave up due to methane production, although this is most often not possible to observe. Due to supersaturation and methane loss it is generally only the top of the methane gradient that can be analyzed accurately. The deeper gradient was simply extrapolated linearly in our model.

Convex downward sulfate profiles due to sulfate consumption would mean lower sulfate flux to the SMT than assumed in the linear FGD model. As a result, the SMT would be located at a shallower depth relative to the FGD than assumed by the FGD model in order to balance the sulfate and methane fluxes. For a similar reason, concave upward methane profiles also mean that the SMT is located at a shallower depth than calculated from a linear methane profile by the FGD model. In Fig. 6 we have plotted the methane fluxes relative to the SMT depth (Table S1), both based on methane data measured in retrieved sediment cores. It should thereby be noted that we did not have sulfate data for all the cores. To be consequent, we therefore used as the SMT in Fig. 6 the depth where the linear methane gradients hit zero, instead of using the depth of equimolar sulfate and methane concentration. As seen by the examples in Fig. 2B and 2C, this had insignificant effect in the cores with the shallowest SMT where a small shift in estimated SMT depth could otherwise have the largest relative effect.

In Fig. 6 we have also plotted the linear function expressed in Eq. (5), which describes the theoretical relationship between methane flux (JCH4) and SMT depth (ZSO4), assuming linear gradients and a 1:1 stoichiometry. The function is seen to mostly overestimate the diffusive methane fluxes. If a 1:0.7 stoichiometry of sulfate vs. methane flux is taken into account the function moves into the data swarm (Fig. 6). If the shallowing of the SMT due to convex-down sulfate profiles and concave-up methane profiles is also taken into account, the function will expectedly provide a good description of the data. However, to show this would require a more complex transport-reaction model of organic matter degradation calculated for each site, which is beyond the objective of this study.

Such an analysis was done by Dale et al. (2009), however, who used a transport-reaction model with two organic carbon pools of different reactivity and with assumed bioirrigation to model sulfate reduction and methanogenesis in gassy sediments of Aarhus Bay in relation to the FGD. The modeled methane fluxes to the SMT were compared to the methane fluxes calculated from methane measured in sediment cores from 19 stations in a gassy area. The measured fluxes were in average 88% of the modeled fluxes with an r2 of 0.77. The model derived an inverse relationship between the methane flux and the FGD (JCH4 = a × FGD-1.05), which is very similar to the relation shown in Fig. 4B (Eq. (9)). This comparison shows that the simpler, linear FGD model in our study provided a similarly good fit to the calculated methane fluxes although our model assumptions are only generalized simplifications.

4.3. Methane flux budget Within the studied area of the Bornholm Basin, we estimate that 2,857 km2 or 67% has insignificant dissolved methane concentrations (<10 µM) throughout the Holocene mud layer (Table 1). In the periphery of the area with methane-enriched sediments the HML thickness was 4-8 m. Methane fluxes were here <1 mmol m-2 d-1 (Fig. 5B) and accounted for only a fifth of the entire methane flux in the basin. The area with free gas in the sediment, 943 km2, was largely identical to the area with HML thickness >8 m. For comparison, the area with free gas in the Arkona Basin to the west of Bornholm, analyzed by similar seismic surveys, yet at lower spatial resolution, is fifty percent larger, 1,500 km2 (Thiessen et al., 2006). As the Arkona Basin is shallower than the Bornholm basin, 40-50 m, bubble formation starts there at a 35-50% lower methane concentration.

The total subsurface methane fluxes in the Bornholm Basin, estimated from the HML and the FGD models, are 14.4 × 105 and 12.5 × 105 mol CH4 d-1 or 17.3 and 15.0 metric ton C d-1, respectively. Thus, the FGD-modeled methane flux for the basin is 13% lower than the HML-modeled flux. Two factors may lead to this difference: a) the methane flux in the non-gassy sediments with an HML thickness of <8 m are not included in the FGD budget (Fig. 5, Table 1). This means that the FGD model misses 19% of the methane flux included in the HML model. These 19% should be added to the FGD budget to cover the entire area; b) organoclastic sulfate reduction in the SMT is not taken into account by the assumption of a 1:1 flux ratio of sulfate and methane. If a flux ratio of 1:0.7 is used instead, then the FGD model fluxes will decrease by 8%. The revised methane flux budget calculated from the FGD model thereby corresponds within a few percent to the HML budget, which confirms the agreement between the two models. The methane fluxes calculated from Bornholm Basin sediment cores were 0.13-0.31 mmol m-2 d-1 in methane-accumulating but non-gassy sediments and 0.57-3.70 mmol m-2 d-1 in gassy sediments. These fluxes are in a similar range as determined for other Baltic Sea regions, e.g. 0.12-0.43 mmol m-2 d-1 in the Arkona Basin (Mogollon et al. (2012), 0.15-0.40 mmol m-2 d-1 in the Gotland Deep (Piker et al., 1998), 0.03-0.13 mmol m-2 d-1 in the Eckernförde Bay (Whiticar, 2002), and 0.05-0.15 mmol m-2 d-1 in the Aarhus Bay (Flury et al., 2016).

Egger et al. (2018) developed a global seafloor map of methane fluxes to the SMT with 10×10 km resolution using algorithms based on sedimentation rate, water depth, distance to coast, and seasurface particulate organic carbon. The algorithms were calibrated with methane data from 740

marine sediment cores, including 58 cores from the Baltic Sea. The authors calculated global mean methane fluxes for the inner continental shelf of 0.86 mmol m-2 d-1 at 0-10 m water depth and 0.27 mmol m-2 d-1 at 10-50 m water depth. For the entire Bornholm Basin they calculated a methane flux of 12 × 105 mol CH4 d-1 (14 metric ton C d-1), i.e. very similar to our estimate. Although the study by Egger et al. (2018) also included 12 cores from the Bornholm Basin, this was not the reason for the good correspondence but rather that the global algorithms applied were also representative for this basin.

4.4. Free gas in the sediment

Accumulation of free gas may lead to upwards migration of the buoyant gas bubbles. In a model study of Arkona Basin sediments (Mogollon et al., 2012), buoyant interstitial gas movement was assumed not to take place, which led to an extremely high modeled gas content of up to 8% by volume. In a rather similar sediment from Aarhus Bay, Dale et al. (2009) concluded from modelling that gas migration did take place. Flury et al. (2016), on the other hand, found that the calculated methane flux in Aarhus Bay could be explained by molecular diffusion without gas migration. We conclude that it is difficult from modeling alone to determine whether gas bubbles rise in the gassy sediment and then dissolve again as they reach the FGD. However, our models are independent of whether or not methane bubbles are rising up through the gassy sediment beneath the FGD. The FGD model only assumes that gas bubbles at the FGD are in equilibrium with the partial pressure of methane in the ambient pore water at that depth. The observation that FGD-modeled methane fluxes are higher in the hotspots than the HML-modeled fluxes (Fig. 5E) is therefore not related to bubble migration.

High production rates and accumulation of methane may lead to ebullition of gas across the sediment surface (e.g., Kipphut and Martens, 1982; Martínez-Carreño and Garcia-Gil, 2013). Gas escape structures such as pockmarks were not observed in the Bornholm Basin during acoustic mapping with swath bathymetry or with sediment echosounder. Acoustic plumes as indication of rising gas bubbles in the water column were also not observed during our five expeditions. Elevated methane concentrations of 10-50 nM have been measured in the Bornholm Basin with highest values at midwater depths of 50-70 m rather than in the bottom water (Gülzow et al., 2014). The methane apparently originated from bottom water of the shallower Arkona Basin, which enters the Bornholm Basin by frequent intrusions, rather than from the bottom of the Bornholm Basin. We have therefore

not included advective gas flux across the sediment-water interface in the Bornholm Basin methane budget.

Although the dissolved methane in Bornholm Basin sediments clearly originates from production in the Holocene mud (e.g. Fig. 2B and 3B) it is a question whether the gas bubbles might have a deeper source. Tóth et al. (2014b, 2015) made a detailed bubble analysis in a transect through a gassy area just northeast of Bornholm where the FGD ranged from 2 to 7 mbsf (red arrow in Fig. 1C). By the use of multiple acoustic frequency bands they determined a resonance peak at 3-4 kHz corresponding to a gas bubble radius of 2-3 mm, whereby most bubbles were probably smaller. Based on a compressible fluid model the authors concluded that the gas was confined to the Holocene mud and was not present in the underlying Lateglacial and Postglacial deposits. This shows that the methane is exclusively Holocene of age and biogenic of origin and does not have an older and deeper source. The underlying Lateglacial to Postglacial clay is rather a sink for methane, in particular where the HML thickness is <6 m and relic sulfate still remains in the limnic clay (Fig. 2B and 3B).

Tóth et al. (2015) calculated that the average free gas content in the studied mud layer was only 0.046% by volume at in situ pressure, and that the total, depth-integrated gas volume corresponded to 2.8 mol CH4 m-2. This is only 2% of all the methane dissolved in the pore water in the sediment column according to our data. As the mean methane flux in the same gassy area was 0.5-2 mmol m-2 d-1 (Fig. 5B and D), the methane reservoir contained in gas bubbles corresponded to 4-15 years of methane production. Since the Holocene mud at 2-7 m depth is thousands of years old and has very low rate of methanogenesis (Beulig et al., 2018), this suggests that the gas bubbles constitute an insignificant and non-dynamic pool in the overall methane budget.

5. Conclusions

Half of the global flux of particulate organic carbon to the seafloor takes place in coastal waters and on the inner shelf at water depths of less than 50 m (Dunne et al., 2007). This is also where nearly two thirds of the global methane production in marine sediments takes place (Egger et al., 2018). These coastal regions represent the marine ecosystems that are most exposed to eutrophication, oxygen depletion and global warming (Diaz and Rosenberg, 2008; Middelburg and Levin, 2009). Each of these factors may enhance organic mineralization and methane production rates in the

seabed (Best et al., 2006). It is therefore critical to understand the magnitude of methane production in marine sediments and the present and future controls on the methane fluxes.

By a combination of seismo-acoustic transects and methane analyses in multiple sediment cores we mapped methane fluxes in an entire sediment basin and quantified the total methane production in the seabed. Two independent approaches to model these methane fluxes yielded similar total flux data and rather similar geographic distributions of the fluxes. One approach is based on acoustic tracking of shallow gas. The other is seismo-acoustic tracking of Holocene sediment deposits, an approach which, to our knowledge, is here used for the first time to quantify and map the flux of methane to the SMT. The results show that remote detection of shallow gas and sediment stratigraphy by seismo-acoustic methods provides information that is indispensable for the interpretation and extrapolation of sediment core data. Thus, even within an apparently uniform sediment basin the sediment core data demonstrated extreme spatial heterogeneity in methane production, ranging from essentially zero methane to high methane and to very shallow methane gas bubbles. In spite of the high number of coring sites, it would therefore not have possible to develop an accurate budget of the total methane flux in the Bornholm Basin based only on sediment core data.

The Baltic Sea is an intracontinental sea at rather high latitude (54-66°N) where the last GlacialInterglacial transition left a distinct sediment stratigraphy that could be used for our HML modeling of methane fluxes. Along all global ocean margins, however, continental shelves were dry land for tens of thousands of years around the Last Glacial Maximum, where the global eustatic sea-level reached a low-stand 120 m below the present sea-level (e.g., Fleming et al., 1998). During this period weathering and erosion shaped distinct Glacial landscapes on the inner shelf. In many protected shelf areas Holocene marine sedimentation has since then covered the Glacial landscapes with a thick mud blanket. Such Holocene mud is a major source of marine biogenic methane today and is critical for the future development of methane dynamics in ocean margin sediments. We believe that an HML model approach, as developed in this study and possibly combined with an FGD model approach, will have wide applicability for the mapping and quantification of methane production in coastal and continental shelf sediments.

Acknowledgements

We thank Jeanette Pedersen and Karina Bomholt Oest for assistance with the analytical work and André Pellerin for graphical assistance. We thank all colleagues and crew who organized and participated in the five research expeditions and assisted by the sampling and analyses of sediment core material. During the research projects, the seismic data acquisition was made in a close cooperation with R. Endler, The Leibniz Institute for Baltic Sea Research (IOW), whom we also thank for providing additional archive data. We thank the four reviewers for their very constructive and helpful recommendations of how to improve the original manuscript. The expeditions were carried out in association with the research projects, METROL (http://metrol.mpibremen.de/; coordinators B. B. Jørgensen and C. Borowski) and BALTIC GAS (http://www.balticgas.au.dk/; coordinators B. B. Jørgensen and H. Fossing). The study was supported by the two projects, METROL (EU 5th FP, grant # EVK3-CT-2002-00080) and BALTIC GAS (EU 7th FP, BONUS program, grant # 217246), by an ERC Advanced Grant to BBJ (MICROENERGY, EU 7th FP, grant # 294200), by the Danish Center for Marine Research (AUBO16 Cruise), and by the Danish National Research Foundation (DNRF grant # 104).

References

Andrén E., Andrén T., Sohlenius G. (2000) The Holocene history of the south-western Baltic Sea as reflected in a sediment core from the Bornholm Basin. Boreas 29, 233–250. Andrén T., Björck S., Andrén E., Conley D., Zillén L. and Anjar J. (2011) The development of the Baltic Sea Basin during the last 130 ka. In The Baltic Sea Basin (eds. J. Harff, S. Björck and P. Hoth). Springer, Berlin, pp. 75–97. Andrén T., Jørgensen B. B., Cotterill C., Green, S., and Expedition 347 Scientists (2015a) “Baltic sea paleoenvironment,” in Proc. IODP, 347: College Station, TX (Integrated Ocean Drilling Program). Available online at: http://publications.iodp.org/proceedings/347/347title.htm Andrén T., Jørgensen B. B., Cotterill C., Green, S., and Expedition 347 Scientists (2015b) Site M0065. In Andrén, T., Jørgensen, B.B., Cotterill, C., Green, S., and the Expedition 347 Scientists, Proc. IODP, 347: College Station, TX (Integrated Ocean Drilling Program). doi:10.2204/iodp.proc.347.109.2015 Best A. I., Richardson M. D., Boudreau B. P., Judo A. G. Leifer I. et al. (2006) Shallow seabed methane gas could pose coastal hazard. EOS Trans. Amer. Geophys. Union 87, 213-220. Beulig F., Røy H., Glombitza C. and Jørgensen B. B. (2018) Control on rate and pathway of anaerobic organic carbon degradation in the seabed. Proc. Natl. Acad. Sci. (USA) 115, 367-372. Björck S. (1995) A review of the history of the Baltic Sea, 13.0-8.0 ka BP. Quatern. Int. 27, 19-40. Borowski W. S., Paull C. K. and Ussler W. (1996) Marine pore-water sulfate profiles indicate in situ methane flux from underlying gas hydrate. Geology 24, 655-658. Borowski W. S., Paull C. K. and Ussler W. (1999) Global and local variations of interstitial sulfate gradients in deep-water, continental margin sediments: Sensitivity to underlying methane and gas hydrates. Mar. Geol. 159, 131-154. Boudreau B. P. (1996) The diffusive tortuosity of fine-grained unlithified sediments. Geochim. Cosmochim. Acta 16, 3139-3142. Brook E. J., Sowers T. and Orchardo J. (1996) Rapid variations in atmospheric methane concentration during the past 110,000 years. Science 273, 1087-1091. Carstensen J., Andersen J. H., Gustafsson B. G. and Conley D. J. (2014a) Deoxygenation of the Baltic Sea during the last century. Proc. Natl. Acad. Sci. (USA) 111, 5628-5633. Carstensen J., Conley D. J., Bonsdorff E., Gustafsson B. G., Hietanen S. et al. (2014b) Hypoxia in the Baltic Sea: Biogeochemical cycles, benthic fauna, and management. Ambio 43, 26-36. Canfield D. E., Kristensen E. and Thamdrup B. (2005) Aquatic Geomicrobiology. Elsevier, San Diego, California.

Christoffersen P. L., Christiansen C., Jensen J. B., Leipe T. and Hille S. (2007) Depositional conditions and organic matter distribution in the Bornholm Basin, Baltic Sea. Geo-Mar. Lett. 27, 325-338. Conley D. J., Björck S., Bonsdorff E., Carstensen J., Destouni G. et al. (2009) Hypoxia-related processes in the Baltic Sea. Env. Sci. Tech. 43, 3412-3420. Dale A. W., Aguilera D. R., Regnier P., Fossing H., Knab N. J. and Jørgensen B. B. (2008) Seasonal dynamics of the depth and rate of anaerobic oxidation of methane in Aarhus Bay (Denmark) sediments. J. Mar. Res. 66, 127-155. Dale A. W., Regnier P., van Cappellen P., Fossing H., Jensen J. B. and Jørgensen B. B. (2009) Remote quantification of methane fluxes in gassy marine sediments through seismic survey. Geology 37, 235-238. Diaz, R. J. and Rosenberg R. (2008) Spreading dead zones and consequences for marine ecosystems. Science 321, 926-929. Duan Z. H., Møller N., Greenberg J. and Weare J. H. (1992) The prediction of methane solubility in natural waters to high ionic strength from 0 to 250°C and from 0 to 1600 bar. Geochim. Cosmochim. Acta 56, 1451-1460. Dunne J. P., Sarmiento J. L. and Gnanadesikan A. (2016) A synthesis of global particle export from the surface ocean and cycling through the ocean interior and on the seafloor. Global Biogeochem. Cycles 21, #4006. Egger M., Hagens M., Sapart C.J., Dijkstra N., van Helmond N. A. G. M. et al. (2017) Iron oxide reduction in methane-rich deep Baltic Sea sediments. Geochim. Cosmochim. Acta 207, 256-276 Egger M., Riedinger N., Mogollón, J. and Jørgensen B. B. (2018) Global diffusive fluxes of methane in marine sediments. Nature Geosci. 11, 421-425. Fleischer P., Orsi T. H., Richardson M. D. and Anderson A. L. (2001) Distribution of free gas in marine sediments: a global overview. Geo-Mar. Lett. 21, 103-122. Fleming K., Johnston P., Zwartz D., Yokoyama Y., Lambeck K. and Chappell J. (1998) Refining the eustatic sea-level curve since the Last Glacial Maximum using far- and intermediate-field sites. Earth Planet. Sci. Lett. 163, 327-342. Flury S., Røy H., Dale A. W., Fossing H., Tóth Z., Spiess V., Jensen J. B. and Jørgensen B. B. (2016) Controls on subsurface methane fluxes and shallow gas formation in Baltic Sea sediment (Aarhus Bay, Denmark). Geochim. Cosmochim. Acta 188, 297-309. Gülzow W., Gräwe U., Kedzior S., Schmale O. and Rehder G. (2014) Seasonal variation of methane in the water column of Arkona and Bornholm Basin, western Baltic Sea. J. Mar. Syst. 139, 332-347.

Gustafsson B. and Westman P. (2002) On the causes for salinity variations in the Baltic Sea during the last 8500 years. Palaeogeography 17, 1-14. Holmkvist L., Kamyshny Jr. A., Brüchert V., Ferdelman T. and Jørgensen B. B. (2014) Sulfidization of lacustrine glacial clay upon Holocene marine transgression (Arkona Basin, Baltic Sea). Geochim. Cosmochim. Acta 142: 75-94. Iversen N. and Jørgensen B. B. (1993) Diffusion coefficients of sulfate and methane in marine sediments: influence of porosity. Geochim. Cosmochim. Acta 57, 571-578. Jensen J. B., Bennike O., Witkowski A., Lemke W. and Kuijpers A. (1997) The Baltic Ice Lake in the southwestern Baltic: sequence-, chrono- and biostratigraphy, Boreas 26, 217-236. Jensen J. B. and Endler R. (2012) Methane distribution in Holocene marine sediments in the Bornholm Basin, southern Scandinavia. Geol. Surv. Denmark Greenland Bull. 26, 21-24. Jensen J. B., Moros M., Endler R. and IODP Expedition 347 members (2017) The Bornholm Basin, southern Scandinavia: a complex history from Late Cretaceous structural developments to recent sedimentation. Boreas 46, 3-17. doi 10.1111/bor.12194 Jørgensen B. B. and Parkes R. J. (2010) Role of sulfate reduction and methane production by organic carbon degradation in eutrophic fjord sediments (Limfjorden, Denmark). Limnol. Oceanogr. 55, 1338-1352. Kipphut G. W. and Martens C. S. (1982) Biogeochemical cycling in an organic-rich coastal marine basin. 3. Dissolved gas transport in methane-saturated sediments. Geochim. Cosmochim. Acta 46, 2049-2060. Knab N. J., Cragg B. B., Borowski C., Parkes R. J., Pancost R. and Jørgensen B. B. (2008) Anaerobic oxidation of methane (AOM) in marine sediments from the Skagerrak (Denmark): I. Geochemical and microbiological analyses. Geochim. Cosmochim. Acta 72, 2868–2879. Kögler F. C. and Larsen B. (1979) West Bornholm Basin in the Baltic Sea – geological structure and Quaternary sediments. Boreas 8, 1-22. Kortekaas M. (2007) Post-glacial history of sea-level and environmental change in the southern Baltic Sea. Ph. D. thesis, Lund University, 34 pp. Kvenvolden, K. (1988) Methane hydrates and global climate. Global Biogeochem. Cycles 2, 221229. doi.org/10.1029/GB002i003p00221 Lafuerza S., Sultan N., Canals M., Frigola J., Berné S., Jouet G., Galavazi M. and Sierro F. J. (2009) Overpressure within upper continental slope sediments from CPTU data, Gulf of Lion, NW Mediterranean Sea. Int. J. Earth Sci. 98, 751-768. Laier T. and Jensen J. B. (2007) Shallow gas depth-contour map of the Skagerrak-western Baltic Sea region. Geo-Mar. Lett. 27, 127-141.

Leipe, T., Tauber F., Vallius H., Virtasalo J., Uścinowicz S., Kowalski N., Hille S., Lindgren S. and Myllyvirta T. (2010) Particulate organic carbon (POC) in surface sediments of the Baltic Sea. Geo-Mar Lett. 31, 175-188. Leppakoski E. (1969) Transitory return of benthic fauna of Bornholm Basin after extermination by oxygen insufficiency. Cah. Biol. Mar. 10, 163. Martínez-Carreño N. and Garcia-Gil S. (2013) The Holocene gas system of the Ría de Vigo (NW Spain): Factors controlling the location of gas accumulations, seeps and pockmarks. Mar. Geol. 344, 82-100. Meischner D. and Rumohr J. (1974) A light-weight high-momentum gravity corer for subaqueous sediments. Senckenb. Marit. 6, 105-117. Middelburg J. J. (1989) A simple rate model for organic matter decomposition in marine sediments. Geochim. Cosmochim. Acta 53, 1577-1581. Middelburg J. J. and Levin L. A. (2009) Coastal hypoxia and sediment biogeochemistry. Biogeosciences 6, 1273-1293. Mogollón J. M., Dale A. W., Fossing H. and Regnier P. (2012) Timescales for the development of methanogenesis and free gas layers in recently-deposited sediments of Arkona Basin (Baltic Sea). Biogeosciences 9, 1915-1933. Mogollón J. M., Dale A. W., Jensen J. B., Schlueter M. and Regnier P. (2013) A method for the calculation of anaerobic oxidation of methane rates across regional scales: an example from the Belt Seas and The Sound (North Sea-Baltic Sea transition). Geo-Mar. Lett. 33, 299-310. Moros M., Lemke W., Kuijpers A., Endler R., Jensen J. B., Bennike O., and Gingele F. (2002) Regressions and transgressions of the Baltic basin reflected by a new high-resolution deglacial and Postglacial lithostratigraphy for Arkona Basin sediments (western Baltic Sea). Boreas 31, 151-162. Niewöhner C., Hensen C., Kasten S., Zabel M. and Schulz H. D. (1998) Deep sulfate reduction completely mediated by anaerobic methane oxidation in sediments of the upwelling area off Namibia. Geochim. Cosmochim. Acta 62, 455-464. Parkes R. J., Cragg B. A., Banning N., Brock F., Webster G., Fry J. C., Hornibrook E., Pancost R. D., Kelly S., Knab N., Jørgensen B. B., Rinna J. and Weightman A. J. (2007) Biogeochemistry and biodiversity of methane cycling in subsurface marine sediments (Skagerak, Denmark). Environ. Microbiol. 9, 1146-1161. Piker L., Schmaljohann R. and Imhoff J. F. (1998) Dissimilatory sulfate reduction and methane production in Gotland Deep sediments (Baltic Sea) during a transition period from oxic to anoxic bottom water (1993-1996). Aquat. Microb. Ecol. 14, 183-193.

Povinec P. P., du Bois P. B., Kershaw P. J., Nies H., Scotto P. (2003) Temporal and spatial trends in the distribution of Cs-137 in surface waters of Northern European Seas—a record of 40 years of investigations. Deep-Sea Res II 50, 2785–2801 Rak D. and Wieczorek P. (2012) Variability of temperature and salinity over the last decade in selected regions of the southern Baltic Sea. Oceanologia 54, 339-354. Reeburgh W. S. (2007) Oceanic methane biogeochemistry. Chem. Rev. 107, 486-513. Regnier P., Dale A., Arndt S., LaRowe D., Mogollón J. and Cappellen P. V. (2011) Quantitative analysis of anaerobic oxidation of methane (AOM) in marine sediments: A modeling perspective. Earth-Sci. Rev. 106, 105-130. Riedinger N., Brunner B., Krastel S., Arnold G. L., Wehrmann L. M. et al. (2017) Sulfur cycling in an iron oxide-dominated, dynamic marine depositional system: The Argentine Continental Margin. Front. Earth Sci. 5, #33. Schulz H. D. (2006) Quantification of early diagenesis: dissolved constituents in marine pore water. In Marine Geochemistry (eds. H. D. Schulz and M. Zabel), second ed. Springer, Berlin, pp. 75–124. Schulz H. D., Dahmke A., Schinzel U., Wallmann K. and Zabel M. (1994) Early diagenetic processes, fluxes, and reaction rates in sediments of the South Atlantic. Geochim. Cosmochim. Acta 58, 2041-2060. Schneider von Deimling J., Weinrebe W. Tóth Zs., Fossing H., Endler R., Rehder G. and Spiess V. (2013) A low frequency multibeam assessment: Spatial mapping of shallow gas by enhanced penetration and angular response anomaly. Mar. Petrol. Geol. 44, 217-222. Shakhova N., Semiletov I., Salyuk A., Yusupov V., Kosmach D. and Gustafsson Ö. (2010) Extensive methane venting to the atmosphere from sediments of the East Siberian Arctic shelf. Science 327, 1246-1250. Thiessen O., Schmidt M., Theilen F., Schmitt M. and Klein G. (2006) Methane formation and distribution of acoustic turbidity in organic-rich surface sediments in the Arkona Basin, Baltic Sea. Cont. Shelf Res. 26, 2469-2483. Tóth Zs., Spiess V. and Jensen J.B. (2014a) Seismo-acoustic signatures of shallow free gas in the Bornholm Basin, Baltic Sea. Cont. Shelf Res. 88, 228-239. Tóth Zs., Spiess V., Mogollón J. M. and Jensen J. B. (2014b) Estimating the free gas content in Baltic Sea sediments using compressional wave velocity from marine seismic data. J. Geophys. Res.: Solid Earth 119, 8577-8593. Tóth Zs., Spiess V. and Keil H. (2015) Frequency dependence in seismoacoustic imaging of shallow free gas due to gas bubble resonance. J. Geophys. Res.: Solid Earth 120, 8056-8072.

Väli G., Meier H. E. M. and Elken J. (2013) Simulated halocline variability in the Baltic Sea and its implact on hypoxia during 1961-2007. J. Geophys. Res.: Oceans 118, 6982-7000. Wallmann K., Pinero E., Burwicz E., Haeckel M., Hensen C., Dale A. and Ruepke L. (2012) The global inventory of methane hydrate in marine sediments: a theoretical approach. Energies 5, 2449-2498. Whiticar M. J. (2002) Diagenetic relationships of methanogenesis, nutrients, acoustic turbidity, pockmarks and freshwater seepages in Eckernförde bay. Mar. Geol. 182, 29-53. Wunderlich J. and Müller S. (2003) High-resolution sub-bottom profiling using parametric acoustics. Internat. Ocean Syst. 7, 6-11. Yamamoto S., Alcauskas J. B. and Crozier T. E. (1976) Solubility of of methane in distilled water and seawater. J. Chem. Eng. Data 21, 78-80.

Depth (m)

Area (km2)

Areal flux (mmol m-2 d-1)

Total flux (10 mol CH4 d-1)

Total flux (ton C d-1)

HML 0-4 4-6 6-8 8-10 10-12 12-14 14-16 16-18 18-20 >20 Total

2857 407 313 189 134 144 123 57 19 18 4261

0.00 0.21 0.60 0.98 1.36 1.74 2.13 2.51 2.89 4.90

0.00 0.87 1.87 1.85 1.82 2.51 2.62 1.42 0.55 0.88 14.39

0.00 1.05 2.24 2.22 2.18 3.01 3.15 1.70 0.66 1.06 17.27

FGD 0.0-0.5 0.5-1.0 1.0-2.0 2.0-3.0 3.0-4.0 4.0-5.0 >5.0 Total

68 190 179 127 138 96 144 943

5.19 2.36 1.18 0.69 0.49 0.38 0.30

3.53 4.48 2.11 0.88 0.68 0.36 0.42 12.46

4.24 5.37 2.53 1.05 0.81 0.44 0.51 14.95

5

Table 1. Methane flux in the Bornholm Basin sediment calculated by the Holocene Mud Layer (HML) model and the Free Gas Depth (FGD) model. The table shows: depth intervals of the HML and the FGD models used to calculate methane fluxes, total area of each depth interval, methane flux per unit area, and total flux in the Bornholm Basin expressed in mol CH4 and in metric ton Carbon per day.

Figure legends

Figure 1. A) Location of the Bornholm Basin (BB) and Arkona Basin (AB) in the southwestern Baltic Sea. B) Bathymetric map of the Bornholm Basin with seismo-acoustic lines. Numbers on bathymetric contour lines show water depth in m. The red line indicates the seismic transect shown in Fig. 3. C) Sediment coring stations for methane measurements during five expeditions. Red arrow indicates coring stations located on acoustic line GeoB10-044 (cf. Tóth et al., 2014b). D) Methane fluxes at individual stations calculated from sediment core data. (Map: Baltic Sea Hydrographic Commission, 2013, Baltic Sea Bathymetry Database version 0.9.3. Downloaded from http://data.bshc.pro/ on 30.11.2015).

Figure 2. Three types of sulfate and methane distributions in cores from the SW Bornholm Basin taken during cruises AUBO16 and PO392. A) Sulfate minimum and no methane accumulation. Data were combined from a gravity core and Rumohr core at Site AUBO16-BB04 (Beulig et al., 2018). B) Data combined from gravity core and Rumohr core PO392-374200 taken outside of gassy area (see Fig. 3A). The thickness of the Holocene mud layer (HML) is indicated by a double arrow. C) Data from Rumohr core PO392-374190 taken within a gassy area (see Fig. 3A). The linear gradients fitted to the sulfate and methane data define a sulfate-methane transition (SMT) at 0.18 m depth and a free gas depth (FGD) at 0.60 m below which gas bubbles occurred. The thick vertical gray line shows in situ saturation of methane. Black circles: sulfate, open circles: methane. For other lines and symbols, see text. Note the different depth scales.

Figure 3. A) Seismo-acoustic transect using Innomar sediment echo sounder (10 kHz) (red line in Fig. 1B and C) across the Bornholm Basin in direction S to N (left to right). The colored horizons within the seismic profile mark the boundaries between the deep Baltic Ice Lake clay, the Yoldia Sea clay, the Ancylus Lake clay, and the Littorina Sea mud. Core sampling positions and core depths are indicated by red bars. Vertical broken lines and "seis" numbers indicate the seismo-acoustic lines used to compile the figure. B) Methane concentration profiles in the sediment at the same coring positions (indicated by station/core number). All methane panels have the same scale in mM CH4. C) The upward methane flux to the SMT calculated for each core.

Figure 4. A) Methane fluxes calculated from the measured methane concentrations in each core (“Core CH4 flux”) versus the acoustically measured thickness of the Holocene mud layer (HML). The symbols distinguish cores with low methane <10 µM (black circles), with only dissolved

methane (gray circles), or with dissolved methane and in situ gas bubbles (open circles). The regression line is calculated only for data with HML >4 m. B) Methane fluxes calculated by the FGD model for each of the gassy sites plotted against the reciprocal, acoustically determined free gas depth, FGD-acoustic. C) CH4 fluxes calculated from individual cores in the gassy areas plotted against CH4 fluxes calculated from the FGD model at the same sites. D) Free gas depths calculated by the FGD model and free gas depths measured acoustically at individual coring sites, both plotted as 1/FGD (see text).

Figure 5. Map of the Bornholm Basin showing the geographic distributions of: A) thickness of the marine Holocene mud layer (HML, m); B) methane flux calculated from the HML model; C) free gas depth (FGD, mbsf); D) methane flux calculated from the FGD model; E) difference between methane fluxes calculated from the HML model and the FGD model; purple: HML flux = FGD flux (±2); light blue: HML flux < FGD flux; orange-yellow: HML flux > FGD flux. All flux units are in mmol CH4 m-2 day-1.

Figure 6. CH4 fluxes versus depth of the sulfate-methane transition (SMT) determined in each core. The curves show the theoretical FGD model relationship based on 1:1 relationship (black curve) or a 1:0.7 relationship (gray curve) between J SO4 and JCH4 (see text).