Towards accounting for dissolved iron speciation in global ocean models

Abstract. The trace metal iron (Fe) is now routinely included in state-of-the-art ocean general circulation and biogeochemistry models (OGCBMs) because of its key role as a limiting nutrient in regions of the world ocean important for carbon cycling and air-sea CO2 exchange. However, the complexities of the seawater Fe cycle, which impact its speciation and bioavailability, are simplified in such OGCBMs due to gaps in understanding and to avoid high computational costs. In a similar fashion to inorganic carbon speciation, we outline a means by which the complex speciation of Fe can be included in global OGCBMs in a reasonably cost-effective manner. We construct an Fe speciation model based on hypothesised relationships between rate constants and environmental variables (temperature, light, oxygen, pH, salinity) and assumptions regarding the binding strengths of Fe complexing organic ligands and test hypotheses regarding their distributions. As a result, we find that the global distribution of different Fe species is tightly controlled by spatio-temporal environmental variability and the distribution of Fe binding ligands. Impacts on bioavailable Fe are highly sensitive to assumptions regarding which Fe species are bioavailable and how those species vary in space and time. When forced by representations of future ocean circulation and climate we find large changes to the speciation of Fe governed by pH mediated changes to redox kinetics. We speculate that these changes may exert selective pressure on phytoplankton Fe uptake strategies in the future ocean. In future work, more information on the sources and sinks of ocean Fe ligands, their bioavailability, the cycling of colloidal Fe species and kinetics of Fe-surface coordination reactions would be invaluable. We hope our modeling approach can provide a means by which new observations of Fe speciation can be tested against hypotheses of the processes present in governing the ocean Fe cycle in an integrated sense

more organic ligands (Wu et al., 2001). Using electrochemical techniques it has been shown that over large parts of the world ocean >99 % of dFe is actually complexed to organic ligands of typically unknown provenance (e.g., Gledhill and van den Berg, 1994;Van den Berg, 1995;Rue and Bruland, 1995;Wu and Luther, 1995;Boye et al., 2003;. This is important as it prevents the precipitation/scavenging of free inorganic Fe(III) to solid forms, which are effectively lost from the dissolved pool and bioavailable Fe species. Fe ligands thus effectively increase both the solubility and residence time of dFe in the ocean, as well as exerting a control on its bioavailability.
It has been shown that phytoplankton can access organically bound Fe (e.g. Maldonado and Price, 1999), but not all forms of complexed Fe are equally bioavailable to different groups of phytoplankton (Hutchins et al. 1999). Mechanisms for the uptake of organically bound Fe typically involves the reduction of the organically bound Fe by either specific reductases at the cell surface (Maldonado and Price 2001;Salmon et al., 2006;Maldonado et al., 2006), which is a mechanism well known from terrestrial plants (Moog and Brüggemann, 1994), or by excretion of reactive oxygen species (Shaked et al., 2005). Due to the diffusive loss of reduced Fe away from the cell, each of these mechanisms inevitably leads to some loss of Fe (Völker and Wolf-Gladrow 1999). Alternatively, phytoplankton can assimilate organically complexed Fe via specific uptake mechanisms (Boukhalfa and Price, 2002), which do not make the whole Fe pool bioavailable. All these mechanisms to access organically complexed Fe seem costly compared to the uptake of inorganic Fe species, which can be assimilated relatively straightforwardly (Hudson and Morel 1990;Morel et al., 2008) and means inorganic Fe can be thought to be the most bioavailable Fe fraction. Organically bound, as well as inorganic colloidal, Fe(III) can also be photoreduced in the presence of light to produce Fe(II) (e.g., Barbeau et al., 2003;Croot et al., 2008). As such, the speciation, residence time and bioavailability of Fe in the ocean depend on a suite of processes that are themselves highly sensitive to the environmental conditions of the ocean.
In the context of its complex speciation and cycling, dFe is treated very simply in "state-of-the-art" OGCBMs, with only a single dFe pool represented and ligand complexation accounted for assuming a single ligand of uniform concentration (e.g., Parekh et al., 2004;Aumont and Bopp, 2006;Moore and Braucher, 2008;Galbraith et al., 2010). Spatiotemporal variability in Fe speciation, cycling and bioavailability is therefore ignored. Alongside the lack of constraints from observations, this is mostly due to the prohibitive computational cost of simulating rapid Fe cycle reactions at the global scale. Three-dimensional regional models have modeled the Fe cycle in a prognostic fashion for Fe-limited waters and noted the potential role of environmental variability in governing the supply of Fe to phytoplankton (Tagliabue and Arrigo, 2006). Similar models have also been employed in a one-dimensional framework at time series sites in the subtropical and tropical Atlantic Ocean (Weber et al., 2005(Weber et al., , 2007Ye et al., 2009). Recently, Tagliabue et al. (2009) included the first order impact of light and temperature on Fe speciation in a 3-D OGCBM and suggested that the role of environmental variability in Fe speciation could be important in governing the residence time and bioavailability of dFe in the ocean.
Over the coming century, the ocean is predicted to undergo a great deal of environmental change, especially the Fe-limited Southern Ocean. It is likely that temperatures will rise, stratification will increase, light levels will increase, pH will fall (due to the uptake of anthropogenic CO 2 ) and reduced sea ice will extend the growing season. All of these changes might impact upon the speciation of Fe and some experimental evidence from mesocosm experiments indeed suggests "acidification" induces changes to Fe(II) levels (Breitbarth et al., 2010), while laboratory experiments using synthetic ligands lead to modifications to Fe bioavailability that depend upon the type of chelator considered (Shi et al., 2010). As it stands, even including the first order impact of light and temperature on the marine Fe cycle (as per Tagliabue et al., 2009) will not resolve the matrix of parallel changes resulting from climate change that will impact Fe cycle rate processes (e.g., oxidation rates, photoreduction rates), the dFe concentration itself and the concentration of ligands. To address these questions we require a tool that can resolve the speciation of Fe in a semi-prognostic manner at the global scale in a 'cost effective' manner. In this study, we outline a new approach that permits the 'semi-prognostic' modeling of Fe speciation at the global scale using an analytical approach similar to that typically employed for inorganic carbon speciation. We then use this model to speculate how Fe speciation might respond to the climate associated with an atmospheric CO 2 concentration of ∼1000 ppm as an illustration of how our Fe speciation model can be applied.

Theoretical framework
Our approach rests on the concept of "fast" and "slow" Fe cycle reactions that assumes, similar to modules that compute inorganic carbon speciation in OGCBMs, that there are a subset of "fast" Fe speciation reactions that approach equilibrium within the time step of the model (normally around one to two hours). For example, the chemical reactions that govern dFe speciation (oxidation, photoreduction, the formation and dissociation of Fe-ligand complexes etc.) are assumed to be "fast" reactions. Examples of "slow" reactions that would need to be computed prognostically by the OGCBM include scavenging of free inorganic Fe(III) onto particles or uptake and recycling of Fe by biology. We assume dFe speciation to be a 'fast' problem and therefore well suited to similar analytical approaches as have been successfully employed in Biogeosciences, 8, 3025-3039, 2011 www.biogeosciences.net/8/3025/2011/ OGCBMs that seek to compute inorganic carbon speciation for air-sea CO 2 exchange or pH calculations. The full equations of the dFe model used here and their sensitivity to various parameter modifications are presented in Tagliabue and Arrigo (2006) and Tagliabue et al. (2009). The state variables of the model are the free concentrations of Fe(II) (Fe(II) ), Fe(III) (Fe(III) ), Fe(III) bound to the weak non-bioavailable ligand (FeL W ) and the strong bioavailable ligand (FeL S ), solid Fe(III) (Fe p , which could be thought to include colloidal Fe(III)), the total dFe concentration (Fe T ), the uncomplexed weak (L W ) and strong (L S ) ligands and the total concentration of L W (L WT ) and L S (L ST ). We disregard here complexation of Fe(II) by organic ligands. Ferrous iron complexes have been suggested to be possibly responsible for the long residence time of Fe(II) in the SOIREE iron fertilization experiment (Croot et al., 2001), but have only been demonstrated in riverine or coastal waters with high fulvic acid concentrations Sulzberger, 1996 andRose andWaite, 2003). However, such Fe(II) specific ligands may be difficult to identify using current techniques when they are at low abundance (e.g., Croot et al., 2007Croot et al., , 2008. Rate constants required by the model are the oxidation of Fe(II) (k ox , which is a function of temperature, pH, salinity and oxygen concentrations), photoreduction of FeL W (k phW , which is a function of irradiance as per Tagliabue et al., 2009) and FeL S (k phS ), the formation of FeL W (k lW ) and FeL S (k lS ), the dissociation of FeL W (k bW ) and FeL S (k bS ), the precipitation of Fe(III) to Fe P (k pcp ) and the remineralization of Fe P (k r ). Re-arranging the differential equations for the Fe species results in the following four governing equations: Additional constraints are that the concentrations of Fe T , L WT and L ST must be conserved over the fast timescale: In order to solve the model analytically first requires a rearrangement of Eqs. (3 and 4) to yield: and These equations are then inserted into Eq. (5) to result in: where a = 1 + k pcp /k r , b = 1 + k phW /k ox , and c = 1+k phS /k ox . From Eqs. (6 and 7) it follows that the free ligand concentrations are: Equation (12) can now be used in combination with Eq.
(2) to produce: and solved for FeL S : where K S = (k bS + k phS )/k lS . Equation (14) is then combined with Eq. (10) to solve for the concentration of FeL W : Inserting Eq. (15) into Eq. (1) permits us to obtain an equation for Fe(III) which, after simplification and sorting into powers of Fe(III) , yields a third order polynomial solution for the concentration of Fe(III) : where K W = (k bW + k phW )/k lW . Equation (16) can be solved analytically or iteratively and has three solutions, but only one is positive and thus a realizable Fe(III) concentration. Therefore by first solving Eq. (16) for the Fe(III) concentration, one can then proceed to solve for the Fe P concentration (Eq. 9), FeL S concentration (Eq. 14), FeL W concentration (Eq. 15), and finally the Fe(II) concentration (Eq. 8). Thus for a given set of rate constants, which are either fixed or vary as a function of environmental variables, and the concentrations of Fe T , L WT , L ST , the procedure outlined above analytically solves for the concentrations of the 5 Fe species (Fe(II) , Fe(III) , FeL W , FeL S , and Fe P ) at considerably less computational expense than a prognostic solution.

Modeling framework and experiments
We decided to include our analytical solution for dFe speciation within the PISCES OGCBM ( 2006), since this model has been widely used for ocean biogeochemistry and climate applications, including some addressing Fe speciation (e.g., Tagliabue et al., 2009). Firstly, the analytical solution of Eq. (16) was solved iteratively in each grid cell of the model at each time step to yield the Fe(III) concentration. From this, the concentrations of all other Fe species can then be computed. The analytical solution uses properties that are either provided by the PISCES model (Fe T , L WT , L ST ) or computed in each grid cell and at each time step from variables simulated by the PISCES model (k ox , k phW , k phS , k lW , k lS , k bW , k bS , k pcp and k r ). For example, k ox will vary in space and time as a function of the temperature, pH, oxygen concentration, and salinity, following the equation of Santana-Casiano et al. (2005), while k phW and k phS will vary with depth and season following available radiation each particular model grid cell (all other rate constants are initially fixed in space and time). The total dFe pool is also modified each time step by phytoplankton uptake and remineralisation and all other source -sink terms for dFe traditionally included in the PISCES model (see: Aumont and Bopp, 2006 for a full list of Fe equations). We find that the calculated Fe T computed from the sum of all species calculated analytically is generally less than ±1 % in error relative to the dFe tracer prognostically simulated by PISCES (which is an input to the speciation solution) and only reaches a maximum of ±5 % error in a few isolated grid cells (below 75 m and 150 m the error is less than ±1 % and ±0.1 %, respectively). This demonstrates that our procedure has an acceptable error in calculating Fe speciation, especially considering the global nature of its application and the necessity to retain a degree of computational efficiency.

Rate constants
Values for the rate constants are taken from the published literature and, apart from the examples detailed here, are identical to those described by Tagliabue et al. (2009). For this study, we used the k ox equation (s −1 ) as described by Santana-Casiano et al. (2005, personal communication), which is a function of temperature, salinity and pH: log 10 k ox = 35.407 − 6.7109 pH + 0.5342 pH 2 Where Tk is the temperature in • K, pH is the pH (free scale) and S is salinity. The realised rate of Fe(II) oxidation (k ox ) is then modified by the oxygen (mol L −1 ) concentration (Santana-Casiano and Gonzalez-Davila, personal communication, 2010) using The kinetic characteristics of L W are assumed to be similar to Phaeophytin-type ligands and rate constants are taken from Witter et al. (2000), with a log conditional stability (log(k lW /k bW ) of 11.00 M −1 . L S is assumed to have the kinetic characteristics of dessferroxamine B-type ligands (Witter et al., 2000) with a log conditional stability of 12.12 M −1 .
In the absence of other information, the kinetic characteristics of L W and L S are fixed in space and time, and are within the range of measurements made in situ for "strong" and "weak" ligands (e.g., Rue and Bruland, 1995;Boye et al., 2003Boye et al., , 2006Cullen et al., 2006). Initially we define "bioavailable" dFe (bFe) as the sum of Fe(II) , Fe(III) and FeL S .

Parameterisation of Fe binding ligands
Most OGCBMs assume that the concentration of dFe binding ligands is fixed at between 0.6 and 1 nM and only assume one fully bioavailable ligand is present (e.g., Aumont and Bopp, 2006;Moore and Braucher, 2009;Tagliabue et al., 2009;Galbraith et al., 2010). Nevertheless, there is ample experimental evidence of at least two ligand classes and highly variable concentrations (e.g., Buck andBruland, 2007 Hunter and. While parameterizing the sources and sinks of two ligand classes in an OGCBM is perhaps out of reach at this moment (it has been done for a one-dimensional model, Ye et al., 2009), there is some data showing an relationship between ligands and dissolved organic carbon (DOC) concentrations (Wagener et al., 2008, see also: Hiemstra and van Riemsdijk, 2006). We therefore decided to use the relationship from the observations of Wagener et al. (2008) to permit us to have ligand concentrations that vary as a function of total DOC concentrations (DOC TOT , in µmol L −1 ) that are already prognostically simulated by PISCES).
PISCES includes a semi-labile DOC pool as a prognostic tracer (Aumont et al., 2001) and we therefore assume a constant refractory DOC pool of 40 µM to arrive at a total DOC concentration (DOC TOT ). Sources of DOC (and thus sources of ligands) in our model are exudation during photosythesis, zooplankton grazing, disaggregation of particles etc, with DOC lost due to bacterial activity and aggregation. Observations show ligand concentration minima of 0.4 nM at 40 µM DOC (Wagener et al., 2008). So following the philosophy of Hunter and Boyd (2007), we set the minimum L W concentration to be 0.4 nM at 40 µM DOC (representing a "refractory" weak ligand pool), and for DOC concentrations greater than 40 µM we portioned two thirds of the "extra" ligand into L S and one third into L W . In this fashion we account for the active production of strong ligands (L S ) by the euphotic zone biotic community, as well as subsurface production of weak ligands (L w ) in a relatively simple fashion. The spatial distribution of L T using our DOC-linked parameterization is shown in Fig. 1

Model experiments
We decided to use our Fe speciation OGCBM to conduct some illustrative model experiments. We firstly simulated Fe chemistry for the "present" climate using an atmospheric CO 2 level of 368.87 ppm (corresponding to observations from the year 2000), an ocean circulation from NEMO that arises from atmospheric re-analysis products (Aumont et al., 2008) and depart from a simulation conducted from 1860-2000 forced by atmospheric CO 2 observations (to ensure a correct ocean pH). Initially, we used the parameterisation of the Fe cycle as described as above, but we also conducted some illustrative sensitivity tests to examine assumptions regarding the nature of the variability associated with the Fe binding ligand pool. Finally, in order to appraise the possible impact of climate change on Fe speciation, we used 2 representations of ocean circulation, as well as initialization files for ocean biogeochemistry (to include the requisite DIC and pH changes), from the IPSL-CM5 coupled model at atmospheric CO 2 levels of 298.06 ppm and 1086.64 pmm (from a transient coupled simulation from pre-industrial CO 2 levels to 4 × CO 2 ) and conducted 10 year simulations with the Fe speciation analytical solution included in PISCES. Our model experiments are summarized in Table 1. In calculating the bFe concentrations in the following, we test two different assumptions, the first assumes that phytoplankton are not relient on free inorganic species, and thus "bFe" also includes organically complexed Fe. The other assumption assumes the contrary, that only free inorganic Fe(II) and Fe(III) are present as "bioavailable" species and this bFe = Fe(II) + Fe(III) . In a sense, this approach allows us to appraise the degree of Fe limitation experienced by only relying on free inorganic Fe species and not making the cellular investment to assimilate organically complexed Fe directly.  Fig. 2a and b) and over most of the ocean range between 0 and 100 pM. The annual mean (seasonal variability in fFe(II) is discussed below) proportion of the dFe pool present as Fe(II) (fFe(II)) ranges from 0 to around 30 % and is maximal at high latitudes (10-30 %), moderate in upwelling regions (3-4 %) and very low in the tropical oceans (<1 %) (Fig. 2c. For example, fFe(II) increases as one moves south in the Southern Ocean from <5 % near South Africa to ∼30 % at around 55 • S, with similar degree of change in the high latitude northern Oceans (Fig. 2c). In general, the latitude-longitude variability in fFe(II) is tightly linked to the variability in k ox for the surface ocean (Fig. 2d). Irradiance governs the depth distribution of fFe(II), with Fe(II) only making up an appreciable fraction of dFe at depths shallower than ∼100 m (Fig. 3a). An exception to this are suboxic zones, wherein the reduction in k ox results in Fe(II) levels >100 pM (Fig. 3b, 200-300 m). Organically complexed Fe(III) (FeL S and FeL W ) makes up almost 100 % of the dFe pool over most of the global ocean, declining slightly to around 85 % in polar waters where Fe(II) is greater due to supply from photoreduction and reduced oxidation rates.

Biovailable Fe
bFe shows variability that is linked to photochemistry, organic complexation and irradiance, as well as being highly sensitive to which species are assumed to be bioavailable. If bFe is assumed to encompass Fe(II), Fe(III) and Fe(III)L S then the proportion of the dFe pool present as bFe (fbFe) varies between 50 to 90 % in surface waters (when annually averaged, Fig. 4a). Variability in fbFe at the surface is positively related to irradiance (due to greater photoproduction of Fe ) and the total ligand concentration (due to reduced losses as Fe P ) and is negatively related to the dFe concentration (due to over-saturation of ligands and loss as Fe P ) (Fig. 4b, c and d). fbFe declines with depth due to the reduced irradiance and lower ligand concentrations at depth (at least in our DOC-based ligand parameterization). If we assume that only Fe(II) and Fe(III) are assumed to make up bFe, then the reduction in fbFe is striking (Fig. 5). Only at high latitudes, where Fe(II) is greatest (Fig. 2a), can fbFe approach even 25 % of the dFe pool and over large parts of the ocean fbFe is <5 % of dFe (Fig. 5). This is due to the reduced residence  time for Fe (the sum of Fe(II) and Fe(III)) away from polar waters and would imply that phytoplankton reliant on Fe would be chronically Fe limited in these waters (Tagliabue et al., 2009). Accordingly, fbFe is strongly and negatively related to spatial variability in Fe(II) oxidation rates when bFe = Fe(II)+Fe(III), and thus generally tracks variability in fFe(II).

Importance of seasonality
Seasonality plays an important role in Fe speciation, especially at high latitudes where there are large changes in environmental variables (temperature, irradiance etc, Tagliabue and Arrigo, 2006). During the winter-spring transition in the high latitude northern and southern hemispheres we suggest an increase in fFe(II) that is maximal in October-December in the Southern Ocean and June-July in the North Atlantic and sub-Arctic Pacific (Fig. 6a). Similarly, fbFe also increases from winter to spring as mixed layers shallow and irradiance levels are increased (Fig. 6b). For the Fe-limited Southern Ocean, fbFe increases from ∼50 % in winter (due to sea-ice or very deep winter mixed layers) to ∼80 % by spring when waters are ice-free and characterized by welllit stratified surface waters. Parallel to the increasing Fe(II), the organically complexed fraction of dFe declines between winter and spring.

Comparison with observations
The most obvious and widespread dataset with which to compare the model is the simulated dFe concentration (from the sum of the Fe species computed by our speciation model). This compares well to a new database of ∼13 000 dFe measurements (Tagliabue et al., 2011, R = 0.52 and 0.54 for the entire water column and 0-50 m, respectively), but although this shows our speciation model does not give unrealistic dFe concentrations in general, the good statistical reproduction of dFe (that arises from our speciation model) probably more reflects the successful simulation of dFe in PISCES. Fe speciation measurements are obviously rarer than those for dFe, but one candidate to compare our speciation model to is Fe(II). At high latitudes, the model appears to do a reasonable job, with surface Fe(II) concentrations of ∼25 pM south of New Zealand comparing well with a range of 19-46 pM from Croot et al. (2007) and modeled values of ∼30-60 pM from the western sub-Arctic Pacific within the range of ∼20-40 pM from Roy et al. (2008), with Fe(II) observations making up as much as 50 % of the dFe pool (modeled values are 30-40 %). In addition, the depth profile from Roy et al. (2008) is also relatively well reproduced by the model, except for the reduced attenuation of Fe(II) with depth in the model (Fig. 7). Earlier Southern Ocean observations of 0-45 pM using a towed fish (Bowie et al., 2002) are also reasonably well reproduced by the model. . It is noteworthy, that the observed latitudinal trends in both Fe(II) and fFe(II) were only significant for daytime stations (Fe(II)) and for both daytime and all stations (fFe(II)), but never when only night-time stations were considered (Sarthou et al., 2011). In the eastern North Atlantic, the onshore-offshore trend (from >250 to 100-150 pM) observed by Boye et al. (2003) is similar to that in the model and observed offshore values (∼100 pM) are, in general, only slightly underestimated by the model (although the limit of detection in Boye et al. (2003) was 100 pM). This is probably due to the onshore-offshore trend in dFe concentrations, since we do not include a specific source of Fe(II) at the margin (there is a margin source of dFe in PISCES). High modeled Fe(II) levels in the Baltic Sea agree well with measurements of Breitbarth et al. (2009). In general, this suggests that the model is reasonably well reproducing the dominant processes governing the Fe(II) distribution in higher latitude Atlantic, Pacific and Southern Oceans. It is worth noting that Fe(II) as a transient species shows significant variability on short time and space scales (e.g., Croot et al., 2007;Roy et al., 2009) and it is for this reason that we evaluated the general trends in the modelled Fe(II) concentrations rather than comparing the model and observations "point by point" (where a model error in mixed layer depth, temperature, light etc could drive a deviation from observations that is not due to the speciation model per se).
The model does a poorer job when compared to the comprehensive lower latitude Pacific Ocean measurements of Hansard et al. (2009) along ∼30 • N (line P02) and ∼152 • W (line P16N), with observed values of >30 pM (as high as >100 pM in some places) greatly underestimated by the model (generally <5 pM). This could be due to either differences in methods to other studies, errors in the modeled dFe field (which is closely linked to absolute Fe(II) concentrations, Fig. 2b), processes missing in our speciation model,  and 100 pM and are therefore very low, relative to the Fe(II) concentrations measured by Hansard et al. (2009), which are of the same order. This suggests it would be very difficult to achieve any appreciable Fe(II) in this region whilst modeled dFe values remain so low. Additional sensitivity tests focused on producing more Fe(II) in this region (drastically reducing k ox or increasing photoreduction, not shown) do not permit any appreciable accumulation of Fe(II) therein and result in unrealistically high Fe(II) concentrations in the high latitudes. Fe(II) might also be underestimated in the lower latitude ocean because the diurnal cycle in irradiance is not included in our OGCBM. Another important issue to bear in mind is that we compare point measurements to the monthly mean model output. Models and observations (e.g., Bowie et al., 2002;Tagliabue and Arrigo, 2006;Croot et al., 2007Croot et al., , 2008Roy et al., 2008;Breitbarth et al., 2009;Ye et al., 2009;Sarthou et al., 2011) show a high degree of variability in Fe(II) in response to changing environmental conditions (especially solar radiation on the diel cycle), although it is noticeable that Hansard et al. (2009) note no diurnal cycle in their observations, in contrast to other studies. Our model may not capture these "extreme" events that are highly specific to the time and location of each precise sample. This is because despite a 1.5 h timestep, we do not include the diurnal cycle and compare monthly mean modeled Fe(II) to point measurements. Therefore, we must conclude that a combination of a very low modeled dFe concentration, lack of high frequency variability and perhaps also an impact of acidified samples and errors in our formulated Fe cycle (see below) precludes a good reproduction of the reported Fe(II) data in the subtropical Pacific Ocean (Fe(II) levels are higher in the tropical Atlantic due to greater dust input of dFe). However we do note the increased Fe(II) in suboxic zones (Fig. 3b), in the eastern tropical Pacific in particular, that compare well to measured increases in Fe(II) at low oxygen levels (e.g., Hopkinson and Barbeau, 2007). As regards the degree of organic complexation, our results of virtually 100 % complexation of dFe agrees with all available observations (e.g., Boye et al., 2003;Boye et al., 2006;Buck and Bruland, 2007) and lesser complexation where Fe inputs are high is in accord with the findings from an artificial Fe enrichment experiment (Boye et al., 2005).
Overall, our speciation model can be seen to do a much better job at higher latitudes, rather than lower latitudes without significant Fe inputs (where Fe(II) levels appear too low).

Sensitivity tests
If we assume that Fe binding ligands are fixed in space and time, then we find that Fe speciation and cycling is modified. For example, fixing ligands at 0.6 nM results in a lower ligand concentration over much of the ocean than from our DOC-based parameterization (Fig. 1). A lower ligand concentration unsurprisingly results in a reduction in the proportion of the total dFe pool that is organically complexed (Fig. 8a). This then impacts the total dFe pool, which is reduced (due to greater losses as Fe P ), especially in regions of high Fe inputs (beneath zones of dust deposition and near coasts). The impact of fixed (and generally lower) ligand concentrations upon bFe depends on the assumed make up of the bFe pool. If bFe = Fe(II)+Fe(III)+FeL S , then bFe declines if ligands are fixed (Fig. 8b), due to reduced stabilization of bFe by L S (as its concentration is reduced). On the other hand, if bFe is assumed to be only made up of Fe(II) and Fe(III), then assuming lower and fixed ligand concentrations actually increases bFe (albeit from very low levels, Fig. 8c), particularly in areas of high Fe input, due to lesser complexation by organic ligands (which are assumed inaccessible to phytoplankton in this formulation of bFe). Thus the nature of the ocean ligand pool has impacts upon the general speciation of Fe, as well as its residence time and bioavailability. We reiterate that for phytoplankton that can access organically complexed Fe, bFe declines when ligands are fixed (at generally lower levels than measured), while for phytoplankton reliant on inorganic Fe, bFe increases when ligands are fixed since more Fe is in inorganic forms.

Fe speciation at four times CO 2
At atmospheric CO 2 levels of approximately 1000 ppm the environmental properties of the ocean are unsurprisingly greatly modified. In general, and similar to previous studies with fully coupled climate-OGCBMs (e.g., Steinacher et al., 2010), the surface ocean is warmer, more stratified (reduced mixed layer depth) and has a lower pH. In the Southern Ocean, sea ice coverage is also reduced which lengthens the growing season. Our objective here is not to comprehensively analyze these aspects (this is for other more focused papers), but to examine how Fe speciation changes using our analytical approach. Turning firstly to Fe(II), we find large increases in fFe(II) (chosen to remove the effect of climate on absolute dFe concentrations) due to climate change that are maximal in the high latitude oceans (Figure). fFe(II) increases by as much as >40 % in the high latitude Southern and Northern Oceans (Fig. 9a) and must be responding to changes to oxidation and photoreduction rates. However, a closer inspection reveals that the largest changes in fFe(II) occur where photoreduction rates were not significantly changed and that there is a very close relationship between the predicted changes to fFe(II) and oxidation rates (Fig. 9b). This suggests that the impact of reduced pH is overriding the impact of greater temperature to yield a net reduction in the future oxidation rate of Fe(II), especially in the high latitude oceans. Similar accumulations of Fe(II) at lower pH were obtained in mesocosm experiments by Breitbarth et al. (2010), where pH changes impacting oxidation rates were also found to be the dominant effect.
Modifications to the Fe(II) concentrations due to high CO 2 induced changes to oxidation rates have implications for the speciation and bioavailability of Fe. Firstly, the greater proportion of the dFe pool present as Fe(II) reduces the amount of Fe that is complexed by organic ligands, but since this Fe is instead retained as Fe(II) species (rather than Fe(III)), there is not a great impact on losses of dFe as Fe P . As we noted for the modern climate, the impacts on fbFe depend upon the assumptions regarding the bFe pool. If bFe = Fe(II) + Fe(III) + FeL S (i.e., organically complexed Fe is available), then climate change and ocean acidification have only a modest impact on fbFe, with fbFe increasing by <10 % in the Fe-limited Southern Ocean or by up to 20 % in the Arctic (Fig. 9c). This is because we assume all bFe (Fe(II), Fe(III), and FeL S ) species to be similarly bioavailable. On the other hand, if only Fe(II) and Fe(III) are assumed bioavailable (i.e., organically complexed Fe is not available), then large increases in fbFe, which parallel those in fFe(II), are found in the Fe-limited Southern Ocean (Fig. 9d) as more dFe now remains in the Fe(II) state, relative to organically complexed Fe(III) pool, due to the reduced oxidation rates. An implication of this result is that climate/acidification induced changes to Fe speciation (that primarily result from pH changes) might make phytoplankton that rely only on Fe(II) and Fe(III) more competitive in the future. It seems that free inorganic Fe species are more readily available than those organically complexed (Morel et al., 2008) and it may be that the cellular investment necessary to access organically complexed Fe (see e.g., Maldonado et al., 2006) would become less advantageous as a result of ocean acidification if fFe(II) were to increase due to pH changes. As an illustration, assuming a general half saturation constant (Kµ) for growth as a function of Fe of 0.05 nM and defining waters as nominally "Fe limited" when the bFe concentration <Kµ, we find that the "Fe limited" area of Southern Ocean surface waters (south of 40 • S) either changes insignificantly in response to climate change and ocean acidification (∼0 %) for phytoplankton that access organically complexed Fe or declines greatly (−17 %) if only Fe(II) and Fe(III) are bioavailable (i.e., organically complexed Fe not available). In other words, the "Fe limited area" of the Southern Ocean for phytoplankton that do not access organically complexed Fe (and rely on free inorganic Fe) would be reduced by 17 % in the future due to a shift in the balance between free inorganic and organically complexed Fe. This illustrates the potential advantage that might accrue for phytoplankton species that eschew the cellular investment necessary to access organically complexed Fe in the future "acidified" ocean, thereby making species that rely on free inorganic Fe potentially more competitive. We speculate that this potential lesser reliance on organically complexed Fe by the phytoplankton community (as a result of adaptation or species shifts) and the physiological trade offs that result could have implications for future growth rates and carbon fixation.

Improvements to the speciation model
While our Fe speciation model is complex, relative to contemporary treatments of Fe cycling in global OGCBMs, there are a number of simplifications and processes that could be included/tested in the future. For example, processes such as Fe(III) reduction as mediated by superoxide, the direct photoreduction of Fe(III) and Fe(III) colloids, as well as including the role of superoxide and hydrogen peroxide in the oxidation of Fe(II) could also be important in governing the spatio-temporal variability in Fe(II) concentrations. The Fe speciation model of Ye et al. (2009) is a good candidate model with which to explore the potential importance of such processes. However, as this model is even more complex than our Fe speciation model, it cannot be solved analytically anymore. Nevertheless, an iterative numerical solution is possible and leads to vast savings in computational time compared to solving the full kinetic equations in the one-dimensional setting by Ye et al. (2009) (Völker et al., 2011). A good candidate addition to the current speciation model that would not overcomplicate the analytical solution might be Fe(II) ligands, which could assist in reproducing the relatively high Fe(II) suggested in the low latitudes (e.g., Hansard et al., 2009). The presence of Fe(II) ligands has been noted in rainwater (e.g., Willey et al., 2008) and suggested in the open ocean as well (e.g., Croot et al., 2001) and their presence may assist in stabilising Fe(II) for a number of hours. It is not difficult to include an Fe(II)binding ligand in the analytical solution presented here, but this would require more information on its specific binding strength, as well as its sources and concentration. In terms of dFe, a better understanding of colloidal Fe species (I.e, those in the 0.02-0.2 µm size fraction) and in particular the nature of organic complexation (e.g., Wu et al., 2001;Boye et al., 2010), as well as the degree and timescales of cycling back to soluble Fe species would help us understand how to separate this potentiall important Fe species between "slow" and "fast" timescales. Specific sources of Fe species could also be important, with observational studies also suggesting that continental margins and organic matter remineralisation can supply Fe(II) (Boye et al., 2006;Sarthou et al., 2011), likely stabilised to some degree. However, including and appraising (for example) Fe(II) sources is not possible in our analytical approach and would be better tested by fully prognostic, and thus necessarily regional, Fe speciation models. Nevertheless, our approach of separating the 'fast' and "slow" Fe cycle reactions will permit us to test many question regarding the controls on Fe speciation in the global ocean during future studies such as the presence and cycling of Fe ligands, including perhaps those specific to Fe(II), cycling of colloidal Fe and questions regarding Fe bioavailability (see below).

Modeling Fe binding ligands
We have shown that variability in ocean Fe binding ligands exert a critical control on the speciation of dFe and, as such, on the residence time and bioavailability of dFe.
Here we have used the semi-labile DOC pool as simulated by PISCES, alongside a relationship derived from field observations (Wagener et al., 2008) to allow ligands to vary in our model. However, while this is likely to be a cost effective improvement upon a fixed uniform ligand concentration (as currently employed in other OGCBMs), it will be important to carefully compare the distribution with widespread in situ ligand data to better constrain its viability. For example, while our DOC-linked parameterization will account for ligand production in surface waters, remineralisation (Boyd et al., 2010;Ibisanmi et al., 2011) is not a source of ligands and it absence in our parameterization may hinder the reproduction of subsurface peaks in ligand concentrations (Boye et al., 2010;Ibisanmi et al., 2011;Thuroczy et al., 2011). In the future it may be useful to include a prognostic simulation of the production of weak ligands from organic matter breakdown, as well as the production of strong ligands by the biota, possibly mediated by Fe stress as per Ye et al. (2009), which could be then compared against ligand profiles. A prognostic ligand model would need to be simulated as part of the "slow" Fe cycle reactions and thus require a long model spin up (order of hundreds to thousands of years) in order to correctly simulate deep water concentrations. But this would be feasible in a 3-D global OGCBM, since it would only require the addition of two new tracers (L W and L S ). To that end, more data on the concentrations (especially their profiles, e.g., Boye et al., 2010;Thuroczy et al., 2011;Ibisanmi et al., 2011) and binding strengths of ocean binding ligands will prove invaluable.

Impact of climate and pH on Fe speciation
While our Fe speciation model has difficulties in reproducing the Fe(II) concentrations measured by Hansard et al. (2009, notwithstanding metholodogical issues), our model does a good job in the regions we predict to be impacted by climate change and ocean acidification (the high latitudes). It will be necessary to carefully understand the reasons behind the low modeled dFe concentrations in the tropical Pacific as this certainly restricts the accumulation of Fe(II) therein. Our model suggests that ocean acidification, rather than climate, is likely to exert the strongest control on the evolution of Fe speciation (independent of any effects on absolute dFe concentrations) over the coming century through its mediation of redox kinetics. The greater fraction of Fe(II) we simulate agrees with results from mesocosm experiments using natural seawater with bubbled CO 2 (Breitbarth et al., 2010). Laboratory results using synthetic ligands have shown that the complexation of Fe by ligands could also change with ocean acidification as a function of the degree of protonation of a given ligand, with increases, decreases and no change in complexation possible (Shi et al., 2010). If the in situ ocean ligand pool can be better characterized, and perhaps connected to different production pathways (sensu Hunter and Boyd, 2007), then we could test the combined impact of climate and pH on Fe redox speciation and ligand complexation/cycling in the future. For example, if remineralization is an important source of Fe binding ligands (Boyd et al., 2010;Ibisanmi et al., 2011), then it may be important to understand climate impacts on ligand production rates, as well as their vertical supply from the subsurface to the surface ocean, in concert with effects on Fe redox speciation and potential "direct" effects of pH on organic complexation.. Nevertheless, we note that understanding the ultimate impact on the biota will critically depend on the assumptions regarding the nature of the in situ bFe pool. Reducing the uncertainties associated with the costs and benefits of phytoplankton assimilating different Fe species is crucial.

Modeling Fe bioavailability
Modeled bFe concentrations are highly sensitive to environmental variability and what Fe species are assumed to be available. In the future, it might be worthwhile to parameterize specific accessibilities of different Fe species to phytoplankton. For example, recent kinetic models (e.g., Völker and Wolf-Gladrow, 1999;Shaked et al., 2005;Salmon et al., 2006;Morel et al., 2008) could be included in our model to more mechanistically treat the bioavailability of the different Fe species we simulate. Or we could assume different accessibilities of strong and weakly complexed Fe to different phytoplankton functional types following the ideas of Hutchins et al. (1999). The subsequent impact of changes in Fe speciation on primary productivity could then be assessed. In doing so, it would also be important to add more detail to the formulation of the phytoplankton Fe quota (e.g., Flynn, 2003;Buitenhuis and Geider, 2010) so that the impact of environmental changes on the number of photosynthetic units, Fe concentrations, nitrate reductase, respiration rates as well as different adaptive physiological strategies can also feedback on the phytoplankton demand for Fe (e.g., Raven, 1988;Raven et al., 1999;Strzepek and Harrison, 2004). Assumptions regarding the availability of specific Fe species could then be tested to explore how the impact of environmental variability on Fe speciation might feedback onto viable phytoplankton Fe uptake strategies. In addition, innovative experimental procedures involving characteristic Fe binding ligands and keystone phytoplankton species are critical in understanding the degree to which organically complexed Fe is bioavailable, relative to free inorganic Fe and under what environmental conditions, as well as the cellular trade offs that  result (if any). In doing so, we will be better able to evaluate the impact of climate driven changes in Fe speciation (on multiple timescales) on phytoplankton productivity and growth rates.

Conclusions
Using an analogy with the computation of inorganic carbon speciation in OGCBMs, we outline a means by which Fe speciation can be solved analytically in a cost effective manner in global models. Our approach rests on the division of the Fe cycle into "fast" and "slow" reactions and permits us to simulate 3D Fe speciation using a global OGCBM. We use our model to show that the distribution of different Fe species is tightly controlled by the dFe concentration, the distribution and concentrations of Fe-binding ligands and environmental variables (temperature, light, oxygen and pH). When compared directly to measurements of Fe(II), our model does a reasonable job of reproducing observations in the high latitude oceans (notwithstanding the variability of measured Fe(II) on short time and space scales), but appears to systematically underestimates Fe(II) in the low latitude Pacific Ocean (although these may be overestimated). This could result from errors in the modeled dFe field, the absence of the diurnal cycle and high frequency variability, or missing processes from our Fe cycle model (such as Fe(II) binding ligands or specific Fe(II) sources). Using our model under future climate suggests that climate change and, in particular, ocean acidification will impact Fe cycling and speciation, especially in the Fe limited Southern Ocean. We predict significant increases in Fe(II) due to acidification, which could reduce the "Fe limited area" of the Southern Ocean by ∼20 % for species that rely solely on assimilating on inorganic Fe. We speculate that a dFe pool that has an increased 'free' inorganic component might exert a selective pressure on competitive Fe uptake strategies exhibited by phytoplankton in the future ocean. Finally, our 'analytical solution' approach can be used as a framework within which to test our understanding of Fe speciation at the global scale in future studies if desired. As such, it can provide a means by which new measurements of Fe speciation can be evaluated against hypotheses regarding the quantitative formulation of the processes at play in the ocean's Fe cycle.