Model

A transient diagenetic model was developed simulating the cycling of 20 chemical species in the sediment (Table 1). The model domain is represented by a one-dimensional grid of 550 cells that captures the interval from the sediment-water interface (SWI) to a depth of 55 cm. Key environmental parameters, such as the porosity, temperature and transport characteristics of the sediment were either taken from literature or constrained from the model (Table 2). The model accounts for biogeochemical reactions (Table 3) and transport processes through burial and bioturbation, as well as molecular/ionic diffusion for dissolved species (Soetaert et al. 1996; Wang and Van Cappellen 1996; Boudreau 1997). The distribution of solid and dissolved species is described by Eqs. 2 and 3, respectively:

(2)

(3)

where C_i and C_j denote concentrations of solid and aqueous species, respectively. ϕ represents the sediment porosity, D_b the bioturbation coefficient, D_j the diffusion coefficients of dissolved species in the sediment, and v and u the advective velocities of the solid and the aqueous phases, respectively (Table 2). The R term denotes the net source or sink of a chemical species from the sum of chemical reactions it participates in.

Table 1. Chemical species included in the model

Species	Notation	Type
Organic mattera		Solid
Oxygen	O2	Solute
Nitrate		Solute
Sulfate		Solute
Iron oxideb	Fe	Solid
Methane	CH4	Solute
Iron	Fe2+	Solute
Ammonium		Solute
Hydrogen sulfidec	ΣH2S	Solute
Elemental sulfur	S0	Solid
Iron monosulfide	FeS	Solid
Pyrite	FeS2	Solid
Phosphate		Solute
Dissolved inorganic carbon	DIC	Solute
Iron-bound phosphorusd	FeP	Solid
Siderite	FeCO3	Solid

• ^a Superscript “α,” “β,” and “γ” denote highly reactive, less reactive, and refractory organic matter, respectively.
• ^b Superscript “α,” “β,” and “γ” denote an increase in the degree of crystallinity.
• ^c ΣH₂S includes H₂S, HS⁻, and S²⁻.
• ^d Iron-bound P accounts for both Fe(OH)₃-bound P (FeP_ox) and P in vivianite (FeP_viv), i.e., FeP = FeP_ox + FeP_viv.

To account for sediment compaction, the porosity is described by

(4)

where x is the depth in the sediment, ϕ₀ the porosity at the top of the sediment, ϕ_∞ the porosity at depth in the sediment, and γ the porosity attenuation factor (e-folding distance, see Table 2).

Table 2. Environmental parameters used in the model

Description	Symbol	Value or expression	Source
Porosity at surface	ϕ0	0.94	b
Porosity at depth	ϕ∞	0.89	b
Porosity e-folding distance	γ	5 cm	a
Sediment density	ρ	2.65 g cm−3	c
Temperature	T	4 °C	b
Salinity	S	6	b
Advective velocity (cm y−1) of solids at depth	v∞		—
Bioturbation coefficient at SWI	Db0	2.7 cm2 y−1	c
Mixed layer depth	ζ	4 cm	a

• Sources: (a) Model constrained; (b) Egger et al. (2015a); (c) Reed et al. (2011). Note: At this site, lugworms were absent (Norkko, personal communication), making bioirrigation a negligible process.

Diffusion coefficients D_j in the sediment were calculated as a function of temperature and salinity (Table 2) following Boudreau (1997) and include a correction for the tortuosity θ in the porous medium (Boudreau 1996), where θ is calculated as follows

(5)

The coefficient D_b was used to account for transport of solids and solutes caused by mixing of the sediment through macrofaunal activity, i.e., bioturbation, and was implemented as

(6)

where D_b0 denotes the bioturbation coefficient near the SWI and ζ the mixed layer depth (Table 2). Following Meysman et al. (2005), the advective velocities v and u of solids and solutes, respectively, caused by sediment accumulation were described by depth-dependent functions to account for changes in the porosity.

The reaction network in the model describes the breakdown of OM, the reoxidation of reduced metabolites, and various mineral formation and dissolution reactions (Table 3). OM degradation occurs through respiratory reactions in which O₂, Fe(OH)₃, , and serve as terminal electron acceptors and through methanogenesis. Following the multi-G approach (Jørgensen 1978; Westrich and Berner 1984), three state variables for OM were used to distinguish between highly reactive (α), less reactive (β), and nonreactive (γ). OM includes organic carbon (C), organic nitrogen (N), and organic phosphorus (P) in a C : N : P ratio of 106 : 15.4 : 1. For Fe(OH)₃ an α, β and γ phase were included to denote variations in crystallinity, which affect the reactivity of Fe(OH)₃ toward OM and ΣH₂S (Tables 3 and 4). The α phase reacts with both OM and ΣH₂S, the β phase persists in the sediment past OM degradation but is susceptible to sulfide scavenging, while the γ phase is not reactive at all. The segregated reactivity of the different Fe(OH)₃ forms was necessary to allow Fe(OH)₃ to persist past the zone of organoclastic Fe(OH)₃ reduction (avoiding reaction with OM) and to be reduced at further depth creating the observed Fe²⁺ accumulation. The model includes oxidation of CH₄ coupled to O₂, , and Fe(OH)₃, where CH₄ is assumed to react with both the α and β phase of Fe(OH)₃. Although MnO₂ has also been shown to be a thermodynamically favorable electron acceptor for AOM (Beal et al. 2009), this study focuses on the potential role of Fe(OH)₃ for AOM only due to the relatively low sedimentary manganese content below the SMTZ at site US5B ( < 40 µmol/µmolg⁻¹) when compared to iron (Egger et al. 2015b) and the likely presence of most of the manganese in the form of manganese carbonate at the relevant depths.

Table 3. Reaction pathways and stoichiometries implemented in the model

Primary redox reactionsa
OMα,β + a O2→ a CO2 + b  + c  + a H2O	R1
OMα,β + 0.8a → a CO2 + b  + c  + 0.6a H2O + 0.4a N2 + 0.8a OH−	R2
OMα,β + 4a  + 4aχα FePox→ a CO2 + b  + (c + 4aχα)  + 3a H2O + 4a Fe2+ + 8a OH−	R3
OMα,β + 0.5a → a CO2 + b  + c  + 0.5a H2S + a OH−	R4
OMα,β→ 0.5a CO2 + b  + c  + 0.5a CH4	R5
Other reactions
2 O2 +   + 2 →  + 2 CO2 + 3 H2O	R6
O2 + 4 Fe2+ + 8  + 2 H2O + 4χα → 4  + 4χα FePox + 8 CO2	R7
2 O2 + FeS→  + Fe2+	R8
7 O2 + 2 FeS2 + 4 OH−→ 4  + 2 Fe2+ + 2 H2O	R9
2 O2 + H2S + 2 →  + 2 CO2 + 2 H2O	R10
2 O2 + CH4→ CO2 + 2 H2O	R11
2  + 2χα,β FeP + H2S + 4 CO2→ 2 Fe2+ + 2χα,β  + S0 + 4  + 2 H2O	R12
Fe2+ + H2S + 2 → FeS + 2 CO2 + 2 H2O	R13
+ CH4 + CO2→ 2  + H2S	R14
CH4 + 8  + 8χα,β FePox→  + 8 Fe2+ + 8χα,β  + 6 H2O + 15 OH−	R15
4 S0 + 2 H2O + 2 OH−→ 3 H2S +	R16
FeS + S0→ FeS2	R17
+ (χα − χβ) FePox→ Fe  + (χα − χβ)	R18
3 Fe2+ + 2 → 2 FePviv	R19
Fe2+ +   + OH−→ FeCO3 + H2O	R20
FeS + H2S→ FeS2 + H2	R21
+ (χα,βmax − χα,β) →  + FePox	R22

• ^a Organic matter (OM) is of the form (CH₂O)_a( )_b( )_c. The “a,” “b,” and “c” refer to the C : N : P ratio in organic matter. χ^α,β,γ denote the P : Fe ratio of , which were estimated based on the sedimentary profiles of and FeP_ox (χ^α = 0.18, χ^β = 0.1 and χ^γ = 0.06).

Table 4. Reaction equations implemented in the model

• The superscript i denotes the highly reactive α and less reactive β fractions of organic matter; superscript j denotes the less crystalline α and more crystalline β fractions of iron oxyhydroxides.

Dissolved inorganic carbon (DIC) and porewater phosphate ( ) were included in the model to allow formation of the reduced Fe minerals siderite and vivianite, which act as key sinks for porewater Fe²⁺ (Tables 3 and 4). DIC was calculated as the sum of the carbon in CO₂ and that was produced or consumed in the modeled reactions (Table 3). Dissolved was allowed to adsorb to existing (see R22, Table 3). It was assumed that the initial P:Fe ratios (χ of arriving at the sediment surface (χ^α = 0.18, χ^β = 0.1) were lower than an estimated maximum binding capacity of under high porewater concentrations ( = 0.25, = 0.14). The boundary conditions at the SWI were specified as variable, time-dependent fluxes for total sediment accumulation and inputs of OM and Fe(OH)₃ (see model application section), and as constant fluxes for all other solids and fixed concentrations for dissolved species (Table 5). For all species a zero gradient boundary condition was set at the bottom of the model domain. Bottom-water O₂ and concentrations were imposed, while the concentrations of reduced species were set to zero. All CH₄ was assumed to be produced within the model domain.

Table 5. Boundary conditions in the model. See Fig. 1 for the time-dependent OM and Fe(OH)₃ fluxes at the SWI. Flux of OM in units of C. The flux of FeCO₃ was adjusted to the time-dependent sedimentation rate to allow for a constant background concentration of FeCO₃ of 120 µmol g⁻¹ with depth in the sediment

Solids	Flux at SWI
	0 mol m−2 y−1
	0 mol m−2 y−1
	0 mol m−2 y−1
	0.1 mol m−2 y−1
	0.02–0.76 mol m−2 y−1

Solutes	Bottom-water concentration
	0.18 mol m−3
	10 × 10−3 mol m−3
	4.8 mol m−3
	0 mol m−3
	0 mol m−3
	0 mol m−3
	10 × 10−3 mol m−3
CDIC	0.15 mol m−3

• For all chemical species a zero-gradient boundary condition was specified at the bottom of the model domain.

The model code was written in R using the marelac (Soetaert et al. 2010) geochemical dataset package and the ReacTran (Soetaert and Meysman 2012) package to calculate the transport in porous media. The ordinary differential equations were solved with the lsoda integrator algorithm (Hindmarsh 1983; Petzold 1983).

Model application and sensitivity analyses

The model was applied to a site in the Bothnian Sea, which is a brackish coastal basin with an average bottom water salinity of 5–6. The site considered, US5B, is located in the deepest part of the basin at 214 m water depth and is characterized by fine-grained, organic-rich sediments. The recent upward displacement of the SMTZ at this site and subsequent fixation of its position near the sediment surface reflects strong temporal changes in OM input over the past decades that have been attributed to a period of eutrophication followed by a recovery phase (Slomp et al. 2013; Egger et al. 2015a, 2015b). These observations are in accordance with long-term monitoring studies based on inorganic nutrients, water transparency and chlorophyll a in the Bothnian Sea that reveal substantial changes in anthropogenic loading of land-derived nutrients over recent decades (Fleming-Lehtinen et al. 2008; Savchuk et al. 2008; Fleming-Lehtinen and Laamanen 2012). Sedimentation rates at this site have previously been estimated to be high (up to 1.9 cm y⁻¹) and have been shown to vary in space and time (Mattila et al. 2006; Egger et al. 2015b). Repeated sampling campaigns indicate that the depth of the peak in solid phase sulfur associated with the SMTZ did not change position relative to the sediment-water interface between 2008 and 2012, suggesting a recent decline in sedimentation rate (Egger et al. 2015b).

The model was calibrated to depth profiles of sediment organic carbon, porosity, solid sulfur (including FeS and FeS₂), Fe(OH)₃, siderite, iron-bound phosphorus (P) (which represents both P associated with Fe(OH)₃ (FeP_ox) and P in vivianite (FeP_viv)) and porewater profiles of CH₄, , , ΣH₂S, Fe²⁺, and that were collected at site US5B in August 2012 (Egger et al. 2015a, 2015b). Reaction parameters were mostly taken from the literature (Table 6). Where parameters were not available or no fit to the data could be obtained with existing parameter ranges, reaction parameters were constrained using the model. The relevant parameters are , k₈, k₁₀, k₁₂, k₁₄, k₁₅, k₁₆, and k₁₇ (Table 6) and thus include the rate constant of AOM coupled to Fe(OH)₃ (k₁₀). Rate constants for the reaction of ΣH₂S with and differed slightly ( and , Table 6). This is in accordance with differences in reactivity of iron oxide minerals of varying crystallinity to ΣH₂S as reported in the literature (Canfield et al. 1992; Poulton et al. 2004).

Table 6. Reaction parameters used in the model

Parameter	Value	Units	Source	Range given by source (values used by source)
	1.62	y−1	a	1.62
	0.0086	y−1	a	0.0086
	20	µM	b	1–30, (20)
	4	µM	b	4–80, (2)
	65	µmol g−1	b	65–100, (65, 100)
	1.6	mM	b	1.6, (1.6)
	5000	mM−1 y−1	b	10,000, (5000)
	140,000	mM−1 y−1	b	140,000, (140,000)
	300	mM−1 y−1	b	300, (300)
	1	mM−1 y−1	b	1
	160	mM−1 y−1	b	≥160, (160)
	107	mM−1 y−1	b	107, (107)
	8	mM−1 y−1	b	≤100, (8)
	5.6	mM−1 y−1	b	≤100, (8)
k8	1000	mM−1 y−1		Model constrained
k9	120	mM−1 y−1	b	10, (10)
k10	0.0074	mM−1 y−1		Model constrained
k11	3	y−1	c	3, (3)
k12	0.001	mM−1 y−1		Model constrained
k13	0.6	y−1	c	0.6, (0.6)
k14	0.51	mM−1 y−1		Model constrained
k15	0. 39	mM−1 y−1		Model constrained
k16	0.001	mM−1 y−1		Model constrained
k17	0.047	mM−1 y−1		Model constrained

• Sources: (a) Reed et al. (2011); (b) Wang and Van Cappellen (1996); (c) Berg et al. (2003).

The model was run to steady state for 110 y and temporal changes in OM and Fe(OH)₃ loading and sediment accumulation rates were then adjusted to fit the data (Fig. 1). The resulting scenario for this site is as follows: Sedimentation rates were initially high at 0.50 g cm⁻² y⁻¹ around 1970, reached a maximum of 0.63 g cm⁻² y⁻¹ between the years 1992 and 2002 and then rapidly decreased during the last 10 y of the simulation (0.03 g cm⁻² y⁻¹). Before the onset of anthropogenic eutrophication, the influx of organic carbon was 6.2 mol m⁻² y⁻¹ (Fig. 1), which falls within the range of 1.1 to 8.2 mol m⁻² y⁻¹ estimated previously for Bothnian Sea sediments (Algesten et al. 2006). The OM loading then began to increase, followed by a large OM pulse during the late 1990s and early 2000s and lower inputs during the last 10 y (Fig. 1). Fe(OH)₃ fluxes increased during the late 1990s and early 2000s and subsequently decreased during the last 10 y (Fig. 1), consistent with changes in OM loading.

The transient evolution of OM loading reported here compares well to qualitative trends reported in previous studies on eutrophication in the Bothnian Sea (Fleming-Lehtinen et al. 2008; Fleming-Lehtinen and Laamanen 2012). No quantitative data are available to further constrain our scenario. However, we note that, based on the location of our study site in a deep basin, it is likely that major shifts in lateral input of OM and FeOH₃ and sediment focusing play a key role in controlling the strong temporal changes. Eutrophication-induced temporal changes in lateral transfer of iron from shallow areas to deep basins have been well-described for other Baltic Sea basins (Fehr et al. 2010; Jilbert and Slomp 2013; Lenz et al. 2015).

The transient scenario described above (“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 setting the OM and Fe(OH)₃ loadings, the O₂ and bottom-water concentrations, and the rate parameter of Fe-AOM (k₁₀, Table 6) to 80 and 120% of the original values (specified in the baseline scenario). To assess the role of environmental forcings in more detail, the OM loading was also changed simultaneously with the Fe(OH)₃ loading and bottom-water concentration.

Coastal eutrophication and Fe-AOM

The model results presented in this study confirm that increased inputs of OM to the sediment can induce an upward vertical migration of the SMTZ in coastal sediments, as suggested previously for this site (Slomp et al. 2013; Egger et al. 2015a, 2015b). Our model results show that increased rates of reduction and methanogenesis linked to enhanced inputs of OM (Fig. 1) can initiate a relatively rapid upward movement of the SMTZ (Fig. 3) thus bringing it close to the sediment surface. Reduced inputs of OM over the last 10 y have now brought the upward displacement of the SMTZ to a halt. The fixation of the SMTZ at its current position has resulted in the accumulation of authigenic FeS_x minerals in a relatively narrow zone over the past decade (Fig. 3).

High sediment accumulation rates and Fe(OH)₃ loading have allowed a significant portion of the Fe(OH)₃ to become buried below the newly established SMTZ. The preservation of reducible Fe(OH)₃ below the zone of reduction has stimulated Fe-AOM, which has led to the build up of dissolved Fe²⁺ in the porewater below the SMTZ. Our results support earlier suggestions by Egger et al. (2015a) that anthropogenic eutrophication of coastal systems may increase the importance of Fe-AOM in coastal environments by inducing nonsteady state diagenesis.

Controls on the Fe-AOM reaction

Besides the quest for conclusive evidence that AOM is truly coupled to Fe(OH)₃ reduction in natural environments and for the microorganisms carrying out the reaction, there is also the question what processes could potentially control in-situ rates of Fe-AOM in sediments. While laboratory studies suggest that the availability Fe(OH)₃ can boost the rate of AOM (Beal et al. 2009; Sivan et al. 2014; Egger et al. 2015a) and that the addition of multiple electron acceptors may decrease the AOM rate (Segarra et al. 2013), the potential response of Fe-AOM to variations in environmental conditions in the field has so far not yet been described.

Here, we calibrated a transient sedimentary model to data from a coastal brackish basin where various lines of evidence suggest that Fe-AOM occurs (Egger et al. 2015a). The model assumes a direct coupling between Fe(OH)₃ reduction and CH₄ oxidation as proposed by Beal et al. (2009) (see Eq. 1), rather than indirect iron stimulated SO₄-AOM (Sivan et al. 2014) because ΣH₂S and are undetectable below the SMTZ. Starting from this assumption, variations in key environmental forcings can be implemented in the model to discern potential controls on Fe-AOM.

The model assumes three electron acceptors which can react with the available CH₄: O₂, , and Fe(OH)₃. Our results show that Fe-AOM is mostly restricted to parts of the sediment where other electron acceptors than Fe(OH)₃ are absent. The oxidation of CH₄ coupled to O₂ is by far the most favorable kinetically (k₅ in comparison to k₉ and k₁₀, Table 6), but is constrained to the upper millimeters of sediment where O₂ is present. The rate constant for SO₄-AOM is five orders of magnitude lower than that of aerobic CH₄ oxidation, but occurs over a larger depth interval because is present throughout the upper 10 cm. In the model, the kinetics of Fe-AOM follow the same bimolecular rate law as the other CH₄ oxidation pathways. CH₄ in the top is oxidized by O₂ and , while the Fe-AOM reaction is only optimal at depths where is depleted (compare Figs. 3, 4).

Enhanced bottom-water concentrations significantly reduce Fe-AOM rates as less CH₄ is available for Fe-AOM. This is due to enhanced organoclastic reduction, which suppresses methanogenesis. For example, a three-fold increase of the concentration (corresponding to 14.5 mM ) leads to a CH₄ oxidation rate that is only 46% of that in the baseline scenario and a Fe-AOM rate near zero (Table 7).

The availability of Fe(OH)₃ for AOM depends on the fraction of the sedimentary Fe(OH)₃ that can react with OM and ΣH₂S. A higher OM loading leads to a steeper gradient (Fig. 6a) and thus more diffusion of into the sediment. This enhances ΣH₂S production and subsequent Fe(OH)₃ scavenging by ΣH₂S. Therefore, systems that receive high reactive OM fluxes but no associated high Fe(OH)₃ flux, will not be suitable environments for Fe-AOM in the surface sediment. The accumulation of Fe²⁺ can only occur when the Fe(OH)₃ loading is high enough to constrain the ΣH₂S to the upper part of the sediment. Figure 8a shows that almost no Fe-AOM was observed in simulations with a 50% lower Fe(OH)₃ loading.

The impact of changes in Fe(OH)₃ loading on OM degradation is generally limited as the organoclastic reduction rate is much larger (Table 7). As the methanogenesis rate is dependent on the total sum of the other OM degradation pathways (Table 4), organoclastic Fe(OH)₃ reduction has a fairly small impact on the in-situ CH₄ production. In line with Egger et al. (2015a), we find that Fe-AOM only accounts for a small fraction of the total CH₄ oxidation rate, and therefore has a relatively small effect on the CH₄ profile.

Signatures of Fe-AOM

As previously outlined, direct identification of the Fe-AOM process in natural environments is difficult because of the lack of unequivocal geochemical signatures. It is therefore useful to test how the geochemical profiles of the coupled iron-sulfur-carbon cycles can serve as signatures for the presence of Fe-AOM. By simulating a range of environments, we have constructed a group of theoretical sedimentary settings that we can use to assess how the presence of Fe-AOM is reflected in the geochemical profiles of its associated species.

The usefulness of CH₄ profiles for identifying Fe-AOM is limited, as the simulation outcomes show that CH₄ is not very sensitive to different intensities of Fe-AOM. In the simulation with a 20% higher Fe(OH)₃ loading CH₄ concentrations were slightly lower than obtained in calculations with the baseline scenario. This was not only caused by the increased Fe-AOM rates, but also by the higher rates of organoclastic Fe(OH)₃ reduction (Table 7). The latter aspect affects the partitioning of OM degradation among electron acceptors and slightly suppresses methanogenesis. Fe-AOM leads to the production of large quantities of Fe²⁺, which can lead to the accumulation of Fe²⁺ at depth. However, to interpret the impact of Fe-AOM on Fe²⁺ profiles, the sinks of Fe²⁺ need to be taken into account (Fig. 7). Most dissolved Fe²⁺ is removed through FeS_x formation at a rate that is dependent on multiple factors. The availability of ΣH₂S is controlled by the reduction rate, which in turn is dependent on the bottom-water concentration and the OM loading. The depth and breadth of the Fe-AOM front is consequential, as Fe²⁺ produced higher up in the sediment is more likely to react with ΣH₂S. Also a significant part of the Fe²⁺ is stored in carbonate and phosphate minerals, i.e., siderite and vivianite (Fig. 7). In the simulations for the baseline scenario, Fe²⁺ does not accumulate at depth before 1982 (Fig. 3), while the total Fe-AOM rates over depth are significant (Fig. 4). Therefore, Fe²⁺ concentration profiles are not a robust signature of Fe-AOM.

Conversely, Fe²⁺ accumulation at depth, may signal the occurrence of Fe-AOM, as it can explain the concurrent production of CH₄ and Fe²⁺. Such a simultaneous production is hard to explain otherwise if the assumption is correct that methanogens are generally outcompeted for common substrates by organoclastic Fe-reducers (Egger et al. 2015a and references therein). In the simulations, Fe-AOM is the only process that can produce Fe²⁺ at depth. The other Fe²⁺ producing reactions are dependent on either O₂, ∑H₂S, or the alpha phase of Fe(OH)₃, which are only available in the upper few centimeters of sediment (see Figs. 3, 4). To achieve deep Fe²⁺ production through Fe-AOM in our scenario, the Fe(OH)₃ needs to penetrate to sufficient depth. This requires a fraction of the Fe(OH)₃ pool to be unreactive toward OM, which is in line with the characteristics of the Fe phase. Mass balance calculations show that the Fe phase is almost exclusively involved in the Fe-AOM reaction (Fig. 7).

Environmental settings

Our results suggest that Fe-AOM is promoted in brackish coastal marine sediments with relatively high OM fluxes and Fe(OH)₃ loadings. Estuarine environments such as Chesapeake Bay, Long Island Sound and the Scheldt estuary all exhibit these conditions. Thus they are all potential candidate sites for investigating the occurrence and scale of AOM coupled to metal oxide reduction. From our sensitivity analyses we see that higher bottom-water concentrations significantly reduce the probability of the occurrence and the intensity of Fe-AOM (Fig. 8). This indicates that the occurrence of this process in surface sediments is, in principle, restricted to low salinity environments. However, Fe-AOM can also occur in deep-sea sediments with typical bottom-water concentrations for marine settings (∼28 mM) if there is a low sulfate zone at depth, as shown by Riedinger et al. (2014). At their site, CH₄ diffuses up from greater depth and Fe-AOM is believed to occur at 6–8 m depth below the SWI. In our model simulations, higher bottom-water concentrations reduce Fe-AOM rates primarily by suppressing methanogenesis. However, higher bottom-water concentrations will be less consequential for Fe-AOM when Fe(OH)₃ is buried to greater depth (outside our model domain) under oligotrophic conditions. At depth, the CH₄ can be formed in situ when the system shifts from oligotrophic to eutrophic conditions, which leads to an upward movement of the SMTZ and methanogenic zone. There can also be an ex situ deep CH₄ source, which explains why, for example, cold seep environments are also key settings for Fe-AOM (Beal et al. 2009; Wankel et al. 2012).

A necessary precondition for Fe-AOM, is the presence of sufficient amounts of Fe(OH)₃ to react with the available CH₄, as there is a 8 : 1 stoichiometric ratio of Fe(OH)₃ to CH₄ in the net reaction. We show that relatively high fluxes of Fe(OH)₃ are necessary to provide enough material for the reaction to occur. Furthermore, the simulations indicate that Fe-AOM is favored in the presence of less reactive iron phases that can penetrate deeper into the sediment. Settings, which receive high amounts of Fe(OH)₃, such as coastal systems near large rivers (Poulton and Raiswell 2002), are therefore also candidates for hosting the reaction.

The impact of OM loading on Fe-AOM is complicated as it affects the depth penetration of , the reduction rate, the rate of methanogenesis, and the persistence of Fe(OH)₃ at depth. The sensitivity analyses show that a wide range of OM loadings may lead to substantial Fe-AOM rates (Fig. 8), but that too much and too little OM can have a detrimental effect on the process.

In summary, the optimal environments for the Fe-AOM process are sedimentary settings which are depositional and receive sufficient OM loading (much of which can be refractory), receiving high influxes of Fe(OH)₃ (preferably of a more refractory nature), and which have sufficiently low salinity. Coastal areas, with brackish water columns and a connection to rivers are thus ideal environments for the Fe-AOM process. In more marine settings, Fe-AOM is far less likely, but may still occur, provided there is an ex-situ source of CH₄ and enough Fe(OH)₃ to fuel the reaction.

Sediment-water exchange fluxes of Fe²⁺

Continued shuttling of iron through repeated reduction-oxidation cycles and lateral iron transport on oxic and hypoxic continental shelves depends on the replenishment of iron from near-coastal regions. Our results suggest that, besides ∑H₂S-driven and organoclastic Fe(OH)₃ reduction, Fe-AOM could contribute to benthic release of Fe²⁺ and thus to the supply of iron from estuarine sediments to continental shelves. In the model result for 2012, the Fe²⁺ efflux to the overlying water is negligible. However, the results for the sensitivity analyses show that increased OM loadings can enhance exchange fluxes of Fe²⁺ to the overlying water (Fig. 8), with values of up to 0.04 mol m⁻² y⁻¹ being reached. These Fe²⁺ fluxes fall within the low end of the typical range estimated for continental margin sediments overlain by oxic waters (0 to ca. 0.12 mol m⁻² y⁻¹; e.g., Severmann et al. 2010; Dale et al. 2015). The increased release of Fe²⁺ on increased OM loading is the result of a shift in the zones of Fe-AOM and organoclastic Fe(OH)₃ reduction toward the SWI (Fig. 4). Increased bottom-water concentrations hinder Fe²⁺ effluxes due to lower rates of Fe(OH)₃ reduction (Fig. 8) and more sulfidic scavenging of Fe²⁺. Very low bottom-water concentrations and high OM loadings may lead to Fe²⁺ effluxes > 0.12 mol m⁻² y⁻¹, which suggests that Fe-AOM has the potential to greatly impact iron dynamics in low-salinity coastal environments.