The fingerprint of climate variability on the surface ocean cycling of iron and its isotopes

The essential micronutrient iron (Fe) limits phytoplankton growth when dissolved Fe (dFe) concentrations are too low to meet biological demands. However, many of the processes that remove, supply, or transform Fe are poorly constrained, which limits our ability to predict how ocean productivity responds to ongoing and future changes in climate. In recent years, isotopic signatures (δ⁵⁶Fe) of Fe have increasingly been used to gain insight into the ocean Fe cycle, as distinct δ⁵⁶Fe endmembers of external Fe sources and δ⁵⁶Fe fractionation during processes such as Fe uptake by phytoplankton can leave a characteristic imprint on dFe signatures (δ⁵⁶Fe_diss). However, given the relative novelty of these measurements, the temporal scale of δ⁵⁶Fe_diss observations is limited. Thus, it is unclear how the changes in ocean physics and biogeochemistry associated with ongoing or future climate change will affect δ⁵⁶Fe_diss on interannual to decadal timescales. To explore the response of δ⁵⁶Fe_diss to such climate variability, we conducted a suite of experiments with a global ocean model with active δ⁵⁶Fe cycling under two climate scenarios. The first scenario is based on an atmospheric reanalysis and includes recent climate variability (1958–2021), whereas the second comes from a historical and high-emissions climate change simulation to 2100. We find that under recent climatic conditions (1975–2021), interannual δ⁵⁶Fe_diss variability is highest in the tropical Pacific due to circulation and productivity changes related to the El Niño–Southern Oscillation (ENSO), which alter both endmember and uptake fractionation effects on δ⁵⁶Fe_diss by redistributing dFe from different external sources and shifting nutrient limitation patterns. While the tropical Pacific will remain a hotspot of δ⁵⁶Fe_diss variability in the future, the most substantial end-of-century δ⁵⁶Fe_diss changes will occur in the Southern Hemisphere at middle to high latitudes. These arise from uptake fractionation effects due to shifts in nutrient limitation. Based on these strong responses to climate variability, ongoing measurements of δ⁵⁶Fe_diss may help diagnose changes in external Fe supply and ocean nutrient limitation.

1 Introduction

The micronutrient iron (Fe) is thought to control primary productivity in large parts of the global ocean where limited supply and/or rapid removal keep dissolved Fe (dFe) concentrations low. However, our understanding of the ocean Fe cycle is still limited, as it involves a multitude of internal cycling processes and various supply mechanisms, both of which are often poorly constrained (Boyd and Ellwood, 2010; Tagliabue et al., 2017). One way to help disentangle the web of processes and sources is to measure the isotopic signatures of Fe (namely ${\mathit{\delta }}^{\mathrm{56}}\mathrm{Fe}\left[\mathrm{‰}\right]=\left[{\left({}^{\mathrm{56}}\mathrm{Fe}{/}^{\mathrm{54}}\mathrm{Fe}\right)}_{\mathrm{sample}}/{\left({}^{\mathrm{56}}\mathrm{Fe}{/}^{\mathrm{54}}\mathrm{Fe}\right)}_{\mathrm{standard}}-\mathrm{1}\right]\cdot \mathrm{1000}$), as some Fe cycle processes distinctively alter δ⁵⁶Fe via fractionation and many external sources of Fe have characteristic δ⁵⁶Fe endmember signatures. Thus, δ⁵⁶Fe has been used to study various aspects of the ocean cycle. For instance, changes in δ⁵⁶Fe of the dissolved and particulate Fe pools have been observed for different phytoplankton bloom stages, indicating that phytoplankton may preferentially take up light Fe so that the δ⁵⁶Fe of the dFe pool (δ⁵⁶Fe_diss) becomes increasingly heavy (Ellwood et al., 2015). Consequently, heavy δ⁵⁶Fe_diss observed in other low dFe surface ocean systems were suggested to be due to ongoing phytoplankton uptake and biological recycling, possibly combined with fractionation during the complexation of Fe by organic ligands (Ellwood et al., 2020; Sieber et al., 2021). On the other hand, the distinct δ⁵⁶Fe endmembers of external sources have been used in mass balance approaches (of δ⁵⁶Fe_diss and dFe) to estimate the relative importance of each source in supplying new Fe (Conway and John, 2014; Pinedo-González et al., 2020). Such approaches take advantage of the wide range of source δ⁵⁶Fe endmembers from exceedingly light δ⁵⁶Fe associated with reductive sediments (as low as −3.3‰; Homoky et al., 2009; Severmann et al., 2010) and aeolian Fe from anthropogenic emissions (as low as −4.0 ‰, Kurisu et al., 2019) to crustal δ⁵⁶Fe (ca. 0.1 ‰) for dFe supplied by non-reductive sedimentary processes or dust deposition (Conway et al., 2019; Homoky et al., 2013; Radic et al., 2011; Waeles et al., 2007).

Overall, δ⁵⁶Fe_diss is likely determined by a combination of δ⁵⁶Fe fractionation during internal cycling and supply of new Fe from external sources with characteristic δ⁵⁶Fe endmembers (König et al., 2021, 2022). Both variable source endmembers and fractionation during phytoplankton uptake and organic complexation are needed for a global Fe isotope model to reasonably reproduce δ⁵⁶Fe_diss observations (König et al., 2021). However, most δ⁵⁶Fe_diss observations and the majority of dFe concentration measurements have been obtained from single occupations of stations as part of GEOTRACES transects and from a few process studies, one of which focuses on seasonal changes (Ellwood et al., 2015). Consequently, there is only limited observation-based information available about the temporal variability in δ⁵⁶Fe_diss (and Fe cycling in general) and its response to changes in ocean physics driven by climate variability. So far, only one study has reported repeat δ⁵⁶Fe_diss measurements from re-occupations of three different stations in the Atlantic (Conway et al., 2016). While at each station the shape of the δ⁵⁶Fe_diss profiles in this work was similar for both occupations, some discrepancies were detected for the station near Cabo Verde and the upper part of the profile in Cape Basin near South Africa. Variable δ⁵⁶Fe_diss near Cabo Verde were suggested to be related to changes in the relative contribution of dFe from different external sources, namely Fe supply by reductive sediments and dust dissolution. In the Cape Basin, changes in surface currents were thought to be responsible, as the differences in upper ocean δ⁵⁶Fe_diss and dFe coincided with changes in temperature, oxygen, and salinity. These observations suggest that δ⁵⁶Fe_diss can vary over interannual timescales, especially in the upper ocean. Substantial temporal variability in local Fe sources and cycling has also been observed for repeat measurements of Fe concentration at time series stations in the North Pacific, Atlantic, and Mediterranean (Bonnet and Guieu, 2006; Fitzsimmons et al., 2015; Nishioka et al., 2001; Schallenberg et al., 2015; Sedwick et al., 2005, 2020), but no parallel Fe isotope measurements were made.

Modelling work suggests that seasonal variability in dFe supply from aeolian deposition and winter mixing can induce variability in surface ocean δ⁵⁶Fe_diss both directly and indirectly (König et al., 2022). Direct effects occur due to the variable δ⁵⁶Fe endmembers of each source, and indirect effects also arise as differences in Fe supply can change the degree of δ⁵⁶Fe fractionation during phytoplankton uptake and/or complexation by organic ligands. The change in the degree of uptake fractionation is suggested to be greatest in areas where Fe was the limiting nutrient of primary production. However, while aeolian Fe deposition varied among the years in this study, changes in ocean physics driven by climate variability were not accounted for, as the same (seasonally variable) physical forcing fields were applied each year. We would expect climate-driven changes in ocean circulation to induce similar, if not larger, changes in δ⁵⁶Fe_diss than variable external Fe supply, as such changes can redistribute dFe and other nutrients more significantly. Changes in nutrient availability and other parameters (e.g. temperature) can then lead to shifts in local upper ocean biogeochemistry, potentially altering the effect of fractionating processes on δ⁵⁶Fe_diss, for instance during phytoplankton uptake (Ellwood et al., 2015). On interannual timescales, climate-related changes in ocean physics are mainly driven by climate variations such as the Southern Annular Mode or the El Niño–Southern Oscillation (ENSO), which has been shown to impact Fe cycling in the California Current System (Johnson et al., 1999). On longer timescales, global warming will likely lead to changes in ocean physics and, consequently, biogeochemistry, where the extent of such effects depends on current and future CO₂ emissions (Bindoff et al., 2019; Kwiatkowski et al., 2020), but their effects on Fe and Fe isotopes remain poorly constrained.

To test how surface ocean δ⁵⁶Fe_diss responds to interannual and decadal climate variability, we set up two sets of model experiments. In a set of hindcast simulations, variability in ocean physics was driven by the inherent variability in an atmospheric reanalysis product, covering the years 1958 to 2021. In a second set of experiments, we used outputs from a 249-year climate change simulation (1852–2100), which includes internal climate variability as well as effects of long-term global warming. These experiments allowed us to detect hotspots of surface ocean δ⁵⁶Fe_diss variability in the present and future climate and identify the driving forces behind such variability.

2 Methods

2.1 Experimental design

We conducted both sets of experiments (referred to as “hindcast” and “climate change”) by applying offline physical forcing fields to a version of the PISCES biogeochemical ocean model with active δ⁵⁶Fe cycling, which has been used previously to study the effect of internal cycling and (variable) external δ⁵⁶Fe endmembers on δ⁵⁶Fe_diss (König et al., 2021, 2022). This model is based on PISCES-v2 (Aumont et al., 2015) which represents two phytoplankton, two zooplankton, and two particle size classes (with variable particle reactivity; Aumont et al., 2017), five nutrients (NO₃, NH₄, PO₄, Fe, and Si), oxygen, and the carbonate system. Its Fe cycle representation is comparably complex (Tagliabue et al., 2016), and includes a prognostic Fe-binding ligand tracer (Völker and Tagliabue, 2015) and Fe input from four external sources (rivers, hydrothermal vents, dust deposition, sediments) plus sea ice (source or sink; Aumont et al., 2015). In the δ⁵⁶Fe version of this model, each Fe tracer is split into a heavy (⁵⁶Fe) and light (⁵⁴Fe) pool (König et al., 2021), as illustrated for heavy and light dFe in Text S1 in the Supplement. Distinct δ⁵⁶Fe endmembers are applied to each of the four external Fe sources and δ⁵⁶Fe fractionation is applied to two of the internal cycling processes: complexation by organic ligands (preference for heavy ⁵⁶Fe) and uptake by phytoplankton (preference for light ⁵⁴Fe; König et al., 2021; Text S1).

For the hindcast model experiments, we applied 5 d physical forcing fields obtained from coupled NEMO-SI³–PISCES simulations forced by the JRA-55 atmospheric reanalysis (Tsujino et al., 2018), as described in Buchanan and Tagliabue (2021), extended to 2021. We conducted three repeat cycles of the 64-year period (1830–2021). As changes in climate in the later years of the cycle may impact the earlier years of the last repeat cycle, we focus on the last 46 years of each simulation (1975–2021).

To investigate the impact of climate change on δ⁵⁶Fe_diss, we set up climate change simulations using offline forcing fields from the IPSL-CM5A climate model (Dufresne et al., 2013), as described in previous studies (e.g. Richon and Tagliabue, 2021). For each experiment (i.e. standard simulation and sensitivity tests), we first set up pre-industrial (PI) control simulations from 1801 to 2100 using forcings from an IPSL-CM5A experiment with fixed pre-industrial atmospheric CO₂. We then set up climate change experiments from 1852 to 2100, initialized using the 1851 PI control output and forced by fields from IPSL-CM5A experiments with atmospheric CO₂ concentrations following historical pathways until 2005 and switching to the high-emissions RCP8.5 scenario thereafter.

2.2 Sensitivity experiments and isolating δ⁵⁶Fe_diss components

We conducted an equivalent set of “standard” δ⁵⁶Fe model simulations and sensitivity tests for both hindcast and climate change physical settings (Table 1). For the standard simulations, we used the δ⁵⁶Fe model set-up described in König et al. (2021), whereas for the sensitivity experiments we turned off either uptake or complexation fractionation by setting their respective fractionation factors to 1. This allowed us to estimate the extent of δ⁵⁶Fe_diss variability that is caused by changes in fractionation and source endmember effects, using the following set of calculations.

Table 1 Overview of experiments.

Firstly, we “split up” δ⁵⁶Fe_diss into uptake fractionation (δ⁵⁶Fe_UF), complexation fractionation (δ⁵⁶Fe_CF), and endmember (δ⁵⁶Fe_EM) contributions, i.e. the effect which each fractionation process and the source endmembers have on “overall” δ⁵⁶Fe_diss. In the case of δ⁵⁶Fe_UF and δ⁵⁶Fe_CF, we calculated their respective contributions using the δ⁵⁶Fe_diss of the standard experiment and those obtained from sensitivity experiments where either the uptake fractionation factor (“noUF”) or complexation fractionation factor (“noCF”) was set to 1 (Eqs. 1 and 2).

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{\delta }}^{\mathrm{56}}{\mathrm{Fe}}_{\mathrm{UF}}={\mathit{\delta }}^{\mathrm{56}}{\mathrm{Fe}}_{\mathrm{diss}}-{\mathit{\delta }}^{\mathrm{56}}{\mathrm{Fe}}_{\mathrm{diss},\mathrm{noUF}}\text{(2)}& {\displaystyle }{\mathit{\delta }}^{\mathrm{56}}{\mathrm{Fe}}_{\mathrm{CF}}={\mathit{\delta }}^{\mathrm{56}}{\mathrm{Fe}}_{\mathrm{diss}}-{\mathit{\delta }}^{\mathrm{56}}{\mathrm{Fe}}_{\mathrm{diss},\mathrm{noCF}}\end{array}$$

Thanks to the additive nature of fractionation and endmember effects in our model, which we confirmed for the hindcast experiments (Fig. S1 in the Supplement), the endmember effect δ⁵⁶Fe_EM could be calculated by subtracting the two fractionation effects from δ⁵⁶Fe_diss (Eq. 3).

$$\begin{array}{}\text{(3)}& {\mathit{\delta }}^{\mathrm{56}}{\mathrm{Fe}}_{\mathrm{EM}}={\mathit{\delta }}^{\mathrm{56}}{\mathrm{Fe}}_{\mathrm{diss}}-{\mathit{\delta }}^{\mathrm{56}}{\mathrm{Fe}}_{\mathrm{UF}}-{\mathit{\delta }}^{\mathrm{56}}{\mathrm{Fe}}_{\mathrm{CF}}\end{array}$$

We also set up an experiment (“neuSED”) where the sediment δ⁵⁶Fe endmember was set to 0 ‰ (only for the hindcast physical setting) to study endmember effects in areas where the surface ocean sediment endmember of the model (−1 ‰; König et al., 2021) is likely too light (Sect. 4.1).

2.3 Determining interannual variability in δ⁵⁶Fe_diss and its components

To quantify the interannual variability in δ⁵⁶Fe_diss and its three components (δ⁵⁶Fe_UF, δ⁵⁶Fe_UF, and δ⁵⁶Fe), we first applied a 12-month running mean boxcar filter to each value (which were based on monthly model output) to smooth out seasonal effects. We then calculated the interannual variance (VAR) and standard deviation (SD) for δ⁵⁶Fe_diss and for each component over the respective time periods of interest (1975–2021 for the hindcast and 2006–2100 for the climate change experiments). Note that we focus here on the surface ocean (0–10 m), as δ⁵⁶Fe_diss variability beneath the surface is often driven by small lateral or vertical movement in steep δ⁵⁶Fe_diss gradients, for instance, due to shifts in mixed layer depth. To compare the interannual variability in each of the three components to δ⁵⁶Fe_diss variability, we calculated the “relative variability”, i.e. the ratio of their VAR to that of δ⁵⁶Fe_diss, where subscript i = UF, CF, or EM (Eq. 4).

$$\begin{array}{}\text{(4)}& {\mathit{\delta }}^{\mathrm{56}}{\mathrm{Fe}}_i\text{relative variability}={\displaystyle \frac{\mathrm{VAR}\left({\mathit{\delta }}^{\mathrm{56}}{\mathrm{Fe}}_i\right)}{\mathrm{VAR}\left({\mathit{\delta }}^{\mathrm{56}}{\mathrm{Fe}}_{\mathrm{diss}}\right)}}\end{array}$$

Note that the sum of the three ratios is less (greater) than 1 in areas with positive (negative) covariance between two or three of the components, i.e. where their responses to climate variability are reinforcing (opposing) each other.

3 Results & discussion

3.1 Hotspots of present and future δ⁵⁶Fe_diss variability in the surface ocean

3.1.1 Areas of high δ⁵⁶Fe_diss variability in the present climate (1975–2021)

In the 1975–2021 period of the hindcast simulations, interannual surface ocean δ⁵⁶Fe_diss variability is highest in the tropical Pacific due to multiple drivers. This is illustrated by maxima in the interannual δ⁵⁶Fe_diss SD, with additional areas of localized elevated variability present around 40^∘ S (Fig. 1a). Interestingly, there appears to be no obvious connection between the interannual δ⁵⁶Fe_diss SD and the interannual SD of dFe concentration and primary productivity (Fig. S2), indicating that (absolute) changes in these values are not (solely) driving δ⁵⁶Fe_diss variability. Comparing this δ⁵⁶Fe_diss SD to the average δ⁵⁶Fe_diss values (Fig. 1b) indicates that the δ⁵⁶Fe_diss SD is often elevated in areas with steep horizontal δ⁵⁶Fe_diss gradients, with shifts in circulation causing local changes in δ⁵⁶Fe_diss. Elevated variability in δ⁵⁶Fe_UF and δ⁵⁶Fe_EM (Fig. 1c, d) suggests that the high δ⁵⁶Fe_diss variability in the tropical Pacific during our simulations is due to changes in both uptake fractionation and the circulation of external source δ⁵⁶Fe endmembers. At high latitudes, a combination of changes in uptake and, to a lesser extent, complexation fractionation appear to drive δ⁵⁶Fe_diss variability (Fig. 1c, e), although variability is generally lower in these areas (Fig. 1a). The elevated δ⁵⁶Fe_diss SD at around 40^∘ S is caused by horizontal shifts in areas with heavy δ⁵⁶Fe_diss likely related to strong (accumulated) uptake fractionation effects (Fig. 1c), as discussed in Sect. 4.2. Note that covariance between the components can be substantial for some regions, especially between δ⁵⁶Fe_UF and δ⁵⁶Fe_EM (Fig. 1f), which can increase or decrease the δ⁵⁶Fe_diss variability, depending on location (Fig. S3). Furthermore, δ⁵⁶Fe_diss variability (unsurprisingly) increases when seasonal effects are accounted for (by calculating δ⁵⁶Fe_diss SD from monthly model outputs without applying a 12-month running mean; Fig. S4). Seasonal effects are thereby largest at high latitudes mostly due to uptake fractionation effects, especially in the Southern Hemisphere (south of ca. 40^∘ S). However, a combination of endmember changes due to ocean circulation and uptake fractionation remain important in the tropical Pacific (Fig. S4).

Figure 1 Surface ocean δ⁵⁶Fe_diss interannual variability and its drivers in the present climate (1975–2021). Surface ocean (0–10 m) (a) interannual δ⁵⁶Fe_diss SD (‰) and (b) average δ⁵⁶Fe_diss (‰) of the hindcast standard experiment. Ratios of interannual variability in (c) δ⁵⁶Fe_UF, (d) δ⁵⁶Fe_EM, and (e) δ⁵⁶Fe_CF to interannual δ⁵⁶Fe_diss variability (Eq. 4); drivers are denoted as “endmembers” (EM), “uptake fractionation” (UF), “complexation fractionation” (CF), and various combinations thereafter. (f) Covariance between δ⁵⁶Fe_UF and δ⁵⁶Fe_EM (12-month running mean).

3.1.2 Future δ⁵⁶Fe_diss variability (2006–2100) due to climate change

The major changes in ocean circulation and physical conditions in our high-emissions climate change experiments lead to substantial δ⁵⁶Fe_diss variability over the next century. The interannual δ⁵⁶Fe_diss SD over the 2006–2100 period of our climate change simulation is considerably higher than for the parallel PI control simulation without any climate change effects (and higher than for the hindcast simulation, Fig. 2a and b). Over the shorter period of the hindcast experiments (1975–2021), elevated δ⁵⁶Fe_diss SD is mainly due to temporal variability around a mean δ⁵⁶Fe_diss value, whereas for the climate change experiments, elevated δ⁵⁶Fe_diss SD is also related to a change in the mean δ⁵⁶Fe_diss over the next century (Sect. 3.3). The δ⁵⁶Fe_diss SD of the PI control thereby resembles that of the standard hindcast experiment (Fig. 1a; note difference in scale), with some discrepancies due to the differences in model circulation.

Figure 2 Surface ocean δ⁵⁶Fe_diss interannual variability under future climate scenarios (2006–2100). Surface ocean (0–10 m) interannual SD δ⁵⁶Fe_diss (‰) for (a) climate change and (b) PI control experiments, and (c) the difference between the two (i.e. climate change SD δ⁵⁶Fe_diss − PI control SD δ⁵⁶Fe_diss; ‰). Ratios of interannual variability in (d) δ⁵⁶Fe_UF, (e) δ⁵⁶Fe_CF, and (f) δ⁵⁶Fe_EM to interannual δ⁵⁶Fe_diss variability (Eq. 4).

In the tropical Pacific, where interannual δ⁵⁶Fe_diss SD was highest over the period covered by the hindcast experiment (1975–2021; Fig. 1a), variability remains high under future climate change conditions (Fig. 2a). However, while the magnitude of δ⁵⁶Fe_diss SD is similar for the climate change experiment, there are some regional differences. δ⁵⁶Fe_diss SD is elevated in the eastern equatorial Pacific compared to the PI control but lower in the western part and subtropical regions (Fig. 2c). Generally, a combination of both endmember (δ⁵⁶Fe_EM) and uptake fractionation (δ⁵⁶Fe_UF) effects appears to be responsible for the δ⁵⁶Fe_diss variability in the tropical Pacific (Fig. 2d, f), similar to the hindcast experiment (Fig. 1f); however, in the eastern equatorial Pacific where δ⁵⁶Fe_diss variability increased, δ⁵⁶Fe_EM appears to be the controlling factor (see also Sect. 3.2.3).

A notable increase in interannual δ⁵⁶Fe_diss variability due to climate change occurs in the Southern Hemisphere. In the region between ca. 30 and 50^∘ S (Fig. 2a, b), where δ⁵⁶Fe_diss SD was already elevated under present climate conditions (hindcast simulation; Fig. 1a), δ⁵⁶Fe_diss SD more than doubles due to climate change (Fig. 2c). While some of this variability seems to be related to endmember effects (δ⁵⁶Fe_EM), changes in uptake fractionation (δ⁵⁶Fe_UF) appear to dominate the future δ⁵⁶Fe_diss variability in most of the Southern Hemisphere at middle to high latitudes (Fig. 2d). This is due to changes in δ⁵⁶Fe_UF between 2006 and 2100 caused by the impact of climate change on primary production (Sect. 3.3).

3.2 δ⁵⁶Fe_diss variability in the tropical Pacific

3.2.1 Current (1975–2021) δ⁵⁶Fe_diss variability in the tropical Pacific driven by ENSO variability

In the tropical Pacific, where δ⁵⁶Fe_diss variability in the hindcast simulation is highest (Fig. 1), there is a clear link between surface ocean δ⁵⁶Fe_diss and changes in sea surface temperature (SST) associated with different ENSO phases. This can be illustrated by focussing on the equatorial region between 5^∘ S and 5^∘ N (Fig. 3). As the hindcast simulations were forced by an atmospheric reanalysis product (see Sect. 2), the modelled changes in SST agree very well with observations in the equatorial Pacific, including the timing and extent of the El Niño and La Niña transitions (Fig. S5).

Figure 3 Variability in SST as well as δ⁵⁶Fe_diss and its drivers in the equatorial Pacific from the hindcast experiment. Time series (1975–2021) of monthly mean surface ocean (a) SST anomaly (^∘C) and Ocean Niño Index (ONI; red: El Niño, grey: La Niña), (b) SST (^∘C), (c) δ⁵⁶Fe_diss (‰), (d) δ⁵⁶Fe_EM (‰), and (e) δ⁵⁶Fe_UF (‰), averaged from 5^∘ N to 5^∘ S. The SST anomaly was calculated by subtracting the 1975–2021 average SST from the monthly outputs. The Ocean Niño Index uses the same key as the SST anomaly and was calculated from smoothed SST anomalies (3-month running mean) of the ENSO 3.4 region (120–170^∘ W, 5^∘ N–5^∘ S).

In neutral or weak El Niño/La Niña phases, when SST anomalies and Ocean Niño Index (ONI) are less than ca. ± 0.5 ^∘C (Fig. 3a), δ⁵⁶Fe_diss is light ($<-\mathrm{0.2}$ ‰) close to the continental margins of Papua New Guinea (PNG) in the west and South America in the east and heavier in the offshore regions of the open ocean (up to >1 ‰; Fig. 3c). However, during an El Niño event, the area of light surface ocean δ⁵⁶Fe_diss close to the PNG margin expands eastward (Fig. 3c), with similar timing and extent as the characteristic warming signal observed in the central and eastern equatorial Pacific (Fig. 3a, b). While these changes are visible in all El Niño years, they are most prominent for the very strong El Niño events of 1982/1983, 1997/1998, and 2015/2016. During these three extreme events (dark red ONI), surface ocean δ⁵⁶Fe_diss becomes anomalously light in most areas along the Equator, except in the easternmost part (around 120^∘ W), where it is heavier (Fig. 3c). For weaker El Niños (light red ONI) and those where the maximum SST anomalies occur in the centre of the equatorial Pacific, the light δ⁵⁶Fe_diss anomaly is restricted to the western part of the equatorial Pacific (Fig. 3c). Conversely, La Niñas events (grey ONI) are associated with a δ⁵⁶Fe_diss decrease in the eastern part of the basin, albeit to a lesser degree than for strong El Niños, as the light δ⁵⁶Fe_diss area close to the South American margin expands westward (Fig. 3c).

3.2.2 Mechanisms behind ENSO-driven δ⁵⁶Fe_diss variability

Splitting δ⁵⁶Fe_diss into its main drivers in the tropical Pacific, δ⁵⁶Fe_EM and δ⁵⁶Fe_UF, shows that δ⁵⁶Fe_diss variability during ENSO cycles is due to a combination of endmember and uptake fractionation effects (Fig. 3d, e). These δ⁵⁶Fe_EM and δ⁵⁶Fe_UF changes are due to circulation changes which affect the redistribution of Fe from external sources and lead to changes in primary productivity and, consequently, Fe uptake rates. We illustrate this with a schematic of the main mechanisms involved (Fig. 4), and use the time period 1995–1996 as an example, as it includes a very strong El Niño followed by a strong La Niña (Fig. 5). During El Niño events, such as in 1997, the weakened trade winds lead to a reversal in the upper ocean currents from predominantly westward flow to eastward flow in most equatorial regions (Figs. 4a, b; 5h). This reversal has far-reaching effects on upper ocean biogeochemistry and Fe cycling, such as higher dFe concentrations in the west, lower concentrations in the east (Fig. 5d), and a general decrease in the ratio of Fe uptake to dFe concentration (Fig. 5e), Fe limitation (Fig. 5f), and primary productivity (Fig. 5g). These changes, in turn, impact δ⁵⁶Fe_EM (Fig. 5b) and δ⁵⁶Fe_UF (Fig. 5c), as discussed below. Eventually, as the El Niño event breaks down, the incoming La Niña leads to westward currents across the basin (Fig. 4c), which are even stronger than in the neutral phase that precedes El Niño (Fig. 4a), as exemplified by the transition in 1998 (Fig. 5). Note that the same mechanisms also operate during El Niños or La Niñas that are weaker than those of 1997/1998 but to a lesser extent, leading to reduced changes in the parameters in question (Fig. S6). Furthermore, similarly decreasing trends in dFe concentration and primary productivity as simulated by our model for the equatorial Pacific (Fig. 5) have also been observed in the California Current System during the 1997/1998 El Niño, where coastal upwelling was also suppressed (Johnson et al., 1999).

Figure 4 Overview of processes behind δ⁵⁶Fe_diss changes in the equatorial Pacific during an El Niño/La Niña cycle. Note that the example δ⁵⁶Fe_diss (‰), dFe concentration (µmol m⁻³), and primary productivity (molC m⁻³ yr⁻¹) curves were taken from the same months as in Fig. 5, i.e. July 1996 for neutral, December 1997 for El Niño, and December 1998 for La Niña phases, using average values from 5^∘ S to 5^∘ N. The dFe and macronutrient limitation plots were taken from the same months and indicate if dFe (red) or macronutrients are limiting phytoplankton growth in (parts of) the 5^∘ S to 5^∘ N region.

Figure 5 Mechanisms behind surface ocean δ⁵⁶Fe_diss changes in the equatorial Pacific from the hindcast experiment during a strong El Niño/La Niña cycle. Upper half: time series (1995–1999; from hindcast experiments) of monthly mean surface ocean (a) δ⁵⁶Fe_diss (‰), (b) δ⁵⁶Fe_UF (‰), (c) δ⁵⁶Fe_UF (‰), (d) dFe concentration (µmol m⁻³), (e) ratio of Fe uptake to dFe concentration, (f) Fe limitation, (g) primary productivity (molC m⁻³ yr⁻¹), and (d) upper ocean (0–50 m average) eastward currents (m s⁻¹), averaged from 5^∘ N to 5^∘ S. The Ocean Niño Index is included on the left (red: El Niño, grey: La Niña). Lower half: map view of the same parameters in July 1996 (neutral), December 1997 (strong El Niño), December 1998 (strong La Niña).

The reversal of upper ocean currents during El Niño/La Niña cycles has a direct impact on δ⁵⁶Fe_EM. In neutral phases before the on-set of El Niño (e.g. July 1996), waters with elevated dFe concentration are upwelled at the South American margin and transported westward (Fig. 5d, h). The δ⁵⁶Fe endmember signature of this dFe is light, as it is mostly sourced from reducing sediments at the South American margin (Fig. 5b). As the currents reverse during an El Niño event such as in 1997, the upwelling and transport at the eastern margin is reduced and instead, high dFe waters with light δ⁵⁶Fe_EM from margin sediments in the PNG region are preferentially transported eastwards (Fig. 5b, d). Thus, δ⁵⁶Fe_EM is light in the west but remains near a crustal signature (i.e. ca. 0.1 ‰) in the eastern tropical Pacific, where the main source of external dFe is now dust deposition. Finally, as the El Niño conditions collapse and La Niña develops (as in 1998), strong westward currents lead to a similar, but larger, light δ⁵⁶Fe_EM anomaly from the South American margin as in the neutral phase, which, in the case of the strong El Niño/La Niña cycle of 1997/1998, extends across the basin (Fig. 5b). Note that the sedimentary Fe input from the PNG region may not be as isotopically light as suggested by our model so that the endmember effects in the western part of the basin may be reduced, as discussed in Sect. 4.1. Also, the sedimentary Fe supply from the PNG region to the equatorial undercurrent (EUC) may be underestimated in our model compared to some (but not all) of the available regional observations (Fig. S7). As the EUC supplies Fe to the surface ocean (e.g. Coale et al., 1996; Kaupp et al., 2011), this may cause additional impacts on surface ocean dFe concentrations and δ⁵⁶Fe_diss as well as their response to climate variability, which affects the strength of the EUC (e.g. Firing et al., 1983; Johnson et al., 2000).

The impact of circulation changes associated with the El Niño/La Niña cycle on δ⁵⁶Fe_UF depends on the ratio of Fe uptake to dFe concentrations. Impacts of El Niño/La Niña on uptake fractionation integrate the changes in primary production, limiting nutrients, and, consequently, Fe uptake rates by phytoplankton. Uptake fractionation becomes strongest in our model when Fe uptake rates are high relative to dFe concentrations (Fig. 5e). This leads to heavy δ⁵⁶Fe_UF. When the ratio of Fe uptake to dFe concentration is high, the system is usually Fe limited (Fig. 5f). During the neutral phases such as in summer 1996, such “high Fe uptake, low dFe concentration” conditions with heavy δ⁵⁶Fe_UF are prevalent in the central and eastern tropical Pacific (Fig. 5c, e), which are relatively productive and strongly Fe limited (Fig. 5f, g). As El Niño develops, the decreased upwelling and westward transport of macronutrients (e.g. nitrate) from the South American margin leads to a decrease in productivity and a switch from Fe to nitrogen limitation in the central and eastern Pacific (Fig. 5f, g). When combined with the input of dFe from the western margin, this shift of nutrient limitation leads to low ratios of Fe uptake to dFe concentration (Fig. 5e), and therefore lower uptake fractionation effects (i.e. relatively light δ⁵⁶Fe_UF) across the basin. During the following La Niña, dFe concentrations are relatively high in the Fe-limited eastern and central Pacific (Fig. 5d). These elevated dFe concentrations are due to decreased Fe uptake during the nitrogen-limited El Niño phase and the increased westward transport of dFe from the South American margin during La Niña. Consequently, the ratio of Fe uptake to dFe concentration is relatively low and uptake fractionation effects moderate (Fig. 5c, e), despite high primary productivity fuelled by the increased upwelling of nitrate and dFe (Fig. 5g).

3.2.3 Future (2006–2100) δ⁵⁶Fe_diss variability in the tropical Pacific driven by climate change

The climate change-induced decrease in surface ocean δ⁵⁶Fe_diss variability in the western equatorial Pacific and increase in the eastern equatorial Pacific over the next century, shown in Fig. 2, is driven by decreased upwelling of cold water in the east (Figs. S8, 6, 7). This reduces the input and westward transport of both macronutrients and isotopically light sedimentary dFe from the South American margin to the surface ocean and decreases primary production. These changes are all similar to those occurring during El Niño events from our hindcast simulations. Consequently, the western half of the equatorial Pacific becomes nitrogen limited and primary productivity is low, leading to roughly constant uptake fractionation, i.e. invariable δ⁵⁶Fe_UF. Conversely, the effects of different ENSO phases are now concentrated in the eastern equatorial Pacific, where variability is high for both δ⁵⁶Fe_EM and δ⁵⁶Fe_UF, due to changes in upwelling of light dFe and variability in the degree of Fe limitation (or switches to nitrogen limitation) respectively.

Importantly, circulation patterns during the historical period (1852–2005) in the climate change simulations are not as realistic as for the hindcast experiments, which were forced by an atmospheric reanalysis. For instance, the cold SST anomaly at the Equator is too strong and extends too far west in the climate change simulations, which is a common bias in such climate models (e.g. Planton et al., 2021). This “cold tongue bias” is caused by upwelling in the eastern Pacific which is too strong in the historical era of climate change experiment. This is also evident when comparing the equatorial SST of the climate change and hindcast simulations and is likely also responsible for the lower δ⁵⁶Fe_EM variability in the western tropical Pacific in the climate change simulation (Fig. S8). In addition, there is some disagreement between observations and the ENSO cycles produced by the climate change model (in the historical period of 1982–2005), for instance, regarding seasonal timing and zonal pattern (Bellenger et al., 2014; Planton et al., 2021). As models which better agree with historical ENSO characteristics predict the frequency of extreme ENSO events to increase (Cai et al., 2015, 2021), future variability in δ⁵⁶Fe_diss may be higher than predicted by our climate change experiment.

While the physical model shows some biases in terms of SST and ENSO characteristics, the modelled primary productivity changes in the tropical Pacific in response to changing SST during ENSO cycles agree well with available constraints from SST observations and satellite estimates of primary production (Kwiatkowski et al., 2017; Tagliabue et al., 2020). This strengthens our confidence in the effect of ENSO cycles on δ⁵⁶Fe_UF, as it indicates that the impact of the different ENSO phases on the model biogeochemistry is realistic. It should be noted that the physical setting of the climate change experiments is based on model simulations using the high-emissions RCP8.5 scenario. As smaller changes in biogeochemical and physical conditions are predicted for lower emission scenarios (Dufresne et al., 2013; Kwiatkowski et al., 2020), we would expect the δ⁵⁶Fe_diss changes also to be reduced in parallel.

3.3 Impact of climate change on surface ocean δ⁵⁶Fe_diss in the Southern Hemisphere

Future climate change-induced alterations to δ⁵⁶Fe_diss become largest in the Southern Hemisphere. The largest changes occur in a region between about 30 and 50^∘ S (Figs. 2, 6a), where stronger uptake fractionation effects (i.e. heavier δ⁵⁶Fe_UF; Fig. 6b) lead to increasingly heavy δ⁵⁶Fe_diss in the Pacific and Indian oceans (Fig. 6a). These changes in uptake fractionation are ultimately caused by higher primary production rates (Fig. 7c), which lead to stronger Fe limitation (Fig. 7b) and thus higher δ⁵⁶Fe_UF (Fig. 6b). The increases in δ⁵⁶Fe_UF between 30 and 50^∘ S are reinforced by an increase in δ⁵⁶Fe_EM in the western parts of each basin (Fig. 6c) due to a southward shift of low dFe areas with relatively heavy δ⁵⁶Fe_EM (Figs. 6c, 7a). Conversely, the δ⁵⁶Fe_diss decrease at lower latitudes is due to lower primary productivity and Fe limitation (Fig. 7b, c), and thus lighter δ⁵⁶Fe_UF (Fig. 6b). A notable decrease in δ⁵⁶Fe_diss thereby occurs just north of the heavy δ⁵⁶Fe_diss area due to a decrease in Fe limitation as the region of high, Fe-limited primary productivity moves southward. Similarly, the strong decrease in δ⁵⁶Fe_UF (and thus δ⁵⁶Fe_diss) in the South Atlantic is due to a decrease in Fe limitation and, in some parts, a shift to macronutrient limitation. Many of the changes in dFe and NO₃ concentration, primary productivity, and Fe limitation (Fig. 6) that drive δ⁵⁶Fe_diss changes over the next century were also observed in the only other available study on the impact of climate change on the ocean Fe cycle (Misumi et al., 2014). For instance, these authors also report increased future primary productivity at around 40^∘ S, which is Fe limited in the Indian and Pacific but no longer in the Atlantic (Misumi et al., 2014). The slight difference in depth range and time periods used in Fig. 6 in this study relative to those reported by Misumi et al. (2014) has little impact on the trends of each parameter. One major difference between the two simulations concerns primary productivity in the eastern tropical Pacific, which is predicted to decrease in our model, due to a switch from Fe to NO₃ limitation, but increases in the simulation of Misumi et al. (2014), in which phytoplankton remain mostly Fe limited (see also Sect. 4.2 and Tagliabue et al., 2020).

Figure 6 Climate change effect on δ⁵⁶Fe_diss and its drivers. Present decade (2011–2020 average) and end-of-century (2091–2100 average) climate change simulation values of (a) δ⁵⁶Fe_diss (‰), (b) δ⁵⁶Fe_UF (‰), (c) δ⁵⁶Fe_EM (‰), and (d) δ⁵⁶Fe_UF (‰). The difference in future minus present values is shown in the bottom row.

Figure 7 Climate change effect on underlying mechanisms of δ⁵⁶Fe_diss drivers. Present decade (2011–2020 average) and end-of-century (2091–2100 average) climate change simulation values of (a) dFe concentration (µmol m⁻³), (b) Fe limitation, (c) primary production (molC m⁻³ yr⁻¹), and (d) nitrate (mmol m⁻³). The difference in future minus present values is shown in the bottom row.

The δ⁵⁶Fe_diss decrease in some higher-latitude areas (around ca. 60^∘ S) is driven by a higher concentration of dFe, likely due to redistribution of “new” dFe from external sources. This is consistent with the slight decrease in δ⁵⁶Fe_EM and the lower δ⁵⁶Fe_CF, which is generally characteristic of “younger” water (Fig. 6c, d; König et al., 2021). Such increases in dFe concentrations also lower the δ⁵⁶Fe_UF effect (Fig. 6b), as the Fe uptake to dFe concentration ratio has decreased.

4 Wider implications

As discussed above, Fe isotopes have been used in the past to track both external Fe sources and internal cycling. Here, we demonstrated how both historical and future surface ocean δ⁵⁶Fe_diss variability arises from a combination of the redistribution of “new” Fe from external sources and changes to internal cycling, especially Fe uptake by phytoplankton under Fe limiting conditions. This raises the question as to whether observed variations in δ⁵⁶Fe_diss can be used to infer alterations in external Fe distribution and Fe uptake and limitation.

4.1 Using δ⁵⁶Fe_diss to track changes in input of new Fe from external sources

Variability in surface ocean δ⁵⁶Fe_diss in the tropical Pacific is a great example of how δ⁵⁶Fe_diss may track new Fe input from external sources. Here, δ⁵⁶Fe_EM is likely responsible for a substantial fraction of δ⁵⁶Fe_diss variability (Figs. 1f, 2d) and tracks ENSO-related circulation changes that redistribute Fe zonally from different external sources (Sect. 3.2). Specifically, the eastward transport of isotopically light Fe from PNG sediments during El Niño events and the westward transport of light Fe from Peru sediments are discernible as light δ⁵⁶Fe_EM (and δ⁵⁶Fe_diss) anomalies in the respective regions (Fig. 3c, d). Simultaneous eastward and westward shifts in the open ocean region that receives most dust Fe input are also evident based on their more crustal δ⁵⁶Fe_EM (Fig. 3d).

Tracking external Fe input using δ⁵⁶Fe_diss is sensitive to the choice of δ⁵⁶Fe endmembers for the different (potential) new Fe sources. This has been discussed previously in a range of studies (e.g. Conway and John, 2014; König et al., 2022; Pinedo-González et al., 2020). For instance, in the tropical Pacific our model includes major new Fe inputs from continental margins, with this sedimentary-sourced Fe being isotopically light, with a δ⁵⁶Fe endmember of −1 ‰ in the upper ca. 400 m (as it is assumed to be released by reductive dissolution at these depths independent of location; König et al., 2021). Input of such light δ⁵⁶Fe_diss appears to be realistic for the eastern margin, as light δ⁵⁶Fe_diss have been observed in the Peru upwelling region (Chever et al., 2015; Fitzsimmons et al., 2016; John et al., 2018). However, isotopically heavier dFe has been observed in the western and central tropical Pacific and has been suggested to be released from PNG sediments via non-reductive processes with crustal (ca. 0.1 ‰) δ⁵⁶Fe endmembers (Labatut et al., 2014; Radic et al., 2011). A sensitivity test with neutral (0 ‰) sediment δ⁵⁶Fe endmember shows that, without any light sedimentary Fe input, δ⁵⁶Fe_EM variability and, consequently, δ⁵⁶Fe_diss variability are reduced in the surface ocean (Fig. S9). This consequently increases the dominance of fractionation effects (mainly δ⁵⁶Fe_UF; Fig. S9). Thus, in areas such as the western tropical Pacific and areas further east which receive Fe from the EUC (Sect. 3.2.2), δ⁵⁶Fe_EM variability may be reduced, as the modelled sediment Fe input from the PNG margin is too light. A greater extent of non-reductive Fe input would mean that changes in new Fe input or redistribution by ocean circulation would be harder to detect using δ⁵⁶Fe_diss. This is especially true in areas where fractionation effects are strong and thereby conceal changes in δ⁵⁶Fe_EM. On the other hand, δ⁵⁶Fe_diss could be useful to track new Fe input and redistribution in areas such as the eastern tropical Pacific, where the δ⁵⁶Fe endmember of new Fe is observed to be light, and where δ⁵⁶Fe_EM is responsible for a substantial part of δ⁵⁶Fe_diss variability.

Another complicating factor for tracking Fe input changes with δ⁵⁶Fe_diss is when fractionation and endmember effects overlap each other. Most fractionating processes lead to heavier δ⁵⁶Fe_diss in our model (except for colloidal pumping; König et al., 2021), whereas endmember effects are comparably light, even for dust Fe or Fe from non-reductive sediments with crustal δ⁵⁶Fe. Thus, if fractionation effects (δ⁵⁶Fe_UF and/or δ⁵⁶Fe_CF) do not vary in unison with δ⁵⁶Fe_EM, changes in Fe input may be hidden (as illustrated for the North Pacific in Fig. S3). However, where variability in fractionation and endmember effects occur simultaneously, fractionation effects may reinforce δ⁵⁶Fe_EM variability (Fig. S3). This is mostly the case in the tropical Pacific, especially at the Equator (Fig. S6), with some discrepancies in the central Pacific south of the Equator (Fig. S3). Here, changes in the redistribution of external Fe could even be deducted from δ⁵⁶Fe_UF alone.

In summary, it is necessary to have a good grasp on the relevant δ⁵⁶Fe endmembers when inferring changes in external Fe input or redistribution from δ⁵⁶Fe_diss, as well as understanding variability in δ⁵⁶Fe_UF and δ⁵⁶Fe_CF relative to δ⁵⁶Fe_EM (Fig. S3). Such changes are therefore easiest to spot where the external Fe source in question is known to provide isotopically light Fe and where fractionation effects are weak and/or invariable in time.

4.2 Using δ⁵⁶Fe_diss to track changes in Fe limitation and recycled Fe

Areas of strong Fe limitation are generally associated with heavy δ⁵⁶Fe_diss driven by strong uptake fractionation effects (i.e. heavy δ⁵⁶Fe_UF), often reinforced by relatively heavy δ⁵⁶Fe_EM due to limited input of light Fe from external sources (Sect. 3.2, 3.3). This suggests that heavy δ⁵⁶Fe_diss could be a potentially useful indicator of Fe limitation, which opens up the possibility of detecting or even monitoring changes in Fe limitation using δ⁵⁶Fe_diss. Such information would be valuable, since Fe is thought to limit primary production in large parts of the global ocean (Moore et al., 2013). However, the specific extent of Fe limitation is unknown and may vary over seasonal, interannual, and decadal timescales. Poorly constrained climate change-driven changes in Fe limitation have been demonstrated to explain large parts of the uncertainty in model projections of net primary production (NPP) in the tropical Pacific, where future NPP is predicted to depend decisively on whether or not phytoplankton switch from Fe to NO₃ limitation (Tagliabue et al., 2020; see also Sect. 3.3).

Outputs from our hindcast experiments show that the connection between heavy surface ocean δ⁵⁶Fe_diss and Fe limitation generally holds and that seasonal changes in Fe limitation are tracked by δ⁵⁶Fe_UF and, consequently, δ⁵⁶Fe_diss (Fig. S10). As discussed for the tropical Pacific (Sect. 3.2), this is due to the fact that uptake fractionation is strongest (i.e. δ⁵⁶Fe_UF heaviest) where the ratio of Fe uptake to dFe concentration is highest; this ratio also determines the degree of Fe limitation (where Fe is the limiting nutrient; Figs. 3, S5). However, in some regions, namely in the south and south-eastern Pacific, δ⁵⁶Fe_UF (and δ⁵⁶Fe_diss) can be heavy despite a relatively moderate ratio of Fe uptake to dFe concentration (Fig. S10), and δ⁵⁶Fe_UF variability appears to be anti-correlated to Fe limitation patterns (Fig. S3). While we cannot rule out that the heavy δ⁵⁶Fe_UF is caused by model artefacts, it may also be related to the ocean circulation pattern of these areas. In contrast to the equatorial Pacific, where ENSO dynamics supply water with light δ⁵⁶Fe_diss every few years, much of the water in the southeast Pacific originates in the seasonally Fe-limited Southern Ocean and contains Fe with relatively heavy δ⁵⁶Fe_diss (Fig. S10). Here, the heavy δ⁵⁶Fe_diss therefore might not (only) be caused by local processes but could be mostly a “legacy” signature. It may also indicate that this dFe has been continuously recycled by phytoplankton on its way to the south-east Pacific. Such recycling would leave dFe isotopically heavy given the preferential uptake of light Fe in our model (Sect. 2), as parts of the (isotopically lighter) phytoplankton Fe are removed from the surface ocean via particle settling. When deriving the local degree of Fe limitation from δ⁵⁶Fe_diss, it is therefore important to consider the “background” δ⁵⁶Fe_diss based on the origin of dFe and previous processing. Moreover, local δ⁵⁶Fe_EM effects should also be considered, especially if they are variable in time and space.

As discussed above, the connection between heavy δ⁵⁶Fe_diss and Fe limitation is a consequence of the uptake fractionation parametrization of our model. For this parametrization we use a constant fractionation factor (α = 0.9995) for Fe uptake by both phytoplankton classes, based on observations from an Fe-limited Southern Ocean eddy (Ellwood et al., 2020). Thus, the uptake fractionation strength is independent of factors such as dFe concentration, temperature, or species composition. This is in contrast to uptake of other nutrients such as ammonium, for which fractionation was found to decrease when concentrations were low (Sigman and Fripiat, 2019). While it is conceivable that δ⁵⁶Fe fractionation during phytoplankton uptake is impacted by dFe availability and/or other parameters, observational studies on the topic are limited. Exceedingly heavy surface ocean δ⁵⁶Fe_diss values have been observed in multiple, likely Fe-limited, systems (e.g. in the eastern tropical Pacific, Southern Ocean, South Atlantic; Fig. S10); for some of these systems, the heavy δ⁵⁶Fe_diss has indeed been attributed to continuous biological processing (i.e. uptake and recycling) combined with fractionation during organic complexation (Ellwood et al., 2020; Sieber et al., 2021). This supports the concentration-independent uptake fractionation of our model since heavy surface ocean δ⁵⁶Fe_diss would not emerge in model simulations if uptake fractionation was set to decrease for low dFe. Nevertheless, more observational or, ideally, experimental data are necessary to determine which factors may impact uptake fractionation and to find the appropriate fractionation factor parametrizations.

Overall, heavy δ⁵⁶Fe_diss could be a useful indicator of Fe limitation if complicating factors such as past processing and variable δ⁵⁶Fe_EM are taken into account. Fe limitation derived from δ⁵⁶Fe_diss observations could thereby complement other measures of nutrient limitation, such as the (scarce) limitation data from incubation experiments (Moore et al., 2013) or the limitation inferred from nutrient deficiencies (Moore, 2016). Our modelling results suggest that changes in Fe limitation can induce strong seasonal variability in δ⁵⁶Fe_diss in some locations (this study, König et al., 2022) but also variability over interannual and decadal scales. Thus, it would be worthwhile to observe δ⁵⁶Fe_diss changes over seasonal and interannual timescales in the form of a δ⁵⁶Fe_diss time series, ideally in a place with variable (degrees of) Fe limitation. A possible candidate location could be the Southern Ocean Time Series (140^∘ E, 47^∘ S), which is Fe limited in summer (Boyd et al., 2001; Sedwick et al., 1999). Here, our model predicts substantial seasonal variability in surface ocean δ⁵⁶Fe_diss (between ca. 0.1 ‰ and 0.7 ‰); available observations are within this range (Barrett et al., 2021; Ellwood et al., 2020). Moreover, based on predicted present and future variability in Fe limitation and, consequently, δ⁵⁶Fe_diss in our model (Figs. 1–2), the equatorial and south-eastern Pacific could be potential locations for future studies that explore the interannual changes in iron limitation.

5 Conclusion

Simulations of a global ocean model with an active δ⁵⁶Fe cycle and variable climate forcings show that surface ocean δ⁵⁶Fe_diss responds distinctly to the changes in ocean physics and biogeochemistry triggered by natural climate variability and long-term global warming effects. Changes in δ⁵⁶Fe_diss thereby integrate both variations in external Fe supply and redistribution and shifts in upper ocean biogeochemistry (especially presence and degree of Fe limitation), which alter δ⁵⁶Fe endmember and uptake fractionation effects respectively. We therefore suggest that regular observations of surface ocean δ⁵⁶Fe_diss as part of a long-term time series could be useful to track climate-driven changes to upper ocean Fe supply or the degree of Fe limitation in phytoplankton. As endmember and fractionation effects can overlap, δ⁵⁶Fe_diss and dFe concentration data should ideally be accompanied by ancillary measurements to constrain changes in upper ocean circulation or primary productivity. In areas with high seasonal variability in δ⁵⁶Fe_diss (e.g. due to strong seasonality in primary productivity), a sub-annual sampling interval should be considered.