Constraining Global Marine Iron Sources and Ligand-Mediated Scavenging Fluxes With GEOTRACES Dissolved Iron Measurements in an Ocean Biogeochemical Model

Iron is a key micronutrient controlling phytoplankton growth in vast regions of the global ocean. Despite its importance, uncertainties high regarding external iron source fluxes and internal cycling on a global scale. In this study, we used a global dissolved iron data set, including GEOTRACES measurements, to constrain source and scavenging fluxes in the marine iron component of a global ocean biogeochemical model. Our model simulations tested three key uncertainties: source inputs of atmospheric soluble iron deposition (varying from 1.4 to 3.4 Gmol/yr), reductive sedimentary iron release (14–117 Gmol/yr), and compared a variable ligand parameterization to a constant distribution. In each simulation, scavenging rates were tuned to reproduce the observed global mean iron inventory for consistency. The variable ligand parameterization improved the global model-data misfit the most, suggesting that heterotrophic bacteria are an important source of ligands to the ocean. Model simulations containing high source fluxes of atmospheric soluble iron deposition (3.4 Gmol/ yr) and reductive sedimentary iron release (114 Gmol/yr) further improved the model most notably in the surface ocean. High scavenging rates were then required to maintain the iron inventory resulting in relatively short surface and global ocean residence times of 0.83 and 7.5 years, respectively. The model simulates a tight spatial coupling between source inputs and scavenging rates, which may be too strong due to underrepresented ligands near source inputs, contributing to large uncertainties when constraining individual fluxes with dissolved iron concentrations. Model biases remain high and are discussed to help improve global marine iron cycle models. ,

Another key aspect of marine iron models is the representation of ligands that organically bind DFe and thereby prevent it from being scavenged to sinking particulates. Some models still prescribe a globally constant ligand concentration typically at 1 nM, while others account for ligand distributions via a parameterization or directly simulating ligands as a prognostic tracer. Ligands are thought to be produced by microbes as a by-product during the production of organic matter (Gledhill & Buck, 2012), including by heterotrophic siderophores that flourish when systems become iron stressed (Bundy et al., 2018). This has led modelers to predict ligand concentrations by assuming they are produced during the production of organic matter (e.g., Völker & Tagliabue, 2015) or by prescribing a relationship to other organic tracers such as dissolved organic matter (DOM) and apparent oxygen utilization (AOU) (e.g., Misumi et al., 2013;Pham & Ito, 2018;Tagliabue & Völker, 2011).
The uncertainties associated with external source fluxes and scavenging represent key gaps in understanding the global marine iron cycle. This hampers accurate estimates of the DFe budget, residence time and, consequently, its sensitivity to environmental perturbations and climate change. While the rapidly increasing amount of DFe measurements is improving our knowledge of the distribution and inventory of dissolved iron in the ocean, constraining external fluxes has proved to be more difficult. As a result, the range of residence times estimated by the current global marine iron cycle models ranges from less than a decade to multiple centuries (Tagliabue et al., 2016), which limits our ability to confidently predict the impact of changes to the marine iron cycle on productivity in a future ocean. Observational estimates fall within a similar range (Johnson et al., 1997), noting that more recent studies estimate much shorter residence times in the upper ocean (∼10 days-4 years) (Croot et al., 2004;Sarthou et al., 2003) depending on the local dynamics, iron pools considered, and source inputs in different regions (Black et al., 2020).
In this study, we use a global marine DFe data set to constrain the iron cycle fluxes in a global marine biogeochemical model. We analyze model sensitivity simulations that focus on three key uncertainties: varying source fluxes of (a) atmospheric soluble iron deposition and (b) reductive sedimentary iron release, as well as the role of a (c) variable ligand distribution on DFe distribution and scavenging rates. The resulting DFe concentrations in each model simulation are evaluated against observations to determine the most realistic marine iron cycle fluxes among the model scenarios.

Model Description
We used the UVic Earth System Climate Model (Weaver et al., 2001) version 2.9 (Eby et al., 2009). In the following section, we provide a general overview of the model components then focus on improvements made to the marine iron cycle in this study, whereas other modifications applied to all model simulations are described in the supplementary information.

Physical Model
The physical ocean-atmosphere-sea ice model includes a three-dimensional (1.8° × 3.6°, 19 vertical levels) general circulation model of the ocean (Modular Ocean Model 2) with parameterizations such as diffusive mixing along and across isopycnals and eddy-induced tracer advection (Gent & McWilliams, 1990). The physical configuration is based on Somes et al. (2017) and includes parameterizations such as computation of tidally induced diapycnal mixing over rough topography on the sub-grid scale (Schmittner & Egbert, 2014), anisotropic viscosity (Large et al., 2001;Somes et al., 2010), and enhanced zonal isopycnal mixing schemes in the tropics to better represent zonal equatorial undercurrents (Getzlaff & Dietze, 2013). A two-dimensional, single level energy-moisture balance atmosphere and a dynamic-thermodynamic sea ice model are used, forced with prescribed monthly climatological winds (Kalnay et al., 1996) and constant ice sheets (Peltier, 2004).

Marine Biogeochemical Model
The updated marine ecosystem-biogeochemical model coupled within the ocean circulation model is based on the Model of Ocean Biogeochemistry and Isotopes ( Figure S1). It combines latest features from previous studies focusing on the nitrogen cycle (Somes & Oschlies, 2015), iron cycle , and carbon chemistry (Kvale et al., 2015), and is also constrained by isotope systems of 13 C and 15 N (Schmittner & Somes, 2016) (not shown here). Our model experiments were simulated for over 5,000 years under pre-industrial boundary conditions as they approached their quasi steady state.

Base Configuration
The marine iron model configuration is based on the previous UVic Kiel Marine Biogeochemistry Model (KMBM) , including improvements implemented in Muglia et al. (2017) (Figure 1). The marine iron model includes explicit tracers for DFe and particulate iron (PFe). All phytoplankton grow with a constant elemental stoichiometry ratio of iron relative to nitrogen. The sources of DFe to the ocean are atmospheric soluble deposition (Luo et al., 2008), reductive dissolution and release from sediments (Elrod et al., 2004;Moore & Braucher, 2008), and hydrothermal fluxes (Tagliabue et al., 2010) (Table 2, Figure 2). The ligand concentration determines the fraction of DFe that is organically complexed and thus unavailable for scavenging, whereas the remaining free DFe (DFe') pool can be scavenged to PFe, which then sinks and remineralizes at the same rate as POM (Table S1). In the base simulation #1, ligands are prescribed to be globally constant at 1 nM as in previous iterations of the model. This simulation is given the name SrcLow_LigCon to reflect its differences (i.e., low source inputs of atmospheric soluble deposition and reductive sedimentary iron release, and constant ligand distribution) from further changes made to the marine iron model in this study (see subsections below and Tables 1 and 2).

Scavenging
The formulation for scavenging and partitioning of free and organically complexed DFe is based on from previous model parameterizations Galbraith et al., 2010). Scavenging of DFe' to PFe occurs via two mechanisms in the model: (a) absorption onto POM following (Honeyman et al., 1988;Parekh et al., 2004) which depends only on DFe' and the inorganic scavenging rate constant (kFe prp ) following the scheme of Galbraith et al. (2010). This inorganic scavenging term primarily represents colloidal aggregation into larger, sinking particles as well as lithogenic scavenging not explicitly accounted for in our model. Here we use a non-linear formulation for inorganic scavenging following Galbraith et al. (2010) which was designed to account for high lithogenic scavenging rates to better reproduce DFe where atmospheric deposition is high (e.g., tropical and subtropical North Atlantic) (Pham & Ito, 2019;Ye & Völker, 2017). Note that we included a slightly higher non-linear exponent (2.) compared to Galbraith et al., 2010 (1.5) that better reproduced DFe in high atmospheric deposition areas in our model. This difference may be related to the fact that Galbraith et al., 2010 model included higher phytoplankton iron quotas when DFe is high which further reduces DFe in that model, whereas our model formulation assumes constant iron stoichiometry due to high uncertainties associated with this process. Thus, our model performed better with higher scavenging rates to reduce the overestimation of DFe in these high deposition areas.
In each model simulation, the scavenging rate constants (kFe org , kFe prp ) were manually tuned so that each simulation contains a nearly identical global iron inventory with an average global DFe concentration of 0.7 ± 0.03 nM ( Table 2). The inorganic scavenging rate constant was adjusted until the model reproduced the mean observed DFe concentration in the ocean interior since it is the dominant form of scavenging there, whereas the POM scavenging rate constant was adjusted to reproduce declining DFe concentrations toward the surface ocean ( Figure 4). The globally integrated rates of the different scavenging processes are shown in Table 2, vertically integrated rates from high and low source input simulations in Figure 2, and total basin-scale averages in Figure 4.

Ligand Parameterization
In the base model configuration, a constant ligand concentration of 1 nM is applied globally, and thus has LigCon in its model name (see Table 1). However, the distribution of ligands in the real ocean is variable (e.g., Völker & Tagliabue, 2015). Since iron-binding ligands are thought to be produced during the production of organic matter (Gledhill & Buck, 2012), which might explain why DOM and AOU may qualitatively reflect some observed ligand concentration patterns (Misumi et al., 2013;Pham & Ito, 2018;Tagliabue & Völker, 2011). However, a first global model-data comparison with ligands simulated as prognostic tracers found ligand distributions difficult to constrain with available observations and is further complicated by large variations in binding strength of different types of ligands (Völker & Tagliabue, 2015). Therefore, to maintain computational efficiency, we pragmatically chose to implement ligand concentrations as a function of existing tracers rather than include additional prognostic tracers.
SOMES ET AL.
where α (0.015 nmol ligand/(mmol O 2 m −3 ) 0.8 ) and β (0.21 nmol ligand/(mmol DON m −3 ) 0.8 ) are generic parameters that determine ligand concentration associated with the tracers AOU and DON, respectively. The parameters α and β were chosen so that the global ligand mean concentration remained at 1 nM, consistent with simulation #1 with constant ligands, but now reflects changes in their spatial distribution ( Figure 3). Model simulations with this variable ligand parameterization (simulations #2-5, see Table 1) have LigVar in their respective model simulation name.
Although we follow previous studies for the variable ligand parameterization (Misumi et al., 2013;Pham & Ito, 2018;Tagliabue & Völker, 2011), a few notable changes have been made in our version. Since AOU can be negative in the surface ocean due to dissolved oxygen supersaturation, we applied a minimum ligand concentration of 0.5 nM. Previous ligand parameterizations have also applied minimum ligand concentrations to account for ligands associated with more refractory forms of DOM not explicitly included in our model (Aumont et al., 2015;Tagliabue & Völker, 2011). We also applied an exponential parameter (0.8) to the AOU and DON terms, which reduces ligands associated to these tracers particularly when their concentrations are high. This helped the model from overestimating DFe concentrations when AOU and DON concentrations are at their highest concentrations in the model.

Reductive Sedimentary Iron Release Parameterization
The base model version uses reductive sedimentary iron release based on the Moore and Braucher (2008) implementation of Elrod et al. (2004), where the Fe flux from the sediments (Fe sed ) is determined by the sedimentary iron release rate ( FeSed = 0.27 μmol Fe/mmol C ox m −2 d −1 ), and organic carbon oxidation (C ox ) in the sediments. The base model version uses the DFe flux rate from Nickelsen et al. (2015) that is lower than suggested by Elrod et al. (2004) (0.72 μmol Fe mmol C ox −1 m −2 d −1 ). Since this formulation yields lower global rates of this source input in the model compared with other implemented sedimentary functions included in this study (described below), model simulations with this sedimentary iron release implementation (#1-2) contain the name SrcLow, noting they also include a low source input of atmospheric soluble iron deposition (see Section 2.3.5 below).
We also implemented the sedimentary iron release function proposed by Dale et al. (2015), who compiled a global data set of sedimentary DFe fluxes to constrain their model estimate. While it has a strong depend-SOMES ET AL.  Black et al. (2020) by including the upper 250 m and account for sinking particulate iron out of this layer as the sink flux. Since our particulate iron pool includes both biogenic (i.e., produced during primary production) and authigenic (i.e., produced by scavenging) iron in the model, this model residence time is comparable to their mean dissolved, biogenic + authigenic estimate, which ranges from 0.1 to 4 years depending on location.
where γ FeSedMax is the maximum flux under steady-state conditions, and bwO 2 is dissolved oxygen concentration in bottom waters interacting with the sediments.
We test two scenarios with the Dale et al.  Tables 1 and 2). This reduced value was chosen to test a global sedimentary DFe flux approximately halfway in between SedHigh and SedLow since their fluxes differ by a large amount. Note that the SedMid simulation does not produce a significantly different spatial distribution compared to SedHigh.

Atmospheric Soluble Iron Deposition
We applied the atmospheric soluble iron deposition mask from Luo et al. (2008)    deposition models that explicitly accounted for the soluble iron deposition rather than assuming a constant solubility from total deposition.
Another estimate we test in this study applies the average flux from four recent atmospheric soluble iron deposition models (Myriokefalitakis et al., 2018). The intermodel average global soluble deposition rate is 3.4 Gmol yr −1 with similar patterns to Luo et al. (2008) but higher rates most notably in the North Atlantic. This simulation with high atmospheric soluble iron deposition (AtmHigh; Figure 2) is applied to the simulation with high sedimentary release and variable ligands and is therefore named Atm + SedHigh_LigVar.

Global Dissolved Iron Data Set
The DFe database used in this study is a collection of observations from both GEOTRACES Intermediate Data Product 2017 (7,520 points; Schlitzer et al., 2018) and prior observations compiled by Tagliabue et al. (2012) (12,371 points). Note that we excluded 37 measurements (19 from GEOTRACES, 18 from prior) with high DFe concentrations between 10 and 216 nM mainly from locations with high hydrothermal activities, but also some near-shore settings (e.g., Laptev Sea, Bristol Bay, Peruvian coastal waters near urban area of Trujillo) and around small islands not resolved in the model (e.g., Kerguelen, Indonesian and Coronation), and thus the data set used here contains concentrations up to 10 nM. We then interpolated the data onto the UVic model grid using the PyFerret SCAT2GRIDGAUSS function developed by NOAA's Pacific Marine Environmental Laboratory, which is a Gaussian interpolation function based on Kessler and McCreary (1993). This gridded data was used for the model-data comparison (Figures 3-7) and to calculate model-data statistical metrics (i.e., correlation coefficient, (uncorrected) standard deviation, and root-mean-squared error) (Figure 8). It covers 5,917 grid points since many observations overlap and thus are averaged on corresponding grid points. Since we compare to annual model results, we interpolated all observations onto the grid and thus temporal aspects and variability of the data are not taken into account or investigated in this study.
Model-data misfit statistical metrics are sensitive to unresolved outlier concentrations and spatial extent of the data interpolation onto the model grid. However, these aspects do not affect which simulations best reproduce the global data set according to statistical metrics. This is illustrated by comparing metrics calculated from all observations (triangles) to only GEOTRACES (circles) in Figure 8. The statistical metrics slightly improve when comparing against only GEOTRACES observations, with the only exception being root-mean-squared error for model simulation #1 in the surface ocean, but the relative improvements in the model simulations are nearly identical. The arbitrary exclusion concentration threshold of 10 nM was chosen as a balance between including as many observations as possible while still being able to calculate useful statistical metrics that are not dominated by these outlier concentrations.

Variable Ligand Distribution
The simulation with constant ligands does not reproduce the major basin-scale features of the observed DFe distribution, despite that its globally averaged depth profile is generally consistent with observations ( Figure 4c). Most notably, simulations with constant ligands significantly overestimate the DFe in the interior Southern Ocean (Figure 4o), a critical ocean basin for Fe-limited phytoplankton growth. LigCon thus overestimates supply of DFe via upwelling, and underestimates Fe limitation of phytoplankton growth, which is a key deficiency in the base configuration and previous model versions (e.g., Muglia et al. (2017)). They also underestimate DFe in intermediate waters in the Indian and Pacific Ocean (Figures 4k and 5b), which we have averaged together since they have similar deep ocean biogeochemical tracer profiles relative to the global average ( Figure S1).
The simulations with variable ligand concentrations (#2-5; LigVar) better reproduce the ocean interior distribution of DFe ( Figure 5). This is primarily due to the AOU dependence of the variable ligand parameterization that mainly determines  Southern and Indian-Pacific Oceans, which contain relatively low and high values of AOU and thus ligand concentrations, respectively, according to our parameterization (see Figures 3 and S1). Lower ligand concentrations in the Southern Ocean enhances scavenging causing lower DFe concentrations, with the opposite effect occurring in the Indian-Pacific Ocean, resulting in better reproduced observations in both basins. Therefore, the interior DFe distribution with the variable ligand parameterization is better partitioned with respect to observations (Figures 4 and 5) and improves the global model-data misfit by 9.2% when averaging across our three metrics (i.e., correlation coefficient (R), normalized standard deviation (nSTD), and SOMES ET AL.
The concentration of semi-refractory DON largely determines ligand concentrations in the surface ocean ( Figure 3a). DON concentrations are higher around the high productivity regimes in the low latitudes with generally decreasing values toward higher latitudes (Somes & Oschlies, 2015) ( Figure S2). This pattern is reflected in the surface DFe distribution that shows the same latitudinal trend in the variable ligand model (Figures 6c-6d). While this meridional DFe pattern better reproduces low DFe concentrations in the open Southern Ocean, it creates larger model-data biases on high latitude continental shelves in the Bering Sea, Weddell Sea, and European shelf seas (Figures 6a-6d, 7c and 7e). This shows that while the overall variable ligand effect significantly improves the global DFe distribution (Figure 7), model-data biases in some regions (e.g., high latitude continental shelf seas) still increase, which contributes to a smaller average metric improvement (3.9%) in the surface layer compared to the global ocean.

Sedimentary Iron Release
The simulations with low sedimentary source inputs (#1-2 SrcLow) provide a relatively poor fit to observed DFe concentrations according to the statistical metrics ( Figure 8). They fail to reproduce the high DFe concentrations near continental margins (Figures 6 and 7), suggesting higher sedimentary release rates are necessary to explain these features. The simulated DFe distribution also lacks the strong spatial gradient toward depleted concentrations in many open ocean regions in the observations. These overly smooth gradients in SOMES ET AL.

10.1029/2021GB006948
13 of 20 SrcLow are the result of low sedimentary release rates and subsequent low scavenging rates that are then required to reproduce the global mean DFe inventory, resulting in a relatively long global mean residence time of 35 years among our simulations ( Table 2).
The simulations with higher sedimentary release rates (Figure 2e) produce higher DFe concentrations in continental shelf seas (Figures 6 and 7), particularly where bottom water oxygen is low in the low latitudes. The simulations applying high-end sedimentary Fe release rates (SedHigh) modestly outperformed simulations assuming lower rates across all calculated statistical metrics (Figure 8) on average by 3.3% in the global ocean and slightly higher by 3.8% in the surface layer, with the intermediate release rate scenario SedMid performed between SedLow and SedHigh. Therefore, our model-data analysis suggests that high-end estimates for global reductive sedimentary iron release rates are the most realistic.
One region that was notably improved by high sedimentary release rates was the low latitude margins near oxygen deficient zones (ODZs) (Figures 6 and 7). Observations there in both the eastern tropical South Pacific off Peru (Figure 7a), eastern tropical South Atlantic off Namibia (Figure 7d), and northern Indian Ocean show high DFe concentrations that are best reproduced in SedHigh scenarios. Since SedHigh simulations also contain high scavenging rates, they better reproduce the lowest DFe concentrations in the offshore open ocean locations as well.
The high DFe concentrations on high latitude continental shelf systems (Figures 6, 7c and 7e) are not improved in SedHigh_LigVar due to the interactions with ligands and scavenging. Decreasing surface ligand concentrations toward high latitude systems ( Figure 3) allow scavenging to compensate the additional sediment-derived DFe more efficiently, in contrast to low latitude systems near ODZs (e.g., Tropical Pacific) that contain higher ligands allowing DFe to be retained in the water column. This causes the simulation with constant ligands to retain slightly higher DFe compared to simulations with variable ligands in high latitude continental shelf systems (e.g., Bering Sea ( Figure 7c) and European Shelf Seas (Figure 7e)), despite that these simulations with variable ligands include much higher sedimentary release rates there (e.g., Se-dHigh_LigVar, Figure 2). This demonstrates that more efficient scavenging rates associated with low ligands can overcompensate the high sedimentary release rates in determining DFe concentrations in the model.

Atmospheric Soluble Deposition
The two soluble atmospheric deposition scenarios tested here predict similar spatial depositional patterns ( Figure 2), with the more recent GESAMP intermodel average (Myriokefalitakis et al., 2018) providing a significantly higher global deposition rate (3.4 Gmol yr −1 ) relative to the low estimate from Luo et al. (2008) (1.4 Gmol yr −1 ). These enhanced rates cause higher DFe concentrations mainly from the Saharan dust plume in subtropical North Atlantic, but also to a lesser degree in the Arabian Sea and North Pacific (Figures 6g, 6h, 7c and 7f). The impact of including higher soluble deposition only slightly improves the global model-data statistical metrics by 0.7% globally and 1.5% in the surface layer, making it difficult to determine the most realistic rates based on our model-data DFe comparison alone.

High Scavenging Effect
In model simulations with high source fluxes (e.g., #5 Atm + SedHigh_LigVar), higher scavenging rates are necessary to maintain a realistic global DFe inventory (Tables 1 and 2, Figures 2h-2i and 3). Scavenging is thus more efficient at reducing DFe concentrations in the high source flux simulations. In regions far away from the source fluxes, particularly in the deep ocean and open Southern Ocean (e.g., see Figure 6), the model simulations with higher source fluxes actually contain lower DFe because the enhanced scavenging outweighs the source fluxes in these areas (Figure 4). Lower DFe concentrations in these deep and open ocean regions better reproduce observations further improving the model-data misfit metrics (Figure 8). The combined effects of high atmospheric and sedimentary source inputs, which also includes highest scavenging rates, contributed to the largest improvement in the surface ocean across our metrics (5.5% improvement relative to SrcLow_LigVar).

Model-Data Constraints and Uncertainties
The variable ligand parameterization improved the model's ability to reproduce the global distribution of DFe observations the most. This is most evident in the interior ocean due to AOU dependency of this parameterization. Since ligands are produced when dissolved oxygen is consumed during the respiration of POM via heterotrophic microbes in the variable ligand parameterization, their concentrations reach maximum values in old Pacific intermediate waters (Figure 3). High ligands reduce scavenging that causes the model to better reproduce high observed DFe concentrations there (Figures 4k and 5c), a feature that has also been demonstrated in other models (Misumi et al., 2013;Pham & Ito, 2018;Frants et al., 2016). This model improvement suggests that ligand production by heterotrophic bacteria is a key mechanism maintaining the global marine iron cycle.
The model simulations that include higher source inputs and scavenging rates show a subtle but continuous improvement in the model-data misfit metrics particularly in the surface ocean (Figure 8). This is in contrast to the intermodel comparison study of Tagliabue et al. (2016), which showed no clear relationship between model performance and source inputs, as well as an inverse modeling study of Pasquier and Holzer (2017), which could not find an optimal solution among their large set of model simulations varying source inputs. However, Pasquier and Holzer (2017) only tested relatively low sedimentary release rates (up to 22 Gmol/yr compared to 117 Gmol/yr in this study) and also did not include an oxygen dependency that has a strong influence in our parameterization. Our analysis emphasizes that future modeling studies should test these important factors associated with reductive sedimentary DFe release that contributed to the model improvements in this study.
The ligand and high sedimentary DFe release effects have similar impacts on DFe spatial distributions making it difficult to constrain their individual impacts with DFe concentrations alone. This spatial overlap is most pronounced above ODZs in the eastern tropical Pacific, eastern tropical Atlantic, and Northern Indian Ocean (Figure 6). This spatial covariance occurs because when AOU is high, bottom water oxygen is typically low. Therefore, DFe concentrations are enhanced both by reduced scavenging due to high ligands where AOU is high, as well as by higher sedimentary DFe release rates where bottom water oxygen is low. Future studies should examine the integrative DFe cycling in these systems (e.g., sedimentary release and scavenging rates, ligand concentrations) to give additional insights on individual processes contributions to total DFe concentrations.
Despite high sedimentary release rates, the SedHigh model simulations still underestimate DFe on most continental shelf systems (Figure 7). The poorly resolved coastal dynamics in our coarse resolution circulation model is likely a key model deficiency preventing the model from representing many coastal dynamics where sedimentary DFe fluxes are high. Coarse resolution models underestimate coastal upwelling and the nutrient input on narrow shelf systems that drive productivity. This bias causes underestimated POM production as well as overestimated dissolved bottom water oxygen concentrations, both of which contribute to underestimating reductive sedimentary DFe release rates on coastal shelf systems.
Further complicating matters are interactions between sedimentary DFe release rates, ligands, and scavenging. For example, our SedHigh_LigVar model simulation releases significantly higher DFe on high latitude shelves (Figures 2e-2f). However, only a small part of this DFe remains in the dissolved pool since scavenging efficiently converts it to particulate iron that eventually sinks back to the sediments (Figures 2h-2i). Therefore, our model underestimation of DFe concentrations remains despite high DFe release rates. This strong spatial coupling between source and scavenging fluxes has also been demonstrated in other modeling studies Pasquier & Holzer, 2017), which also found that this tight spatial coupling significantly contributes to the difficulty in constraining source inputs. The exclusion of riverine inputs that may also directly include ligands could also contribute to overly efficient scavenging resulting in underestimated DFe. If our ligand parameterization predicted higher concentrations on these high latitude shelf systems, which has been indicated by ligand observations (Völker & Tagliabue, 2015), this would prevent rapid scavenging of DFe released from sediments and better reproduce observations.
Sedimentary DFe release rates may still be underestimated even in our high release scenario. Note that our highest tested global sedimentary release rate (117 Gmol yr −1 ) was not the highest from the marine iron model intercomparison (up to 194 Gmol yr −1 ) (Tagliabue et al., 2016), and every model scenario tested here with increased source fluxes improved the model-data misfit metrics (Figure 8). Potentially important sedimentary processes not included in the model are non-reductive dissolution and release from reactive sediments in tectonically active or volcanic regions (Conway & John, 2014;Homoky et al., 2013) and sedimentary colloidal production (Homoky et al., 2021), which could further contribute to higher total sedimentary DFe release rates that may improve the model-data misfit.
An important limitation of applying these empirical functions of reductive sedimentary DFe release (e.g., Dale et al., 2015;Elrod et al., 2004) in global models is that total iron balance within the sediments is not explicitly accounted for. Thus, these parameterizations can potentially represent an unlimited long-term supply of DFe to the ocean which is unrealistic. This simplification can be justified because many important sources of particulate Fe to the sediment are not yet included in the model, for example, atmospheric and riverine input of lithogenic material and in situ production at volcanic islands or active margins, which provide DFe for release. Also note that the Dale et al. (2015) parameterization applied in the SedHigh simulations sets a maximum rate determined under steady-state conditions which caps potentially unrealistic high release rates. While this simplification is likely not a significant deficiency in steady-state model simulations presented here, this should be considered in transient simulations with substantial enhancement of sedimentary DFe fluxes.
Atmospheric deposition often occurs at high rates over continental shelves (e.g., North Pacific, Patagonia) and ODZs (e.g., Arabian Sea), again making it difficult to constrain individual processes driving DFe concentrations when multiple processes act together in close spatial proximity. For example, our high atmospheric soluble deposition scenario helps reproduce high DFe concentrations in the Arabian Sea ( Figure 7f). However, our model underestimates the extent of the Arabian Sea ODZ which could be the real cause driving high DFe concentrations there via high sedimentary DFe release, reduced scavenging, and/or enhanced redox cycling (Moffett et al., 2007). Instead the model ODZ is mostly misplaced to the Bay of Bengal, where higher simulated DFe there in the model better matches observations within the real ODZ in the Arabian Sea (see star symbols in Figure 7f).
The model simulations do not resolve the high variance of the observations which is reflected in the underestimated standard deviation (Figures 4 and 8). This occurs everywhere in the ocean and is most pronounced in the Southern Ocean due to it containing very low DFe in the open ocean but also high concentrations near islands, continental margins, and hydrothermal vents (Figures 4-6). Although not a focus of this study, the model was not able to reproduce the full spatial extent of high DFe concentrations near hydrothermal vents at mid-ocean depths (Figures 4 and 5), despite that this source is included (Table 2). Previous modeling studies were only able to reproduce this high DFe extent when assuming that the hydrothermal vents were also a significant source of ligands Resing et al., 2015) or included stabilization via reversible scavenging (Roshan et al., 2020), both of which we have not accounted for in our model. This emphasizes that future model versions should include all important ligands and scavenging dynamics to better represent their importance in marine iron models, but that a more robust global database of ligand concentrations including their binding strength would be required (Völker & Tagliabue, 2015).
High variance in the global data set may not reflect mean climatological conditions simulated by the preindustrial steady-state model results given the highly dynamic nature of DFe cycling particularly in the surface ocean with short residence times (Black et al., 2020). The spatial and temporal sparsity of the data set likely contribute to high variance as well. But note that the standard deviation was significantly improved in our best model simulation with variable ligands and high source/scavenging fluxes (Atm + SedHigh_LigVar; see Figures 4, 8b and 8e) suggesting that a model with low residence times can better reproduce the high variance and strong gradients in the DFe observations. Since most DFe observations have been collected in recent decades, there could already be a significant anthropogenic impact (e.g., enhanced deoxygenation, atmospheric/riverine pollutants) on the global marine iron cycle not included in these model simulations, especially if the marine DFe residence time operates on decadal timescales or less. Future additions and expansion to the global DFe data set as well as comparison with transient model simulations at the same period of data collection will improve uncertainties in future model-data analyses.

A Global Marine Iron Cycle With a Residence Time Under a Decade?
Our model simulations testing various external source fluxes in the global marine iron cycle result in global average residence times ranging from 7.5 to 36 years. The simulation that best reproduces the observations (Atm + SedHigh_LigVar) has the lowest residence time (global: 7.5 years; surface ocean: 0.83 years) among our model experiments. This low-end residence time is caused in large part due to the high source fluxes, with the reductive sedimentary release being the most important with the highest global rate in our simulations. These high source fluxes need to be compensated by efficient scavenging and subsequent removal via burial in the sediments to reproduce the distribution and global mean inventory in DFe observations, a model feature that was also found in other modeling studies (e.g., see Frants et al., 2016;Pasquier & Holzer, 2017). This is in general agreement with observational studies focusing on the surface layer (Black et al., 2020;Sarthou et al., 2003). For example, Black et al. (2020) estimated similar residence times throughout the global surface ocean (0-250 m) for DFe ranging from approximately 1 month to 4 years depending on the region and specific iron pools considered, although noting that the uncertainties remain large (i.e., equal or greater than the absolute value of the estimate in each region). These generally low surface residence times are captured in our model simulations that range from 0.83 to 3.12 years (Table 2). However, residence times of individual molecules and regions can further vary depending on the local coupling of source inputs, scavenging efficiency, and regeneration (e.g., Holzer et al. (2016); Pasquier and Holzer (2018); Tagliabue et al., 2019). For instance, DFe in the ocean interior is more stable and controlled by the amount of ligands that reduces scavenging and removal to the sediments via sinking particulates, contributing to the longer global residence times.

Marine Iron Flux Impacts on Global Ocean Biogeochemistry
An interesting feature of the model simulations is that there is surprisingly little change to globally averaged marine productivity and export production (Table 3). This occurs in large part in the model because scavenging was also increased in high sedimentary iron release scenarios, and thus much of the additional DFe fluxes from the sediments is efficiently scavenged to particulate iron that sinks back to the sediments before it can be transported to the surface ocean where it may stimulate additional productivity. This general impact was also found in a model study using a previous iteration of the model version used here but comparing different complexities of the marine iron configurations (Yao et al., 2019) as well as other inverse modeling studies (Pasquier & Holzer, 2017. However, it must be noted that all of these model studies, including this one, only evaluated steady-state simulations in which uncertain parameters were manually tuned or optimized to best reproduce observations. Therefore, they are not necessarily indicative to how the iron dynamics in the model may respond to and impact marine productivity in externally forced transient scenarios. There is a notable decrease in marine productivity and export production in the Southern Ocean among our model simulations with better representations of the global iron distribution (Table 3). The variable ligand parameterization predicts less ligands in the Southern Ocean (Figure 3), which allows higher scavenging to SOMES ET AL.  reduce DFe that better reproduces observations. Furthermore, since external iron sources in the Southern Ocean are small (Figures 2 and 4m), the enhanced scavenging in the high source flux simulations removes more DFe than source fluxes add to the Southern Ocean. Therefore, DFe levels further decrease in the Southern Ocean (Figures 4o and 6) in the high source flux scenarios. The high scavenging in our best model simulation with variable ligands and high source fluxes (Atm + SedHigh_LigVar) reduces DFe, marine productivity and resulting oxygen consumption during remineralization of POM, thereby increasing dissolved oxygen concentrations at depth. This effect is significant enough to increase average global dissolved oxygen concentrations by 8% in the model because water masses formed in the Southern Ocean contribute to much of the global deep ocean (Table 3). This emphasizes the importance of simulating a robust global marine iron cycle most importantly in the Southern Ocean.

Conclusions
In this study we tested various rates of atmospheric soluble deposition, reductive sedimentary release, and variable ligand distributions within a marine iron component in a global ocean biogeochemical model. The simulations that best reproduce the global DFe observations include highest tested source fluxes and a variable ligand parameterization. The most striking feature in the global DFe observations that supports this hypothesis is the strong gradients that often occur with high concentrations near source fluxes and low concentrations in adjacent open ocean regions. This high source flux/scavenging iron cycling regime causes a relatively short residence times of less than a decade in the global oceans and less than a year in the surface ocean. The short residence time implies that the global marine iron cycle is highly sensitive to environmental perturbations in the Anthropocene and geological past. Uncertainties remain high due to model parameterizations of complex, poorly understood, and often intertwined processes (e.g., ligand production and subsequent control on scavenging near source inputs) and the sparsity of DFe and ligand measurements throughout the global ocean. Nevertheless, our model-data analysis suggests the marine iron cycle operates with high global source inputs and scavenging rates and low residence times compare to most previous estimates.