Timescales for the development of methanogenesis and free gas layers in recently-deposited sediments of Arkona Basin (Baltic Sea)

Abstract. Arkona Basin (southwestern Baltic Sea) is a seasonally-hypoxic basin characterized by the presence of free methane gas in its youngest organic-rich muddy stratum. Through the use of reactive transport models, this study tracks the development of the methane geochemistry in Arkona Basin as this muddy sediment became deposited during the last 8 kyr. Four cores are modeled each pertaining to a unique geochemical scenario according to their respective contemporary geochemical profiles. Ultimately the thickness of the muddy sediment and the flux of particulate organic carbon are crucial in determining the advent of both methanogenesis and free methane gas, the timescales over which methanogenesis takes over as a dominant reaction pathway for organic matter degradation, and the timescales required for free methane gas to form.


Introduction
Methane, a potent greenhouse gas, is ubiquitously present within marine sediments in dissolved, gaseous and hydrate form. A broad array of evidence suggests that the masses and fluxes of methane in seafloor sediments can vary significantly over time. For instance, there are many features indicative of gas expulsion (e.g. pockmarks, Nelson et al., 1979;Hovland et al., 1984), or remnants of methanogenesis where none exists today (e.g. 13 C-enriched carbonate, Malone et al., 2002). Several studies and models have focused on the sources and sinks of methane in the slope (Davie and Buffett, 2001;Jørgensen et al., 2001;Haeckel et al., 2004;Jørgensen et al., 2004) and on the shelf (Dale et al., 2008a,b;Mogollón et al., 2009;Mogollón et al., 2011;Regnier et al., 2011). The focus of this study, however, is to trace back the development of the methane cycle in shelf sediments of the Baltic Sea. Except for high-latitude regions, shelf systems are generally not favorable for gas hydrate accumulation; they are also much more sensitive to fluctuating conditions in the water-column such as temperature or bottom-water hypoxia. It is thus expected that their dynamics might differ significantly from sediments on the slope.
In shelf sediments underlying shallow and hypoxic water columns, organic matter degradation may completely exhaust sulphate within the first meters of the sediment, allowing for extensive methane production (Jørgensen and Kasten, 2006). In many instances the formation of free methane gas occurs in these settings due to the low solubility concentrations resulting from the shallow water depths (around 10 mM for 50 m depth). The fate of this gas varies depending on the type of sediment where it occurs . Nevertheless, geophysical surveys have revealed widespread areas of acoustic blanking in marine sediments caused by a substantial portion of free gas that is either trapped or moving at slow velocities in the sediment (Wever and Fiedler, 1995;Laier and Jensen, 2007).
The methane produced in sediments can be consumed aerobically in the presence of oxygen, and/or more commonly, by anaerobic oxidation of methane (AOM), a process by which microorganisms oxidize methane utilizing sulphate as the terminal electron acceptor (Iversen and Jørgensen, 1985;Boetius et al., 2000). AOM occurs in a narrow zone of the sediment, termed the sulphate-methane transition zone (SMTZ), and in passive sediments, effectively hinders dissolved methane escape from the sediment. Furthermore, AOM causes methane undersaturation which may prevent free gas from forming near the sediment water interface (SWI), and, may even enhance its dissolution (Dale et al., 2009b;Mogollón et al., 2009;Mogollón et al., 2011;Regnier et al., 2011).
In shelf sediments, present-day biogenic methanogenesis and AOM rates have been quantified in the laboratory (e.g. Crill and Martens, 1983;Treude et al., 2005;Parkes et al., 2007;Knab et al., 2008) and in modeling studies (e.g. Regnier et al., 2011, and references therein), but have generally been confined to the top sediment layers (0-3 m depth). Furthermore, few modeling studies (e.g. Dale et al., 2008b;Arndt et al., 2009) have focused on the long-term geochemical impact of AOM and none have attempted to reconstruct the evolution of methane turnover rates as a result of long-term changes in environmental conditions over millennial timescales. In this context, by simulating the sedimentary history of the methane cycle since its inception, the required time for both the development of a methanogenic zone and eventually a gas forming zone can be estimated. These estimates can then be validated by comparing contemporary measured and simulated profiles. In this study we use methane, sulphate, particulate organic carbon (POC) concentrations and sulphate reduction rates to constrain the modeled reaction rates. The proposed approach could easily be extrapolated to other shelf settings, for instance in cases where pre-Holocene conditions are defined by a hiatus that can be used as a marker for initial conditions.
The aforementioned studies have found that depthintegrated AOM and methanogenesis rates in shallow seas vary by several orders of magnitude, ranging from 0.01 to 10 mol m −2 yr −1 (Crill and Martens, 1983;Treude et al., 2005;Knab et al., 2008;Regnier et al., 2011). While gassy sites have been shown to display greater depth-integrated rates of AOM , published studies so far have not shown discernible differences in methanogenesis rates between gassy and non-gassy sites. This lack of variation may indicate that free gas formation is highly influenced by the methanogenic history of the sediment, and not strictly controlled by higher methanogenesis rates at gassy sites.
The Baltic Sea represents a highly-productive continental shelf setting with a large-scale estuarine-type circulation. Depth-integrated methanogenesis rates as high as 1 mol m −2 yr −1 have been reported in the underlying sediments (Schmaljohann, 1996;Treude et al., 2005). Furthermore, free gas is commonly found in several of the deep water areas in the southern Baltic, such as Arkona and Bornholm Basins, which were created by glacial scouring during the last ice age (Jensen, 1995). In these basins, stratification of the water column due to both dense saline waters entering from the North Sea and an increase in river runoff from deglaciation, coupled to a rise in primary productivity (Bianchi et al., 2000) have led to the deposition and burial of an organic-rich sediment commonly referred to in the literature as Gyttja Clay, Baltic Sapropel, or Holocene organicrich mud (HORM); the term used in this study. At these locations, the degradation of organic carbon within these HORMs becomes the driving force for biogenic methane formation due to the absence of deep methane sources. The aim of this study is to hindcast and quantify the methane cycle since seawater first intruded into the Baltic Sea. Thus, we develop a reactive transport model forced by transient boundary conditions to explore the time-dependent geochemical dynamics within sediments characterized by high organic matter accumulation rates. The model spans the entire depositional history of the HORMs. A lacustrine sediment with no organic carbon or sulphate and methane in the porewater represents the initial condition for the model. The model is calibrated using POC, methane, sulphate, and sulphate reduction rates measured in contemporary sediments. We elucidate the history of the methane cycle at four study sites within Arkona Basin and establish the time and site-dependent magnitude of POC degradation rates, interpret the transient features observed in both methane and sulphate profiles, investigate the methane dynamics toward the deeper part of the methanogenic zone, and determine the time and length scales required to initiate methanogenesis and free gas formation.

Arkona Basin
Arkona Basin is located north of Rügen Island (Germany), west of Bornholm Island (Denmark), and south of mainland Sweden, (Fig. 1). Modern day bottom water salinity in Arkona Basin oscillates from 16 to 20 (Omstedt and Axell, 1998). It is a seasonally-hypoxic basin with a maximum water depth of circa 50 m (Jensen et al., 1999;Moros et al., 2002;Thießen et al., 2006). The sediments exhibit a welldefined gaseous horizon toward the geographical center of the basin (Thießen et al., 2006;Laier and Jensen, 2007) and an organic-rich fluffy layer (5-10 dry wt %) is present at the sediment surface (Schulz and Emeis, 2000).
Significant Holocene post glacial sea-level variations have led to the flooding of Arkona Basin's western barrier (the Darss Sill), producing 4-5 distinct stratigraphical stages (Björck, 1995). (1) The Baltic Ice Lake stage beginning roughly 13.5 kyr before present (BP) as a result of deglaciation (Jensen, 1995;Jensen et al., 1997). During this stage, evidence points toward probable drainage of Arkona Basin waters through Øresund ( Fig. 1) (Jensen, 1995). At circa 10 kyr BP, when major drainage events through south-central Sweden lowered the sea level in the Baltic Sea by 25 m, most of Arkona Basin became exposed (Björck, 1995), although the northeastern part of Arkona Basin may contain shoreline deposits attributed to the (2) Yoldia Sea (Björck, 1995; Lemke (1998). Superimposed is the spatial extent of free gas according to Jensen (2007) -black dashed line, andThießen et al. (2006), white line. Kortekaas et al., 2007). Further deglaciation in the Baltic Sea led to a transgression at around 9.5 kyr BP in Arkona Basin, which was cut off from the North Sea, forming the freshwater (3) Ancylus Lake stage (Björck, 1995;Sohlenius et al., 2001). Later regressions in Arkona Basin during the Ancylus Lake Stage around 8.9 kyr coupled to rising sea level in the Skagerrak/Kattegat led to the pulsated invasion of saline waters into Arkona Basin, and a progressive shift towards the brackish (4) Mastogloia Sea Stage around 8.5 kyr BP (Jensen et al., 1999;Witkowski et al., 2005). The continued sea level rise in the Kattegat since then gave way to the (5) Littorina Sea stage (8 kyr BP until present), during which Arkona Basin became brackish. The exact timing of Arkona Basin's stages is still debated due to the range of 14 C dates recorded by different studies and the unknown reservoir age that can be used to correct measured marine 14 C dates (Sohlenius et al., 1996;Gustafsson and Westman, 2002;Moros et al., 2002;Thießen et al., 2006;Kortekaas et al., 2007). In Arkona Basin, the Ancy- Jensen et al. (1999) Witkowski et al. (2005  lus Lake -Mastogloia Sea transition is not easily discernible in the sediment record and thus requires a detailed diatom stratigraphic analysis that can demonstrate the invasion of brackish-water species (Jensen et al., 1999;Witkowski et al., 2005). The Mastogloia Sea -Littorina Sea transition, on the other hand, can be sedimentologically and geochemically observed, since it is characterized by an increase in the concentration of particulate organic carbon (POC) as well as a change from light-dark laminated gray sediments to a dark grayish-green sediment with decreasing sediment depth (Thießen et al., 2006;Kortekaas et al., 2007), marking the beginning of the HORM deposition.
In this study, we assume that the initial HORM deposition in the southwestern Baltic Sea began 8.0 kyr BP (Jensen, 1995;Jensen et al., 1997Jensen et al., , 1999Witkowski et al., 2005). This date is consistent with the averaged date for the Littorina Sea Stage in the Baltic Sea given by Gustafsson and Westman (2002) if one considers this event in Arkona Basin preceded the averaged date for the Baltic Sea by 500 yr. Figure 2 shows a general overview of the two most recent stages of the Baltic Sea and the chronostratigraphy employed in the model.

Sample collection
Arkona Basin sediments were sampled September 2004 at Stations A1, A3, A5, and A7 ( Fig. 1) to the depths of 129 cm, 330 cm, 49 cm, and 437 cm, respectively, aboard the R/V Gunnar Thorson (Copenhagen). Different coring techniques were applied. Undisturbed surface sediments were sampled at Stations A1 and A7 by use of a Rumohr Lot corer whereas a multi-corer sampling device was deployed at Station A5. Gravity coring was also performed at Stations A1, A3, and A7; nevertheless, these gravity cores may be subject to surface sediment loss (up to 40 cm) as experienced in other expeditions (Dale et al., 2008a;Knab et al., 2008). Thus concentration gradients observed in these cores were aligned with similar observations in the Rumohr Lot cores to obtain "undisturbed" concentration profiles. Samples were analyzed for porosity, methane, sulphate, POC concentrations, and, at Station A3, sulphate reduction rates. Methane was measured using gas chromatography upon sediment sampling, while sulphate was measured using non-suppressed anion exchange chromatography, and POC was measured in acidwashed samples through combustion in a CHN-analyzer (Dale et al., 2008a). Sulphate reduction rates were measured by tracing isotopically-labeled sulphate in incubated sediments (Jørgensen, 1978;Fossing et al., 2000). No groundwater seepage was observed at any of the study sites.

Model setup
A one-dimensional (depth, z) model that tracks the deposition and decomposition of POC (chemically represented by carbohydrate, CH 2 O) with time (t) has been developed. It represents a simplified version of the model in Mogollón et al. (2011) such that compaction is fitted through an exponential porosity function as opposed to calculated based on the effective stress-porosity relation of the sediment. POC degradation is assumed to take place through organoclastic sulphate reduction (Reaction R1) and methanogenesis (Reaction R2), according to the following net stoichiometric reactions: Methane is present in the dissolved and/or gaseous forms with possible exchange between the two phases: Methane is also consumed anaerobically in the presence of sulphate by methanotrophs: The reaction network used in the model considers these four processes and the corresponding chemical species: POC (% dry weight), dissolved methane (mM porewater), sulphate (mM porewater) and free methane gas (gas volume fraction φ g , -). In addition, chloride is simulated as a non-reactive tracer and as a proxy for salinity (see below). Porosity, solid and aqueous phase velocities in Arkona Basin are invariant in time due to the assumption of steady state compaction and the absence of historical markers for impressed fluid flow. However, the biogeochemical reactions and the free gas phase dynamics were modeled as transient phenomena in order to capture the sedimentary geochemical evolution during the HORM deposition. Onedimensional conservation equations for POC concentrations (Eq. 1), aqueous species, C i (i = CH 4 , SO 2− 4 , Cl −1 ) (Eq. 2) and the free gas phase (φ g ) (Eq. 3) were described as follows (Boudreau, 1997;Mogollón et al., 2009): Note that the dependence of each variable on either depth from the SWI (| z ), time (| t ), or both (| z,t ), is explicitly indicated. v s , v a , v g are the burial velocity of solids, the aqueous phase velocity and the gas phase velocity, respectively, and φ is the porosity of the sediment. D * i is the effective diffusion coefficient which is the molecular diffusion constant corrected for tortuosity (Boudreau, 1997). is the universal gas constant, and R gas is the molar rate of free methane gas formation. Biologically induced mixing processes (bioturbation and bioirrigation) and sediment resuspension due to storm events were neglected here since they are essentially decoupled to the much deeper processes of methanogenesis, free gas formation and AOM (Dale et al., 2008a). Besides, Station A5, where the shallowest SMTZ is observed, is characterized by an exponentially decreasing POC profile indicative of a marginal influence of bioturbation.
Steady-state compaction of the sediment with an exponentially decreasing porosity with depth is described as: where β is the depth-attenuation coefficient, v s | 0 is the burial velocity of solids at the SWI, φ| 0 is the porosity at the SWI, and φ| ∞ is the porosity at great depth. Porosity changes due to the presence of a gas phase were assumed to be negligible. With the assumption of steady-state compaction (Eq. 5), no biologically induced sediment mixing, and an exponential, time-invariant porosity profile (Eq. 4), the age of a solid particle undergoing burial (inert or reactive) at any given depth (α| z ) can be derived by integrating v s | z = dz/dα| z and solving for α| z , yielding: where α 0 is the age of the solid particle at deposition. The burial velocity of solids at the SWI for the different stations can be determined similarly to Eq. (6): where Z is the sediment layer thickness and A is the time since first deposition of the HORMs. The burial velocity of solids at each station was determined based on the HORM thickness (Table 3), which was obtained from Lemke (1998). Pore-water advection v a may be influenced by both compaction and externally impressed flows: φ| ∞ is the externally impressed flow below the compacted layer. Due to lack of evidence for externally impressed fluid flow in Arkona Basin sediments, q is assumed to be zero here, which means that both the aqueous phase and the solid phase are buried at the same velocity in fully compacted sediments. Seep formation was not detected at any of the four study sites (or in Arkona Basin as a whole), and thus free gas was assumed to undergo burial in the sediment through: where v g is the free gas advection. The above equation assumes that no gas escapes the sediment. The gas phase can thus be considered as a maximum possible value (see "Results and discussion"). The reaction rate for POC degradation (R POC ) was modeled using a continuum model with a gamma distribution for POC reactivity (Boudreau and Ruddick, 1991): where Q R POC is a parameter that regulates the temperature dependency for organic matter decay and T ref is the reference temperature (Table 2). a and ν are adjustable parameters related to the initial distribution of organic matter reactivities. In this respect, a is the average life-span of the fastest decaying components of POC, that is, high a values indicate that the most labile components degrade more slowly and vice-versa. ν controls the distribution of the most refractive organic matter such that, if ν is low, organic matter is dominated by refractive components. Assuming extracellular hydrolysis to be the rate limiting step controlling the redox geochemistry, we constrain these variables using the methane, sulphate, POC and sulphate reduction rate profiles provided for Station A3. The reactive continuum model has been applied to a range of field data by Boudreau and Ruddick (1991) and a comparison to other models of organic matter degradation can found in Regnier et al. (2011).
The reaction term for sulphate R SO 2− 4 includes both organoclastic sulphate reduction and methanotrophic sulfate reduction AOM (such that the rate of sulphate reduction includes both organoclastic sulphate reduction and methanotrophic sulphate reduction), where AOM is expressed using a bimolecular rate law: where, ρ s is the solid-phase density, n C is the molar mass of carbon, k bi is the bimolecular rate constant for AOM, and K SO 2− 4 is the half-saturation constant for sulphate. f SO 2− 4 | z,t is a factor which controls the partitioning rate of organic matter degradation between organoclastic sulfate reduction and methanogenesis, such that, . Q R AOM is a factor that regulates the temperature dependency of the AOM rate reaction (Table 2).
The reaction rate for methane, R CH 4 , includes methanogenesis, AOM, and free gas formation/dissolution: Gas formation is assumed to be diffusion controlled, whereas gas dissolution is both interface and diffusion controlled (Mogollón et al., 2009). The molar rate of gas formation (R gas ) is thus described as follows: where k int and k dif describe interface-controlled and diffusion-controlled mass transfer respectively. n is the bubble density (assumed constant through the core), c diss is a kinetic mass transfer coefficient, C * CH 4 is the methane saturation concentration, and c λ is the diffusive length boundary constant. Note that given the pressure, temperature and salinity regime in Arkona Bain, methane hydrate cannot form and is thus excluded from the model. The variable notation and  Tables 1 and 2, respectively. The methane saturation concentration, C * CH 4 , depends on the local temperature, pressure and salinity (S = 0.03 + 1.805×C Cl − ×35.45×1.005×10 −3 , where S is unitless and C Cl − is in mM) conditions, and is estimated based on a previously determined algorithm , which is applicable to the conditions observed in Arkona Basin since the beginning of the Littorina Sea Stage (P 3-8 bar, T 273-285 K, S 5-30).
Temperature is modeled explicitly assuming that no production or consumption of heat occurs within the sediment. The weak dependence of the thermal diffusivity on temperature is also neglected . Heat transport through sediments is thus described as (Woodside and Messmer, 1961): where c, ρ, and k represent the specific heat capacity, density, and thermal conductivity for the solid ("s"), aqueous ("a") and gas ("g") phases, and n CH 4 represents the methane molar mass. The first term on the right hand side describes heat conduction and the second term heat advection. Note that, for integration purposes, in Eq. (14) the denominator terms on the right hand side were considered time independent.
Pressure is calculated according to the following equation: where P atm is the atmospheric pressure, g is the gravitational constant, H SL is the water depth, and ρ a is the density of the aqueous phase.

Initial and boundary conditions
Initial conditions for the simulations were based on the assumption that any organic matter deposited during the limnic Ancylus Lake and the Mastogloia Sea stage is highly unreactive. The present cores did not penetrate into the Ancylus Lake stage sediments. Nevertheless, other cores from Arkona Basin have shown that Ancylus Lake sediments have a light gray color typical of oxygen exposure (Kortekaas et al., 2007), have a significantly smaller organic matter content with respect to the Littorina Sea stage (Sohlenius et al., 2001), and, in several intervals, are even characterized by organic-free sands (Jensen et al., 1999). Thus, at the beginning of the simulation (8.5 kyr BP) the sediment is assumed void of reactive POC and methane. It is also assumed that sulphate and chloride were absent during the Ancylus Lake stage, and that the brackish water intrusion during the Mastogloia Sea Stage (8.5-8.0 kyr BP) can be simulated according to a linear increase in concentrations of both variables (e.g. Gustafsson and Westman, 2002). The initial temperature profile was assumed to be constant, with an assigned value equal to the SWI temperature at 8500 yr BP (Fig. 3, upper panel). The POC flux was adjusted concurrent with its reactivity to fit the contemporary organic matter profiles and was assumed constant in time. This assumption is justified based on the observed exponential decrease in the present-day POC data for the Littorina Sea stage (over the last 8 kyr) sediments.
The sea level at any given time (H SL | t ) is: where H SL0 is the current sea level and H SL is the sea level change (in m). Sea level rise at Arkona Basin during the past 8500 yr was mathematically formulated to fit the shoreline displacement curve established for the Holocene by Bennike and Jensen (1998). The following regression was used (Fig. 3, middle panel): where t * is time in yr BP. Boundary conditions at the SWI for sulphate were set as time dependent variables that were adjusted according to paleorecords of the Baltic Sea (Fig. 3, middle panel). The boundary conditions at the SWI for chloride (in mM) were assumed to be 19.3 times those of sulfate (in mM). Variations in both chloride and sulphate in the Arkona Basin are attributed to the salinity gradients observed in bottom waters at Arkona Basin, and are thus related to the migration of dense bottom water currents migrating from the Kattegat and into the Baltic Sea (west to east). To a lesser extent, in Arkona Basin salinity variations may also be due to the mixing with saline-poor waters (∼8) of the Pomeranian Bay, located southeast of Arkona Basin. Chloride and  sulphate were assumed to have fluctuated according to the salinity trends from Gustafsson and Westman (2002). This idealized curve is based on diatom and mollusk assemblages for the entire Baltic Sea, yielding values that are somewhat lower that those currently measured in the bottom waters of Arkona Basin. The temporal trend of these variations, however, should not be significantly different throughout the Baltic. In this respect, we apply the trend in Gustafsson and Westman (2002) by assuming that the historical relative deviations from present values at any location are consistent throughout Arkona Basin. For example, at 5.5 kyr BP the salinity was 1.7 times modern values according to Gustafsson and Westman (2002), thus, we use a 1.7 times correc-tion with respect to the modern salinity values at each station during that period. We also applied the same technique to the bottom water sulphate concentrations, assuming a constant bottom-water chlorinity/sulphate molar ratio of 19.3. Dissolved methane is assumed to have remained constant at 0 mM just above the SWI.
The boundary conditions at the lower limit (L) of the model were assumed to be a no-flux boundary for dissolved species, and L was placed at a depth greater than the diffusive length scale over 8000 yr (L > 4D CH 4 T , T = 8000 yr), that is at L > 28 m. The length of the model sediment column was 30 m in all simulations.  Seppä et al., 2005), in the relative shoreline displacement (e.g. the sea-level drop) with respect to modern-day levels (middle panel -adapted from Bennike and Jensen, 1998), and in the bottom-water sulphate concentration imposed for all 4 cores (bottom panel, adapted from Gustafsson and Westman, 2002) in Arkona Basin during the past 8.5 kyr. Note that the temperature in the top panel represents the inter-annual signal over which a seasonal signal with a temperature amplitude of 4 K was superimposed. The vertical dashed line represents the time (8 kyr BP) when the Holocene organic-rich mud deposition begins.
Millennial changes in temperature at the SWI during the Littorina sea stage were taken from the record of Seppä et al. (2005) for a Swedish lake after filtering out minor oscillations (<1 K) (Fig. 3, top panel). A seasonal forcing (following Dale et al., 2008a) where the temperature amplitude (4 K) was set to the third quartile bottom-water temperature range observed in Omstedt and Axell (1998) was also applied to the top boundary condition for temperature. Note that the geothermal gradient is ignored, which at the bottom boundary would produce an increase in temperature of circa 0.7 K. We therefore assume that the bottom boundary for temperature also shifts according to the trend in Fig. 3, since heat diffusion is most likely able to keep up with the gradual change of a couple of kelvin over millennial scales.

Organic matter deposition rate model
The depositional rate model is based on steady state compaction and the assumption that the accumulation of organic rich sediment during the last 8 kyr occurred under a constant POC flux to the sediment surface. The age of the sediment was thus calculated using Eq. (7) which correctly integrates the effects of compaction into the sedimentation rate.
Not accounting for porosity changes, or not accounting for the integrated effect of the porosity changes in the velocity field can lead to significant errors in calculating the age of solid particles and/or the age of the sediment at any given depth. For example, Fig. 4 compares the measured and calculated sediment age using the porosity data of Kortekaas et al. (2007) according to four different algorithms. In the figure, the filled squares with associated error bars represent measured optically stimulated luminescence ages in Kortekaas et al. (2007). The best-fit solid line represents the predicted age model based on Eq. (6), which properly takes into account the changing burial velocity of solids v s | z . This correct formulation is compared with alternate (incorrect) formulations which include using the depth-dependent velocity (α = z/v s | z , dashed line), which ignores the cumulative effects of compaction. Furthermore, using the burial velocity of solids at the sediment water interface (α = z/v s | 0 , dash dot line) or the burial velocity of solids at great depth (α = z/v s | ∞ , dash dot dot line) completely ignore the effects of compaction.
The fit between the simulated age profiles and the measured age profiles suggests that the assumption of steady state compaction with an exponentially-decreasing porosity profile holds for this temporal interval in Arkona Basin sediments. Although age data for the simulated cores in this study  Kortekaas et al., 2007). Filled squares with associated error bars represent measured ages in Kortekaas et al. (2007). The best-fit solid line represents the predicted age model based on Eq. (6), which properly takes into account the changing burial velocity of solids v s | z . This correct formulation is compared with alternate (incorrect) formulations which include using the depth-dependent velocity without taking to integrated effects into account (α = z/v s | z , dashed line), using the burial velocity of solids at the sediment water interface (α = z/v s | 0 , dash dot line), and using the burial velocity of solids at great depth (α = z/v s | ∞ , dash dot dot line). (b) Porosity measurements (filled squares) and simulated porosity profile (line) for the same core. The porosity profile was used for the above calculations.Note that the maximum depth shown coincides with the HORM thickness.
is lacking, using the HORM thickness map of Lemke (1998) and Eq. (6), the site-specific sedimentation rates can be extracted. The HORM thickness throughout the Arkona Basin varies widely, from less than 1 m to more than 10 m (Fig. 1), and reflects the paleobathymetry since the Mastoglioa Sea stage (Moros et al., 2002). The HORM thicknesses and the fitted porosity profiles at the 4 stations (Fig. 5c, f, i, l), were subsequently used to calcualte the burial velocity of solids at the SWI (v s | 0 ) for the 4 stations can be determined using Eq. (7). This yields values of 0.0480 cm yr −1 , 0.0740 cm yr −1 , 0.190 cm yr −1 and 0.215 cm yr −1 at the A1, A3, A7, and A5 sites, respectively (Table 3). The 4-fold variation in the thickness of the HORM sequence at Arkona Basin translates into v s | 0 values with similar relative variability. The porosity profiles used to calculate v s | 0 were fitted to the data assuming that the porosity at great depth was the same for all sites (φ| ∞ = 0.7). The depth-attenuation coefficient (β) and the porosity at the SWI (φ| 0 ) were varied simultaneously to produce the best visual fit.

POC deposition and degradation
Values for a and ν in the reactive continuum model were initially determined for Station A3, the location with the most complete set of measured data (methane, sulphate, sulphate reduction rate, POC, and porosity profiles). Although the model could be further constrained using DIC, alkalinity, 13 C-methane, and 13 C-dissolved inorganic carbon (e.g. Martens et al., 1999), these dataset were not available for the current stations. The determined values for a and ν were then applied to the remaining stations, based on the findings from Schulz and Emeis (2000), who argued that the organic matter types deposited throughout the basin were similar. The bestfit was obtained using a = 5 yr, indicating a relatively fast decay of the reactive pools, and ν = 0.135, indicating that POC degradation at Arkona Basin has a high apparent order of reaction. This relatively low a value for a basin characterized by high sedimentation rates agrees with the findings of Boudreau and Ruddick (1991), where a decreases with an increase in the sedimentation rate. Our ν value of 0.135 falls in the range reported by the same authors, and compares favorably with ν = 0.125; the value extracted for the POC degradation rate experiments of Westrich (1983). A sensitivity analysis reveals that POC profiles respond strongly to the a parameter only near the SWI where the most reactive components are being degraded (Fig. 6). At greater depth, the POC profiles then shift with respect to the baseline profile, with higher POC concentrations when a value increases (indicating a longer degradation time for the fast-reacting components) and vice-versa (indicating that a significant amount of the fastest-reacting components was consumed earlier). ν, which controls the apparent order of the reaction (Eq. 10), influences the shape of the POC profile throughout the sediment. Higher ν values indicate lower apparent orders of reaction and thus a lower POC degradation dependence on the POC amount, whilst the converse is true for low ν values. Although variations for a and ν can influence the POC profiles, the sulphate reduction rates are barely affected, except within the first few centimeters of the sediment. Note that in Fig. 6, variations in the a and ν parameters were performed over different ranges (400 % and 50 %, respectively). The variability in POC concentrations thus mainly reflects changes in v, especially in the methanic zone where the labile components are no longer present. POC deposition fluxes over the last 8 kyr reveal high spatial variability (0.620, 1.02, 1.80, and 2.70 mol C m −2 yr −1 for Stations A1, A3, A7, and A5 respectively, Table 3). Given that at the four stations the POC reactivity is assumed constant and the POC concentrations asymptote at similar values (around 2 wt %), for this particular model the POC flux scales quantitatively to the burial velocity of solids.

Contemporary geochemical concentrations and turnover rates
Simulated sulphate, methane, and organic carbon concentrations as well as methane gas volume fractions are plotted for all stations along with the corresponding measured data where available (Fig. 5, left panels). Total sulphate reduction rate and AOM rate profiles (Fig. 5, mid panels), and porosities (Fig. 5, right panels) are also shown. The measured data (Fig. 5) provide insight into the contemporary methane dynamics taking place in the upper portions of the HORM. The sampling gear deployed (Rumohr Lot and gravity cores) did not penetrate through the base of the mud at any of the stations. The profiles reflect the relative flux of POC to the seafloor and the depth over which POC mineralization occurs. At the station with the thinnest mud layer (∼169 cm, A1, Lemke, 1998) and lowest POC flux (0.620 mol C m −2 yr −1 ), the rate of organoclastic sulphate reduction was low enough to allow sulphate to persist down to and through the base of the mud layer and, consequently, methane was not detected (Fig. 5a). With the increasing mud thickness at Stations A3, A5, and A7 (Lemke, 1998), sulphate is completely consumed in the SMTZ at 100, 70, and 50 cm, respectively (Fig. 5d, g, j). This allowed dissolved methane to accumulate according to the characteristic sigmoid pattern from zero above the SMTZ until ambient saturation concentrations (e.g. Martens et al., 1999;Dale et al., 2008a;Knab et al., 2008). Note that much of the dissolved methane at Station A7 degassed upon core recovery, leaving a residual pool of with concentration of ca. 2 mM. The shoaling of the SMTZ across the basin is a consequence of both thicker HORM sequences and higher POC fluxes, such that sulphate is removed by both "top down" consump-tion due to organoclastic sulphate reduction and by "bottom up" consumption due to AOM, driven by upward diffusion of methane generated in the deeper portions of the HORM. To illustrate this effect, sulphate reduction rates measured at Station A3 (Fig. 5e) show the typical concave down profile from organic carbon degradation with a peak due to methane-based sulphate reduction indicating the SMTZ (e.g. Jørgensen and Kasten, 2006).
The simulated profiles not only corroborate these observations from the measured data, but also can be used to predict what happens deeper in the sediment. At Station A3, a deep zone of methanogenesis stretches from below the SMTZ to the bottom of the HORM layer, with a maximum methane concentration of ca. 2 mM at 320 cm. Although methanogenesis is active down to the deepest extent of the HORM, POC degradation rates are low at the bottom of the HORM sequence and dissolved methane does not reach saturation levels (ca. 10 mM). Similarly, at Station A7, the methane concentration gradually builds up to 8 mM at 400 cm, roughly 3 mM below the saturation limit. At Station A5 on the other hand (Fig. 5j), which is located within the region of gassy sediments (Thießen et al., 2006;Laier and Jensen, 2007), the simulations reveal that methane exceeds its saturation value of about 11 mM at ca. 280 cm, forming a gas phase with a maximum volume fraction equal to 8 % of total sediment (circa 11 % of the pore space).
At stations where dissolved methane accumulates (A3, A5, A7), the model further predicts that methane concentrations level off near the bottom of the HORM and decrease as they diffuse into the limnic sediment. This net downward diffusion of dissolved methane suggests that the Holocene limnic clay (the Ancylus Lake sediment) acts as a sink for methane. This unusual feature results from the fact that during the early marine transgressions of the Mastogloia Sea stage (8.5-8.0 kyr BP) sulphate diffused into the limnic sediment, and currently acts as the oxidizing agent for downward-diffusing methane (see below). Although geochemical pore water observations in Arkona Basin are mostly restricted to the Littorina Sea stage and thus this feature remains uncorroborated, similarly decreasing methane concentrations below the pre-Littorina sediments in the Bornholm Basin do provide independent support for the model results (Fossing, unpublished data). Simulated present-day values for depth-integrated methanogenesis rates ( MET) at Stations A3, A7, and A5 (0.0431 0.152, 0.431 mol m −2 yr −1 , respectively) are within the range of previously reported values at shallow gassy marine sediments (0.01-10 mol m −2 yr −1 , Table 4). The increasing MET trend from Stations A1, A3, A7 is consistent with the increasing HORM thickness (and therefore increasing sedimentation rates and POC fluxes) among these stations. In comparison, at Cape Lookout Bight, where ebullition has been extensively observed during the summer months (Martens and Klump, 1980), and the sedimentation rates reach values of 10.3 cm yr −1 (Martens et al., 1998), MET has been measured to be at least an order of magnitude higher.

Dynamics of methanogenesis and free gas formation over the Littorina Sea stage
The variation in contemporary geochemical profiles between the 4 study sites is a consequence of the lateral (spatial) heterogeneity caused by the difference in POC fluxes and sedimentation rates. Furthermore, the redox dynamics are also time dependent as illustrated by 6 transient time slices taken from core A5 (Fig. 7). The first time slice (a), at time 8 kyr BP, reflects the conditions at the onset of organic matter deposition. The second time slice (b) illustrates the formation of a sulphate minimum as a result of the transient evolution of organoclastic sulphate reduction (analogous to present-day A1). The third time slice (c) illustrates the formation of a methanogenic zone bounded by two SMTZs, where sulphate diffuses from both the SWI and from below the methanogenic zone. From there on, the two SMTZs move apart as a result of AOM at both redox boundaries. The fourth panel (d) shows further development of the system with a fast-migrating downward methane front and a slowmigrating upward methane front, hindered by downwarddiffusing sulfate (analogous to present-day A3). The fifth panel (e) illustrates an advanced state of the reconstructed history with a thick methanogenic zone, a shallow SMTZ and a methane concentration that is just below gas saturation (analogous to present-day A7). The last panel (f) represents the present-day redox zonation at Station A5, where a permanent gas horizon has formed.

Depth (cm)
Time ( Focusing more on the temporal geochemistry at each station, Fig. 8 shows the transient evolution of the POC, sulphate, methane and gas volume fraction in Stations A1, A3, A7, A5 over the past 8 kyr. Although organic matter deposition at the beginning of the Littorina Sea stage leads to organoclastic sulphate reduction in the upper reaches of the sediment, the combined effects of increasing sulphate concentrations at the SWI and fast diffusion rates allow sulphate to penetrate deep into the sediment. With the progressive accumulation of HORMs and an increase in the thickness of the organoclastic sulphate reduction zone, however, sulphate penetration eventually becomes abated in the limnic sediment. At Station A1 where the POC flux is lowest, sulphate penetrates deep into the modeled domain and remains in the sediment up to present day (Fig. 8). Organoclastic sulphate reduction leads to a sulphate minimum zone during the first century of deposition (see below). Temporal sulphate variations at the SWI are evident from Fig. 8, yet, these changes do not lead to significant fluctuations in the depth-integrated sulphate reduction rates after 1 kyr (Fig. 9, A1) since sulphate is not limiting for POC degradation. Ultimately, the thickness of the HORM at Station A1 is too small to allow methanogenesis.
At Station A3, the development of a methanogenic zone at around 250 cm depth begins at 2.2 kyr BP (Fig. 8) once sulphate is completely reduced. The onset of methanogenesis supplies an additional electron donor for sulphate reduction through AOM. AOM not only takes place at the top of the methanic interval with sulphate diffusing from the SWI, but also at the bottom of the methane zone due to the upward diffusion of relic sulphate below the methanic zone. A similar pattern also occurs at A5 and A7 and is illustrated in detail in Fig. 7 (Station A5) as an example. The expansion of the methanic zone leads to the divergence of two SMTZs or methane fronts (Fig. 9, A3). The shoaling of the top SMTZ (SMTZ1) is due to AOM, which results in a steepening of the sulphate gradient. Since the relic sulphate is not replenished, the deepening of the lower SMTZ (SMTZ2) occurs at a faster rate (Fig. 9, A3). After the onset of methanogenesis and AOM, the shallow SMTZ progressively shoals. Methane concentrations reach 1.8 mM, well below the 10 mM saturation concentration (Figs. 5g, 8). Although double SMTZs have not usually been detected in  Yearly-averaged temporal development (8 kyr BP until present) of depth-integrated rates of sulphate reduction (green), AOM (red) and the SMTZ depth (black) at all stations. Concomitantly with the extension of the methanogenic zone, methane concentrations increase and methane diffuses increasingly closer to the sediment surface and the bottom of the Holocene organic rich mud creating two SMTZs where methane is consumed through AOM. At A5, the free gas depth (FGD) is indicated (blue). Note that the left y-axis governs the rates (green, red) while the right y-axis governs the depths (black, blue). medium sized cores (<5 m), geochemical profiles of long cores (>10 m) collected in the Saanich inlet (Whiticar and Elvert, 2001) show similar bell-shaped methane profiles bordered by sulphate which indicate double AOM fronts. The conditions at Saanich Inlet are similar to those of Arkona Basin such that the lower SMTZ is most likely due to remnant sulphate present in Pleistocene glacial marine clays before the deposition of an organic rich unit. Data from Knab et al. (2008) for Station S11 in the Skagerrak also reveal a drop in the methane concentrations and a slight increase in the sulphate concentrations within the deepest sediment layers, which may indicate methane consumption by an old sulphate source. At Station A7 methanogenesis started 5500 yr BP at 250 cm-depth (Figs. 8 and 9, A7) and thus earlier than at A3, which allows for more methane buildup in the porewater. After 2500 yr of methanogenesis, the methane concentration is within 3 mM of the methane saturation concentration (circa 11 mM).
Due to the higher POC deposition flux, methanogenesis at Station A5 begins at circa 6400 yr BP at 180 cm depth, leading to a methane build up that exceeds the in-situ methane solubility concentration (A5, Figs. 8, and 9). At 3000 yr BP free methane gas phase starts to form at 360 cm. While seasonal variations in bottom-water temperature lead to a recurrent oscillation between methane supersaturation and undersaturation, a free gas phase is present throughout the year and the free gas layer expands as methanogenesis proceeds (not shown).
The contemporary rates of methanogenesis at Arkona Basin are plotted against POC fluxes in Fig. 10. Here, the simulation results are placed in the context of contemporary geochemical zonations in Arkona Basin; that is, zones featuring no methanogenesis, methanogenesis and AOM, methanogenesis, AOM and an ephemeral gas phase, and finally, methanogenesis, AOM and a permanent gas phase. The approximated boundaries between different stages of redox evolution are also delineated, based on simulations using an averaged porosity profile and interpolated burial velocity of solids and POC fluxes. It is important to note that rising temperatures shift the ephemeral gas phase boundary leftwards, which may increase the seasonal gas inventories and increase the potential for gas ebullition and methane escape.

Time scale for methane buildup at Arkona Basin
The accumulation of the most reactive POC pools during the initial stages of the HORM deposition lead to a sharp increase of depth-integrated POC degradation rates in the first centuries at all stations (Fig. 11). The depth-integrated POC degradation rates become larger with time as the HORM accumulates, but the rate increase becomes smaller with time since it reflects the degradation of the more refractory POC pools. Although the burial of the most refractory components provides little additional POC substrate, the ongo- ing sulphate reduction in sediments with HORMs exceeding 1 m is sufficient to deplete sulphate, eventually leading to methanogenesis. Since the HORM accumulation proceeds at the kyr timescale, a sulphate-free sedimentary interval requires 5800, 2500, and 1500 yr at Stations A3, A7, and A5, respectively. Prognostic simulations at A1 reveal that it will take another 500 yr before methanogenesis occurs at this site.
Once methanogenesis starts, it influences directly the sulphate dynamics through AOM. This process contributes toward a rapid shoaling of the SMTZ depth ( Fig. 9, A3, A7, A5), which increases both the thickness of the methanogenic zone and the downward flux of sulphate through diffusion. As a result, the percentage contribution of AOM toward sulphate reduction is commensurate with the depth of the SMTZ. Stations A3, A7, and A5 have a present-day SMTZ depth of 100 cm, 60 cm, and 40 cm, as well as a 13 %, 26 %, and 47 % AOM contribution toward sulphate reduction, respectively. Methanogenesis becomes an increasingly important POC degradation pathway with time as both the SMTZ becomes shallower and the HORM becomes thicker. Nevertheless, while the thickness of the methanic zone may exceed that of organoclastic sulphate reducing zone 10 fold (e.g. present day values at A5 of ∼60 cm vs. ∼600 cm, respectively) organoclastic sulphate reduction still contributes to the major portion of the anaerobic POC degradation pathway.
Thus, the shoaling of the SMTZ depth coupled to the accretion of HORMs deeper into the sediments means that the relative contribution of methanogenesis toward total POC degradation becomes progressively more important with time. The stations with thicker HORMs, higher POC fluxes, and shallower SMTZs exhibit a higher percentage methanogenesis (Fig. 11). The imbalance between depth-integrated methanogenesis and depth-integrated AOM rates (in both SMTZs) dictates the depth-integrated rate of methane accumulation. In Arkona Basin, the net depth-integrated methane accumulation is 10-15 % of the depth-integrated methanogenesis rates. Therefore, methane accumulation in the sediment increases steadily resulting in a concentration increase and concomitant methane diffusion toward both SMTZs, located near the sediment surface and below the HORM, respectively. The top SMTZ depth stabilizes after the initial shoaling when the downward diffusing sulphate flux from the SWI is balanced by the depth-integrated sulphate reduction rate. Downcore, however, the methane front continues to propagate downward at a rate of approximately 3.5 m kyr −1 , 6 m kyr −1 and 7 m kyr −1 at Stations A3, A7, and A5, respectively, since there is no sulphate replenishment in the pre-Holocene sediments that can stabilize the downward methane front. At Station A5 simulations predict that free methane gas has been present over the last 3.7 kyr, that is, 2.2 kyr after the first appearance of dissolved methane. This result indicates that although gas formation likely requires high methanogenesis rates, it also requires a sufficient time span in order for methane to build up to saturation concentrations. The time required for the methane concentration to reach the respective saturation concentration will thus be highly variable, depending on the local methanogenesis rates, as well as the pressure, temperature, and salinity conditions which govern the methane saturation concentration.
Once supersaturation in the sediment is reached, further increases in dissolved methane are abated by the production of free methane gas . Within the free gas zone, free gas production thus helps maintain a dissolved methane profile which closely follows the methane saturation curve, which is relatively invariant compared to the methane profile above the FGD (Fig. 7, final panel). In this respect, upward methane diffusion out of the gas zone is substantially reduced and most methane generated in the sediment ends up in the gas phase. Presently at Station A5, a total of 0.042 mol m −2 yr −1 go toward gas production. If this value is averaged over the gas forming zone (3.1 m), and translated into percentage gas according to local pressure (ca. 5.8 bar) and temperature conditions (278 K), then the contemporary depth-averaged gas production amounts to 0.0054 % yr −1 . It is important to note, however, that relic instances of gas ebullition and/or buoyant interstitial gas movement in the sediment (processes not taken into account in the present model) may lead to a lowering of the calculated methane gas accumulation rate. In the case of gas ebullition, the gas accumulation rate would be lowered by the amount of gas that has escaped unto the overlying water column, whereas in the case of interstitial gas movement, these rates would diminish by the rate of gas advection, and dissolution in the AOMinduced methane undersaturated zone.

Conclusions
Arkona Basin is a seasonally-hypoxic marine system in the southwestern Baltic Sea characterized by an organic-rich Holocene mud layer which increases in thickness toward the center of the depression. Areas with the thickest mud layer have free methane gas in the deeper parts of the muddy sediment. Millennial-scale simulations capturing the deposition of the organic-rich mud since the onset of the present brackish stage of the Baltic Sea reveal that the flux of organic carbon and the burial velocity of solids ultimately determine the propensity toward methanogenesis and free gas formation. They furthermore coincide with previous studies (e.g. Moros et al., 2002) in indicating that sedimentation rates have been higher toward the geographical gas-rich center of the basin. Our simulations show that methanogenesis only occurs in the basin where the mud thicknesses exceeds 3 m. In addition, at 45 m water depth, the methane solubility exceeds 10 mM, and thus, both high burial velocities (>0.16 cm yr −1 ) and high organic matter fluxes (>1.9 mol C m −2 yr −1 ) are required for gas to form. Although the organic matter reactivity may have slight variations throughout the basin, it is the cross-basin heterogeneity in sedimentation rate and organic matter flux which drives the formation of gas.
At present, between 0 and 45 % of the total anaerobic degradation rate of organic matter is due to methanogenesis, with the remainder attributable to sulphate reduction. These high methanogenic fractions may be related to the ongoing dynamics of the methanogenesis zone. Once sulphate is consumed to levels low enough to allow methane to accumulate, the methanogenic zone rapidly expands, channeling more POC degradation through fermentation pathways. As AOM increases because of the upward shift in the SMTZ, organoclastic sulphate reduction rates decrease and more sediment is exposed to degradation by methanogenesis. The spreading of this zone, however, eventually becomes limited as AOM and organoclastic sulphate reduction are balanced by increased bottom-water sulphate diffusion induced by the sharpening of the sulphate gradient, thus reducing the rate at which the SMTZ approaches the sediment surface. The time required to reach this quasi-steady state is around 3-4 kyr.
The simulations can also be used to elucidate spatiotemporal information that cannot be analyzed from measured data alone, such as methane diffusion down through the base of the Holocene mud into organic-poor freshwater clays, and the time-scale required for both methane and free gas formation since the onset of high organic matter depositional fluxes. Eventually, methane production and consumption rates can be scaled to parameters such as the sedimentation rate and the organic matter flux in order to spatially resolve and quantify the regional methane turnover rates in a given location. The results from these models can also be compared to similar models developed for the slope and the deep sea in order to further decode the global dynamics of methane formation.