Inverse-model estimates of the ocean ’ s coupled phosphorus , silicon , and iron cycles

The ocean’s nutrient cycles are important for the carbon balance of the climate system and for shaping the ocean’s distribution of dissolved elements. Dissolved iron (dFe) is a key limiting micronutrient, but iron scavenging is observationally poorly constrained, leading to large uncertainties in the external sources of iron and hence in the state of the marine iron cycle. Here we build a steady-state model of the ocean’s coupled phosphorus, silicon, and iron cycles embedded in a dataassimilated steady-state global ocean circulation. The model includes the redissolution of scavenged iron, parameterization of subgrid topography, and small, large, and diatom phytoplankton functional classes. Phytoplankton concentrations are implicitly represented in the parameterization of biological nutrient utilization through an equilibrium logistic model. Our formulation thus has only three coupled nutrient tracers, the three-dimensional distributions of which are found using a Newton solver. The very efficient numerics allow us to use the model in inverse mode to objectively constrain many biogeochemical parameters by minimizing the mismatch between modeled and observed nutrient and phytoplankton concentrations. Iron source and sink parameters cannot jointly be optimized because of local compensation between regeneration, recycling, and scavenging. We therefore consider a family of possible state estimates corresponding to a wide range of external iron source strengths. All state estimates have a similar mismatch with the observed nutrient concentrations and very similar large-scale dFe distributions. However, the relative contributions of aeolian, sedimentary, and hydrothermal iron to the total dFe concentration differ widely depending on the sources. Both the magnitude and pattern of the phosphorus and opal exports are well constrained, with global values of 8.1 ± 0.3TmolPyr−1 (or, in carbon units, 10.3 ± 0.4 PgCyr−1) and 171.± 3.TmolSiyr−1. We diagnose the phosphorus and opal exports supported by aeolian, sedimentary, and hydrothermal iron. The geographic patterns of the export supported by each iron type are well constrained across the family of state estimates. Sedimentary-iron-supported export is important in shelf and large-scale upwelling regions, while hydrothermal iron contributes to export mostly in the Southern Ocean. The fraction of the global export supported by a given iron type varies systematically with its fractional contribution to the total iron source. Aeolian iron is most efficient in supporting export in the sense that its fractional contribution to export exceeds its fractional contribution to the total source. Per source-injected molecule, aeolian iron supports 3.1 ± 0.8 times more phosphorus export and 2.0 ± 0.5 times more opal export than the other iron types. Conversely, per injected molecule, sedimentary and hydrothermal iron support 2.3 ± 0.6 and 4. ± 2. times less phosphorus export, and 1.9± 0.5 and 2.± 1. times less opal export than the other iron types.


Introduction
The ocean's nutrient cycles control the primary productivity of the global marine ecosystem and the ocean's biological carbon pump, which are crucial components of the global carbon cycle that regulate atmospheric CO 2 concentrations.The nutrient cycling of the ocean is governed by the interplay Published by Copernicus Publications on behalf of the European Geosciences Union.
of the ocean's advective-diffusive circulation, biological utilization, biogenic particle transport, and the external sources and sinks of nutrients.The cycles of macro-and micronutrients are coupled through co-limitation on biological uptake and through the scavenging of micronutrients such as iron by sinking organic matter.
We focus on dissolved iron (dFe) as a key micronutrient because of its well-documented fundamental role in primary production (e.g., Boyd and Ellwood, 2010).Indeed, dFe was suggested to limit oceanic phytoplankton growth as early as the 1930s (e.g., Gran, 1931;Hart, 1934).Since then, numerous studies have reported that iron deficiency limits productivity over vast areas, particularly over the high-nutrient lowchlorophyll (HNLC) regions like the Southern Ocean (e.g., de Baar et al., 1995;Lundry et al., 1997;Martin and Fitzwater, 1988;Boyd et al., 2007;Boyd and Ellwood, 2010).Martin (1990) went as far as to suggest that perturbations in the iron cycle played a crucial role in past climate fluctuations.More recently, iron-enrichment field experiments (e.g., Boyd et al., 2007) and model simulations (e.g., Nickelsen and Oschlies, 2015) have demonstrated the importance of iron for the global biological pump.
We model phosphate (PO 4 ) because it is essential to the metabolism of all living organisms (e.g., Smith, 1984;Howarth, 1988), which allows all biological production to be keyed to phosphate utilization (e.g., Kwon and Primeau, 2008;Primeau et al., 2013;Holzer and Primeau, 2013).Silicic acid (Si(OH) 4 ) was considered because of the importance of diatoms in marine ecosystems (e.g., Nelson et al., 1995;Buesseler, 1998;Moore et al., 2004;Brzezinski et al., 2011b) and because the pronounced silicon trapping of the Southern Ocean (e.g., Holzer et al., 2014) might be sensitive to iron availability.
With a changing climate, we expect not only changes in the ocean circulation, but also changes in the winds, hydrological cycle, and land use, and hence in the aeolian iron supply.To understand how such changes impact the ocean's nutrient cycles, it is necessary to model the coupling between the nutrients mechanistically.While global biogeochemistry models have been used extensively for this purpose (e.g., Tagliabue et al., 2014;Nickelsen and Oschlies, 2015), none of these models have been objectively constrained by the available observations.Here, we formulate an inverse model of the ocean's coupled macronutrient and iron cycles embedded in a data-assimilated global circulation.The biogeochemical parameters of the model are determined by objectively minimizing the mismatch with observed nutrient and phytoplankton concentrations.To ensure the optimization problem remains tractable with a reasonable computational burden, we formulate a model of intermediate complexity.
The intercomparison of iron models by Tagliabue et al. (2016) showed that current models contain significant uncertainties.Despite the fact that the models have iron source strengths that range over nearly 2 orders of magnitude, all models can be tuned to roughly the same mean dFe concentration with an intermodel variance of only 27 %.This is due to essentially unconstrained scavenging rates so that models are free to employ different scavenging strengths to balance the sources at roughly comparable dFe concentrations.All the models of the intercomparison are prognostic forward models that are computationally too expensive for exploring the biogeochemical parameter space systematically or for computing the sensitivity with respect to multiple parameters (e.g., Kwon and Primeau, 2006).Here we aim to close this gap in our ability to constrain the iron cycle objectively by formulating a numerically highly efficient model of the iron cycle that is mechanistically coupled to key macronutrient cycles and that is embedded in a data-assimilated global circulation.
We build on the simple iron-only inverse model of Frants et al. (2016), for which biogeochemical parameters were optimized and which allowed novel diagnostics to be computed such as the mean iron age and rigorous source attribution of dFe (Holzer et al., 2016).Consistent with the findings of the iron-model intercomparison, Frants et al. (2016) showed that current dFe observations cannot constrain the iron sources because of local compensation between sources and sinks.Frants et al. (2016) therefore explored a family of state estimates corresponding to a range of aeolian source strengths, all of which are consistent with the currently available dFe observations.
However, in the study of Frants et al. (2016) the iron cycle could not interact with a prescribed phosphate cycle, and the silicon cycle was not considered at all.Moreover, for the family of state estimates of Frants et al. (2016) only the aeolian source was varied and it is unclear if the wide range of hydrothermal and sediment iron sources found in the literature is consistent with the dFe observations.The model of Frants et al. (2016) is therefore unsuitable for investigating the effects of iron-source perturbations on biological production or for exploring how much export is supported by the different iron sources.Here, we overcome these shortcomings by explicitly coupling the iron, phosphorus, and silicon cycles through their mutual co-limitations so that the macronutrients can respond to changes in dFe.We furthermore refine the modeling of the sedimentary iron source, improve the representation of iron scavenging to include redissolution, and model three phytoplankton functional classes with concentrations that are derived from a steady-state logistic equation (e.g., Dunne et al., 2005).Through these advances we are able to produce, for the first time, a family of data-constrained state estimates of the coupled Fe-P-Si cycles for a wide range of not only aeolian, but also hydrothermal and sedimentary sources.We find that these state estimates are roughly equally consistent with the observed macronutrient and dFe concentrations regardless of the iron source strengths.Analysis of our family of state estimates shows that the uncertainty in the iron sources stems not only from compensation between overlapping sources and sinks, but also from the ability of the different iron source types (aeolian, hydrothermal, and sedimentary) to compensate for each other despite their different spatial distributions.
We use our inverse-model estimates of the coupled Fe-P-Si cycles to address an important open question about the marine iron cycle: what are the relative contributions of the different iron sources to supporting the world ocean's export production?While there have been perturbation experiments with forward models (Tagliabue et al., 2009(Tagliabue et al., , 2010(Tagliabue et al., , 2014) ) where one type of source (e.g., hydrothermal or sedimentary) was shut down to assess the importance of dFe from the missing source, such experiments cannot quantify the true contribution of hydrothermal or sedimentary iron to biological production because of the nonlinearity of the iron cycle (Holzer et al., 2016).Moreover, these numerical experiments were conducted with definite but highly uncertain choices of the iron sources, and the models were not objectively constrained by the observations.Thus, in addition to presenting the first inverse model of the coupled Fe-P-Si cycles, we address the following key scientific questions: 1. How well can the modeled dFe, PO 4 , and Si(OH) 4 concentrations be fitted to observations for widely differing iron sources, and are there limits on the iron source strengths that are consistent with the observed dFe concentrations?
2. What are the nutrient limitation patterns that emerge from the data-constrained estimates of the coupled nutrient cycles, given that direct observational data on these patterns are very sparse?
3. How well constrained are the phosphorus and opal exports of optimized state estimates with widely different iron sources?
4. What fractions of the phosphorus and opal exports are supported by aeolian, hydrothermal, and sedimentary iron, and how do these fractions vary with the ironsource strengths?
In the following, we detail the model formulation in Sect. 2 and the optimization strategy in Sect.3. In Sect.4, we quantify the fidelity of the family of state estimates to the nutrient observations.We examine nutrient limitation in Sect.5, export production in Sect.6, and iron-attributed export in Sect.7. Caveats are discussed in Sect.8 and we summarize and conclude in Sect.9.

Biogeochemical model
We distinguish three phytoplankton functional groups, nondiatom small and large phytoplankton as well as diatoms, with a nominal separation between small and large at a cell diameter of 2 µm.We denote the molar PO 4 uptake rate per unit volume of each class by U c , where the subscript c ∈ {lrg, sml, dia} identifies functional class.The uptake rates U c are only nonzero in the model's upper 73.4 m (two layers), the model's euphotic zone.
We consider the three nutrients PO 4 , Si(OH) 4 , and dFe and denote their concentrations by χ i , with i ∈ {P, Si, Fe}.We write the steady-state tracer equations for these concentrations by keying all biological production to the uptake of phosphate U c as follows: T χ Fe = c (S Fe c − 1)R Fe : P U c (3) + (S s,POP − 1)J POP + (S s,bSi − 1)J bSi − J dst + s A + s S + s H .
In Eqs. ( 1)-(3), T is the advection-eddy-diffusion operator, the operators S i c model the biogenic transport and remineralization of nutrient i taken up by functional class c, and the operators S s,POP and S s,bSi model the particle transport of scavenged iron and its partial redissolution at depth as the scavenging particles remineralize or dissolve (details in Sect.2.2).The iron scavenging rates per unit volume are J POP for scavenging by particulate organic phosphorus (POP), J bSi for scavenging by opal particles, and J dst for scavenging by mineral dust (details in Sect.2.4.2).The terms s A , s S , and s H are the aeolian, sediment, and hydrothermal iron sources (details in Sect.2.4.1).The factors R Si : P and R Fe : P are the stoichiometric uptake ratios that allow us to key all production to phosphorus.These ratios are functions of the nutrient concentrations as described in Sect.2.3.3.
The terms proportional to γ g in Eqs.
(1)-(2) fix the global mean phosphate and silicic acid concentrations through weak relaxation to their observed global means χ obs P and χ obs Si .This is necessary because the phosphorus and silicon cycles have no external sources and sinks to set the global mean in steady state.(For phosphate and silicic acid, external sources, e.g., riverine input, and loss to sediment burial are neglected.)We choose the restoring timescale γ −1 g = 10 6 years ("geological" restoring); there is no sensitivity to the precise value of γ g .
Equations ( 1)-( 3) are coupled via the uptake of phosphate, U c , which depends on the concentrations of all three nutrients, via the iron scavenging, which depends on the export fluxes of organic matter and opal, and via the sedimentary release of dFe, which is keyed to the flux of organic matter onto the sediments (Elrod et al., 2004), as discussed in detail below.

Circulation
We use the data-assimilated, steady (nonseasonal) circulation of Primeau et al. (2013), which has a horizontal resolution of 2 The circulation is also constrained dynamically, and the data assimilation used the wind-stress climatology of Trenberth et al. (1989) and specified horizontal and vertical viscosities of 5 × 10 4 and 10 −4 m 2 s −1 , respectively.The circulation's advective-diffusive transport operator has fixed horizontal and vertical eddy diffusivities of 10 3 and 10 −5 m 2 s −1 , respectively.We emphasize that the circulation effectively provides a ventilation-weighted transport because it has been optimized against PO 4 and the ventilation tracers CFC-11 and radiocarbon.The steady circulation thus does not bias estimates of preformed nutrients in the way an annual-average circulation would.

Biogenic transport
Organic matter sinks as POP, dissolves, and remineralizes at depth.Inverse models of the phosphorus cycle (Primeau et al., 2013;Holzer and Primeau, 2013;Pasquier and Holzer, 2016) suggest that dissolved organic phosphorus (DOP) represents a relatively small fraction of the total dissolved phosphorus that we neglect here (no DOP tracer) for simplicity and numerical efficiency.We note that our estimates of phosphorus export effectively capture the export due to DOP, despite DOP not being explicitly represented as a separate tracer.This is because the optimization of the biogeochemical parameters minimizes the mismatch with the observed PO 4 distribution, which in the real ocean is determined by the remineralization of all organic phosphorus, including DOP.
Because the particle transport is much faster than the fluid transport across a grid box, we approximate particle transport and remineralization, which acts as an interior source of nutrients, as instantaneous.We model this process for each phytoplankton functional class by the "source" operator, S P c , which reassigns a "detrital" fraction f c of the uptake rate to a remineralization rate throughout the water column, while a fraction 1 − f c remineralizes in situ where the uptake occurred.We therefore express S P c in terms of a biogenic redistribution operator B P such that where it is understood that multiplication by the field f c precedes the action of B P .(The operator B P does not have a functional class subscript c because it redistributes a unit uptake with the same profile regardless of functional class.)We assume that the remineralization of organic matter releases dFe and phosphate in the same ratio with which they were taken up.Therefore, S Fe c = S P c .The values of f c , which set the export efficiency of each class, are not directly constrainable by the data used here.The nutrient concentrations constrain only the total export (i.e., summed over all classes), while the phytoplankton concentrations constrain the uptake, but not the export, of each class.We therefore use the optimized detrital fractions from the work of Dunne et al. (2005): the detrital f c fractions are modeled as decreasing with temperature T so that f c = f 0 c e −k f T , with k f = 0.032 • C −1 independent of class, f 0 sml = 0.14, f 0 lrg = 0.74, and we assign f dia = f lrg .The large and diatom classes are thus ∼ 5 times more efficient at exporting organic matter than the small class.We acknowledge that Richardson and Jackson (2007) suggested that small and large phytoplankton have similar export efficiencies.However, their very sparse data do not provide strong evidence that f sml = f lrg , only that their values are uncertain.Indeed, our state estimates using the f c values of Dunne et al. (2005) are consistent with the data presented by Richardson and Jackson (2007).Plots of the fractional uptake of each class versus its fractional export (not shown here) are broadly consistent with Fig. 1b of Richardson and Jackson (2007).
Following Najjar et al. (1992), we assume that the detrital production rate is fluxed as POP through the base of the euphotic zone at z e = 73.4m with φ POP The redistribution operator B Si similarly injects silicic acid into the aphotic water column with the divergence of the opal flux, φ bSi , which attenuates because of temperature-dependent opal dissolution following Gnanadesikan (1999) and Holzer et al. (2014).For each latitude and longitude, φ bSi is computed as the solution to ∂ z φ bSi (z) = −(κ max Si /w Si ) exp(−T E /T (z))φ bSi (z), with the boundary condition φ bSi (z e ) = 0 z e R Si : P f dia U dia dz.We use T E = 11, 481 K as Gnanadesikan (1999) and the same detrital fraction f dia for the opal export and diatom POP export.The parameter combination κ max Si /w Si has nearly the same value as determined by Holzer et al. (2014), but was re-optimized here for a simple restoring-type model that takes subgrid topography into account (see Appendix B).
The scavenging operators S s,POP and S s,bSi act on J POP and J bSi to redistribute a fraction of the iron scavenged at every layer throughout the water column below.In terms of the corresponding redistribution operators, we write where the fractions f POP and f bSi were both fixed at 0.9 (see Appendix D).The operators B s,POP and B s,bSi in effect "recycle" scavenged iron.They are very similar to B P and B Si but in addition to distributing scavenged iron from the euphotic zone to the aphotic zone, they also redistribute the scaveng-ing rates of every aphotic layer to a source of redissolving iron with the divergence of the scavenging particle fluxes.The flux of scavenged iron into the bottom is assumed to be lost forever so that there would be iron loss even for 100 % efficient recycling of scavenged iron.(For details see Appendices A and B.) To compute accurate particle fluxes for constructing all S operators, we take subgrid topography into account (as done by Moore and Braucher, 2008), using the high-resolution ETOPO2V2c data set (National Geophysical Data Center, 2006).This is done by calculating for each grid box the fractional area occupied by the subgrid topography, which is also the fraction of the particle flux intercepted by the subgrid topography.For each grid box, the fraction of the flux intercepted is instantly remineralized and redissolved for S i c or buried in the sediments for S s,POP and S s,bSi (details in Appendix B).

Uptake rates
The PO 4 uptake rate at a given point is a function of the local temperature T , irradiance I , and nutrient concentrations.The uptake rate for functional class c is calculated as the product of its phytoplankton concentration, p c , and its specific growth rate, µ c , as where τ c is the timescale for growth, p max c is the phytoplankton concentration under ideal conditions, and F I,c and F N,c are dimensionless factors in the interval [0, 1) that represent light and nutrient limitation, respectively, as defined below.We derive Eq. ( 6) similarly to Dunne et al. (2005) and Galbraith et al. (2010) as follows.
First, p c is calculated diagnostically by assuming steady state between growth and mortality, which avoids the need to carry explicit plankton concentration tracers.This is justified by the coarse resolution of our model, which implies transport timescales across a grid box that are much larger than the typical timescales for phytoplankton growth.Based on Dunne et al. (2005)'s mortality formulation, p c can be modeled by a logistic equation where the p c /p * c fraction scales the specific mortality rate λ, and p * c has been referred to as the "pivotal" population density (e.g., Dunne et al., 2005;Galbraith et al., 2010).Equation (7) has a nontrivial steady state, given by We assume that all phytoplankton classes share the same specific mortality rate λ, which depends only on temperature.
For simplicity, we follow Galbraith et al. (2010) and approximate the T dependence of λ to be identical to that of the growth rate, which was determined by Eppley (1972) to be of the form e κT .Thus, λ = λ 0 e κT , where λ 0 is a constant.Our formulation differs from that of Dunne et al. (2005) and Galbraith et al. (2010), who raise the ratio p c /p * c to a power α = 1 or α = 1/3 to differentiate between small and large phytoplankton classes.Here, we instead differentiate between classes by assigning them different half-saturation rates and maximum uptake-rate constants similarly to the work of Matsumoto et al. (2008) (details in Sect. 2.3.1 and 2.3.3).
We model the specific growth rate µ c as multiplicatively co-limited (Saito et al., 2008) by temperature, light, and nutrients: where τ c is the growth timescale at 0 • C under ideal conditions and the temperature dependence e κT (Eppley, 1972) is identical to that used for the mortality rate (e.g., Galbraith et al., 2010).To group parameters for more efficient optimization, we define p max c = p * c /(λ 0 τ c ), so that diagnostic Eq. ( 8) for the phytoplankton concentration becomes Substituting Eq. ( 9) and Eq. ( 10) into U c = µ c p c gives Eq. ( 6), which is similar to the uptake formulations of Doney et al. (2006) and Matsumoto et al. (2008).We note that in the Sea of Japan the model's circulation produces unrealistic nutrient trapping, likely due to underresolved currents.We therefore set the specific growth rate in the Sea of Japan to zero, effectively removing it from the computational domain of the biogeochemical model.

Nutrient limitation
We model the limitation of functional class c by nutrient i by a Monod function (Monod, 1942) of the concentration, where k i c is the half-saturation constant that determines the scale on which the concentration influences uptake (because only diatoms take up silicon k Si lrg = 0 and k Si sml = 0.) For the co-limitation of all three nutrients, we use the type-I multiplicative form (Saito et al., 2008) We chose the Monod model over the arguably more realistic quota model (e.g., Flynn, 2003) Droop, 2009).However, here we prefer the smoothness of the multiplicative formulation because differentiability is a theoretical requirement for Newton's method to converge (e.g., Kelley, 2003a).A product of PO 4 , dFe, and irradiance Monod terms was also used by Parekh et al. (2005) and Dutkiewicz et al. (2006) in their coupled phosphorus-iron model.

Light limitation
We prescribe the irradiance I and model light limitation using a simple Monod factor with half-saturation constant k I,c for class c (e.g., Doney et al., 2006).We use an annual mean I derived from photosynthetically active radiation (PAR) measured over the period 2002-2015 by the Modis Aqua satellite (NASA Goddard Space Flight Center, 2014).The surface PAR at location (x, y), denoted by I 0 (x, y), was converted to W m −2 using 2.77 × 10 18 quanta s −1 W −1 (Morel and Smith, 1974).Irradiance is modeled as exponentially attenuated with depth z so that with k −1 w = 25 m.

Elemental uptake ratios
Because we key all biological production to PO 4 utilization, we must specify the Fe : P and Si : P elemental uptake ratios for the iron and silicon cycles.The Fe : P uptake ratio, R Fe : P , is known to increase and saturate with increasing dFe concentration (e.g., Sunda and Huntsman, 1997).We follow Galbraith et al. (2010) and model the dFe dependence as a simple Monod term where R Fe : P 0 is the maximal Fe : P uptake ratio.As noted by Galbraith et al. (2010), this formulation ignores the effects of light limitation suggested by several studies (e.g., Sunda and Huntsman, 1997;Strzepek et al., 2012).In principle, R Fe : P 0 and k Fe : P could be different for different functional classes.However, constraining class-dependent Fe quotas is beyond the scope of what is possible with our inverse model: different values of R Fe : P 0 and k Fe : P for each class would directly compensate for one another in the global Fe export.
Equation ( 14) effectively encodes a minimum iron requirement of zero.Thus, for very low dFe concentrations, R Fe : P would fall below a realistic cell quota.However, this has no mechanistic consequence, because for such low dFe concentrations there is essentially no uptake in our formulation.This is because R Fe : P is proportional to a dFe Monod term, while P uptake is proportional to the square of a dFe Monod term.Thus, as dFe becomes small, the uptake goes to zero faster than R Fe : P itself.Simply put, this means that when the Fe : P ratio is unrealistically small, it does not matter because there is no P or Fe uptake.When we introduced a nonzero minimum for R Fe : P it tended to be optimized to zero, which means that a simple Monod factor suffices to capture the dFe dependence of the Fe : P ratio.
The Monod formulation (14) does capture luxury iron uptake (e.g., Marchetti et al., 2009a) when the half-saturation constant of Eq. ( 14) exceeds the half-saturation constant of the iron limitation in Eq. ( 11), as made explicit by Galbraith et al. (2010).This is the case for our optimized value of k Fe : P so that phytoplankton has the luxury to increase its iron uptake with increasing dFe concentration even when iron is not limiting.
Our representation of the R Si : P uptake ratio takes into consideration field studies and iron enrichment experiments, which have indicated that in HNLC and upwelling regions iron limitation leads to increased diatom silicification, i.e., increased cellular Si : N and Si : P ratios (e.g., Takeda, 1998;Hutchins and Bruland, 1998;Franck et al., 2000;Brzezinski et al., 2003).However, there is no literature consensus on a mechanistic formulation of the iron dependence of silicic acid uptake.For example, Matsumoto et al. (2013) assume that the Si : N uptake ratio is inversely proportional to the dFe concentration (capped at a minimum), while Jin et al. (2006) assume the Si : N ratio to depend only on the Si(OH) 4 concentration.Others suggest that the dFe concentration only impacts the diatom growth rate and not the cellular Si : C ratio, while the Si(OH) 4 concentration impacts the cellular Si : C ratio and not growth rate (e.g., Marchetti et al., 2009b;Brzezinski et al., 2011a).Here, we chose to retain the effects of increased silicification due to iron limitation and the impact of high Si(OH) 4 concentration on silicification (Brzezinski 2016, personal communication).We model these effects with the formulation The ratio involving χ Fe produces increased silicification when iron is deficient, while the Monod term for χ Si produces increased silicification in silicon-replete environments: if χ Fe → 0 and χ Si k Si Si : P , then R Si : P → R Si m , while if χ Fe k Fe Si : P or χ Si → 0, then R Si : P → R Si 0 .The minimum and maximum Si : P ratios R Si 0 and R Si m , as well as the constants k Fe Si : P and k Si Si : P were tuned rather than fully optimized to achieve the observation-based fractional uptake of each functional class (see Sect. 3.3 on optimization for details).(Plots of the experimental data that show increased silicification under conditions of low dFe can be seen in Fig. 6 of Franck et al., 2000 andin Fig. 7 of Brzezinski et al., 2003.)2.4 Iron model

Iron sources
The aeolian source, s A , is based on the spatial pattern of the surface flux of atmospheric soluble iron of Luo et al. (2008), obtained from an atmospheric model for current climate conditions that includes size-partitioned mineral dust, biomass burning, and industrial emissions.Because the global strength σ A ≡ s A (r) d 3 r of the aeolian source is highly uncertain (e.g., Tagliabue et al., 2016), we scale the global amplitude of this pattern to an initial guess for σ A that is then refined in our final optimization step (see Sect. 3.3 and 3.4 below).We note that the model of Luo et al. (2008) estimates a soluble aeolian iron flux into the ocean of σ A ∼ 6 Gmol yr −1 .(The global source strength of iron type k is defined as σ k ≡ s k (r) d 3 r.) The sedimentary source, s S , has the pattern of the POP flux reaching the sediments (Elrod et al., 2004;Frants et al., 2016) and accounts for both resolved and subgrid topography.The amplitude of this pattern is the global sediment iron source strength σ S , which is an optimized parameter.The dependence of the sediment redox reaction on dissolved oxygen (e.g., Galbraith et al., 2010) is ignored here for simplicity and to avoid carrying oxygen as another tracer.Unlike in the model of Frants et al. (2016), the phosphorus cycle and POP flux are not prescribed but coupled to the iron and silicon cycles as described above.
To model the hydrothermal source, s H , we use the 3 He source pattern of the Ocean-Carbon Cycle Model Intercomparison Project (OCMIP) protocol (Dutay et al., 2004), and jointly optimize the hydrothermal iron source strengths σ H,ATL , σ H,PAC , σ H,IND , and σ H,SO of the Atlantic, Pacific, Indian, and Southern Ocean ridge systems, as in the work of Frants et al. (2016).

Iron sinks
Dissolved iron can be either chelated by ligands or "free".We assume that scavenging acts only on the concentration χ Fe of free iron so that chelation by ligands protects dFe from being scavenged.Scavenging is modeled as a first-order process (e.g., Aumont et al., 2015) so that the scavenging rate is proportional to the product of χ Fe and the concentration of the scavenging particles χ j , for j ∈ {POP, bSi, dst}, the three types of particles considered.For each particle type, the scavenging rate per unit volume is thus modeled as where the constants κ POP scv , κ bSi scv , and κ dst scv are optimizable parameters.
To compute the concentration of the scavenging particles, we use the fact that the flux divergences generated by the biogenic transport operators must be balanced by local rem-ineralization or dissolution rates; that is, and Although we use the nominal values of κ P and κ max Si listed in Table 1, note that these constants only enter the scavenging rates (Eq.16) through the combinations κ POP scv /κ P and κ bSi scv /κ Si , where κ POP scv and κ bSi scv are optimized.The concentration of dust particles is modeled as vertically uniform due to sinking mineral dust that does not dissolve or resuspend from sediments (e.g., Moore and Braucher, 2008).We use the geographic pattern of the dust mass flux into the ocean provided by Luo et al. (2008), which we convert to a particle concentration using a nominal sinking speed of w dst = 50 m day −1 .The exact value of w dst does not matter because the dust scavenging rate depends only on κ dst scv /w dst , and κ dst scv is optimized.The key control on shaping the free-iron concentration, and hence the scavenging, is the ligand concentration L. Chemical equilibrium between ligands, total dFe, and free iron determines χ Fe as a quadratic function of the (total) dFe concentration (see, e.g., Frants et al., 2016).We used the same ligand stability constant of K L = 8 × 10 10 kg (mol Lig) −1 as Frants et al. ( 2016).The ligand concentration itself is modeled to have a uniform background value L b that can be enhanced in old waters (Misumi et al., 2013) and in hydrothermal plumes (e.g., Bennett et al., 2008;Hawkes et al., 2013), similar to the formulation of Frants et al. (2016).Specifically, we use where L H and L sw are the elevated hydrothermal and aged "seawater" ligand concentrations, which are modeled as follows.The hydrothermal ligand plumes are computed from the source-sink balance where v is a mask that is unity for grid boxes containing vent sites (taken from the OCMIP 3 He source from the study of Dutay et al., 2004) and zero elsewhere.The timescale τ v = 1 s clamps the ligand concentration to L v at the vents, and the timescale τ b controls the plume spread by setting the rate with which L H decays away from the vents.The ligand concentration L sw is enhanced in old waters according to where (r) is the ideal mean water age (easily computed for our model).We choose max = 1600 years following Frants As is the case for most iron models, there is no need to explicitly represent the chemical precipitation of dFe.This is because in most formulations the scavenging rates increase rapidly when dFe exceeds a certain threshold.For our model this threshold is set by the ligand concentration L: when dFe concentrations exceed L, χ Fe and hence the scavenging rates rise rapidly.
3 Numerical method, parameter optimization, and family of state estimates

Steady-state solution
All three-dimensional fields (e.g., the concentrations χ i ) are discretized on our model grid and organized into column vectors (length n ∼ 2 × 10 5 at our resolution).Linear operators such as T , S i c , and S s,j are correspondingly organized into n × n sparse matrices.The steady-state tracer Eqs. ( 1)-( 3) then become a 3n × 3n system of equations that are nonlinear because of the iron scavenging and the co-limitation of the PO 4 uptake.
The 3n × 3n system is solved efficiently using Newton's method (e.g., Kelley, 2003a, b).Convergence of the Newton method depends on the initial guess for the solution and is not guaranteed.For the initial guess of χ P and χ Si we use the annual mean fields of the World Ocean Atlas (WOA13, Garcia et al., 2014) interpolated to our grid, and for the initial guess of χ Fe we use the dFe fields estimated by Frants et al. (2016).The Newton solver typically converges to numerical precision in ∼ 10 iterations.

Cost function
We optimize the model parameters by systematically minimizing a quadratic cost function of the mismatch between modeled and observed fields.For PO 4 and Si(OH) 4 , for which gridded climatologies are available, we define the weights based on the grid box volumes, organized into vector v, as , and where we have normalized the weights by the total ocean volume V and the squared global mean observed concentrations.This nondimensionalizes the quadratic cost terms and scales them to the same order of magnitude.For dFe, for which only sparse observations are available, we also define weights w Fe based on grid box volumes, but observations that are part of a vertical profile receive additional weight as detailed in Appendix C. With diagonal weight matrix W i = diag(w i ) for nutrient i, its cost for the mismatch with observations is then given by where δχ i ≡ χ i − χ obs i .For χ obs P and χ obs Si we use WOA13 fields interpolated to our grid, and for χ obs Fe we used the GEO-TRACES intermediate data product (Mawji et al., 2015) and the data set compiled by Tagliabue et al. (2012).
The cost terms for the nutrient mismatch do not provide a strong constraint on the relative sizes of the phytoplankton classes because the nutrients are determined by their combined export.We therefore include additional terms in our cost function that constrain the phytoplankton concentrations p c to the recent satellite derived estimates of Kostadinov et al. (2009).These estimates provide phytoplankton concentrations for picophytoplankton (0.5-2 µm in diameter), nanophytoplankton (2-20 µm), and microphytoplankton (20-50 µm), which we identify with our small, large, and diatom functional classes.We use the entire mission composite data set as the satellite climatology (Kostadinov et al., 2016).
Because of the large dynamic range of the phytoplankton concentrations, we consider mismatches in the log of the concentrations; that is, where p obs c is introduced to limit the logarithm where the phytoplankton concentration falls to zero, and p 0 = 1 mg C m −3 nondimensionalizes the argument of the logarithm.For each class, we construct normalized weight vectors: where V eup is the global euphotic volume.
Organizing mismatches into vectors and weights into diagonal matrices, we calculate the cost for the phytoplankton concentration mismatch as and combine the costs for the nutrient and plankton mismatches into the total cost: which we minimize to constrain our model parameters by the available observations.In Eq. ( 26) the ω weights were chosen such that the four cost terms contribute roughly equally to the total cost for a typical member of our family of state estimates.This was achieved with (ω P , ω Si , ω Fe , ω plk ) = (1, 0.47, 0.044, 0.30), the smaller weight for dFe reflecting its larger root mean square (rms) mismatch and hence much larger cost.

Optimization strategy
Our model has ∼ 50 biogeochemical parameters that can in principle be determined through objective optimization given appropriate observational data.However, even with perfect data, some parameters can compensate for others (e.g., two parameters appearing as a ratio) so that not all parameters are independent.Other parameters cannot be optimized because the mismatch with available nutrient and phytoplankton data is not sensitive to their value.In practice, it therefore is not possible to optimize all parameters, and care is needed to optimize only those parameters that independently shape the nutrient and phytoplankton concentrations.
The parameters associated with the remineralization of phosphate and dissolution of opal are well constrained by the high-quality climatologies of PO 4 and Si(OH) 4 .However, the iron cycle is relatively poorly constrained because the dFe data are much more sparse in both time and space, and estimates of the iron sources range over 2 orders of magnitude (e.g., Tagliabue et al., 2016).Moreover, the ligand field that determines the scavengable free iron is highly uncertain.Given these challenges, the recent inverse model of the iron cycle by Frants et al. ( 2016) considered a family of state estimates for a range of external source strengths, an approach we will follow here for our coupled model.
Another key consideration is computational cost.Even with the numerically efficient Newton solver, optimization typically requires hundreds of solutions of Eqs. ( 1)-( 3) per optimized parameter.We therefore optimized no more than 13 parameters at a time.We acknowledge that the minimum attained by sequentially optimizing groups of independent parameters is generally different from jointly optimizing all independent parameters, but computational and practical considerations demanded a sequential approach.We justify this a posteriori by the fact that we are able to achieve fits to the observed nutrient concentration fields with rms mismatches similar to those of other recent data-constrained models (e.g., Primeau et al., 2013;Holzer et al., 2014;Frants et al., 2016).Given these considerations, we adopted the following strategy: i. Parameters that are measurable and considered wellknown, as well as parameters that are unconstrainable by our cost function or with values that are not critical, because they are strongly compensated for by other parameters, were assigned values from the literature as collected in Table 1.The considerations that entered our choice of prescribed parameters are detailed in Appendix D.
ii.The parameters that set the phosphate remineralization and opal dissolution profiles were optimized by minimizing the mismatch with PO 4 and Si(OH) 4 concentration data from the WOA13 using separate singlenutrient models.For the Si cycle, we used the model of Holzer et al. (2014) and verified that the opal sinking speed parameter w Si was not affected by the inclusion of subgrid topography (Appendix B).For the P cycle, we used a similar conditional restoring model without POP but with subgrid topography, and optimized the Martin exponent b.The resulting values of w Si and b (Table 1) were held fixed for all optimizations of the coupled nutrient cycling model.
iii.The remaining parameters were optimized using our coupled model.We first assign initial values for all these parameters and then sequentially update these initial values by optimizing subsets of parameters as detailed in Appendix D. Both initial and final optimized parameter values are collected in Table 2. (For the parameters of the iron cycle, Table 2 gives the values of our typical state estimate and the range across a family of state estimates, the different members of which have different external iron sources.)

Family of state estimates
Figure 1 shows the quality of the fit to nutrient and phytoplankton data for all our optimized state estimates, which span a wide range of source strengths.For ease of presentation, state estimates are divided at σ H = 1 Gmol Fe yr −1 into low and high hydrothermal cases, with σ H ranging from 0.073 to 11. Gmol yr −1 .For high σ H , we focused on correspondingly higher aeolian and sedimentary source regimes.Source-parameter space was not explored uniformly because (i) the final step of our optimization adjusted our initial choice of sources, and because (ii) some source choices produced spurious numerical difficulties for the Newton solver.All state estimates fit the macronutrient fields about equally well, but the overall quality of fit as quantified by the square root of the quadratic mismatch ("total cost", top panels of Fig. 1) gets systematically worse with increasing aeolian source strength, σ A , especially for high hydrothermal sources.This worsening fit for high σ A is reflected in the mismatch of all three nutrients.The state estimates with a total cost that exceeds the smallest misfit by ∼ 5 % (corresponding to a total cost above 15.4) are discarded from our family.This essentially eliminates state estimates with σ A 22 Gmol yr −1 and σ H 5 Gmol yr −1 (black crosses in Fig. 1).(If we include the "crossed-out" states in plots below that show scatter across the family of state estimates, the visual impact is virtually imperceptible.)While it is clear from Fig. 1 that high-σ A states are less likely, we hasten to add that the cost threshold for inclusion in the family is arbitrary as we do not have a formal error covariance to convert the cost into a likelihood.
For a small fraction of our state estimates, the optimization pushed the maximum possible Fe : P ratio R Fe : P 0 to nearzero values.These cases are unrealistic because zero R Fe : P 0 means significant P uptake and export are maintained without Fe uptake.We therefore also exclude cases for which the optimized R Fe : P 0 falls below 0.5 mmol Fe (mol P) −1 from our family of state estimates.Removing these unphysical outliers has negligible visual impact on plots that show the entire family of state estimates.In terms of total cost, there is little sensitivity to the strength of the sedimentary source -scavenging can be optimized for a sedimentary source ranging over 2 orders of magnitude for an overall similar quality of fit.For low σ H , there are small opposing rms mismatches for PO 4 and dFe, with a slightly better PO 4 fit for higher sedimentary source and a slightly better dFe fit for lower sedimentary source, although the variation in the mismatch is less that 1 % of the global mean concentrations.
While the mismatch for dFe is substantial at ∼ 45 % of the global mean dFe concentration, the smallest dFe mismatch occurs when all three sources are low.The dFe mismatch rapidly increases with σ A , consistent with the findings of the much simpler model of Frants et al. (2016).The overall cost and the mismatch for each nutrient are insensitive to the strength of the hydrothermal source.
While Fig. 1 shows some variations with the source strengths in the overall quality of the fit, it is clear that the iron sources and scavenging sinks are poorly constrained by the available nutrient and phytoplankton observational data.Given the uncertainties in the sources and the small cost differential between family members, it is not appropriate to single out the state with the lowest numerical cost as being the most realistic.We therefore use the entire family of state estimates below to assess the robustness of our results and to elucidate the systematic variations of the carbon and opal exports with the fractional source of each iron type.The uncertainty in the value of any metric is assigned from its spread across the family.
As a typical representative of our family of state estimates, for which we plot results below, we selected the state with (σ A , σ S , σ H ) = (5.3,1.7, 0.9) Gmol Fe yr −1 .This state is typical in that it lies at the mode of the distribution of overall rms misfit values and, for most quantities, tends to lie in the middle of the range across the family.
We emphasize that the variations across the family of state estimates explored here are variations of the fully optimized biogeochemical states.These variations cannot be used to infer the system's response to dFe perturbations for which the other biogeochemical parameters would not change.Such perturbations, which are of great interest in themselves, are beyond the scope of this paper and will be examined in a separate publication.

Fidelity to observations
We now examine in more detail how well our estimates match the observations against which they were optimized.We focus here on the nutrient fields, which contribute the bulk of the cost function.Each phytoplankton field contributes only ∼ 10 % to the cost function and serves primarily to differentiate between the uptake of small and large phytoplankton -a comparison between estimated and observed phytoplankton concentrations is provided in Appendix E.Where there is little variation across the family, we focus on our typical state estimate.For iron-related quantities that have by construction significant spread across the family, we will focus on the systematic variations of the optimized states with the dFe sources.

Nutrient concentrations
The nutrient concentrations are well constrained for all members of our family of state estimates.We quantify the overall fit of the modeled nutrient concentrations in terms of the joint probability density function (pdf) of the modeled and observed concentrations.This joint pdf may be thought of as the binned scatter plot of the modeled versus observed values for all grid boxes.The binning for each nutrient was weighted using the associated cost-function weights, w i .These joint pdfs are shown in Fig. 2 for all three nutrients for our typical state estimate.Both the PO 4 and Si(OH) 4 pdfs fall close to the 1 : 1 line, showing high fidelity to observations.For PO 4 the cost-weighted rms error is ∼ 5 % of its global mean of 2.17 µM.In comparison, Primeau et al. (2013) achieved an rms mismatch of 3 % by jointly optimizing the uptake rate of each grid box with the circulation.Silicic acid has a larger rms mismatch of ∼ 12 % relative to its global mean of 89.1 µM.This is similar to the 13 % rms error reported by Holzer et al. (2014), who used the same circulation but a much simpler model of the silicon cycle.
The global mean dFe concentration is well constrained within the narrow range of 0.56-0.68nM across the family of state estimates.For iron, the joint probability is by necessity computed using only those grid boxes that contain dFe observations.The scatter from the 1 : 1 line is much larger than for the macronutrients with a substantial rms mismatch of 0.29 nM or 44 % of the mean.This mismatch is comparable to that of other models (e.g., Tagliabue et al., 2016).Compared to the simpler model of Frants et al. (2016), the joint pdf shows that our dFe field has a wider, more realistic dynamic range.We note that, while Frants et al. ( 2016) report an rms mismatch of only 0.19 nM, they also employed different weights for the model-observation mismatch.If we recompute the rms mismatch of the optimized dFe field of Frants et al. ( 2016) using the weights of this work, we also obtain a 0.29 nM mismatch.
The relatively large mismatch for dFe not only quantifies model deficiencies, but to a large degree also reflects the fact that we are comparing snapshot observations against a steady-state coarse-resolution model.The dFe observations  2012) and from the GEOTRACES Intermediate Data Product 2014 (Mawji et al., 2015).have difficult-to-quantify temporal and spatial sampling biases, and dFe being a trace element is sensitive to episodic events in the aeolian source (e.g., Croot et al., 2004), and possibly to internal episodic events such as submarine volcanism (e.g., Massoth et al., 1995).

dFe profiles
To quantify the spatial structure of the dFe mismatch, we examine vertical profiles for each basin.For both model and observations, we only use the grid boxes that contain observations and average horizontally over the basins using the cost weights w Fe .The resulting profiles are shown in Fig. 3.The family of model profiles generally overlaps with the observational uncertainties.The estimates are particularly close to the observations near the surface.In the abyssal oceans, the spread in the family of profiles is larger.This spread is in part a reflection of the weights in our cost function.Most dFe observations are available in the upper ocean, implying a small variance of the mean concentration and hence large weights, while deep observations tend to be sparser with smaller weights (for details on the weights see Appendix C).
Figure 3 also shows systematic biases in the inferred dFe concentrations.Biases are particularly strong in the Pacific, where the observations tend to be underpredicted by as much as ∼ 0.3 nM above ∼ 1500 m and overpredicted by ∼ 0.2 nM below ∼ 2000 m depth.The typical estimated Pacific profiles are too linear in the upper 1500 m, with vertical gradients that are too weak above ∼ 300 m and too strong below ∼ 1000 m.In the Atlantic, a smaller low bias of ∼ 0.15 nM can be seen between ∼ 500 and ∼ 1300 m depth.
These biases could be due to deficiencies in our model such as oversimplified ligand parameterization, but one must also keep in mind that there are hard-to-quantify biases in the observations.The observations are too sparse to form a reliable climatology, and it is remarkable that we can fit the available observations as well as we do.The larger biases in the Pacific could well be due to the absence of Pacific transects in the GEOTRACES Intermediate Data Product 2014 (Mawji et al., 2015), which means that mismatches in the Pacific incur a relatively smaller penalty in our cost function.
Figures 1d, 2, and 3 are appropriate quantitative comparisons between estimated and observed dFe, given that essentially raw bottle data are compared with a coarse-resolution steady-state model.For completeness, Appendix G also compares the main transects included in the GEOTRACES Intermediate Data Product 2014 with our typical state estimate.For a given phytoplankton functional class, different nutrients are known to limit biological production in different parts of the ocean (e.g., Moore et al., 2001).These geographic limitation patterns are a fundamental fingerprint of upper-ocean ecosystem dynamics.Knowledge of the limitation patterns is important for understanding how the global nutrient cycles operate in the current climate and for assessing possible future changes of the global ocean ecosystem.
Limiting nutrients can be determined observationally (e.g., Moore et al., 2013), and from biogeochemical models (e.g., Moore et al., 2004).Here, we estimate the limitation patterns from our optimized inverse-model state estimates.In our model, the biological uptake of each functional class (Eq.6) is limited through F N,c , the product defined in Eq. ( 11) of three Monod terms, one for each nutrient.We define the deficiency D i c of functional class c in nutrient i as the complement of the corresponding Monod factor, i.e., as . We deem nutrient i to be "limiting" class c if D i c > 0.5 or, equivalently, if χ i < k i c ; i.e., if the nutrient concentration falls below its half-saturation value for uptake.
To display the pattern of the nutrient limitations, we could use the fact that we have three nutrients with which to define an RGB color as (D P c , D Si c , D Fe c ).However, because the resulting colors vary continuously, it is hard to quantify the resulting patterns.We therefore define the limiting RGB color as (L P c , L Si c , L Fe c ), where L i c = 1 if D i c > 0.5 and L i c = 0 otherwise.This partitions the RGB color cube into eight colors that define the eight limitation regimes shown in Fig. 4. Specifically, the color is black (0, 0, 0) if all nutrients are available in sufficient quantities so that none are deemed limiting, white (1, 1, 1) if all three nutrients are limiting, red (1, 0, 0) if only dFe is limiting, green (0, 1, 0) if only Si(OH) 4 is limiting, and blue (0, 0, 1) if only PO 4 is limiting.The remaining three possibilities correspond to two nutrients being co-limiting: magenta (1, 0, 1) if dFe and PO 4 are colimiting, cyan (1, 1, 0) if PO 4 and Si(OH) 4 are co-limiting, and yellow (0, 1, 1) if Si(OH) 4 and dFe are co-limiting.
Figure 4 shows the limitation patterns of all three phytoplankton classes.The large and diatom classes have similar patterns of iron limitation in the Southern Ocean, eastern tropical Pacific, and North Pacific.For both classes, the Indian Ocean and North Atlantic are largely PO 4 limited.The subtropical gyres of the Indian Ocean and North Atlantic are PO 4 limited for the large class and PO 4 -Si(OH) 4 colimited for diatoms.The differences between the large and diatom classes stem from the Si(OH) 4 requirement of diatoms.Because the large class requires no silicic acid, its limitation map has no areas where all three nutrients are limiting (white).The subtropical gyres of the Pacific and South Atlantic are dFe and PO 4 co-limited for the large class, while for diatoms the center of these gyres are limited in all three nutrients.For diatoms, the eastern margins of the Pacific sub-tropical gyres show Si-Fe co-limitation (yellow).Only a few grid boxes in the Arctic are solely limited by silicic acid (green).The completely nutrient replete regions of the Arctic and Weddell Sea reflect low biological utilization driven in our model by light limitation through the prescribed irradiance field, I .
The small phytoplankton class shows a much simpler pattern.Limitation occurs primarily in the subtropical oceans with small patches of iron limitation also in the Southern Ocean and tropical Pacific.Iron limitation dominates the subtropical South Pacific, while PO 4 limitation occurs primarily in the subtropical gyres of the southern Indian Ocean and North Atlantic.The rest of the ocean is largely nutrient replete for the small class.The limitation patterns are robust across our family of state estimates, with areas of each type of limitation generally varying by 5 % or less.
The general features seen in Fig. 4 broadly agree with the observational data (in situ and bottle nutrient addition experiments) reported by Moore et al. (2013).Similar to our estimates, the observations show Fe limitation in the Southern Ocean, subpolar North Pacific, and eastern tropical Pacific.The observations also indicate Fe limitation in the North Atlantic, which for our state estimates is also present in small patches in the western subpolar North Atlantic and becomes slightly more pronounced for the states with a higher total iron source.Moore et al. (2013) report Si limitation in the Pacific sector of the Southern Ocean at its northern boundary, where silicic acid concentration sharply decreases.This is consistent with our yellow region of joint Si and Fe limitation along the eastern edge of the Pacific subtropical gyres.Consistent with our estimates, the observations show PO 4 limitation in the North Atlantic subtropical gyre and in the equatorial Atlantic.
Our limitation patterns can also be compared to those calculated for summer conditions in the BEC model of Moore et al. (2004).However, it must be kept in mind that (i) the BEC model has a different circulation and a different formulation of biogeochemical cycles (e.g., explicitly representing the nitrogen cycle and diazotrophs) and that (ii) Moore et al. (2004) define limitation in terms of the minimum Monod factor, while we use a threshold of 1/2 for the Monod factors and jointly consider three Monod terms to define the type of limitation.For diatoms, the Fe limitation pattern reported by Moore et al. (2004)  with our definitions there is very little PO 4 limitation in the Pacific for the small class, which is iron limited or nutrient replete in most of the Pacific.In addition, we note that the annual-mean nature of our estimates is another possible reason for the differences.

Export production
A key metric of the nutrient cycles is their export production, which determines the strength of the biological pump (e.g., Pasquier and Holzer, 2016).Export production is not directly available from satellite measurements, but observationally constrained estimates are easily calculated from our inverse model.The phosphorus export flux, P , is simply the flux of organic phosphorus into the aphotic zone that is remineralized there, which we compute using the operators S P c (sinking and remineralization) as For plotting, we convert P to a carbon export flux using a constant C : P ratio of 106 : 1.We acknowledge that this simple unit conversion underestimates the true carbon export, because we do not explicitly represent DOC.Semi-labile DOC has a longer typical lifetime than DOP, resulting in an effectively larger C : P ratio for dissolved organic matter (DOM) than for particulate organic matter.Using the dataassimilated phosphorus cycle of Primeau et al. (2013), which explicitly carries both PO 4 and DOP, and applying a C : P ratio of 225 : 1 for DOP (as determined by the DOM OPT simulation of Letscher et al. (2015) for semi-labile DOM), we estimate that the simple unit conversion of the POP export (Eq.27) underestimates the carbon export by ∼ 12 %.We similarly calculate the opal export as and the iron export associated with the remineralization of organic matter as where the vertical integrals are over the model aphotic zone (bottom to 73.4 m depth).
Figure 5a shows a map of the phosphorus export flux (converted to carbon units), together with its zonal integral for each member of our family of state estimates.The spatial pattern shows some differences with the estimate of Primeau et al. (2013) (blue curve in Fig. 5a).Our estimate of the carbon export has 1.5-2 times larger tropical and highlatitude peaks, but is closer to the satellite-derived estimates of Dunne et al. (2007).Our estimate also has sharper meridional gradients, which can also be seen in satellite-derived estimates of production (e.g., Frants et al., 2016).Our globally integrated phosphorus export of 7.5-8.6Tmol P yr −1 (9.5-11.Pg C yr −1 ) is also larger than the 7.4 ± 2.5 Pg C yr −1 estimate of Primeau et al. (2013) 2013) consider the phosphorus cycle in isolation and optimize a single spatially varying uptake timescale for each grid box, while we explicitly represent three phytoplankton functional classes with different, optimized globally uniform uptake timescales, τ c .We note that, if we use the same growth timescale for each phytoplankton class, our model's phosphorus export remains close to that of Primeau et al. (2013).
Our estimates of export production compare well with the satellite-based estimates of 9.7-12.Pg C yr −1 by Gnanadesikan et al. (2004).Our estimates also lie within the wide range of 9-28 Pg C yr −1 of the Ocean-Carbon Cycle Model Intercomparison Project 2 (OCMIP-2, Najjar et al., 2007), and compare well the OCMIP-2 mean particle export of 13 ± 3 Pg C yr −1 .(Because our model does not carry DOM, we compare our phosphorus export in carbon units to both the export production and the particle export reported by OCMIP-2.) Our estimates of phosphorus export in the subtropical gyres compare well to the POP exports of Letscher et al. (2016) in spite of the fact that we do not explicitly represent DOP.Using the same masks to define the subtropical gyres, we estimate a global mean subtropical phosphorus export of 10. ± 1. mmol P m −2 yr −1 (mean and standard deviation across our family of estimates), while the estimate of Letscher et al. (2016) is 10. ± 2. mmol P m −2 yr −1 .This underscores the fact that our inverse model implicitly captures the effects of DOP lateral transport and utilization, which Letscher et al. (2016) estimate to contribute (29 ± 9) % of the subtropical phosphorus export.
Figure 5b shows a map of the opal export, together with its zonal integral.As expected, opal export is most pronounced at high latitudes, particularly in the Southern Ocean.In spite of our relatively complex formulation of silicic acid utilization in terms of co-limitations, the spatial pattern of the opal export and its global total of 164.-177.Tmol Si yr −1 compare well with the estimates by Holzer et al. (2014) (171 ± 31 Tmol Si yr −1 ).Other estimates of the global opal export range from 69 to 185 Tmol Si yr −1 (e.g., Moore et al., 2004;Sarmiento et al., 2007;Heinze et al., 2003).
There is very little spread in the carbon and opal export productions across our family of state estimates as can be seen by the tightly clustered zonal integrals plotted in grey in Fig. 5a,b.This shows that the carbon and opal exports are well constrained despite the wide range of iron inputs.
Figure 5c shows a map of the iron export associated with organic matter, but not including the iron export carried by scavenging particles.The phosphorus and iron exports have broadly similar patterns, with differences that reflect variations in the local Fe : P uptake ratio.In the irondeficient Southern Ocean, the Fe : P ratio is smaller than its global mean, which results in Southern Ocean iron export that is less efficient than for phosphorus (for iron, the peak Southern Ocean export relative to the tropical peak is lower than for phosphorus).As expected from the widely varying iron source strengths across our family of estimates, the globally integrated iron export spans a wide range of 0.87-5.6Gmol Fe yr −1 .However, the geographic pattern of the iron export is robust across the family: the zonally integrated iron exports normalized by their global integrals collapse onto a well-defined cluster of curves.The spread in the thus normalized iron export is similar to the spread in the (unnormalized) phosphorus export, but slightly larger due to variations in the Fe : P ratio.
All export fields of Fig. 5 show near-zero export in the Weddell Sea, in contrast to what restoring-type models tend to show.For example, the opal export estimated by Holzer et al. (2014) has a local maximum in the Weddell Sea.The Weddell Sea minimum here is due to near-zero satellite measurements of PAR in this region.This may well be an artifact of the satellite data, for which the irradiance in the Weddell Sea varies substantially depending on which years are averaged.

Iron cycle
Here we document some of the key features of the iron cycle as constrained by our inverse model.Certain features such as the dFe concentration field are robustly constrained by the observations regardless of iron source strengths, while other features, such as the relative importance of hydrothermal iron, vary systematically with the source strengths.

Iron sources and sinks
The pattern of the aeolian source is identical for all family members because we only vary its global source strength, σ A .The sediment source is keyed to export production, which is well constrained across the family of state estimates.Therefore, the sedimentary iron source patterns are very similar across all state estimates, with only the global strength σ S varying among state estimates.The initial hydrothermal pattern is set by the OCMIP 3 He source (Dutay et al., 2004), but for total hydrothermal sources larger than ∼ 0.5 Gmol yr −1 , the optimized contributions from each basin changed substantially.Across our family of estimates the mean and standard deviations of the percentage contributions from each basin to the total hydrothermal source are (15 ± 9) % for the Atlantic, (52 ± 6) % for the Pacific, (16 ± 2) % for the Indian Ocean, and (17 ± 3) % for the Southern Ocean (south of 40 • S).For reference, the vertically integrated iron sources of our typical state estimate are plotted in Appendix H.
Because of the small variations in the source patterns, the patterns of the vertically integrated total sinks of dFe, dz (1 − S P s )J POP , dz (1 − S Si s )J bSi , and dz J dst , also vary little across the family of state estimates (see Appendix H for plots of these quantities for our typical state estimate).Note that these sinks balance the total source exactly in steady state.For our typical state, scavenging by POP and opal each account for about half of the total iron sink.The patterns of POP and opal scavenging are determined by the phosphorus and opal exports and by the concentration of free iron.Consequently, the POP scavenging sink is strongest in the tropics, and the opal scavenging sink is strongest in the Southern Ocean.The sink due to mineral dust scavenging reflects the pattern of the aeolian dust input modulated by the free-iron concentration.For our family of state estimates, the sink due to dust scavenging is essentially negligible, being about 3 orders of magnitude smaller than the POP and opal scavenging.
We note that the partition of scavenging among the different particle types cannot be inferred robustly from our in-verse model.This is because the nutrient and phytoplankton data used do not provide separate constraints on the scavenging by each particle type, only on the total amount of scavenging.Moreover, scavenging by one particle type can be compensated for by another type because of overlap in their spatial patterns.However, the partition among particle types does vary systematically across our family of estimates.Scavenging by dust is negligible for all state estimates, while the fraction scavenged by POP ranges from ∼ 10 % for the lowest iron sources to saturation near 100 % for the highest iron sources.(The complementary fraction is due to opal scavenging.)

dFe Concentration and source attribution
Figure 6 shows the typical state's zonally averaged dFe concentration for each basin and for the global ocean.For each zonal average, Fig. 6 also shows the corresponding profile of horizontally averaged dFe for each member of our family.The profiles are tightly clustered showing that the large-scale features of the dFe field are well constrained despite the large variations of the iron sources.The inverse model fits the observed dFe field for widely different sources by adjusting the corresponding scavenging.While these adjustments keep the total dFe field close to the observations, the relative contributions from the aeolian, sediment, and hydrothermal sources are unconstrained and can vary widely.
We calculate dFe concentrations due to each source following Holzer et al. (2016) by replacing the dFe concentration tracer Eq. ( 3) by an equivalent linear diagnostic system that has the same solution.This linear system, corresponding to a given solution of the full nonlinear system, is obtained by replacing the iron uptake and scavenging by linear operators.Specifically, the dFe uptake R Fe : P U c is replaced with L U,c χ Fe and the scavenging rate J j with L J,j χ Fe , where the linear operators, organized into matrix form, are simply specified from the uptake and scavenging rates of the nonlinear solution as L U,c ≡ diag(R Fe : P U c /χ Fe ) and L J,j ≡ diag(J j /χ Fe ).The dFe concentration, χ k Fe , due to source s k (with k ∈ {A, S, H}) is then computed by replacing the total source in the linear equivalent system with s k .
Figure 6 also shows the profiles of the individual source components of dFe, color coded according to the fractional strength of the aeolian source.In contrast to the profiles of the total dFe, these individual source components vary widely across the family of state estimates, but such that the total concentration χ Fe = χ A Fe + χ S Fe + χ H Fe is relatively tightly constrained.For example, for low aeolian sources (yellow profiles in Fig. 6), the low concentration of aeolian iron χ A Fe is largely compensated for by a larger sediment contribution χ S Fe .The concentrations of hydrothermal iron vary less systematically with the aeolian source, but all family members have very similarly shaped hydrothermal dFe profiles.The amplitudes of the hydrothermal dFe profiles can be seen to vary by roughly an order of magnitude across the majority of www.biogeosciences.net/14/4125/2017/Biogeosciences, 14, 4125-4159, 2017 states, effectively fine-tuning the total dFe concentration to be as close to the observations as possible.
7.3 Iron-type attributed export

Phosphorus export
We quantify the contribution of each iron type to the export production as follows.In our formulation, nonzero dFe is necessary for nonzero phosphate uptake U c .The uptake U c (r) inferred at point r is supported by the dFe concentration at r, which is a mixture of aeolian, sedimentary, and hydrothermal dFe.Thus, the uptake supported by iron type k is given by (χ k Fe (r)/χ Fe (r)) U c (r); that is, the local uptake supported by dFe of type k must be in proportion to the concentration fraction χ k Fe /χ Fe .(Note that k χ k Fe /χ Fe = 1.)For a given nonlinear solution, the phosphorus export production supported by iron type k, denoted by P k , is therefore calculated by replacing the uptake U c in Eq. ( 27) with (χ k Fe /χ Fe ) U c (r).
While the total export production is well constrained regardless of the iron source strengths, the production sup-ported by a given iron type varies substantially with the magnitude of the corresponding source.(Summing over the three iron types yields the well-constrained total.)However, regardless of the source amplitudes, the geographic patterns of the export supported by each iron type is similar across the entire family of state estimates.
Figure 7 shows P k ≡ P k / P k , which is the export flux supported by iron type k normalized by the global mean export P k .The P k patterns are plotted for our typical state estimate, together with zonal averages of P k for all family members.The zonally averaged patterns can be seen to differ little among family members.Even for the pattern of the export supported by hydrothermal dFe, which varies most, all family members share the broadly similar features of peak export in the Southern Ocean, with secondary peaks in the tropics and in the Northern Hemisphere subpolar oceans.
Figure 7a shows that aeolian iron supports export primarily in the tropics and in the subpolar oceans.The tropics receive direct input of fresh aeolian iron, while the subpolar oceans receive upwelling regenerated iron (Holzer et al., 2016).The aeolian-iron-supported export pattern is very sim- ilar to the pattern of the total export flux shown in Fig. 5a.(Note that here we plot zonal averages, while Fig. 5 shows zonal integrals.)For sedimentary dFe to support export, it must be transported from the ocean bottom into the euphotic zone.Consequently, the pattern of the export supported by sedimentary dFe (Fig. 7b) is dominated by regions of upwelling in the tropical and subpolar oceans and by regions of shallow depth (both resolved and subgrid), where there is high organic matter flux, such as the seas around Indonesia.The pattern of export supported by hydrothermal dFe (Fig. 7c) is dominated by the Southern Ocean, where most of the density classes into which hydrothermal fluid is injected outcrop.Secondary regions of hydrothermal-iron-supported export are associated with upwelling in the tropics and in the subpolar oceans of the Northern Hemisphere.
Underscoring the similar source distribution of hydrothermal dFe and mantle 3 He, the pattern of hydrothermal-iron- Figure 8. Percent global phosphorus export (equivalently carbon export) supported by each iron type (aeolian, sedimentary, hydrothermal) versus the corresponding fractional source of that iron type.The superposed lines are least-squares fits to theoretical relationships with fixed relative export-support efficiencies.(See text for details.) supported export production is similar to the pattern with which mantle 3 He outgasses to the atmosphere (e.g., Holzer et al., 2017).We do not expect an exact correspondence in the patterns because hydrothermal dFe is subject to scavenging losses, while 3 He is not, and our ratio of hydrothermal dFe source to mantle 3 He source is different for different basins.(The ranges of the ratio of the optimized hydrothermal iron source to the mantle 3 He source across the family were 0.00087-3.3,0.098-8.1,0.2-15., and 0.025-2.8 in units of Mmol Fe (mol 3 He) −1 , for the Atlantic, Pacific, Indian, and Southern Ocean basins, respectively.) While the total phosphorus export is well constrained and varies little across our family of state estimates, the magnitude of a given iron-type-supported export production varies systematically with the corresponding fractional source strength.To quantify these systematic variations, Fig. 8 plots the fraction φ P k ≡ P k / P of the globally averaged iron-type-k-supported export to the total global export as a function of the corresponding fractional global iron source strengths σ k ≡ σ k /σ tot , where σ tot ≡ k σ k .Note that, if a given source strength σ k constitutes 100 % of the total, then it must support 100 % of the export and that if σ k = 0 then it supports 0 % of the export; i.e., the relationship between φ P k and σ k must pass through the points (0, 0) % and (100, 100) %.
Figure 8 shows that, depending on the state, aeolian dFe supports ∼ 20-100 % of the global export, with the low end of the range corresponding to an aeolian source of only ∼ 5 % of the total source.(We did not explore lower fractional aeolian sources.)Sedimentary iron supports ∼ 0-80 % of the global export, with the high end of the range corresponding to a sediment source as high as ∼ 90 % of the total source.B. Pasquier and M. Holzer: Inverse-model of coupled Fe-P-Si cycles Hydrothermal iron supports the least export, ranging from ∼ 0-18 % for hydrothermal sources as large as ∼ 45 % of the total source.
The key point of Fig. 8 is that aeolian iron can be considered to be the most efficient type of iron for supporting export production: for a given fraction of the total source, the fraction of export supported by aeolian iron is larger (i.e., the aeolian points all lie above the 1 : 1 line by as much as 30 %).In other words, per source-injected dFe molecule, aeolian iron supports more export than the other iron types.Sedimentary and hydrothermal dFe make fractional contributions to export that are less than their fractional sources (the sedimentary and hydrothermal points lie below the 1 : 1 line by as much as ∼ 20 %).
The fact that the scatter plots of φ k versus σ k in Fig. 8 are reasonably compact suggests a simple underlying relationship.If we define a given source type's efficiency in supporting phosphorus export by P k ≡ φ P k / σ k , we see from Fig. 8  that P k varies with σ k .However, one might expect the efficiency of source type k relative to the efficiency of the other sources to be more constant: this relative export-support efficiency is controlled by dFe transport pathways and by scavenging, which is in turn controlled by the well-constrained organic-matter export.To investigate this possibility, we note that the export-support efficiency of the sources other than s k is given by P k ≡ (1 − φ P k )/(1 − σ k ) so that the exportsupport efficiency of source s k relative to the other sources is e P k = P k / P k .If e P k is constant, then it follows algebraically that φ P k = e P k σ k /[1+(e P k −1) σ k ]; i.e., the relationship between φ P k and σ k is determined by the single parameter e P k .Note that e P k is the slope of this theoretical φ k versus σ k relationship at the origin.Nonlinear least-squares fits of this functional form to the ( φ P k , σ k ) pairs of our family of states approximate the scatter plots well (lines in Fig. 8) and result in relative exportsupport efficiencies of e P A = 3.1 ± 0.8, e P S = 0.4 ± 0.2, and e P H = 0.3 ± 0.1, where the uncertainty for each source type is the standard deviation of the corresponding residuals.Thus, per source-injected molecule, aeolian iron supports 3.1 ± 0.8 times more phosphorus export than the other sources, while sedimentary and hydrothermal iron support 1/e P S = 2.3 ± 0.6 and 1/e P H = 4. ± 2. times less export than the other sources.The ability of aeolian iron to make disproportionately large contributions to supporting organic-matter export, quantified here by a relative export-support efficiency greater than unity, is presumably due to fresh aeolian iron being directly injected into the euphotic zone.The less-than-unity relative export-support efficiencies of sedimentary and hydrothermal iron reflect the fact that iron from benthic sources is generally subject to scavenging before it even reaches the euphotic zone.Because most large sedimentary sources are relatively shallow, a typical sedimentary dFe molecule will undergo less scavenging en route to the euphotic zone than a typical hydrothermal dFe molecule, which is quantified here by the higher relative export-support efficiency of sedimen- tary iron.These arguments are supported by the fact that, if we calculate φ P k only for the Southern Ocean, where the aeolian source is small and most aeolian iron is supplied as upwelled regenerated iron (Holzer et al., 2016), then all ( φ P k , σ k ) pairs lie closer to the 1 : 1 line.(In terms of relative export-support efficiencies, e P A is reduced to 2.0 ± 0.5, while e P S and e P H are increased to 0.5 ± 0.1 and 0.6 ± 0.2, respectively.)

Opal export
The opal export supported by each iron type can be calculated analogously, and the corresponding geographic patterns are shown in Fig. 9. Similar to the total opal export (Fig. 5b), the patterns of the opal export supported by each iron type emphasize regions with high diatom concentrations, namely the Southern Ocean and subpolar North Pacific and North Atlantic where there is upwelling and/or vertical mixing.
Figure 10.Percent global opal export supported by each iron type (aeolian, sedimentary, hydrothermal) versus the corresponding fractional source of that iron type.Lines represent fits to theoretical curves with fixed relative export-support efficiencies.(See text for details.) Aeolian-iron-supported opal export (Fig. 9a) is large in the Southern Ocean, but most pronounced in the subpolar North Pacific, where both diatom production is significant and aeolian input is high downwind from Asia's deserts.While tropical opal export is of secondary importance, the tropics are most pronounced for aeolian-supported export, again because of the direct source there.The pattern of sedimentaryiron-supported opal export (Fig. 9b) is broadly similar to that for aeolian dFe, but weaker in the tropics.The pattern of hydrothermal-iron-supported opal export (Fig. 9c) is dominated by the Southern Ocean, where diatom production is high and where most hydrothermal iron upwells.The patterns of iron-type-supported opal export have tightly clustered zonal means.
The amplitude of the opal-export patterns varies systematically with the iron source strength as summarized in Fig. 10, which shows the fractional iron-type-supported opal export, φ Si k ≡ Si k / Si , as a function of the corresponding fractional dFe source, σ k .While aeolian dFe is also the most efficient iron type for supporting opal export, aeolian dFe is less efficient for opal export than for phosphorus export (the aeolian points fall closer to the 1 : 1 line).Conversely, sedimentary iron is slightly more efficient in supporting opal export than in supporting phosphorus export, and hydrothermal dFe is only slightly less efficient than sedimentary dFe.
Similar to our preceding analysis of phosphorus export, we define the relative opal export-support efficiencies by Nonlinear least-squares fits give relative opal exportsupport efficiencies of e Si A = 2.3 ± 0.5, e Si S = 0.5 ± 0.2, and e Si H = 0.5 ± 0.2.Per source-injected molecule, aeolian iron is thus 2.3 ± 0.5 times more efficient in supporting opal export than the other sources, while sedimentary and hydrothermal iron are 1.9±0.5 and 2.±1.times less efficient, respectively, than the other sources.
The lower efficiency of aeolian iron for supporting opal export is consistent with the fact that opal export occurs primarily in the Southern Ocean, where direct aeolian input is small.Similarly, the greater efficiency of sedimentary and hydrothermal iron is consistent with the bulk of the opal export occurring in the upwelling regions of the Southern Ocean where access to deep iron sources is greatest.This is supported by the fact that plots of φ P k and φ Si k for the Southern Ocean only versus σ k (not shown) are both nearly identical to the plot of fractional global opal export of Fig. 10.(The relative Southern-Ocean-only opal export-support efficiencies are e Si A = 1.9 ± 0.5, e Si S = 0.5 ± 0.2, and e Si H = 0.7 ± 0.2.) 8 Discussion and caveats Our approach has a number of limitations that should be kept in mind.Most importantly, inverse-model estimates are only as good as the data used to constrain them.The dFe observations are too sparse in space and time to construct a gridded annual mean climatology like those available for PO 4 and Si(OH) 4 .We averaged the available dFe data to minimize observational biases, but in many places observations are only available for one time of year and likely contain seasonal biases.Other biases are likely introduced when dFe measurements alias episodic source events such mineral dust downwind from the major deserts (e.g., Croot et al., 2004;Johnson et al., 2010).In the near future, GEOTRACES will release an expanded data product that will include Pacific transects that were not available for the Intermediate Data Product 2014 used here.The additional dFe observations will help constrain the hydrothermal sources, particularly the strength of the Pacific source relative to that of the other basins.Important nonnutrient observational fields for our inverse model are the satellite-measured photosynthetically active radiation (PAR) and ocean-color-derived estimates of the size-partitioned phytoplankton concentrations.Small-scale features of the PAR field, e.g., in the Weddell Sea where ice and cloud cover play a role, are uncertain with the PAR for different time averages showing different features.The satellite-based estimates of phytoplankton concentrations also carry unquantified uncertainties due to a number of assumptions (Kostadinov et al., 2016).In our inverse model, these estimates provide constraints on how biological production is partitioned among the different functional classes.The unquantified uncertainties warrant re-evaluation as independent satellite-derived estimates become available in the future.
Most biogeochemical parameters are determined through objective optimization against available observations, but the construction of the cost function and the choice of which parameters are optimized and which are prescribed are necessarily subjective.For example, choosing a different set of B. Pasquier and M. Holzer: Inverse-model of coupled Fe-P-Si cycles weights (ω P , ω Fe , ω Si , ω plk ) to combine the four terms of the cost function would result in different optimal parameters.Similarly, assigning greater importance to dFe data measured as part of a vertical profile introduces another arbitrary weight.As for any nonlinear least-squares problem, it is also important to recognize that any minimum of the cost function found numerically is not guaranteed to be the global minimum and it is always possible that a better fit exists for a different set of parameters.By the same token, depending on the choice of initial state, the optimizer may find a local minimum that has grossly unrealistic features and must be rejected.
The uncertainty in key metrics (e.g., global phosphorus export) was quantified in terms of their spread across our family of state estimates and in terms of systematic variations with the iron source strengths.While our efficient numerics allow us to easily determine the linear sensitivities of any metric with respect to all parameters (from which one can also estimate uncertainty), we did not do so here because the spread in the metric across the family is more relevant.Given the large set of parameters x j and several interesting metrics M i , a detailed investigation of all the sensitivities ∂M i /∂x j evaluated at the optimal states is beyond the scope of this study.In principle, one can estimate the uncertainty of the optimal parameters themselves using a Bayesian framework (e.g., Teng et al., 2014).However, this requires the construction of suitable covariances and is also beyond the scope of this study.
A key limitation of our approach is that seasonality is ignored and we use a steady circulation.This circulation is constructed so that its transport reproduces the annual-mean observed temperature, salinity, CFC-11, radiocarbon, and PO 4 fields with minimal error.The circulation is hence not a simple average, but an effective ventilation-weighted mean.However, we acknowledge that effects due to the seasonal covariance of biological production and circulation cannot be captured.
Our model of the macronutrient cycles makes a number of simplifying approximations.We ignore external inputs of silicic acid and therefore also neglect permanent burial of opal in sediments.While this approximation has been shown to have a negligible impact on particle fluxes (Sarmiento et al., 2007), we acknowledge that our estimates will miss features such as silicic acid plumes due to crustal fluid venting (Johnson et al., 2006).The uncertainty of the silicon cycle that is most difficult to quantify stems from our simple parameterization of opal dissolution, which does not account for partial frustule protection by decaying organic material or the effect of digestion by zooplankton.Another key uncertainty lies in our parameterization of the Si : P uptake ratio, particularly its dependence on dFe.While our empirical formulation captures known dependencies qualitatively, a first-principle derivation based on cell biology is currently lacking.These remarks apply equally to the Fe : P uptake ratio.
Although our model of the iron cycle includes an explicit representation of the redissolution of scavenged iron, effects of subgrid topography, and dynamic coupling to the phosphorus and silicon cycles, and is thus much more complex and mechanistic than the iron model of Frants et al. (2016), it was still necessary to make a number of simplifying approximations.Specifically, we do not model ligands dynamically, we ignore colloidal iron (e.g., Fitzsimmons and Boyle, 2014), and do not represent some iron sources that may be locally important, such as input from icebergs (Klunder et al., 2011(Klunder et al., , 2014)).We also assume that PO 4 and dFe are remineralized with the same Martin curve and in the same ratio in which they were utilized.The recent work by Twining et al. (2014) suggests that sinking diatoms release phosphorus higher in the water column than iron, but we do not have sufficient information to model these effects.Given the large uncertainties in the external iron sources, the neglected details are likely of second order for estimating the large-scale dFe concentration.
Other uncertainties concern the phosphorus cycle to which the uptake of the other elements is keyed.While the optimized phosphate fields have the smallest misfit with observations, our model of the phosphorus cycle makes several simplifying approximations.The Martin exponent is assumed to be globally uniform although in reality it almost certainly varies spatially (Weber et al., 2016), potentially leading to underestimated gradients in our model.To avoid carrying an additional tracer, we approximated DOP to have zero lifetime.In reality, DOP has a wide range of lifetimes, and the lifetime of semi-labile DOP is typically assumed to be a fraction of a year (e.g., Primeau et al., 2013).However, the neglect of DOP is unlikely to seriously affect our estimates.DOP represents only a tiny fraction (less that 1 %) of the total phosphorus pool (e.g., Pasquier and Holzer, 2016), and by using a Martin exponent optimized for a restoring model without DOP, we were able to match PO 4 concentrations to within 5 % of the observations.We emphasize that the carbon export reported here was simply our estimate of the phosphorus export converted to carbon units.No effort was made to compute a more realistic carbon export such as could be achieved with an explicit representation of the carbon cycle (which would require additional tracers and was numerically too expensive) and the C : P export ratio was treated as globally uniform.While a globally uniform export ratio is acceptable for a unit conversion, the true C : P export ratio is now known to vary spatially (e.g., Galbraith and Martiny, 2015;Teng et al., 2014).Although using the regionally varying C : P ratios would be more realistic, we find that variations in C : P have only modest effects on the globally integrated carbon export: (i) applying the variable P : C relation of Galbraith and Martiny (2015) to the phosphorus export of our typical state gives a carbon export of 8.5 ± 0.4 Pg C yr −1 or 9.4 ± 0.9 Pg C yr −1 when we use their log-binned parameter values.(ii) Applying the regional C : P inverse-model estimates of Teng et al. (2014) gives a carbon export of 10. ± 2. Pg C yr −1 .Both this estimate and the one based on the log-binned regression agree within their uncertainties with our simple unitconversion value of 10.3 ± 0.4 Pg C yr −1 .

Summary and conclusions
We have formulated a steady-state model of the coupled phosphorus, silicon, and iron cycles that is embedded in a steady data-assimilated global circulation.The model is of intermediate complexity and couples the nutrient cycles through co-limitations on biological uptake and through the scavenging of iron by organic particles.The concentrations of the small, large, and diatom phytoplankton functional classes are calculated diagnostically, which avoids the need for plankton concentration tracers.We explicitly represent iron scavenging by POP, opal, and mineral-dust particles, and the redissolution of POP-and opal-scavenged iron.Subgrid topography is parameterized for the sedimentary iron sources and intercepts all vertical fluxes.The relative simplicity of the biogeochemical model and the matrix formulation of the steady-state advective-diffusive transport afford highly efficient numerics.Steady-state solutions are readily found using a Newton solver, which permits the model to be used in inverse mode to constrain many of the biogeochemical parameters through objective optimization.The optimization minimizes the mismatch with the observed nutrient concentrations and with satellite-derived estimates of phytoplankton concentrations.
Our estimates of the macronutrient concentrations closely match the observational WOA13 climatology with volumeweighted rms errors of 5 % for phosphate and 12 % for silicic acid relative to the global mean.The modeled dFe concentration has a larger cost-weighted rms mismatch of ∼ 45 % relative to the global mean.However, the cost-weighted basinaveraged vertical dFe profiles for the Atlantic and Southern Ocean generally lie within the observational uncertainties.The Pacific dFe profiles show systematic biases, in part because there are relatively few dFe observations for the Pacific, with no Pacific transect in the GEOTRACES Intermediate Data Product 2014.The fractional global biomass of each phytoplankton functional class lies within 7 % of the observation-based estimates by Kostadinov et al. (2009).
Given that even the order of magnitude of the iron sources is uncertain, we produced a family of state estimates with a wide range of iron source strengths.Because different iron source strengths are compensated for by optimally adjusting the scavenging parameters, each family member fits the observations with roughly the same fidelity.This means that the available observed dFe and phytoplankton concentrations by themselves are insufficient to constrain the sources.This conclusion can also be gleaned from the model intercomparison of Tagliabue et al. (2016) and was reached using an inverse model by Frants et al. (2016), who also considered a family of state estimates.However, while Frants et al. ( 2016) varied only the aeolian source, our estimates here explore a range of sedimentary and hydrothermal source strengths in addition to a much wider range of aeolian source strengths.
We partitioned the dFe concentration field into its aeolian, sedimentary, and hydrothermal components without perturbing the system using the approach of Holzer et al. (2016).While the individual source components vary widely depending on the source strengths, we find that the total dFe concentration is well constrained.Variations in the aeolian component are compensated for primarily by sedimentary dFe.Both the compensations between different iron types and between effective sources and sinks suggest that a more dense sampling of the ocean's dFe field by future measurement campaigns may not provide the information necessary for constraining the iron sources.The required information may ultimately have to come from better direct quantification of the source and/or scavenging processes themselves.
Nutrient limitation patterns were defined by jointly considering whether the PO 4 , Si(OH) 4 , and dFe concentrations fell below their half-saturation values for uptake.Iron limitation was thus deemed to occur where only dFe fell below its half-saturation value, phosphate-iron co-limitation occurred where both PO 4 and dFe fell below their half-saturation values, and so on.The resulting limitation patterns are robust across our family of state estimates and broadly consistent with direct observations (Moore et al., 2013) and with alternatively defined limitation patterns in the BEC model (Moore et al., 2004).The large and diatom functional classes show iron limitation in the Southern Ocean, eastern tropical Pacific, and subpolar North Pacific, with PO 4 -dFe and (for diatoms) PO 4 -Si(OH) 4 -dFe co-limitations in the Pacific and South Atlantic subtropical gyres.The Indian Ocean, tropical Atlantic, and North Atlantic are largely iron replete (i.e., not limited in the sense defined) with PO 4 limitation and for diatoms PO 4 -Si(OH) 4 co-limitation.
The export productions of phosphorus and opal are well constrained across our family of state estimates, in terms of both pattern and magnitude.Because we model three phytoplankton functional classes with distinct optimized uptake timescales, our phosphorus export (expressed in carbon units) of 9.5-11.Pg C yr −1 is ∼ 40 % larger than that estimated by Primeau et al. (2013) and closer in spatial pattern to the satellite-based estimates of Dunne et al. (2007).The opal export of 164.-177.Tmol Si yr −1 overlaps with the estimate of Holzer et al. (2014), who used a simple restoringtype model of the silicon cycle uncoupled from other nutrients.
We quantified the role of the iron cycle in shaping the phosphorus and opal export productions.We find that each iron source type (aeolian, sedimentary, hydrothermal) supports phosphorus and opal exports with a distinct geographic pattern that is robust across the family of state estimates.The export pattern supported by a given iron type reflects the nature of its source.Sedimentary and hydrothermal iron B. Pasquier and M. Holzer: Inverse-model of coupled Fe-P-Si cycles support phosphorus export that is dominantly shaped by the large-scale patterns of upwelling, which brings these iron types to the surface.Aeolian iron supports export that is shaped by both the pattern of direct aeolian input and by large-scale upwelling, which brings regenerated as well as scavenged and redissolved aeolian dFe back into the euphotic zone.For opal export, the signature of each iron type is qualitatively similar, but compared to phosphorus export, the opal export patterns tend to be weaker in the tropics and stronger at high latitudes, especially in the Southern Ocean where diatom concentrations and silicon trapping are strongest.
The fraction of the globally integrated export supported by a given iron type varies systematically with its fractional global source.These variations quantify the export-support efficiency of each iron type per source-injected molecule.Aeolian iron is most efficient and supports a fraction of the global export that is larger than its fractional source, while sedimentary and hydrothermal iron are less efficient, supporting fractions of the global export that are less than their fractional sources.This is because dFe from deeper sources is more likely to be scavenged en route to the euphotic zone.The relative export-support efficiency of each iron type is robust across our family of state estimates.Per source-injected molecule, aeolian iron supports 3.1 ± 0.8 times more phosphorus export and 2.3 ± 0.5 times more opal export than the other iron types.Conversely, sedimentary and hydrothermal iron are respectively 2.3 ± 0.6 and 4. ± 2. times less efficient in supporting phosphorus export and 1.9 ± 0.5 and 2. ± 1. times less efficient in supporting opal export than the other iron types.
Our optimized model is ideally suited for investigating the response of the global ocean ecosystem to a variety of biogeochemical perturbations.In the future, we will report on the model's response to perturbations in the iron supply and on a more comprehensive analysis of the detailed workings of the iron cycle.
Data availability.The temperature, phosphate, and silicic acid data used in this study are available from the World Ocean Atlas v2 2013 (www.nodc.noaa.gov/OC5/woa13/woa13data.html).The dissolved iron data used in this study, including the data set of Tagliabue et al. (2012), are available from GEOTRACES (www.bodc.ac.uk/geotraces/data).The satellite estimates of the concentrations of picophytoplankton, nanophytoplankton, and microphytoplankton are available from the PANGAEA data repository (doi:10.1594/PANGAEA.859005).The yearly irradiance data from NASA's MODIS Aqua PAR are available from the OceanColor website (https://oceancolor.gsfc.nasa.gov/).B. Pasquier and M. Holzer: Inverse-model of coupled Fe-P-Si cycles both data sets and remove dFe observations above 2.71 nM, which probably correspond to transient states with short timescales that cannot be captured by our steady-state model.In order to compensate the fact that most dFe observations are close to the surface, we give more weight to observations that are part of a "profile".(A dFe observation is deemed to belong to a "profile" if there are 10 or more observations at the same latitude and longitude, and if one of those was recorded deeper than 2000 m.)Because the dFe observations do not sample the seasonal cycle uniformly, we adopt an approach similar to Frants et al. (2016) to reduce potential sampling bias when we interpolate the data to our model grid: if multiple dFe observations lie in the same grid cell, we first take the seasonal averages, which we then average again to estimate the annual mean.
As in Eq. ( 22) for PO 4 and Si(OH) 4 , we use volume weights to evaluate the dFe concentration mismatch with the observations.However, because not all model cells contain dFe observations, we define a dFe-specific vector of grid box volumes, v all Fe , which has nonzero elements only for grid boxes that contain at least one dFe observation.We also define a dFe "profile-specific" vector, v pro Fe , which is nonzero only for grid boxes that contain "profile" observations.The corresponding weights are defined by where V all Fe is the total volume of grid cells which contain a dFe observation, and V pro Fe the total volume of grid cells containing "profile" observations.We define the total dFe weight vector, w Fe , for the mismatch with observations in Eq. ( 23 where we give extra weight to the "profile" observations.The 1 : 4 ratio was manually adjusted until "profile" observations were deemed to have sufficiently strong influence on the state estimates.We also tried different approaches for weighting the model-observation dFe mismatch, including the use of inverse variances (Frants et al., 2016), but found no significant difference in our results.
Appendix D: Optimization strategy details

D1 Prescribed parameters
The following considerations determined which parameters were not optimized and how their values were chosen.The recyclable fractions of POP and opal scavenging, f POP and f bSi , compensate with the maximum Fe : P uptake ratio, R Fe : P 0 , and thus were prescribed at 90 % (Moore and Braucher, 2008).(This compensation results from the biological iron pump having almost the same effect as the combination of scavenging and recycling.)Similarly, the detrital fractions, f 0 c , which set the particle export ratio, are directly compensated by all the other parameters in the uptake formulation.We therefore followed Dunne et al. (2005) and assigned their "small" detrital fraction to f 0 sml and their "large" detrital fraction to both f 0 lrg and f 0 dia .When trying to optimize the silicon half-saturation rate k Si dia , starting from a value of 1 mmol m −3 (e.g., Matsumoto et al., 2013), we found that the optimal value always remained within a few percent of this initial value.This is in part due to the fact that in regions of high diatom concentration the Monod term for silicic acid is near saturation so that there is little sensitivity to the precise value of k Si dia .Moreover, there appears to be consistency across the literature that k Si dia = 1 mmol m −3 .We therefore simply fixed k Si dia at this value for numerical efficiency.

D2 Choice of initial parameter values
We first chose an initial set of values for the remaining parameters as collected in Table 2.The parameters of the iron cycle were taken from of the typical state estimate of Frants et al. ( 2016) except for the half-saturation constant of the Fe : P ratio, which was taken from the work of Galbraith et al. (2010), and the scavenging-rate parameters.The initial parameters for POP and opal scavenging, κ POP scv and κ bSI scv , were determined so that the globally integrated scavenging of each process was initially ∼ 5 Gmol Fe yr −1 (the typical total source/sink strength reported by Frants et al., 2016).The initial value of the dust scavenging rate parameter, κ POP dst , was chosen so that the sink due to dust scavenging was ∼ 10 % of the total sink of the initial state.
The initial irradiance half-saturation constants were taken from the work of Doney et al. (2006).The initial uptake half-saturation constants k i c were taken from the work of Matsumoto et al. (2013).The uptake timescales τ c were set to an initial value of 6 days and optimized subject to the constraint τ sml ≥ τ lrg ≥ τ dia .The initial values of the maximum phytoplankton concentrations were calculated as p max c = p * /(λ 0 τ c ) using p * = 0.018 mmol P m −3 (Galbraith et al., 2010) and λ −1 0 = 5.26 d (Dunne et al., 2005).The initial values of the parameters of the Si : P ratio were set so that k Fe Si : P and k Si Si : P were on the order of typical dFe and Si(OH) 4 concentrations, while R Si 0 and R Si m were based on corresponding Si : N uptake ratios found in the literature and converted using N : P = 16 : 1.Thus, in terms of Si : N units, R Si 0 was chosen to be of the order of the minimum Si : N uptake ratio used by Matsumoto et al. (2013)   macronutrient availability, is the large-scale transport into the euphotic zone, particularly the transport into iron-limited regions such as the Southern Ocean.We have no reason to think that this large-scale transport is suspect as evidenced by realistic large-scale patterns of production that are robust across our family of states with widely varying iron source strengths.

Appendix H: Iron source and sink patterns
Figure H1 shows the vertically integrated sources of dFe with a logarithmic color scale.The aeolian soluble iron deposition pattern is identical to that of the study of Luo et al. (2008), albeit limited to the oceans.The tropical Atlantic close to the Sahara, the Arabian Sea, and the Bay of Bengal are the regions of largest aeolian iron deposition.The hydrothermal iron sources follow the mid-ocean ridges with the pattern of the OCMIP protocol, but are independently scaled for the Atlantic, Pacific, Indian, and Southern Ocean basins.Sedimentary iron is more intense where export production is large and in areas where oceans are shallower, because in both cases, a large flux of organic matter (POP in our model) reaches the sediment.The subgrid topography plays a significant role in the pattern of sedimentary iron, in particular for coastal regions and large underwater plateaus, e.g., near the Kerguelen Islands or the Falkland Islands.Because of unrealistic circulation features in the Sea of Japan, we zero all sources there consistently with zeroing out production in the Sea of Japan.
Figure H2 shows the vertically integrated sinks that balance the sources of Fig. H1.The scavenging due to sinking mineral dust particles is about 3 orders of magnitude smaller than the sink due to organic and opal particle scavenging and could be neglected without changing our estimates appreciably.Although the pattern of the scavenging sinks has significant local variations among our family of state estimates, the zonally averaged pattern (vertically integrated sink normalized by its global mean) is broadly similar across the family.
c (z e ) = 0 z e f c U c dz, and that the POP flux attenuates with depth according to the Martin power law φ POP c (z) = φ POP c (z e )(z/z e ) −b due to remineralization in the aphotic zone.The operator B P therefore injects PO 4 with the divergence of φ POP c into the aphotic water column.The flux into the ocean bottom is remineralized in the lowest grid box as in the work of Primeau et al. (2013).The exponent b was determined to be b = 0.82 using a restoring-type phosphate-only model.(Most parameters were optimized for the full coupled model -for details of our optimization strategy see Sect.3.3.) and M. Holzer: Inverse-model of coupled Fe-P-Si cycles et al. (2016), and L max , τ b , L v , and L b are optimizable parameters.

B.Figure 1 .
Figure 1.Total cost metric and rms mismatch of the nutrient concentrations as a function of the aeolian, hydrothermal, and sedimentary iron source strengths (σ A , σ S , σ H ) plotted for all our optimized state estimates.State estimates with a total nondimensional cost that exceeds 15.4 are indicated by black crosses and were excluded from our family.Plots on the left show state estimates for which σ H < 1 Gmol yr −1 , while for plots on the right σ H ≥ 1 Gmol yr −1 .(a) Square root of the total cost expressed as a nominal percentage representative of the mean rms mismatch of the nutrient and phytoplankton concentrations.(b) Rms mismatch of the PO 4 concentration as a percentage of the global mean PO 4 concentration.(c) As (b) but for Si(OH) 4 .(d) As (b) but for dFe.(e) The value of the hydrothermal source σ H for each family member.

Figure 2 .
Figure2.Joint distribution of the cost-weighted observed and modeled concentrations of PO 4 , Si(OH) 4 , and dFe.The percentiles of the cumulative distribution are defined such that x % of the distribution lies outside the x-percentile contour.Large percentiles thus correspond to high densities.For PO 4 and Si(OH) 4 , the WOA13 observations were interpolated to our model grid.The dFe observations were interpolated to our model grid from the data compilation ofTagliabue et al. (2012) and from the GEOTRACES Intermediate Data Product 2014(Mawji et al., 2015).

Figure 3 .
Figure3.Basin-wide, cost-weighted average profiles of the (red) observed and (grey) modeled dFe concentrations for the Atlantic and Pacific oceans (both north of 40 • S), and the Southern Ocean (south of 40 • S).The profiles of our typical state estimate are shown in black.The error bars represent the combined standard error associated with the spatial standard deviation from the basin-mean profile and the observational standard deviation for each grid box.These were added in weighted quadrature using the weights for dFe mismatch from our cost function.

B
. Pasquier and M. Holzer: Inverse-model of coupled Fe-P-Si cycles 5 Limiting nutrients

Figure 4 .
Figure 4.The patterns of limiting nutrients for each phytoplankton functional class.The color cube at the bottom right shows the eight possible limitation regimes of our inverse model: red corresponds to dFe limitation, blue to P limitation, and green to Si limitation.Cyan, yellow, and magenta correspond to co-limitations of P and Si, dFe and Si, and P and dFe, respectively.White corresponds to co-limitation of all three nutrients, while black indicates no limitation.(See text for the definitions of the deficiencies D P c , D Si c , and D Fe c of the cube axes.)

Figure 5 .
Figure5.Local export production for each nutrient (maps on the left) and its zonal integral (curves on the right).Maps are shown for our typical state estimate, while we plot the zonal integral of each family member (scaled for dFe) in grey and the typical state estimate in black.(a) Phosphorus export, expressed in carbon units using C : P = 106 : 1.The blue zonal integral is the export production estimate ofPrimeau et al. (2013).(b) Opal export, where the blue zonal integral is the estimate ofHolzer et al. (2014).(c) Iron export, with its zonal integrals expressed as a percentage of the global iron export.The globally integrated exports for the typical estimate are indicated in the plot titles together with their ranges across the family of estimates in parentheses.

Figure 6 .
Figure6.Estimates of the dFe concentration in each basin (ATL, PAC, IND) and globally (GBL).The zonal averages in latitude-depth space on the left show the total dFe field of our typical state estimate.The corresponding horizontally averaged profiles of total dFe are shown in grey for each family member and in black for the typical state estimate.The three columns of plots on the right show the source-partitioned dFe profiles of each family member.The individual source-partitioned profiles are color coded according to the percent contribution of the aeolian source to the total iron source, with our typical state estimate shown in black.

Figure 7 .
Figure 7. Phosphorus export supported by each iron type [aeolian (a), sedimentary (b), hydrothermal (c)] normalized by its global mean.The maps show our typical state estimate, while zonal averages of the normalized phosphorus export are shown for each family member in grey, with the typical state estimate in black.

Figure 9 .
Figure 9. Opal export supported by each iron type [aeolian (a), sedimentary (b), hydrothermal (c)] normalized by its global mean.The map is for our typical state estimate, while zonal averages of the normalized opal export are shown for each family member in grey, with the typical state estimate in black.

Figure F1 .
Figure F1.Local export production (maps on left) and its zonal integral (curves on the right) expressed in carbon units (using C : P = 106 : 1).Maps are shown for our typical state estimate, while we plot the zonal integral of each family member in grey and the typical state estimate in black.The export productions are plotted for each phytoplankton functional class: small (top plots, a), large (middle plots, b), and diatom (bottom plots, c).Note the different color scale for the small class.The globally integrated exports for the typical estimate are indicated in the plot titles together with their ranges across the family of estimates in parentheses.

Figure H1 .Figure H2 .
Figure G1.Dissolved iron concentrations of the typical state estimate (contours) compared to the GEOTRACES Intermediate Data Product 2014 (dots).The abscissa runs south to north or west to east along the transects (map).

Table 1 .
Parameters that were prescribed from the literature or that were separately optimized in a submodel.

Table 2 .
Optimized parameters and range across family of state estimates. .
and R Si m was