Assessing the potential for non-turbulent methane escape from the East Siberian Arctic Shelf

The East Siberian Arctic Shelf (ESAS) hosts large yet poorly quantified reservoirs of subsea permafrost and associated gas hydrates. It has been suggested that the global-warming induced thawing and dissociation of these reservoirs is currently releasing methane (CH₄) to the shallow coastal ocean and ultimately the atmosphere. However, a major unknown in assessing the contribution of this CH₄ flux to the global CH₄ cycle and its climate feedbacks is the fate of CH₄ as it migrates towards the sediment–water interface. In marine sediments, (an)aerobic oxidation reactions generally act as a very efficient methane sink. However, a number of environmental conditions can reduce the efficiency of this biofilter. Here, we used a reaction-transport model to assess the efficiency of the benthic methane filter and, thus, the potential for benthic methane escape across a wide range of environmental conditions that could be encountered on the East Siberian Arctic Shelf. Results show that, under steady-state conditions, anaerobic oxidation of methane (AOM) acts as an efficient biofilter. However, high CH₄ escape is simulated for rapidly accumulating and/or active sediments and can be further enhanced by the presence of organic matter with intermediate reactivity and/or intense local transport processes, such as bioirrigation. In addition, in active settings, the sudden onset of CH₄ flux triggered by, for instance, permafrost thaw or hydrate destabilization can also drive a high non-turbulent methane escape of up to 19 µmol CH₄ cm⁻² yr⁻¹ during a transient, multi-decadal period. This “window of opportunity” arises due to delayed response of the resident microbial community to suddenly changing CH₄ fluxes. A first-order estimate of non-turbulent, benthic methane efflux from the Laptev Sea is derived as well. We find that, under present-day conditions, non-turbulent methane efflux from Laptev Sea sediments does not exceed 1 Gg CH₄ yr⁻¹. As a consequence, we conclude that previously published estimates of ocean–atmosphere CH₄ fluxes from the ESAS cannot be supported by non-turbulent, benthic methane escape.

1 Introduction

The Siberian shelf represents the largest shelf on Earth (∼3 millions km² Wegner et al., 2015) and spreads from the Kara Sea to the Laptev, East Siberian and Chuckhi seas. The East Siberian Arctic Shelf (ESAS) corresponds to the broad area beneath the shallow (∼45 m water depth, James et al., 2016) Laptev and East Siberian Arctic Sea (Romanovskii et al., 2004; Shakhova et al., 2010a) and represents the largest region on the Siberian shelf (Romanovskii et al., 2005), covering about 25 % of the total Arctic shelf (Shakhova et al., 2010a).

Although similar in many aspects to other shelf environments, a distinguishing feature of the ESAS is the presence of subsea permafrost and associated gas hydrates buried in the sediment (Sloan and Koh, 2007; Romanovskii et al., 2005). Subsea permafrost is a terrestrial relict that mainly formed during glacial periods, when retreating sea levels (with a minimum of 120 m below the current level around the Last Glacial Maximum) exposed Arctic shelves (Fairbanks, 1989; Bauch et al., 2001). Under these conditions, permafrost aggraded on the shelf and was subsequently submersed (Romanovskii and Hubberten, 2001; Romanovskii et al., 2005) by rising sea levels during the Holocene sea transgression (12–5 kyr BP) (Bauch et al., 2001). Gas hydrates are solid, methane-concentrated formations in which a gas molecule is trapped in a cage of water molecules (Ruppel and Kessler, 2017). They are thermodynamically stable under specific temperature–pressure–salinity in the ocean floor, including areas beneath the subsea permafrost (Sloan and Koh, 2007).

Little is known about the total amount of carbon stored in subsea permafrost or its original extent. Published estimates of carbon reservoir sizes diverge by orders of magnitude. For instance Shakhova et al. (2010a) estimate that 1175 Pg C are locked in subsea permafrost on the ESAS alone, while McGuire et al. (2009) calculate that, across the entire Arctic shelf, 9.4 Pg C reside in upper sediments and 1.5–49 Pg C (2–65 Pg CH₄) in methane gas hydrates. Thus, the size of the Arctic subsea permafrost reservoir, its spatial distribution, as well as its biogeochemical and physical characteristics remain poorly known.

These knowledge gaps are critical as climate change is amplified in polar regions. The Arctic is currently warming at a rate twice as fast as the global mean (Trenberth et al., 2007; Bekryaev et al., 2010; Jeffries and Richter-Menge, 2012; Christensen et al., 2013). Recent observations indicate that bottom water temperatures in the coastal and inner shelf regions of the ESAS (water depth <30 m, Dmitrenko et al., 2011) are rising, while the central shelf sea may be subject to intense episodic warming (Janout et al., 2016). The increasing influx of warmer Atlantic water into the Arctic Ocean – the so-called Atlantification (Polyakov et al., 2017; Barton et al., 2018) – will not only further enhance this warming, but will also influence circulation and salinity patterns on the shelf (Carmack et al., 1995; Zhang et al., 1998; Biastoch et al., 2011). At the same time, it has been long recognized that the Arctic is a potential hotspot for methane emissions. Extensive methane gas bubbling has been observed in the Laptev Sea and has been directly linked to these environmental changes (Shakhova et al., 2010b, 2014). Shakhova et al. (2014) suggest that warming induced subsea permafrost thaw and hydrate destabilization may support methane emissions of up to 17 Tg CH₄ yr⁻¹ from the ESAS alone. Projected change in temperature (Shakhova et al., 2017, 2019) due to climate change is expected to further destabilize Arctic subsea permafrost and gas hydrate reservoirs and might thus enhance further methane emissions (Piechura and Walczowski, 1995; Westbrook et al., 2009; Reagan and Moridis, 2009; Biastoch et al., 2011; Hunter et al., 2013; Drake et al., 2015; Ruppel and Kessler, 2017). However, a number of recent studies have questioned the significance of subsea permafrost thaw and hydrate destabilization for methane efflux from Arctic sediment (Thornton et al., 2016; Ruppel and Kessler, 2017), for methane concentrations in Arctic Ocean waters (Overduin et al., 2015; Sapart et al., 2017) and, ultimately, for methane emissions from the Arctic waters (Ruppel and Kessler, 2017; Sparrow et al., 2018). Thus, the contribution of subsea permafrost thaw and gas hydrate destabilization to methane emissions from the warming Arctic shelf and, ultimately, methane–climate feedbacks remains poorly quantified (James et al., 2016; Saunois et al., 2016). As a consequence, it has not received much attention in the recent IPCC special report (Masson-Delmotte et al., 2018). At present, a major unknown is the strength of methane sinks in Arctic sediments and waters and their influence on methane emissions (Ruppel and Kessler, 2017). Therefore, improved assessments of the present and future climate impact of permafrost thaw and hydrate destabilization require not only a better knowledge of the Arctic subsea permafrost and hydrate distribution, reservoir size and characteristics, but also a better quantitative understanding of Arctic methane sinks.

In marine sediments, upward migrating methane is generally efficiently consumed by the anaerobic oxidation of methane (AOM) and, to a lesser extent, the aerobic oxidation of methane (AeOM) (Hinrichs and Boetius, 2002; Reeburgh, 2007; Knittel and Boetius, 2009). Although the exact AOM process has not been fully understood yet (James et al., 2016; McGlynn et al., 2015; Milucka et al., 2012; Wegener et al., 2015; Dean et al., 2018), it is thought that AOM is mediated by a consortium of methane-oxidizing Archaea and sulfate-reducing bacteria (Boetius et al., 2000) according to the reaction (Hinrichs and Boetius, 2002; Krüger et al., 2003)

$$\begin{array}{}\text{(1)}& {\mathrm{CH}}_{\mathrm{4}}+{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\to {\mathrm{HCO}}_{\mathrm{3}}^{-}+{\mathrm{HS}}^{-}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}.\end{array}$$

A recent assessment indicates that, in global sediments, around 45–61 Tg CH₄ yr⁻¹ (Egger et al., 2018) are consumed by AOM, thus significantly reducing previously published estimates of 320–360 Pg CH₄ yr⁻¹ (Hinrichs and Boetius, 2002; Reeburgh, 2007).

AOM generally acts as a particularly efficient biofilter for upward migrating methane and oxidizes up to 100 % of the methane flux coming from below (e.g. Regnier et al., 2011). However, a number of environmental conditions can reduce the efficiency of this AOM biofilter, allowing methane to escape from the sediment (Iversen and Jorgensen, 1985; Piker et al., 1998; Jørgensen et al., 2001; Treude et al., 2005; Knab et al., 2008; Dale et al., 2008c; Thang et al., 2013; Egger et al., 2016). It has been shown that, in particular, high sedimentation rates (Egger et al., 2016), slow microbial growth (Dale et al., 2006, 2008c) or the accumulation of free gas can promote methane efflux from the sediment. These findings are particularly relevant for potential methane escape from Arctic shelf sediments. The Siberian shelf is the largest sedimentary basin in the world (Gramberg et al., 1983) and shelf areas close to the large Arctic rivers reveal sedimentation rates than can be up to 5 times faster than rates that are typically observed in the ocean (Leifer et al., 2017). In addition, the Arctic shelf is subject to large seasonal, as well as climate-induced long-term changes in environmental conditions, namely ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration in seawater and availability of CH₄ in the sediments coming from deeper strata. These factors may influence the efficiency of the AOM biofilter through their effect on microbial biomass dynamics. Finally, observations from the ESAS also indicate that methane gas accumulates in the sediments. When free gas pockets grow enough, methane tends to migrate upwards along pathways with higher permeability or where fractures occur (Yakushev, 1989; Boudreau et al., 2005; Wright et al., 2008; Shakhova et al., 2014, 2015, 2017; Leifer et al., 2017) and might even crack the sediments themselves (O'Connor et al., 2010; Overduin et al., 2016; Yao et al., 2019; Stranne et al., 2019). However, despite a wealth of AOM-related research, a holistic, quantitative evaluation of the most important environmental controls on the efficiency of the AOM biofilter and its impact on methane escape from marine sediments is currently lacking. Thus our ability to understand and quantify AOM sink in ESAS sediments and thus the climate impact of subsea permafrost thaw and gas hydrate destabilization is seriously compromised.

Therefore, we here use a one-dimensional reaction-transport model approach to understand and quantify the efficiency of the AOM biofilter and its influence on the potential benthic release of methane in response to a plausible range of upward migrating dissolved methane fluxes from thawing permafrost and/or dissociating methane gas hydrates on the warming ESAS. The developed model accounts for the most pertinent primary and secondary redox processes, as well as mineral precipitation, methane gas formation and fast equilibrium reactions. Both active sites (characterized by an upward water flow) and passive sites (without an upward water flow) are investigated. We limit our model analysis to non-turbulent methane efflux, because methane in gaseous form is not directly accessible for the AOM community. As a consequence, free gas bubbles are less prone to be consumed by AOM and methane gas either is trapped in the sediments or rapidly migrates upwards through cracks, faults or fractures (Boudreau, 2012), bypassing the AOM biofilter.

The model is forced with a set of boundary conditions that are broadly representative of the conditions potentially encountered on the ESAS. It is applied to conduct a comprehensive one-at-a-time, steady-state sensitivity study over the entire plausible range of (1) sedimentation rates, (2) active fluid flow velocities, (3) AOM rate constants, (4) organic matter reactivity and (5) non-local transport activity encountered on the ESAS. In addition, we also evaluated the influence of (1) seasonal variability and (2) idealized, projected climate change on the efficiency of the AOM biofilter and potential non-turbulent methane escape from the ESAS under transient conditions. For this purpose, the model is extended by adopting an explicit description of AOM biomass dynamics and a bioenergetic rate law for AOM (Dale et al., 2006, 2008c, b). Finally, the results of all sensitivity study runs are used to identify the most important controls on methane efflux and derive a transfer function that allows establishment of a first-order estimate of methane escape from the ESAS.

The specific aims of this work are thus (1) to identify and quantitatively understand the most important environmental controls on the efficiency of the AOM biofilter, as well as (2) its significance in reducing upward migrating methane fluxes originating from thawing subsea permafrost or destabilizing methane gas hydrates under a plausible range of environmental conditions encountered on the present and future Siberian shelf. Model results are then used to (3) identify environmental conditions (and thus areas on the ESAS) that favour non-turbulent dissolved methane fluxes across the sediment–water interface and (4) derive transfer functions that allow estimation of the potential for non-turbulent CH₄ escape from ESAS sediments, thus providing first-order constraints on the Arctic methane budget.

2 Methods

2.1 BRNS: reaction-transport model

The Biogeochemical Reaction Network Simulator (BRNS) (Regnier et al., 2002; Aguilera et al., 2005; Centler et al., 2010) – an adaptive simulation environment suitable for simulating large, mixed kinetic-equilibrium reaction networks in porous media (e.g. Jourabchi et al., 2005; Thullner et al., 2005; Dale et al., 2009) – is used to quantitatively explore the fluxes and transformations of methane in a sediment column representative for ESAS conditions. For this purpose, we set up a reaction network (Tables S1, S2 in the Supplement), model parameters (Table S6), as well as boundary conditions (Table S7) that cover the conditions encountered on the present-day Siberian shelf.

In the BRNS, the general mass conservation for each solid and dissolved species is described by a set of coupled advection–diffusion–reaction equations in porous media which are solved simultaneously (e.g. Berner, 1980; Boudreau, 1997; note that dependencies on z and t have been omitted for simplicity):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathit{\xi }{C}_i}{\partial t}}& ={\displaystyle \frac{\partial }{\partial z}}\left[(D_i+D_{\mathrm{b},i})\mathit{\xi }{\displaystyle \frac{\partial C_i}{\partial z}}\right]\\ & -{\displaystyle \frac{\partial }{\partial z}}\left(v\mathit{\xi }{C}_i\right)+{\mathit{\alpha }}_i\mathit{\xi }\left(C_i\right(\mathrm{0})-C_i)+{\mathcal{S}}_i.\end{array}\end{array}$$

C_i is the concentration of the species i (mass per porewater volume for dissolved species or mass per solid matrix volume for a solid species); ξ i.e. the porosity ξ=φ for dissolved species and $\mathit{\xi }={\mathit{\phi }}_s=\mathrm{1}-\mathit{\phi }$ for solid species. D_i is the effective diffusion coefficient for species i and is affected by salinity, temperature and tortuosity (see Table S6). D_b denotes the bioturbation coefficient and v is the advective velocity. For solid species v=ω with ω being the burial rate, while the advective velocity for dissolved species is given by the sum of the burial rate and an advective flow velocity, v_up, i.e. $v=\mathit{\omega }+v_{\mathrm{up}}$. A site where v_up≠0 is defined as an active site, while a site with no advective upward water flow is defined as passive. α_i is the bioirrigation coefficient (α_i=0 for solid species) and C_i(0) is the concentration of the species i at the sediment–water interface (SWI). The reaction term 𝒮_i is written as

$$\begin{array}{}\text{(3)}& {\mathcal{S}}_i=\sum _j{\mathit{\lambda }}_{ij}{R}_j,\end{array}$$

where λ_ij are the stoichiometric coefficients of all reaction rates R_j that affect species i.

2.1.1 Transport

The effective diffusion coefficients D_i are determined by correcting the diffusion coefficients in free solution $D_i^{\mathrm{0}}$ (Boudreau, 1997) for tortuosity θ and temperature. Tortuosity is calculated by means of porosity φ according to a modified Weissberg relation (Boudreau, 1997): $\mathit{\theta }=\mathrm{1}-\ln \left({\mathit{\phi }}^{\mathrm{2}}\right)$. Note that the effective diffusion coefficients used in the model neglect pressure effects. Following Dale et al. (2008a), migration of methane gas is simply parameterized via a pseudo-diffusive term, with an apparent gas diffusion coefficient, $D_{{\mathrm{CH}}_{\mathrm{4}}}\left(g\right)$. Bioturbation in the upper decimetres of the sediment is simulated using a diffusive term (e.g. Boudreau, 1986), with a constant bioturbation coefficient, $D_b^{\mathrm{0}}$. The model assumes that bioturbation ceases at the bioturbation depth, z_bio (Boudreau, 1997). Bioirrigation is included in the mass conservation equation as a source or a sink function analogous to a kinetic rate. It is calculated as the product of the irrigation intensity, α (α=0 for all solids), and the difference in concentration of species i relative to the concentration at the SWI, C_i(0). The bioirrigation rate α, is evaluated from the bioirrigation coefficient at the sediment surface (α₀) and the bioirrigation attenuation depth (z_irr) and is given by Eq. (S9) in the Supplement. Porosity is assumed to decrease with depth according to an exponential decay (Athy, 1930):

$$\begin{array}{}\text{(4)}& \mathit{\phi }\left(z\right)={\mathit{\phi }}_{\mathrm{0}}{e}^{-c_{\mathrm{0}}z},\end{array}$$

with φ₀ the porosity at the SWI and c₀ the typical length scale for compaction.

2.1.2 Biogeochemical network

The reaction network implemented here (33 species, 37 reactions) encompasses the most pertinent primary and secondary redox reactions, equilibrium reactions and mineral precipitation and adsorption reactions. A summary of the reactions, their stoichiometry and their rate formulations can be found in Tables S2 and S3. The following section provides a short description of the implemented reaction network, as well as a more detailed description of the reactions that affect the production/consumption of methane. A complete description can be found in the Supplement.

The BRNS model accounts for the degradation of organic matter by aerobic degradation, denitrification, manganese oxide reduction, iron reduction, sulfate reduction and methanogenesis (Table S2). Organic matter degradation is described by means of the reactive continuum model (RCM) (Aris, 1968; Ho and Aris, 1987; Boudreau and Ruddick, 1991) that describes compound-specific reactivities (Tesi et al., 2014) and, thus, captures the widely observed decrease in apparent organic matter reactivity with degradation state. The relative importance of each metabolic pathway is simulated through a series of kinetic limitation terms, reflecting their sequential utilization in the order of their decreasing Gibbs energy yields (Table S1). After all terminal electron acceptors (TEAs) are consumed, the remaining organic matter may be degraded by methanogenesis. The rates of secondary redox reactions (Table S3), are described by bimolecular rate laws (e.g. Wang and Van Cappellen, 1996). Adsorption reactions are considered as fast equilibrium processes (Table S3, R28–R31 in the Supplement). Mineral precipitation rates are simulated according to kinetic–thermodynamic rate laws (Table S3, R16–R24).

As described above, methane is produced during organic matter degradation by methanogens in deeper sediment layers, once all TEAs are depleted (Table S2, R6). If the concentration of dissolved methane exceeds the saturation concentration [CH₄]^*, methane gas forms. The transfer rate of methane between the dissolved and gaseous phase is linearly controlled by the departure of the simulated dissolved methane concentration from the saturation concentration (Haeckel et al., 2004; Hensen and Wallmann, 2005; Tishchenko et al., 2005; Mogollón et al., 2009; Graves et al., 2017). [CH₄]^* is calculated according to Dale et al. (2008a), derived from the formulation proposed by Duan et al. (1992) for which [CH₄]^* depends on in situ salinity, pressure and temperature. Here, we assume that the formed methane gas is inaccessible to microbial activity and hence bypasses anaerobic and/or aerobic oxidation zones. In contrast, dissolved methane can be consumed by anaerobic (AOM) or aerobic oxidation of methane (AeOM). Free gas can re-dissolve into porewater once porewater methane concentration falls below the saturation level and may then become available to methanotrophs. AeOM rate is simply described by a bimolecular rate law (Table S3, R14). The description of AOM depends on the model scenario. For steady-state simulations, we apply a simple bimolecular rate:

$$\begin{array}{}\text{(5)}& {\mathrm{rate}}_{\mathrm{AOM}}=k_{\mathrm{AOM}}\left[{\mathrm{CH}}_{\mathrm{4}}\right]\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right].\end{array}$$

It is the simplest and most commonly used formulation of the AOM rate in reaction-transport models (e.g. Regnier et al., 2011). It accounts for kinetic controls and assumes that, under steady-state conditions, bioenergetic controls are negligible (Dale et al., 2006; Regnier et al., 2011).

For transient model simulations, we apply a bioenergetic rate law in combination with an explicit description of the AOM-performing biomass (Dale et al., 2006, 2008c). It has been shown that the rates of redox reactions, whose energy yield is used by micro-organisms to grow, can be coupled to biomass growth rates via a kinetic Monod term and a thermodynamic Boltzmann term (e.g. Rittmann and VanBriesen, 2019). Hence, the time derivative of AOM-performing biomass (B) can be written as

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}B}{\mathrm{d}t}}={\mathit{\mu }}_{\mathrm{g}}B\cdot F_{\mathrm{K}}\cdot F_{\mathrm{T}}-{\mathit{\mu }}_{\mathrm{d}}{B}^{\mathrm{2}},\end{array}$$

where μ_g is the growth rate and μ_d is the decay rate. F_K is the kinetic constraint given by

$$\begin{array}{}\text{(7)}& F_{\mathrm{K}}={\displaystyle \frac{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}{K_{\mathrm{m}}^{{\mathrm{CH}}_{\mathrm{4}}}+\left[{\mathrm{CH}}_{\mathrm{4}}\right]}}\cdot {\displaystyle \frac{\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]}{K_{\mathrm{m}}^{{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}}+\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]}},\end{array}$$

with $K_{\mathrm{m}}^{{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}}$ half saturation constant of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and $K_{\mathrm{m}}^{{\mathrm{CH}}_{\mathrm{4}}}$ half saturation constant of CH₄, according to a typical Michaelis–Menten formulation for enzymatically catalysed reactions. F_T represent the thermodynamic limitation and is given by

$$\begin{array}{}\text{(8)}& \left\{\begin{array}{ll}\mathrm{1}-\exp \left(\frac{\mathrm{\Delta }{G}_{\mathrm{r}}+\mathrm{\Delta }{G}_{\mathrm{BQ}}}{\mathit{\chi }RT}\right),& \text{if }\frac{\mathrm{\Delta }{G}_{\mathrm{r}}+\mathrm{\Delta }{G}_{\mathrm{BQ}}}{\mathit{\chi }RT}<\mathrm{0}\\ \mathrm{0},& \text{if }\frac{\mathrm{\Delta }{G}_{\mathrm{r}}+\mathrm{\Delta }{G}_{\mathrm{BQ}}}{\mathit{\chi }RT}>\mathrm{0},\end{array}\right.\end{array}$$

where R is the gas constant, T is the absolute temperature, χ is the average number of electrons transferred per reaction per mole of ATP produced (Jin and Bethke, 2005), ΔG_r is the Gibbs free energy of the reaction and ΔG_BQ=20 kJ (mol e⁻)⁻¹ is the minimum energy needed to support synthesis of $\sim \frac{\mathrm{1}}{\mathrm{3}}-\frac{\mathrm{1}}{\mathrm{4}}$ mol ATP (Dale et al., 2008c). In order to be thermodynamically favourable, the total energy ΔG_r+ΔG_BQ has to be negative, meaning that Gibbs free energy provided by the catabolic reaction is sufficient to sustain the microbial biomass growth. ΔG_r is given by

$$\begin{array}{}\text{(9)}& \mathrm{\Delta }{G}_{\mathrm{r}}=\mathrm{\Delta }{G}_{\mathrm{r}}^{\mathrm{0}}+RT\ln \left(\mathit{\gamma }{\displaystyle \frac{\left[{\mathrm{HS}}^{-}\right]\cdot \left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]}{\left[{\mathrm{CH}}_{\mathrm{4}}\right]\cdot \left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]}}\right),\end{array}$$

with $\mathrm{\Delta }{G}_{\mathrm{r}}^{\mathrm{0}}$ the standard free energy of the reaction, the second term the deviation from standard conditions (temperature and reaction quotient) on Gibbs free energy, and γ a parameter representing departure from ideal behaviour.

The link between substrate consumption and microbial growth (anabolism) is given by Dale et al. (2006)

$$\begin{array}{}\text{(10)}& \begin{array}{rl}\mathrm{13.8}\mathrm{SD}\cdot {\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}& +\mathrm{14.3}\mathrm{SD}\cdot {\mathrm{CH}}_{\mathrm{4}}+\mathrm{0.2}\mathrm{SD}\cdot {\mathrm{NH}}_{\mathrm{4}}^{+}\\ & +\mathrm{0.3}\mathrm{SD}\cdot {\mathrm{H}}^{+}\to \mathrm{0.2}\mathrm{B}\\ & +\mathrm{13.3}\mathrm{SD}\cdot {\mathrm{HCO}}_{\mathrm{3}}^{-}+\mathrm{13.8}\mathrm{SD}\cdot {\mathrm{HS}}^{-},\end{array}\end{array}$$

assuming that the cellular composition of the biomass B is equal to C₅H₇O₂N (Bruce and Perry, 2001; Dale et al., 2006, 2008c; Rittmann and McCarty, 2012). SD $=(\mathrm{1}-\mathit{\phi })/\mathit{\phi }$ is the conversion factor between dissolved and solid species, here represented by micro-organisms (which are assumed to be attached to the solid matrix). Catabolism is linked to biomass growth (anabolism) through the growth yield. We apply a yield of 0.0713 (Dale et al., 2006), which falls at the upper end of reported AOM growth yields, i.e. 0.05–0.07 (Dale et al., 2006; Nauhaus et al., 2007).

2.1.3 Boundary conditions

Boundary conditions place the model in its environmental context. For dissolved species, constant bottom water concentrations (Dirichlet boundary conditions) are applied at the sediment–water interface, while a known flux condition (Neumann boundary condition) are applied for solid species. At the lower boundary, a zero gradient flux boundary condition ($\partial C/\partial z=\mathrm{0}$) is considered for all species except methane, for which a Dirichlet condition is specified to account for methane supplied from thawing permafrost and/or dissociating gas hydrates below.

2.2 Model evaluation

To evaluate the performance of the BRNS set-up in capturing the main diagenetic patterns observed in Arctic shelf sediments we run the model for two case study sites in the area of interest: (1) a site offshore Kotelny Island in the central region of the Laptev Sea, north of the Lena River delta (76.171^∘ N, 129.333^∘ E, 56 m water depth) collected during the SWERUS-C3 expedition in summer 2014 (Brüchert et al., 2018; Brüchert, 2020).

Although observations are merely available for the first 22 cm, the first 3 m of sediment are simulated to allow for the full development of the early diagenetic network, thus also accounting for biogeochemical processes (e.g. methanogenesis) in deeper sediment layers that potentially affect biogeochemical cycling in the shallower sediment. Observations at the site indicate the absence of active flow and the advective velocity v_up is thus set to zero. Upper boundary conditions and model parameters are constrained on the basis of the observations reported (Brüchert, 2020) (Table S4). The observed organic carbon profile is imposed in the first 19 cm (Table S5) and organic carbon contents in deeper sediments are calculated on the basis of the reactive continuum model for organic matter degradation (described in Sects. S2 and S3) and the deepest observed value. In addition, the possibility of a source of methane is implemented at the bottom of the modelled sediment column by applying a Dirichlet boundary condition, thus taking into account the possible presence of methane seeping from deep sediments as results of destabilizing gas hydrates/subsea permafrost – a distinguishing feature of the ESAS sediments. The methane boundary condition is determined by model fitting (see below).

When evaluating model performance, particular attention is given to sulfate, methane, ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$), phosphates (${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$) and dissolved inorganic carbon (DIC) depth profiles. While the former two species are of main interest for evaluating simulated AOM dynamics, the remaining three serve as indicators for OM degradation dynamics since they are metabolic byproducts of degradation (see Table S2). Moreover, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ is only affected by nitrification (R7) and adsorption (R28). The latter, although important, acts homogeneously throughout the sediment (considering the slight variation in sediment porosity, LaRowe et al., 2017). It can thus only cause uniform shifts in the $\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]$ profile, but does not affect the overall shape of the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ depth profile. Similarly, ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ is only consumed by fluorapatite precipitation (R22) and adsorption processes (R29 and R31). Fluorapatite precipitation controls maximum dissolved ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ concentrations, while the mineral adsorption process (R29) exerts a homogeneous influence and the interaction with Fe(OH)₃ is expected to be minor and mainly affects ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ within the iron reduction zone.

OM reactivity parameters (a and ν), bottom methane concentration ([CH₄]₋) and reaction rates are varied to find the best fit between observed and simulated profiles. Methane concentrations at the bottom of the model domain can also exceed the saturation concentration $[{\mathrm{CH}}_{\mathrm{4}}{]}^{*}=\mathrm{14}$ mM (estimated according to the value reported in Dale et al., 2008a) to include the possibility of methane in gaseous form.

2.3 Modelling strategy

2.3.1 Steady-state sensitivity analysis

To evaluate the main physical and biogeochemical controls on the efficiency of the AOM biofilter and its impact on non-turbulent methane emission from deep methane sources such as dissociating permafrost and/or disintegrating methane gas hydrates in ESAS sediments, we conduct a comprehensive, steady-state sensitivity study. For this purpose, we design a set of two baseline scenarios that are broadly representative for environmental conditions encountered on the shallow ESAS:

1. a passive case, i.e. v_up=0 cm yr⁻¹;

2. an active case, i.e. with v_up=1 cm yr⁻¹, a value which falls within the range of fluid flow velocities $v_{\mathrm{up}}=\mathrm{0.005}-\mathrm{30}$ cm yr⁻¹ observed across a wide range of different active environments (Regnier et al., 2011).

For both baseline scenarios, we assume a water depth of 30 m, which is slightly shallower than the average water depth of the ESAS ∼45 m (James et al., 2016), since we are here interested in the shallow, near-coastal part of the shelf that potentially hosts large subsea permafrost reservoirs and is most affected by the warming. Temperature is set equal to 0^∘ C and is thus similar to the yearly average of −0.79^∘ C observed in the Laptev Sea at a depth of about 30 m (Dmitrenko et al., 2011). The bioturbation coefficients $D_{\mathrm{b}}^{\mathrm{0}}$ and bioirrigation coefficients α₀ (Thullner et al., 2009) are then derived from global empirical relationships according to Middelburg et al. (1997) and Thullner et al. (2009), respectively. The methane saturation concentration [CH₄]^* is calculated on the basis of the relationship proposed by Dale et al. (2008a) assuming a soil matrix density of 2.41 g cm⁻³. Values of φ₀ and c₀ (see Eq. 4) are determined based on LaRowe et al. (2017). Boundary conditions are reported in Table S7 and informed by observations. They are chosen to be broadly representative of the wider Siberian shelf environment.

Each sensitivity study run is forced with a range of different dissolved [CH₄] concentrations at the lower model boundary, mimicking different methane fluxes from thawing subsea permafrost and/or disintegrating methane gas hydrates at depth. The applied set of methane concentrations at the lower boundary range from zero to the methane gas saturation concentration $[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=$ 0–20–100–330–1169–5455 µM and also include the highest methane concentration that has been to date observed in ESAS cores (Overduin et al., 2015) ($[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{1.169}$ mM).

Table 1 and Table S6 summarize the parameters applied in the baseline simulation and Table S7 provides an overview of the applied upper boundary conditions.

Burwicz et al. (2011)Thullner et al. (2009)Regnier et al. (2011)

Table 1 Model parameters changed in the “one-at-time” sensitivity studies. Reported values are for the baseline simulations.

To assess the influence of environmental conditions on the efficiency of the AOM biofilter and its influence on non-turbulent methane emission from dissociating permafrost and/or disintegrating methane gas hydrates in ESAS sediments a set of five “one-at-time” parameter variation experiments is designed. It encompasses the most important controls on benthic methane cycling (Regnier et al., 2011; Meister et al., 2013; Egger et al., 2018) and parameter variation experiments are performed for both the passive as well as active baseline scenarios.

1. Sedimentation rate ω. The sedimentation rate is varied over 2 orders of magnitude (0.03–0.123–0.17–1.5 cm yr⁻¹). Maximum values are comparable to terrestrial sediment accumulation rates in the Lena River delta (Bolshiyanov et al., 2015), fast marine sedimentation rates during the early Holocene sea transgression (Bauch et al., 2001) and marine accumulation on subsea permafrost deposit in Buor Khaya Bay (∼1.1 cm yr⁻¹, inferred from Overduin et al., 2015), while minimum values are representative of sedimentation rates found in the East Siberian Arctic Sea (Stein et al., 2001 in Levitan and Lavrushin, 2009). The baseline value of ω is calculated based on the empirical global relationship proposed by Burwicz et al. (2011).

2. Active fluid flow v_up. Buoyancy-induced motion (Baker and Osterkamp, 1988), water streams channeled through fault lines or groundwater discharge (Charkin et al., 2017) can cause active fluid flow in Arctic shelf sediments underlain by subsea permafrost or gas hydrates (Judd and Hovland, 2009; Semenov et al., 2019). Therefore, v_up is varied from 0–0.3–0.5–1–3–7–10 cm yr⁻¹. This interval falls in the range of reported upward advective water velocities in marine sediments 0.005–30 cm yr⁻¹ (Regnier et al., 2011).

3. AOM constant k_AOM. Rate constants implicitly account for factors that are not explicitly described in the model and thus tend to show a strong variability between sites. A comprehensive compilation of published model AOM rate constants (Regnier et al., 2011) reveals a variability of over 6 orders of magnitude (10–10⁷ M⁻¹ yr⁻¹). The AOM rate constant k_AOM (Eq. 5) is thus varied over the range $k_{\mathrm{AOM}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{2}}$–5×10³–5×10⁴–5×10⁵–5×10⁶–5×10⁷ M⁻¹ yr⁻¹.

4. Organic matter reactivity (i.e. RCM parameters). Although the apparent OM reactivity is controlled by a combination of two parameters (a and ν), previous studies indicate a less pronounced variability in ν (Arndt et al., 2013; Sales de Freitas, 2018), as well as a strong control of a on the sulfate–methane transition zone (SMTZ) depth (Regnier et al., 2011; Meister et al., 2013). Thus, ν was kept constant, while a was varied over the entire range of previously published values a=0.1–1–10–100–500–1000 years (Arndt et al., 2013). Studies about ESAS organic matter degradation show a reactivity of deposited organic matter which is compatible with the RCM parameters we explored. For instance, Bröder et al. (2016) found a half life for the organic carbon in the East Siberian Arctic Shelf of 19–27 years, which would correspond to a=3.4–4.8 years, with ν=0.125.

5. Bioirrigation coefficient α₀. Bioirrigation activity remains largely unconstrained on the Siberian shelf due to the scarcity of observational data (Teal et al., 2008). However, environmental stressors, such as ice scouring (e.g. Shakhova et al., 2017, and references therein) and trawling, which can dig furrows up to few metres (Shakhova et al., 2017) are detrimental to the local fauna, thus suggesting a low bioirrigation intensity. However, observations from other polar sites indicate that although biological diversity and activity is often low, it might be locally enhanced (Clough et al., 1997). In addition, ice scouring might also enhance non-local transport seasonally. We therefore varied α₀ over the entire range of plausible values: 0–33–66–99.5–120–240 yr⁻¹ (Thullner et al., 2009).

2.3.2 Transient sensitivity study

Dale et al. (2008c) showed that temporally varying environmental conditions may reduce the efficiency of the benthic AOM filter and facilitate methane escape due to the delayed response of the microbial community to changing conditions. Therefore, in addition to the steady-state sensitivity study, we also performed a series of transient simulations to explore the impacts of seasonal and projected climate change on benthic methane effluxes on the ESAS in response to changing upward methane fluxes from dissociating permafrost and/or disintegrating methane gas hydrates. Transient simulations were run with a bioenergetic rate law for AOM (Eq. 6) and an explicit description of AOM biomass. Simulation results from the passive steady-state baseline run with $[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{0}$ mM were used as initial conditions for the transient experiments. Four different transient environmental perturbation scenarios that reflect seasonal (1, 2) as well as idealized future (3, 4) environmental variability on the ESAS were run with three different values of v_up=0–1–5 cm yr⁻¹ over a period of 200 years.

1. Seasonal CH₄: seasonal change in methane supply from permafrost thaw and/or hydrate destabilization. CH₄ concentration at the bottom of the sediment column: null for 6 months, then increasing up to a peak of [CH₄]₋ (20–100–330–1169–5455 µM) for the remaining 6 months of the year and again back to null concentration.

2. Seasonal ${\mathrm{CH}}_{\mathrm{4}}+{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$: seasonal freshening of waters due to riverine discharge and sea ice melt. During winter, higher bottom salinity (Dmitrenko et al., 2011) results in higher sulfate concentration (Dickson and Goyet, 1994), while lower salinities and thus sulfate concentrations characterize the melt season. The bottom boundary condition for methane [CH₄]₋ follows an opposite trend: it is set to zero during the winter months and increases in Arctic summer.

3. Linear CH₄: slow increase in methane supply from permafrost thaw and/or hydrate destabilization. A linear increase in the bottom boundary methane concentration [CH₄]₋ (from 0 up to the peak) over 200 years is applied.

4. Sudden CH₄: abrupt increase in methane supply from permafrost thaw and/or hydrate destabilization. An instantaneous change in bottom boundary methane concentration – from 0 to one of the peak value [CH₄]₋ – is applied.

2.3.3 Analysed output

For each simulation we evaluate the effect of the respective parameter change on

1. the non-turbulent (i.e. not-ebullition-driven) flux of methane from the sediments into the water column;

2. the depth of the SMTZ;

3. the efficiency (η) of the AOM biofilter (see Appendix A for the exact definition of AOM applied here).

In addition, fluxes of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and CH₄ at the SMTZ, the maximum and integrated AOM rate and the Damköhler number (Da) for AOM and methanogenesis are also calculated. The Damköhler number is defined as Eq. (B4) (see Appendix B) and sets the ratio between the typical transport timescale and the typical reaction timescale. If D_a<1, the reaction timescale is longer than transport timescale (i.e. the reaction is slower) and the process is reaction-limited. If D_a>1 the process is transport-limited. Finally, for transient simulations, the integrated AOM-performing biomass (ΣB) was also analysed.

Figure 1 Pore water concentration profiles for O₂, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, CH₄, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ and DIC at site 14-3 on the Laptev Sea (76.171^∘ N, 129.333^∘ E, 56 m water depth). Dots represents the measurements and continuous lines the simulated results. The boundary conditions and model parameters employed in the model are reported in Table S4, the measured organic carbon content in Table S5. For O₂ no measured profile is available.

3 Results and discussion

3.1 Case study: sediment core on the Laptev Sea shelf

Figure 1 compares simulated and observed depth profiles for site 14-3. Cores were retrieved during the SWERUS-C3 campaign (Miller et al., 2017; Brüchert et al., 2018; Brüchert, 2020). Simulation results show an overall good agreement with measurements, but also reveal a slight overestimation of ${\mathrm{NH}}_{\mathrm{4}}^{+}$. Data–model fitting reveals that reconciling simulated and observed CH₄ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ depth profiles requires a diffusive flux of CH₄ through the lower model boundary (i.e. a bottom boundary concentration of $[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{16}$ mM). Neither higher marine OM contents in sediment layers below the first 22 cm for which observations are available, nor higher reactivities results, can satisfactorily reproduce the observed sulfate depletion and observed gradients. Model–data fitting thus not only highlights the important role of AOM in controlling the SMTZ, but also indicates that upward migrating methane from deep, pre-Holocene sources, such as subsea permafrost in the sediment might be an ubiquitous feature on the Siberian shelf.

The simulated ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ and DIC profiles are in good agreement with data, suggesting that the degradation dynamics of marine organic matter and adsorption are well captured by the model – although the maximum concentration of ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ at depth is mostly controlled by the saturation value of $\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]$. The largest discrepancy between data and modelling results is observed for ${\mathrm{NH}}_{\mathrm{4}}^{+}$. Observed ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations first increase to a maximum at about 6 cm depth and the slightly decrease in the lower sediment layers, whereas simulated ${\mathrm{NH}}_{\mathrm{4}}^{+}$ show an asymptotic increase in ${\mathrm{NH}}_{\mathrm{4}}^{+}$ concentrations. The observed ${\mathrm{NH}}_{\mathrm{4}}^{+}$ profile might either indicate changes in OM reactivity and/or characteristics or spatially heterogeneous adsorption/desorption dynamics. Such heterogeneity is not incorporated into the model, and accounting for it would require additional information.

3.2 Main physical and biogeochemical controls on potential non-turbulent methane flux from ESAS sediments

3.2.1 General patterns of methane and sulfate cycling on the ESAS

The comprehensive ensemble of all sensitivity experiments allows exploration of the general patterns of methane and sulfate cycling under a range of environmental conditions that is broadly representative for conditions encountered on the ESAS at present (Fig. 2). Model results confirm that AOM is an efficient sink for the diffusive CH₄ supply from below. For most of the investigated environmental conditions (95 % of the runs), 95 %–99.9 % of the upward diffusing CH₄ is consumed within the SMTZ, resulting in very small or negligible methane effluxes ($\le {\mathrm{10}}^{-\mathrm{2}}$ µmol CH₄ cm⁻² yr⁻¹) from the sediment. If upscaled to the total area of the ESAS ($\sim \mathrm{1.485}\times {\mathrm{10}}^{\mathrm{6}}$ km², Wegener et al., 2015), for which methane outgassing estimates have been published, the smallest simulated non-turbulent methane flux (i.e. $\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{13}}$ µmol cm⁻² yr⁻¹, Fig. 2b) would sum up to a total flux of 2.1 mmol CH₄ yr⁻¹, resulting in a negligible role of non-turbulent, benthic methane fluxes to the Arctic methane budget.

However, results also show that, under a specific set of environmental conditions that lower the efficiency of the AOM biofilter (see detailed discussion below), non-turbulent CH₄ escape from ESAS sediments can reach values of up to 27 µmol CH₄ cm⁻² yr⁻¹. Simulation analysis shows that these high effluxes and, thus, low AOM biofilter efficiencies are generally obtained for environmental conditions that cause a shallow location of the SMTZ (<18 cm) and that they are very sensitive to changes in environmental conditions that would cause a deepening of the SMTZ. For instance, a deepening of the SMTZ from 18 to 26 cm results in a rapid increase in AOM efficiency from 1 % to 98 % (Fig. 2a). Furthermore, results indicate that, for SMTZ depths larger than 26 cm, AOM remains an efficient barrier across the full spectrum of investigated environmental conditions (Fig. 2). The observed link between AOM filter efficiency and SMTZ is reflected in the strong (semilog) linear relationship between methane flux at the SWI and the SMTZ depth (Fig. 2b). Such a relationship reveals the pivotal connections between these two quantities and mirrors the empirically found linear log–log relationship between measured CH₄ fluxes at the SMTZ and the SMTZ depths (Fig. S5 in the Supplement) by Egger et al. (2018). Maximum simulated CH₄ effluxes are thus comparable in magnitude to fluxes reported from other settings potentially sensible for CH₄ emissions. These include mud volcanoes, e.g. in the Gulf of Cadiz: 2.1–40.7 µmol CH₄ cm⁻² yr⁻¹ (Niemann et al., 2006a); Håkon Mosby mud volcano in the Barents Sea: 0.03 µmol CH₄ cm⁻² yr⁻¹ (Niemann et al., 2006b) and coastal settings, e.g. a Dutch coastal reservoir (20–80 µmol CH₄ cm⁻² yr⁻¹, Egger et al. (2016)) or tidal flats (4–800 µmol CH₄ cm⁻² yr⁻¹ Borges and Abril, 2011). Upscaling the highest simulated non-turbulent flux (27.48 µmol CH₄ cm⁻² yr⁻¹) to the ESAS results in a total efflux of 0.408 Tmol CH₄ yr${}^{-\mathrm{1}}=\mathrm{6.52}$ Tg CH₄ yr⁻¹. This value represents an estimated upper limit which, for comparison, equals ∼10 % of global marine seepage at seabed level (Saunois et al., 2016) and is in magnitude similar to the global methane efflux that has been estimated for upper continental slope sediments on a centennial timescale (4.73 Tg CH₄ yr⁻¹, Kretschmer et al., 2015).

Figure 2 Aggregation of all the simulation performed for the “one-at-time” sensitivity study. (a) AOM efficiency versus the depth of the SMTZ. (b) Scatter plot and semi-log fit of the methane flux ($J_{{\mathrm{CH}}_{\mathrm{4}}}$) at the SWI versus SMTZ depth.

Further insights into the general drivers that control methane dynamics in ESAS sediments are provided by Damköhler numbers. Damköhler numbers for simulated methanogenesis (Da_MG) and AOM (Da_AOM) are reported in Fig. S2. Da_MG (purple circles) are <1, span a range of ∼0.0021–0.43 and are thus comparable to previously reported Da_MG of 0.22 for methane gas hydrate bearing sites, such as the Hydrate Ridge and Kithley Canyon (Chatterjee et al., 2011). They reveal that methanogenesis is always slower than methane transport and that CH₄ dynamics driven by methanogenesis are thus reaction-limited. This result is consistent with the fact that methanogenesis rates are merely supported by the slow influx and transport of OM by burial and bioturbation.

In contrast, high Da_AOM values (Da_AOM=32–2.78×10⁵ – Fig. S2, orange circles), show that AOM is transport-limited, suggesting a sensitive role of transport parameters in determining AOM efficiency and in controlling methane flux across the SMTZ and subsequently the SWI.

3.2.2 Environmental controls and mechanisms of methane escape from ESAS sediments

The simulated general patterns of methane and sulfate cycling on the ESAS broadly corroborate previous findings regarding the dominant environmental controls on AOM biofilter efficiency and SMTZ depth (Regnier et al., 2011; Egger et al., 2018; Meister et al., 2013; Winkel et al., 2018). However, they also challenge intuitive views on the factors that favour high CH₄ escape through the SWI. In particular, they highlight the essential link between AOM efficiency and SMTZ depth and the central importance of environmental conditions that control the depth of the SMTZ. In addition, they suggest that transport processes play a dominant role for non-turbulent methane effluxes from ESAS sediments. The following sections explore the role of each of the investigated environmental conditions on methane efflux in more detail and also shed light on the mechanisms behind non-turbulent methane escape from ESAS sediments.

Figure 3 (a) Barplot of the methane flux at the SWI versus ω for passive case (plain style) and active case (pattern style) and the [CH₄]₋ reported in the text. The squared value of ω is the reference value. (b) Semilog plot of methane flux at SWI versus v_up for the different [CH₄]₋ reported in the text. (c) Log–log plot of methane efflux at SWI versus ω for passive case (diamonds) and active case (circle). The log–log fit is also displayed. (d) Log–log plot of SMTZ depth versus ω for passive case (diamonds) and active case (circle) with log–log fit. The red line is the trend found by Egger et al. (2018) (the term log(100) is to take into account unit conversion).

3.2.3 Role of advective transport

Figure 3a illustrates the effects of sedimentation rate ω on the flux of methane across the SWI. For both active (v_up=1 cm yr⁻¹) and passive (v_up=0 cm yr⁻¹) settings, simulated CH₄ effluxes increase as a power with sedimentation rate (log–log linear; see Fig. 3c) from $\mathrm{5.5}\times {\mathrm{10}}^{-\mathrm{15}}$ µmol CH₄ cm⁻² yr⁻¹ for low sedimentation rates (ω=0.03 cm yr⁻¹) to values as high as 27.5 µmol CH₄ cm⁻² yr⁻¹ for high sedimentation rates (ω=1.5 cm yr⁻¹). Accordingly AOM acts as an efficient filter for upward diffusing methane (with η∼100 %, see Fig. S4), in slowly accumulating sediments.

In contrast, the efficiency of the AOM biofilter drops to 50 %–0 % for high sedimentation rates. The main driver behind the simulated high CH₄ fluxes and low AOM efficiencies in these rapidly accumulating sediments, are enhanced methanogenesis rates. High sedimentation rates not only facilitate the supply of organic matter to the methanogenic zone of the sediment, but also reduce residence times in the upper sediment layer, resulting in a lower OM age (see Eqs. S14 and S16 in the Supplement) or degradation state (see Eq. S12) within the methanogenic zone. The enhanced supply of reactive OM to anoxic sediment layers supports higher methanogenesis rates, resulting in higher methane porewater concentrations and an upward shift of the SMTZ.

In addition, the presence of active fluid flow further enhances methane efflux. The CH₄ fluxes from below adds complexity to the overall methane dynamics and this effect is investigated further by contrasting Damköhler numbers for passive and active settings on the shelf. Table S8 shows that for low to intermediate sedimentation rates, Da_AOM values significantly decrease with v_up, indicating that less and less methane consumption occurs within the typical transport timescale τ_T, thus leading to a reduction in AOM biofilter efficiency.

Maximum simulated flux differences between active and passive settings can reach up to 10 orders of magnitude. However, flux differences quickly decrease with increasing sedimentation rates. Rapidly accumulating sediments show almost no difference in efflux between active and passive sites (Fig. 3a). In contrast to sedimentation rates, the mechanism behind the control of v_up on non-turbulent methane efflux is straightforward and self-evident. Active flow enhances the upward transport of CH₄, shifting the SMTZ upwards and, thus, increasing CH₄ concentrations at shallow sediment depths (see Fig. 3d). The apparent paradox of the CH₄ efflux insensitive to fluid flow in fast accumulating sediments can be resolved by examining the dissolved CH₄ depth profiles (Fig. 4). Simulated depth profiles are nearly identical and reveal CH₄ concentrations at or near the saturation concentration. In fast accumulating sediments, high methanogenesis rates result in an over-saturation of porewaters directly below the generally shallow SMTZ. High methanogenesis rates thus support the build-up of methane gas. Methane gas formation also explains why, in these cases, integrated methanogenesis exceeds non-turbulent CH₄ fluxes by up to 6 times. In rapidly accumulating, active and passive sediments, non-turbulent CH₄ fluxes are thus essentially identical. However, active settings will be characterized by the additional build-up of gaseous CH₄ and its potential escape through the sediment–water interface – a process not simulated in the present study.

Figure 4 Porewater profiles in the case of ω=1.5 cm yr⁻¹ for CH₄ (a), ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (b) and gaseous CH₄ (c). Dashed lines are simulation in passive scenario with $[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{0}$ mM, while continuous lines simulations display active scenario with $[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{5.455}$ mM, corresponding to the saturation concentration in the environmental conditions considered for the representative profile.

Model results thus show that the dominant mechanism behind the observed transport-control on non-turbulent CH₄ efflux is an overall increase in CH₄ concentration and an upward shift of the SMTZ rather than an increasing relative contribution of advective transport processes to the total efflux. In fact, a comparison of the different methane transport processes across the SWI (i.e. molecular diffusion, bioturbation-induced diffusion, bioirrigation and advective transport, as explained in Sect. S1 in the Supplement) shows that the relative contribution of both the advection and molecular diffusion flux to the total flux is small and further decreases with increasing v_up (Fig. S3). High non-turbulent methane effluxes in rapidly accumulating and/or active settings are thus largely driven by the non-local irrigation flux (see Sect. 3.2.5 for more details on the role of irrigation). With increasing ω or v_up, the SMTZ shifts upwards, resulting in higher methane concentrations at shallow sediment depths and thereby reinforcing the relative contribution of non-local transport for CH₄ fluxes, as well as lowering the efficiency of the AOM barrier from η∼100 % to η∼78 %. The important role of the SMTZ location as a key control on CH₄ efflux is further confirmed by the observed power-law relationship between the location of the SMTZ and ω (Fig. 3d). This result is qualitatively in agreement with the global compilation of empirical data by Egger et al. (2018), which reveals the same log–log decreasing trend between SMTZ and sedimentation rate. Our results are also consistent with observations from brackish sediments that show that sedimentation rates >10 cm yr⁻¹ give rise to high non-turbulent CH₄ fluxes (20–80 µmol CH₄ cm⁻² yr⁻¹) and a high OM burial efficiency (∼78 %, Egger et al., 2016). Egger and co-workers explained these findings by the slow growth of AOM micro-organisms and the resulting inability of the microbial community to consume all of the CH₄ produced. However, our results show that the same pattern can be observed without having to invoke a low AOM efficiency. Our simulations thus indicate that the rapid burial of reactive organic matter to deeper sediment layers in rapidly accumulating sediments is sufficient to explain high CH₄ effluxes.

3.2.4 Role of organic matter quality

The quality of organic matter deposited onto the sediment exerts an additional control on CH₄ efflux. Figure 5 illustrates the influence of organic matter quality (as a function of OM degradation model parameter a; see Eq. S12) and sedimentation rate ω on non-turbulent methane efflux for both active and passive settings, as well as different methane fluxes from below. Results corroborate the dominant influence of sedimentation rates on methane efflux, while organic matter quality exerts a secondary control. This also means that, in order to assess the main features of possible CH₄ efflux in terms of modelling, capturing the details of organic matter quality is not fundamental. Maximum fluxes are generally simulated for rapidly accumulating sediments ω>0.5 cm yr⁻¹ that receive organic matter of intermediate quality (a=10–100 years).

Figure 5 Flux of methane at the SWI as dependent on a and ω. For $[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{0}$ mM (a, c) and $[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{5.455}$ mM (b, d), and passive (a, b) and active (c, d) case. The circle with pattern corresponds to the baseline simulation.

These findings are in agreement with previously published studies (Regnier et al., 2011; Meister et al., 2013) and can be explained with the fact that high methanogenesis rates require a supply of reactive OM to the methanogenic zone. If organic matter quality is high (a<10 years), methanogenesis becomes substrate limited due to the rapid degradation of organic matter through energetically more favourable degradation pathways in the shallow sediments. In turn, if organic matter quality is low (a>100 years), methanogenesis becomes reactivity limited. And this is especially true for ESAS sediments at depth. Therefore, the ideal combinations of organic matter reactivity and sedimentation rate that result in maximum methane effluxes correspond to conditions characterized by OM that is (i) sufficiently reactive to support enhanced methanogenesis rates and thus an accumulation of CH₄ at depth, but (ii) sufficiently unreactive (in comparison to the burial rate) to escape the complete degradation in non-methanogenic sediments. BRNS outcomes show that the onset of active fluid flow and an enhanced methane supply from below (i.e. higher CH₄ concentration at the lower boundary) increase the CH₄ efflux through the SWI without altering the overall patterns (see Fig. 5a–b versus 5c–d).

3.2.5 Role of non-local transport

The analysis of the influence of bioirrigation on non-turbulent CH₄ efflux from the ESAS (Fig. S12) shows that such a process enhances methane efflux in sediments that are characterized by a shallow SMTZ. However, bioirrigation exerts a limited effect under a range of environmental conditions that favour a deep or shallow SMTZ location, respectively.

In passive settings, changes in bioirrigation coefficient (α₀) exert a limited control on CH₄ effluxes. For most passive scenarios, the SMTZ is located well below the sediment layer affected by bioirrigation (z_irr=3.5 cm, which implies that bioirrigation is strongly suppressed below 15 cm) and, thus, changes in α₀ have no effect on methane efflux.

In contrast to passive settings, active settings reveal a rapid increase in methane efflux with the onset of bioirrigation activity. Methane effluxes first increase by up to 5 orders of magnitude from α₀=0 yr⁻¹ to α₀=5 yr⁻¹, reaching maximum effluxes of ∼0.02 µmol CH₄ cm⁻² yr⁻¹, before remaining almost constant with a further increase in bioirrigation coefficients (up to 240 yr⁻¹). The simulated increase in methane efflux is a direct effect of the transport process itself, which enhances the upward transport of methane accumulating in the upper sediment layers, including layers below the generally shallow SMTZ.

These results are corroborated by the concomitant analysis of CH₄ dynamics over the three-dimensional transport coefficient ω, v_up and α₀ space shown in Fig. 6.

Figure 6 Efflux of methane at the SWI as dependent on v_up and ω for α₀=0 (a), α₀=5 (b), α₀=10 (c) and α₀=33 yr⁻¹ (d). Circles represent simulations outcomes. Results for α₀≠0 yr⁻¹ are almost the same. The lower boundary condition for methane is $[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{1.169}$ mM. A three-dimensional version of the plot is reported in Fig. S13.

A comparison between simulations with α₀=0 yr⁻¹ and α₀≠0 yr⁻¹ (α₀=5, α₀=10 and α₀=33 yr⁻¹) shows that irrigation increases the CH₄ efflux at low to intermediate sedimentation rates and/or high v_up (lower-left corner of the phase space in both plots). However, maximum methane effluxes that are simulated for high sedimentation rates or v_up are almost identical between bioirrigated and non-irrigated sites despite the differences in dominant transport mechanism (diffusion when α₀=0 yr⁻¹; irrigation when α₀≠0 yr⁻¹). Under these conditions (i.e. high v_up and/or high ω), the SMTZ is located close to the SWI. In such settings, non-local transport becomes the dominant transport process in bioirrigated sediments (see Sect. 3.2.3 and Fig. S3) because it weakens concentration gradients near the SWI and, thus, contributes to a substantial reduction in the gradient-driven, diffusive transport terms. As a consequence, simulated CH₄ efflux at the SWI are broadly similar for all of the investigated α₀≠0 yr⁻¹ (Fig. 6b, c, d). It is worth noticing that, independently on the α₀, CH₄ efflux for ω=0.03 and v_up=10 cm yr⁻¹ is ∼1 µmol CH₄ cm⁻² yr⁻¹ – a value almost identical to the one reported in Luff and Wallmann (2003) – 1.4 µmol CH₄ cm⁻² yr⁻¹ – for sediments characterized by v_up=10 and ω=0.0275 cm yr⁻¹.

3.2.6 AOM rate constant

Given its crucial role in AOM biogeochemistry, one would expect a pronounced influence of the kinetic rates constant, k_AOM, on non-turbulent methane effluxes. However, modelling reveals that k_AOM only plays a minor role for non-turbulent methane fluxes across the SWI (see Figs. S14, S15). An increase in k_AOM can reduce methane effluxes from passive shelf sediments by up to 5 orders of magnitude. Still, its effect remains small compared, for instance, to the response to variations in sedimentation rate, which can change methane efflux by up to 14 orders of magnitude. The most important effect of increasing k_AOM is the increasing linearity of the [CH₄] and $\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]$ profiles around the SMTZ and the concurrent narrowing and downward movement of the SMTZ, which can result in a reduction in methane efflux. Simulations thus show that the AOM biofilter and, as a consequence, non-turbulent methane effluxes from sediments, are not affected by the exact value of the kinetic rate constant, at least in the range we analysed. This is in disagreement with results by Dale et al. (2008c), which show that, in dynamic settings subject to large methane fluxes, an increase of 3 orders of magnitude in k_AOM (from 10² to 10⁵ M⁻¹ yr⁻¹) leads to a reduction in steady-state methane fluxes below 10⁻² µmol CH₄ cm⁻² yr⁻¹. However, this discrepancy might be ascribable to the high water flow velocity employed in their simulation (v_up=10 cm yr⁻¹), 10 times higher than the one we considered in our active simulations. Finally, on the shallow ESAS, dissolved methane concentrations are limited by the comparably low gas saturation concentration, resulting in a minor influence of k_AOM on methane fluxes (as the AOM rate is proportional to the CH₄ concentration). An indirect support to our findings regarding the secondary role of k_AOM on the AOM itself comes from Luff and Wallmann (2003). They showed that, as long as it is not null, the actual value of k_AOM is unimportant for the precipitation of authigenic carbonate. Since the authigenic carbonate precipitation is largely driven by alkalinity produced during AOM (e.g. Aloisi et al., 2004; Luff et al., 2005; Karaca et al., 2010; Pierre et al., 2012; Crémière et al., 2016b, a; Meister et al., 2018), the observed independence of precipitation rates from k_AOM bolsters our conclusion.

3.2.7 Summary of steady-state experiments

The results of the steady-state sensitivity study indicate that, under environmental conditions that are broadly representative for the ESAS, low AOM efficiencies and thus high non-turbulent CH₄ effluxes from thawing subsea permafrost and/or dissociating methane gas hydrates (larger than 4 µmol CH₄ cm⁻² yr⁻¹) are promoted by intense advective transport (sedimentation rate ω>1 cm yr⁻¹, active fluid flow v_up>7 cm yr⁻¹). Under these conditions, CH₄ efflux can be further enhanced by moderate OM reactivity (a=10–10² years) and intense non-local transport processes, such as bioirrigation (irrigation constant α₀>0 yr⁻¹). Overall, non-turbulent benthic escape of CH₄ from deep sources appears to be mainly controlled by the concurrent effects of ω, v_up and α₀. In contrast, maximum AOM rates, k_AOM, exert no influence on the AOM filter efficiency.

3.2.8 Geographic pattern and potential for non-turbulent methane emissions from Laptev Sea sediments

One strength of models is that they can provide the explorative means to assess dynamics at spatial/temporal scales that cannot easily be assessed by observations alone. In particular, transfer functions, simple look-up tables and neural networks that are derived from, or trained on, a large ensemble of individual model simulations over a broad range of plausible boundary conditions have been frequently and successfully used to investigate regional and even global dynamics (Gypens et al., 2008; Marquardt et al., 2010; Dale et al., 2015, 2017; Capet et al., 2016; Bowles et al., 2014). Such a quantitative framework for first-order estimates of potential non-turbulent methane escape from ESAS sediments can also be derived from the results of the sensitivity study we ran.

Model results indicate that sedimentation rate exerts the dominant control on benthic escape of methane from thawing subsea permafrost and/or dissociating methane gas hydrates on the ESAS. The functional relationship between sedimentation rate and methane flux across the SWI reported in Fig. 3c thus allows estimation of a potential non-turbulent, benthic methane efflux derived from deep sources for a given sedimentation rate. Thus, if the spatial distributions of these environmental controls on methane efflux are known, a first-order geographical distribution of potential non-turbulent methane escape from the Siberian shelf can be derived. However, the availability of observational data from the Siberian shelf is extremely scarce. Therefore, we here focus on the Laptev Sea – a comparable well-studied part of the Siberian shelf. The Laptev Sea is well known for its subsea permafrost and gas hydrate content and subject to large riverine inputs from the Lena River. To derive a map of sedimentation rates for Laptev Sea shelf sediments, we use published linear sedimentation rates (Table S9) and extrapolate these values to the entire region by applying a simple three-dimensional kriging method (see Fig. 7a), using the International Bathymetric Chart of the Arctic Ocean (IBCAO) (Jakobsson et al., 2012) and employing longitude, latitude and water depth as predictors for ω.

Observations indicate that sedimentation rates are highest (ω=0.45 cm yr⁻¹) close to the mouth of the Lena River and Moustakh Island in the Buor-Khaya Gulf. As a consequence, the vicinity of the river mouth, as well as the area along the shallow bathymetric profile towards the NE of the Lena delta are characterized by comparably high sedimentation rates (ω=0.27–0.42 cm yr⁻¹). The relatively shallow areas (∼10 m deep) around the New Siberian islands reveal intermediate values (ω=0.06–0.12), while minimum sedimentation rates (∼0.002–0.03 cm yr⁻¹) roughly follow the 55 m isobath down to the continental slope at 100 m. Deeper shelf areas are characterized by a more homogeneous distribution of sedimentation rates with values around 0.03–0.06 cm yr⁻¹.

Table 2 Estimated flux of CH₄ at SWI in mol yr⁻¹ for different depth regions of Laptev Sea in a passive (v_up=0 cm yr⁻¹) and active (v_up=1 cm yr⁻¹) case.

Figure 7 (a) Values of the sedimentation rate extrapolated for the whole Laptev Sea via a simple kriging method. The reference values (circles) are the ones reported in Table S9. (b, c) (Log) values of the potential non-turbulent dissolved methane emissions at the SWI considering the relationship presented in Fig. 3c for passive (b) and active (c) cases.

Estimated non-turbulent methane effluxes corresponding to the highest measured sedimentation rates close to the Lena mouth do not exceed $\mathrm{1.57}\times {\mathrm{10}}^{-\mathrm{1}}$ µmol CH₄ cm⁻² yr⁻¹ assuming the presence of active fluid flow and $\mathrm{2.25}\times {\mathrm{10}}^{-\mathrm{5}}$ µmol CH₄ cm⁻² yr⁻¹ for passive settings. These findings are not surprising as steady-state sensitivity results indicate that high CH₄ efflux requires sedimentation rates of ω>1 cm yr⁻¹. The regional non-turbulent CH₄ efflux budget for different depth sections of the Laptev Sea assuming the absence of active fluid flow in Laptev Sea shelf sediments (see Table 2) thus indicates that non-turbulent CH₄ efflux is negligible. Even if we assume the omnipresence of an active fluid flow of v_up=1 cm yr⁻¹, the estimated non-turbulent methane efflux merely sums up to 9.39×10⁶ mol CH₄ yr⁻¹ (∼0.1 Gg CH₄ yr⁻¹) over the entire Laptev Sea area of 527.4×10³ km². Such small effluxes would most likely be subject to further oxidation in the water column, thus limiting any potential impact on atmospheric methane concentrations and climate.

Higher advective fluid flow velocities, intermediate organic matter reactivity and/or a more intense macrobenthic biological activity could increase these estimates of non-turbulent methane escape from the Laptev Sea shelf. Higher advective fluid flow velocities (i.e. v_up>1 cm yr⁻¹), possibly in connection with active seepages, groundwater discharges and fault lines (the latter follow parallel pattern in Laptev Sea (Drachev et al., 1998) on the direction SW–NE from the west of Lena delta up to the little Lyakhovsky and Kotelny islands), could result in methane effluxes of up to 10–10^1.3 µmol CH₄ cm⁻² yr⁻¹ (see Figs. 5 and 6). However, such high fluid flow velocities would be only found locally and would thus merely give rise to a number of methane emission hotspots that would not change the overall non-turbulent methane flux budget. In addition, intermediate organic matter reactivity, in particular in the fast accumulating sediments close to the coastline and the Lena River delta that receive more reactive organic matter from thawing terrestrial permafrost (Wild, 2019) could result in a higher estimated non-turbulent methane escape. However, our sensitivity study results show that OM reactivity merely plays a secondary role, suggesting that changes in OM reactivity would only change efflux by less than an order of magnitude assuming both a=100 years or a=1 years. Changes in bioirrigation intensity would exert merely a limited effect on efflux estimates, as bioirrigation has already been included in the estimate calculations. The absence of bioirrigation, which is known to be patchy in Arctic sediments, could act both in the direction of further reducing (limiting the bioirrigated flux from the sediments) or increasing (by limiting the flux of TEAs from the seawater and therefore oxidation) the estimated non-turbulent methane efflux. Additional physical reworking, such as ice scouring or dredging, may also have such an opposite effect: it could reduce the methane efflux (by enhancing the flux of TEAs into the sediments) but it could also intensify it (by removal of the upper sediment layer).

Model results thus show that, under present-day, steady-state environmental conditions, AOM acts as an efficient biofilter for potential non-turbulent methane fluxes in Laptev Sea sediments. The estimated non-turbulent methane escape from Laptev Sea shelf sediments cannot support previously estimated methane outgassing fluxes of few teragrams of CH₄ yr⁻¹ (Berchet et al., 2016; Thornton et al., 2020) or even tens of teragrams of CH₄ yr⁻¹ (Shakhova et al., 2014). If such outgassing were to be supported by methane efflux from Laptev Sea sediments, it would require the build-up of CH₄ gas reservoirs in Laptev Sea sediments of at least similar or larger size than the evaded amount, as well as the preferential and rapid transport of this CH₄ gas to the atmosphere. Nevertheless, model results also suggest that projected trends of terrestrial permafrost thawing and coastal permafrost degradation (Vonk et al., 2012) might increase the importance of non-turbulent methane escape for the Arctic's methane budget by potentially increasing sedimentation rates through coastal erosion (vast amount of debris and terrigenous material) and increased riverine inputs (Guo et al., 2007); active fluid flow through permafrost and methane gas hydrate degradation (James et al., 2016; Ruppel and Kessler, 2017); organic matter reactivity through an enhanced delivery of more reactive permafrost organic matter (Wild et al., 2019) and/or an enhanced macrobenthic activity through warming and Atlantification. However, the magnitude of these projected environmental changes and thus their effect on non-turbulent methane escape from ESAS sediments is difficult to assess.

3.3 Methane efflux dynamics in response to seasonal and long-term environmental variability

The steady-state sensitivity results reveal that, under steady-state conditions, AOM represents an efficient biofilter for upward migrating methane from thawing permafrost and/or dissociating methane gas hydrates on the ESAS. However, transient dynamics induced by, for instance, seasonally or climate-change-driven variability in environmental conditions, may weaken the efficiency of the AOM biofilter. Therefore, we additionally explore the potential for non-turbulent methane escape from thawing subsea permafrost and/or dissociating methane gas hydrate in ESAS sediments under transient conditions. Table 3 summarizes the maximum simulated non-turbulent methane fluxes for two kinds of environmental change scenarios: seasonal and long-term. With the former, we explore seasonal changes in deep methane flux and seasonal freshening of bottom waters. With the latter instead, we investigate the impacts of a slow linear increase and a sudden maximum increase in deep methane flux (see Sect. 2.3.2).

Table 3 Maximum of methane fluxes (in µmol cm⁻² yr⁻¹) at SWI for the four analysed transient scenarios. Values in round parenthesis indicate the year after the beginning of simulation corresponding to the reported maximum.

Results reveal that the transient response of simulated non-turbulent methane efflux is similar for all environmental scenarios, but instead significantly differs for passive and active sites. In general, passive settings do not allow for significant methane escape (Fig. S17). Although transient methane efflux monotonously increases over the simulated period, it only reaches a maximum value of 0.03–0.05 µmol CH₄ cm⁻² yr⁻¹ (Fig. S17). Similarly, the simulated SMTZ depth merely migrates 11.5–29 cm upwards (Fig. S18). Over the simulated 200 years, the integrated non-turbulent methane escape from passive settings for all environmental change scenarios barely reaches 3–4 µmol CH₄ cm⁻².

In contrast, active settings (i.e. v_up=1 cm yr⁻¹) exhibit an initial increase in CH₄ fluxes to maxima of 0.55–0.83 µmol CH₄ cm⁻² yr⁻¹ over the first 50 years. This growth coincides with a rapid upward shift of the SMTZ by 100 cm. Methane escape then temporarily drops by 17 %–20 % until years 70–75, when it begins to increase again until the end of the simulation. During this second phase (i.e. after the first 50 years), the SMTZ remains stationary. Temporally integrated methane efflux (over 200 years) increases with active fluid flow rate from 66–121 µmol CH₄ cm⁻² for v_up=1 cm yr⁻¹ to $\sim \mathrm{0.95}-\mathrm{1.154}$ mmolCH₄ cm⁻² for v_up=5 cm yr⁻¹. A large fraction of these emissions (30 % and 48 %–87 %, respectively) occurs in the first 100 years after the perturbation.

Model results thus indicate that the exact temporal character of environmental changes does not exert an important influence on non-turbulent methane efflux. Conversely, both microbial growth dynamics and the presence/absence of active fluid flow (Table 3) largely control the transient response to environmental change. The reasons for this are 2-fold. First, the response time of the resident AOM community is longer than the characteristic timescales under investigation, thus smoothing out the impact of environmental perturbations. Second, active fluid flow enhances the impact of the perturbation by triggering a significant upward shift of the SMTZ. In particular, the initial movement of the SMTZ prevents the establishment of an efficient AOM community at the SMTZ: this creates a “window of opportunity” for methane escape. In contrast, the comparably slow and limited movements of the SMTZ in passive settings (Fig. S18) enables the establishment of an AOM community that acts as an efficient biofilter for upward-migrating methane.

The following sections explore the factors that control the creation of such a window of opportunity and discusses the mechanisms behind the simulated methane escape. Given the overall similar transient response of non-turbulent methane fluxes to different environmental scenarios (Figs. S17, S18), we will base the following discussion on scenario 4, namely a step-like CH₄ forcing with v_up=1 cm yr⁻¹ and a specific bottom concentration ($[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{1.169}$ mM). The reason for selecting this scenario is simple. In contrast to the other scenarios, scenario 4 allows for a straightforward definition of the initial and final state, which facilitates the attribution of a typical response timescale for the system.

3.3.1 Window of opportunity

Figure 8 illustrates the temporal evolution of the simulated (a) filter efficiency and AOM rate, (b) CH₄ efflux, (c) SMTZ depth and (d) AOM biomass for transient scenario 4 in the case of v_up=1 cm yr⁻¹ and $[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{1.169}$ mM. The onset of a sudden methane flux from thawing permafrost and/or dissociating methane gas hydrates below the sediment column triggers the rapid movement of the SMTZ. Simulation outputs show that the velocity at which the SMTZ moves upward (v_SMTZ) is solely controlled by v_up, as evident from the constant v_SMTZ∼2.46 cm yr⁻¹ for all the transient scenarios with v_up=1 cm yr⁻¹ (Fig. S18). The initial upwards movement of the SMTZ delays the microbial response since the transient dynamics inhibits the establishment of a resident AOM community sufficiently large to consume upward migrating methane. The AOM rate, and thus the filter efficiency, is controlled by the AOM biomass dynamics (Eq. 6), which in turn is determined by the kinetic (F_K, Eq. 7) and thermodynamic (F_T, Eq. 8) constraints. Figure 10 illustrates the typical depth profiles of the thermodynamic and kinetic terms in the bioenergetic AOM formulation (Eq. 6), as well as their evolution in response to the onset of a sudden methane flux from below. Initially, although kinetically possible (i.e. F_K≠0), AOM is inhibited by thermodynamic constraints (i.e. F_T=0). For v_up=1 cm yr⁻¹ during the first 23 years, AOM biomass thus remains largely constant (Figs. 8d and S21a) and, as a consequence, AOM rate and filter efficiency are zero. In this period, aerobic methane oxidation represents the only barrier to upward diffusing methane. However, this barrier is weak due to the limited availability of oxygen and the competition with aerobic organic matter degradation as well as additional secondary redox reactions that also consume oxygen (see Table S3). As a consequence, CH₄ efflux increases. The initial methane efflux is largely supported by in situ methanogenesis since the advective transport of methane (occurring at $v_{\mathrm{up}}-\mathit{\omega }=\mathrm{0.877}$ cm yr⁻¹, corresponding to 20.17 cm in 23 years) is too slow to allow methane from below 3 m to reach the sediment–water interface.

After the first 23 years, thermodynamic constraints ease and AOM begins to efficiently consume upward migrating methane at the SMTZ by 40 % (Fig. 8a). However, as consumption occurs at the SMTZ (for the specific case at a sediment depth of 100.4 cm), it does not immediately affect the methane efflux at the SWI. The time required for the consumption signal to propagate to the SWI with velocity $\stackrel{\mathrm{‾}}{v}=v_{\mathrm{SMTZ}}+v_{\mathrm{up}}-\mathit{\omega }=\mathrm{3.337}$ cm yr⁻¹ is therefore $\frac{\mathrm{100.4}\mathrm{cm}}{\mathrm{3.337}\mathrm{cm}{\mathrm{yr}}^{-\mathrm{1}}}=\mathrm{30.1}$ years. Consequently, methane efflux further increases. This methane efflux is now also supported from deep sources such as thawing permafrost and/or dissociating methane gas hydrates, which have started to contribute to methane efflux between years 7 and 20 (assuming typical values of v_up reported for active marine sediments of 0.5–5 cm yr⁻¹). Methane efflux typically peaks 2–3 decades after the onset of methane supply. Maximum methane efflux increases with v_up: from 0.5–0.6 µmol CH₄ cm⁻² yr⁻¹ for v_up=1 cm yr⁻¹ to 11–19 µmol CH₄ cm⁻² yr⁻¹ for v_up=5 cm yr⁻¹. However, the duration of this initial “window of opportunity” for methane escape decreases with increasing v_up. In general, simulated maximum methane fluxes fall within the range of previous models applied to different environments (Sommer et al., 2006; Dale et al., 2008c) but do not reach the high values measured in other settings (Linke et al., 2005; Regnier et al., 2011).

After the initial “window of opportunity” (i.e. $\mathrm{23}+\mathrm{30.1}=\mathrm{53.1}$ years), the effect of an efficient methane consumption at the SMTZ starts to reduce the non-turbulent methane efflux at the SWI (Fig. 8b). This reduction lasts until the upward movement of the SMTZ slows down. At this point, the AOM filter efficiency reaches a quasi-stationary level of ∼85 % (as Fig. 8a). Meanwhile, in situ methanogenesis continues to produce methane, which is not entirely consumed by the AOM community that already reached its full capacity. As a consequence, methane fluxes at SWI increase again until a new steady state is reached.

Figure 8 Time evolution over 200 years for the case of an active set-up with v_up=1 cm yr⁻¹ and a step-like methane forcing from below from 0 to $[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{1.169}$ mM. (a) AOM vertically integrated rate (blue) and AOM efficiency (red). (b) CH₄ flux at SWI. (c) SMTZ depth. (d) Vertically integrated biomass (number of cells).

Figure 9 Vertical profiles at the end of transient simulation (after 200 years) with bioenergetic AOM formulation for the case $[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{1.169}$ mM and v_up=5 cm yr⁻¹. (a) Bimolecular product $\left[{\mathrm{CH}}_{\mathrm{4}}\right]\cdot \left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]$. (b) AOM rate according to the bioenergetic formulation (blue) and, for comparison, according to bimolecular formulation used for the steady-state simulations (red). (c) Apparent k_AOM, estimated from Eq. (5).

Figure 10 Vertical profile of F_T, F_K, $F_{\mathrm{Tot}}=F_{\mathrm{K}}\cdot F_{\mathrm{T}}$ and the AOM (scaled to the maximum) for three instants in time: 8 years (a), 19 years (b) and 200 years (c) of simulation, for the case $[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{1.169}$ mM and v_up=5 cm yr⁻¹.

3.3.2 Final new steady state

The final new steady-state value of methane efflux (Fig. 8b, S17 and S19) is generally in good agreement with Dale et al. (2008c), who reported an efflux of the same order of magnitude (3 µmol CH₄ cm⁻² yr⁻¹) for the new steady state at the end of a transient run with v_up=10 cm yr⁻¹ and $[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{70}$ mM. Simulations with v_up=5, ω=0.123 cm yr⁻¹ and $[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{1.169}$ mM (Fig. S19) offer a better understanding of the model. In this case, the final new steady state is about 2 orders of magnitude larger than the efflux of ∼0.1 µmol CH₄ cm⁻² yr⁻¹ simulated in the steady-state simulations, with bimolecular rate law, under identical environmental conditions (inferred from Fig. 6).

The reason for this discrepancy can be clarified by plotting the apparent k_AOM for transient simulations. Such a value is calculated by computing an apparent bimolecular rate constant k_AOM (as in Eq. 5) from the transient bioenergetic simulations for the new final steady state. Results are shown in Fig. 9. Figure 9a illustrates that the concentration product $\left[{\mathrm{CH}}_{\mathrm{4}}\right]\cdot \left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]$ is wider than the AOM rate profile (Fig. 9b, blue curve). Figure 9c also shows that the apparent k_AOM is not uniform: it reaches a maximum value of 109 M⁻¹ yr⁻¹, but remains well below 100 M⁻¹ yr⁻¹ at most depths. Compared to the values typically applied for bimolecular rate laws (i.e. k_AOM=10²–10⁷ M⁻¹ yr⁻¹), these values are rather low and reflect the ongoing thermodynamic limitation of AOM. F_T remains the main constraint on AOM throughout the simulation (Fig. 10c). A more uniform sulfide concentration – [HS⁻] enters in defining F_T – in lower sediments combined with the upward movement of the SMTZ pushes the maximum of F_T upwards, thus limiting the zone where AOM is thermodynamically favourable (∼13 cm deep).

Integrated biomass ΣB ranges from $\sim \mathrm{1.2}\times {\mathrm{10}}^{\mathrm{10}}$ to 3.5×10¹¹ cells cm⁻² (except for simulation with v_up=5 cm yr⁻¹ and $[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{5.455}$ mM, whose ΣB$=\mathrm{1.2}\times {\mathrm{10}}^{\mathrm{12}}$). These values are comparable with AOM biomass reported in Treude et al. (2003) (1.5–1.8×10¹⁰ cells cm⁻²) or with values simulated in Dale et al. (2008c) (3.7×10¹¹ cells cm⁻² for v_up=5 cm yr⁻¹). In addition, the maximum simulated biomass for active settings (0.5–2.5×10¹⁰ cells cm⁻³) agrees well with previously reported values, ranging from 0.27 to 7.4×10¹⁰ cells cm⁻³ (Dale et al., 2008c). Integrated AOM rates (ΣAOM) are instead smaller then previously published rates for shallow, active sites above the shelf break (Boetius et al., 2000; Haese et al., 2003; Luff and Wallmann, 2003; Linke et al., 2005; Wallmann et al., 2006b; Dale et al., 2008c), but comparable to those observed in active sites below the shelf break (Aloisi et al., 2004; Wallmann et al., 2006a; Maher et al., 2006) or in passive settings (Borowski et al., 1996; Martens et al., 1998; Fossing et al., 2000; Jørgensen et al., 2001; Dale et al., 2008c). The discrepancy may be due to different environmental conditions encountered at these sites. For instance, Dale et al. (2008c) applied an advective velocity of v_up=10 cm yr⁻¹ and $[{\mathrm{CH}}_{\mathrm{4}}{]}_{-}=\mathrm{60}$ mM. While differences in v_up affect the ΣAOM, its effect on ΣB is negligible since an efficient AOM microbial filter has to account for at least >10¹⁰ cells cm⁻³ (Lösekann et al., 2007; Knittel and Boetius, 2009).

4 Conclusions

In this study, we evaluate the potential for non-turbulent, benthic methane escape from thawing subsea permafrost and/or dissociating methane gas hydrates in both passive as well as active settings and under a range of environmental conditions that are broadly representative for conditions encountered on the present and future East Siberian Arctic Shelf (ESAS). We identify the most important biogeochemical and physical controls on non-turbulent methane escape from those sediments under steady-state conditions, as well as in response to environmental variability on seasonal and centennial timescales. Based on model results, we derive a simple transfer function that allows establishment of a first-order regional estimate of (not-turbulent) methane efflux and of potential methane consumption in Laptev Sea sediments.

Model results reveal that AOM is an efficient sink for upward migrating, dissolved methane in ESAS sediments. Simulated non-turbulent methane effluxes are negligible for a broad range of environmental conditions under both steady-state and transient conditions. Since AOM is a transport-limited process, transport parameters exert a dominant control on the efficiency of the AOM biofilter and, ultimately, on the methane efflux at the SWI. Both steady-state and transient model results confirm the key role of advective transport (mainly sedimentation and active fluid flow) in supporting methane escape from Arctic shelf sediments. Under steady-state conditions, high methane effluxes (up to 27.5 µmol cm⁻² yr⁻¹) are generally found for sediments that are characterized by high sedimentation rates and/or active fluid flow (sedimentation rate ω>0.7 cm yr⁻¹, active fluid flow v_up>6 cm yr⁻¹). Under these conditions, methane efflux can be further enhanced by intermediate organic matter reactivity (RCM model parameter a=10–10² years) even though the control exerted by organic matter is only secondary with respect to the transport parameters. Finally intense local transport processes, such as bioirrigation (irrigation constant α₀>1 yr⁻¹), do also contribute to larger methane effluxes. Our results indicate therefore that present methane efflux from ESAS sediments can be supported by methane gas escape and non-turbulent CH₄ efflux from rapidly accumulating and/or active sediments (e.g. coastal settings, portions close to river mouths or submarine slumps). In particular, active sites sediments may release methane in response to the onset or increase in permafrost thawing or CH₄ gas hydrate destabilization.

High methane escape (up to 11–19 µmol CH₄ cm⁻² yr⁻¹ corresponding to 2.6–4.5 TgCH₄ yr⁻¹ if upscaled to the ESAS) can occur during a transient period following the onset of methane flux from the deep sediments. Under these conditions, substantial methane escape from sediments requires the presence of active fluid flow that supports a significant and rapid upward migration of the SMTZ in response to the onset of CH₄ flux from below. Such rapid and pronounced movements create a window of opportunity for non-turbulent methane escape by inhibiting the accumulation of AOM-performing biomass within the SMTZ – mainly through thermodynamic constraints – thereby perturbing the efficiency of the AOM biofilter. The magnitude of methane effluxes, as well as the duration of this window of opportunity, is largely controlled by the active flow velocity. In addition, results of transient scenario runs indicated that the characteristic response time of the AOM biofilter is of the order of few decades (20–30 years), thus exceeding seasonal–interannual variability. Consequently, seasonal variation of bottom methane and seawater sulfates exert a negligible effect on methane escape through the sediment–water interface.

AOM generally acts as an efficient biofilter for upward migrating CH₄ under environmental conditions that are representative for the present-day ESAS with potentially important yet unquantified implications for the Arctic ocean's alkalinity budget and, thus, CO₂ fluxes. Our results thus suggest that previously published fluxes estimated from ESAS waters to the atmosphere cannot be supported by non-turbulent methane efflux alone.

A regional upscaling of non-turbulent methane efflux for the Laptev Sea shelf using a model-derived transfer function that relates sedimentation rate and methane efflux merely sums up to ∼0.1 GgCH₄ yr⁻¹. Nevertheless, it also suggests that the evaluation of methane efflux from Siberian shelf sediments should pay particular attention to the dynamic and rapidly changing Arctic coastal areas close to big river mouths, as well as areas that may favour preferential methane gas release (e.g. rapidly eroding coastlines, fault lines or shallow sea floors, i.e. <30 m). In addition, our findings call for more data concerning sedimentation and active fluid flow rates, as well as the reactivity of depositing organic matter and bioirrigation rates in Arctic shelf sediments.

In conclusion, we argue that the evaluation of projected subsea permafrost thaw and/or hydrate destabilization impacts on the Arctic environment requires models that include an explicit description of (1) methane gas, (2) AOM biomass, as well as (3) the entire network of the most pertinent biogeochemical reactions. Such approaches, valid globally for all the shelves underlain by methane reservoirs (e.g. continental slopes), are even more recommended in order to enable a robust quantification of methane escape from the Arctic shelf to the Arctic ocean, settings even more sensible to the rapidly changing environmental conditions. Finally such refined modelling will also help evaluate the impact of subsea permafrost thaw and methane destabilization on Arctic alkalinity and biogeochemical cycling.

Appendix A: AOM efficiency η

If we identify the SMTZ region as the portion of the sediment column where the rate of AOM is 1 % of the maximum, we can define the efficiency of the AOM filter η as

$$\begin{array}{}\text{(A1)}& \mathit{\eta }\left(\mathit{\%}\right)=\left(\mathrm{1}-{\displaystyle \frac{J_{{\mathrm{CH}}_{\mathrm{4}}}^{+}}{J_{{\mathrm{CH}}_{\mathrm{4}}}^{-}}}\right)\cdot \mathrm{100},\end{array}$$

where $J_{{\mathrm{CH}}_{\mathrm{4}}}^{+}$ is the methane flux at the shallowest point where the AOM rate is 1 % of the maximum (upper dashed line in Fig. A1), and $J_{{\mathrm{CH}}_{\mathrm{4}}}^{-}$ is methane flux at the deepest point where the AOM rate is 1 % of the maximum (lower dashed line in Fig. A1).

Figure A1 Typical sediment profile of $\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]$, [CH₄] and AOM rate. Units are mM for concentration and mM yr⁻¹ for rate. The region between the two dashed lines represents the zone where AOM rate is larger than 1 % of it its maximum and defines the sulfate–methane transition zone (SMTZ). The fluxes $J_{{\mathrm{CH}}_{\mathrm{4}}}^{-}$ and $J_{{\mathrm{CH}}_{\mathrm{4}}}^{+}$ are the fluxes used in the definition of η of Eq. (A1).

Appendix B: Damköhler number

The Damköhler number Da is a dimensionless quantity which relates timescales typical of transport processes to timescales typical of chemical reactions. It compares the consumption/production rate with the advective transport and is defined as

$$\begin{array}{}\text{(B1)}& \mathit{\text{Da}}={\mathit{\tau }}_{\mathrm{T}}/{\mathit{\tau }}_{\mathrm{R}},\end{array}$$

where τ_T is the advective timescale and τ_R is the reaction timescale. τ_R is defined as 1∕K_R where K_R is the reaction rate of AOM or methanogenesis. If we call R the reaction rate then K_R reads as

$$\begin{array}{}\text{(B2)}& K_{\mathrm{R}}={\displaystyle \frac{\mathrm{1}}{\mathcal{L}}}\underset{\mathcal{L}}{\int }{\displaystyle \frac{R}{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}}\mathrm{d}z,\end{array}$$

where ℒ is the width where the reaction rate is larger than 1 % of the maximum rate. τ_T is instead defined as

$$\begin{array}{}\text{(B3)}& {\mathit{\tau }}_{\mathrm{T}}={\displaystyle \frac{\mathcal{L}}{|v_{\mathrm{up}}-\mathit{\omega }|}},\end{array}$$

where v_up−ω is the effective advective velocity. Da can be the expressed by

$$\begin{array}{}\text{(B4)}& \mathit{\text{Da}}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{T}}}{{\mathit{\tau }}_{\mathrm{R}}}}={\displaystyle \frac{\mathrm{1}}{|v_{\mathrm{up}}-\mathit{\omega }|}}\underset{\mathcal{L}}{\int }{\displaystyle \frac{R}{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}}\mathrm{d}z.\end{array}$$