Anoxic iron and sulphur cycling in the cGENIE.muffin Earth system model (v0.9.16)

The coupled biogeochemical cycles of iron and sulphur 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 sulphur (’euxinic’) vs. 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 continued, episodically, to develop regions of euxinic or ferruginous conditions, 5 such as associated with past key intervals of organic carbon deposition (e.g. during the Cretaceous) as well as extinction events (e.g. at the Permian/Triassic boundary). A better understanding of the cycling of iron and sulphur 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. 10 Here, we extend the ’muffin’ release of the intermediate-complexity Earth system model cGENIE, to now include an anoxic iron and sulphur cycle, enabling the model to simulate ferruginous and euxinic redox states as well as the precipitation of reduced iron and sulphur minerals (pyrite, siderite, greenalite) and attendant iron and sulphur isotope signatures, which we describe in full. While we cannot make direct model comparison with present-day (oxic) ocean observations, we use an idealized ocean configuration to explore model sensitivity across a selection of key parameters. We also present the spatial patterns of 15 concentrations and δFe isotope signatures of both dissolved and solid-phase Fe species in an anoxic ocean as an example application. Our sensitivity analyses show how the first-order results of the model are relatively robust against the choice default kinetic parameters 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) 20 cycling to the iron cycle. 1 https://doi.org/10.5194/gmd-2020-312 Preprint. Discussion started: 14 November 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
The biogeochemical cycles of iron and sulphur 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 oxidised states, both sulphur (in the form of sulphate; SO 2− 4 ) and iron (in the form of iron (oxyhydr)oxides; F eOOH) 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, 5 stratified lakes) (Thamdrup, 2000;Crowe et al., 2008;Raiswell and Canfield, 2012). In metabolising organic matter, microbial reduction of SO 2− 4 and F eOOH produce reduced sulphide (H 2 S) and ferrous iron (F e 2+ ), respectively. When present at the same location, H 2 S and F e 2+ combine into iron monosulphides and eventually pyrite (F eS 2 ) (Rickard, 1997(Rickard, , 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 2 S vs. F e 2+ , the precipitation of F eS 2 can lead to an anoxic water body becoming either 10 iron-rich ('ferruginous') or sulphide-rich ('euxinic') (Canfield, 1998;Poulton and Canfield, 2011) -H 2 S:F e 2+ ratios greater than 2:1 (the S:F e ratio in F eS 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 sulphur would have dominated the marine redox landscape (Poulton and Canfield, 2011;Raiswell and Canfield, 2012). 15 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 F e 2+ to their oxidised counterparts (Kharecha et al., 2005;Canfield 20 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 sulphur. The ability to simulate in models the evolution of these cycles as well as geochemical distributions in the ocean 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 F e 2+ (Fig. 1), which 25 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 hypothesised to reflect a transition to oxygen-rich bottom waters that prevented the buildup 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 30 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 of oceanic sulphate concentrations, allowing sulphide (produced by microbial sulphate reduction) to remove reduced iron from solution by the precipitation of F eS 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 has 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;Raiswell and Canfield, 2012). After (van de Velde et al., in review). See text for details. 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 et al., 2005;Poulton and Canfield, 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 5 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 of a certain point in time. Furthermore, 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 F e 2+ or H 2 S) are still 10 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, 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 carbon and various solid-Earth derived tracers Colbourn et al., 2013;Adloff et al., 2020). Here, we use the current 'muffin' release that 5 encompasses a range of developments and/or additions in the representation of: temperature-dependent metabolic processes in the ocean (Crichton et al., 2020), 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 weatheringrelevant trace-metals and isotopes (Adloff et al., 2020).
The climate component in cGENIE -C-GOLDSTEIN -has a reduced physics (frictional geostrophic) 2-D energy-moisture 10 balance model of the atmosphere coupled to a 3-D ocean circulation model and a 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 employs a much reduced spatial and temporal resolution relative to most high-resolution ocean general circulation models. While this precludes exploration of very detailed spatial patterns, it 15 does provide a flexibility which is not available in more high-resolution models, and the relatively short run time of cGENIE (around 1 day 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 20 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. 25 For the purpose of this study -implementing and characterising the coupled cycling of iron and sulphur in an anoxic ocean -we adopt a deliberately idealised model configuration. We configure the ocean model on a 18x18 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 30 of land connecting North and South poles (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 idealised boundary conditions of zonally-averaged wind stress and speed, plus a zonally-averaged planetary albedo, all following Vervoort et al. (in review). The solar constant is set to modern (1368.0 W m −2 ).

Continental configuration and climatology
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 manuscript), with the exception of the parameterisation controlling ocean mixing. When using the standard implementation of isoneutral diffusion and eddy-induced advection (Edwards and Marsh, 2005;Marsh 5 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 parameterisation, reverting to the original un-adjusted 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 parameterisations, and found only a minimal difference in the large scale ocean circulation (i.e., a ∼ 2 Sv stronger Atlantic meridional overturning circulation in

The biological carbon pump
In this paper, we adopt a representation of biological export from the surface ocean driven by an implicit (i.e. unresolved) biological community with a highly parameterised 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) as characterised by Tagliabue et al. (2016). The governing equations are summarised below. 5 Biological productivity in the euphotic zone (taken to be the surface layer in the ocean model) is controlled by dissolved phosphate (P O 3− 4 ) and dissolved iron (F e 3+ ) availability, the fractional ice coverage of each grid cell (A), mean ambient light and temperature. With the exception of the parameterisation of the temperature term, the equation for photosynthetic nutrient uptake follows Doney et al. (2006): Rates of photosynthetic nutrient uptake are scaled to ambient dissolved nutrient concentrations (P O 3− 4 and F e 3+ ), according to an optimal uptake timescale (τ bio =63.4 days; Meyer et al., 2016), and converting [F e 3+ ] to the equivalent dissolved phosphate concentration via the F e : P Redfield ratio (red P/F e ). The various limitation terms 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 20 for each limiting component ( 1 nM ) and are as used in Tagliabue et al. (2016). The influence of temperature on biological export is parameterised as: where k T 0 (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 two change per temperature change of 10 • C (Reinhard et al., 2020a). 25 A proportion (ν = 0.66) of P O 3− 4 taken up by biota is partitioned in dissolved organic phosphorus (DOP) while the remainder -as particulate organic phosphorus (POP) -is exported vertically out of the surface ocean. Because the biological configuration used here does not resolve explicit standing plankton biomass, the export flux of POP is always equal to the rate 7 https://doi.org /10.5194/gmd-2020-312 Preprint. Discussion started: 14 November 2020 c Author(s) 2020. CC BY 4.0 License. of P O 3− 4 uptake: where ρ is the density of seawater and h e the thickness of the euphotic zone (=80.84 m). The particulate organic carbon (POC) export flux is calculated using a fixed Redfield ratio as: After export from the surface grid cells, POC is remineralised instantaneously throughout the water column following a Martin -type curve (Martin et al., 1987), with a specified decay constant b (=0.7) ( Fig. 3a): 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 leads to a better reconstruction of the distribution of CaCO 3 in deep-sea sediments (Wilson, pers. comm.).

Oceanic iron and sulphur cycling
The motivation behind this paper to provide a tool that can aid understanding of the key interactions the key interactions between the biological carbon pump, the oxygenation of the atmosphere and ocean, and the marine biogeochemical cycles of iron (Fe) and sulphur (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 only consider a relatively limited number of potential 15 redox states. In this section we briefly describe the general conceptual cycle of Fe and S, as illustrated in Fig. 3b, and in the following subsections discuss the assumptions, reaction equations, and parameters for each cycle. The model equations that describe the biogeochemical reactions (summarised 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 mineralisation of POC, which is produced in the photic zone by 20 photosynthesis, and subsequently sinks through the water column . The sinking flux of POC follows a 'Martin'-type decay (Fig. 3a), which determines the mineralisation rate (R min ) at each depth layer in the water column, and thus R min at depth z is calculated as, Note that mineralisation only occurs below the surface layers; z > 1 in Fig. 2. The mineralisation of POC is coupled to the 25 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 mineralisation scheme includes aerobic respiration (AR, R1 in Table 1), dissimilatory iron reduction (DIR, R2 in Table 1), dissimilatory sulphate 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 (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 relationship with respect to TEA concentration 5 that allows for a non-linear closure of the system: where K 0 are the half-saturation constants for the four primary redox reactions and K i are 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 10 all of the POC mineralisation. When oxygen starts to become depleted, DIR, and subsequently DSR, become the dominant mineralisation pathways. When both iron oxides and sulphate are exhausted, MG represents the final mineralisation pathway.
While the occurrence and parametrisation 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 section 3.1.2. It should be noted that a ('Martin'- Table 1. List of biogeochemical reactions included in reduced Fe-S scheme within the cGENIE model. A distinction is made between mineralisation ("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 here any details for reactions involving the methane cycle, though they are included in the simulations described here. For more information about the parametrisation 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

R1
Aerobic respiration AR CH2O + O2 → CO2 + H2O 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, in any one depth interval in the ocean, exactly the same proportion of organic matter will be degraded. We avoid the alternative here -a fully kinetic set of equations where 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 to implement such a scheme effectively, requires knowledge about the composition of 5 settling organic matter and how its relative reactivity changes with time (Ridgwell, 2011;LaRowe and Van Cappellen, 2011).

R2
As a consequence of organic matter remineralisation in the model, DIR produces ferrous iron (F e 2+ ), which is re-oxidised 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 sulphide (because we do not consider sulphide speciation explicitly, reduced sulphide is defined as ΣH 2 S = H 2 S+HS − ) that comes into contact with oxygen is 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 (see text for details). The values of the kinetic constants are listed in Table 3.

R11
Greenalite precipitation k greenalite e (b greenalite log 10 (SI greenalite )) 1 Two different kinetic constants are used for the reactions between F e 3+ and sulphide (k SM I,d ) and F eOOH and sulphide (k SM I,s ), see text for details. 2 The formation of iron oxide minerals is parametrised according to Parekh et al. (2004), and is explained in more detail in section 3.1.1.
Pyrite, siderite and greenalite are subsequently buried in the sediment, together with any solid FeOOH that has not reacted (the 10 half-life of iron oxides in euxinic waters is in the order of 10s -100s of days, which is comparable to the residence time of a particle in the ocean; Poulton et al., 2004a) (Fig. 3b). 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).

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 the today's ocean volume is oxic; Diaz and Rosenberg, 2008), nor does it contain a coupling between the iron and sulphur cycles (e.g., via the precipitation of F eS 2 ). The 5 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 sulphur cycling would have dominated ocean biogeochemistry (Poulton and Canfield, 2011;Raiswell and Canfield, 2012). In addition to summarising the existing oxic cycle, we expand the model to include key processes operating in anoxic environments, as follows.
3.1.1 The oxic iron cycle 10 In oxygenated waters, iron is predominately present in its oxidised form (F e 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; total dissolved F e 3+ total consists of free F e 3+ f ree and ligand-bound F e 3+ ligand , and only the free fraction can be scavenged (and thus form F eOOH), or taken up by biology to fuel primary productivity (see section 2.2). The scavenging rate of iron is a function of the concentration of free F e 3+ 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 h −1 ) is a scavenging constant calibrated to the modern day distribution of F e (Tagliabue et al., 2016). We assume that the complexation reaction is always in equilibrium, which allows for the calculation of the amount of F e 3+ f ree at each time step by the conservation equations where K F eL sp is the stability constant of the L − F e 3+ complex (1.0 × 10 11 M −1 ; Table 3), and L total is the total amount of ligand.
Oxidised iron is very insoluble in seawater and will rapidly form particulate oxides (F eOOH), driving down equilibrium dissolved F e 3+ abundance to picomolar values without stabilisation by ligands (Liu and Millero, 2002). Thus, F e 3+ f ree will 25 form colloidal or nanoparticulate iron oxides -F eOOH -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 F e 3+ f ree can be thought of as a mix of purely dissolved F e 3+ and a range of colloidal and nanoparticulate F e 3+ phases that do not settle efficiently through the water column. The F e 3+ f ree 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, F e 3+ f ree (that is not stabilised by ligands) is scavenged by settling POC (using a fixed scavenging rate; Eq. 11) and 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 5 Gaillard, 2002;Tagliabue et al., 2016). Marine snow formation is what likely happens today, where 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). When scavenged on POC, F eOOH can be used for DIR when oxygen becomes depleted (see Section 3.1.2).
Note that we also provide the option in the model for F e 3+ f ree to precipitate as a pure F eOOH phase, without being associated 10 with POC. In this case, dissimilatory iron reduction does not occur in the water column since both POC and F eOOH 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', where all F e 3+ f ree that passes the solubility threshold of 0.5 nM (Liu and Millero, 2002) is considered particulate and sinks through the water column. 15 To close the oceanic iron cycle, we have introduced a new dissolved iron species -F e 2+ -any reduced F e 2+ that is mixed into the oxic zone is being rapidly oxidise to F e 3+ total via ferrous iron oxidation (FIO), where the rate equation can be expressed as 20 and where k F IO is a reaction rate constant (k F IO = 0.115 × 10 6 M −1 h −1 ) ( Table 2; Millero et al., 1987a).
In summary, we have extended the basic representation of Parekh et al. (2004), that 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 F e 2+ and F e 3+ total ), with F e 3+ total being split into free (and scavengeable) and ligand-bound phase.

The anoxic iron cycle 25
In an anoxic water column, oxidised forms of iron can be reduced to F e 2+ either via dissimilatory iron reduction or via sulphur-mediated iron reduction. We describe the two processes individually as follows.
Dissimilatory iron reduction in the water column 30 The majority of oxidised iron exists in particulate form due to the low solubility of F e 3+ , which poses a challenge for microorganisms performing dissimilatory iron reduction in the water column. Iron reducers need physical contact with the iron mineral 13 https://doi.org/10.5194/gmd-2020-312 Preprint. Discussion started: 14 November 2020 c Author(s) 2020. CC BY 4.0 License.
to reduce iron (Gorby, 2006), which is difficult when both organic matter and F eOOH are in particulate form and are dispersed separately throughout the aqueous medium. Even in a sediment column, where all particles are packed closely together, DIR generally requires sediment homogenisation by burrowing fauna to become volumetrically important (Thamdrup, 2000;van de Velde and Meysman, 2016). Indeed, in a sulphate-rich, ancient marine brine, F e 3+ was found to be the terminal electron acceptor, but sulphur ultimately acted as a redox shuttle between organic matter and iron oxides (Mikucki et al., 2009). Studies 5 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 mineralisation 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 10 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. 15 To model DIR in the water column, we use the mathematical expression of DIR in marine sediments, where F eOOH is a common electron acceptor for organic matter mineralisation (Thamdrup, 2000). Dissimilatory iron reduction has an overall reaction stoichiometry of, where it is implicitly assumed that every iron oxide particle consists solely of F e 3+ , rather than a mixture of redox states 20 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 µmol cm −3 (Van Cappellen and Wang, 1996;Thamdrup, 2000). This limitation is expressed as in Eq. 10 and follows the conventional limitation-inhibition scheme (Soetaert et al., 1996), where the parameter K 0,F eOOH expresses the concentration at which DIR occurs at half of its maximum rate. In marine sediments, DIR 25 generally occurs before sulphate 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 in Eq. 10, where K i,F eOOH expresses the concentration above which FeOOH inhibits other mineralisation pathways (DSR and MG). In early diagenetic models, this concentration is generally assumed to be identical to the limitation parameter K 0,F eOOH (Van Cappellen and Wang, 1996;Soetaert et al., 1996;30 van de Velde and Meysman, 2016). As a baseline value, we assume that K 0,F eOOH and K i,F eOOH both equal 10 −3 mol kg −1 (∼ 1 µmol cm −3 , 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,F eOOH and K i,F eOOH in our model sensitivity testing (section 4.3).

Sulphur-mediated iron reduction
Oxidised iron in the ocean can also be reduced via sulphur-mediated iron reduction (SMI), which follows the stoichiometry: where the reaction rate can be expressed as (Poulton et al., 2004a): Note that [F e 3+ ] in Eq. 17 can represent dissolved F e 3+ total or solid F eOOH. We assume that the dissolved form is highly reactive with sulphide (a reaction time of ∼ 5 minutes, which is comparable to the reactivity of freshly precipitated hydrous 10 ferric oxide) and has a reaction rate constant of k SM I,d = 2.64 × 10 2 M −0.5 h −1 (Poulton et al., 2004a). The solid form of F eOOH can represent a number of different oxidised iron minerals which are reactive towards sulphide on timescales ranging from hours to hundreds of days (Poulton et al., 2004a). For our baseline simulations, we assume that all F eOOH precipitates as lepidocrocite, which has a reactivity constant of k SM I,s = 1.98 × 10 0 M −0.5 h −1 (Poulton et al., 2004a). Lepidocrocite is less crystalline than the other (non-hydrous) iron oxides and precipitates when the rate of F e 3+ supply is low relative to the rate 15 of precipitation, as is generally the case in natural systems (Crosby et al., 1983). We discuss the model sensitivity to choices of k SM I,s in Section 4.3. The reduction of F e 3+ produces F e 2+ , which can either be re-oxidised when it comes in contact with O 2 (see section 3.1.1) or can form reduced minerals.

Reduced iron mineral formation
Reduced F e 2+ can form complexes with a number of ligands that are common in seawater, including Cl − , SO 2− 4 and bicar-20 bonate (Millero et al., 1995). Under anoxic conditions, however, the most important ligand is free sulphide (Rickard, 2006), which is produced by sulphate reduction. Together, dissolved F e 2+ f ree and the aqueous iron-sulphide complex (F eS aq ) make up ∼ 100 % of the total dissolved iron pool in sulphidic-anoxic seawater (Fig. 4a). The thermodynamic equilibrium can be calculated as 25 which was obtained from the visualMINTEQ database (Gustafsson, 2019), and is based on stability constants calculated by Luther et al. (1996).
Since dissolved F e 2+ f ree and the F eS aq complex are the dominant dissolved forms of reduced ferrous iron in natural waters (in anoxic seawater devoid of sulphide, dissolved F e 2+ f ree still represents ∼80 % of the dissolved F e 2+ total pool; Fig. 4b), we choose to only consider those two species. We do not explicitly model the F eS 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 F e 3+ ; see section 3.1.1). This is important, since F eS aq forms particulate F eS p , which is the precursor of pyrite (F eS 2 ) (Rickard, 1997(Rickard, , 2006. Furthermore, siderite or greenalite can only precipitate from F e 2+ f ree (Tucker, 2001;Tosca et al., 2015), and their precipitation rates are consequently dependent on the dissolved F e 2+ f ree concentration. Before each precipitation reaction, we calculate the concentration of F eS aq , F e 2+ f ree and H 2 S f ree using equation 18. We 5 then compare the F eS aq concentration to the solubility threshold (10 −5.7 M ;Rickard, 2006). If F eS aq surpasses the solubility value, some fraction of it precipitates instantaneously as particulate iron monosulfide (F eS p ) (Suits and Wilkin, 1998). Note that here we consider F eS 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 F eS p and free ΣH 2 S, with the production of free hydrogen gas (as S 2− in ΣH 2 S and F eS p has to be oxidised to S 1− in F eS 2 ; Rickard, 1997), Since hydrogen-gas is highly reactive it is almost instantaneously oxidised, with the reduction of an electron acceptor, like O 2 , SO 2− 4 or CO 2 . Given that pyrite precipitates under sulphidic conditions, oxygen is absent, and sulphate is the most likely electron acceptor  (Rickard, 1997). The rate equation for F eS 2 precipitation can then be written as Any F e 2+ that is not complexed with sulphide (F e 2+ f ree ) can form siderite (F eCO 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 F eCO 3 (OH) 1/2 . Since the thermodynamic data required to calculate the mineral solubility product of AFC is 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),
Another reduced iron mineral potentially important for anoxic, Fe-rich environments is the iron-silicate mineral greenalite (F e 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 20 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) SI greenalite = log 10 IAP K greenalite sp (28) 25 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 sulphur cycle
The oxidised form of dissolved sulphur in the ocean is sulphate (SO 2− 4 ). Under anaerobic conditions, sulphate is used as terminal electron acceptor during organic matter remineralisation: The produced sulphide can then be re-oxidised when it comes into contact with oxygen via canonical sulphide oxidation (CSO) The rate of CSO in cGENIE is dependent on the concentration of the electron donor (ΣH 2 S) and acceptor (O 2 ) and a second order rate constant k CSO = 0.625 × 10 6 M −2 h −1 (Zhang and Millero, 1993): 10 Alternatively, ΣH 2 S can be re-oxidised with F e 3+ (either F e 3+ total or solid F eOOH) during sulphide-mediated iron reduction (as described in reaction 16).
In cGENIE, the eventual sink for sulphur is the precipitation of pyrite (reaction 19) (ignoring sulphide 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.   3) -a sink likely less important on a globally integrated basis during Precambrian time, or during any other period in which 15 ocean [SO 2− 4 ] was relatively low (Grotzinger and Kasting, 1993;Crowe et al., 2014;Fakhraee et al., 2019). Planned future developments to cGENIE will incorporate an explicit gypsum cycle.

Isotope geochemistry
A particularly important application 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 sulphur (Raiswell and Canfield, 2012). As noted 20 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 (F e for example has four stable forms 54 F e, 56 F e, 57 F e and 58 F e) can potentially be used to track physicochemical processes that act to partition stable isotopes according to thermodynamic or 25 kinetic principles. For F e the most abundant isotopes are 54 F e and 56 F e, and deviations in the 56 F e/ 54 F e ratio of Fe-bearing aqueous and mineral phases from that of a reference material can be described using conventional δ-notation:   (Thamdrup, 2000),  (Johnson , 1982) where ( 56 F e 54 F e ) 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. https://doi.org/10.5194/gmd-2020-312 Preprint. Discussion started: 14 November 2020 c Author(s) 2020. CC BY 4.0 License.
To model the isotopic signatures of F e and S, we track the concentrations of the 'bulk' pools (C i ) and the isotope-specific pools ( 56 C i for example, in the case of F e). The isotopic signature of an F e species C i is then calculated as Each individual reaction R k is assigned a fractionation factor 56 R k (in ‰; Table 4), which relates to 56 α R k as Fractionations are then implemented by calculating the rate for the 56 F e pool from the bulk reaction rate R k In this way, we assign a fractionation factor to each of the reactions considered in our model, and are able to track the isotopic 10 signature of each F e and S species. This allows us to simulate the δ 56 F e 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 sulphur reducers (or iron oxidisers) express different fractionation factors, even though the overall reaction remains the same (Gomes and Johnston, 2017;Pellerin et al., 2019). This 15 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/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 20 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 end-member hypotheses for ocean chemistry.
4 Model testing 25 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  (Kaplan and Rittenberg, 1964 (Wiesli et al., 2004) 1 The isotopic fractionation for iron oxide precipitation is driven by the oxidation reaction (FIO). 2 The isotope fraction factor for greenalite precipitation is currently unknown.
or euxinic conditions is essentially absent. Only highly restricted basins such as the Cariaco basin or the Black and Baltic Seas 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. Hence, instead of calibrating our model results on a large observational dataset from the modern oceans, we will show the spatial concentration and isotope features of a hypothetical anoxic ocean (section 4.2), and -where possible 5 -compare our predicted reaction rates to rates obtained from ancient ocean analogues. Subsequently, we broadly illustrate the sensitivity of the model output to the newly introduced parameters of the iron-sulphur cycling (section 4.3). As discussed below, our model is largely robust across a range of values for key parameters, and predicts reaction and process rates that are comparable to those obtained experimentally or observed in modern analogue environments.

Model configuration 10
For testing and characterizing the Fe-S cycle model developments, we configure cGENIE with a single continent that runs from pole to pole ( Fig. 2c; Fig. 5a). Each model experiment is initialised 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 interpretation of the results.
The model is run in 'closed' configuration for all elements (notably C and P ), in which the ocean-atmosphere inventory for each element is always conserved. An exception is made for F e and S, which enter the surface ocean in the form of F e 2+ and SO 2− 4 , and exit the ocean as F eOOH, F eS 2 , F eCO 3 or F e 3 Si 2 O 5 (OH) 4 . Fluxes of F e and S are chosen to be in balance with respect to F eS 2 burial (S:Fe=2:1), and to represent the best estimate of present-day weathering fluxes to the ocean (excluding 5 reprocessing in inner shelf sediment settings) (F F e 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 F e 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 F e 2+ , broadly comparable to a plausible input flux to Proterozoic oceans (15.1 × 10 12 mol yr −1 ; 10 Thompson et al., 2019). This F e 2+ flux in our model setup 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 initialise the model with a 15 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 run the ocean circulation to steady state for 20,000 years. From that spin-up, we then run each individual experiment for 40,000 years and present model output from the 40,000th year of integration. It is important to emphasise that the idealised 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 20 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. 5. Even though we have a symmetric continent, the overturning patterns are skewed to the southern hemisphere, with a strong anticlockwise circulation at around -60°N (Fig. 5c). In contrast, the barotropic streamfunction shows a large degree of symmetry ( Fig. 5d). Since the aim of this manuscript is to describe the newly developed Fe-S chemistry, and not ocean circulation in real 25 or fake worlds, we will not further discuss the emerging patterns of ocean circulation.

Spatial output
In the following section, we briefly discuss the spatial model output and modelled redox cycling for our baseline configuration.
Precipitation fluxes of oxidised F e 3+ are specified according to the scavenging scheme (Eq. 11). We focus here on: (i) three depth slices of particle sinking fluxes of P OC, F eOOH, F eS 2 , and F e 3 Si 2 O 5 (OH) 4 ( Fig. 6; F eCO 3 fluxes were negligible 30 and are not shown); (ii) the average and three vertical longitudinal slices of O 2 , dissolved F e 3+ total , F e 2+ total and dissolved ΣH 2 S (Fig. 7); and (iii) an overview of globally averaged reaction rates and iron mineral burial fluxes (Table 5).

Bulk fluxes and concentrations
The P OC 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. 6a-c). The spatial pattern of solid F eOOH flux (following the scavenging scheme described earlier -see section 3.1.1) matches the P OC flux pattern immediately below the surface layer, with higher values at the equator and poles (light brown shading in Fig. 6d). Most of this F eOOH is reduced in anoxic subsurface layers and 5 the flux declines with depth. However, at the poles where deep convection allows for greater oxygen (and F e 3+ ) penetration into the ocean interior (Fig. 7a-h), sinking particles have more time to continue to scavenge and accumulate F e 3+ and the flux increases with depth (Fig. 6e,f). The maximum F eOOH flux at the seafloor reaches ∼ 300 µmol m −2 d −1 (Fig. 6f), comparable to that observed in typical modern ocean margin sediments (see e.g.; van de Velde and Meysman, 2016). Fluxes of F eS 2 reach their maximum near the seafloor (where a source of deep F e 2+ is supplied via the imposed hydrothermal 10 flux), with maximum fluxes similar to those of F eOOH (Fig. 6g-i). The spatial pattern of F e 3 Si 2 O 5 (OH) 4 formation is very similar to F eS 2 , but the fluxes are several orders of magnitude lower (Fig. 6j-l). Fluxes of F eCO 3 where near-zero everywhere (data not shown), consistent with recent work suggesting that water column precipitation of F eCO 3 is difficult to achieve, even in iron-dominated oceans Tosca et al., 2019).
We additionally find that even in an idealised ocean with a simple symmetrical continental configuration, complex spatial 15 patterns emerge in the Fe-S redox chemistry (Fig. 7). For instance, the poles are more well-ventilated, allowing oxygenated  conditions and persistence of F e 3+ to a few kilometres depth in our benchmark simulation (Fig. 7a-h). The concentrations of F e 3+ are much higher than the nano -picomolar concentrations we observe in the ocean today (Tagliabue et al., 2016) -a consequence of much higher rates of iron delivery to the surface ocean (through upwelling of reduced F e 2+ ), which allows F e 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 compare well to oxic water layers overlying anoxic 5 deep water (Taillefert and Gaillard, 2002;Crowe et al., 2008). Note that this F e 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 F e 2+ and ΣH 2 S (Fig. 7l,p).
Dissolved F e 2+ total reaches higher concentrations in the deeper ocean, largely as a result of deep hydrothermal inputs (Fig. 7i-l), whereas the highest ΣH 2 S concentrations are spatially constrained to to areas of more intense P OC degradation along the 10 equator, just below the oxic zone (Fig. 7m-p).
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 mineralisation (Table   5), consistent with previous observations from ferruginous Lake Matano and the euxinic Black Sea (Konovalov et al., 2006;5 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 6-7 orders of magnitude less important than all other mineralisation pathways (Table 5). A 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., ∼ 300m in Lake Matano versus 5000m in our idealized ocean), there is a much stronger coupling between sedimentary processes (i.e., the recycling of Fe as a benthic flux) 10 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).

Isotope patterns
Figure 8a-e shows the modelled stable Fe isotope patterns for all key Fe-bearing dissolved and solid-phase species (with the exception of F eCO 3 , which is a negligible component in our benchmark simulation). In our model simulations, all dissolved 5 F e that enters the ocean is assigned an isotope signature of 0.0 ‰. This allows us to observe the isotope fractionation of all F e phases relative to the F e that entered the ocean (Fig. 9a). The dissolved iron phases show similar isotopic signatures (∼ 1.1 ‰; Fig. 8a,b; Fig. 9a), which can be explained by the large amount of isotopically light F eS 2 burial (Fig. 9a), which drives δ 56 F e−F e 2+ to heavier values. The isotope signature of buried oxidised iron (F eOOH) is around 1.8 ‰ heavier than the F e that entered the ocean, whereas the major burial fraction of reduced iron (F eS 2 ) has an isotope signature that is ∼ 1 ‰ lighter, 10 which reflects their relative importance as an Fe burial phase (Fig. 8c,e; Fig. 9a). These isotope values are broadly comparable to phase-specific stable Fe isotope observations from ancient sedimentary rocks (Heard and Dauphas, 2020). Similarly, Fig.   8f-g shows 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 ‰. Sulphate is isotopically enriched relative to its input value (∼ 3.3 ‰; Fig. 8f; Fig. 9b), 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 sulphide is isotopically lighter than sulphate (∼ -5 ‰ on a global scale; up to ∼ -15 ‰ locally), while buried F eS 2 expresses an isotope signature of ∼ -3.25 ‰ (Fig. 9b). Our baseline simulation thus suggests that the expressed isotope fractionation between SO 2− 4 and F eS 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 5 strong potential for tracking iron and sulphur isotope signatures for comparison with Earth's rock record.

Sensitivity analysis
We evaluate model output sensitivity to four key parameters of our modelled iron-sulphur cycle (K 0,F eOOH , K i,F eOOH , k SM I,s , k P yP ) using the Elementary Effect Test (EET; Morris, 1991). These parameters were chosen because they are either unconstrained by laboratory experiments (K 0,F eOOH , K i,F eOOH ; Section 3.1.2), represent a complex mixture of different iron 10 minerals with different reactivities (k SM I,s ; Poulton et al., 2004a) or are expected to have a strong and potentially difficult to forecast influence on other reactions (k P yP ; van de Velde et al., in review). Other parameters introduced in the model are either relatively well-constrained by laboratory studies and calibrated on field data (k CSO , k F IO ; Millero et al., 1987a, b;Ridgwell Figure 9. Isotope mass balance for the (a) 56 F e and (b) 34 S systems from our baseline simulation. Purple arrows are influxes, black arrows outfluxes. Bubbles represent the inventory of the whole ocean. Note that in our simulation, the S system is not in steady-state (because F eS2 burial is too low). 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-sulphur cycle (k AF C , b AF C , k greenalite , b greenalite ).
The EET method estimates global sensitivity by calculating the mean of r finite differences ('Elementary Effects') (Pianosi et al., 2016), 5 where x j i represents the j th value of the i th parameter, ∆ j i the variation on the i th parameter, f () is the model output for a given set of parameters and c i is a scaling factor. A higher mean value 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. 10; Pianosi et al., 2016). 10 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 6. We vary the limitation and inhibition constants K 0,F eOOH and K i,F eOOH 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 sulphide (k SM I,s ) is defined as the reactivity of hematite (5.34 10 −3 M −0.5 h −1 ; Poulton et al., 2004a), and the upper bound is taken 15 to be 5 orders of magnitude higher (5.34 10 2 M −0.5 h −1 , which is comparable to freshly precipitated hydrous ferric oxide). We test the pyrite precipitation constant (k P yp ) across a range between 0.3708 10 −3 M −1 h −1 and 0.3708 10 2 M −1 h −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., 2020). Our sensitivity ensemble consists of one hundred individual model experiments, using Latin-Hypercube sampling approach (using the SAFE toolbox; 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 5 refer the interested reader to Campolongo et al. (2011).
The EET analysis suggests that changes in K 0,F eOOH can significantly impact DIR, whereas the inhibition constant K i,F eOOH has a relatively minor impact on model output (Fig. 10). This is expected, as both parameters act only on the mineralisation 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 sulphide (Mikucki et al., 2009). Even though both K 0,F eOOH 10 and K i,F eOOH are the least constrained by experimental results (see section 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 P yP and k SM I,s both exerted more notable impact on model output across a range of diagnostics (Fig. 10). In particular, k P yP had a significant impact across all model output diagnostics analysed here both in terms of model sensitivity and interactivity with other parameters (Fig. 10).
Because pyrite precipitation controls the inventories of both dissolved F e 2+ and dissolved ΣH 2 S, reducing or increasing the 15 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 where we changed a parameter from its baseline value to its lower and upper bound, whilst keeping 20 the other 3 parameters at their baseline values (Table 6). Figure 11 reveals that AR, DSR, mean surface O 2 concentrations and CSO are relatively insensitive to parameter choices (Fig. 11a,c,d,i). Dissimilatory iron reduction is only important when K 0,F eOOH is set at its lowest value, and even then it is three orders of magnitude smaller than AR or SR (Fig. 11a-c). Consistent with the EET analysis (Fig. 10), K 0,F eOOH and K i,F eOOH have no influence on the model output, aside from the magnitude of DIR, which in itself is of less importance 25 as mineralisation pathway ( Fig. 11b and Section 4.2.1). In contrast, changes in the pyrite precipitation constant k P yP impact Figure 10. Sensitivity analysis of the four key-parameters of the iron-sulphur 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, points that plot higher indicate that the parameter is more interactive with other parameters. Data processing was done with the SAFER toolbox of Pianosi et al. (2015).
several model outputs. When pyrite precipitation rates are elevated, F e 2+ and ΣH 2 S are removed from local seawater more rapidly, which results in: 1. a decreased build-up of F e 2+ and ΣH 2 S in the ocean interior (Fig. 11f,g), leading to  Table 6. Panels (a)-(c) and (h)-(n) are in mol yr −1 , panels (d)-(g) are in mol kg −1 .
2. a decrease in both aerobic and anaerobic re-oxidation pathways (FIO, SMI d and SMI s ; Fig. 11h-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. 11a,c) and 4. a decrease in F e 3+ concentrations in the surface ocean (Fig. 11e), and 5 5. less burial of F eOOH, increasing the burial of reduced iron minerals (Fig. 11l-n) The effect on greenalite burial is non-linear (Fig. 11o), which indicates that at higher values of k P yP , pyrite precipitation is competing with greenalite precipitation. Overall, our sensitivity analysis suggests that k P yP 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 well constrained. The fourth parameter, which also influences the model output, is the reactivity parameter of solid iron oxides (k SM I,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, 1992;Poulton et al., 2004a). Therefore, we explore the sensitivity of this parameter in more detail using a range of measured reactivity constants by Poulton et al. (2004a), whilst keeping all other parameters at their baseline values (Fig. 12). Increasing the reactivity of particulate F eOOH 5 1. increases the anaerobic re-oxidation reaction of sulphide with F eOOH (SMI s ) at the expense of the aerobic re-oxidation reaction (CSO) (Fig. 12i,j), which then leads to 2. an increase in the F e 2+ total inventory and a decrease in the ΣH 2 S inventory (Fig. 12f,g), more F e 2+ total leads to 3. more FIO, and thus a higher surface F e 3+ concentration, and more re-oxidation of ΣH 2 S with dissolved F e 3+ (SMI d ; Fig. 12e,h,k) 10 4. Because of the higher reactivity of the F eOOH particles, less F eOOH is buried, and more reduced F e-minerals are buried ( Fig. 12l-o) Although it is clear that changing k SM I,s impacts model output, the overall magnitude of the effect is moderate when compared to changing k P yP . Nevertheless, the choice of k SM I,s is critical. We choose a baseline reactivity value (Table 3) comparable to lepidocrocite for several reasons. Firstly, we assume all F e 3+ that is not scavenged represents a 'colloidal' pool, 15 with a reactivity of similar to that of hydrous ferric oxide. When F e 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 F eOOH 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 section 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, 20 field evidence suggests that F eOOH 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 F eOOH minerals that would resist reduction passing through a sulphidic 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, 25 these observations support the presumption that once F e 3+ is scavenged, its reactivity is less than hydrous ferric oxide but higher than that of goethite, motivating our default k SM I,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 30 to include coupled, anoxic Fe and S biogeochemical cycles. Because direct tuning with actual measured concentrations and . Lines indicate experimental reactivity parameters as presented in (Poulton et al., 2004a). Panels (a)-(c) and (h)-(n) are in mol yr −1 , panels (d)-(g) are in mol kg −1 . Note that the x-axis is logarithmic.
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 parametrisation 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 5 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 10 cycle (see e.g., Dale et al., 2015). In particular, the model does not include a representation of benthic iron reduction and 33 https://doi.org/10.5194/gmd-2020-312 Preprint. Discussion started: 14 November 2020 c Author(s) 2020. CC BY 4.0 License.
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 will 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 5 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, and in particular the scavenging and remobilisation 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 10 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 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 to validate it against observations from modern anoxic systems and the geologic record.

15
Code availability. The specific version used of the cGENIE.muffin model used in this paper is tagged as release v0.9.16, and has been assigned a DOI:10.5281/zenodo.4033262 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.16

20
Configuration files for the specific experiments presented in the paper can be found in the directory: genie-userconfigs/MS/vandeveldeetal.GMD.2020 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 release. A manual, detailing code installation, basic model configuration, plus an extensive series of tutorials covering various aspects of muffin capability, experimental design, and results output and 25 processing, is provided on GitHub. The latex source of the manual, along with a pre-built PDF format version, can be obtained by cloning: https://github.com/derpycode/muffindoc Author contributions. SJV and AR developed the model with input from all authors. CR provided geological context. SJV and DH ran the sensitivity tests and analysed the output. SJV wrote the paper with input from all authors.
Competing interests. The authors declare that they have no conflict of interest.