Gene-Based Modeling of Methane Oxidation in Coastal Sediments: Constraints on the Efficiency of the Microbial Methane Filter

Methane is a powerful greenhouse gas that is produced in large quantities in marine sediments. Microbially mediated oxidation of methane in sediments, when in balance with methane production, prevents the release of methane to the overlying water. Here, we present a gene-based reactive transport model that includes both microbial and geochemical dynamics and use it to investigate whether the rate of growth of methane oxidizers in sediments impacts the efficiency of the microbial methane filter. We focus on iron- and methane-rich coastal sediments and, with the model, show that at our site, up to 10% of all methane removed is oxidized by iron and manganese oxides, with the remainder accounted for by oxygen and sulfate. We demonstrate that the slow growth rate of anaerobic methane-oxidizing microbes limits their ability to respond to transient perturbations, resulting in periodic benthic release of methane. Eutrophication and deoxygenation decrease the efficiency of the microbial methane filter further, thereby enhancing the role of coastal environments as a source of methane to the atmosphere.

(91 kBq) was injected in the syringes.The sediment was incubated for 24 h in the dark under a nitrogen atmosphere, after which it was transferred to 50 mL centrifuge tubes containing 20 mL of deoxygenated 20% zinc acetate to precipitate dissolved H 2 S and inhibit biological activity [Fossing and Jørgensen, 1989;Kallmeyer et al., 2004].Upon analysis, samples were rinsed twice using deoxygenated bottom water (10 mL) and centrifuged for removal of unreacted 35 SO 2 - 4 [Egger et al., 2016b].The reduced S was extracted with an acidic chrome chloride solution (48 h) via the passive diffusion method [Burton et al., 2008].Sulfate S1 reduction rates were calculated by comparing the activity (decays per minute) of the radiolabeled total reduced inorganic sulfur (aTRIS) to the total SO 2 -  4   (aTOT) radiotracer [Kallmeyer et al., 2004].
A.1.2Fe and Mn reduction and NH + 4 production rates Two cores were sliced under a nitrogen atmosphere for 6 sediment depth intervals of 3 cm (0.5-3.5, 3.5-6.5, 10-13, 20-23, 30-33, 40-43 cm) into beakers in which the sediment was subsequently homogenized.Additionally, the surface layer 0-0.5 cm was sampled in a 50 ml greiner.Samples from the beakers were taken at ca. t = 0, 8, 20, 30 and 48 hours under a nitrogen atmosphere.For the top sediment, samples were only taken at t=0 and t=48 because of the limited amount of material.
During sampling, approximately 10 ml of sediment from each beaker was transferred into a 50 ml greiner tube, which was centrifuged at 4500 rpm for 10 minutes to extract porewater.Subsequently the porewater was filtered over 0.45 µm and subsampled for ICP-OES (2 mL, acidified with 20 µL 35% suprapur HCl) and NH 4 (remainder of sample).ICP-OES samples were stored at 4 • C, NH + 4 samples were stored at -20 • C. All beakers and the greiner were stored under a nitrogen atmosphere at in-situ temperature during incubation.Dissolved Fe, Mn and NH + 4 concentrations were determined in the same way as porewater samples.Fe and Mn reduction and NH + 4 production rates were determined by fitting a linear regression to the concentration time series for each beaker/greiner.Only regressions with R 2 >0.5 were considered.

A.1.3 CH 4 production rates
Sediments for incubations to measure potential CH 4 production rates were sliced into intervals of 3 to 4 cm (0-4, 9-12, 21-24, 33-36, 49-52 and 69-72 cm) under a nitrogen atmosphere.The slices were stored anoxically in the dark at 4 • C for one month from sample collection until bottles were assembled.
For this, 5 g of wet sediments was placed in 60 ml-serum bottles, and sulfate-free artificial seawater (ASW) medium pH 7.5 was added to create a 1:1 diluted slurry.The ASW medium was adapted from that of [Kester et al., 1967] to achieve a salinity of 5.3.No trace elements or vitamins were added.
Methane production was monitored via injection of 100 µl-headspace samples into an HP 5890 gas chromatograph with a detection limit of <1 ppm.Each gas sample was measured in triplicate and S2 results were averaged.

A.1.4 Bromide incubations
A sediment core was incubated with the inert tracer bromide to determine bioirrigation rates [Martin and Banta, 1992].Directly after core retrieval, the volume of the overlying water was adjusted to ca.
500 mL and a concentrated bromide solution was added to the overlying water to achieve a concentration of ca. 3 mmol L −1 following [Lenstra et al., 2019].During a two-day incubation at in-situ temperature, the overlying water was kept saturated with O 2 and well-mixed by bubbling with air.After incubation, the top 20 cm of sediment was sliced into intervals of 0.5 to 2 cm and centrifuged for 20 min at 4500 rpm.After centrifugation, the supernatant was filtered (0.45 µm) and stored at 4 • C.
Bromide concentrations were determined with ion chromatography.Subsequently porewater bromide was modeled using a 1-D, nonlocal exchange function [Emerson et al., 1984;Boudreau, 1984].Here, dissolved bromide is described as (Eq.A.1), where ϕ is the depth dependent sediment porosity, z is depth in the sediment (cm) and D s is the sediment diffusion coefficient for bromide.The molecular diffusion coefficient D s at site NB8 was corrected for tortuosity in the porous medium [Boudreau, 1996a] and for the ambient salinity S, temperature T (Table A.4) and pressure using the R package CRAN: marelac [Soetaert et al., 2010], which implements the constitutive relations listed in [Boudreau, 1997].α is the nonlocal bioirrigation function for bromide, which is assumed to be time-invariant.c is the concentration of bromide in the porewater and c 0 is the bromide concentration of the overlying water.In the model, porosity is interpolated linearly for every depth layer in the model (1 layer every mm) between measured points.The bioirrigation function (α) was determined by fitting porewater bromide depth profiles after incubation to equation A.1.

S3
A.2 Reactive transport model

General model description
The reactions in the reactive transport model describe organic matter degradation coupled to different electron acceptors.Reactions are divided in primary redox reactions and other biogeochemical reactions (Table A .6).The succession of oxidants during organic matter degradation [Froelich et al., 1979] is described by means of Monod kinetics (Table A .3;Boudreau [1997]).This means that the oxidants with the highest metabolic free energy yield are preferentially used until they become limiting and the oxidant with the next highest energy yield is used [Berg et al., 2003;Boudreau, 1996b;Wang and Van Cappellen, 1996].Respiratory reactions occur where O 2 , NO - 3 , MnO 2 , Fe(OH) 3 and SO 4 2− serve as electron acceptors, and finally, organic matter is subject to methanogenesis (Lenstra et al. [2018]; Table A.3). Organic matter includes carbon (C) and nitrogen (N) in a C:N ratio of 7.1:1 (Table A .4).
In the model, oxidation of CH 4 is possible with O 2 , NO - 3 , MnO 2 , Fe(OH) 3 and SO 4 2− .Due to strong variations in the incorporation of Mn and other cations in the structure of vivianite [Rothe et al., 2016;Kubeneck et al., 2021], this mineral is not included in the model.Dissolved inorganic carbon in the model is calculated as the sum of the carbon in CO 2 and HCO 2 - 3 , which is produced or consumed by modeled reactions (Table A .6).
The generic mass conservation equations for solids and solutes are described by Eq.A.2 and A.3; where C s is the concentration of solid species (mol L −1 ), C aq is the concentration of dissolved species (mol L −1 ), t is time (yr), φ is the sediment porosity, v and u are the advective velocities of solid and dissolved species (cm yr −1 ), respectively.Variables v and u were described by a depth-dependent S4 function to account for changes in porosity [Meysman et al., 2005].Distance from the sedimentwater interface is z (cm), D ′ is the diffusion coefficient of dissolved species (cm 2 yr −1 ), corrected for tortuosity in the porous medium [Eq. A. 4 Boudreau, 1996c].∑ R s and ∑ R aq are the net reaction rates from the chemical reaction (Table A.3) for solid and dissolved species.
Porosity (φ) is described by Eq.A.5 to account for sediment compaction [Meysman et al., 2005;Reed et al., 2011b], where φ 0 is the porosity at the sediment-water interface, φ ∞ is the porosity at depth and y is the porosity attenuation factor/e-folding distance (Table A.4).The model code was written in R with the use of the marelac geochemical dataset package [Soetaert et al., 2010].To calculate the transport in porous media, the R package Reactran was used [Soetaert and Meysman, 2012].The set of ordinary differential equations was solved numerically with the Lsode integrator algorithm [Petzold, 1983].
Zero gradient boundary conditions were applied to the base of the model domain for all chemical species.The total depth of the model was set to 80 cm (divided into 800 grid cells of 0.1 cm).Reaction parameters were mostly taken from literature or obtained within existing parameter ranges (Table A.7).
If these were not available, or no fit to the data could be obtained with existing ranges, parameters were constrained by fitting the model to the measured data.

A.2.1 Microbial modeling and equations
The model considers 4 different microbial groups that couple the oxidation of CH 4 to O 2 , SO 2 - 4 , Fe oxides and Mn oxides, respectively (Table A.6: R27-R32).Several archaea have been shown to mediate the anaerobic oxidation of CH 4 (AOM), including the ANME-1, ANME-2a-c, Methanoperedenaceae, and the ANME-3 [Leu et al., 2020].It has been suggested that archaea can switch oxidation pathways from SO 2 - 4 to for example Fe and Mn oxides, but this does not hold for all groups.For example ANME-1 can perform S-AOM [Boetius et al., 2000] but is likely not involved in Fe-AOM [Aromokeye et al., 2020].In the model, we therefore decided to include 4 different groups of mi-

S5
crobes that correspond to a particular metabolisms without accounting for switches between different metabolisms.We also did not include cells that facilitate the oxidation of CH 4 coupled to NO x because this is likely a negligible pathway in marine sediments [Jørgensen, 2021].The production of cells at any depth is driven by the release of energy from their reactions, and is proportional to the Gibbs free energy multiplied by the reaction rate.The cell specific reaction rate (mol yr −1 cell −1 ) in given as follows, where, V r is the maximum gene-specific rate (mol cell −1 yr −1 ; Table A where, γ e r is the negative stoichiometric coefficient and ∆G r is the Gibbs free energy of the reaction (kJ mol −1 ).The Gibbs free energy of the reaction is calculated using, Here, ∆G 0 r is the standard Gibbs free energy of the reaction (kJ mol −1 ) that depends on the local temperature and pressure as calculated using the CHNOSZ R package [Dick, 2008].R g is the gas constant (8.3145J K −1 mol −1 ), T is bottom water Temperature (K) and Q r is the reaction quotient (unitless).
Ft is the dimensionless function that varies between 0 (complete kinetic/thermodynamic limitation) and 1 (no kinetic/thermodynamic limitation), where, F is the Faraday constant (96.485 kJ V −1 ), ∆Ψ is the electric potential across the membrane (0.12 mV).Both constants are taken from Reed et al. [2014].

A.2.2 Model parameterization
Because of a lack of reliable information, maximum rate constants of R27-R32 (Table A.6) were calibrated to geochemical depth profiles and measured rate depth profiles.The cell mass used in equation 1 (main paper) is 5 * 10 −13 g cell −1 , similar as in Louca et al. [2016].
Cell death rates (q r ) are constant and are related to their maximum growth rate (Table A

A.2.3 Transient modeling scenario
The model was run to steady state in 200 years.Subsequently, temporal changes were implemented that follow the periodic pulses of organic matter, Fe oxides and Mn oxides that occur ca.every 20 years The transient scenario as described above (i.e.our baseline scenario) is subsequently used in a sensitivity analysis.In this analysis, the effects of variations in single parameters on the system were evaluated by running the scenario again but in this case changing (1) the bottom water salinity (0-25); (2) the bottom water O 2 concentration (0-275 µmol L −1 ); (3) the organic matter deposition (factor S7 0.01-2) and ( 4) the Fe and Mn oxide deposition (factor 0.5-2) from the original values.We note that organic matter in our model consists of both carbon and nitrogen (Table A.3).
A.3 Discussion of CH 4 production and SRR CH 4 production rates Methanogenesis rates below the SMTZ determined in the laboratory were much higher than modeled rates in our RTM (Table A.8).There are three likely explanations for this difference.Firstly, mixing of sediments in slurry incubations can lead to higher reaction rates because of easier access of microbes to substrates than is the case in-situ.Secondly, the rate measurements were not performed at insitu temperature.While the bottom water temperature was ca.2.8 • C, the methanogenesis rates were determined at room temperature, i.e. at ca. 21 • C.This likely led to enhanced rates of microbial activity, including CH 4 production.Thirdly, methanogenesis rates were determined ca. 1 month after sampling.During sample storage, other electron acceptors used in organic matter degradation, such as Fe and Mn oxides, which are known to be present at depth in these sediments, may have been lost, allowing methanogenesis to take over.For a detailed discussion of the problems with determining methanogenesis rates we refer to section 7 in Reeburgh [2007].

A.4 The sensitivity of changes in maximum growth rate and death rate
The model is very sensitive to changes in the maximum growth rate of the microorganisms (varied between 0.5 and 1.5 times the values in the baseline scenario).Porewater depth profiles of SO .When the maximum growth rate increases, the SMTZ is located relatively close to the sediment water interface.We find that for FeOx-and MnOx-ANMEs there is an optimum maximum growth rate, which is similar to the baseline scenario, where their abundances are high.
We found that the our model is not very sensitive to changes in the death rate of the microorgan- Before inc.After inc.

Cell abundance cells cm −3
Cell abundance cells cm −3

Cell abundance cells cm −3
Cell abundance cells cm −3

Cell abundance cells cm −3
Cell abundance cells cm −3

Cell abundance cells cm −3
Cell abundance cells cm −3
are however several reasons why the determined SRR, especially below the SMTZ, might not reflect actual in-situ SRR.Firstly, cryptic sulfur cycling, where sulfide is oxidized by Fe oxides, forms SO 2 - 4 and is subsequently reduced[Holmkvist et al., 2011], is not included in the model.In our model, SO 2 SMTZ.If the measurements would indeed reflect the actual in-situ SRR, the difference between our model result and the measurements could be used to quantify the rate of SO 2 - during coring.Previous work has shown that this can lead to a slight overestimation of the SO 2 - 4 concentration at depth[Pellerin et al., 2018].The low SO2 -  4     concentrations can lead to an overestimation of the SRR.For a detailed discussion of the problems with SRR measurements and their comparison with modeled SRR we refer to section 3 in Jørgensen [2021].

Figure A. 1 :Figure
Figure A.1: Results of bromide tracer incubations at station NB8.Black and red diamonds indicate porewater bromide concentrations prior to (t=0) and after incubation, respectively.The black line indicates modeled porewater bromide concentrations, with diffusion and bioirrigation (α function).

Figure A. 4 :
Figure A.4: Sensitivity analysis where we varied the maximum growth rate for all microorganisms with different factors (i.e.0.5; 0.75; 0.9; 1.1; 1.25 and 1.5).(A) Porewater depth profiles for the different sensitivity analysis.(B) Cell abundances related to porewater depth profiles in figure A. The baseline scenario (orange) is fitted to the measured data and is the same as in Fig. 2.

Figure A. 5 :
Figure A.5: Sensitivity analysis where we varied the death rate for all microorganisms with different factors (i.e.0.5; 0.75; 0.9; 1.1; 1.25 and 1.5).(A) Porewater depth profiles for the different sensitivity analysis.(B) Cell abundances related to porewater depth profiles in figure A. The baseline scenario (orange) is fitted to the measured data and is the same as in Fig. 2.
ratio of organic matter C/N 7.1 mol mol −1 d Sources: (a) Measured; (b) Model constrained; (c) Reed et al. [2011b]; (d) [Lenstra et al., 2018].Table A.5: Boundary conditions of solids and solutes at the sediment-water interface in the model.Time dependent fluxes of OM α,β,γ , Fe(OH) 3 α,β,γ and MnO 2 α,β at the sediment-water interface are shown in figure A.3.For all chemical species a zero-gradient boundary condition was specified at the bottom of the model domain.Reaction equations implemented in the model.

Table A .
3: Reaction pathways and stoichiometries implemented in the model.