Exploring the Iron‐Binding Potential of the Ocean Using a Combined pH and DOC Parameterization

The major part of dissolved iron (DFe) in seawater is bound to organic matter, which prevents iron from adsorptive removal by sinking particles and essentially regulates the residence time of DFe and its availability for marine biota. Characteristics of iron‐binding ligands highly depend on their biological origin and physico‐chemical properties of seawater which may intensely alter under ongoing climate change. To understand environmental controls on the iron binding, we applied a function of pH and dissolved organic carbon (DOC) to describe changes in the binding strength of organic ligands in a global biogeochemical model (REcoM). This function was derived based on calculations using a thermodynamic chemical speciation model NICA. This parameterization considerably improved the modeled DFe distribution, particularly in the surface Pacific and the global mesopelagic and deep waters, compared to our previous model simulations with a constant ligand or one prognostic ligand. This indicates that the organic binding of iron is apparently controlled by seawater pH in addition to its link to organic matter. We tested further the response of this control to environmental changes in a simulation with future pH of a high emission scenario. The response of the binding potential shows a complex pattern in different regions and is regulated by factors that have opposite effects on the binding potential. The relative contributions of these factors can change over time by a continual change of environmental conditions. A dynamic feedback system therefore needs to be considered to predict the marine iron cycle under ongoing climate change.


Introduction
Iron, one of the most biologically important trace metals, controls phytoplankton productivity and thus its CO 2 uptake in a large area of the worlds' ocean, particularly in the high nutrient low chlorophyll regions (Boyd et al., 2007). Furthermore, in tropical oligotrophic waters, iron has a great impact on nutrient consumption and primary productivity via regulation of nitrogen fixation (Mark Moore et al., 2009). Many external sources bring iron into the ocean, for example, dust deposition, rivers, sediments, and hydrothermal vents. Some of these iron inputs vary strongly with climate conditions, which leads to feedbacks between iron supply, marine biogeochemistry, and climate (e.g., Jickells et al., 2005). The iron cycle is therefore highly relevant to climate change.
Iron is the fourth most abundant element in the Earth's crust. The low solubility of iron in the modern oxygenated ocean however limits the oceanic iron inventory and is thus a key factor in determining iron supply to marine phytoplankton. The equilibrium free iron concentration in seawater is around 0.01 μmol m −3 at pH of 8.1 and 25 • C (Liu & Millero, 2002), much lower than the concentration required to support phytoplankton productivity (Sunda & Huntsman, 1995). Appreciably higher in situ dissolved iron (DFe) concentrations are facilitated by complexation with organic matter that acts as a buffer against precipitation and particle adsorption of iron (Gledhill & Buck, 2012;Liu & Millero, 2002). Knowledge of the processes affecting organic complexation is therefore vital for understanding interactions between the marine iron cycle and climate.
Our knowledge of the chemical nature of organic iron-binding ligands, their distribution, and cycling in the ocean has been growing substantially during the last decades. Ligands with high affinity, so-called siderophores, are actively produced by iron stressed microorganisms in order to bind iron for subsequent acquisition by the microbe (Buck et al., 2010;Gledhill et al., 2004;King et al., 2012;Wilhelm & Trick, 1994). Their contribution to the total ligand pool is generally low, but they have been detected in many ocean regions (Boiteau et al., 2016;Bundy et al., 2018;Mawji et al., 2008;Velasquez et al., 2011). Besides siderophores, almost the entire pool of dissolved organic matter (DOM) can bind iron to some extent. A covariance between ligands and dissolved organic carbon (DOC) has been observed during a mesocosms experiment in the Mediterranean Sea (Wagener et al., 2008) as well as an in situ correlation between nitrate/phosphate and Fe solubility (e.g., Schlosser & Croot, 2009). The DOM pool contains various functional groups providing binding sites for iron, for example, amine (as in porphyrin-like ligands), hydroxyl, and carboxyl groups (Zhang et al., 2019). The role of humic substances and exopolysaccharides as organic ligands has been extensively studied (e.g., Hassler et al., 2011;Laglera & van den Berg, 2009). Those compounds could be released by grazing and microbial degradation of organic matter Poorvin et al., 2011;Sato et al., 2007) or have terrestrial origins (Laglera & van den Berg, 2009). Some of the binding sites have higher binding strengths, similar to siderophores, while others are likely much weaker (Gledhill & Buck, 2012).
In addition to binding site composition, environmental conditions also influence the affinity of the organic binding sites for iron. The acid-base chemistry of the binding sites is affected by sea water pH (Gledhill et al., 2015), while their abundance is linked to DOC (Avendaño et al., 2016;Gledhill et al., 2015;Hiemstra & van Riemsdijk, 2006) and DFe concentration (Gledhill & Gerringa, 2017). DFe concentration is also influenced by pH and temperature via controls on solubility (Liu & Millero, 2002), while DFe, pH, and DOC are all influenced by biological productivity. Furthermore, competition between metal ions is also a possible factor affecting binding site affinities (Abualhaija et al., 2015). Therefore, the impact of biological activity on the availability of binding sites and how this changes with environmental conditions are highly complex, and incorporation of these factors into models in order to better predict iron-climate feedbacks is not straightforward.
Thermodynamic chemical speciation models such as the Windermere Humic Acid model (WHAM) (Stockdale et al., 2015(Stockdale et al., , 2016 and the non-ideal competitive adsorption (NICA) model (Kinniburgh et al., 1996;Milne et al., 2003) describe the binding of metals under different environmental conditions. These approaches consider multiple binding sites (WHAM) or competitive interactions between cations and heterogeneous anionic binding sites (NICA) and are thought to better represent the complexity of organic matter binding sites in aqueous environments (Turner et al., 2016). The NICA model is used to describe the binding of metals to a heterogeneous mix of binding sites in natural waters. The model was derived to describe the binding of metals to organic matter in freshwater and soil environments and in these environments electrostatic interactions between charged ions at low ionic strengths are important, so the NICA model is frequently combined with the Donnan model . In seawater, high ionic strengths reduce electrostatic interactions to negligible quantities and the electrostatic component of the model can thus be neglected. The NICA model has been described extensively (e.g., the above refs). Typically the model is based on a bimodial distribution of binding sites, with one mode of sites having acidic properties and described by a NICA affinity constant (NICA logK H + 3.7-4) and the other mode more basic (NICA logK H + 8-10). Since the model is parameterized to account for competitive binding between cations and the hydrogen ion, it can account for the impact of variable pH on metal binding. In freshwater and soil systems, the model has been parameterized for humic and fulvic acids, and a set of generic parameters derived (Milne et al., 2003). Determinations of acid-base properties of dissolved organic matter in marine samples suggest a similar heterogeneous distribution of binding sites (Louis et al., 2009;Muller & Bleie, 2008), and the NICA model has thus been applied to cation binding to marine dissolved organic matter (DOM) (Avendaño et al., 2016;Gledhill & Gerringa, 2017;Hiemstra & van Riemsdijk, 2006;Ulfsbo et al., 2015), although derivation of NICA constants for binding of marine DOM with major cations (H + , Mg 2+ , Ca 2+ ) has yet to be reported.
Global ocean biogeochemical models consider the organic complexation of iron with different degrees of complexities (Tagliabue et al., 2016). While some still assume constant ligand concentrations, others have been trying to capture the spatial variability of organic ligands with different assumptions of their relationship to biological activities. For example, Misumi et al. (2013) linked ligand concentration to oxygen consumption, and Tagliabue and Völker (2011) considered an empirical relationship between total ligand concentration and DOC. In common with other global iron models, we previously described the role of iron-binding ligands in the iron cycle with a constant ligand (Ye & Völker, 2017) and with a prognostic ligand assuming a source of ligands from DOC (Völker & Tagliabue, 2015). Although surface DFe distributions in these models encompassed a similar range to observed DFe concentrations, there were some systematic discrepancies: DFe concentrations are generally too low in the Pacific Ocean, and the models could not reproduce the observed inter-basin pattern in the deep ocean when keeping the range of concentrations close FIGURE 1. Results of fitting iron titration curves using a Langmuir isotherm and assuming two ligands (a); the assumed linear dependence of logk on pH (b) and quadratic dependence on DOC (c). Circles are the iron titration data (a) from the NICA model calculations and derived conditional stability constants for each titration experiment (b and c). In (b) and (c), L 1 is in blue and L 2 in red. Curves of the same color in (b) are simulations at different DOC concentrations and those in (c) are simulations at different pH values.
to measured data (for details, see section 3.1.4). These discrepancies indicate that some processes controlling the spatial variability of ligands and DFe are still likely missing in biogeochemical models of iron.
In this study we make the first steps toward understanding the influence of environmental controls on iron binding in an ocean biogeochemical model. We consider the influence of changes in ocean pH and DOC through simplified equations derived from NICA calculations to parameterize iron binding as a function of pH and DOC. Our study thus takes into account not only a ligand source from DOC but also the environmental control of iron binding by pH. We use the global biogeochemical model in this study to test (1) if the spatial pattern of seawater pH, together with the distribution pattern of DOC, provides more constraints on the spatial variability of organic complexation and improves the model-data agreement of DFe concentration and (2) how the capacity of the ocean to bind iron, DFe concentration, and marine productivity responds to potential climate change driven changes in pH.

Experiments With the NICA Model
We used the NICA model in combination with the speciation program visual Minteq (Visual MINTEQ version 3.0., https://vminteq.lwr.kth.se/) to model iron binding to dissolved organic matter (DOM) over a range of iron concentrations (0.1-2 μmol m −3 ), dissolved organic carbon concentrations (40-100 mmol m −3 ) and pH values (7.2-8.2, free scale). For the NICA model, we applied the general parameters from Milne et al. (2003) for DOM binding with H + , Ca 2+ , and Mg 2+ and used parameters empirically derived in Avendaño et al. (2016) for iron binding. The generic acid-base parameters and binding site densities are thus derived from terrestrial humic substances, since there were no equivalent parameters available for marine DOM at the time of this study. Iron binding parameters were derived from iron titrations of North West European shelf waters undertaken over a range of pH values (Avendaño et al., 2016). It should be noted that the high proton affinity phenolic functional groups dominate iron binding over the pH range explored in this study. The sensitivity of our results to changes in temperature (5-25 • C) has been tested, but since there is very little parameterization for reaction enthalpies in the NICA model as applied in visual MINTEQ, we observed a minimal effect of temperature on iron binding and thus neglected temperature in further calculations. Data calculated with the NICA model are available in the supporting information Table S1.

Application of a Langmuir Isotherm to Titration Data
A Langmuir isotherm with an assumption of two ligands was fitted to the iron titration curves simulated by the NICA model ( Figure 1a). Conditional stability constant (K 1 and K 2 ) and ligand concentration (L tot1 and L tot2 ) were derived for each NICA simulation with fixed values of DOC and pH.
NICA simulations at different pH and DOC concentrations allowed us to derive a dependence of the conditional stability constants and ligand concentrations on pH and DOC. Plotting logK against pH in all experiments revealed a linear relationship between these two parameters ( Figure 1b). In contrast, a polynomial fit was used for the dependence of logK on DOC. A quadratic function of DOC ( Figure 1c) was found to best describe the variations of the simulated titration curves (Figure 1a). The ligand concentrations did not vary greatly with pH or DOC and were thus assumed to remain constant at 0.6 μmol m −3 for L tot1 and 1.7 μmol m −3 for L tot2 . Thus, the functions (Equations 2 and 3) of pH and DOC were derived for the hypothetical conditional stability constants k2 and k1 (or logk 2 and logk 1 , respectively). This fitting result is pure operational for the best reproduction of data from the NICA simulations. It should be emphasized that these constants are not comparable to the thermodynamic equilibrium constants determined from experimental iron-ligand titration curves (e.g., K cond FeL( Fe ′ ) ) and the fixed values of L tot1 and L tot2 were not chosen based on observed ligand concentrations in the ocean.

Parameterization of the pH and DOC Dependency of Fe Complexation in a Global Model
The two equations (Equations 2 and 3) were implemented into a 3D global biogeochemical model REcoM2 (Regulated Ecosystem Model 2) (Hauck et al., 2013) with a complex description of the iron cycle (Ye & Völker, 2017). REcoM2 describes two phytoplankton classes, diatoms and non-diatoms; a generic zooplankton; and one class of organic sinking particles whose sinking speed increases with depth. The iron cycle in the model is driven by two external iron sources: dust deposition and sedimentary input, and by biological uptake and remineralization, and scavenging onto biogenic and lithogenic particles. DFe input flux from the atmosphere is calculated from the daily dust deposition field from Mahowald et al. (2003) assuming 3.5% iron content in dust particles and 2% solubility. Only unbound iron undergoes scavenging and scavenging removal of DFe is proportional to DFe and particle concentrations.
We considered two ligands in the model and used GLODAP pH data Olsen et al., 2016) and Equations 2 and 3 in combination with modeled DOC to calculate logk 1 and logk 2 with constant ligand concentrations. Since our model derives labile DOC, we added a refractory fraction of 40 mmol C m −3 in order to take into account the entire DOC pool. The constants from a to e were determined to be −2·10 −4 , 0.034, −1.67, 24.36, and 2.67, respectively. This parameterization of organic complexation is expected to produce a higher variability of DFe distribution but also much higher DFe concentrations than assuming a constant ligand concentration of 1 μmol m −3 (the total ligand concentration is 0.6 + 1.7 = 2.3 μmol m −3 ), if using the same scavenging rates as in Ye and Völker (2017). To keep modeled DFe close to global observations, the scavenging rate for organic particles was increased to 0.752 (mmol C m −3 ) −1 day −1 which is about 50-fold the rate used in Ye and Völker (2017). More details of the model description of particle dynamics as well as the model setup for the iron cycle can be found in Ye and Völker (2017).
REcoM2 is coupled with MITgcm (Marshall et al., 1997), spanning the latitude range from 80 • N to 80 • S at a zonal resolution of 2 • and a meridional resolution between 0.39 • and 2 • . In the vertical it has 30 layers, increasing in thickness from 10 mat the surface to 500 m below a depth of 3,700 m. All simulations in this study were run for 1,000 years, and the output of the last 10 years was used for analysis.

Experiments With REcoM
A standard run R stand was conducted with the parameters given above and GLODAP present-day pH values. Model sensitivity of the binding potential of organic ligands, DFe concentration, and thus biological production to the pH change was studied in R ph , with the same parameters but future pH values predicted in CMIP5 (Coupled Model Intercomparison Project, https://cmip.llnl.gov/cmip5/index.html) using NorESM1-ME for a high emissions scenario (RCP8.5) (Taylor et al., 2012). The simulation with NorESM1-ME was carried out for the period from 2006 to 2100 with a 900-year spin-up, and the output of the year 2100 was used to run R ph . Details of this pH field and its difference to GLODAP data are described in section 3.3.1.
R stand of this study was compared with a model run with a constant ligand concentration of 1 μmol m −3 (Ye & Völker, 2017) (R constL ) and a run with one prognostic ligand (Völker & Tagliabue, 2015) (R progL ) (section 3.1.4). R constL assumed a conditional stability constant of 10 11.3 μmol −1 m 3 and took additionally into account a scavenging loss of iron by absorption onto lithogenic particles, besides that onto biogenic particles. R progL used the same conditional stability constant and described the cycling of one ligand class with the following sources and sinks: The ligand is released by degradation of particulate organic matter and produced by living organisms (related to DOC production) and is destroyed by microbial and photochemical degradation and removed from water column by biological uptake and formation of aggregates. R progL and R stand have some common features in the biological control of iron binding: Organic complexation of iron is enhanced in regions with high biological productivity, since degradation of a large amount of organic matter releases ligand in R progL and DOC which leads to high binding potential of ligands in R stand . It is however unique in R stand that degradation of particulate organic carbon releases DIC and protons which lower pH and can affect iron binding in more remote regions via transport of low pH water masses. Both R constL and R progL were run for 1,000 years, and the output of the last 10 years was used for comparison. Model-data comparison was made using the collected observations in GEOTRACES IDP2017 (Schlitzer et al., 2018).

Modeled DFe 3.1.1. DFe in Surface Waters
The modeled surface DFe (<100 m) is high in regions receiving most of the dust and in coastal regions and shows an inter-basin gradient from the Atlantic and Indian Ocean, over the Southern Ocean and to the Pacific Ocean ( Figure 2a). In most parts of the ocean, surface DFe is close to observations, except in the subtropical North Atlantic and northern part of Arabian Sea where an overestimation of aeolian iron input in the model seems likely.

DFe in Intermediate Waters
In the intermediate waters (500-1,000 m), the modeled DFe concentrations are clearly higher than those at the surface and particularly high in the North and eastern equatorial Pacific, northern part of Indian Ocean and slightly enhanced in the eastern Atlantic close to the African continent and in some parts of the Southern Ocean ( Figure 2b). Lower concentrations are found in the western Atlantic Ocean and South subtropical gyres in the Pacific and Indian Ocean.
This pattern is basically confirmed by the observations although DFe is apparently too low in the tropical and subtropical North Atlantic Ocean and too high in some regions in the Southern Ocean. These are primarily consequences of the modeled pattern of DOC concentrations. Only one labile DOC fraction is considered in our model with a life time of 10 days and a temperature-dependent degradation. The refractory DOC fraction is assumed to be constant at 40 mmol C m −3 . This results in too high DOC concentrations in the Southern Ocean and too low values in the subtropics, compared with the measured and modeled global DOC distribution (Hansell et al., 2009).
The enhanced DFe in the upwelling regions (e.g., the Mauritanian and Peru upwelling region) in the model is consistent with observations of enhanced DFe in these regions (e.g., Buck et al., 2018;Hatta et al., 2015). These are typically attributed to remineralization processes . Misumi et al. (2013) used a parameterization of ligand concentration based on Apparent Oxygen Utilization (AOU) to capture this feature in their model. Here, we achieved a similar result in R stand . However, in our case the increased DFe is not explained by the release of ligands during remineralization of organic matter (as in Misumi et al., 2013), but by the lower pH in upwelled waters.

DFe in Deep Waters
In the deep ocean (2,000-3,000 m), the model shows generally lower concentrations in the Atlantic Ocean compared to other ocean basins ( Figure 2c) and a weak inter-basin gradient of ∼0.1 μmol m −3 is found between the Atlantic and Pacific Ocean. To compare with the measured DFe concentrations in the deep Atlantic and Pacific Ocean, we calculated the basin mean for all data measured between 2,000 and 3,000 m. Since our model does not take into account hydrothermal sources, some high values (>2 μmol m −3 ) measured along the GP16 cruise in the Pacific (Fitzsimmons et al.,201), GA02 (Rijkenberg et al., 2014), GA03 , and GAc01 (Saito et al., 2013) in the Atlantic Ocean were omitted for the statistical analysis (i.e., basin mean, bias, and root mean square). An inter-basin gradient is also found in the measurements with averages of 0.72 ± 0.20 μmol m −3 in the Atlantic and 0.93 ± 0.33 μmol m −3 in the Pacific, in agreement with the modeled deep DFe pattern.

Model Runs With Different Parameterization of Fe Organic Complexation
We compared three model runs with different parameterization of organic complexation: R stand , R constL , and R progL (for details, see section 2.4). While the three runs were tuned to match the range of observed surface DFe concentrations, there are some significant spacial discrepancies in DFe distributions. In R constL , DFe concentrations in part of the surface ocean were too low and resulted in an overestimation of the extent of iron limitation in the Pacific and Southern Ocean ( Figure S1A). In R progL , DFe was far too high in regions receiving atmospheric dust ( Figure S2A). Furthermore in the deep ocean (Figures S1B, S1C, S2B, and S2C), both R constL and R progL failed to reproduce the observed basin variability in regions which are not affected by hydrothermal vents (Ellwood et al., 2018;Nishioka et al., 2013;Resing et al., 2015). In the deep Atlantic Ocean, R constL produced too low DFe due to strong scavenging removal by lithogenic particles while a clear positive Pacific-Atlantic gradient with higher concentrations in the deep Atlantic Ocean was obtained in R progL . The latter pattern can not be confirmed by observations Rijkenberg et al., 2014).
DFe distribution changed considerably in R stand when compared to R constL and R progL . In surface waters, modeled concentrations of DFe are generally elevated in part of the Pacific and Southern Ocean with higher biological productivity (Figures 3a and 3b), which increases the spatial variability of DFe in these regions. The basin mean in the two ocean basins is 0.15 and 0.27 μmol m −3 , respectively. In contrast, DFe in R constL and R progL is nearly uniform in the Pacific Ocean and in both ocean basins clearly lower than the observations. The bias is then improved from −0.1 in both R constL and R progL to −0.05, μmol m −3 in R stand in the Pacific and from −0.06 in both R constL and R progL to 0.01 μmol m −3 in R stand in the Southern Ocean. This is likely a result of a greater overall potential for iron binding linked to increased DOC.
The too high concentrations of DFe produced in R progL in part of the Atlantic and Indian Ocean under dust plumes are significantly alleviated in the current model ( Figure 3b). Aerosol-derived inputs are more spatially constrained in the current model when compared with R constL , which also considered lithogenic scavenging (Figure 3a). Although the localized impact of dust derived iron is high (even higher than in R constL in the subtropical North Atlantic due to the higher binding potential in the current model), the added iron is not propagated so that overall the concentrations of DFe in dust influenced areas of the surface ocean are lower in R stand . This is explained by the link between ligand binding potential and DOC: since the productivity decreases dramatically when moving from the equator to subtropical gyres, DOC concentrations, and thus ligand binding potential also decrease. The dependency of ligand binding potential on DOC concentration creates a link between growth limitation by macronutrients and the DFe concentration in the gyres.
The differences from R stand to R constL and R progL show similar patterns in the intermediate and deep waters (Figures 3c to 3f). R stand has a clear increase of DFe in both intermediate and deep waters compared to R constL and a reversed inter-basin trend between the Atlantic and Pacific Ocean compared to R progL . Vertically, the basin-wide averaged DFe (not shown) increases with depth in R stand from the surface down to its maximum between 800-1,000 m, as typically observed in field studies (e.g., Boyd & Ellwood, 2010). Below that, DFe remains nearly constant with depth. In both R constL and R progL , the deep DFe maximum is about 200 m shallower. Overall, an improvement of bias has been achieved from 0.1 to −0.05 μmol m −3 in Deep Atlantic and from −0.51 to −0.18 μmol m −3 in Deep Pacific. Furthermore, RMS (root mean square) has been reduced from 0.25 to 0.21 in Deep Atlantic and from 0.62 to 0.39 in Deep Pacific. However, since the measurements of deep DFe concentrations are still scarce and the DFe ranges in the two deep ocean basins are relatively close to each other (see section 3.1.3), this pattern remains to be verified when more data become available.

Interactions With the Binding Potential of Ligands
We have shown that the combination of lithogenic scavenging and two prognostic ligands parameterized by pH and DOC concentration result in improvements to the modeled distribution of DFe relative to our other models. Here we further examine why these changes occur using the concept of the side reaction coefficient FeL , which describes the binding potential of ligands for iron in natural organic matter.

Side Reaction Coefficient FeL
Generally, the side reaction coefficient describes the likelihood that a particular complex will form given the presence of competing ions within a complex mixture (Ringbom, 1963). For the complex FeL, it can be derived from concentration of free ligand L ′ and conditional stability constant K cond FeL(Fe ′ ) or from the formed complexes FeL and free iron concentration Fe ′ (Equation 4). Using FeL , we describe not only the current relationship between ligand and DFe but also the potential of the ocean to keep iron in solution if external iron supply changes.
We use the side reaction coefficient both in order to compare modeled data with observations of iron speciation and in order to compare the potential for different regions of the ocean to bind more iron (i.e., where FeL is higher). We use FeL to compare modeled to observed iron speciation because differences in applied detection windows during analysis can lead to systematic bias in ligand concentrations and logK (Gledhill & Gerringa, 2017).

Modeled Side Reaction Coefficient
We first calculated free iron concentration Fe ′ from modeled DFe, total ligand concentration, k 1 and k 2 , and then FeL using a variant of Equation 4: The side reaction coefficient tends to be high either where organic matter concentrations are high or where DFe inputs are low so that more free binding sites exist in the water. However, FeL is also affected by pH since it depends on changes in the relative binding strengths of organic ligands and hydroxide ions. These features are reflected in the calculation of FeL (Equation 4) via the parameterization of the hypothetical conditional stability constants (k) (Equations 2 and 3) in the model. In surface waters, k is regulated by both DOC and pH; while in deep waters, pH control dominates the spatial distribution of k because DOC concentrations are nearly homogeneous.
At the surface, modeled FeL is low in Fe-replete regions, for example, the tropical and subtropical Atlantic Ocean, Arabian Sea and most coastal regions, and high in Fe-limiting regions, for example, the North and equatorial Pacific and Southern Ocean (Figure 4a). FeL is moderate in the subtropical Pacific gyres and coastal regions. Biological productivity and thus DOC production is extremely low in the Pacific Gyres ( Figure  S3B), explaining the relatively low logk there in the model. Lacking a considerable iron source, the binding potential in these regions remains moderate. In coastal regions, iron supply is high and thus phytoplankton growth is predominately limited by macronutrients. With the flexible stoichiometry, the model can describe the overconsumption of carbon by phytoplankton when growing under strong nitrogen limitation (Schartau et al., 2007). This leads to a higher DOC exudation in regions strongly limited by nitrogen and thus intensifies the DOC production in the coastal regions. High DOC and high DFe work against each other in affecting the binding potential of sea water and result in moderate values of FeL as well.
Below the surface, DFe concentrations increase as a result of remineralization (Figure 2b). DOC varies less in intermediate waters, so FeL is nearly homogeneous around 10 3.5 (Figure 4b). Slightly enhanced values are obtained in the gyres. Due to the low productivity in surface waters, less iron is remineralized by microbial degradation of organic matter. This is reflected by low DFe in intermediate waters (Figure 2b) in the subtropical Atlantic and Indian Gyres, and the South Pacific Gyre. Thus, many binding sites of ligands there remain unoccupied by iron, resulting in a higher binding potential. Furthermore, the constant refractory DOC concentration results in a greater sensitivity of FeL and thus DFe to pH. The significantly lower pH in the North Pacific (7.5-7.9) compared to the other gyres (pH> 8.0, see Figure S3C), results in higher FeL in the North Pacific Gyre. Here we see that reduced pH leads to both higher binding potential and higher DFe concentrations (Figure 2b).
In the deep ocean, the spatial variability of FeL becomes weaker (Figure 4c). Higher FeL values are found in the North Pacific due to lower pH and in some dust regions, for example, the subtropical North Atlantic and Arabian Sea, caused by strong scavenging of iron by small lithogenic particles in the deep ocean.
We compared FeL calculated from model data with that calculated from observed data using Equation 5. The comparison was conducted for all depth layers (Figure 4), but we focus here on the surface layer ( Figure  4a), because conditional stability constants of organic ligands are usually measured after buffering pH to a range between 8.0-8.2 and observed FeL will thus only be valid for this pH range. The calculated values for observations at depth, where pH in fact is generally lower than this range, can thus be biased.
Both the modeled and observed FeL values vary in the range between 10 and 10 4 . In the subtropical North Atlantic Ocean, the observed FeL is clearly higher than the modeled FeL (up to two orders of magnitudes), whereas in the Southern Ocean the majority of the observed FeL is lower than the modeled one. Several reasons related to both the analytical determination of iron binding and to the model parameterization should be considered to explain these discrepancies. The analytical determination of side reaction coefficients can be subject to bias resulting from the conditions and assumptions made during the analysis and interpretation of titration data (Gledhill & Gerringa, 2017;Pižeta et al., 2015). Thus, in the Southern Ocean, conditions applied during the analytical determination of iron complexation may result in lower FeL because the applied detection window is too low to detect very strong iron-binding sites, which will play a more important role when iron concentrations are low. In contrast in the tropical North Atlantic, where DOC concentrations are low, but iron concentrations are relatively high, the applied detection window may be too strong, and thus weaker binding sites would not be detected.
With respect to the model parameterization, we assume a constant solubility of iron from dust particles (2%) which could result in a too high iron supply by dust deposition close to the source regions (Z. , resulting in a too high Fe ′ concentration and thus a too low FeL : Fe ′ ratio ( FeL ) under the Saharan dust plume. Furthermore, the model produces too high DOC concentrations in the Southern Ocean and too low values in the subtropics (section 3.1.2), which will have a concomitant impact on FeL . In addition, the functions implemented in our model to describe the binding strength of organic ligands have been derived by fitting the data from the NICA experiments which assume that all ligand classes vary in proportion to DOC concentrations. This could limit the power of using this parameterization to explain the overall binding potential in surface waters where other types of ligands also play a role. For example siderophores are more influenced by iron availability and microbial community composition than DOC concentrations. Finally, we compare the annual mean of the model with observational data that may contain a seasonal signal.

Binding Potential of the World's Ocean
Using the side reaction coefficient of organic iron-binding ligands ( FeL ), one can explore the relationship between binding potential and DFe and provide a further constraint on iron models in analogy to the application of L* (excess ligand over DFe) (e.g., Boyd & Tagliabue, 2015). Furthermore, we suggest the binding potential of the world's ocean could be a helpful parameter for predicting the response of the marine iron cycle to environmental changes.
The surface ocean can be divided into three domains based on the binding potential FeL in the upper 100 m: a high, medium, and low potential region. We defined a high-potential region by FeL ≥ 1, 000 which covers the largest area of the world ocean. The region with a medium binding potential (defined by 100 ≤ FeL < 1,000) includes the Pacific Gyres, some regions strongly influenced by dust input of iron, and coastal regions. In regions with extremely high iron supply, for example, the subtropical North Atlantic and Arabian Sea, almost all the binding sites of ligands become occupied, leading to a very low potential of ligands (defined by FeL < 100).
Future environmental changes such as rising temperatures, ocean acidification, sea-level rise, and deoxygenation will affect the marine biogeochemical cycles in the coming decades in many ways (Gruber, 2011). . Difference between pH from CMIP5 NorESM1-ME used in R ph and that from GLODAP used in R stand in the upper 100 m (a), between 500-1,000 m (b) and between 2,000-3,000 m (c).
In the following section, we want to use the concept of iron-binding potential to investigate a hitherto unconsidered pathway how ocean acidification affects the marine biogeochemistry, namely, via the response of iron-binding strength of ligands to pH change, and its effect on the marine iron cycle and the CO 2 uptake by marine phytoplankton.

Response of the Iron Cycle to Future pH Change
In a new model simulation R ph , we replaced the GLODAP pH field with the one simulated for CMIP5 using NorESM1-ME for a high emissions scenario (RCP8.5). Other parameters remain the same as in the standard simulation (R stand ). With this setup, we aimed to improve our understanding and prediction of the future marine iron cycle, being aware of some limitations of this approach. pH does not only alter organic complexation of iron but also the inorganic iron solubility and phytoplankton growth, for example, the shell formation of calcifying organisms (e.g., Orr et al., 2005). Since the pH field is prescribed in R ph and only considered to calculate the binding strength of organic ligands, the inorganic solubility of iron and calcification are not affected. Furthermore, there is no feedback between biological DIC consumption and pH.

Response of DFe and NPP
pH predicted in CMIP5 NorESM1-ME decreases strongly at the surface compared to GLODAP data (Figure 5a), whereas deeper in the ocean, the North Pacific and Southern Ocean show elevated values to different extents (Figure 5b and 5c). The increase of pH in the intermediate and deep water is mainly caused by the model-data discrepancy between the simulation for present-day (Year 2005) and GLODAP data (see Text S1 and Figure S4). This is tempered a little by the pH decrease from the present-day to the high emission scenario in NorESM1-ME.
Reduced pH in the surface ocean results in an increase of logk ( Figure S5) and in consequence an increase in DFe everywhere in the surface ocean except for some limited regions in the Pacific Gyres (Figure 6a). Furthermore, higher concentrations of DFe are found in the deep ocean except for the high latitudes in the Pacific Ocean ( Figure S6). Overall the entire DFe inventory is enhanced by approximately 10%.
As a consequence, global net primary production increases from 42.6 to 44.8 Pg C year −1 (5%), with a general increase in the HNLC regions and decrease in the Indian and Atlantic Ocean, and some parts in the Southern Ocean (Figure 6b). The export production increases from 10.8 to 11.2 Pg C year −1 . The largest change is found in the South Pacific Ocean between 30 • S and 40 • S where the production was severely limited by iron in R stand . Small phytoplankton benefit most from this release of iron limitation in R ph and NPP for this population increases by up to fivefold. The relative increase in NPP is less than that of DFe because besides iron, phytoplankton growth is also limited by macronutrients and light. The increase of DFe in Fe-limited regions leads to higher NPP and higher consumption of macronutrients. Overall, this is offset by intensification of growth limitation by macronutrients in some presently Fe-replete regions (e.g., the Atlantic and Indian Ocean) where a decrease of NPP is observed. A similar effect has also been found in oligotrophic regions in the global simulation of in situ Fe fertilization studies (Aumont & Bopp, 2006). Surprisingly, NPP is unchanged or decreased in the Southern Ocean (> 60 • S) despite a marginal increase in DFe. Here the increased response of small phytoplankton and their resultant higher density causes stronger attenuation of light so that diatoms become more light-limited and total NPP in that region thus decreases.

Factors Driving Changes in Ligand Binding Potential FeL
The parameters used in the NICA model for this study result in an increase in organic binding site strength relative to hydroxide complexation as pH decreases (Avendaño et al., 2016). This manifests itself in an increased logk throughout the surface ocean in R ph ( Figure S5) and has an impact on FeL , since FeL is the product of logk and free ligand (Equation 4). Nevertheless, FeL decreases in most of the surface ocean in R ph (Figure 6c). This pattern is determined by the complex interactions between DFe, pH and DOC and NPP. We explain these interactions in this section by illustrating various pathways regulating the different responses of FeL to the pH change (Figure 7). First of all, we should keep in mind that an increase in DFe (and thus Fe ′ ) and an increase in logk (due to decreased pH, possibly further intensified by increased DOC) work in the opposite way to affect FeL . Therefore, the balance between the negative effect by DFe and the positive effect by logk can help us understand the response of FeL .
The green loop in Figure 7 demonstrates a pathway in which the positive effect by logk dominates. Here, significant change in both pH and DOC result in a large increase of FeL . This pathway explains the change in the South Pacific between 30 and 40 • S in the model (Figure 6c). Being far from any significant iron sources, newly bound iron is to a large part consumed by the enhanced productivity (Figure 6b) which results in much higher DOC and thus an increase of logk. Therefore, the binding potential increases.
The red loop shows a different trend: DFe increases as expected from the increase of logk and the negative effect of DFe on FeL dominates. This pathway explains the change in most of the HNLC regions including the Southern Ocean north of 60 • , as well as the large area of the Pacific Gyres where FeL declines by up to one order of magnitude ( Figure 6c). With more DFe becoming available, phytoplankton growth and thus DOC production increased until they are switched from Fe-to N-or Si-limited, preventing further increase of logk. Hence, more binding sites are occupied by iron, leading to a lower free ligand concentration and reduced binding potential.
Further types of pathways are found in regions where phytoplankton growth does not respond positively to increased DFe and thus the effect of DOC on logk can be neglected. The higher logk, which results from a lower pH, enhances DFe but not biological production because of the lower availability of macronutrients or stronger light limitation addressed above. Whether FeL responds positively or negatively in these pathways (Figure 7, blue loops), depends on the strength of the iron source. In the North Atlantic where the iron supply is very high (Figure 7, blue solid loop), much more iron is prevented from scavenging removal by binding with ligands. Thus many binding sites of ligands are occupied and the binding potential FeL decreases. It works in a similar way to most of the HNLC regions (red loop). These two types of regions do not have a similar magnitude of iron supply, but the relative importance of the iron supply to the change in DOC FIGURE 8. Conceptual presentation of the relationship between DFe, DOC and the binding potential of ligands in the model. The shape of the isolines has been derived based on the surface distribution of log along the axes of DOC and DFe concentration in R stand ( Figure S7). The slope of the isolines should qualitatively illustrate the trend and may be larger than in Figure  S7. Red, green, and blue arrows show the response at different stages within a temporal evolution, and grey arrows show the direction of a temporal evolution.
production is similar so that the negative effect of DFe dominates in both pathways.
Another pathway without significant biological response is found if the positive effect by logk (due to pH drop) overpowers the negative one by DFe (Figure 7, blue dashed loop). For example, some regions in the North Atlantic the South Atlantic, parts of the Indian Ocean and the Southern Ocean south of 60 • S have a lower iron supply and similar logk values to the North Atlantic, thus more binding sites remain unoccupied, resulting in an increase of FeL .

Evolution of the Response Type
While these are response types for a specific climate condition, a continual change will lead to changes in which type of response dominates in a specific region over time, particularly when besides pH, the external sources of iron also vary in a future scenario. To illustrate the change of response type over time, we plotted the isolines of log( ) against DFe and DOC concentration in surface waters in Figure 8 and discuss the temporal evolution by moving a focused region (the black dashed cross) in directions with different conditions. The starting point of the black dashed cross could be any position depending on the initial FeL value.
Six stages can be identified during the evolution. We start from Stage I (Figure 8, along the green solid arrow) which presents a positive knock on effect in the strongly Fe-limited regions. Here an increase in DFe due to higher logk (given by lower pH) leads to a higher production and DOC concentration, and further higher logk and FeL (described as the green loop in Figure 7). This knock on effect will continue until phytoplankton growth is shifted from predominantly Fe-limitation to N/Si-limitation or light-limitation. Then the biological production increases more slowly and with it the DOC concentration. If the iron supply remains high or an external source is enhanced (e.g., caused by climate change), we will move on to Stage II which is characterized by the larger effect of DFe relative to that of DOC (described as the red loop in Figure 7). This results in a decrease of FeL ( Figure  8, along the red solid arrow). With a further rising DFe concentration, Stage III of Fe-replete regions will be generated (described as the solid blue loop in Figure 7). Here biological production is not affected any more by the change in DFe and free binding sites of organic ligands are reduced by a further increase of DFe, leading to a decrease of FeL (Figure 8, moving along the blue solid arrow).
The evolution in the other direction is demonstrated with the dotted arrows in Figure 8. This scenario requires a decline of DFe which is not presented in R ph but has been posited under certain future climate scenarios, for example, dust deposition could decrease in some regions in the future (Kok et al., 2018). Starting in a Fe-replete region in Stage III, we now move with the decline of iron supply toward a system (Stage IV) where DFe concentration is still high enough so that phytoplankton growth is hardly affected by the DFe decline, but more binding sites of organic ligands become free, resulting in an increase of FeL (Figure 8, moving along the blue dotted arrow). If the reduction of iron supply further progresses, phytoplankton becomes Fe-limited again and DOC decreases slowly. Then we arrive in Stage V (along the red dotted arrow) where the effect of DFe is stronger than that of DOC, similar to Stage II, and thus FeL increases slowly. With a continuous decrease of DFe, we will end up with Stage VI with a negative knock on effect, since the Fe-limitation gets intensified and DOC drops rapidly, resulting in a decrease of FeL .

Conclusions
In this study a function of pH and DOC has been derived to describe changes in the binding strength of organic ligands based on a heterogeneous distribution of many binding sites of different strengths as described in the NICA model. The function has been applied in the global biogeochemical model REcoM to parameterize the dependence of iron binding on pH and DOC. The model also incorporated scavenging onto both biogenic and lithogenic particles.
In the surface ocean, the modeled binding by organic ligands is determined by both DOC and pH, whereas deeper in the water column, pH alone controls binding since DOC concentration is low and distributed almost homogeneously. DFe in the surface Pacific and Southern Ocean shows a significant increase and higher spatial variability, compared to previous model simulations with a constant ligand or one prognostic ligand released from degradation of organic matter. In the mesopelagic and deep waters, DFe is considerably elevated, particularly in the Pacific and Indian Ocean, and a Pacific-Atlantic gradient has been generated in the deep ocean with higher concentrations in the deep Pacific. We compared model results with GEOTRACES measurements and the model-data agreement of DFe concentration was improved in both the surface and deep ocean, although the deep Pacific-Atlantic gradient remains to be proven as more observational data become available. The improvements in model-data agreement indicates that the organic binding of iron is apparently controlled by sea water pH in addition to its link to organic matter. We tested the response of this control to future environmental changes in a sensitivity study by considering the future pH of a high emission scenario. We examined changes in biological productivity, DFe and the binding potential of organic ligands to future changes in pH. We found a 5% increase in NPP, which resulted from a 10% increase in the DFe inventory. Two factors, the change in binding strength itself (resulted from the parameterization with pH and DOC) and that in DFe (caused by the change in binding strength), act in an opposite way in regulating the ligand binding potential. Hence, how the iron cycle in different ocean regions responds to the future pH mainly depends on the relative contribution of these two factors, and furthermore on Fe-limitation status and the ability of Fe acquisition of different phytoplankton groups. Furthermore, the pathway of response of a specific ocean region to further changes can evolve over time resulting in a dynamic feedback system. However, one needs to be careful with the interpretation of the results of the simulation with future pH. In this study we only examined the impact of pH on primary productivity with respect to organic complexation of iron. Changes in other processes such as phytoplankton uptake of iron, nitrogen fixation and phytoplankton iron demand (Hutchins & Boyd, 2016) were not taken into account. Studies on those processes reported different responses to pH change (e.g., Breitbarth  This study is the first to investigate the effect of pH on organic complexation of iron in a global iron model. Despite some limitations, for example, the pH field is prescribed and only the labile DOC pool is calculated in the model, our results strongly support further biogeochemical modeling studies that include the relationship between organic complexation of iron and sea water pH, particularly for predicting the future iron cycle under ocean acidification. Inclusion of an interactive description of pH and a good representation of present-day DOC concentrations together with more comprehensive understanding of mechanisms controlling phytoplankton growth under Fe-limitation would further improve the predictive power of the models. parameterization of organic complexation. PANGAEA (https://doi. org/10.1594/PANGAEA.915466). This work is supported by the project PalMod (PaleoModeling, BMBF 01LP1505C) and DFG project (Development of a consistent thermodynamic model of trace element-organic matter interactions in the Ocean, GL 807/2-1).