The iron budget in ocean surface waters in the 20th and 21st centuries: projections by the Community Earth System Model version 1

We investigated the simulated iron budget in ocean surface waters in the 1990s and 2090s using the Community Earth System Model version 1 and the Representative Concentration Pathway 8.5 future CO 2 emission scenario. We assumed that exogenous iron inputs did not change during the whole simulation period; thus, iron budget changes were attributed solely to changes in ocean circulation and mixing in response to projected global warming, and the resulting impacts on marine biogeochemistry. The model simulated the major features of ocean circulation and dissolved iron distribution for the present climate. Detailed iron budget analysis revealed that roughly 70% of the iron supplied to surface waters in high-nutrient, low-chlorophyll (HNLC) regions is contributed by ocean circulation and mixing processes, but the dominant supply mechanism differed by region: upwelling in the eastern equatorial Pacific and vertical mixing in the Southern Ocean. For the 2090s, our model projected an increased iron supply to HNLC waters, even though enhanced stratification was predicted to reduce iron entrainment from deeper waters. This unexpected result is attributed largely to changes in gyre-scale circulations that intensified the advective supply of iron to HNLC waters. The simulated primary and export production in the 2090s decreased globally by 6 and 13%, respectively, whereas in the HNLC regions, they increased by 11 and 6%, respectively. Roughly half of the elevated production could be attributed to the intensified iron supply. The projected ocean circulation and mixing changes are consistent with recent observations of responses to the warming climate and with other Coupled Model Intercomparison Project model projections. We conclude that future ocean circulation has the potential to increase iron supply to HNLC waters and will potentially buffer future reductions in ocean productivity.


Introduction
Iron is an essential nutrient for marine phytoplankton growth and limits primary production in the three major highnutrient low-chlorophyll (HNLC) regions of the world's oceans, namely, the eastern equatorial Pacific, the Southern Ocean, and the subarctic North Pacific (Martin and Fitzwater, 1988;Martin et al., 1990Martin et al., , 1994. To reduce the uncertainties of climate projections, it is essential to understand iron cycling in the present ocean and to evaluate changes in response to global warming. The iron budget in ocean surface waters is determined by a combination of exogenous inputs and physical and biogeochemical ocean processes (e.g., Fung et al., 2000). Various processes are responsible for the exogenous inputs: deposition of iron-bearing dust (Duce and Tindale, 1991;Tegen and Fung, 1994;Mahowald et al., 1999;Jickells et al., 2005); reduction and resuspension of sedimentary iron (Johnson et al., 1999;Elrod et al., 2004;Moore et al., 2004;Lam et al., 2006;Nishioka et al., 2007;Lam and Bishop, 2008;Moore and Braucher, 2008;Tagliabue et al., 2009;Misumi et al., 2011); and fluvial and hydrothermal inputs (Wetz et al., 2006;Tagliabue et al., 2010;Klunder et al., 2011;Nishioka et al., 2013;Saito et al., 2013). Most sedimentary and hydrothermal iron is supplied from continental margins and the ocean abyss to open ocean surface waters by physical ocean processes. Because of the relatively short residence time of iron in seawater, favorable ocean bathymetry, circulation, or mixing conditions are required for iron to be transported from sedimentary or hydrothermal sources to surface waters of the open ocean; such conditions may occur, for example, downstream of islands (Blain et al., 2001) or where ocean currents flow swiftly (de Baar et al., 1995;Mackey et al., 2002;Slemons et al., 2009;Misumi et al., 2011). Complexation with organic ligands greatly lengthens the residence time of iron in seawater (e.g., Misumi et al., 2011), so the concurrent release of iron and organic ligands from hydrothermal vents (Bennett et al., 2008;Toner et al., 2009) may facilitate long-distance iron transport from a hydrothermal source. Most iron in rivers is thought to flocculate and be precipitated in estuaries, because seawater cations neutralize the negatively charged iron-bearing colloids (Boyle et al., 1977). The release of riverine organic ligands may also contribute to the distribution of iron in shelf regions (Nakayama et al., 2011;Nishimura et al., 2012;Klunder et al., 2012).
Iron in surface waters is utilized by phytoplankton and is also subject to particle scavenging; it is recycled in the euphotic layer, and a part of the iron is exported to deeper waters as sinking particles. A fraction of the exported iron is then remineralized or desorbed from sinking particles and resupplied to surface waters by upwelling and mixing. In the climatological mean state, iron supplied by exogenous sources and by physical ocean processes is balanced by iron removal via biogeochemical processes.
In the HNLC regions, physical ocean processes play an important role in supplying new iron (not regenerated within the euphotic layer) to surface waters where exogenous sources are relatively limited owing to geographic distance from continental land masses. For example, in the eastern equatorial Pacific, physical ocean processes account for most of the supply of new iron to surface waters (Coale et al., 1996;Fung et al., 2000). The iron in this region derives mainly from sediments around Papua New Guinea and is transported eastward via the Equatorial Undercurrent (EUC) and then supplied to surface waters by upwelling (Mackey et al., 2002;Slemons et al., 2010;Kaupp et al., 2011). The amount of iron supplied is thus sensitive to circulation changes in the western equatorial Pacific. Ryan et al. (2006) suggested that the amount of iron transported is linked to the occurrence of El Niño events, because El Niño events regulate the intensity of the New Guinea Coastal Undercurrent (NGCU), which is the primary source of iron-enriched waters to the EUC.
Various observational and modeling studies have projected changes in meridional overturning and gyre-scale circulations and water column stratification in response to global warming. The Intergovernmental Panel on Climate Change (IPCC) Fourth Assessment Report (AR4) concluded that warmer climatic conditions enhance surface water stratification (Meehl et al., 2007), which can reduce the amount of iron supplied by vertical mixing (Steinacher et al., 2010). In the western equatorial Pacific, Coupled Model Intercomparison Project (CMIP) phase 3 (CMIP3) models have projected a substantial decrease of the South Equatorial Counter Current (SECC) flow (Ganachaud et al., 2013). Several studies have shown that increased radiative forcing causes southern subtropical gyres to spin-up and shift southward (Beal et al., 2011;Meijers et al., 2012). Such circulation and mixing changes are likely to alter the iron budget in HNLC waters.
In addition to the direct effect of physical ocean processes, the physical iron supply can also be modified via changes in the dissolved iron distribution. Reduced entrainment of macronutrients as a result of enhanced stratification can reduce biological iron utilization and particle scavenging in macronutrient-limited areas, and surplus iron can be transported laterally to downstream regions. This "downstream" effect can be important, because stratification preferentially reduces macronutrient concentrations in surface waters relative to iron, the supply of which depends largely on exogenous sources.
In this study, we used a CMIP phase 5 (CMIP5) model to investigate the iron budget in ocean surface waters, with a particular focus on HNLC waters. We evaluate future iron cycling changes and their impact on ocean productivity under one global warming scenario on the basis of our understanding of the iron budget of the present ocean. The model projections show that in the 2090s changes in ocean circulation will intensify the amount of iron supplied to HNLC waters, despite enhanced stratification of those waters. Thus, future oceanic iron-cycle changes will potentially enhance K. Misumi et al.: The iron budget in ocean surface waters in the 20th and 21st centuries 35 ocean productivity in HNLC waters and thus partly compensate for decreases in ocean productivity caused by reduced macronutrient entrainment.
In section 2, we describe the climate model used and the analytical method. We present the simulated results in Sect. 3. In Sect. 4, we compare the simulated iron cycle for the 1990s with other estimations and discuss future ocean iron budget changes in HNLC waters and their potential impact on ocean productivity. We summarize our conclusions in the final section.

Model
We used the Community Earth System Model (CESM1; Hurrell et al., 2013; available online at http://www.cesm.ucar. edu), a CMIP5 model developed for contributing to the IPCC Fifth Assessment Report (AR5). We used the CESM1-BGC model configuration (Lindsay et al., 2014); this configuration is mostly the same as the default configuration of the previously released Community Climate System Model version 4 (CCSM4; Gent et al., 2011). Different from the default CCSM4 configuration, the CESM1-BGC configuration includes a marine ecosystem module, three-dimensional atmospheric CO 2 tracers, and interactive coupling of the oceanic components to each other and to the land model. In addition, the CESM1-BGC configuration calculates the absorption of short-wave radiation in the upper ocean using the simulated chlorophyll field, whereas the default CCSM4 configuration calculates it by using a climatology based on satellite observations. The simulated physical states differ little, however, between the CESM1-BGC and CCSM4 configurations; thus, we refer readers to recently published CCSM4 papers for details (e.g., Gent et al., 2011;Danabasoglu et al., 2012). CESM1 consists of four major component models, the atmospheric, land, ocean, and sea-ice components, which are coupled interactively by a flux coupler (CPL7; Craig et al., 2011). The atmospheric and land components are the Community Atmospheric Model version 4 (CAM4; Neale et al., 2013) and Community Land Model version 4 (Lawrence et al., 2012), respectively. Both models adopt a uniform horizontal resolution of 1.25 • × 0.9 • , and the CAM4 has 26 vertical layers. The ocean component is the Parallel Ocean Program version 2 (POP2; Smith et al., 2010;Danabasoglu et al., 2012), and the sea-ice component is the Community Ice Code version 4 (Hunke and Lipscomb, 2008;Holland et al., 2012). These components adopt a non-uniform horizontal grid on which the northern pole is displaced onto Greenland. This grid has a zonal resolution of 1.125 • and a meridional resolution that varies from 0.27 • around the Equator, and the maximum resolution is about 0.64 • located in the northwestern Pacific. The POP2 has 60 vertical levels, which are 10 m thick in the upper 150 m; below 150 m their thickness increases with depth.
The BEC model considers two external iron sources: aeolian dust  and seafloor sediments (Moore and Braucher, 2008). Aeolian dust is assumed to contain a constant fraction of iron (3.5 wt %), and a small fraction, 2 %, of the dust iron is assumed to dissolve instantaneously at the sea surface. The model also considers the dissolution of dust iron as it sinks in the water column (Moore and Braucher, 2008). The BEC model uses the monthly mean climatology of aeolian dust and the annual mean climatology of the sedimentary iron flux. We assumed that these fluxes remained constant during the simulation period (from 1850 to 2100). Riverine and hydrothermal iron inputs and iron incorporated in sea ice or icebergs are not considered in this version of the BEC model.
The BEC model considers three pools of iron: dissolved iron, semilabile dissolved organic iron (DOFe), and particulate iron (PFe) pools. The dissolved iron pool represents total "dissolved" iron, including soluble and colloidal iron, iron complexed with organic ligands and free iron. Dissolved iron is assumed to be the only iron pool utilized by phytoplankton, and it is subject to particle scavenging. The ironscavenging rate is calculated by considering ambient sinking particle fluxes (particulate organic matter, dust, biogenic silica, and calcium carbonate) and, implicitly, the effect of complexation with binding ligands (for details see Moore and Braucher, 2008;and Misumi et al., 2011). Although in preliminary tests, we introduced spatially variable ligand concentrations into the BEC model ), a conventional approach with a homogeneous ligand concentration (0.6 nM) was used for this simulation. Iron in phytoplankton was routed to each iron pool via excretion, mortality, and sloppy feeding; the proportion allocated to each pool depended on the plankton class. For example, more iron is transferred to PFe via mortality and zooplankton grazing of diatoms than is transferred via the mortality and zooplankton grazing of smaller phytoplankton classes. DOFe is transported by physical ocean processes and is remineralized on a timescale of 100 days. PFe sinks instantaneously and is remineralized at depth in the same grid where it originates. The remineralization length scale is prescribed based on Armstrong et al. (2002) (for details see also Moore et al., 2004;Moore et al., 2008). Most scavenged iron (90 %) is routed to PFe, and the remainder is removed from the ocean, supposing transformation to an insoluble form and loss to the sediments. A change in phytoplankton community composition, therefore, affects iron-recycling rates in surface waters. For example, a shift of the community composition away from diatoms toward smaller phytoplankton results in a decrease in the amount of iron routed to PFe via mortality and zooplankton grazing, as well as a decrease in iron scavenging because of the reduced production of biogenic sinking particles.
In this study, we base our findings on results from two simulation experiments, a 20th century simulation and a 21st century simulation, using prescribed anthropogenic CO 2 emissions. Lindsay et al. (2014) have documented the simulation for the 1850-2005 period, called the 20C PROG experiment, and have described the initial conditions, spin-up procedure, and transient forcing for the 20th century simulation. For the simulation from 2005 to 2100, the model was run with the Representative Concentration Pathway 8.5 (RCP8.5) CO 2 emission scenario (Moss et al., 2010. This scenario assumes a nominal anthropogenic forcing of 8.5 W m −2 by 2100. The simulated data are available online at the CMIP5 website (http://cmip-pcmdi.llnl. gov/cmip5/). Moore et al. (2013) and Long et al. (2013) have documented the model skill in simulating ecology and biogeochemical tracers averaged over the 1990s in comparison with observational data sets. In this paper, we focus on the ocean iron budget averaged over two periods, the 1990s and the 2090s. This work complements the broader ecosystem analysis for these two time periods by Moore et al. (2013).

Iron budget analysis
The BEC model calculates the net dissolved iron tendency (TEND) or time rate of change of the dissolved iron concentration in each grid box as the sum of the following terms: where the terms represent simulated dissolved iron tendencies due to advection by large-scale flows (ADV); isopycnal, parameterized eddy mixing and surface boundary-layer mixing (MIX); biogeochemical processes (BGC) and external forcings (FRC). PHYS, that is, the sum of the ADV and MIX terms, represents the net physical ocean processes. The ADV term is subdivided into horizontal (ADVh) and vertical (ADVv) components. The difference between the 2090s and the 1990s is denoted by putting a "δ" before the term name. In this study, we consider the iron budget in the upper 100 m of the ocean, unless otherwise noted.
For each specific region, we define the iron budget amplitude (Amp) as follows: where terms i represents the individual iron budget terms (Eq. 1), and dV represents a specific region. The Amp ratio (rAmp) is defined as follows: where Amp 1990s represents the decadal average of Amp values in the 1990s. The ratio represents a magnitude of change in iron budget relative to the 1990s. Changes in advective iron supply rates can be caused by both the direct effect of a circulation change and the indirect effect of a change in the dissolved iron distribution. To diagnose the relative importance of the direct and indirect effects, we decompose the advection term. The advection term can be written as the iron flux convergence: where v represents three-dimensional velocity and Fe represents the dissolved iron concentration. The advection term in the 1990s can be written as where the box bracket [X] decade indicates the decadal monthly mean, the angle brackets X indicate the annual mean, and the prime superscript X denotes deviation from the decadal monthly mean of the monthly mean output of variable X. The difference in the advection term between the 2090s and the 1990s is, thus, The sum of the covariance terms is small in regions where this analysis was applied; thus, we neglect them. The first term of the last line in Eq. (6) can be rewritten as Then, by substituting Eq. (7) into Eq. (6), we can decompose δADV into three terms: The first term on the right-hand side represents the advective iron supply change due to circulation change, given the iron distribution in the 1990s, and the second term represents the iron supply change due to the iron distribution change, given the circulation in the 1990s. The last term represents the non-linear effect. By evaluating the relative importance of the first and second terms we can elucidate which processes are dominantly responsible for the net advective iron supply change. In the following, we denote the horizontal flux convergence as −∇ · (uFe) and the vertical flux convergence as −∂(wFe)/∂z.

Overview of the simulated physical and biogeochemical fields
The atmospheric CO 2 concentration in the 1990s was simulated to be 376 ppm, a value somewhat higher than the observed concentration in the 1990s (359.82 ppm; Ed Dlugokencky and Pieter Tans, NOAA/ESRL http://www.esrl. noaa.gov/gmd/ccgg/trends/). The simulated concentration increased to 1076 ppm by the 2090s (Fig. 1a). The globally averaged simulated sea surface temperature increased from 18.7 to 21.7 • C from the 1990s to the 2090s (Fig. 1b).
The simulated primary and export productions decreased from 56.1 to 52.9 Gt C yr −1 and from 8.07 to 7.00 Gt C yr −1 from the 1990s to the 2090s, respectively (Fig. 1c, d; Table 1). The simulated marine primary and export productions in the 1990s fell within the range of previous estimates (36.5-67 Gt C yr −1 for primary production; Longhurst et al., 1995;Antoine et al., 1996;Behrenfeld and Falkowski, 1997;Behrenfeld et al., 2005; 10 ± 3 Gt C yr −1 for export production; Sarmiento and Gruber, 2006). Although the globally averaged surface nitrate concentration was projected to decrease from 7.58 to 6.27 µM from the 1990s to the 2090s (Fig. 1e), the surface dissolved iron concentration was projected to increase from 0.69 to 0.75 nM (Fig. 1f). The model reproduced known basin-scale barotropic (vertically integrated) stream function (BSF) distributions well, but the results for the 1990s showed some biases in detail (Fig. 2a). The simulated Antarctic Circumpolar Current (ACC) transport through the Drake Passage was 170 Sv, roughly 25 % larger than the observational estimate (137 ± 8 Sv) (Cunningham et al., 2003). The simulated South Atlantic subtropical gyre transport was up to 58 Sv, which is more than double observed values (Large and Danabasoglu, 2006;Danabasoglu et al., 2012). In the 2090s, the overall patterns of the simulated BSF distributions are similar to those in the 1990s (Fig. 2b), but with moderate differences in some regions (Fig. 2c). In the Southern Hemisphere oceans south of 20 • S, the alternation of positive, negative, and positive values in the meridional direction imply that the southern subtropical gyres and the ACC shift southward. For example, the zonal average of the zero-contour latitude of the BSF in the Southern Ocean was projected to shift from 44.8 to 46.4 • S between the 1990s and the 2090s. The South Atlantic subtropical gyre transport was projected to increase by 40 % to 80 Sv in the 2090s.  The model simulated large-scale patterns of primary and export production rates compare well with satellite-based estimations (Fig. 3a, c; cf., Sarmiento and Gruber, 2006), though the model somewhat overestimated those rates in the subtropical oceans and underestimated them in the subarctic North Pacific. The simulated surface nitrate concentrations were also higher in the subtropical oceans and lower in the subarctic North Pacific ( Fig. 3e; see also Moore et al., 2013).
The simulated surface dissolved iron concentrations were overestimated, especially in the Indian and the equatorial Atlantic oceans, and they were underestimated in the Arctic Ocean (Figs. 3g, 4a, b). The overestimation is likely due to the underestimation of iron removal, because aggregation removal is not considered in our model , and the underestimation can be attributed to the fact that riverine inputs of neither iron nor organic ligands are considered (Nakayama et al., 2011;Nishimura et al., 2012;Klunder et al., 2012). The simulated iron concentrations were modestly overestimated in the Kuroshio extension region; this bias is likely attributable to the misplaced Kuroshio Current in our model (Misumi et al., 2011). In other areas, the model biases were relatively small, generally less than ±0.2 nM. Point-by-point comparison of surface dissolved iron concentrations between field observations and simulated data showed large dispersion (Fig. 4c), but the basin-scale distribution was reasonably simulated (Fig. 4d). The model successfully simulated the known HNLC regions (i.e., ironlimited areas for phytoplankton growth), in the eastern equatorial Pacific, the Southern Ocean and the subarctic North Pacific, though it underestimated the areal extent of the subarctic North Pacific HNLC region ( Fig. 5; see also . Although the globally integrated primary and export production rates decreased in the RCP8.5 simulation (Fig. 1c,  d), in iron-limited areas in particular, these rates were simulated to increase in the 2090s (Fig. 3b, d; Table 1). This result implies that changes in iron cycling fuel production in HNLC regions. Moore et al. (2013) has briefly described this response, and we describe it in greater detail. Surface nitrate concentrations were projected to decrease almost everywhere globally in the 2090s as a result of decreased entrainment from deeper waters due mainly to enhanced stratification ( Fig. 3f; see also Moore et al., 2013). The projected reduction of surface nitrate concentrations in nitrate-limited areas also led to a decrease in the biological uptake of iron and particle scavenging and resulted in increased surface dissolved iron concentrations in non-HNLC regions in the 2090s (Fig. 3h). The areal extent of the iron-limited areas was projected to shrink by the 2090s (Fig. 5).

Simulated oceanic iron budget in global surface waters
The values for the total TEND term (see Sect. 2.2) in the 1990s were very small compared with the values of its individual component terms (Fig. 6), a result that implies that surface water iron concentrations are in a quasi-steady state.
Simulated FRC values were positive in the global ocean, and they were especially large in the equatorial Atlantic, the northern Indian Ocean, and near coasts ( Fig. 6g) because of iron supplied from aeolian dust and continental margin sediments. Simulated PHYS values, which represent iron supplied via three-dimensional advection and mixing, were positive in most regions (Fig. 6c). Negative PHYS values mean that large external iron inputs exceed the amount of locally consumed iron, which is represented by the BGC term, and that the surplus iron is exported to adjacent areas or to deeper waters. Negative BGC values mean that iron is exported from surface waters by the biological pump and by particle scavenging (Fig. 6e). Globally, in the 1990s, the total amount of iron supplied to the upper 100 m accounted for by the FRC term was 40.8 Gmol yr −1 ; i.e., 8.29 Gmol yr −1 from aeolian dust and 32.5 Gmol yr −1 from sediments ( Table 2). The PHYS term accounted for only 4.09 Gmol yr −1 , roughly 10 % of total iron supply to the upper 100 m. Although the global integral was small, physical ocean processes were considered to play a key role in the surface iron budget because the PHYS term values exceeded the FRC term values in 73 % of the global ocean areas (Fig. 6c, g). The decomposition of the PHYS term shows that surface boundary-layer mixing was the most important physical process supplying iron to surface waters (60 % of the total PHYS value). The BGC term accounted for the removal of 44.9 Gmol yr −1 of iron from surface waters, an amount that balances the total iron inputs accounted for by the FRC and PHYS terms. The spatial patterns of the PHYS and BGC term differences between the 2090s and 1990s are similar to each other but with opposite sign (Fig. 6d, f). At a global scale, the net iron supply accounted for by the PHYS term decreased to 3.83 Gmol yr −1 in the 2090s (Table 2).
It should be noted that the PHYS values integrated globally over the upper 100 m (4.09 and 3.83 Gmol yr −1 for the 1990s and 2090s, respectively) represent the net vertical sup- include transport of iron from any origin. However, because the PHYS values are negative (flux divergence) in regions where physical processes receive iron from the other processes (iron inputs from the external sources, represented in the FRC term; and iron regenerated from organic matters and inorganic particles, included in the BGC term), the PHYS values integrated over a specific domain represents iron transport by net physical processes.

Iron supply processes to HNLC waters
The total iron supply to HNLC waters was 2.33 Gmol yr −1 in the 1990s, roughly two-thirds of which was accounted for by the PHYS term (1.60 Gmol yr −1 ; Table 2). Positive δPHYS values in HNLC areas (Fig. 6d) indicate elevated amounts of iron supplied via physical processes in the 2090s compared with the 1990s, whereas the corresponding negative BGC term values indicate intensified export of iron by the bio-logical pump and particle scavenging (Fig. 6f). Although the globally integrated PHYS term decreased in the 2090s, the value integrated over the HNLC areas increased from 1.60 to 1.81 Gmol yr −1 (Table 2).
We defined two HNLC regions for analysis by dividing the simulated iron-limited areas for diatoms in the 1990s according to latitude (Fig. 5c): the eastern equatorial Pacific (EE-QPAC), between 30 • S and 30 • N; and the Southern Ocean (SO) south of 30 • S. We excluded the subarctic North Pacific from our analysis because the simulated iron concentrations and limitation showed a relatively large bias (Fig. 4b) and the projected change in this region was more subtle.

Eastern equatorial Pacific
In the EEQPAC, in the 1990s, external iron inputs (FRC) are small (20 %) (Fig. 7a), and most iron is supplied by the ADVv component (65 %), that is, by equatorial upwelling. Much of this iron is transported by EUC and originates in the western equatorial Pacific . Large vertical iron fluxes from subsurface to surface waters are seen in the eastern equatorial Pacific (Fig. 8g). The vertical flux divergence (negative values) in the subsurface water (Fig. 8g) is supported by the zonal supply of iron (positive values) from the western region (Fig. 8c). A large zonal iron flux divergence at around 150 • E (Fig. 8c) corresponds to a large meridional iron flux convergence (Fig. 8e). The simulated horizontal iron transport for the 1990s shows large northward iron fluxes along the coast of Papua New Guinea (Fig. 9a), where iron inputs from shelf sediments (Fig. 9e) cause the divergence of these fluxes. The northward fluxes converge meridionally at around 150 • E in the equatorial region (Fig. 9a); as a result, a large meridional flux convergence is seen at around 150 • E (Fig. 8e). Thus, the iron supplied to the eastern equatorial Pacific originates mainly in shelf sediments of Papua New Guinea and is transported via the EUC.
In the 2090s, the simulated result showed increased ADVv and larger negative BGC values (Fig. 7b, c), the implication being that an elevated iron supply by vertical advection intensifies biological iron utilization and particle scavenging. The difference of the iron flux between the 2090s and the 1990s along the Equator suggests intensification of iron transport by a process that looks very similar to that seen in the 1990s (Fig. 8b, d, f, h). The elevated eastward iron transport (Fig. 8d) is triggered by a larger meridional flux convergence in the 2090s, observed at around 150 • E (Fig. 8f);   Fig. 4. Spatial distributions of the observed dissolved iron concentration (a) and the difference between the simulated results for the 1990s and the observed data (b) in the upper 100 m. The observed data were compiled by Tagliabue et al. (2012), with subsequent subsetting (see the "ALL-COAST-OUT" subset of Misumi et al., 2013). Scatter plots comparing the observed data with the simulated results in the 1990s in the upper 100 m (c). The simulated results are those for the same month sub-sampled at the grid points nearest to the measured locations. N, number of data; R 2 , coefficient of determination; and RMSE, root mean square error. The same plot but of basin-scale averages (d). Horizontal and vertical bars represent the standard deviations of the observed and simulated data, respectively, for the individual basins. The basin averages and the standard deviations of the simulated results were calculated using the sub-sampled data shown in (c). The background colors in (a) show the regions used for calculating basin averages. this increased meridional flux convergence can be attributed to the combined effect of circulation and production changes in this region (Fig. 9b, d; the mechanisms are described in detail in Sect. 3.4). The rAmp value in EEQPAC was 0.253 (Fig. 7c); thus, the simulated iron input to and output from the surface water of EEQPAC increases by roughly 25 % relative to that in the 1990s.

Southern Ocean
In the SO, in the 1990s, the largest positive value is seen in the MIX term (Fig. 7d); further decomposition of the MIX term revealed that surface boundary-layer mixing is the dominant contributor (55 % of the total MIX input). The second highest contributor to the iron supply in the SO is the FRC term (Fig. 7d). Although the net contribution to the iron budget is small, the ADVh term plays an important role in distributing iron over a wide area of the SO (Fig. 10a). The ad-vection fluxes diverge near the coast; this divergence indicates that the ocean currents along the coast entrain iron from external sources and offshore-transported iron by horizontal mixing (Fig. 10c). The advection fluxes penetrate the open ocean where the currents separate from the coasts (Fig. 10a). Some of the fluxes that cross the SO boundary converge downstream; these fluxes represent the advective iron supply to the SO. Large horizontal advection fluxes and their convergence are observed particularly in the Brazil-Malvinas confluence zone. Both the energetic current and the plentiful iron supply from external sources in the upstream region contribute to this extensive transport. Weakly positive ADVv values are broadly observed in the open ocean region (Fig. 10b); the values represent the iron supply by the mean flow of the meridional overturning circulation.
A modest reduction of the MIX term is seen in the 2090s for the SO (Fig. 7e, f). This reduction is mainly due to the shoaling of the mixed layer depth, and its impact is particularly large near the coast of Antarctica (Fig. 11c). The simulated ADVh increases in the 2090s (Fig. 7e, f), and this increase exceeds the reduction of the iron supply by the MIX term; therefore, the elevated iron supply to the SO in the 2090s is due to horizontal advection. The simulated iron input to/output from the surface water of the SO increases by roughly 7 % relative to that in the 1990s (Fig. 7f). The simulated spatial distribution of δADVh values shows intensification of iron transport in the Brazil-Malvinas confluence zone (Fig. 11a).

Mechanisms causing elevated advective iron transport to HNLC waters
Our model projected elevated advective iron transport in the 2090s both to the eastern equatorial Pacific and to the At-lantic sector of the Southern Ocean. This intensified advective iron transport can be explained either by the direct effect of a circulation change or by an indirect effect, driven by a change in the dissolved iron distribution, or by an interaction between direct and indirect effects. Here, we examine mechanisms intensifying the advective iron transport by using Eq. (8) to decompose δADV.

Eastern equatorial Pacific
The simulated δADV in the eastern equatorial Pacific was driven by changes in the western region. The simulated spatial pattern of δADVh in the western region (Fig. 12a) is similar to that of δPHYS (Fig. 9b), the implication being that the changes in the horizontal advection dominate the total physical iron supply change. The somewhat patchy structure of  δADVh compared to that of δPHYS is due to the exclusion of the mixing term, which modulates the concentration gradient developed by the advection term. The simulated positive and negative dissolved iron concentration anomalies (differences between the 2090s and the 1990s, green and magenta lines in Fig. 12) correspond well to the δ [u] [Fe] 1990s flux convergence and divergence, respectively (Fig. 12b). This argues that the circulation change is the primary cause of the dissolved iron concentration anomalies. The westward δ [u] [Fe] 1990s fluxes along with the flux divergence around 7 • S east of Papua New Guinea show the weakening of eastward iron transport by the SECC (Fig. 9a); the attenuation of the SECC intensifies northward iron transport along the coast and causes the positive and negative concentration anomalies in the western equatorial Pacific. We find eastward iron fluxes along the Equator east of 160 • E in the [u] 1990s δ [Fe] component (Fig. 12c), but not in the δ [u] [Fe] 1990s component (Fig. 12b); therefore, the intensified eastward iron transport along the Equator is driven by the developed west-east concentration gradient and not by a change in the EUC velocity. Decomposition of the vertical component (δADVv) using Eq. (8) shows the dominant contribution of the [w] 1990s δ [Fe] component of the iron supply to the surface water (Fig. 13). This result indicates that an increased vertical concentration gradient caused by the penetration of the positive iron concentration anomaly causes the increased iron supply to the surface waters of the eastern equatorial Pacific.
The influence of the production decrease in the western equatorial Pacific (Fig. 3b) on the development of the westeast concentration gradient is relatively small. We diagnosed the dissolved iron tendency caused by the production change by multiplying the net source/sink term for dissolved inorganic carbon by the Fe : C ratio (6.0 µmol Fe (mol C) −1 ), that is, the maximum iron uptake ratio for small phytoplankton and diatoms used in the BEC model (Fig. 14a). Weak positive tendencies overlap the positive concentration anomaly, but the values are much smaller than those caused by δADVh ( Fig. 12; note that the color scale in Fig. 14a is an order of magnitude smaller than that used in Fig. 12).
In summary, weakening of the SECC leads to a positive dissolved iron concentration anomaly in the western equatorial Pacific, and the resultant west-east concentration gradient promotes iron transport to the eastern region via the EUC. Penetration of higher dissolved-iron concentrations to subsurface waters enlarges the vertical concentration gradient along the Equator, and it intensifies vertical iron advection to the surface water in the eastern equatorial Pacific. Therefore, the direct effect of the circulation change and the associated concentration gradient change intensifies iron advection to the eastern equatorial Pacific.

Southern Ocean
In the Atlantic sector of the Southern Ocean, we found an intensified iron supply by horizontal advection in the 2090s. The simulated negative and positive dissolved iron  1990s (a, c, e and g) and δADV (b, d, f and h).
concentration anomalies in the subtropical region and the Brazil-Malvinas confluence zone correspond well to the iron flux divergence and convergence by δADVh (Fig. 15a); thus, the change in horizontal advection is mainly responsible for the concentration anomalies. But the fact that the decomposition of the δADVh term shows complicated spatial patterns in each component (Fig. 15b-d) suggests that there is an interaction among the components.
In the subtropical region, the δADVh flux divergence is explained well by the δ [u] [Fe] 1990s component, and the spatial pattern of the divergence corresponds well to the negative concentration anomaly (Fig. 15a, b). We also find large southward δ [u] [Fe] 1990s fluxes in the Brazil Current region. These results mean that the spin-up and southward shift of the South Atlantic subtropical gyre causes more iron to be exported to the downstream region via the Brazil Current. After the Brazil Current separates from the coast, meandering δ [u] causes alternating negative and positive −∇ · δ [u] [Fe] 1990s values where the direction of δ [u] is northward and southward, respectively (Fig. 15b), because [Fe] 1990s increases northward in this region (Fig. 3g). A similar spatial pattern but with opposite sign is observed in the −∇ · [u] 1990s δ [Fe] values in the Brazil-Malvinas confluence zone (Fig. 15c); this result implies that the  [u] 1990s δ [Fe] component modulates the concentration gradients caused by the circulation change. Therefore, changes in the velocity field are the main driver of iron transport from the subtropical region to the tongue-shaped region of the positive concentration anomaly.
East of 10 • W, however, the [u] 1990s δ [Fe] component shows a large flux convergence of the positive concentration anomaly (Fig. 15c), whereas the δ [u] [Fe] 1990s component does not (Fig. 15b). This result indicates that the concentration gradient contributes mainly to the transport of iron to the region further downstream. Because the positive concentration anomaly penetrates an iron-limited area, a downward iron concentration gradient is developed in the flow direction as a result of enhanced biological iron uptake in the downstream region. A modest influence of the non-linear component is observed locally (Fig. 13d), but the large-scale pattern is explained well by the interaction between the δ [u] [Fe] 1990s and [u] 1990s δ [Fe] components. Taken together, these results indicate that the intensification and southward shift of the South Atlantic subtropical gyre and the associated concentration gradient change causes the enhanced horizontal iron advection to the Atlantic sector of the Southern Ocean.
Reduced biological production in the subtropical region contributes locally to the intensification of advective iron flux transport ( Fig. 14b; note that the color scale is half that of Fig. 15), but the overall contribution is much smaller than that of the circulation change. If the "downstream" effect were dominant, then a positive dissolved iron concentration anomaly would be developed in the subtropical region and [u] 1990s δ [Fe] fluxes would transport iron from the subtropical region to the Brazil-Malvinas confluence zone. Given that the simulated dissolved iron concentrations decrease in the subtropical region and the direction of the simulated [u] 1990s δ [Fe] fluxes is opposite to the flow direction of the Brazil Current (Fig. 15c), the production decrease in the subtropical region makes only a minor contribution to the enhanced iron transport to the Southern Ocean in the 2090s.

Simulated iron transport processes and projected changes under future climate conditions
In the eastern equatorial Pacific, we found that most iron originates from shelf sediments around Papua New Guinea. Sedimentary iron is entrained by the western boundary current and then transported by the EUC to the eastern part of the region. Previous observational studies have well documented this west-to-east iron transport (Mackey et al., 2002;Ryan et al., 2006;Slemons et al., 2010) and we also briefly described such transport in our earlier study . The observational studies suggest that the NGCU plays a major role in this transport, but our model failed to simulate the NGCU because the Vitiaz Strait, the pathway of the NGCU between New Guinea and New Britain Island, is closed in our coarse-resolution model. Instead, a western boundary current flows northward along the eastern coast of New Britain Island and transports sedimentary iron to the EUC. Considering that our model simulates the EUC transport well (Danabasoglu et al., 2012), the simulated west-to-east iron transport processes seem realistic in this region, though the detailed transport path is different. Our model projected a weakening of the SECC under future climate conditions that results in an iron flux convergence in the western equatorial Pacific and an intensification of iron transport to the EUC and, thus, to the eastern equatorial Pacific. Consistent with our results, the multimodel ensemble of CMIP3 simulations has also projected substantial weakening of the SECC (Ganachaud et al., 2013). Consequently, it is likely that future ocean circulation changes will elevate the iron supply to surface waters of the eastern equatorial Pacific.  and 1990s (a-c). Boyd et al. (2012) summarized various iron supply mechanisms to the surface waters of the Southern Ocean. They suggest fluxes of the new iron supply ranging from O(1) µmol Fe m −2 yr −1 for the vertical diffusive supply to O(100) µmol Fe m −2 yr −1 for iron from sediment resuspension in waters less than 1000 m deep. If we convert the simulated spatial maps of iron supply to surface waters (Fig. 6c, g) to units of µmol Fe m −2 yr −1 by multiplying by 100 m, we find that the simulated physical iron supply is within the range of their estimate. Boyd et al. (2012) suggested that the largest iron supply to surface water of the Southern Ocean was from sediment resuspension, whereas our model showed that the largest contribution was from surface boundarylayer mixing around Antarctica. The continental shelf around Antarctica is steep and narrow; thus, most sedimentary iron is supplied below 100 m and is transported to the upper 100 m via surface boundary-layer mixing. Thus, the dominant iron supply process simulated in this study is consistent with the estimate of Boyd et al. (2012).
Despite the predicted enhanced stratification, our model projected an intensified iron supply to the Southern Ocean, based on basin-scale integration of the results. This intensification is driven by the spin-up and southward shift of the South Atlantic subtropical gyre, which enables more ironrich waters to intrude into the iron-limited area via the Brazil Current. Recent field data have shown a tendency for the subtropical front to shift southward as the climate warms (Beal et al., 2011, and references therein), and CMIP5 simulations commonly project the spin-up and southward shift of southern subtropical gyres under future climate conditions (Meijers et al., 2012). Thus, advective iron transport to the Southern Ocean is likely to intensify in the future climate. Reducing model biases in simulations of the southern subtropical gyres would make it possible to better quantify the impact of these changes.
Although the projected circulation change increased iron supply, the surface nitrate concentrations decreased both in the eastern equatorial Pacific and Southern Ocean in the 2090s (Fig. 3f). There are at least two reasons that could explain the different model response between nitrate and iron. First, water transported to the iron-limited area is deficient in nitrate relative to iron because the circulation change  transports water from the nitrate-limited area to iron-limited area. Second, intensification of western boundary currents enhanced entrainment of sedimentary iron near the coasts, while it does not for nitrate. Thus, future ocean circulation changes will affect individual nutrients differently owing to differences in nutrient limitations, and their source and sink.

Relative importance of the direct effect of a circulation change and the "downstream" effect in intensifying advective iron transport
Our model results show a decrease of nitrate concentrations in surface waters (Fig. 3f). Enhanced nitrate limitation reduces production in nitrate-limited areas (Fig. 3b, d) and promotes shifts in the phytoplankton community composition away from diatoms toward small phytoplankton . These changes reduce biological utilization of iron, decrease scavenging by biogenic particles, and elevate recycling rates of iron in the upper ocean. The combined effect increases dissolved iron concentrations by more than 0.1 nM in some open ocean regions (Fig. 3h). We previously inferred that the elevated iron concentrations contributed substantially to the elevated advective iron transport to adjacent HNLC areas . The detailed flux analysis presented here, however, revealed that this "downstream" effect is rather small compared with the Misumi et al., Fig. 13 Fig. 13. The same as Fig. 8h (a). Decomposition of δADVv using Eq. (8) (b-d).
direct effect of the circulation change in the equatorial Pacific and Atlantic sector of the Southern Ocean. The "downstream" effect has been documented by several modeling studies (Dutkiewicz et al., 2005;Moore et al., 2006;Tagliabue et al., 2008;Krishnamurthy et al., 2009). These studies argued the "downstream" effect in sensitivity experiments changing aeolian iron supply without considering changes in ocean dynamics, while this study focused more on the effects of changes in ocean dynamics. Changes in aeolian iron supply alter horizontal iron distributions in surface waters modifying lateral gradients and increasing the "downstream" effects. In this study, we demonstrated that the direct effect of the circulation changes will likely elevate iron supply to HNLC waters in the future. The combined effects of the direct circulation change influence and the "downstream" effects, and their relative importance should be investigated in future work.

Impact of iron cycle changes on primary and export production rates
An intensified iron supply to HNLC waters in the 2090s contributes to elevated primary and export production rates in the HNLC regions. The simulated primary and export production rates integrated over the HNLC regions increase by 1.8 Gt C yr −1 and 0.15 Gt C yr −1 , respectively, despite the decreased rate of production at a global scale (Table 1). It is important to note, however, that even in HNLC regions, factors other than iron availability, such as improvement of light and temperature conditions under warmer climatic conditions, may also contribute to elevated production rates (e.g., Steinacher et al., 2010).
To identify the contribution of the iron cycle changes to the changes in production, we calculated linear regression slopes and coefficients of determination (R 2 ) between the simulated temporal changes of dissolved iron supply to the upper 100 m accounted for by physical processes and primary and export production rates (Fig. 16). Where the production rates are notably enhanced by the iron supply increase, the regression slopes should be positive, and the R 2 values should be large. Here, we define the IRON-DRIVEN AREA as the area where the regression slopes are positive and R 2 values are 0.5 or more (dotted regions in Fig. 16). This provides a conservative estimate of where the HNLC production increases are mainly Misumi et al., Fig. 15   Fig. 15. The same as Fig. 12, but for the Brazil-Malvinas confluence zone. The maximum vector length is set to 20 µmol m −3 yr −1 , and advection vectors with magnitudes of less than 2 µmol m −3 yr −1 are omitted. Contour lines are the same as in Fig. 14b, but the zero contour is omitted. The gray dotted line represents the border of the simulated HNLC region in the 1990s. driven by increased iron inputs driven by changes in ocean circulation and mixing.
Although the areal extent of the IRON-DRIVEN AREA is not very large, it corresponds to those areas where the largest changes in production and export between the 1990s and the 2090s are simulated (Fig. 3b, d). Increases in primary and export production rates integrated over the IRON-DRIVEN AREA are 0.74 Gt C yr −1 and 0.08 Gt C yr −1 , respectively. Thus, in our simulation at least half of the total HNLC changes in production and export are caused by intensified iron inputs due to ocean circulation/mixing changes. The actual contribution of the ocean circulation/mixing changes would be higher, as there are extensive HNLC areas with positive regressions but R 2 values less than 0.5 (Fig. 16).
A previous ocean biogeochemical model intercomparison study showed that four models, one of which was our previous version (CCSM3), projected decreases in globally integrated primary and export production rates (Steinacher et al., 2010). Among the four models, only CCSM3 projected an increase in production in the eastern equatorial Pacific. This  increase of production in the eastern equatorial Pacific is also seen in our CESM1 result (Fig. 3b, d). We demonstrated that the elevated production is attributed to intensified sedimentary iron transport triggered by weakening of the SECC. Because two of the four models (MPIM and CSM1.4) did not consider sedimentary iron sources, they could not represent this response. The remaining model (IPSL) did consider sedimentary iron sources but nevertheless projected a decrease in production in the eastern equatorial Pacific (Steinacher et al., 2010). Comparing the simulated iron transport process under current climate and projected circulation changes between the models will provide insights into the different model responses.

Other factors changing the iron cycle in the future ocean
In this study, we investigated future oceanic iron-cycle changes by focusing on changes in ocean dynamics. Other factors, including changes in ocean biogeochemistry and in external iron inputs, may also alter the marine iron cycle in the future. An increase in anthropogenic carbon reduces ocean pH (Caldeira and Wickett, 2003;Orr et al., 2005;Doney et al., 2009c;Doney 2010), and ocean acidification changes iron speciation, but it is not known whether it increases or decreases bioavailable iron. A reduction in pH increases the thermodynamic solubility of Fe(III) near a pH of 8 (Kuma et al., 1996;Liu and Millero, 2002) and slows the oxidation rate of Fe(II) (Millero et al., 1987;Breitbarth et al., 2010). A four times atmospheric CO 2 experiment considering effects of dissolved oxygen, temperature, light and pH on iron speciation demonstrated that decrease in pH increases the fraction of Fe(II) and the pH effect overrides the others (Tagliabue and Völker, 2011). Ocean acidification also increases the relative reactivity of Fe(III) to Cu, which may relieve iron stress in phytoplankton (Hirose, 2007). In contrast, Shi et al. (2010) pointed out that ocean acidification may decrease the bioavailability of iron by elevating the binding strength of chelators with acidic binding groups. A re-cent modeling study considering iron speciation has demonstrated that changes in bioavailable iron in a more acidic ocean strongly depend on which speciation of iron we assume to be bioavailable (Tagliabue and Völker, 2011). Increases in seawater temperature and enhanced stratification decrease the dissolved oxygen content in the ocean interior (e.g., Keeling et al., 2010). Several studies have observed high dissolved iron concentrations, especially Fe(II), in oxygen minimum zones (OMZs) (Hopkinson and Barbeau, 2007;Moffett et al., 2007;Kondo and Moffett, 2013), and higher iron solubility and active biological reduction are regarded as plausible causes of the observed high dissolvediron concentrations. Thus, expansion of OMZs may increase dissolved iron concentrations in subsurface waters and hence increase iron concentrations in upwelling waters. Moreover, enhanced stratification may increase the amount of iron-binding ligands in subsurface waters by reducing photodegradation . Increases in seawater temperature and changes in the composition of phytoplankton communities will also affect both sources and sinks of iron-binding ligands.
To evaluate the net changes in iron biogeochemistry in response to global warming, these processes must be quantified and introduced together into numerical models. Most current-generation models simulating global-scale ocean biogeochemical cycling assume only a single "dissolved" iron pool and consider complexation with a single ligand having a uniform concentration Aumont and Bopp, 2006;Moore and Braucher, 2008;Galbraith et al., 2010). Recent efforts to incorporate iron speciation (Tagliabue and Völker, 2011) and iron-binding ligand dynamics (Tagliabue and Völker, 2011;Misumi et al., 2013) form the basis for the next generation of oceanic iron-cycle models.
The investigation of changes in exogenous iron inputs under a future warmer climate is also important. Model simulations have projected both an increase and a decrease in emissions of mineral aerosols under future climate states Tegen et al., 2004). A decrease in cloud water pH can increase soluble iron (Luo et al., 2008;Mahowald et al., 2009), and one observational study suggested that a large iron deposition event was facilitated by acidic fog (Iwamoto et al., 2011). Luo et al. (2008) pointed out that combustion sources make a large contribution to soluble iron; changes in anthropogenic iron sources may thus perturb marine productivity in the future (Krishnamurthy et al., 2009).
Because sea ice and icebergs play an important role in offshore iron transport in the Southern Ocean (Smith et al., 2007;Lannuzel et al., 2008;Lancelot et al., 2009), a reduction of their production rates will change the spatial pattern of the iron supply. Recent studies have suggested that ice sheets and melt water play a role in supplying iron to the ocean (Raiswell et al. 2008;Wadham et al. 2013;Bhatia et al., 2013); the inclusion of such sources in numerical models (e.g., Death et al., 2013) is potentially important because the high latitudes are experiencing rapid warming.

Conclusions
We investigated the iron budget in the upper ocean during the 1990s and the 2090s by analyzing the simulated results of Community Earth System Model version 1. Although there are some notable biases in the simulated results, comparisons with available field data indicate that the model captured large-scale features of oceanic circulation and dissolved iron distributions. Our results predict that physical transport strongly impacts the supply of iron to surface waters in HNLC regions; thus, changes in physical ocean processes will strongly impact future iron budgets in HNLC waters.
Although it is commonly assumed that enhanced ocean stratification in a future warmer climate will reduce nutrient entrainment from deeper waters, we found that intensified large-scale gyre circulation and shifts in their location may elevate the iron supply to HNLC waters. The simulated oceanic circulation changes in the equatorial Pacific and the Southern Ocean are generally consistent with recent observations under warming climate conditions and other simulated results projected by CMIP3 and CMIP5 models. Therefore, it is likely that future ocean circulation changes will elevate the iron supply to HNLC waters. We also demonstrated that changes in the iron supply due to physical ocean processes elevate primary and export production rates and have the potential to partially buffer any future production decrease caused by enhanced stratification. Better understanding of the responses of large-scale circulation to a warming climate is essential to make projections of future iron supply changes more robust. Future oceanic iron-cycle models should also take into account iron speciation and ligand dynamics. Investigations of changes in various exogenous iron sources under a future climate are also essential.