Iron and sulfur cycling in the cGENIE.muffin Earth system model (v0.9.21)

The coupled biogeochemical cycles of iron and sulfur are central to the long-term biogeochemical evolution of Earth's oceans. For instance, before the development of a persistently oxygenated deep ocean, the ocean interior likely alternated between states buffered by reduced sulfur (“euxinic”) and buffered by reduced iron (“ferruginous”), with important implications for the cycles and hence bioavailability of dissolved iron (and phosphate). Even after atmospheric oxygen concentrations rose to modern-like values, the ocean episodically continued to develop regions of euxinic or ferruginous conditions, such as those associated with past key intervals of organic carbon deposition (e.g. during the Cretaceous) and extinction events (e.g. at the Permian–Triassic boundary). A better understanding of the cycling of iron and sulfur in an anoxic ocean, how geochemical patterns in the ocean relate to the available spatially heterogeneous geological observations, and quantification of the feedback strengths between nutrient cycling, biological productivity, and ocean redox requires a spatially resolved representation of ocean circulation together with an extended set of (bio)geochemical reactions.

Here, we extend the “muffin” release of the intermediate-complexity Earth system model cGENIE to now include an anoxic iron and sulfur cycle (expanding the existing oxic iron and sulfur cycles), enabling the model to simulate ferruginous and euxinic redox states as well as the precipitation of reduced iron and sulfur minerals (pyrite, siderite, greenalite) and attendant iron and sulfur isotope signatures, which we describe in full. Because tests against present-day (oxic) ocean iron cycling exercises only a small part of the new code, we use an idealized ocean configuration to explore model sensitivity across a selection of key parameters. We also present the spatial patterns of concentrations and δ⁵⁶Fe and δ³⁴S isotope signatures of both dissolved and solid-phase Fe and S species in an anoxic ocean as an example application. Our sensitivity analyses show that the first-order results of the model are relatively robust against the choice of kinetic parameter values within the Fe–S system and that simulated concentrations and reaction rates are comparable to those observed in process analogues for ancient oceans (i.e. anoxic lakes). Future model developments will address sedimentary recycling and benthic iron fluxes back to the water column, together with the coupling of nutrient (in particular phosphate) cycling to the iron cycle.

1 Introduction

The biogeochemical cycles of iron and sulfur are tightly coupled in the marine environment and play fundamental roles in the evolution and functioning of the Earth system (Raiswell and Canfield, 2012). In their main oxidized states, both sulfur (in the form of sulfate; SO${}_{\mathrm{4}}^{\mathrm{2}-}$) and iron (in the form of iron (oxyhydr)oxides; FeOOH) are important electron acceptors in the oxidation of organic matter in anoxic environments such as marine sediments or oxygen-deficient water columns (e.g. Black Sea, stratified lakes) (Thamdrup, 2000; Crowe et al., 2008; Raiswell and Canfield, 2012). In metabolizing organic matter, microbial reduction of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and FeOOH produces reduced sulfide (H₂S) and ferrous iron (Fe²⁺), respectively. When present at the same location, H₂S and Fe²⁺ combine into iron monosulfides and eventually pyrite (FeS₂) (Rickard, 1997, 2006), the burial of which couples the short-term surface cycles of Fe and S with their long-term geological cycles (Berner, 1989). Depending on the relative ocean inventory of H₂S vs. Fe²⁺, the precipitation of FeS₂ can lead to an anoxic water body becoming either iron-rich (“ferruginous”) or sulfide-rich (“euxinic”) (Canfield, 1998; Poulton and Canfield, 2011), with ${\mathrm{H}}_{\mathrm{2}}\mathrm{S}:{\mathrm{Fe}}^{\mathrm{2}+}$ ratios greater than 2:1 (the S:Fe ratio in FeS₂) promoting euxinic conditions and lower ratios promoting ferruginous conditions (Poulton and Canfield, 2011). For most of Earth's history, the ocean interior is thought to have been predominantly anoxic (Fig. 1; Lyons et al., 2014), which implies that reduced forms of iron and sulfur would have dominated the marine redox landscape (Poulton and Canfield, 2011; Raiswell and Canfield, 2012).

Whether an anoxic water body becomes ferruginous or euxinic can have significant impacts on the availability of nutrients, with ferruginous conditions potentially leading to phosphate limitation and euxinic conditions potentially leading to depletion of key biological trace elements (Van Cappellen and Ingall, 1996; Bjerrum and Canfield, 2002; Reinhard et al., 2013, 2017; Wallmann et al., 2019). Furthermore, before the advent of oxygenic photosynthesis, the productivity of marine ecosystems was likely, at least partly, fuelled by the oxidation of H₂S or Fe²⁺ to their oxidized counterparts (Kharecha et al., 2005; Canfield et al., 2006; Ozaki et al., 2018; Thompson et al., 2019). As a result, the long-term evolution of the Earth system and the structure of marine ecosystems are closely tied to the evolution of the biogeochemical cycles of iron and sulfur. The ability to simulate the evolution of these cycles and geochemical distributions in the ocean in models hence becomes key to better understanding the early evolution of microbial ecosystems.

During the early stages of Earth's history (i.e. Archean to mid-Proterozoic), the ocean was rich in Fe²⁺ (Fig. 1), which enabled the extensive deposition of banded iron formations (BIFs) and ferruginous shales during those eons (Bekker et al., 2010; Planavsky et al., 2011; Sperling et al., 2015; Konhauser et al., 2017). BIFs are rare in sediments deposited after 1.8 billion years ago, which was initially hypothesized to reflect a transition to oxygen-rich bottom waters that prevented the build-up of soluble iron by removing it as insoluble iron oxides (Cloud, 1972; Holland, 1984). In a seminal paper, Canfield (1998) suggested that the abundance of atmospheric oxygen during the Proterozoic (which was at most ∼10 % of today; Canfield and Teske, 1996; Lyons et al., 2014) was too low to oxygenate the deeper waters of the ocean. Instead, he proposed that the disappearance of BIFs was driven by an increase in oceanic sulfate concentrations, allowing sulfide (produced by microbial sulfate reduction) to remove reduced iron from solution by the precipitation of FeS₂ (Canfield, 1998). Hence, this hypothesis implied that the deep ocean was euxinic for most of the Proterozoic. Since then, a large wealth of geochemical proxy data have been collected, aided by the development of a sequential Fe extraction scheme that helps to differentiate between oxic, euxinic, and ferruginous conditions by determining seven operationally defined iron mineral pools (Poulton and Canfield, 2005, 2011). These paleo-reconstructions of Archean and Proterozoic deposits helped to further refine the spatial and temporal history of ocean redox. In contrast to the earlier view, the Proterozoic eon appears to have been dominated by ferruginous conditions (Canfield et al., 2008; Planavsky et al., 2011; Poulton and Canfield, 2011; Guilbaud et al., 2015), likely interspersed with euxinic excursions along the shelf (Canfield et al., 2008; Poulton et al., 2010; Raiswell and Canfield, 2012) (Fig. 1). However, to date, these redox landscapes have been largely qualitative, since field-based observations are restricted by the limited number of available unaltered deposits from a certain point in time and sample acquisition for a depth transect is an enormously labour-intensive task (see e.g. Poulton et al., 2010). Because of these limitations, both the spatial redox pattern and the exact nature of ferruginous versus euxinic conditions (i.e. the concentrations of Fe²⁺ or H₂S) are still relatively unconstrained. The uncertainty with respect to the ocean redox conditions propagates into hypotheses on nutrient availability and ecosystem evolution (see e.g. Reinhard et al., 2017).

Figure 1 First-order evolution of Earth's ocean redox landscape. Insets show conceptual models of the spatial redox structure of the ocean at certain points in Earth's history. Based on Poulton et al. (2010), Poulton and Canfield (2011), and Raiswell and Canfield (2012). After van de Velde et al. (2020b). See text for details.

Whether the ocean interior is ferruginous or euxinic can significantly impact the supply of essential nutrients to surface ocean environments and thus nutrient availability for primary producers (Van Cappellen and Ingall, 1996; Bjerrum and Canfield, 2002; Guilbaud et al., 2020). For example, iron oxides (FeOOH) can efficiently scavenge phosphate (${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$) from a water column, which has been suggested to limit oceanic phosphate availability in the past (Berner, 1973; Van Cappellen and Ingall, 1996; Bjerrum and Canfield, 2002; Jones et al., 2015). In contrast, a euxinic water column could potentially induce trace metal limitations by titrating out essential trace elements like iron or molybdenum (Reinhard et al., 2013; Wallmann et al., 2019). Indeed, lower nutrient availability – specifically phosphate – has been suggested to have limited primary productivity for the majority of the Proterozoic (Bjerrum and Canfield, 2002; Reinhard et al., 2017; Ozaki et al., 2019). At the same time, nutrient availability is thought to have been critical in shaping and driving eukaryote evolution and proliferation (Brocks et al., 2017; Reinhard et al., 2020b). However, many of the most significant innovations in the history of Earth's biosphere have likely occurred in specific sites of the ocean that did not reflect the average state of the ocean as a whole (Nisbet and Sleep, 2001). Hence, the ability to reconstruct iron and sulfur cycling in a spatially explicit way is critical for exploring the relationships between biospheric evolution and changes in ocean redox.

None of the currently available suite of global models that explicitly represent biogeochemical cycling under low-oxygen marine environmental conditions include an extensive treatment of biogeochemical iron cycling (Ozaki et al., 2011; Laakso and Schrag, 2014; Hülse et al., 2017; Lenton et al., 2018; Reinhard et al., 2020a). Consequently, many essential interactions between the iron cycle and other elemental cycles (e.g. sulfur burial via pyrite precipitation) are abstracted and parameterized using techniques that may or may not be mechanistically robust across a range of scenarios. Moreover, most of the ocean models used to simulate the early stages of Earth's history are essentially box or one-dimensional ocean models (Ozaki et al., 2011; Laakso and Schrag, 2014; Lenton et al., 2018) and are unable to resolve the spatial patterns necessary to start contrasting with the geological record (although some two- or three-box models attempt to distinguish between the coastal zone and the open ocean; Laakso and Schrag, 2019; Thompson et al., 2019; Alcott et al., 2019). Indeed, previous simulations of past oceans with three-dimensional ocean models have already indicated the importance of spatial patterns in ocean redox (e.g. Olson et al., 2013), specifically when considering habitability for complex eukaryotic life (Reinhard et al., 2016, 2020b).

In this paper, we present the development of a coupled anoxic oceanic iron and sulfur cycle embedded within the “muffin” release of the carbon-centric Grid-ENabled Integrated Earth system model, cGENIE (note that the oxic cycle of both Fe and S already exists in previous versions; described in Tagliabue et al., 2016, and Ridgwell et al., 2007, respectively). The aim is to extend the functionality of cGENIE into regimes in which the ocean interior is pervasively anoxic, including much of Earth's Precambrian history and periods of significant perturbation to ocean redox during the Paleozoic and Mesozoic (during so-called “ocean anoxic events”; OAEs) when coupled iron–sulfur cycling dominated biogeochemical interactions in the ocean interior. Our extension explicitly accounts for the formation, burial, and isotopic compositions of key mineral phases used in paleo-environmental reconstructions – specifically iron oxide (FeOOH), siderite (FeCO₃), greenalite (Fe₃Si₂O₅(OH)₄), and pyrite (FeS₂) (Poulton and Canfield, 2005, 2011) – which allows quantitative comparison with available data (e.g. Rouxel et al., 2005; Heard and Dauphas, 2020). In the next section (Sect. 2), we briefly discuss the cGENIE model framework, focusing on the features that are most relevant for our purpose of modelling pervasively anoxic oceans. Section 3 describes the included iron–sulfur reactions and the reasoning behind the chosen parameterizations. In Sect. 4 we discuss a modern configuration, an example configuration, and a series of sensitivity experiments for the chosen parameterization. Finally, in Sect. 5 we discuss model limitations and potential future developments.

2 The cGENIE.muffin Earth system model framework

cGEnIE is an Earth system model of intermediate complexity (EMIC) which comprises a modular framework that incorporates different components of the Earth system, including ocean circulation and biogeochemical cycling, ocean–atmosphere and ocean–sediment exchange, and the long-term (geological) cycle of carbon and various solid-Earth-derived tracers (Ridgwell et al., 2007; Ridgwell and Hargreaves, 2007; Colbourn et al., 2013; Adloff et al., 2020). Here, we use the current “muffin” release that encompasses a range of developments and/or additions in the representation of temperature-dependent metabolic processes in the ocean (Crichton et al., 2021), ocean–atmosphere cycling of methane (Reinhard et al., 2020a), marine ecosystems (Ward et al., 2018), organic matter preservation and burial in marine sediments (Hülse et al., 2018), and geological cycles of weathering-relevant trace metals and isotopes (Adloff et al., 2020).

The climate component in cGENIE – C-GOLDSTEIN – consists of a 2-D energy–moisture balance model of the atmosphere coupled to a reduced physics (frictional geostrophic) 3-D ocean circulation plus dynamic–thermodynamic sea ice model (see Edwards and Marsh, 2005; Marsh et al., 2011, for full descriptions). In addition to the simplified atmosphere, to further facilitate the simulation of a relatively large number of interacting gaseous, dissolved, and solid tracers across atmosphere, ocean, and marine sediment (and land surface), cGENIE can be configured with a highly reduced spatial and temporal resolution relative to most high-resolution ocean general circulation models (the default temporal resolution, which we use here, requires 48 time steps per year in solving ocean circulation). While this precludes exploration of very detailed spatial patterns, it does provide a flexibility which is not available in more high-resolution models, and the relatively short runtime of cGENIE (around 1 d per 10 000 model years on a single CPU core) allows us to run many different model experiments and carry out comprehensive parameter sweeps and sensitivity analysis (and hence parameter tuning and model calibration). This aligns with our ultimate aim here, which is to explore ocean biogeochemistry during periods of the Mesozoic (>65 million years ago), Paleozoic (>250 million years ago), and Precambrian (>540 million years ago). Most of these changes likely occurred on timescales exceeding 10 000 years, and at present we have virtually no information with respect to seasonality or detailed spatial variability for many of these intervals (see e.g. Poulton et al., 2010; Guilbaud et al., 2015). In addition, key boundary conditions and parameter values for these periods of Earth's history are often poorly constrained, necessitating large model ensembles in order to adequately assess the robustness of any given result.

2.1 Continental configuration and climatology

For the purpose of this study – implementing and characterizing the coupled cycling of iron and sulfur in an anoxic ocean – we adopt a deliberately idealized model configuration. We configure the ocean model on a 18×18 equal-area horizontal grid with 16 logarithmically spaced z-coordinate levels (Fig. 2). The horizontal grid is uniform in longitude (20^∘ resolution) and uniform in the sine of latitude (∼3.2^∘ at the Equator to 19.2^∘ near the poles) (Fig. 2a). The layer thickness in the vertical grid increases from 80.8 m at the surface to 765 m at the deepest layer (Fig. 2b). We adopt a “ridge world” set-up, with a thin strip of land connecting the North and South Pole (Fig. 2c), following Ferreira et al. (2010), which creates a single ocean basin with no circumpolar current (a little akin to the plate configuration prevailing during the late Permian). We apply idealized boundary conditions of zonally averaged wind stress and speed, plus a zonally averaged planetary albedo, all following Vervoort et al. (2021) (Fig. 3a–c). The solar constant is set to modern (1368.0 W m⁻²). It is worth noting that our model set-up is entirely abyssal, while in reality the continental slopes and shelves are important for the biogeochemical cycling of Fe and S. We have chosen our idealized bathymetry and continental configuration as it generates a relatively simple ocean circulation and thus facilitates the interpretation of model output (see below). Our choice allows us to clearly illustrate the dependence of model output on the choice of parameters, whereas more elaborate continental configurations would introduce more spatial complexity and potentially obscure model sensitivity to parameter selection.

The physics parameters controlling the model climatology all follow the 16-level ocean-based configuration assessed by Cao et al. (2009) (Table S1 in their paper), with the exception of the parameterization controlling ocean mixing. When using the standard implementation of isoneutral diffusion and eddy-induced advection (Edwards and Marsh, 2005; Marsh et al., 2011), sharp vertical redox gradients simulated under extreme redox and high dissolved iron conditions resulted in unacceptably negative tracer concentrations at depth, particularly for dissolved iron. We hence disabled this parameterization, reverting to the original unadjusted horizontal + vertical diffusion physics configuration of Edwards and Shepherd (2002). We tested the modern cGENIE configuration of Cao et al. (2009) for both ocean mixing parameterizations and found only a minimal difference in the large-scale ocean circulation (i.e. a ∼2 Sv stronger Atlantic meridional overturning circulation in the unadjusted parameterization).

Figure 2 Schematic depiction of the cGENIE grid used for our simulations. The grid is $\mathrm{18}\times \mathrm{18}\times \mathrm{16}$, uniform in longitude, and uniform in the sine of latitude. The layer thickness in the vertical grid increases from 80.8 m at the surface to 765 m at the deepest layer. (a) Top view of the cGENIE grid with an indication of the layer numbers. (b) Side view of the cGENIE grid with an indication of the layer numbers. (c) Schematic picture of the continental configuration. “Ridge world” has one continent which runs from pole to pole and extends from 0 to 20^∘ E (i.e. the width of a longitudinal grid cell). The ocean is 5000 m deep everywhere.

Figure 3 (a) Albedo, (b) topography, (c) sea surface temperature and wind stress, and (d, e) ocean circulation patterns for all model runs. (d) Barotropic streamfunction. (e) Overturning streamfunction. Positive values indicate clockwise circulation.

2.2 The biological carbon pump

In this paper, we adapt a representation of biological export from the surface ocean driven by an implicit (i.e. unresolved) biological community with a highly parameterized uptake of nutrients in the photic zone. The description of the basic scheme can be found in Ridgwell et al. (2007), although we use the specific configuration (including iron co-limitation) following Tagliabue et al. (2016). The governing equations are summarized below.

Biological productivity in the euphotic zone (taken to be the surface layer in the ocean model) is controlled by dissolved phosphate (${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$) and dissolved iron (Fe, which we will define later) availability, the fractional ice coverage of each grid cell (A), mean ambient light, and temperature. With the exception of the parameterization of the temperature term, the equation for photosynthetic nutrient uptake (M h⁻¹) follows Doney et al. (2006).

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathrm{\Lambda }& =F_I\cdot F_T\cdot min(F_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}},F_{{\mathrm{Fe}}^{\mathrm{3}+}})\cdot (\mathrm{1}-A)\\ & \cdot {\displaystyle \frac{min\left(\right[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}],{\mathrm{red}}_{\mathrm{P}/\mathrm{Fe}}\cdot [\mathrm{Fe}\left]\right)}{{\mathit{\tau }}_{\mathrm{bio}}}}\end{array}\end{array}$$

Rates of photosynthetic nutrient uptake are scaled to ambient dissolved nutrient concentrations (${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ and Fe³⁺) according to an optimal uptake timescale (τ_bio=1521.6 h; Meyer et al., 2016) and converting [Fe] to the equivalent dissolved phosphate concentration via the Fe:P Redfield ratio (${\mathrm{red}}_{\mathrm{P}/\mathrm{Fe}}$). The various limitation terms (all unitless) are

$$\begin{array}{}\text{(2)}& {\displaystyle }{F}_I={\displaystyle \frac{I}{I+{\mathit{\kappa }}_I}},\text{(3)}& {\displaystyle }{F}_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}}={\displaystyle \frac{\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]}{\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]+{\mathit{\kappa }}_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}}}},\text{(4)}& {\displaystyle }{F}_{{\mathrm{Fe}}^{\mathrm{3}+}}={\displaystyle \frac{\left[\mathrm{Fe}\right]}{\left[\mathrm{Fe}\right]+{\mathit{\kappa }}_{\mathrm{Fe}}}},\end{array}$$

where shortwave irradiance I is averaged over the entire mixed layer and is assumed to decay exponentially from the sea surface with a length scale of 20 m, as per Doney et al. (2006). The κ terms in each equation represent half-saturation constants for each limiting component (κ_I=40 W m⁻², ${\mathit{\kappa }}_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}}$=0.1 µM, ${\mathit{\kappa }}_{{\mathrm{Fe}}^{\mathrm{3}+}}=\mathrm{0.1}$ nM) and are as used in Tagliabue et al. (2016). The influence of temperature on biological export is parameterized as

$$\begin{array}{}\text{(5)}& F_T=k_{T\mathrm{0}}\cdot \exp \left({\displaystyle \frac{T}{k_{eT}}}\right),\end{array}$$

where k_T0 (0.59) is a scaling constant, k_eT (15.8 ^∘C) the e-folding temperature, and T is the in situ temperature (^∘C). The scaling constants give rise to an approximately factor of 2 change per temperature change of 10 ^∘C (Reinhard et al., 2020a).

A proportion (ν=0.66) of ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ taken up by biota is partitioned into dissolved organic phosphorus (DOP), while the remainder – as particulate organic phosphorus (POP) – is exported vertically out of the surface ocean. The value of ν has been assigned following the assumptions of the OCMIP-2 protocol (Najjar and Orr, 1999) (there is also an option for enacting temperature-dependent partitioning – and remineralization – of DOP, which presents an alternative to the fixed DOP–POP partitioning used here; Crichton et al., 2021). Because the biological configuration used here does not resolve explicit standing plankton biomass, the export flux of POP ($\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) is always equal to the rate of ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ uptake:

$$\begin{array}{}\text{(6)}& F_{z=h_{\mathrm{e}}}^{\mathrm{POP}}=\underset{h_{\mathrm{e}}}{\overset{\mathrm{0}}{\int }}\mathit{\rho }(\mathrm{1}-\mathit{\nu })\mathrm{\Lambda }\mathrm{d}z,\end{array}$$

where ρ is the density of seawater (kg m⁻³) and h_e the thickness of the euphotic zone (80.84 m). The particulate organic carbon (POC) export flux ($\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) is calculated using a fixed Redfield ratio as

$$\begin{array}{}\text{(7)}& F_{z=h_{\mathrm{e}}}^{\mathrm{POC}}=\mathrm{106}{F}_{z=h_{\mathrm{e}}}^{\mathrm{POP}}.\end{array}$$

After export from the surface grid cells, POC is remineralized instantaneously throughout the water column following a Martin-type curve (Martin et al., 1987), with a specified decay constant b (0.7) (Fig. 4a):

$$\begin{array}{}\text{(8)}& F_z^{\mathrm{POC}}=F_{z=h_{\mathrm{e}}}(z/h_{\mathrm{e}}{)}^{-b}.\end{array}$$

The b value of 0.7 is slightly higher than the global average of 0.6 estimated based on modern observations by Henson et al. (2011), but it leads to a better reconstruction of the distribution of CaCO₃ in deep-sea sediments (Jamie Wilson, personal communication).

3 Oceanic iron and sulfur cycling

The motivation behind this paper is to provide a tool that can aid understanding of the key interactions between the biological carbon pump, the oxygenation of the atmosphere and ocean, and the marine biogeochemical cycles of iron (Fe) and sulfur (S). In taking the first step towards this end, we construct a parsimonious model of ocean biogeochemical cycling based on a simplified speciation scheme for both Fe and S, and we only consider a relatively limited number of potential redox states. In this section we briefly describe the general conceptual cycle of Fe and S, as illustrated in Fig. 4b, and in the following subsections discuss the assumptions, reaction equations, and parameters for each cycle. The model equations that describe the biogeochemical reactions (summarized in Table 1) are given in Table 2, and the kinetic constants, their units, and default values can be found in Table 3.

Figure 4 (a) Example of the Martin curve for vertical organic matter fluxes in the ocean. (b) Conceptual description of the iron–sulfur cycle implemented in the model. For simplicity, only FeOOH (oxidized iron) and Fe²⁺ (reduced iron) are depicted, but the model includes different iron species in both the oxic and anoxic states. Interaction between the methane cycle and sulfur cycle is omitted for clarity. See text for details.

Table 1 List of biogeochemical reactions included in the reduced Fe–S scheme within the cGENIE model. A distinction is made between mineralization (“primary redox reactions”), re-oxidation of reduced products (“secondary redox reactions”), and precipitation reactions. The associated kinetic expressions are listed in Table 2. Note that we do not include any details for reactions involving the methane cycle here, though they are included in the simulations described here. For more information about the parameterization and reactions of the methane cycle in cGENIE, we refer the reader to Reinhard et al. (2020a).

^a DIR can only occur with solid-phase FeOOH, whereas SMI can occur with dissolved Fe³⁺ and solid-phase FeOOH. See text for more details. ^b The calculation of FeS_p occurs implicitly from the concentrations of Fe²⁺ and ΣH₂S as well as the solubility of FeS_aq; see Sect. 3.1.3 for more details.

Redox cycling in the ocean (and sediments) is driven by the mineralization of POC, which is produced in the photic zone by photosynthesis and subsequently sinks through the water column (Ridgwell et al., 2007). A Martin-type decay curve of organic matter flux with depth is prescribed in the model such that regardless of the relative availability of different electron acceptors, the fraction of organic matter that will be degraded (relative to the amount of organic matter produced in the photic zone) is depth-dependent. We avoid the alternative here – a fully kinetic set of equations wherein each electron acceptor is associated with a different rate of degradation – partly because of the additional set of poorly constrained (kinetic rate constant) parameters that would be required and partly because implementing such a scheme effectively requires knowledge about the composition of settling organic matter and how its relative reactivity changes with time (Ridgwell, 2011; LaRowe and Van Cappellen, 2011). Thus, the mineralization rate is dependent on the magnitude of the POC flux, which follows a Martin-type decay (Fig. 4a) that determines the mineralization rate (R_min) at each depth layer in the water column, and thus R_min (M h⁻¹) at depth z is calculated as

$$\begin{array}{}\text{(9)}& R_{\min }\left(z\right)={\displaystyle \frac{F_{\mathrm{POC}}(z-\mathrm{1})-F_{\mathrm{POC}}\left(z\right)}{\mathit{\rho }\mathrm{\Delta }z}},\end{array}$$

where Δz is the thickness of the grid cell at depth z (in metres). Note that mineralization of particulate organic matter only occurs below the surface layers; z>1 in Fig. 2. The mineralization of POC is coupled to the reduction of a terminal electron acceptor (TEA). These TEAs are used according to decreasing energy yield (Froelich et al., 1979), and relative consumption rates are scaled with TEA concentration and the local abundance of inhibitory substances (i.e. a more energy-yielding TEA). Our mineralization scheme includes aerobic respiration (AR, R1 in Table 1), dissimilatory iron reduction (DIR, R2 in Table 1), dissimilatory sulfate reduction (DSR, R3 in Table 1), and methanogenesis (MG, R4 in Table 1). Of these, AR and DSR existed in the original cGENIE biogeochemical framework (Ridgwell et al., 2007), while methanogenesis has been more recently added (Reinhard et al., 2020a). In the current paper, we ignore nitrate reduction (Monteiro et al., 2012).

The rate of TEA consumption is represented by a Michaelis–Menten type of relationship with respect to TEA concentration that allows for a non-linear closure of the system.

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{f}_{\mathrm{AR}}& ={\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{K_{\mathrm{0},{\mathrm{O}}_{\mathrm{2}}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\\ f_{\mathrm{DIR}}& ={\displaystyle \frac{\left[\mathrm{FeOOH}\right]}{K_{\mathrm{0},\mathrm{FeOOH}}+\left[\mathrm{FeOOH}\right]}}{\displaystyle \frac{K_{\mathrm{i},{\mathrm{O}}_{\mathrm{2}}}}{K_{\mathrm{i},{\mathrm{O}}_{\mathrm{2}}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\\ f_{\mathrm{DSR}}& ={\displaystyle \frac{\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]}{K_{\mathrm{0},{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}}+\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]}}\\ & \cdot {\displaystyle \frac{K_{\mathrm{i},\mathrm{FeOOH}}}{K_{\mathrm{i},\mathrm{FeOOH}}+\left[\mathrm{FeOOH}\right]}}{\displaystyle \frac{K_{\mathrm{i},{\mathrm{O}}_{\mathrm{2}}}}{K_{\mathrm{i},{\mathrm{O}}_{\mathrm{2}}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\\ f_{\mathrm{MG}}& ={\displaystyle \frac{K_{\mathrm{i},{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}}}{K_{\mathrm{i},{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}}+\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]}}\\ & \cdot {\displaystyle \frac{K_{\mathrm{i},\mathrm{FeOOH}}}{K_{\mathrm{i},\mathrm{FeOOH}}+\left[\mathrm{FeOOH}\right]}}{\displaystyle \frac{K_{\mathrm{i},{\mathrm{O}}_{\mathrm{2}}}}{K_{\mathrm{i},{\mathrm{O}}_{\mathrm{2}}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\end{array}\end{array}$$

Here, K₀ represents the half-saturation constants for the four primary redox reactions and K_i represents the inhibition constants that act on the less energetic redox reaction (Table 3). The consequence of this scheme is that in the oxic zone of the water column, aerobic respiration is responsible for nearly all of the POC mineralization. When oxygen starts to become depleted, DIR and subsequently DSR become the dominant mineralization pathways. When both iron oxides and sulfate are exhausted, MG represents the final mineralization pathway. While the occurrence and parameterization of DSR and MG in the water column are relatively straightforward, the possibility of DIR in the water column is more complex and is discussed in more detail in Sect. 3.1.2.

Table 2 List of kinetic rate expressions for the reactions included in the cGENIE model. All expressions are based on standard kinetic formulations in biogeochemical models. The values of the kinetic constants are listed in Table 3.

^a Two different kinetic constants are used for the reactions between Fe³⁺ and sulfide (k_SMI,d) and between FeOOH and sulfide (k_SMI,s); see text for details. ^b The formation of iron oxide minerals is parameterized according to Parekh et al. (2004) and is explained in more detail in Sect. 3.1.1. References: [1] Soetaert et al. (1996), [2] van de Velde and Meysman (2016), [3] Millero et al. (1987a), [4] Dale et al. (2015), [5] Ridgwell et al. (2007), [6] Zhang and Millero (1993), [7] Poulton et al. (2004), [8] Parekh et al. (2004), [9] Tagliabue et al. (2016), [10] Rickard (1997), [11] Rickard (2006), [12] Jiang and Tosca (2019), [13] Tosca et al. (2015), [14] Rasmussen et al. (2015)

As a consequence of organic matter remineralization in the model, DIR produces ferrous iron (Fe²⁺), which is re-oxidized when it comes into contact with oxygen (e.g. via upwelling of the reduced compound or downwelling of oxygen) via ferrous iron oxidation (FIO, R5 in Table 1) (Millero et al., 1987a). Similarly, any reduced sulfide (because we do not consider sulfide speciation explicitly, reduced sulfide is defined as $\mathrm{\Sigma }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}={\mathrm{H}}_{\mathrm{2}}\mathrm{S}+{\mathrm{HS}}^{-}$) that comes into contact with oxygen is re-oxidized via canonical sulfur oxidation (CSO, R6 in Table 1) (Millero et al., 1987b). Additionally, H₂S is oxidized by reaction with oxidized iron via sulfide-mediated iron reduction (SMI, R7 in Table 1) (Fig. 4b; Canfield et al., 1992; Poulton et al., 2004; Mikucki et al., 2009). The oxidized form of iron, ferric iron (Fe³⁺), will precipitate out as iron oxide (FeOOH) minerals (IrP, R8 in Table 1).

When Fe²⁺ and H₂S are simultaneously present in the water column they form dissolved FeS (FeS_aq). Once FeS_aq surpasses a solubility threshold of ∼2 µM (Rickard, 2006), particulate FeS_p can form, which further reacts with H₂S to form the mineral pyrite (FeS₂) via pyrite precipitation (PyP, R9 in Table 1) (Rickard, 1997). Alternatively, when Fe²⁺ accumulates past its saturation state with bicarbonate, siderite (FeCO₃) forms via siderite precipitation (SiP, R10 in Table 1) (Jimenez-Lopez et al., 2004; Jiang and Tosca, 2019). Finally, greenalite (Fe₃Si₂O₅(OH)₄) precipitates when Fe²⁺ and dissolved silica (SiO₂) are saturated with respect to greenalite precipitation (GrP, R11 in Table 1) (Tosca et al., 2015; Rasmussen et al., 2015). Pyrite, siderite, and greenalite are subsequently buried in the sediment, together with any solid FeOOH that has not reacted (the half-life of iron oxides in euxinic waters is of the order of tens to hundreds of days, which is comparable to the residence time of a particle in the ocean; Poulton et al., 2004) (Fig. 4b). These four solid-iron phases are the main burial phases for reactive iron and form the basis of the Fe-speciation proxy used to reconstruct local redox conditions in past oceans (Poulton and Canfield, 2011). It should be noted, however, that some of these phases can undergo transformations to other phases in the sediment after deposition (such as greenalite to magnetite). In our current model set-up, no sedimentary processes are included, but future developments will address the sedimentary part of the Fe and S cycle.

Many of the iron and sulfur reactions can go to completion much faster than the biogeochemical time step ($\mathrm{1}/\mathrm{24}$ year in the example conceptual configuration) under typical modern or paleo-geochemical conditions. This can create negative tracer concentrations because transport by ocean circulation acts concurrently on the tracer field in our numerical scheme. We hence place a limit on reactant consumption for any reaction in which a reactant would be completely consumed in a single time step. Specifically, we prescribe a maximum timescale for all geochemical reactions, for which in this paper we chose a value of 45 d.

3.1 The iron cycle

The cGENIE model already included the representation of a simplified iron cycle designed to account for iron limitation of biological productivity at the ocean surface (Tagliabue et al., 2016). However, because its initial use was to model the present-day oceanic iron cycle, it does not contain an anoxic iron cycle (as 95 % of today's ocean volume is oxic; Diaz and Rosenberg, 2008) or any coupling between the iron and sulfur cycles (e.g. via the precipitation of FeS₂). The absence of an anoxic iron cycle limits the suitability of cGENIE for simulating low-oxygen worlds in which the ocean interior is pervasively anoxic and iron and sulfur cycling would have dominated ocean biogeochemistry (Poulton and Canfield, 2011; Raiswell and Canfield, 2012). In addition to summarizing the existing oxic cycle, we expand the model to include key processes operating in anoxic environments.

3.1.1 The oxic iron cycle

In oxygenated waters, iron is predominately present in its oxidized form (Fe³⁺). The oxic iron cycle in cGENIE follows Parekh et al. (2004) and has been recently updated with a revised set of parameters calibrated based on the present-day iron cycle (Tagliabue et al., 2016). Briefly, in the previous scheme total dissolved Fe³⁺ consists of free ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{3}+}$ and ligand-bound ${\mathrm{Fe}}_{\mathrm{ligand}}^{\mathrm{3}+}$, with only the “free” fraction being subject to scavenging on particulate organic carbon (POC) particles (whereas both free and ligand-bound dissolved iron phases are assumed to be bioavailable and fuel primary productivity; see Sect. 2.2). When scavenged on POC, FeOOH sinks as a “marine snow” particle through the water column. This scavenging mechanism is commonly used in reactive-transport models of ferruginous lakes or the modern ocean (see e.g. Taillefert and Gaillard, 2002; Tagliabue et al., 2016). Marine snow formation is what likely happens today, whereby relatively high production of organic matter drives the transport of chemically heterogeneous aggregates (containing iron oxides or barite; Dehairs et al., 1990; Balzano et al., 2009). The scavenging rate of iron is a function of the concentration of ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{3}+}$ together with the magnitude of the POC flux from the grid cell, i.e.

$$\begin{array}{}\text{(11)}& R_{\mathrm{IrP}}=k_{\mathrm{scav}}\left[{\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{3}+}\right]F_{\mathrm{POC}}\left(z\right),\end{array}$$

where k_scav ($\mathrm{1.43}\times {\mathrm{10}}^{-\mathrm{6}}$ mol⁻¹ m²) is a scavenging constant calibrated to the modern-day distribution of Fe (Tagliabue et al., 2016). We assume that the complexation reaction is always in equilibrium, which allows for the calculation of the amount of ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{3}+}$ at each time step by the conservation equations.

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\left[{\mathrm{Fe}}^{\mathrm{3}+}\right]& =\left[{\mathrm{Fe}}_{\mathrm{ligand}}^{\mathrm{3}+}\right]+\left[{\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{3}+}\right]\\ \left[L_{\mathrm{total}}\right]& =\left[{\mathrm{Fe}}_{\mathrm{ligand}}^{\mathrm{3}+}\right]+\left[L_{\mathrm{free}}\right]\\ K_{\mathrm{sp}}^{\mathrm{FeL}}& ={\displaystyle \frac{\left[{\mathrm{Fe}}_{\mathrm{ligand}}^{\mathrm{3}+}\right]}{\left[{\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{3}+}\right]\left[L_{\mathrm{free}}\right]}}\end{array}\end{array}$$

Here, $K_{\mathrm{sp}}^{\mathrm{FeL}}$ is the stability constant of the $L-{\mathrm{Fe}}^{\mathrm{3}+}$ complex (1.0×10¹¹ M⁻¹; Table 3), L_free is ligand unassociated with Fe³⁺, and L_total is the total amount of ligand.

Oxidized iron is highly insoluble in seawater and will rapidly form particulate oxides (FeOOH), driving down equilibrium dissolved Fe³⁺ abundance to picomolar values without stabilization by ligands (Liu and Millero, 2002). Thus, ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{3}+}$ will form colloidal or nanoparticulate iron oxides – FeOOH – which can then adsorb on other particles or aggregate to bigger particles and sink through the water column (Raiswell and Canfield, 2012). In cGENIE, the pool of ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{3}+}$ can be thought of as a mix of purely dissolved Fe³⁺ and a range of colloidal and nanoparticulate Fe³⁺ phases that do not settle efficiently through the water column. The ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{3}+}$ pool can then further react or scavenge onto particles that allow it to more effectively settle through the water column.

In our new formulation of the Fe cycle in cGENIE, ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{3}+}$ (that is not stabilized by ligands) is still the phase that is scavenged by settling POC (using a fixed scavenging rate; Eq. 11), but it is explicitly assumed to be in the form of FeOOH and can be used for DIR when oxygen becomes depleted (see Sect. 3.1.2). We also provide the option in the model for ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{3}+}$ to precipitate as a pure FeOOH phase, without being associated with POC. In this case, dissimilatory iron reduction does not occur in the water column since both POC and FeOOH are not associated within the same particle. This situation would be more appropriate in the case of high amounts of iron oxidation and lower production of organic matter (which was potentially the case in the Archean; Thompson et al., 2019). This configuration is by default implemented using a numerical “cut-off”, at which all ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{3}+}$ that passes the solubility threshold of 0.5 nM (Liu and Millero, 2002) is considered particulate and sinks through the water column.

To close the oceanic iron cycle, we have introduced a new dissolved iron species, Fe²⁺. Any reduced Fe²⁺ that is mixed into the oxic zone is being rapidly oxidize to Fe³⁺ via ferrous iron oxidation (FIO):

$$\begin{array}{}\text{(13)}& {\mathrm{Fe}}^{\mathrm{2}+}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\mathrm{O}}_{\mathrm{2}}+{\mathrm{H}}^{+}\to {\mathrm{Fe}}^{\mathrm{3}+}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{H}}_{\mathrm{2}}\mathrm{O},\end{array}$$

where the rate equation can be expressed as

$$\begin{array}{}\text{(14)}& R_{\mathrm{FIO}}=k_{\mathrm{FIO}}\left[{\mathrm{Fe}}^{\mathrm{2}+}\right]\left[{\mathrm{O}}_{\mathrm{2}}\right]\end{array}$$

and where k_FIO is a reaction rate constant ($k_{\mathrm{FIO}}=\mathrm{0.115}\times {\mathrm{10}}^{\mathrm{6}}$ ${\mathrm{M}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$) (Table 2; Millero et al., 1987a). In the calculation of the nutrient limitation of biological production at the surface, Fe²⁺ is also assumed to be bioavailable, meaning that iron limitation is calculated from the total dissolved iron pool, equal to ${\mathrm{Fe}}^{\mathrm{2}+}+{\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{3}+}+{\mathrm{Fe}}_{\mathrm{ligand}}^{\mathrm{3}+}$.

In summary, we have extended the basic representation of Parekh et al. (2004), which considered only a generic free “dissolved” iron (that could be scavenged) and a ligand-bound iron phase, to now distinguish between the oxidation states of iron (and the tracers Fe²⁺ and Fe³⁺), with Fe³⁺ being split (as previously) into a free (and scavengeable) and ligand-bound phase.

3.1.2 The anoxic iron cycle

In an anoxic water column, oxidized forms of iron can be reduced to Fe²⁺ either via dissimilatory iron reduction or via sulfur-mediated iron reduction. We describe the two processes individually as follows.

Dissimilatory iron reduction in the water column

The majority of oxidized iron exists in particulate form due to the low solubility of Fe³⁺, which poses a challenge for micro-organisms performing dissimilatory iron reduction in the water column. Iron reducers need physical contact with the iron mineral to reduce iron (Gorby, 2006), which is difficult when both organic matter and FeOOH are in particulate form and are dispersed separately throughout the aqueous medium. Even in a sediment column in which all particles are packed closely together, DIR generally requires sediment homogenization by burrowing fauna to become volumetrically important (Thamdrup, 2000; van de Velde and Meysman, 2016). Indeed, in a sulfate-rich, ancient marine brine, Fe³⁺ was found to be the terminal electron acceptor, but sulfur ultimately acted as a redox shuttle between organic matter and iron oxides (Mikucki et al., 2009). Studies of anoxic Lake Pavin suggest that most of the iron reduction occurs in the sediment rather than the water column (Michard et al., 1994; Cosmidis et al., 2014), although the same studies indicated that manganese oxide reduction (which has a very similar mechanism and inhibition concentration as DIR; Lovley, 1991; Thamdrup, 2000) does occur in the water column. However, there is direct evidence for iron reduction occurring coupled to organic matter mineralization in the water column in a range of other anoxic lakes, such as Lake Sammamish (Washington, USA; Balistrieri et al., 1992), Lake Cadagno (Switzerland; Berg et al., 2016), Lake Matano (Indonesia; Crowe et al., 2008), and Paul Lake (Michigan, USA; Taillefert and Gaillard, 2002). These studies suggested that DIR was coupled to either the oxidation of dissolved organic matter (Crowe et al., 2008) or reduction of iron in aggregates with organic matter. The latter has been experimentally shown to occur (Balzano et al., 2009). We consider this to be strong evidence that DIR can be coupled to POC oxidation in a water column, especially in the ocean where the residence time of particles is considerably longer than in a much shallower lake environment.

To model DIR in the water column, we use the mathematical expression of DIR in marine sediments, where FeOOH is a common electron acceptor for organic matter mineralization (Thamdrup, 2000). Dissimilatory iron reduction has an overall reaction stoichiometry of

$$\begin{array}{}\text{(15)}& {\mathrm{CH}}_{\mathrm{2}}\mathrm{O}+\mathrm{4}\mathrm{FeOOH}+\mathrm{7}{\mathrm{H}}^{+}\to {\mathrm{HCO}}_{\mathrm{3}}^{-}+\mathrm{4}{\mathrm{Fe}}^{\mathrm{2}+}+\mathrm{6}{\mathrm{H}}_{\mathrm{2}}\mathrm{O},\end{array}$$

where it is implicitly assumed that every iron oxide particle consists solely of Fe³⁺ rather than a mixture of redox states, as is sometimes encountered at the interface of ferruginous lakes (Zegeye et al., 2012). This assumption greatly simplifies the reaction scheme and parameter set. Dissimilatory iron reduction generally becomes limited when concentrations of FeOOH drop below $\sim \mathrm{30}\times {\mathrm{10}}^{-\mathrm{3}}$ M (Van Cappellen and Wang, 1996; Thamdrup, 2000). This limitation is expressed as $\frac{\left[\mathrm{FeOOH}\right]}{K_{\mathrm{0},\mathrm{FeOOH}}+\left[\mathrm{FeOOH}\right]}$ in Eq. (10) and follows the conventional limitation–inhibition scheme (Soetaert et al., 1996), whereby the parameter K_0,FeOOH expresses the concentration at which DIR occurs at half of its maximum rate (the maximum rate is set by R_min; see Eq. 10 and Table 2). In marine sediments, DIR generally occurs before sulfate reduction since it is a more energy-yielding electron acceptor (Thamdrup, 2000), although the wide range of reactivities for different iron oxide minerals often leads to an overlap of DIR and DSR zones (Postma and Jakobsen, 1996). This limitation is expressed as $\frac{\left[K_{\mathrm{i},\mathrm{FeOOH}}\right]}{K_{\mathrm{i},\mathrm{FeOOH}}+\left[\mathrm{FeOOH}\right]}$ in Eq. (10), where K_i,FeOOH expresses the concentration above which FeOOH inhibits other mineralization pathways (DSR and MG). In early diagenetic models, this concentration is generally assumed to be identical to the limitation parameter K_0,FeOOH (Van Cappellen and Wang, 1996; Soetaert et al., 1996; van de Velde and Meysman, 2016). As a baseline value, we assume that K_0,FeOOH and K_i,FeOOH both equal 10⁻³ M (comparable to the inhibition concentration in marine sediments), but as explained above, these parameter values are subject to a high degree of uncertainty. Therefore, we include parameters K_0,FeOOH and K_i,FeOOH in our model sensitivity testing (Sect. 4.3).

Sulfur-mediated iron reduction

Oxidized iron in the ocean can also be reduced via sulfur-mediated iron reduction (SMI), which follows the stoichiometry (Poulton et al., 2004)

$$\begin{array}{}\text{(16)}& \mathrm{\Sigma }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}+\mathrm{2}{\mathrm{Fe}}^{\mathrm{3}+}\to {\mathrm{S}}^{\mathrm{0}}+\mathrm{2}{\mathrm{Fe}}^{\mathrm{2}+}+\mathrm{2}{\mathrm{H}}^{+}.\end{array}$$

We are not explicitly modelling elemental sulfur (S⁰) but assume that it becomes quantitatively disproportionated into H₂S and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (Finster et al., 1998),

$$\begin{array}{}\text{(17)}& {\mathrm{S}}^{\mathrm{0}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\to {\displaystyle \frac{\mathrm{3}}{\mathrm{4}}}\mathrm{\Sigma }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{H}}^{+},\end{array}$$

which then leads to the overall stoichiometry

$$\begin{array}{}\text{(18)}& \mathrm{\Sigma }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}+\mathrm{8}{\mathrm{Fe}}^{\mathrm{3}+}+\mathrm{4}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\to {\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}+\mathrm{8}{\mathrm{Fe}}^{\mathrm{2}+}+\mathrm{10}{\mathrm{H}}^{+}.\end{array}$$

Note that [Fe³⁺] in Eq. (19) can represent dissolved Fe³⁺ or solid FeOOH. The assumption of quantitative disproportionation implies that pyrite precipitation is not closely coupled to the reaction between ΣH₂S and Fe³⁺ but only occurs via precipitation of FeS_p with ΣH₂S (see Sect. 3.1.3). Laboratory experiments have shown that this assumption is valid for aquatic systems in which Fe³⁺ is not in excess with respect to ΣH₂S (Wan et al., 2017), which is the case for most modern marine systems. We contend that this is also a valid assumption for water-column chemistry for most of Earth's history, as rapid settling of oxidized particulate FeOOH through the water column would prevent high concentrations of Fe³⁺ (in the water column). To achieve an excess of Fe³⁺ over ΣH₂S, the ΣH₂S concentrations would have to be even lower, leading to very negligible rates of iron reduction. However, this reaction pathway would likely become important in the sediment of a low-sulfate ocean (i.e. periods of Archean time).

The kinetic rate for Reaction (18) can then be expressed as (Poulton et al., 2004)

$$\begin{array}{}\text{(19)}& R_{\mathrm{SMI}}=k_{\mathrm{SMI},\mathrm{d}/\mathrm{s}}\left[\mathrm{\Sigma }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}{]}^{\mathrm{0.5}}\right[{\mathrm{Fe}}^{\mathrm{3}+}].\end{array}$$

We assume that the dissolved form is highly reactive with sulfide (a reaction time of ∼5 min, which is comparable to the reactivity of freshly precipitated hydrous ferric oxide) and has a reaction rate constant of $k_{\mathrm{SMI},\mathrm{d}}=\mathrm{2.64}\times {\mathrm{10}}^{\mathrm{2}}$ ${\mathrm{M}}^{-\mathrm{0.5}}{\mathrm{h}}^{-\mathrm{1}}$ (Poulton et al., 2004). The solid form of FeOOH can represent a number of different oxidized iron minerals which are reactive towards sulfide on timescales ranging from hours to hundreds of days (Poulton et al., 2004). For our baseline simulations, we assume that all FeOOH precipitates as lepidocrocite, which has a reactivity constant of $k_{\mathrm{SMI},\mathrm{s}}=\mathrm{1.98}\times {\mathrm{10}}^{\mathrm{0}}$ ${\mathrm{M}}^{-\mathrm{0.5}}{\mathrm{h}}^{-\mathrm{1}}$ (Poulton et al., 2004). Lepidocrocite is less crystalline than the other (non-hydrous) iron oxides and precipitates when the rate of Fe³⁺ supply is low relative to the rate of precipitation, as is generally the case in natural systems (Crosby et al., 1983). We discuss the model sensitivity to choices of k_SMI,s in Sect. 4.3. The reduction of Fe³⁺ produces Fe²⁺, which can either be re-oxidized when it comes in contact with O₂ (see Sect. 3.1.1) or can form reduced minerals.

Both reduction and oxidation reactions of dissolved iron, even with dissolved O₂ and ΣH₂S in nanomolar (nM) concentrations, can proceed very rapidly and are hence subject to the geochemical reaction rate limitation described earlier. Furthermore, because these two reactions form a coupled oxidation–reduction pair, we limit the fastest reaction according to a 45 d timescale but limit the slower one in proportion to their relative unmodified reactions rates. We thereby simulate an equilibrium partitioning between Fe²⁺ and Fe³⁺ according to ambient dissolved oxygen and sulfide concentrations.

3.1.3 Reduced iron mineral formation

Reduced Fe²⁺ can form complexes with a number of inorganic ligands that are common in seawater, including Cl⁻, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, and bicarbonate (Millero et al., 1995). Under anoxic conditions, however, the most important inorganic ligand is free sulfide (Rickard, 2006), which is produced by sulfate reduction. Together, dissolved ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{2}+}$ and the aqueous iron-sulfide complex (FeS_aq) make up ∼100 % of the total dissolved iron pool in sulfidic–anoxic seawater (Fig. 5a). The thermodynamic equilibrium can be calculated as

$$\begin{array}{}\text{(20)}& K_{\mathrm{sp}}^{{\mathrm{FeS}}_{\mathrm{aq}}}={\displaystyle \frac{\left[{\mathrm{Fe}}^{\mathrm{2}+}\right]\left[\mathrm{\Sigma }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}\right]}{\left[{\mathrm{FeS}}_{\mathrm{aq}}\right]}}={\mathrm{10}}^{-\mathrm{5.08}}\mathrm{M},\end{array}$$

which was obtained from the visualMINTEQ database (Gustafsson, 2019) and is based on stability constants calculated by Luther et al. (1996).

Since dissolved ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{2}+}$ and the FeS_aq complex are the dominant dissolved forms of reduced ferrous iron in natural waters (in anoxic seawater devoid of sulfide, dissolved ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{2}+}$ still represents ∼80 % of the dissolved Fe²⁺ pool; Fig. 5b), we choose to only consider those two species. We do not explicitly model the FeS_aq complex but calculate the thermodynamic equilibrium before each reaction proceeds, implicitly assuming that it is reached much faster than any of the kinetic reactions in our reaction set (as is commonly done for the complexation of Fe³⁺; see Sect. 3.1.1). This is important, since FeS_aq forms particulate FeS_p, which is the precursor of pyrite (FeS₂) (Rickard, 1997, 2006). Furthermore, siderite or greenalite can only precipitate from ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{2}+}$ (Tucker, 2001; Tosca et al., 2015), and their precipitation rates are consequently dependent on the dissolved ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{2}+}$ concentration.

Figure 5 Relative importance of different dissolved Fe complexes, obtained by visualMINTEQ (Gustafsson, 2019). Concentrations of other constituents are chosen to represent average seawater, with the dissolved organic carbon (DOC) concentration from Sharp et al. (1995).

Before each precipitation reaction, we calculate the concentration of FeS_aq, ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{2}+}$, and H₂S_free using Eq. (20). We then compare the FeS_aq concentration to the solubility threshold (10^−5.7 M; Rickard, 2006). If FeS_aq surpasses the solubility value, the fraction above the solubility concentration precipitates instantaneously as particulate iron monosulfide (FeS_p) (Suits and Wilkin, 1998). Note that here we consider FeS_p solubility to be independent of pH, which should be valid for seawater pH above 7 (Rickard, 2006). Pyrite can then form locally via reaction between FeS_p and free ΣH₂S, with the production of free hydrogen gas (as S²⁻ in ΣH₂S and FeS_p has to be oxidized to S¹⁻ in FeS₂; Rickard, 1997):

$$\begin{array}{}\text{(21)}& {\mathrm{FeS}}_p+\mathrm{\Sigma }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}\to {\mathrm{FeS}}_{\mathrm{2}}+{\mathrm{H}}_{\mathrm{2}}.\end{array}$$

Since hydrogen gas is highly reactive, it is almost instantaneously oxidized with the reduction of an electron acceptor like O₂, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, or CO₂. Given that pyrite precipitates under sulfidic conditions, oxygen is absent, and sulfate is the most likely electron acceptor.

$$\begin{array}{}\text{(22)}& {\mathrm{H}}_{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{H}}^{+}\to {\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\mathrm{\Sigma }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\end{array}$$

To preserve the redox balance while avoiding the need to model an additional state variable (H₂), we can combine Reactions (21) and (22) to

$$\begin{array}{}\text{(23)}& {\mathrm{FeS}}_p+{\displaystyle \frac{\mathrm{3}}{\mathrm{4}}}\mathrm{\Sigma }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{H}}^{+}\to {\mathrm{FeS}}_{\mathrm{2}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}.\end{array}$$

By using Eq. (23) we implicitly assume that the hydrogen gas produced during pyrite precipitation immediately reacts with ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$. Since pyrite formation is not an equilibrium reaction (pyrite minerals are stable in seawater), the reaction of pyrite is described as a kinetic reaction with a second-order dependency on [FeS_p] and [ΣH₂S], with $k_{\mathrm{PyP}}=\mathrm{0.3708}\times {\mathrm{10}}^{\mathrm{0}}$ ${\mathrm{M}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$ (Rickard, 1997). The rate equation for FeS₂ precipitation can then be written as

$$\begin{array}{}\text{(24)}& R_{\mathrm{PyP}}=k_{\mathrm{PyP}}\left[{\mathrm{FeS}}_p\right]\left[\mathrm{\Sigma }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}\right].\end{array}$$

Any Fe²⁺ that is not complexed with sulfide (${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{2}+}$) can form siderite (FeCO₃), an iron–carbonate mineral (Tucker, 2001).

$$\begin{array}{}\text{(25)}& {\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{2}+}+{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\to {\mathrm{FeCO}}_{\mathrm{3}}\end{array}$$

Recent experiments by Jiang and Tosca (2019) have shown that the precipitation of iron–carbonate phases is controlled by the formation of amorphous Fe carbonate (AFC), which has the stoichiometric formula ${\mathrm{FeCO}}_{\mathrm{3}}(\mathrm{OH}{)}_{\mathrm{1}/\mathrm{2}}$. Since the thermodynamic data required to calculate the mineral solubility product of AFC are currently lacking, Jiang and Tosca (2019) have determined the rate of AFC precipitation based on the ion activity product (IAP):

$$\begin{array}{}\text{(26)}& \begin{array}{rl}{\mathrm{IAP}}_{\mathrm{AFC}}& ={\mathit{\gamma }}_{{\mathrm{Fe}}^{\mathrm{2}+}}\left[{\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{2}+}\right]{\mathit{\gamma }}_{{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]\\ & \cdot \left({\mathit{\gamma }}_{{\mathrm{OH}}^{-}}\right[{\mathrm{OH}}^{-}]{)}^{\mathrm{0.5}},\end{array}\end{array}$$

where γ represents the activity coefficient of a given ion in seawater (Table 3). The rate-limiting step in the precipitation reaction is the spontaneous nucleation of AFC, and hence the empirically derived rate of AFC precipitation follows an exponential function that describes the nucleation rate from water (Steefel and Van Cappellen, 1990; Jiang and Tosca, 2019):

$$\begin{array}{}\text{(27)}& R_{\mathrm{SiP}}=k_{\mathrm{AFC}}{e}^{b_{\mathrm{AFC}}{\log }_{\mathrm{10}}{\mathrm{IAP}}_{\mathrm{AFC}}},\end{array}$$

where k_AFC is $\mathrm{1.963}\times {\mathrm{10}}^{-\mathrm{14}}$ M h⁻¹ and b_AFC equals 9.042 (Table 3).

Another reduced iron mineral potentially important for anoxic, Fe-rich environments is the iron–silicate mineral greenalite (Fe₃Si₂O₅(OH)₄) (Tosca et al., 2015).

$$\begin{array}{}\text{(28)}& \mathrm{3}{\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{2}+}+\mathrm{2}{\mathrm{SiO}}_{\mathrm{2}}\left(\mathrm{aq}\right)+\mathrm{5}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\to {\mathrm{Fe}}_{\mathrm{3}}{\mathrm{Si}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}(\mathrm{OH}{)}_{\mathrm{4}}+\mathrm{6}{\mathrm{H}}^{+}\end{array}$$

The precipitation rate of greenalite is dependent on its degree of supersaturation (i.e. its saturation index), which can be calculated based on its IAP:

$$\begin{array}{}\text{(29)}& {\text{IAP}}_{\mathrm{greenalite}}={\displaystyle \frac{\left({\mathit{\gamma }}_{{\mathrm{Fe}}^{\mathrm{2}+}}\right[{\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{2}+}\left]{)}^{\mathrm{3}}\right({\mathit{\gamma }}_{{\mathrm{SiO}}_{\mathrm{2}}}\left[{\mathrm{SiO}}_{\mathrm{2}}\right]{)}^{\mathrm{2}}}{\left({\mathit{\gamma }}_{{\mathrm{H}}^{+}}\right[{\mathrm{H}}^{+}]{)}^{\mathrm{6}}}},\end{array}$$

where γ represents the activity coefficient of an ion in seawater (Table 3) and $K_{\mathrm{sp}}^{\mathrm{greenalite}}$ its solubility product (Rasmussen et al., 2015):

$$\begin{array}{}\text{(30)}& {\text{SI}}_{\mathrm{greenalite}}={\log }_{\mathrm{10}}\left({\displaystyle \frac{\text{IAP}}{K_{\mathrm{sp}}^{\mathrm{greenalite}}}}\right),\end{array}$$

where $K_{\mathrm{sp}}^{\mathrm{greenalite}}$ equals 3.98×10²⁷ M⁻¹ (Tosca et al., 2015). The rate equation follows an exponential function that describes the nucleation rate from water (as in siderite, see above; Rasmussen et al., 2015; Jiang and Tosca, 2019):

$$\begin{array}{}\text{(31)}& R_{\mathrm{GrP}}=k_{\mathrm{greenalite}}{e}^{b_{\mathrm{greenalite}}{\log }_{\mathrm{10}}{\text{SI}}_{\mathrm{greenalite}}},\end{array}$$

where k_greenalite is $\mathrm{6.996}\times {\mathrm{10}}^{-\mathrm{13}}$ M h⁻¹ and b_greenalite equals 1.856 (Table 3).

3.2 The sulfur cycle

The oxidized form of dissolved sulfur in the ocean is sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$). Under anaerobic conditions, sulfate is used as a terminal electron acceptor during organic matter remineralization:

$$\begin{array}{}\text{(32)}& {\mathrm{CH}}_{\mathrm{2}}\mathrm{O}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}+{\mathrm{H}}^{+}\to {\mathrm{CO}}_{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Sigma }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}.\end{array}$$

The produced sulfide can then be re-oxidized when it comes into contact with oxygen via canonical sulfide oxidation (CSO).

$$\begin{array}{}\text{(33)}& \mathrm{\Sigma }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}+\mathrm{2}{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}+\mathrm{2}{\mathrm{H}}^{+}\end{array}$$

The rate of CSO in cGENIE is dependent on the concentration of the electron donor (ΣH₂S) and acceptor (O₂) as well as a second-order rate constant $k_{\mathrm{CSO}}=\mathrm{0.625}\times {\mathrm{10}}^{\mathrm{6}}$ ${\mathrm{M}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ (Zhang and Millero, 1993):

$$\begin{array}{}\text{(34)}& R_{\mathrm{CSO}}=k_{\mathrm{CSO}}\left[\mathrm{\Sigma }{\mathrm{H}}_{\mathrm{2}}\mathrm{S}\right][{\mathrm{O}}_{\mathrm{2}}{]}^{\mathrm{2}}.\end{array}$$

Alternatively, ΣH₂S can be re-oxidized with Fe³⁺ (either Fe³⁺ or solid FeOOH) during sulfide-mediated iron reduction (as described in Reaction 18).

In cGENIE, the eventual sink for sulfur is the precipitation of pyrite (Reaction 21) (ignoring sulfide reacting with organic matter – see Hülse et al., 2019). We do not currently include precipitation of gypsum (CaSO₄) in our model description (Fig. 4b). Gypsum is an evaporite mineral that precipitates during regional and episodic events of supersaturation and was likely a less important sulfur sink on a globally integrated basis during Precambrian time or during any other period in which ocean [${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$] was relatively low (Grotzinger and Kasting, 1993; Crowe et al., 2014; Fakhraee et al., 2019). Indeed, there is still some debate as to the time-integrated impact of sulfate evaporites on the steady-state global sulfur cycle even during more recent periods of Earth's history (Halevy et al., 2012; Canfield, 2013). However, due to its episodic nature, gypsum could play an important role as a sulfate source during transient events (for example during events of enhanced weathering of a gypsum-rich source; Shields et al., 2019). Planned future developments to cGENIE will incorporate an explicit gypsum cycle, which would allow us to use the cGENIE.muffin model to investigate transient events of enhanced sulfate delivery to the ocean (see e.g. Shields et al., 2019).

Table 3 List of kinetic constants for the reactions included in the cGENIE model. Note that the units are not identical to the units used in the user configuration files for the cGENIE model. A table with the parameter units converted to cGEnIE units is included in Appendix A.

[1] Ridgwell et al. (2007), [2] Thamdrup (2000), [3] Olson et al. (2016), [4] Zhang and Millero (1993), [5] Millero et al. (1987a), [6] Poulton et al. (2004), [7] Ridgwell and Death (2021), [8] Luther et al. (1996), [9] Rickard (1997), [10] Jiang and Tosca (2019), [11] Tosca et al. (2015), [12] Rasmussen et al. (2015), [13] Marion et al. (2011), [14] (following the Davies equation), [15] Johnson (1982).

3.3 Isotope geochemistry

A particularly important application of having a representation of an anoxic Fe–S cycle in cGENIE is exploring ocean redox landscapes during Precambrian time when ocean biogeochemical cycling was dominated by iron and sulfur (Raiswell and Canfield, 2012). As noted above, one way of comparing our model output to available data is the explicit simulation of the burial phases of iron, which allows comparison to the often-used Fe proxy (Poulton and Canfield, 2011). A different and independent constraint potentially exists in the form of Fe or S isotopes (see e.g. Beard et al., 1999; Gomes and Johnston, 2017; van de Velde et al., 2018).

Any chemical element with multiple stable isotopes (Fe, for example, has four stable forms: ⁵⁴Fe, ⁵⁶Fe, ⁵⁷Fe, and ⁵⁸Fe) can potentially be used to track physicochemical processes that act to partition stable isotopes according to thermodynamic or kinetic principles. For Fe the most abundant isotopes are ⁵⁴Fe and ⁵⁶Fe, and deviations in the ${}^{\mathrm{56}}\mathrm{Fe}{/}^{\mathrm{54}}\mathrm{Fe}$ ratio of Fe-bearing aqueous and mineral phases from that of a reference material can be described using conventional δ notation:

$$\begin{array}{}\text{(35)}& {\mathit{\delta }}^{\mathrm{56}}\mathrm{Fe}=\left({\displaystyle \frac{\left(\frac{{}^{\mathrm{56}}\mathrm{Fe}}{{}^{\mathrm{54}}\mathrm{Fe}}\right)}{{\left(\frac{{}^{\mathrm{56}}\mathrm{Fe}}{{}^{\mathrm{54}}\mathrm{Fe}}\right)}_{\mathrm{ref}}}}-\mathrm{1.0}\right)\times \mathrm{1000},\end{array}$$

where $(\frac{{}^{\mathrm{56}}\mathrm{Fe}}{{}^{\mathrm{54}}\mathrm{Fe}}{)}_{\mathrm{ref}}$ is the isotope ratio of a standard reference material (IRMM-14). Any geochemical reaction, be it biotic (mediated by micro-organisms) or abiotic, can induce isotopic fractionation between co-occurring Fe- or S-bearing phases.

To model the isotopic signatures of Fe and S, we track the concentrations of the “bulk” pools (C_i) and the isotope-specific pools (⁵⁶C_i, for example, in the case of Fe). The isotopic signature of an Fe species C_i is then calculated as

$$\begin{array}{}\text{(36)}& {{\mathit{\delta }}^{\mathrm{56}}\mathrm{Fe}}_{C_i}=\left({\displaystyle \frac{\left(\frac{{}^{\mathrm{56}}{C}_i}{C_i{-}^{\mathrm{56}}{C}_i}\right)}{{\left(\frac{{}^{\mathrm{56}}\mathrm{Fe}}{{}^{\mathrm{54}}\mathrm{Fe}}\right)}_{\mathrm{ref}}}}-\mathrm{1.0}\right)\times \mathrm{1000}.\end{array}$$

Each individual reaction R_k is assigned a fractionation factor ${}^{\mathrm{56}}{\mathit{ϵ}}_{R_k}$ (‰; Table 4), which relates to ${}^{\mathrm{56}}{\mathit{\alpha }}_{R_k}$ as

$$\begin{array}{}\text{(37)}& {}^{\mathrm{56}}{\mathit{\alpha }}_{R_k}=\mathrm{1}+{\displaystyle \frac{{}^{\mathrm{56}}{\mathit{ϵ}}_{R_k}}{\mathrm{1000}}}.\end{array}$$

Isotope fractionation is then implemented by calculating an isotope-specific reaction rate (for the ⁵⁶Fe pool) from the bulk reaction rate R_k:

$$\begin{array}{}\text{(38)}& {}^{\mathrm{56}}{R}_k={\displaystyle \frac{{}^{\mathrm{56}}{\mathit{\alpha }}_{R_k}^{\mathrm{56}}{r}_{C_i}}{\mathrm{1}{+}^{\mathrm{56}}{\mathit{\alpha }}_{R_k}^{\mathrm{56}}{r}_{C_i}}}{R}_k,\end{array}$$

where

$$\begin{array}{}\text{(39)}& {}^{\mathrm{56}}{r}_{C_i}={\displaystyle \frac{{}^{\mathrm{56}}{C}_i}{C_i{-}^{\mathrm{56}}{C}_i}}.\end{array}$$

In this way, we assign a fractionation factor to each of the reactions considered in our model and are able to track the isotopic signature of each Fe and S species. This allows us to simulate the δ⁵⁶Fe and δ³⁴S values of dissolved species and solid mineral phases, both of which can potentially be compared to observations from the geological record.

The obvious limitation of assigning a constant fractionation factor to each reaction is that this is incapable of fully capturing natural isotopic variability. For instance, different microbial strains of sulfur reducers (or iron oxidizers) express different fractionation factors, even though the overall reaction remains the same (Gomes and Johnston, 2017; Pellerin et al., 2019). This is reflected in the often broad range of isotopic fractionation factors found in the literature for a given process (Table 4). Other factors influencing microbial fractionation are local environmental conditions, such as electron donor type and availability (Wing and Halevy, 2014; Pellerin et al., 2018) or evolutionary adaptation (Pellerin et al., 2015). Aside from biologically mediated transformations, kinetic effects associated with abiotic aqueous reactions and precipitation of solid phases also affect the fractionation that is eventually recorded in the end product. All these factors can make interpretation of isotope fractionations very complex in natural settings, and it is thus highly unlikely that any particular model simulation will be able to exactly reproduce observed isotope records. Nevertheless, the scheme employed here should be able to discern first-order observations from the geologic record and in some cases could potentially be used to rule out particular endmember hypotheses for ocean chemistry.

Table 4 Isotope fractionation factors.

[1] Kaplan and Rittenberg (1964), [2] Detmers et al. (2001), [3] Sim et al. (2011), [4] Gomes and Johnston (2017), [5] Pellerin et al. (2019), [6] Poser et al. (2014), [7] Fry et al. (1988), [8] Bottcher et al. (1998), [9] Price and Shieh (1979), [10] Beard et al. (1999), [11] Beard et al. (2003), [12] Johnson et al. (2004), [13] Bullen et al. (2001), [14] Crosby et al. (2007), [15] Rolison et al. (2018), [16] Guilbaud et al. (2011), [17] Wiesli et al. (2004).
^a The isotopic fractionation for iron oxide precipitation is driven by the oxidation reaction (FIO). ^b The isotope fraction factor for greenalite precipitation is currently unknown.

4 Model testing

We evaluate our model in two ways. We first test whether our extended model code is able to reproduce the Fe cycle of the contemporary oxygenated ocean by comparing it with the previously validated standard (oxic) Fe cycle in cGEnIE (Tagliabue et al., 2016). Secondly, we assess how the model behaves in an ocean which is predominantly anoxic. For contrasting with observations, we are severely limited because the modern ocean is largely well-oxygenated, with only a few oxygen-minimum zones near highly productive margins such as the Indian Ocean and the Peruvian margin (Keeling et al., 2010). However, even in these regions dissolved oxygen concentrations rarely reach zero, and the development of ferruginous or euxinic conditions is essentially absent. Only highly restricted basins such as the Cariaco Basin or the Black Sea and Baltic Sea can develop euxinic conditions, but these conditions arise as a result of local circulation within silled or enclosed basins and are not likely to be representative of an anoxic open-ocean setting. As a result, we lack observations to which we can directly calibrate our model. However, as our model development consisted of the implementation of well-established, mechanistic biogeochemical reactions with relatively well-defined kinetic rates, our model should simulate realistic rates and concentrations without extensive calibration of model parameters. We illustrate this by showing the spatial concentration and isotope features of a hypothetical anoxic ocean (Sect. 4.2) and – where possible – compare our predicted reaction rates to rates obtained from anoxic process analogues for the ancient oceans. Subsequently, we broadly illustrate the sensitivity of the model output to the newly introduced parameters of the iron–sulfur cycling (Sect. 4.3). As discussed below, our model is largely robust across a range of values for key parameters, and it predicts reaction and process rates that are comparable to those obtained experimentally or observed in modern analogue environments.

4.1 The contemporary ocean

We test our new iron scheme in comparison to the results of the oxic-only scheme presented by Tagliabue et al. (2016), which was based on the modern cGENIE configuration of Cao et al. (2009). For clarity, we carry out this test stepwise in two stages in order to separate out the consequences of resolving the different oxidation states of dissolved iron from the various solid-iron reactions.

1. Firstly, we substitute the original three tracers carried in the ocean – free dissolved iron (which we can equate here to ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{3}+}$), ligand-bound iron (${\mathrm{Fe}}_{\mathrm{ligand}}^{\mathrm{3}+}$), and free ligand (L_free) – for our new tracer scheme of Fe³⁺ (total = free + ligand-bound), Fe²⁺, and TDL (total dissolved ligand). Note that in the original scheme only two tracers were actually needed – total dissolved iron (Fe³⁺) and total dissolved ligand (L) – as the equilibrium between ${\mathrm{Fe}}_{\mathrm{free}}^{\mathrm{3}+}$, ${\mathrm{Fe}}_{\mathrm{ligand}}^{\mathrm{3}+}$, and L_free was re-calculated each time step. Hence, our new tracer scheme is composed of these same two primary tracers (Fe³⁺, L), plus Fe²⁺. We then enable the two-way oxidation–reduction reactions: ${\mathrm{Fe}}^{\mathrm{2}+}\to {\mathrm{Fe}}^{\mathrm{3}+}$ and ${\mathrm{Fe}}^{\mathrm{3}+}\to {\mathrm{Fe}}^{\mathrm{2}+}$ (Eqs. 13 and 18). We retain the same scavenging scheme as before, with iron scavenged by POM as Fe³⁺ and not subject to DIR.

2. The step is as per (1), but we now represent scavenged iron as FeOOH and enable DIR (Eq. 15). We also enable SMI of solid FeOOH with ΣH₂S (Eq. 18) and pyrite precipitation (Eq. 23). Note that the biogeochemical scheme used in Tagliabue et al. (2016) already included sulfate reduction and sulfide oxidation (DSR and CSO; Eqs. 32 and 33).

We find only minor differences between the old cGEnIE Fe cycling scheme (which only includes an oxic Fe cycle – see above) and the step (1) configuration of our new scheme (Fig. 6). Because the old scheme did not include reduced Fe²⁺, we compare total dissolved Fe (TDFe), which is the sum of all dissolved Fe phases (i.e. [TDFe] = [Fe³⁺] + [Fe²⁺]). Adding Fe²⁺ as a tracer, which is released during POC oxidation under anoxic conditions, leads to a slight increase in [TDFe] in oxygen-minimum zones (OMZs) (Fig. 6b) because in our new scheme Fe²⁺ is assumed not to be subject to scavenging. Hence, reducing Fe³⁺ to Fe²⁺ reduces iron loss due to scavenging and enhanced TDFe in OMZs. Consequently, there is a slightly enhanced supply of Fe to the surface waters via upwelling, leading to marginally higher surface [TDFe] (∼ 0.02 nM), in particular above the Peruvian OMZ (Fig. 6e). In step 2, by including a full Fe redox cycle (i.e. the reduction of FeOOH coupled to POC oxidation or sulfide oxidation) this effect is slightly intensified, with DIR and SMI leading to an additional reduction of Fe³⁺ (as FeOOH) and release of Fe²⁺ (which again escapes scavenging). (Fig. 6c, f).

Figure 6 Comparison of total dissolved Fe concentrations ([TDFe] = [Fe³⁺] + [Fe²⁺]) between (a, d) the old cGEnIE scheme (Tagliabue et al., 2016), (b, e) the new scavenging scheme, and (c, f) the full new scheme. (a–c) Zonally averaged TDFe concentrations, (d–f) surface TDFe concentrations (all units: nM). (b, c, e, f) Positive values indicate that the new scheme has higher concentrations, and negative values indicate that the old scheme has higher concentrations.

The small differences in [TDFe] translate to a negligible increase in global export production (Table 5) and indicates that the differences in performance between the old and new cGEnIE Fe cycle scheme are relatively trivial in a modern and nearly fully oxygenated ocean. For both the new and the old scheme, the sink of iron is scavenging of its oxidized form and no FeS₂ is formed (Table 5). Pyrite will only precipitate once FeS has passed its solubility threshold, which is ∼2 µM, and such concentrations are not reached in the modern ocean. Because the new scheme is more mechanistic, as it allows for reduction of Fe³⁺ in anoxic conditions, we believe our new extension will also be beneficial for studies of the modern iron cycle – in particular when addressing topics such as future ocean deoxygenation or intervals of the more recent geological past characterized by much more extensive OMZs than present.

Tagliabue et al. (2016)

Table 5 Comparison of selected model output for the old Fe scheme and the new Fe scheme.

n/a: not applicable.

4.2 An anoxic ocean

For testing and characterizing the Fe–S cycle model developments under anoxic conditions, we configure cGENIE with a single continent that runs from pole to pole (Figs. 2c, 3a). Each model experiment is initialized from a homogeneous and static ocean, with an imposed constant atmospheric O₂ concentration of 0.1 PAL (present atmospheric level, i.e. 21 000 ppm) and an atmospheric CO₂ concentration of ∼16 PAL (5337.6 ppm). These atmospheric boundary conditions are chosen as broadly plausible for the Precambrian Earth system and deliberately preclude the formation of sea ice in order to simplify the resulting spatial tracer patterns in the ocean and their mechanistic interpretation. The model is run in a “closed” configuration for all elements (notably C and P), in which the ocean–atmosphere inventory for each element is always conserved. We chose to keep a closed configuration for C and P because this allows us to fix important boundary conditions (such as pCO₂ and productivity) and look at the emerging ocean redox state under these conditions. An exception is made for Fe and S, which enter the surface ocean in the form of Fe²⁺ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and exit the ocean as FeOOH, FeS₂, FeCO₃, or Fe₃Si₂O₅(OH)₄. Fluxes of Fe and S are chosen to be in balance with respect to FeS₂ burial ($\mathrm{S}:\mathrm{Fe}=\mathrm{2}:\mathrm{1}$) and to represent the best estimate of present-day weathering fluxes to the ocean (excluding reprocessing in inner-shelf sediment settings) ($F_{{\mathrm{Fe}}^{\mathrm{2}+}}=\mathrm{1.3}\times {\mathrm{10}}^{\mathrm{12}}$ mol yr⁻¹, $F_{{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}}=\mathrm{2.6}\times {\mathrm{10}}^{\mathrm{12}}$ mol yr⁻¹; Poulton and Raiswell, 2002; Raiswell and Canfield, 2012). Hydrothermal systems represent another potentially important flux of Fe to the ocean (Tagliabue et al., 2010; Conway and John, 2014; Lough et al., 2019), and this flux is likely to be elevated when ocean chemistry is pervasively anoxic and relatively low in ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (Kump and Seyfried, 2005). Our simulations therefore also include a hydrothermal flux of Fe²⁺, broadly comparable to a plausible input flux to Proterozoic oceans (15.1×10¹² mol yr⁻¹; Thompson et al., 2019). This Fe²⁺ flux in our model set-up is equally distributed along a “hydrothermal ridge” located in the middle of the ocean (a straight line from k=9, l=3 to k=9, l=15; Fig. 2a). In order to balance the S:Fe flux ratio at 2:1 (see above), we must also specify a uniform surface flux of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ of 30.2×10¹² mol yr⁻¹. This is an extremely high S flux and is unlikely to be realistic on long timescales. However, it allows us to quickly diagnose spatial patterns in reducing Fe–S cycling at steady state, without a priori introducing a bias towards ferruginous or euxinic redox states. We initialize the model with an average ligand concentration of 1 nM (comparable to the modern ocean) and a semi-arbitrary marine phosphate inventory that is 50 % of today in order to reduce marine primary productivity (as productivity in the Precambrian ocean was likely lower than today; Crockford et al., 2018; Ozaki et al., 2019), and we run the ocean circulation to steady state for 20 000 years for each individual experiment and present model output from the 20 000th year of integration. It is important to emphasize that the idealized configuration we implement here is not meant to represent any specific period or event in Earth's history, but is rather meant to serve as a broadly plausible and computationally efficient set of boundary conditions for testing the extended model code.

The sea surface temperature and ocean circulation generated by our configuration of cGENIE are shown in Fig. 3. Sea surface temperatures vary from ∼40 ^∘C at the Equator to <10 ^∘C at the poles (Fig. 3c). The barotropic streamfunction shows a large degree of symmetry (Fig. 3d), whereas the overturning patterns are skewed to the Southern Hemisphere, with a strong anticlockwise circulation at around −60^∘ N (Fig. 3e). The overturning streamfunction shows strong upwelling at the Equator and deepwater mixing at the poles (Fig. 3e). The wind stress at the Equator drives surface waters towards the west, which will lead to stronger upwelling on the west side of the continent. The redox patterns discussed below will reflect these two main features: (i) upwelling at the Equator, in particular at the west side of the continent, and (ii) deepwater mixing at the poles.

In the following section, we briefly discuss the spatial model output and modelled redox cycling for our baseline configuration. Precipitation fluxes of oxidized Fe³⁺ are specified according to the scavenging scheme (Eq. 11). We focus here on (i) three depth slices of particle sinking fluxes of POC, FeOOH, FeS₂, and Fe₃Si₂O₅(OH)₄ (Fig. 7; FeCO₃ fluxes were negligible and are not shown); (ii) the average and three vertical longitudinal slices of O₂, dissolved Fe³⁺, Fe²⁺, and dissolved ΣH₂S (Fig. 8); and (iii) an overview of globally averaged reaction rates and iron mineral burial fluxes (Table 6).

The POC flux decreases with depth in the model from a maximum of ∼10 $\mathrm{mmol}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ immediately below the surface layer to near-zero values at the seafloor (Fig. 7a–c). The spatial pattern of solid-FeOOH flux (following the scavenging scheme described earlier – see Sect. 3.1.1) matches the POC flux pattern immediately below the surface layer, with higher values at the Equator and poles (dark blue shading in Fig. 7d). Most of this FeOOH is reduced in anoxic subsurface layers, and the flux declines with depth. However, at the poles where deep convection allows for greater oxygen (and Fe³⁺) penetration into the ocean interior (Fig. 8a–h), sinking particles have more time to continue to scavenge and accumulate Fe³⁺, and the flux increases with depth (Fig. 7e, f). The maximum FeOOH flux at the seafloor reaches ∼300 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (Fig. 7f), comparable to that observed in typical modern ocean margin sediments (see e.g. van de Velde and Meysman, 2016). Fluxes of FeS₂ reach their maximum near the seafloor (where a source of deep Fe²⁺ is supplied via the imposed hydrothermal flux), with maximum fluxes similar to those of FeOOH (Fig. 7g–i). The spatial pattern of Fe₃Si₂O₅(OH)₄ formation is very similar to FeS₂, but the fluxes are several orders of magnitude lower (Fig. 7j–l). Fluxes of FeCO₃ were near zero everywhere (data not shown), consistent with recent work suggesting that water-column precipitation of FeCO₃ is difficult to achieve, even in iron-dominated oceans (Jiang and Tosca, 2019; Tosca et al., 2019). The negligible FeCO₃ fluxes are at odds with Fe-speciation data for Precambrian rocks showing that an important fraction of sedimentary Fe consists of reduced non-sulfurized Fe minerals (Sperling et al., 2015). Our results indicate that these minerals are most likely formed during sedimentary diagenesis, emphasizing the potential importance of processes below the sediment–water interface in structuring Fe-speciation signals in Earth's rock record. Future development will thus include a representation of sedimentary Fe cycling in the sedimentary module OMEN-SED (Hülse et al., 2018).

Figure 7 Model output of (a)–(c) particulate organic carbon flux (J_POC), (d)–(f) solid-iron oxide flux (J_FeOOH), (g)–(i) pyrite flux ($J_{{\mathrm{FeS}}_{\mathrm{2}}}$), and (j)–(l) greenalite flux ($J_{{\mathrm{Fe}}_{\mathrm{3}}{\mathrm{Si}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}(\mathrm{OH}{)}_{\mathrm{4}}}$) for an ocean with 50 % of the modern phosphate inventory. J_POC is given in $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, J_FeOOH and $J_{{\mathrm{FeS}}_{\mathrm{2}}}$ are given in $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, and $J_{{\mathrm{Fe}}_{\mathrm{3}}{\mathrm{Si}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}(\mathrm{OH}{)}_{\mathrm{4}}}$ is given in $\mathrm{nmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. The siderite flux ($J_{{\mathrm{FeCO}}_{\mathrm{3}}}$) was near zero everywhere and is not shown here.

We additionally find that even in an idealized ocean with a simple symmetrical continental configuration, complex spatial patterns emerge in the Fe–S redox chemistry (Fig. 8). For instance, the poles are more well-ventilated, allowing oxygenated conditions and persistence of Fe³⁺ to a few kilometres of depth in our benchmark simulation (Fig. 8a–h). The concentrations of Fe³⁺ are much higher than the nanomolar and picomolar concentrations we observe in the ocean today (Tagliabue et al., 2016), which is a consequence of much higher rates of iron delivery to the surface ocean (through upwelling of reduced Fe²⁺) that allows Fe³⁺ to accumulate to higher concentrations before it is eventually scavenged. Our model predicts concentrations in the hundreds of nanomolar range (up to micromolar at the oxic–anoxic interface), which compares well to oxic water layers overlying anoxic deep water (Taillefert and Gaillard, 2002; Crowe et al., 2008). Note that this Fe³⁺ is likely in colloidal or nanoparticulate form and not truly dissolved. At eastward latitudes, deep convection at the poles is less intense, and upwelling on the eastward edge of the ocean leads to higher export production, which subsequently leads to build-up of reduced Fe²⁺ and ΣH₂S (Fig. 8l and p). Dissolved Fe²⁺ reaches higher concentrations in the deeper ocean, largely as a result of deep hydrothermal inputs (Fig. 8i–l), whereas the highest ΣH₂S concentrations are spatially constrained to areas of more intense POC degradation along the Equator, just below the oxic zone (Fig. 8m–p).

Figure 8 Model output of (a)–(d) dissolved oxygen (O₂), (e)–(h) dissolved ferric iron (Fe³⁺), (i)–(l) total dissolved ferrous iron (Fe²⁺), and (m)–(p) total dissolved sulfide (ΣH₂S) for an ocean with 50 % of the modern phosphate inventory (all concentrations: µM).

Because pervasive anoxia is not present in modern open-ocean environments (see above), evaluation of the realism of our globally integrated reaction rates must rely on comparison to modern process analogues for ancient oceans (e.g. anoxic lake and restricted marine systems). Though we consider this a valid approach, it must be borne in mind that the transport processes in particular are very different in stratified lacustrine and marine systems, and there are reasons to assume that the biogeochemical dynamics of these systems will not strictly map onto pervasively anoxic open-ocean environments. However, though the absolute values are not directly comparable, we can qualitatively compare our results with rates derived from anoxic lake systems. For instance, our default model suggests that DIR is a negligible contributor to POC mineralization (Table 6), consistent with previous observations from ferruginous Lake Matano and the euxinic Black Sea (Konovalov et al., 2006; Crowe et al., 2011). However, Crowe et al. (2011) found that rates of DIR were roughly an order magnitude lower than those of methanogenesis, whereas we find that DIR is 5 orders of magnitude less important than all other mineralization pathways (Table 6). One possible reason for this discrepancy is the importance of sediment recycling in the natural lake system, which our model does not currently represent. Due to the shallower water column (e.g. ∼300 m in Lake Matano versus 5000 m in our idealized ocean), there is a much stronger coupling between sedimentary processes (i.e. the recycling of Fe as a benthic flux) and water-column processes. Indeed, when the total rate of DIR is corrected for iron recycling occurring in sediments within Lake Matano the importance of DIR decreases by several orders of magnitude (Crowe et al., 2008). Additionally, cGENIE treats POC as a concentration assumed to be uniformly dispersed throughout a grid cell, meaning that POC and FeOOH concentrators used in reaction calculation will not take into account the reality of particulate matter being locally aggregated and highly concentrated as it sinks.

Table 6 Globally integrated reaction rates and burial fluxes for our benchmark simulation.

Our model predicts a globally integrated FIO rate of $\sim \mathrm{0.86}\times {\mathrm{10}}^{-\mathrm{7}}$ $\mathrm{mol}\mathrm{Fe}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (Table 6), which is lower than rates estimated for anoxic lakes (12–$\mathrm{51}\times {\mathrm{10}}^{-\mathrm{6}}$ $\mathrm{mol}\mathrm{Fe}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$; Crowe et al., 2008; Walter et al., 2014). However, rates measured near the oxycline ($\mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}}$ $\mathrm{mol}\mathrm{Fe}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) are comparable to rates measured in anoxic lake systems. Interestingly, more than half of the ΣH₂S is re-oxidized via SMI (Table 6), which indicates that Fe is able to act as a relatively efficient intermediate between ΣH₂S oxidation and O₂ reduction. This also occurs in the Black Sea, where metal oxides are responsible for ∼60 % of the sulfide re-oxidation (Konovalov et al., 2006).

Isotope patterns

Figure 9a–e show the modelled stable Fe isotope patterns for all key Fe-bearing dissolved and solid-phase species (with the exception of FeCO₃, which is a negligible component in our benchmark simulation). In our model simulations, all dissolved Fe that enters the ocean is assigned an isotope signature of 0.0 ‰. This allows us to observe the isotope fractionation of all Fe phases relative to the Fe that entered the ocean (Fig. 10a). The dissolved iron phases show similar isotopic signatures (∼1.1 ‰; Figs. 9a, b, 10a), which can be explained by the large amount of isotopically light FeS₂ burial (Fig. 10a) that drives $\mathit{\delta }{}^{\mathrm{56}}\mathrm{Fe}-{\mathrm{Fe}}^{\mathrm{2}+}$ to heavier values. Note that the very heavy ${\mathit{\delta }}^{\mathrm{56}}\mathrm{Fe}-{\mathrm{Fe}}^{\mathrm{3}+}$ values only occur in the ocean interior, where the concentration of Fe³⁺ is virtually zero (Fig. 8e–h). The isotope signature of buried oxidized iron (FeOOH) is around 1.8 ‰ heavier than the Fe that entered the ocean, whereas the major burial fraction of reduced iron (FeS₂) has an isotope signature that is ∼1 ‰ lighter, which reflects their relative importance as an Fe burial phase (Figs. 9c, e, 10a). These isotope values are broadly comparable to phase-specific stable Fe isotope observations from ancient sedimentary rocks (Heard and Dauphas, 2020). Similarly, Fig. 9f–g show the modelled stable S isotope patterns for all key S-bearing dissolved and solid-phase species.

All dissolved S that enters the ocean is assigned an isotope signature of 0.0 ‰. Sulfate is isotopically enriched relative to its input value (∼3.5 ‰; Figs. 9f, 10b), and this difference compares well to the geological record (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ is ∼5 ‰ heavier than its input value; Canfield and Farquhar, 2009). Additionally, free sulfide is isotopically lighter than sulfate ($\sim -\mathrm{5}$ ‰ on a global scale; up to $\sim -\mathrm{15}$ ‰ locally), while buried FeS₂ expresses an isotope signature of $\sim -\mathrm{2.38}$ ‰ (Fig. 10b). Our baseline simulation thus suggests that the expressed isotope fractionation between ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and FeS₂ is only around 6.5 ‰, which is roughly consistent with what is observed in the geological record (Canfield and Farquhar, 2009). Therefore, we conclude that our model has strong potential for tracking iron and sulfur isotope signatures for comparison with Earth's rock record.

Figure 9 (a–e) Modelled stable Fe isotope signatures of key dissolved and solid-phase Fe species. Shown are zonally averaged values for (a) dissolved ferric Fe (${\mathit{\delta }}^{\mathrm{56}}\mathrm{Fe}-{\mathrm{Fe}}^{\mathrm{3}+}$) and (b) dissolved ferrous iron (${\mathit{\delta }}^{\mathrm{56}}\mathrm{Fe}-{\mathrm{Fe}}^{\mathrm{2}+}$), as well as the isotope compositions at the seafloor of (c) iron oxides (δ⁵⁶Fe−FeOOH), (d) greenalite (δ⁵⁶Fe−Fe₃Si₂O₅(OH)₄), and (e) pyrite (δ⁵⁶Fe−FeS₂) for an ocean with 50 % of the modern phosphate inventory. All values are in per mille relative to IRMM-14. (f–g) Modelled stable S isotope signatures of key dissolved and solid-phase S species. Shown are zonally averaged values for (f) sulfate (${\mathit{\delta }}^{\mathrm{34}}\mathrm{S}-{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) and (g) dissolved free sulfide (δ³⁴S−ΣH₂S), as well as (h) the isotope composition at the seafloor of pyrite (δ³⁴S−FeS₂). All values are in per mille relative to VCDT.

Figure 10 Isotope mass balance for the (a) ⁵⁶Fe and (b) ³⁴S systems from our baseline simulation. Purple arrows are influxes, and black arrows are outfluxes. Bubbles represent the whole ocean inventory. Note that in our simulation, the S system is not in steady state (because FeS₂ burial is too low).

4.3 Sensitivity analysis

We evaluate model output sensitivity to four key parameters of our modelled iron–sulfur cycle (K_0,FeOOH, K_i,FeOOH, k_SMI,s, k_PyP) using the elementary effect test (EET; Morris, 1991). These parameters were chosen because they are either unconstrained by laboratory experiments (K_0,FeOOH, K_i,FeOOH; Sect. 3.1.2), represent a complex mixture of different iron minerals with different reactivities (k_SMI,s; Poulton et al., 2004), or are expected to have a strong and potentially difficult to forecast influence on other reactions (k_PyP; van de Velde et al., 2020b). Other parameters introduced in the model are either relatively well-constrained by laboratory studies and calibrated on field data (k_CSO, k_FIO; Millero et al., 1987a, b; Ridgwell et al., 2007), have been calibrated extensively in previous work (k_scav; Tagliabue et al., 2016), or are likely to be of secondary importance to the iron–sulfur cycle (k_AFC, b_AFC, k_greenalite, b_greenalite).

The EET method estimates global sensitivity by calculating the mean of r finite differences (“elementary effects”) (Pianosi et al., 2016):

$$\begin{array}{}\text{(40)}& \begin{array}{rl}{S}_i& ={\displaystyle \frac{\mathrm{1}}{r}}\sum _{j=\mathrm{1}}^r{\text{EE}}^j={\displaystyle \frac{\mathrm{1}}{r}}\sum _{j=\mathrm{1}}^r\\ & \cdot {\displaystyle \frac{f({\overline{x}}_{\mathrm{1}}^j,\mathrm{\dots },{\overline{x}}_i^j+{\mathrm{\Delta }}_i^j,\mathrm{\dots },{\overline{x}}_M^j)-f({\overline{x}}_{\mathrm{1}}^j,\mathrm{\dots },{\overline{x}}_i^j,\mathrm{\dots },{\overline{x}}_M^j)}{{\mathrm{\Delta }}_i^j}}{c}_i,\end{array}\end{array}$$

where $x_i^j$ represents the jth value of the ith parameter, ${\mathrm{\Delta }}_i^j$ represents the variation of the ith parameter, f() is the model output for a given set of parameters, and c_i is a scaling factor. A higher mean value of S_i indicates that a given model output is more sensitive to variations in parameter i. The standard deviation can also be calculated, with a high standard deviation indicating that a parameter interacts with others because its sensitivity changes across the variability space (inlay in Fig. 11; Pianosi et al., 2016).

Table 7 Sensitivity range and baseline values of each parameter tested. Note that the units are not identical to the units used in the user configuration files for the cGENIE model. A table with the parameter units converted to cGEnIE units is included in Appendix A.

We use the EET method, as implemented within the Sensitivity Analysis For Everyone (SAFE) toolbox (Pianosi et al., 2015), to investigate the four chosen model parameters across the ranges specified in Table 7. We vary the limitation and inhibition constants K_0,FeOOH and K_i,FeOOH over 6 orders of magnitude around the default (baseline) value due to the high uncertainty associated with these parameters. The lower bound of the reaction constant of solid iron oxide with sulfide (k_SMI,s) is defined as the reactivity of hematite ($\mathrm{5.34}\times {\mathrm{10}}^{-\mathrm{3}}$ ${\mathrm{M}}^{-\mathrm{0.5}}{\mathrm{h}}^{-\mathrm{1}}$; Poulton et al., 2004), and the upper bound is taken to be 5 orders of magnitude higher (5.34×10² ${\mathrm{M}}^{-\mathrm{0.5}}{\mathrm{h}}^{-\mathrm{1}}$, which is comparable to freshly precipitated hydrous ferric oxide). We test the pyrite precipitation constant (k_Pyp) across a range between $\mathrm{0.3708}\times {\mathrm{10}}^{-\mathrm{3}}$ ${\mathrm{M}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$ and 0.3708×10² ${\mathrm{M}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$, which corresponds to the range of kinetic constants commonly used in diagenetic models (see e.g. Van Cappellen and Wang, 1996; Meysman et al., 2003; Dale et al., 2009; van de Velde et al., 2020a). Our sensitivity ensemble consists of 100 individual model experiments using a Latin hypercube sampling approach (using the SAFE toolbox; Pianosi et al., 2015) to select random starting points $x^j(j=\mathrm{1},\mathrm{\dots };r)$ and parameter variations Δ_i. For more information on the sampling strategy we refer the interested reader to Campolongo et al. (2011).

Figure 11 Sensitivity analysis of the four key parameters of the iron–sulfur cycle for a range of different outputs. Inset shows how to read the graph; points that plot more to the right indicate that the specific output is more sensitive to changes in parameter values, and points that plot higher indicate that the parameter is more interactive with other parameters. Data processing was done with the SAFE toolbox of Pianosi et al. (2015).

The EET analysis suggests that changes in K_0,FeOOH can significantly impact DIR, whereas the inhibition constant K_i,FeOOH has a relatively minor impact on model output (Fig. 11). This is expected, as both parameters act only on the mineralization pathways (Eq. 10). It is also consistent with some of the literature on anoxic systems, which suggests that in many cases the majority of iron reduction in the water column is coupled to sulfide (Mikucki et al., 2009). Even though both K_0,FeOOH and K_i,FeOOH are the least constrained by experimental results (see Sect. 3.1.2), the EET analysis indicates that they are relatively unimportant for the overall model output despite their impact on DIR. In contrast, k_PyP and k_SMI,s both exerted a more notable impact on model output across a range of diagnostics (Fig. 11). In particular, k_PyP had a significant impact across all model output diagnostics analysed here in terms of both model sensitivity and interactivity with other parameters (Fig. 11). Because pyrite precipitation controls the inventories of both dissolved Fe²⁺ and dissolved ΣH₂S, reducing or increasing the kinetic precipitation parameter will affect the ambient concentrations, re-oxidation pathways, and eventual mineral burial for all phases across the Fe–S system.

Unfortunately, an elementary effect test gives no quantitative metric to evaluate the magnitude with which a parameter affects overall model outcome. Therefore, to illustrate the quantitative impact of the possible parameter choices, we ran a separate set of experiments in which we changed a parameter from its baseline value to its lower and upper bound whilst keeping the other three parameters at their baseline values (Table 7).

Figure 12 Global model output for the three different values of the four key parameters (K_0,FeOOH, K_i,FeOOH, k_PyP, and k_SMI,s). Low, baseline, and high values are given in Table 7. Panels (a)–(c) and (h)–(n) are in moles per year (mol yr⁻¹), and panels (d)–(g) are in moles per kilogramme (mol kg⁻¹).

Figure 12 reveals that AR, DSR, mean surface O₂ concentrations, and CSO are relatively insensitive to parameter choices (Fig. 12a, c, d, i). Dissimilatory iron reduction is only important when K_0,FeOOH is set at its lowest value, and even then it is 3 orders of magnitude smaller than AR or SR (Fig. 12a–c). Consistent with the EET analysis (Fig. 11), K_0,FeOOH and K_i,FeOOH have no influence on the model output aside from the magnitude of DIR, which in itself is of less importance than the mineralization pathway (Fig. 12b and Sect. 4.2). In contrast, changes in the pyrite precipitation constant k_PyP impact several model outputs. When pyrite precipitation rates are elevated, Fe²⁺ and ΣH₂S are removed from local seawater more rapidly, which results in the following:

1. a decreased build-up of Fe²⁺ and ΣH₂S in the ocean interior (Fig. 12f, g), leading to

2. a decrease in both aerobic and anaerobic re-oxidation pathways (FIO, SMI_d, and SMI_s; Fig. 12h–k), which then results in

3. more O₂ available for aerobic respiration at the expense of less thermodynamically favourable electron acceptors (i.e. AR increases and DSR decreases; Fig. 12a, c),

4. a decrease in Fe³⁺ concentrations in the surface ocean (Fig. 12e), and

5. less burial of FeOOH, increasing the burial of reduced iron minerals (Fig. 12l–n).

The effect on greenalite burial is non-linear (Fig. 12o), which indicates that at higher values of k_PyP, pyrite precipitation is competing with greenalite precipitation. Overall, our sensitivity analysis suggests that k_PyP is an important parameter for the model output and should be chosen with care. Fortunately, pyrite precipitation has been well-studied in laboratory experiments (Rickard, 1997, 2006), with the result that our baseline value for this parameter is relatively well-constrained.

Figure 13 Global model output for a range of k_SMI,s values (units are ${\mathrm{M}}^{-\mathrm{1}/\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$). Lines indicate experimental reactivity parameters as presented in Poulton et al. (2004). Panels (a)–(c) and (h)–(n) are in moles per year (mol yr⁻¹), and panels (d)–(g) are in moles per kilogramme (mol kg⁻¹). Note that the x axis is logarithmic.

The fourth parameter, which also influences the model output, is the reactivity parameter of solid iron oxides (k_SMI,s). Here, the parameter choice is more complex. Laboratory experiments have shown that different iron oxide minerals exhibit a wide reactivity range (spanning several orders of magnitude) (Canfield et al., 1992; Poulton et al., 2004). Therefore, we explore the sensitivity of this parameter in more detail using a range of measured reactivity constants by Poulton et al. (2004) whilst keeping all other parameters at their baseline values (Fig. 13). Increasing the reactivity of particulate FeOOH does the following:

1. increases the anaerobic re-oxidation reaction of sulfide with FeOOH (SMI_s) at the expense of the aerobic re-oxidation reaction (CSO) (Fig. 13i, j), which then leads to

2. an increase in the Fe²⁺ inventory and a decrease in the ΣH₂S inventory (Fig. 13f, g); more Fe²⁺ leads to

3. more FIO and thus a higher surface Fe³⁺ concentration and more re-oxidation of ΣH₂S with dissolved Fe³⁺ (SMI_d; Fig. 13e, h, k).

4. Because of the higher reactivity of the FeOOH particles, less FeOOH is buried, and more reduced Fe minerals are buried (Fig. 13l–o).

Although it is clear that changing k_SMI,s impacts model output, the overall magnitude of the effect is moderate when compared to changing k_PyP. Nevertheless, the choice of k_SMI,s is critical. We choose a baseline reactivity value (Table 3) comparable to lepidocrocite for several reasons. Firstly, we assume all Fe³⁺ that is not scavenged represents a “colloidal” pool, with a reactivity similar to that of hydrous ferric oxide. When Fe³⁺ becomes scavenged (and is thus in solid state), it has likely undergone some ageing, and it will be less reactive than hydrous ferric oxide. Secondly, any FeOOH that does not react in the water column will end up in the sediment and will, at least in part, be recycled back to the water column (even under oxic conditions; Dale et al., 2015). Our model currently lacks a sedimentary iron cycle (see Sect. 5) and would thus tend to underestimate the overall importance of the iron cycle were we to select a reactivity constant that is too low. Finally, field evidence suggests that FeOOH that is freshly precipitated is highly reactive (Picard et al., 2015; Beam et al., 2018), and thus iron precipitating from the surface ocean is expected to react on relatively short timescales. Any FeOOH minerals that would resist reduction passing through a sulfidic water column are likely unreactive and are thus presumably inert on early diagenetic timescales. Indeed, iron oxide minerals in sediments underlying euxinic water columns tend to show no depth trend, indicating that very little iron reduction is occurring at depth in such systems (see e.g. Xiong et al., 2019). Taken together, these observations support the presumption that once Fe³⁺ is scavenged, its reactivity is less than hydrous ferric oxide but higher than that of goethite, motivating our default k_SMI,s value. Nevertheless, the value used for this parameter should be considered carefully depending on the model application and assumed boundary conditions.

5 Outlook and conclusions

The principal aim of this paper is to provide a detailed description of our extension of the cGENIE biogeochemistry module to include coupled anoxic Fe and S biogeochemical cycles. Because direct tuning with actual measured concentrations and rates is not possible, we have relied heavily on kinetic constants and solubility values extracted from laboratory incubations. While care should be exercised in the application of these kinetic constants to reactions under in situ conditions, our sensitivity analysis indicates that our key model results are robust across a wide range of possible parameter values. In addition, our proposed baseline parameterization yields reaction rates, concentrations, burial fluxes, and stable Fe isotope compositions that broadly compare well to both field measurements of process analogues for ancient ocean systems (i.e. anoxic lakes) and observations from the geologic record. Therefore, we believe that our model description of the anoxic Fe–S cycle is a valuable tool and an important step forward in simulating ocean redox landscapes during periods of Earth's history in which the ocean interior was pervasively anoxic. However, below we highlight some important limitations to our current model architecture and also give some examples of possible future developments.

Most notably, our model is currently unable to resolve any sedimentary processes that would contribute to the global iron cycle (see e.g. Dale et al., 2015). In particular, the model does not include a representation of benthic iron reduction and recycling back into the water column. Building on the improvements in the biogeochemistry of cGENIE described here, in the future we plan to extend the organic-matter-enabled sediment component of cGENIE (OMEN-SED; Hülse et al., 2018) to include an explicit representation of the benthic iron cycle. We anticipate that this will both improve the realism of tracer fields within the ocean interior and make comparisons between predicted sedimentary signals and observations from Earth's sedimentary rock record more accurate and robust. Second, there are likely to be important mechanistic links between the biogeochemistry of Fe and S within the ocean and the local and global recycling and bioavailability of key nutrient species for the biosphere (Bjerrum and Canfield, 2002; Laakso and Schrag, 2014; Jones et al., 2015; Reinhard et al., 2017). Future work will thus also focus on explicitly coupling the anoxic Fe and S biogeochemistry to the phosphorus (P) and nitrogen (N) cycles, in particular the scavenging and remobilization of P under different redox states and the impact of dissolved Fe availability on nitrogen fixation. The modularity of cGENIE also allows the substitution of an explicit plankton ecological model (ECOGEM) for the default biological export scheme (Ward et al., 2018), enabling the exploration of feedbacks between marine ecosystems, nutrient availability, and ocean redox conditions (Reinhard et al., 2020b). Lastly, future work will seek to include other redox-sensitive proxies and bio-essential elements such as molybdenum, uranium, or vanadium (Tribovillard, 2006) within the model code, which will further extend the applicability of our model and help validate it against observations from modern anoxic systems and the geologic record.

Appendix A

Table A1 List of kinetic constants for the reactions included in the cGENIE model – converted to cGEnIE units.

Table A2 Sensitivity range and baseline values of each parameter tested – converted to cGEnIE units.