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

. The East Siberian Arctic Shelf (ESAS) hosts large, yet poorly quantiﬁed 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 4 ) to the shallow coastal ocean and ultimately the atmosphere. However, a major unknown in assessing the contribution of this CH 4 ﬂux to the global CH 4 cycle and its climate feedbacks is the fate of CH 4 as it migrates towards the sediment-water interface. In marine sediments, (an)aerobic oxidation reactions generally act as a very efﬁcient 5 methane sink. Yet, a number of environmental conditions can reduce the efﬁciency of this bioﬁlter. Here, we used a reaction-transport model to assess the efﬁciency of the benthic methane ﬁlter 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 efﬁcient bioﬁlter. Yet, high CH 4 escape is simulated for rapidly accumulating and/or active sediments and can be further enhanced by the presence of organic matter with 10 intermediate reactivity and/or intense local transport processes, such as bioirrigation. In addition, in active settings, the sudden onset of CH 4 ﬂux triggered by, for instance, permafrost thaw or hydrate destabilization can also drives a high non-turbulent methane escape of up to 19 µ mol CH 4 cm − 2 yr − 1 during a transient, multi-decadal period. This "window of opportunity" arises due to delayed response of the resident microbial community to suddenly changing CH 4 ﬂuxes. A ﬁrst-order estimate of non-turbulent, benthic methane efﬂux from the Laptev Sea is derived as well. We ﬁnd that, under present day conditions, 15 non-turbulent methane efﬂux from Laptev Sea sediments does not exceed 1 Gg CH 4 yr − 1 . As a consequence, we conclude that previously published estimates of ocean-atmosphere CH 4 ﬂuxes from the ESAS cannot be supported by non-turbulent, benthic methane escape. To assess the inﬂuence of environmental conditions on the efﬁciency of the AOM bioﬁlter and its inﬂuence on non-turbulent methane emission from dissociating permafrost and/or disintegrating methane gas hydrates in ESAS sediments a set of ﬁve "one-at-time” parameter variation experiments is designed. It encompasses the most important controls on benthic methane cycling (Regnier Meister 2013; Egger and parameter variation experiments are performed for both the passive as well as active baseline scenario: 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 establishing a ﬁrst-order regional estimate of (not-turbulent) methane efﬂux and of potential methane consumption in Laptev Sea sediments. Model results reveal that AOM is an efﬁcient sink for upward migrating, dissolved methane in ESAS sediments. Simulated non-turbulent methane efﬂuxes 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 efﬁciency of the AOM bioﬁlter and, ultimately, on the methane efﬂux at the SWI. Both steady state and transient model results conﬁrm the key role of advective transport (mainly sedimentation and active ﬂuid ﬂow) in supporting methane escape from Arctic shelf sediments. Under steady state conditions, high methane efﬂuxes (up to 27.5 µ mol cm − 2 yr − 1 ) are generally found for sedi- 10 ments that are characterized by high sedimentation rates and/or active ﬂuid ﬂow (sedimentation rate ω > 0 . 7 cm yr − 1 , active ﬂuid ﬂow v up > 6 cm yr − 1 ). Under these conditions, methane efﬂux can be further enhanced by intermediate organic matter reactivity (RCM model parameter a = 10 − 10 2 yr) 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 α 0 > 1 yr − 1 ), do also contribute to larger methane efﬂuxes. Our results indicate therefore that present methane efﬂux from ESAS sediments can be supported by methane gas escape and non-turbulent CH 4 efﬂux 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 of permafrost thawing or CH 4 gas hydrate


Introduction
The Siberian Shelf represents the largest shelf on Earth (∼ 3 millions km 2 Wegner et al. (2015)) and spreads from the Kara 20 Sea to the Laptev, the East Siberian and the Chuckhi Sea. 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 Jr 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,5 2001; Romanovskii et al., 2005) by rising sea levels during the Holocene sea transgression (12 and 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 Jr and Koh, 2007).
Little is known about the total amount of carbon stored in subsea permafrost or its original extentPublished estimates of 10 carbon reservoir sizes diverge by orders of magnitude. For instance Shakhova et al. (2010a) estimate that 1175 PgC are locked in subsea permafrost on the ESAS alone, while McGuire et al. (2009) calculate that, across the entire Arctic shelf, 9.4 PgC reside in upper sediments and 1.5 − 49 PgC (2 − 65 PgCH 4 ) 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 15 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 20 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(Shakhova et al., , 2014. Shakhova et al. (2014) suggest that warming induced subsea permafrost thaw and hydrate destabilization may support methane emissions of up to 17 TgCH 4 yr −1 from the ESAS alone. Projected change in temperature (Shakhova et al., , 2019 due to climate 25 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., 30 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 35 (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 Arctic subsea permafrost and hydrates distribution, reservoir size and characteristics, but also a better quantitative understanding of Arctic methane sinks.
It has been shown that, in particular, high sedimentation rates (Egger et al., 2016), slow microbial growth (Dale et al., 2006(Dale et al., , 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 20 (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 longterm, changes in environmental conditions, namely SO 2− 4 concentration in sea water and availability of CH 4 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 25 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., 2014Shakhova et al., , 2015Shakhova et al., , 2017Leifer 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 30 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 1-D 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 wariming ESAS shelf.
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 5 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 en-10 countered 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 15 biomass dynamics and a bioenergetic rate law for AOM (Dale et al., 2006(Dale et al., , 2008c. 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 establishing 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 20 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) identifying 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 estimating the potential for non-turbulent CH 4 escape from ESAS sediments, thus providing first order constraints on the Arctic methane budget.

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 30 transformations of methane in a sediment column representative for ESAS conditions. For this purpose, we set-up a reaction network (Table S1, S2), 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 a set of coupled advectiondiffusion-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): C i is the concentration of the species i (mass per porewater volume for dissolved species or mass per solid matrix volume for 5 a solid species); ξ i.e. the porosity ξ = ϕ for dissolved species and ξ = ϕ s = 1−ϕ 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 = ω + v 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 10 coefficient (α i = 0 for solid species) and C i (0, t) is the concentration of the species i at the Sediment-Water Interface (SWI).
The reaction term S i is written as: where λ ij are the stoichiometric coefficients of all reaction rates R j that affect species i.

15
The effective diffusion coefficients D i are determined by correcting the diffusion coefficients in free solution D 0 i (Boudreau, 1997) for tortuosity θ and temperature. Tortuosity is calculated by means of porosity ϕ according to a modified Weissberg relation (Boudreau, 1997): θ = 1 − ln(ϕ 2 ). 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 CH4 (g). Bioturbation in the upper decimeters of the sediment is simulated using a diffusive 20 term (e.g., Boudreau (1986)), with a constant bioturbation coefficient, D 0 b . 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 (α 0 ) and the bioirrigation attenuation depth (z irr ) and is given by eq.

25
S9. Porosity is assumed to decrease with depth according to an exponential decay (Athy, 1930): with ϕ 0 : porosity at the Sediment-Water Interface (SWI) and c 0 : typical length scale for compaction.

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 Table S2 and Table 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 produc-5 tion/consumption of methane. A complete description can be found in the supplementary information.
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 10 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-R30). Mineral precipitation rates are simulated 15 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 4 ] * 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 20 and Wallmann, 2005;Tishchenko et al., 2005;Mogollón et al., 2009;Graves et al., 2017). [CH 4 ] * is calculated according to Dale et al. (2008a), derived from the formulation proposed by Duan et al. (1992) for which [CH 4 ] * 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 fall 25 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: It is the simplest and most commonly used formulation of the AOM rate in reaction-transport models (e.g. Regnier et al.
For transient model simulations, we apply a bioenergetic rate law in combination with an explicit description of the AOMperforming biomass (Dale et al., 2006(Dale et al., , 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: where µ g is the growth rate and µ d is the decay rate. F K is the kinetic constraint given by: with K where R is the gas constant, T is the absolute temperature, χ is the average number of electrons transferred per reaction per 10 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 − ) −1 is the minimum energy needed to support synthesis of ∼ 1 3 − 1 4 mol ATP (Dale et al., 2008c). In order to be thermodynamically favorable 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 15 with ∆G 0 r : standard free energy of the reaction, the second term: deviations from standard conditions (temperature and reaction quotient) on Gibbs free energy and γ: a parameter representing departure from ideal beahviour.

Boundary conditions 25
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 (∂C/∂z = 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.

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 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  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 Sections 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 15 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 (NH + 4 ), phosphates (PO 3− 4 ) and dissolved inorganic carbon (DIC) depth profiles. While the former two species are of main interest for evaluating 20 simulated AOM dynamics, the remaining three serve as indicators for OM degradation dynamics since they are metabolic byproducts of degradation (see Table S2). Moreover NH + 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 [NH + 4 ] profile, but does not affect the overall shape of the NH + 4 depth profile. Similarly, PO 3− 4 is only consumed by fluorapatite precipitation (R22) and adsorption processes (R29 and R31). 25 Fluorapatite precipitation controls maximum dissolved PO 3− 4 concentrations, while the mineral adsorption process (R29) exerts a homogeneous influence and the interaction with Fe(OH) 3 is expected to be minor and mainly affects PO 3− 4 within the iron reduction zone. To evaluate the main physical and biogeochemical controls on the efficiency of the AOM biofilter and its impact on nonturbulent 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 5 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 −1 ; 2. an active case, i.e. with v up = 1 cm yr −1 , a value which falls within the range of fluid flow velocities v up = 0.005 − 30 cm yr −1 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  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 4 ] concentrations at the lower model boundary, mimicking different methane fluxes from thawing subsea permafrost and/or disintegrating methane gas hydrates at depth. The  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. 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 5 both the passive as well as active baseline scenario: 1. Sedimentation rate ω. The sedimentation rate is varied over two orders of magnitude (0.03 − 0.123 − 0.17 − 1.5 cm yr −1 ). 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 −1 , inferred from Overduin et al. (2015)), 10 while minimum values are representative of sedimentation rates found in the East Siberian Arctic Sea (Stein et al. (2001) in Levitan and Lavrushin (2009)  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 15 by subsea permafrost or gas hydrates (Judd and Hovland, 2009;Semenov et al., 2019). Therefore, v up is varied from k AOM (eq. 5) is thus varied over the range k AOM = 5 · 10 2 − 5 · 10 3 − 5 · 10 4 − 5 · 10 5 − 5 · 10 6 − 5 · 10 7 M −1 yr −1 .
4. Organic matter reactivity (i.e. RCM parameter). 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 SMTZ depth (Regnier et al., 2011;Meister et al., 2013). Thus, ν was kept 25 constant, while a was varied over the entire range of previously published values a = 0.1 − 1 − 10 − 100 − 500 − 1000 yr (Arndt et al., 2013). Studies about ESAS organic matter degradation shows a reactivity of deposited organic matter which is compatible with the RCM parameter we explored. For instance, Bröder et al. (2016) found an half life for the organic carbon in the East Siberian Arctic Shelf of 19 − 27 yr, which would correspond to an a = 3.4 − 4.8 yr, with ν = 0.125. 5 5. Bioirrigation coefficient α 0 . 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 meters  are detrimental to the local fauna, thus suggesting a low bioirrigation intensity. Yet, 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, 10 ice scouring might also enhance non-local transport seasonally. We therefore, varied α 0 over the entire range of plausible Thullner et al., 2009). 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, 15 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  2. Seasonal CH 4 + SO 2− 4 : 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 sulphate concentrations characterise the melt season. The bottom boundary condition for methane [CH 4 ] − follows an opposite trend: it is set to zero during the winter months and increases in Arctic summer.

Transient Sensitivity Study
3. Linear CH 4 : slow increase in methane supply from permafrost thaw and/or hydrate destabilization. A linear increase of the bottom boundary methane concentration [CH 4 ] − (from 0 up to the peak) over 200 years is applied. 4. Sudden CH 4 : abrupt increase of methane supply from permafrost thaw and/or hydrate destabilization. An instantaneous change of bottom boundary methane concentration -from 0 to one of the peak value [CH 4 ] − -is applied.

Analyzed 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; 5 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 SO 2− 4 and CH 4 at the SMTZ, the maximum and integrated AOM rate and the Damköhler number (D a ) for AOM and methanogenesis are also calculated. Damköhler number is defined as eq. B4 (see Appendix B) and sets the ratio between the typical transport time-scale and the typical reaction time-scale. If D a < 1, the reaction time-scale is longer than 10 transport time-scale (i.e. the reaction is slower) and the process is reaction-limited. If D a > 1 the process is transport-limited.
3 Results and discussion 3.1 Case study: sediment core on the Laptev Sea shelf  Table S4, the measured organic carbon content in Table S5. For O2 no measured profile is available.
dynamics. Such heterogenity is not incorporated in the model and accounting for such a heterogeneity would require additional information.
3.2 Main physical and biogeochemical controls on potential non-turbulent methane flux from ESAS sediments

General patterns of methane and sulfate cycling on the ESAS
The comprehensive ensemble of all sensitivity experiments allows exploring the general patterns of methane and sulfate cycling 5 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 4 supply from below. For most of the investigated environmental conditions (95% of the runs), 95-99.9% of the upward diffusing CH 4 is consumed within the SMTZ, resulting in very small or negligible methane effluxes (≤ 10 −2 µmolCH 4 cm −2 yr −1 ) from the sediment. If upscaled to the total area of the ESAS (∼ 1.485 · 10 6 km 2 , Wegener et al. (2015)), for which methane outgassing estimates have been published, the smallest simulated non-turbulent methane flux (i.e. 1.4 · 10 −13 µmol cm −2 yr −1 , Fig. 2.b) would sum up to a total flux of 2.1 mmolCH 4 yr −1 , resulting in a negligible role of non-turbulent, benthic methane fluxes to the Arctic methane budget.
Yet, 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 4 escape from ESAS sediments can reach values of up to 27 µmolCH 4 cm −2 yr −1 . Simulation analysis shows that these high effluxes and, thus, low AOM biofilter efficiencies are generally obtained for 5 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. 2.a). 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 10 methane flux at the SWI and the SMTZ depth ( Fig.2.b). Such a relationship reveals the pivotal connections between these two quantities and mirrors the empirically found linear log-log relationship between measured CH 4 fluxes at the SMTZ and the SMTZ depths (  (2011)). Upscaling the highest simulated non-turbulent flux (27.48 µmolCH 4 cm −2 yr −1 ) to the ESAS results in a total efflux of 0.408 TmolCH 4 yr −1 = 6.52 TgCH 4 yr −1 . 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 20 slope sediments on a centennial timescale (4.73 TgCH 4 yr −1 , Kretschmer et al. (2015)).
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 (D a M G ) and AOM (D a AOM ) are reported in Fig. S2. D a M G (purple circles) are < 1 , span a range of ∼ 0.0021 − 0.43 and are thus comparable to previously reported D a M G of 0.22 for methane gas hydrate bearing sites, such as Hydrate Ridge and Kithley Canyon (Chatterjee et al., 2011). They reveal that 25 methanogenesis is always slower than methane transport and that CH 4 dynamics driven by methanogenesis are thus reactionlimited. 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 D a AOM values (D a AOM =32-2.78 · 10 5 - 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 30 SMTZ and subsequently the SWI.

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;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 (JCH4) at the SWI versus SMTZ depth. Meister et al., 2013;Winkel et al., 2018). Yet, they also challenge intuitive views on the factors that favour high CH 4 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. They also shed light on the mechanisms behind 5 non-turbulent methane escape from ESAS sediments.  fig. 3.c) from 5.5 · 10 −15 µmolCH 4 cm −2 yr −1 for low sedimentation rates (ω = 0.03 cm yr −1 ) to values as high 10 as 27.5 µmolCH 4 cm −2 yr −1 for high sedimentation rates (ω = 1.5 cm yr −1 ). 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 4 fluxes and low AOM efficiencies in these rapidly accumulating sediments, are enhanced methanogenesis rates. High sedimentation rates facilitate not only the supply of organic matter to the methanogenic zone of the sediment, but also reduce residence times in the 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). upper sediment layer, resulting in a lower OM age (see eq. S14, S16)/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.

Role of advective transport
In addition, the presence of active fluid flow further enhances methane efflux. The CH 4 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 5 settings on the shelf. Table S8 shows that for low to intermediate sedimentation rates, D a AOM values significantly decrease with v up , indicating that less and less methane consumption occurs within the typical transport time scale τ 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. Yet, flux differences quickly decrease with increasing sedimentation rates. Rapidly accumulating sediments show almost no difference 10 in efflux between active and passive sites (Fig. 3.a). 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 4 , shifting the SMTZ upwards and, thus, increasing CH 4 concentrations at shallow sediment depths (see Fig. 3.d). The apparent paradox of the CH 4 efflux insensitive to fluid flow in fast accumulating sediments can be resolved by examining the dissolved CH 4 depth profiles (Fig. 4). Simulated depth profiles are nearly identical and reveal CH 4 concentrations at or near the 15 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 exceed non-turbulent CH 4 fluxes up to 6 times. In rapidly accumulating, active and passive sediments, non-turbulent CH 4 fluxes are thus essentially identical. However, active settings will be characterised by the additional build-up of gaseous CH 4 and its potential escape through the sediment-water 20 interface-a process not simulated in the present study.
Model results thus show that the dominant mechanism behind the observed transport-control on non-turbulent CH 4 efflux is an overall increase in CH 4 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 section 25 S1) 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 section 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 4 fluxes, as well as lowering the efficiency of the AOM barrier 30 from η ∼ 100% to η ∼ 78%. The important role of the SMTZ location as a key control on CH 4 efflux is further confirmed by the observed exponential relationship between the location of the SMTZ and ω (Fig. 3.d). 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 −1 give rise to high non-turbulent CH 4 fluxes (20 − 80 µmolCH 4 cm −2 yr −1 ) and a high OM burial efficiency (∼ 78%, Egger et al. (2016)). Egger and co-workers explained these findings by the slow growth of AOM microorganisms and the resulting inability of the microbial community to consume all of the CH 4 produced. Yet, 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 4 effluxes.

Role of organic matter quality
The quality of organic matter deposited onto the sediment exerts an additional control on CH 4 efflux. Fig. 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 10 control. This also means that, in order to assess the main features of possible CH 4 efflux in terms of modeling, capturing the details of organic matter quality is not fundamental. Maximum fluxes are generally simulated for rapidly accumulating sediments ω > 0.5 cm yr −1 that receive organic matter of intermediate quality (a = 10 − 100 yr). 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 yr), 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 yr), 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 4 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 5 supply from below (i.e. higher CH 4 concentration at the lower boundary) increase the CH 4 efflux through the SWI without altering the overall patterns (see Fig. 5.a-b vs 5.c-d).

Role of non-local transport
The analysis of the influence of bioirrigation on non-turbulent CH 4 efflux from the ESAS (Fig. S12) shows that such a process enhances methane efflux in sediments that are characterized by a shallow SMTZ. Yet, bioirrigation exerts a limited effect under 10 a range of environmental conditions that favour a deep or shallow SMTZ location respectively.
In passive settings, changes in bioirrigation coefficient (α 0 ) exert a limited control on CH 4 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 α 0 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. 15 Methane effluxes first increase by up to 5 orders of magnitude from α 0 =0 yr −1 to α 0 =5 yr −1 , reaching maximum effluxes of ∼ 0.02 µmolCH 4 cm −2 yr −1 , before remaining almost constant with a further increase in bioirrigation coefficients (up to 240 yr −1 ). 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 4 dynamics over the 3-dimensional transport coefficient ω, 20 v up and α 0 space shown in Fig.6.
A comparison between simulations with α 0 = 0 yr −1 and α 0 = 0 yr −1 (α 0 = 5 yr −1 , α 0 = 10 yr −1 and α 0 = 33 yr −1 ) shows that irrigation increases the CH 4 efflux at low to intermediate sedimentation rates and/or high v up (lower-left corner of the phase space in both plots). Yet, 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 25 α 0 = 0 yr −1 ; irrigation when α 0 = 0 yr −1 ). 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 section 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 4 efflux at the SWI are broadly similar for all of the investigated α 0 = 0 yr −1 (Fig. 6.b,c,d). It is worth noticing that, independently on the α 0 , CH 4 efflux for ω = 0.03 cm

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, modeling reveals that k AOM only plays a minor role for non-turbulent methane fluxes across the SWI (see Fig. 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 5 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 4 ] and [SO 2− 4 ] 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 analyzed. This is in disagreement with results by Dale et al. (2008c), which show that, in dynamic settings 10 subject to large methane fluxes, an increase of 3 orders of magnitude in k AOM (from 10 2 M −1 yr −1 to 10 5 M −1 yr −1 ) leads to a reduction in steady state methane fluxes below 10 −2 µmolCH 4 cm −2 yr −1 . However this discrepancy might be ascribable to the high water flow velocity employed in their simulation (v up = 10 cm yr −1 ), ten 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 15 CH 4 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 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.

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 4 effluxes from thawing subsea permafrost and/or dissociating methane gas hydrates (larger than 4 µmolCH 4 cm −2 yr −1 ) are promoted by intense advective transport (sedimentation rate ω > 1 cm yr −1 , active fluid flow v up > 7 cm yr −1 ). Under these conditions, CH 4 efflux can be further enhanced by moderate 25 OM reactivity (a = 10 − 10 2 yr) and intense non-local transport processes, such as bioirrigation (irrigation constant α 0 > 0 yr −1 ). Overall, non-turbulent benthic escape of CH 4 from deep sources appears to be mainly controlled by the concurrent effects of ω, v up and α 0 . In contrast, maximum AOM rates, k AOM , exert no influence on the AOM filter efficiency.

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;Capet et al., 2016;Dale et al., 2017;Bowles et al., 2014). Such a quantitative framework in which first-order estimates of potential non-turbulent methane escape from ESAS sediments can also be derived from the results of the model sensitivity study.
Model results indicate that sedimentation rate exerts the dominant control on benthic escape of methane from thawing subsea 5 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. 3.c thus allows estimating 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 10 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 3D kriging method (see Fig. 7.a), using the International Bathymetric Chart of Arctic Ocean (IBCAO) (Jakobsson et al., 2012) and employing longitude, latitude and water depth as predictors for ω.  Table 2. Estimated flux of CH4 at SWI in mol yr −1 for different depth regions of Laptev Sea in a passive (vup = 0 cm yr −1 ) and active (vup = 1 cm yr −1 ) case.
vup Region (water depth, area) 0 1 0 − 10 m, 7.7 · 10 4 km 2 6.5 8.9 · 10 5 10 − 80 m, 4.5 · 10 5 km 2 296.2 8.5 · 10 6 Estimated non-turbulent methane effluxes corresponding to the highest measured sedimentation rates close to the Lena mouth do not exceed 1.57 · 10 −1 µmolCH 4 cm −2 yr −1 assuming the presence of active fluid flow and 2.25 · 10 −5 µmolCH 4 cm −2 yr −1 for passive settings. These findings are not surprising as steady state sensitivity results indicate that high CH 4 efflux 25 requires sedimentation rates of ω > 1 cm yr −1 . The regional non-turbulent CH 4 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 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. Bottom (Log) Values of the potential non-turbulent dissolved methane emissions at the SWI considering the relationship presented in Fig. 3.c for passive (b) and active (c) cases.
non-turbulent CH 4 efflux is negligible. Even if we assume the omnipresence of an active fluid flow of v up = 1 cm yr −1 , the estimated non-turbulent methane efflux merely sums up to 9.39 · 10 6 molCH 4 /yr (∼ 0.1 GgCH 4 /yr) over the entire Laptev Sea area of 527.4·10 3 km 2 . 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 −1 ), possibly in connection with active seepages, groundwater discharges and fault lines (the latter 5 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 island), could result in methane effluxes of up to 10 − 10 1.3 µmolCH 4 cm −2 yr −1 (see Fig. 5 and Fig. 6). However, such high fluid flow velocities would be only found locally and would thus merely give rise to a number of methane emission hot spots 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 10 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 yr or a = 1 yr. 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 15 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 20 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 4 yr −1 (Berchet et al., 2016;Thornton et al., 2020) or even tens of teragrams of CH 4 yr −1 (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 4 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 25 CH 4 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 30 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.

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. Yet, 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 5 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 Section 2.3.2).
10 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 for all environmental change scenarios barely reaches 3-4 µmolCH 4 cm −2 .
In contrast, active settings (i.e. v up = 1 cm yr −1 ) exhibit an initial increase in CH 4 fluxes to maxima of 0.55-0.83 µmolCH 4 cm −2 yr −1 over the first 50 years. This growth coincides with a rapid upward shift of the SMTZ by 100 cm. Methane escape   (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. Fig. 10 illustrates the 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). During the first 23 years, AOM biomass thus remains largely constant ( Fig. S21.a) 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 5 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 4 efflux increases.
The initial methane efflux is largely supported by in situ methanogenesis since the advective transport of methane (occurring at v up − ω = 0.877 cm yr −1 , corresponding to 20.17 cm in 23 years) is too slow to allow methane from below 3 m to reach the sediment-water interface. 10 After the first 23 years, thermodynamic constraints ease and AOM begins to efficiently consume upward migrating methane at the SMTZ by 40% (Fig. 8.a). 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 velocityv = v SM T Z +v up −ω = 3.337 cm yr −1 is therefore 100.4 cm 3.337 cm yr −1 = 30.1 yr. Consequently, methane efflux further increases. This methane efflux is now also supported from deep sources such as thawing permafrost and/or dissociating 15 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 −1 ). 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 µmolCH 4 cm −2 yr −1 for v up = 1 cm yr −1 to 11 − 19 µmolCH 4 cm −2 yr −1 for v up = 5 cm yr −1 . Yet, 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 20 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. 23 + 30.1 = 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. 8.b). 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. 8.a). 25 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.

Final new steady state
The final new steady state value of methane efflux ( Fig. 10.b, S17 and S19) is generally in good agreement with Dale et al.  two orders of magnitude larger than the efflux of ∼ 0.1 µmolCH 4 cm −2 yr −1 simulated in the steady state simulations, with bimolecular rate law, under identical environmental conditions (inferred from Fig. 6). AOM rate according to the bioenergetic formulation (blue) and, for comparison, according to bimolecular formulation used for the steady-state simulations (red). c. Apparent kAOM , estimated from eq. 5.
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. Panel 9.a illustrates that the concentration product [CH 4 ] · [SO 2− 4 ] is wider than the AOM rate profile (panel 9.b, blue curve). Fig. 9.c also shows that the apparent k AOM is not uniform: it reaches a maximum value of 109 M −1 yr −1 , but remains well below 100 M −1 yr −1 at most depths. Compared to the values typically 5 applied for bimolecular rate laws (i.e. k AOM = 10 2 − 10 7 M −1 yr −1 ), 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. 10.c). 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). 10 Integrated biomass ΣB ranges from ∼ 1.2 · 10 10 to 3.5 · 10 11 cells cm −2 (except for simulation with v up = 5 cm yr −1 and [CH 4 ] − = 5.455 mM, whose ΣB= 1.2·10 12 ). These values are comparable with AOM biomass reported in Treude et al. (2003) (1.5 − 1.8 · 10 10 cells cm −2 ) or with values simulated in Dale et al. (2008c) (3.7 · 10 11 cells cm −2 for v up = 5 cm yr −1 ). In addition, the maximum simulated biomass for active settings (0.5 − 2.5 · 10 10 cells cm −3 ) agrees well with previously reported values, ranging from 0.27 to 7.4 · 10 10 cells cm −3 (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 5 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 −1 and [CH 4 ] − = 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 10 cells cm −3 (Lösekann et al., 2007;Knittel and Boetius, 2009).

10
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 establishing 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 5 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 −2 yr −1 ) are generally found for sedi-10 ments that are characterized by high sedimentation rates and/or active fluid flow (sedimentation rate ω > 0.7 cm yr −1 , active fluid flow v up > 6 cm yr −1 ). Under these conditions, methane efflux can be further enhanced by intermediate organic matter reactivity (RCM model parameter a = 10 − 10 2 yr) 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 α 0 > 1 yr −1 ), do also contribute to larger methane effluxes. Our results indicate therefore that present methane efflux from ESAS 15 sediments can be supported by methane gas escape and non-turbulent CH 4 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 of permafrost thawing or CH 4 gas hydrate destabilization.
High methane escape (up to 11-19 µmolCH 4 cm −2 yr −1 corresponding to 2.6-4.5 TgCH 4 yr −1 if upscaled to the ESAS) can occur during a transient period following the onset of methane flux from the deep sediments. Under these conditions, 20 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 4 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 25 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 sea water sulfates exert a negligible effect on methane escape through the sediment-water interface.
AOM generally acts as an efficient biofilter for upward migrating CH 4 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, 30 thus, CO 2 fluxes. Our results thus suggest that previously published fluxes estimated from ESAS waters to 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 4 yr −1 . 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 favor 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 5 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 modeling will also help evaluate the impact of subsea permafrost thaw and methane destabilization on 10 Arctic alkalinity and biogeochemical cycling. The Damköhler number D a is a dimensionless quantity which relates time scales typical of transport processes to time scales typical of chemical reactions. It compares the consumption/production rate with the advective transport and is defined as 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: where L is the width where the reaction rate is larger than 1% of the maximum rate. τ T is instead defined as where v up − ω is the effective advective velocity. D a can be the expressed by: Dale, A. W., Aguilera, D., Regnier, P., Fossing, H., Knab, N., and Jørgensen, B. B