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 intermediatecomplexity 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 δ56Fe and δ34S 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.

Abstract. 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 intermediatecomplexity 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 δ 56 Fe and δ 34 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.

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 2− 4 ) 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 or-ganic matter, microbial reduction of SO 2− 4 and FeOOH produces reduced sulfide (H 2 S) and ferrous iron (Fe 2+ ), respectively. When present at the same location, H 2 S and Fe 2+ combine into iron monosulfides and eventually pyrite (FeS 2 ) (Rickard, 1997(Rickard, , 2006, the burial of which couples the shortterm surface cycles of Fe and S with their long-term geological cycles (Berner, 1989). Depending on the relative ocean inventory of H 2 S vs. Fe 2+ , the precipitation of FeS 2 can lead to an anoxic water body becoming either iron-rich ("ferruginous") or sulfide-rich ("euxinic") (Canfield, 1998;Poulton and Canfield, 2011), with H 2 S : Fe 2+ ratios greater than 2 : 1 (the S : Fe ratio in FeS 2 ) 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., 2013Reinhard et al., , 2017Wallmann 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 2 S or Fe 2+ 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 2+ (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 2 (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 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 2+ or H 2 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).
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 (PO 3− 4 ) 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- 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. 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 onedimensional 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 .
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, andRidgwell 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 3 ), greenalite (Fe 3 Si 2 O 5 (OH) 4 ), and pyrite (FeS 2 ) 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, oceanatmosphere and ocean-sediment exchange, and the longterm (geological) cycle of carbon and various solid-Earthderived tracers 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 , 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 dynamicthermodynamic 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.

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 −2 ). 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).

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 , 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 (PO 3− 4 ) 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 −1 ) follows Doney et al. (2006).
Rates of photosynthetic nutrient uptake are scaled to ambient dissolved nutrient concentrations (PO 3− 4 and Fe 3+ ) 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 (red P/Fe ). The various limitation terms (all unitless) are 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 −2 , κ PO 3− 4 =0.1 µM, κ Fe 3+ = 0.1 nM) and are as used in Tagliabue et al. (2016). The influence of temperature on biological export is parameterized as where k T 0 (0.59) is a scaling constant, k eT (15.8 • C) the efolding 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 . A proportion (ν = 0.66) of PO 3− 4 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 (mol m −2 h −1 ) is always equal to the rate of PO 3− 4 uptake: where ρ is the density of seawater (kg m −3 ) and h e the thickness of the euphotic zone (80.84 m). The particulate organic carbon (POC) export flux (mol m −2 h −1 ) is calculated using a fixed Redfield ratio as 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): 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 3 in deep-sea sediments (Jamie Wilson, personal communication).

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 bio- logical 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. 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 . 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 Martintype decay (Fig. 4a) that determines the mineralization rate (R min ) at each depth layer in the water column, and thus R min 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 , while methanogenesis has been more recently added ). 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.  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).

Reaction number Reaction name Abbreviation Equation
Primary redox reactions Here, K 0 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. As a consequence of organic matter remineralization in the model, DIR produces ferrous iron (Fe 2+ ), which is reoxidized 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 H 2 S = H 2 S+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 2 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 3+ ), will precipitate out as iron oxide (FeOOH) minerals (IrP, R8 in Table 1).
When Fe 2+ and H 2 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 2 S to form the mineral pyrite (FeS 2 ) via pyrite precipitation (PyP, R9 in Table 1) (Rickard, 1997). Alternatively, when Fe 2+ accumulates past its saturation state with bicarbonate, siderite (FeCO 3 ) forms via siderite precipitation (SiP, R10 in Table 1) (Jimenez-Lopez et al., 2004;Jiang and Tosca, 2019). Finally, greenalite (Fe 3 Si 2 O 5 (OH) 4 ) precipitates when Fe 2+ and dis-solved silica (SiO 2 ) 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 solidiron 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 (1/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.

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 2 ). 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.

The oxic iron cycle
In oxygenated waters, iron is predominately present in its oxidized form (Fe 3+ ). 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 3+ consists of free 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.

Reaction number Reaction name
Kinetic expression Reference Fe 3+ free and ligand-bound Fe 3+ ligand , 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 Fe 3+ free together with the magnitude of the POC flux from the grid cell, i.e.
where k scav (1.43 × 10 −6 mol −1 m 2 ) 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 Fe 3+ free at each time step by the conservation equations.
Here, K FeL sp is the stability constant of the L − Fe 3+ complex (1.0 × 10 11 M −1 ; Table 3), L free is ligand unassociated with Fe 3+ , 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 3+ abundance to picomolar values without stabilization by ligands (Liu and Millero, 2002). Thus, Fe 3+ free 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 Fe 3+ free can be thought of as a mix of purely dissolved Fe 3+ and a range of colloidal and nanoparticulate Fe 3+ phases that do not settle efficiently through the water column. The Fe 3+ free 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, Fe 3+ free (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 Fe 3+ free 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 Fe 3+ free 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 2+ . Any reduced Fe 2+ that is mixed into the oxic zone is being rapidly oxidize to Fe 3+ via ferrous iron oxidation (FIO): where the rate equation can be expressed as and where k FIO is a reaction rate constant (k FIO = 0.115 × 10 6 M −1 h −1 ) ( Table 2; Millero et al., 1987a). In the calculation of the nutrient limitation of biological production at the surface, Fe 2+ is also assumed to be bioavailable, meaning that iron limitation is calculated from the total dissolved iron pool, equal to Fe 2+ + Fe 3+ free + Fe 3+ ligand . 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 ligandbound iron phase, to now distinguish between the oxidation states of iron (and the tracers Fe 2+ and Fe 3+ ), with Fe 3+ being split (as previously) into a free (and scavengeable) and ligand-bound phase.

The anoxic iron cycle
In an anoxic water column, oxidized forms of iron can be reduced to Fe 2+ 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 3+ , 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 3+ 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 (15) where it is implicitly assumed that every iron oxide particle consists solely of Fe 3+ 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 ∼ 30 × 10 −3 M (Van Cappellen and Wang, 1996;Thamdrup, 2000). This limitation is expressed as FeOOH] 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 FeOOH] 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 −3 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 sulfurmediated iron reduction (SMI), which follows the stoichiom-etry (Poulton et al., 2004) We are not explicitly modelling elemental sulfur (S 0 ) but assume that it becomes quantitatively disproportionated into H 2 S and SO 2− 4 (Finster et al., 1998), which then leads to the overall stoichiometry  , 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 3+ (in the water column). To achieve an excess of Fe 3+ over H 2 S, the H 2 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) 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 SMI,d = 2.64 × 10 2 M −0.5 h −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 SMI,s = 1.98 × 10 0 M −0.5 h −1 (Poulton et al., 2004). Lepidocrocite is less crystalline than the other (nonhydrous) iron oxides and precipitates when the rate of Fe 3+ 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 3+ produces Fe 2+ , which can either be reoxidized when it comes in contact with O 2 (see Sect. 3.1.1) or can form reduced minerals.
Both reduction and oxidation reactions of dissolved iron, even with dissolved O 2 and H 2 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 2+ and Fe 3+ according to ambient dissolved oxygen and sulfide concentrations.

Reduced iron mineral formation
Reduced Fe 2+ can form complexes with a number of inorganic ligands that are common in seawater, including Cl − , SO 2− 4 , 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 Fe 2+ free and the aqueous ironsulfide 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 which was obtained from the visualMINTEQ database (Gustafsson, 2019) and is based on stability constants calculated by Luther et al. (1996). Since dissolved Fe 2+ free and the FeS aq complex are the dominant dissolved forms of reduced ferrous iron in natural waters (in anoxic seawater devoid of sulfide, dissolved Fe 2+ free still represents ∼ 80 % of the dissolved Fe 2+ 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 3+ ; see Sect. 3.1.1). This is important, since FeS aq forms particulate FeS p , which is the precursor of pyrite (FeS 2 ) (Rickard, 1997(Rickard, , 2006. Furthermore, siderite or greenalite can only precipitate from Fe 2+ free (Tucker, 2001;Tosca et al., 2015), and their precipitation rates are consequently dependent on the dissolved Fe 2+ free concentration. Before each precipitation reaction, we calculate the concentration of FeS aq , Fe 2+ free , and H 2 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 2 S, with the production of free hydrogen gas (as S 2− in H 2 S 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). and FeS p has to be oxidized to S 1− in FeS 2 ; Rickard, 1997): Since hydrogen gas is highly reactive, it is almost instantaneously oxidized with the reduction of an electron acceptor like O 2 , SO 2− 4 , or CO 2 . Given that pyrite precipitates under sulfidic conditions, oxygen is absent, and sulfate is the most likely electron acceptor.
To preserve the redox balance while avoiding the need to model an additional state variable (H 2 ), we can combine Reactions (21) and (22) By using Eq. (23) we implicitly assume that the hydrogen gas produced during pyrite precipitation immediately reacts with SO 2− 4 . 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 2 S], with k PyP = 0.3708 × 10 0 M −1 h −1 (Rickard, 1997). The rate equation for FeS 2 precipitation can then be written as Any Fe 2+ that is not complexed with sulfide (Fe 2+ free ) can form siderite (FeCO 3 ), an iron-carbonate mineral (Tucker, 2001).
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 FeCO 3 (OH) 1/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): 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): where k AFC is 1.963 × 10 −14 M h −1 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 3 Si 2 O 5 (OH) 4 ) (Tosca et al., 2015).
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: where γ represents the activity coefficient of an ion in seawater (Table 3) and K greenalite sp its solubility product (Rasmussen et al., 2015): where K greenalite sp equals 3.98×10 27 M −1 (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): where k greenalite is 6.996 × 10 −13 M h −1 and b greenalite equals 1.856 (Table 3).

The sulfur cycle
The oxidized form of dissolved sulfur in the ocean is sulfate (SO 2− 4 ). Under anaerobic conditions, sulfate is used as a terminal electron acceptor during organic matter remineralization: The produced sulfide can then be re-oxidized when it comes into contact with oxygen via canonical sulfide oxidation (CSO).
The rate of CSO in cGENIE is dependent on the concentration of the electron donor ( H 2 S) and acceptor (O 2 ) as well as a second-order rate constant k CSO = 0.625 × 10 6 M −2 h −1 (Zhang and Millero, 1993): Alternatively, H 2 S can be re-oxidized with Fe 3+ (either Fe 3+ 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 4 ) 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 [SO 2− 4 ] 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 cGE-NIE 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).

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: 54 Fe, 56 Fe, 57 Fe, and 58 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 54 Fe and 56 Fe, and deviations in the 56 Fe/ 54 Fe ratio of Fe-bearing aqueous and mineral phases from that of a reference material can be described using conventional δ notation: where ( 56 Fe 54 Fe ) 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 isotopespecific pools ( 56 C i , for example, in the case of Fe). The isotopic signature of an Fe species C i is then calculated as Each individual reaction R k is assigned a fractionation factor 56 R k (‰; Table 4), which relates to 56 α R k as Isotope fractionation is then implemented by calculating an isotope-specific reaction rate (for the 56 Fe pool) from the bulk reaction rate R k : where  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 δ 56 Fe and δ 34 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 (Pel-lerin 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.

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  (1964) 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.

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 Fe 3+ free ), ligand-bound iron (Fe 3+ ligand ), and free ligand (L free ) -for our new tracer scheme of Fe 3+ (total = free + ligand-bound), Fe 2+ , and TDL (total dissolved ligand). Note that in the original scheme only two tracers were actually needed -total dissolved iron (Fe 3+ ) and total dissolved ligand (L) -as the equilibrium between Fe 3+ free , Fe 3+ ligand , and L free was recalculated each time step. Hence, our new tracer scheme is composed of these same two primary tracers (Fe 3+ , L), plus Fe 2+ . We then enable the two-way oxidationreduction reactions: Fe 2+ → Fe 3+ and Fe 3+ → Fe 2+ (Eqs. 13 and 18). We retain the same scavenging scheme as before, with iron scavenged by POM as Fe 3+ 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 2 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 2+ , we compare total dissolved Fe (TDFe), which is the sum of all dissolved Fe phases (i.e. [TDFe] = [Fe 3+ ] + [Fe 2+ ]). Adding Fe 2+ 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 2+ is assumed not to be subject to scavenging. Hence, reducing Fe 3+ to Fe 2+ 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 3+ (as FeOOH) and release of Fe 2+ (which again escapes scavenging). (Fig. 6c, f). 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 2 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 3+ 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.

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 2 concentration of 0.1 PAL (present atmospheric level, i.e. 21 000 ppm) and an atmospheric CO 2 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 2 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 2+ and SO 2− 4 and exit the ocean as FeOOH, FeS 2 , FeCO 3 , or Fe 3 Si 2 O 5 (OH) 4 . Fluxes of Fe and S are chosen to be in balance with respect to FeS 2 burial (S : Fe = 2 : 1) and to represent the best estimate of present-day weathering fluxes to the ocean (excluding reprocessing in inner-shelf sediment settings) (F Fe 2+ = 1.3 × 10 12 mol yr −1 , F SO 2− 4 = 2.6 × 10 12 mol yr −1 ; 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 SO 2− 4 (Kump and Seyfried, 2005). Our simulations therefore also include a hydrothermal flux of Fe 2+ , broadly comparable to a plausible input flux to Proterozoic oceans (15.1 × 10 12 mol yr −1 ; Thompson et al., 2019). This Fe 2+ 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 SO 2− 4 of 30.2 × 10 12 mol yr −1 . 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  FeOOH burial mol Fe yr −1 1.814 × 10 9 1.814 × 10 9 1.815 × 10 9 FeS 2 burial mol Fe yr −1 n/a n/a 0.0 n/a: not applicable. (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 3+ are spec-ified according to the scavenging scheme (Eq. 11). We focus here on (i) three depth slices of particle sinking fluxes of POC, FeOOH, FeS 2 , and Fe 3 Si 2 O 5 (OH) 4 ( Fig. 7; FeCO 3 fluxes were negligible and are not shown); (ii) the average and three vertical longitudinal slices of O 2 , dissolved Fe 3+ , Fe 2+ , and dissolved H 2 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 mmol C m −2 d −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 3+ ) penetration into the ocean interior (Fig. 8a-h), sinking particles have more time to continue to scavenge and accumulate Fe 3+ , and the flux increases with depth (Fig. 7e,  f). The maximum FeOOH flux at the seafloor reaches ∼ 300 µmol m −2 d −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 2 reach their maximum near the seafloor (where a source of deep Fe 2+ is supplied via the imposed hydrothermal flux), with maximum fluxes similar to those of FeOOH (Fig. 7g-i). The spatial pattern of Fe 3 Si 2 O 5 (OH) 4 formation is very similar to FeS 2 , but the fluxes are several orders of magnitude lower (Fig. 7jl). Fluxes of FeCO 3 were near zero everywhere (data not shown), consistent with recent work suggesting that watercolumn precipitation of FeCO 3 is difficult to achieve, even in iron-dominated oceans . The negligible FeCO 3 fluxes are at odds with Fespeciation data for Precambrian rocks showing that an important fraction of sedimentary Fe consists of reduced nonsulfurized 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).
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 3+ to a few kilometres of depth in our benchmark simulation (Fig. 8a-h). The concentrations of Fe 3+ 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 2+ ) that allows Fe 3+ 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 3+ 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 2+ and H 2 S ( Fig. 8l and p). Dissolved Fe 2+ reaches higher concentrations in the deeper ocean, largely as a result of deep hydrothermal inputs (Fig. 8i-l), whereas the highest H 2 S concentrations are spatially constrained to areas of more intense POC degradation along the Equator, just below the oxic zone ( Fig. 8m-p).
Because pervasive anoxia is not present in modern openocean 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.
Our model predicts a globally integrated FIO rate of ∼ 0.86 × 10 −7 mol Fe kg −1 yr −1 (Table 6), which is lower than rates estimated for anoxic lakes (12-51 × 10 −6 mol Fe kg −1 yr −1 ; Crowe et al., 2008;Walter et al., 2014). However, rates measured near the oxycline (5 × 10 −6 mol Fe kg −1 yr −1 ) are comparable to rates measured in anoxic lake systems. Interestingly, more than half of the H 2 S is re-oxidized via SMI (Table 6), which indicates that Fe is able to act as a relatively efficient intermediate between H 2 S oxidation and O 2 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 3 , 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 2 burial ( Fig. 10a) that drives δ 56 Fe − Fe 2+ to heavier values. Note that the very heavy δ 56 Fe − Fe 3+ values only occur in the ocean interior, where the concentration of Fe 3+ 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 2 ) 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 (SO 2− 4 is ∼ 5 ‰ heavier than its input value; Canfield and Farquhar, 2009). Additionally, free sulfide is isotopically lighter than sulfate (∼ −5 ‰ on a global scale; up to ∼ −15 ‰ locally), while buried FeS 2 expresses an isotope signature of ∼ −2.38 ‰ (Fig. 10b). Our baseline simulation thus suggests that the expressed isotope fractionation between SO 2− 4 and FeS 2 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.

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;, 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): where x j i represents the j th value of the ith parameter, j i 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). 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 (5.34 × 10 −3 M −0.5 h −1 ; Poulton et al., 2004), and the upper bound is taken to be 5 orders of magnitude higher (5.34 ×  Pianosi et al., 2015) to select random starting points x j (j = 1, . . .; r) and parameter variations i . For more information on the sampling strategy we refer the interested reader to Campolongo et al. (2011).
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 2+ and dissolved H 2 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 reveals that AR, DSR, mean surface O 2 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 2+ and H 2 S are removed from local seawater more rapidly, which results in the following: 1. a decreased build-up of Fe 2+ and H 2 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 2 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 3+ 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(Rickard, , 2006, with the result that our baseline value for this parameter is relatively wellconstrained. 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 mag- 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). nitude) (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 2+ inventory and a decrease in the H 2 S inventory (Fig. 13f, g); more Fe 2+ leads to 3. more FIO and thus a higher surface Fe 3+ concentration and more re-oxidation of H 2 S with dissolved Fe 3+ (SMI d ; Fig. 13e, h, k).  Table 7. Panels (a)-(c) and (h)-(n) are in moles per year (mol yr −1 ), and panels (d)-(g) are in moles per kilogramme (mol kg −1 ).
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 3+ that is not scavenged represents a "colloidal" pool, with a reactivity similar to that of hydrous ferric oxide. When Fe 3+ 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 3+ 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.

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 de-scribed here, in the future we plan to extend the organicmatter-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 be-tween marine ecosystems, nutrient availability, and ocean redox conditions . Lastly, future work will seek to include other redox-sensitive proxies and bioessential 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.
Code availability. The code for the version of the "muffin" release of the cGENIE Earth system model used in this paper is tagged as v0.9.21 and is assigned a DOI: https://doi.org/10.5281/zenodo.4651390 (Ridgwell et al., 2021b). The code is hosted on GitHub and can be obtained by cloning https://github.com/derpycode/cgenie.muffin, changing the directory to cgenie.muffin, and then checking out the specific release: $ git checkout v0.9.21. Configuration files for the specific experiments presented in the paper can be found in the following directory.
genie-userconfigs/MS/vandeveldeetal.GMD.2021 Details of the experiments, plus the command line needed to run each one, are given in the readme.txt file in that directory. All other configuration files and boundary conditions are provided as part of the code release. A manual detailing code installation, basic model configuration, tutorials covering various aspects of model configuration and experimental design, plus results output and processing is assigned a DOI: https://doi.org/10.5281/zenodo.4651394 (Ridgwell et al., 2021a). The LaTeX source of the manual, along with a pre-built PDF format version, can be obtained by cloning the following.
Author contributions. SJvdV and AR developed the model with input from all authors. CTR provided geological context. SJV and DH ran the sensitivity tests and analysed the output. SJvdV wrote the paper with input from all authors.
Competing interests. The authors declare that they have no conflict of interest.
Financial support. This research has been supported by the Belgian American Educational Foundation (postdoctoral fellowship), the National Aeronautics and Space Administration (Nexus for Exoplanet System Science), the NASA Astrobiology Institute (NASA postdoctoral programme fellowship, administered by the Universities Space Research Association under contract with NASA), and the Simons Foundation (postdoctoral fellowship in marine microbial ecology, award 653829).
Review statement. This paper was edited by Andrew Yool and reviewed by three anonymous referees.