Articles | Volume 19, issue 7
Research article
28 Jul 2023
Research article |  | 28 Jul 2023

Impact of iron fertilisation on atmospheric CO2 during the last glaciation

Himadri Saini, Katrin J. Meissner, Laurie Menviel, and Karin Kvale

While several processes have been identified to explain the decrease in atmospheric CO2 during glaciations, a better quantification of the contribution of each of these processes is needed. For example, enhanced aeolian iron input into the ocean during glacial times has been suggested to drive a 5 to 28 ppm atmospheric CO2 decrease. Here, we constrain this contribution by performing a set of sensitivity experiments with different aeolian iron input patterns and iron solubility factors under boundary conditions corresponding to 70 000 years before present (70 ka), a time period characterised by the first observed peak in glacial dust flux. We show that the decrease in CO2 as a function of Southern Ocean iron input follows an exponential decay relationship. This exponential decay response arises due to the saturation of the biological pump efficiency and levels out at ∼21 ppm in our simulations. We show that the changes in atmospheric CO2 are more sensitive to the solubility of iron in the ocean than the regional distribution of the iron fluxes. If surface water iron solubility is considered constant through time, we find a CO2 drawdown of ∼4 to ∼8 ppm. However, there is evidence that iron solubility was higher during glacial times. A best estimate of solubility changing from 1 % during interglacials to 3 % to 5 % under glacial conditions yields a ∼9 to 11 ppm CO2 decrease at 70 ka, while a plausible range of CO2 drawdown between 4 to 16 ppm is obtained using the wider but possible range of 1 % to 10 %. This would account for ∼12 %–50 % of the reconstructed decrease in atmospheric CO2 (∼32 ppm) between 71 and 64 ka. We further find that in our simulations the decrease in atmospheric CO2 concentration is solely driven by iron fluxes south of the Antarctic polar front, while iron fertilisation elsewhere plays a negligible role.

1 Introduction

CO2 drawdown during the Last Glacial Period occurred in multiple steps before reaching a minimum level of ∼190 ppm at the Last Glacial Maximum (LGM, ∼21 000 years ago) (EPICA Community Members et al.2004; Ahn and Brook2008; Lüthi et al.2008; Bereiter et al.2012, 2015). One of the last major drops in atmospheric CO2 concentration occurred at the beginning of Marine Isotope Stage (MIS) 4 (71 000–59 000 years ago, ka hereafter) when CO2 decreased from 233±8 ppm to 201±4 between 71 and 64 ka (Fig. 1b). This period also coincided with the minimum summer insolation at high northern latitudes (∼70 ka, Berger1978), which led to the most pronounced episode of ice sheet growth in the Northern Hemisphere (NH) over the glaciation (Bassinot et al.1994; Petit et al.1999; Grant et al.2012). Even though NH ice sheet extent and volume during MIS4 and MIS3 are not well constrained, numerical modelling experiments and observational records have shown that the extent of the Cordilleran and Scandinavian ice sheets was larger during MIS4 than during the LGM, whereas the Laurentide and Greenland ice sheets were smaller (Lambeck et al.2010; Kleman et al.2013; Batchelor et al.2019). As a result of the glaciation occurring during the early part of MIS4, it is suggested that the global sea level dropped to about 80 m below the present-day sea level (Waelbroeck et al.2002; Grant et al.2012; Batchelor et al.2019; De Deckker et al.2019). Global sea surface temperatures (SSTs) also show a rapid decline of ∼1C in the early MIS4 transition and an overall ∼1.5 drop across MIS4, before rising again to pre-MIS4 levels around 59 ka (Kohfeld and Chase2017).

Figure 1Time series of (a) Antarctic temperature anomalies from present day (C) (Jouzel et al.2007), (b) atmospheric CO2 concentration (ppm) (Bereiter et al.2015), and (c) dust flux (mg m−2 yr−1) (Lambert et al.2012) as recorded in EPICA DOME C ice core. (d) Iron dust accumulation rates (mg m−2 yr−1) from ODP Site 177-1090 (Atlantic) (Martínez-García et al.2011). Iron (%) records in the southern Pacific from Lamy et al. (2014) at sites (e) PS75-076 and (f) PS75-059. (g) Atmospheric δ13CO2 (‰) as recorded in a composite of Antarctic ice cores (purple line, Eggleston et al.2016) and high-resolution records from Taylor Glacier ice cores (dashed cyan line, Menking et al.2022). Shaded areas represent marine isotope stages 2 (MIS2: 27–19 ka) and 4 (MIS4: 71–59 ka).


Several mechanisms have been put forward to explain the drawdown of atmospheric CO2 during glaciations, including the ∼32 ppm drop in CO2 during the MIS4 transition. These mechanisms include higher solubility of CO2 in colder ocean waters (Heinze et al.1991; Kucera et al.2005; Williams and Follows2011; Khatiwala et al.2019), higher carbon sequestration associated with a weaker Atlantic meridional overturning circulation (AMOC) and thus lower ventilation rates during glacial periods (Sigman and Boyle2000; Toggweiler2008; Sigman et al.2010; Watson et al.2015; Jaccard et al.2016; Yu et al.2016; Menviel et al.2017; Kohfeld and Chase2017), and the expansion of sea ice cover leading to more stratified Southern Ocean (SO) waters and smaller air–sea gas exchange (Francois et al.1997; Stephens and Keeling2000; Ferrari et al.2014). However, it has been argued that sea ice expansion does not correlate well with the timing of CO2 drawdown at 70 ka (Kohfeld and Chase2017) and that the major drivers during this transition might in fact be due to a shallower AMOC (Piotrowski et al.2005; Thornalley et al.2013; Yu et al.2016) or a more efficient biological pump due to enhanced iron input to the ocean (Brovkin et al.2012; Menviel et al.2012; Anderson et al.2014; Lamy et al.2014; Martínez-García et al.2014). Interestingly, high-resolution δ13CO2 records from Antarctic ice cores (Menking et al.2022) display a 0.5 ‰ decrease centred at 70.5 ka, followed by a 0.7 ‰ increase (Fig. 1g), indicating a complex set of processes impacting atmospheric CO2 at the MIS5–4 transition. While surface ocean cooling could explain the concurrent CO2 and δ13CO2 decrease, the δ13CO2 increase in the second part of the transition would be consistent with a greater efficiency of the biological pump and increased storage of respired carbon in the deep ocean (Menviel et al.2015; Eggleston et al.2016; Menking et al.2022).

The efficiency of the biological pump is partly dependent on the relative abundance of different marine phytoplankton communities, which further depends on the availability of both macro- and micro-nutrients (Kvale et al.2015a; Saini et al.2021). Martin (1990) suggested that the micro-nutrient iron plays a crucial role as a limiting nutrient in phytoplankton growth, and thus put forward the hypothesis that iron fertilisation, resulting from an increase in atmospheric dust during glacial times, might have increased marine net primary production (NPP) in the Southern Ocean, leading to a decrease in atmospheric CO2 concentration. Antarctic ice core records indeed show peaks in dust fluxes that coincide with lower Antarctic temperatures and lower CO2 levels (Fig. 1) during MIS4 and MIS2 (27–19 ka) (Wolff et al.2006, 2010; Lambert et al.2008, 2012; Martínez-García et al.2011, 2014; Lamy et al.2014). These peaks in dust flux are concurrent with increased export production (EP) in the sub-Antarctic zone (SAZ) of the Southern Ocean during the LGM (Kohfeld et al.2005, 2013; Martínez-García et al.2014), as well as during MIS4 (Lamy et al.2014; Martínez-García et al.2014; Thöle et al.2019; Amsler et al.2022). On the other hand, palaeoceanographic data from the Antarctic zone (AZ) suggest a decrease in EP at the LGM and MIS4 (Kohfeld et al.2005, 2013; Anderson et al.2009; Jaccard et al.2013; Studer et al.2015; Thöle et al.2019; Amsler et al.2022; Weber et al.2022).

Numerous modelling studies have explored the changes in Southern Ocean biogeochemistry due to enhanced aeolian iron input either under pre-industrial (PI) (Aumont and Bopp2006; Tagliabue et al.2009a, 2014) or LGM conditions (Aumont et al.2003; Bopp et al.2003; Oka et al.2011; Lambert et al.2015; Muglia et al.2017; Khatiwala et al.2019; Yamamoto et al.2019; Lambert et al.2021; Saini et al.2021) but not under MIS4 conditions. An exception is Menviel et al. (2012), who simulated a ∼12 ppm CO2 drop as a result of enhanced aeolian iron fluxes into the Southern Ocean at 70 ka in a transient simulation of the last glacial cycle.

In this study, we constrain the impact of enhanced aeolian iron input on atmospheric CO2 concentration during the MIS4 transition. We perform a set of sensitivity experiments with different aeolian iron flux masks and different iron solubilities under 70 ka boundary conditions using a model of intermediate complexity that includes an ocean general circulation model (OGCM) coupled to a recently developed complex ecosystem model (Kvale et al.2021; Saini et al.2021). The ecosystem model used here includes calcifying and silicifying plankton and an iron cycle, which makes it a well-suited model for studying CO2 uptake in the Southern Ocean.

2 Methods

2.1 Model description

This study uses version 2.9 of the University of Victoria Earth System Climate Model (UVic ESCM), which consists of a sea ice model (Semtner1976; Hibler1979; Hunke and Dukowicz1997) coupled to the ocean general circulation model MOM2 (Pacanowski1995), with 19 ocean depth layers and a spatial resolution of 3.6 by 1.8. It also includes a dynamic vegetation model (Meissner et al.2003), a land surface scheme (Meissner et al.2003), a 2-D atmospheric energy moisture balance model integrated vertically (Fanning and Weaver1996), and a sediment model (Archer and Maier-Reimer1994; Meissner et al.2012). The model is forced with seasonally varying wind stress and wind fields (Kalnay et al.1996) and seasonal variations in solar insolation at the top of the atmosphere. Full physical and structural descriptions of the model can be found in Weaver et al. (2001), Meissner et al. (2003), Eby et al. (2009), and Mengis et al. (2020).

The marine carbon cycle is represented by the newly developed Kiel Marine Biogeochemistry Model version 3 (KMBM3) (Kvale et al.2021), which is based on the Nutrient Phytoplankton Detritus Zooplankton model of Schmittner et al. (2005) and Keller et al. (2012). There are four different classes for phytoplankton in this ecosystem model. Three of which include specifically characterised plankton such as diazotrophs that can fix nitrogen, coccolithophores that produce CaCO3 shells, and diatoms that produce opal. The fourth class is for the rest of the types of plankton, which are mostly located in the low-latitude regions. The model also includes prognostic CaCO3 and silica tracers (Kvale et al.2015b, a, 2021) and dissolved nitrate, phosphate, iron, and silica as nutrients. The model also incorporates an iron cycle (Nickelsen et al.2015), including hydrothermal sources. The ecosystem model is described in detail in Kvale et al. (2021).

2.2 Experimental design

A pre-industrial control simulation (PI-control) is integrated under PI boundary conditions, including an atmospheric CO2 concentration of 283.86 ppm and PI iron (pife) and silica (pisi) fluxes. The PI aeolian dust fluxes are based on the PI-BASE dust flux simulation of Mahowald et al. (2006). To obtain iron and silica fluxes, we multiply the dust flux with the percentage distribution of iron and silica in dust (Zhang et al.2015), respectively. The solubility factor for iron in the PI-control simulation is set to 1 %.

All the other experiments are run under 70 ka boundary conditions, including orbital parameters corresponding to the year 70 ka according to Berger (1978), and a continental ice sheet topography generated by an offline ice sheet simulation performed with IcIES (Ice sheet model for Integrated Earth system Studies, Abe-Ouchi et al.2013). The global mean ocean alkalinity was adjusted on the basis of sea level differences. For the control simulation (70ka-control), the model is run into equilibrium with an atmospheric CO2 concentration of 222.5 ppm (Bereiter et al.2015), after which the simulation is integrated for an additional 200 years with prognostic atmospheric CO2; 70ka-control is forced with pife (with 1 % iron solubility) and pisi fluxes.

From this control simulation, we branch off a suite of sensitivity experiments with prognostic atmospheric CO2, using two different glacial iron dust flux estimates. The first iron flux, lambfe, is obtained from the LGM dust estimates of Lambert et al. (2015). This dust flux estimate is calculated by performing a global 2-D interpolation on unevenly distributed LGM dust flux records, most of which are collated in the DIRTMAP (Dust Indicators and Records of Terrestrial and Marine Paleoenvironments) database (Kohfeld and Harrison2001; Maher et al.2010). Their interpolation method assumes that the aerosol concentration in the air decreases exponentially from the source distance. The second iron flux, glacfe, is derived from a dust flux obtained with a model (LGMglac.a, Ohgaito et al.2018). This model includes glaciogenic dust sources in addition to the usual desert dust sources and assumes dry and unvegetated regions as dust sources. It then calculates global dust transport and deposition. Both glacfe and lambfe are then obtained in this study by mapping the iron percentage on dust (Zhang et al.2015) as mentioned above. The aeolian iron fluxes based on glacfe and lambfe dust patterns compared to the PI iron flux (pife) and assuming a 1% solubility are shown in Fig. 2, while the full dust deposition maps for the PI period and the two LGM reconstructions are available in the Supplement of Saini et al. (2021).

Figure 2Aeolian iron dust flux (µmol m−2 yr−1) anomalies for (a) lambfe1% minus pife1% and (b) glacfe1% minus pife1%.

The iron delivered to the ocean surface via aeolian dust deposition needs to be dissolved before it can become available to phytoplankton. The solubility factor of iron is therefore an important parameter that will impact the magnitude of iron fertilisation. Iron solubility varies depending on both provenance and speciation. For example, Baker et al. (2006) find that solubilities in samples of Saharan dust (median of 1.7 %) are significantly lower than solubilities in aerosols from other source regions feeding into the Atlantic Ocean (median of 5.2 %). Schroth et al. (2009) find a significant difference in iron solubility between arid soils (1 %), glacial flour (2 %–3 %) and, of lesser importance to our study, oil combustion products (77 %–81 %). In the present-day Southern Ocean, observations of iron solubilities vary between 0.2 % and 48 % (Ito et al.2019), with the higher values being hypothesised to be caused by pyrogenic iron. The majority of the observational data in the present-day Southern Ocean lie between 1 % and 12 % (Ito et al.2019). Quantifying iron solubilities during past climates is even more challenging. Conway et al. (2015) estimate iron solubility at the Last Glacial Maximum (LGM) based on dust particles in the EPICA Dome C and Berkner Island ice cores. They find that iron solubility was very variable during the LGM interval at Dome C (1 %–42 %), with a mean of 10 % and a median of 6 %. On the other hand, iron solubility at Berkner Island was 1 %–5 % at the LGM, and on average ∼3 % between 23 and 50 ka. The lower iron solubility at Berkner Island is thought to be due to the aerosol composition. Due to its proximity to South America, the solubility found at Berkner island might better represent large-scale aeolian dust deposited in the Southern Ocean. Iron is more bioavailable in dust that originates from physically weathered bedrock than from chemically weathered bedrock (Shoenfelt et al.2017). The analysis of sub-Antarctic marine sediment cores further suggests that aeolian iron was 15 to 20 times more bioavailable during glacial periods than during the current interglacial (Shoenfelt et al.2018).

In this study, we test a global mean solubility range of 1 % to 20 % (Table 1). Based on the studies mentioned above, we define the plausible range of average iron solubility to be 1 % to 10 %, with an increased likelihood at 3 %–5 % during colder climate episodes. During warmer climates, such as the PI period or at the MIS4 onset, we define the most likely solubility factor as 1 %.

Table 1Aeolian iron fluxes for each experiment.

Download Print Version | Download XLSX

However, due to the uncertainties associated with present-day iron solubilities, we perform additional sensitivity experiments under 70 ka boundary conditions using the PI iron dust mask with iron solubility varying between 3 % and 20 %. This approach allows us to estimate the minimum change in CO2 due to glacial dust fluxes, assuming no change in solubility over time. The corresponding CO2 changes can be calculated by taking the difference between CO2 changes achieved with the full experiments (i.e. changing masks and solubilities) and the CO2 changes achieved by only changing solubility. This approach was validated by performing two additional 70 ka equilibrium experiments with the pife mask and an iron solubility of 3 % and 10 % from which we branched off simulations with the lambfe mask and constant solubility (not shown). The resulting CO2 drawdown in these experiments was the same as if it was calculated as the difference between the full experiments and solubility-only experiments. The integrated aeolian iron input for all experiments is listed in Table 1.

We perform five additional sensitivity experiments (Fig. A1) to better quantify the contribution of the Southern Ocean to the total CO2 drawdown. In these experiments, the aeolian iron flux in the Southern Ocean follows the lambfe mask while the PI aeolian iron fluxes are applied outside of the Southern Ocean. In four of these experiments (lambfe10%-30S, lambfe10%-40S, lambfe10%-50S, lambfe10%-60S), the lambfe mask with 10 % solubility (Fig. A2) is applied south of 30, 40, 50, and 60 S, respectively. In the fifth sensitivity experiment (lambfe50%-47S), the aeolian iron input south of 47 S follows the lambfe mask with a solubility factor of 50 % (Fig. A1f), which is equivalent to 23.05 Gmol yr−1 iron input in the Southern Ocean south of 47 S and provides an upper limit on the potential CO2 drawdown.

2.3 Carbon decomposition

To better understand the changes in ocean carbon, we decompose the simulated dissolved inorganic carbon (DIC) into its three major components: respired organic carbon (Creg), DIC generated by dissolution of calcium carbonate (CCaCO3), and preformed carbon (Cpref) as described below.

Figure 3(a) Annual mean SST anomalies (C, shading) at 70 ka compared to the PI period and 15 % annual mean sea-ice concentration for the PI period (green line) and 70 ka (red line). The grey shading on land shows the extent of continental ice sheets at 70 ka (Abe-Ouchi et al.2013). (b) Export production anomalies between 70 ka and the PI period at 177.5 m depth (g C m−2 yr−1). Overlaid dashed contours represent the modern sub-Antarctic front (SAF) in black and Antarctic polar front (APF) in green based on Sokolov and Rintoul (2009). The global overturning streamfunction (Sv) at (c) 70 ka and (d) PI.

Creg is calculated based on the remineralised phosphate in the ocean (Preg) and the carbon to phosphate stoichiometric ratio (RC/P=106):

(1) Δ C reg = Δ P reg × R C / P .

Preg is determined based on apparent oxygen utilisation (AOU) and the phosphate-to-oxygen stoichiometric ratio (RP/O2=1/160), ΔPreg=ΔAOU×RP/O2, where AOU is calculated as the difference between the saturated oxygen concentration (O2sat) and the dissolved oxygen concentration in the ocean (O2): AOU=O2sat-O2.

CCaCO3 is calculated as follows:

(2) Δ C CaCO 3 = 0.5 Δ ALK + R N / P × Δ P tot ,

where RN/P=16. Finally, Cpref is calculated as the remainder as follows:

(3) Δ C pref = Δ DIC - Δ C reg - Δ C CaCO 3 .

To further assess the efficiency of the biological pump, we calculate global P* (Ito and Follows2005), defined as follows:

(4) P * = P reg / P total ,

where Ptotal is the total phosphate content in the ocean.

3 Results

3.1 Simulated ocean conditions at 70 ka

In this section, we describe the simulated physical and biological conditions at 70 ka, at the onset of MIS4 (Fig. 3). In our 70ka-control simulation, the globally averaged ocean temperature is ∼2.8C, about 0.6 C lower than in PI-control but 0.9 C higher than in a LGM simulation integrated with the same model (Saini et al.2021). The globally averaged annual mean SST is 17.1 C, 0.8 lower than in the PI-control and 1.1 higher than in the LGM simulation. The strongest ocean surface cooling (−1.45C) at 70 ka with respect to our PI-control is simulated north of 40 N in the Atlantic and Pacific oceans, while the annual mean SSTs in the SAZ and AZ are ∼0.8 and 0.4 C lower than in the PI-control, respectively (Fig. 3a). At 70 ka, the annual mean Southern Ocean sea ice edge is situated at 55 S in the southern Atlantic Ocean and the southern Indian Ocean, ∼1 further north than in the PI-control simulation, while the change is insignificant in the southern Pacific sector (sea ice edge at ∼62 S). The simulated AMOC strength in the 70ka-control experiment is ∼14 Sv, compared to ∼17 Sv in the PI-control, without significant changes in its depth (Fig. A3). Furthermore, there are no significant changes in simulated Antarctic Bottom Water (AABW) formation (Fig. 3c and d).

The colder conditions, more extensive sea-ice cover, and changes in the deep-ocean convection impact marine productivity (Saini et al.2021). Globally, NPP decreases by ∼9.7 %, while EP decreases by ∼3 % in our 70ka-control simulation compared to PI-control (both integrated with the pife mask and 1 % iron solubility, Table A1). The simulated 18 % EP increase at 70 ka in the northern Pacific (Fig. 3b) can be attributed to a greater diatom and coccolithophore abundance in that region (Fig. A4a and b) resulting from higher nutrient availability. This EP increase in the northern Pacific is, however, inconsistent with the biogenic Ba (Jaccard et al.2005) and δ15N (Gebhardt et al.2008) records from the sub-arctic northern Pacific that suggest a decrease in EP at MIS4 onset. The northern Atlantic shows a complex pattern of anomalies due to changes in the strength and location of deep-ocean convection, which result in an overall decrease in NPP and EP by 16 % and 9 %, respectively.

Within the AZ (south of the APF), diatoms decrease by 14 % in the Pacific, while they increase by 3 % in the Atlantic and Indian sectors (Fig. A4a). The decrease in diatoms in the Pacific sector of the AZ leads to a competitive growth advantage for coccolithophores, which increase by 6.5 % south of the polar front, thus leading to a poleward shift in coccolithophores in the Pacific sector (Fig. A4b). On the other hand, coccolithophores decrease in the Atlantic and Indian sectors of the AZ by 15 % and 8 %, respectively. Diatoms increase by 22 % in the SAZ (north of the SAF) while there are no significant changes in coccolithophores abundance. As a result of the changes in these two plankton species, EP increases by 1.3 % in the SAZ and decreases by 14 % in the AZ (Fig. 3b). The total Southern Ocean (south of 30 S) EP and NPP decrease by 2 % and 7.5 %, respectively, at 70 ka.

3.2 Impact of changes in iron dust flux on atmospheric CO2

Changing the iron flux masks from PI to glacial conditions at 70 ka leads to a 3.8 to 8.3 ppm drop in atmospheric CO2 concentration if we assume that the mean iron solubility remains unchanged (Table 2). Interestingly, for solubilities of 3 % and higher, the drawdown is nearly constant, regardless of the glacial dust flux mask and regardless of the solubility (7.3±1 ppm).

Table 2ΔCO2 simulated by changing the iron mask from PI to glacial data (lambfe and glacfe) at 70 ka but keeping the solubility constant.

Download Print Version | Download XLSX

However, the solubility of iron was likely higher during cold than during warm periods (see Sect. 2.2). We will therefore discuss experiments that switch from a PI iron mask with 1 % solubility to glacial iron masks with higher solubilities from hereon. These experiments show a ∼9 to 19 ppm drop in atmospheric CO2 concentration (Fig. 4a and b, Table 3). These changes in atmospheric CO2 impact the physical conditions of the ocean only slightly in our model. As such, the global overturning circulation, mean SST, and sea ice extent are similar in most of the sensitivity experiments described below.

Figure 4(a) Equilibrated atmospheric CO2 concentration (ppm) as a function of the globally integrated aeolian iron flux into the ocean (Gmol yr−1). (b) Atmospheric ΔCO2 concentration (ppm) as a function of changes in aeolian iron flux into the Southern Ocean south of 47 S. The grey shading represents the range of likely glacial iron solubilities (1 %–10 %) and the associated change in CO2 concentration (−4 to −16 ppm), while the orange shading represents our best estimate of change in CO2 (−9 to −11 ppm) for a solubility of 1 % during warm periods and 3 %–5 % during colder periods. Note that lambfe10%, lambfe10%-30S, and lambfe10%-40S are overlapping. The black curve represents the best fit and suggests a maximum CO2 drawdown of ∼21 ppm due to Southern Ocean iron fertilisation. (c) Global P* as a function of aeolian iron flux into the Southern Ocean south of 47 S.(d) Globally integrated carbon reservoirs (GtC) as a function of aeolian iron flux into the Southern Ocean south of 47 S. The shading represents ocean carbon (blue), atmospheric carbon (red), and terrestrial carbon (green).


Table 3ΔCO2 for the sensitivity experiments compared to 70 ka-control, globally averaged P* values, and NPP and EP values integrated over different regions of the Southern Ocean. Percentage changes from the 70ka-control experiment are provided in brackets.

Download Print Version | Download XLSX

As the global and Southern Ocean iron input increases, the atmospheric CO2 concentration decreases, but it does not do so linearly, as the efficiency of the iron fertilisation weakens with the increasing iron availability (Fig. 4a and b). The patterns of aeolian iron deposition in the Southern Ocean are different for the two masks used here. Both reconstructions show an increase in the southern Atlantic sector compared to PI fluxes, but they differ over the Pacific and Indian sectors, including large regions where dust deposition decreases compared to PI conditions (Fig. 2). It is therefore interesting that our experiments with similar total iron input into the Southern Ocean but different regional patterns (e.g. glacfe20% and lambfe20%, Fig. 4b, pink star and black circle) display similar efficiency in drawing down CO2 regardless of the dust mask used. The glacfe mask is, however, slightly more efficient in drawing down atmospheric CO2. The iron flux into the southern Atlantic sector is higher in the glacfe mask compared to the lambfe mask, increasing iron concentrations near one of our major convection sites in the Weddell Sea. Bottom water formation in this region leads to greater mixing with deeper layers, replenishing surface nutrient concentrations and resulting in higher export production (Marinov et al.2006). Global P* is ∼1.5 % higher in our glacfe experiments than in our lambfe experiments, resulting in a slightly larger CO2 drawdown.

Figure 5Export production anomalies (g C m−2 yr−1) at 177.5 m depth for (a) lambfe3% and (b) glacfe3% compared to 70ka-control (pife1%). Overlaid dashed contours represents modern SAF in black and APF in green based on Sokolov and Rintoul (2009).

The sensitivity experiments provide additional information on the role of the Southern Ocean in our simulations. We find that experiment lambfe10% is as efficient in drawing down CO2 as the experiments where iron fluxes were only enhanced south of 30 or 40 S (lambfe10%-30S, lambfe10%-40S; Fig. 4a and b)). This indicates that in our experiments iron fertilisation in the Southern Ocean south of 40 S is mainly responsible for the atmospheric CO2 drawdown, while fertilisation elsewhere plays a negligible role. On the other hand, a 14 and 8.5 ppm CO2 decrease is simulated in lambfe10%-50S and lambfe10%-60S, respectively, compared to 15.3 ppm in lambfe10%.

Our experiments show that the ecosystem response north and south of ∼47 S (which is roughly the modern SAF) are of opposite sign (Table 3). While an increase in iron availability leads to an increase in NPP, diatom and coccolithophore abundance, and consequently EP south of the SAF (Figs. 5 and 6), the regions north of the SAF show a significant decline in all of these parameters. In our model, the EP increase occurs south of the SAF, where the surface nutrient concentrations are high (Fig. A5a). Enhanced aeolian iron input leads to a more efficient use of nutrients in that region and results in a reduction in northward transport of nutrients north of the polar front (Fig. A6), which in turn leads to a decrease in EP in the SAZ.

Figure 6Zonally integrated anomalies of (a) EP (g C m−1 yr−1) at 177.5 m, (b) NPP (g C m−1 yr−1), (c) diatoms (g C m−1) and (d) coccolithophores (g C m−1) compared to 70ka-control. Please note the multiplication factor allocated to each panel.


The relationship between atmospheric CO2 and the iron flux south of 47 S is as follows (Fig. 4b):

(5) Δ CO 2 = a 1 - b × e - c × ( SOfe ) d ,

where ΔCO2 is the atmospheric CO2 anomaly in parts per million compared to 70ka-control, SOfe denotes the total iron input into the Southern Ocean south of 47 S (in Gmol yr−1), a is −21.51 ppm, b is 1.0, c is 0.487 (Gmol yr−1)−1, and d is 0.66.

We also note that the CO2 drawdown capacity levels out for high iron forcing. Both lambfe and glacfe masks with 20% iron solubility lead to a ∼19 ppm CO2 decrease. In our extreme experiment lambfe50%-47S, with an iron input south of 47 S almost 2.5 times higher than in lambfe20% and glacfe20%, only an additional 2 ppm CO2 drawdown is simulated (Fig. 4b). Our results therefore suggest that iron fertilisation could have led to a maximum CO2 decrease of ∼19–21 ppm at 70 ka.

The enhanced iron flux leads to an increase in Preg in the Southern Ocean, which is subsequently transported northward at depth by AABW. As a result, the globally averaged P*, and therefore the overall efficiency of the soft tissue pump, increases with Southern Ocean aeolian iron input (Fig. 4c). As ocean circulation and ocean temperatures change only slightly in our experiments, the relationship between P* and the resulting CO2 drawdown is almost linear (Fig. A7) (Ito and Follows2005).

However, as the Southern Ocean iron flux increases, the overall iron fertilisation efficiency reduces. While export production increases south of 47 S, thus using nutrients more efficiently (Fig. A8c) and reducing the nutrient advection north of the SAF (Fig. A8b), nitrate limitation increases in the SAZ (Figs. A9b, d and A10a, b), leading to a decrease in export production. For example, in experiment lambfe20%, a 67 % decrease in nitrate concentration is simulated between 40 and 47 S (Fig. A8b). At the same time, silicate limitation decreases in the SAZ (Figs. A9e and A10c) due to a southward shift in both coccolithophores and diatoms and a decrease in diatoms between 40 and 55 S (Fig. A11). This further enhances silica availability in this region, consequently leading to a decrease in silicate limitation. Therefore, the total biological pump efficiency, represented here by changes in P* (Fig. 4c) and Southern Ocean EP to NPP ratio (Fig. A8a), saturate at high iron values due to nitrate limitation north of 50 S in the Southern Ocean. The global P* in lambfe20% and in lambfe50%-47S (Fig. 4c) equal 0.63 and 0.65, respectively, suggesting a maximum efficiency of 65 % in our experimental set-up.

The oceanic carbon reservoir in our sensitivity experiments increases between 20 and 130 Gt C depending on the iron flux scenario (Fig. 4d, blue shade). The simulated decrease in atmospheric CO2, equivalent to 8–45 Gt C (Fig. 4d, red shading), leads to a decrease in surface air temperatures, as well as regional changes in precipitation and soil moisture. In addition, the lower atmospheric CO2 concentration also reduces photosynthesis and consequently litter fall. The direct and indirect effects of a lower atmospheric CO2 concentration result in a terrestrial carbon decrease of 16 to 88 Gt C (Fig. 4d, green shade), out of which 8 to 45 Gt C is from terrestrial vegetation while 8 to 43 Gt C is from soil carbon.

3.3 Impact of aeolian iron input on the distribution of oceanic carbon

As mentioned in the previous section, our sensitivity experiments do not show significant changes in the physical ocean conditions. However, enhanced iron input significantly impacts marine ecosystems. To understand the resulting ocean carbon changes, we start by describing changes in the simulated ecosystems for one of our sensitivity experiments, lambfe3%, compared to 70ka-control. We choose lambfe3% (Fig. A2) because earlier research shows that the lambfe Southern Ocean dust mask provides a better fit with available glacial dust flux proxy records than the glacfe mask (Saini et al.2021). Furthermore, a solubility factor of 3 % corresponds to a likely estimate of iron solubility at 70 ka (see Sect. 2) and has also been used in previous studies (Tagliabue et al.2014; Lambert et al.2015; Muglia et al.2017; Yamamoto et al.2019). All lambfe experiments with different solubility factors show similar anomaly patterns, only the magnitude of the changes varies as a function of solubility (Fig. A11).

A ∼9 ppm drop in atmospheric CO2 is simulated in experiment lambfe3% with a global increase in EP by 1.3 % and a global decrease in NPP by 6.7 % (Tables 3 and A1). Iron increase in the lambfe3% experiment leads to a 37 % increase in diatoms and an 88 % increase in coccolithophores in the AZ (Fig. A4c and d). The simulated increased EP in the AZ (Fig. 5a) leads to greater nutrient utilisation south of the APF (Fig. A8c). In contrast, because of lower nutrient availability, diatom and coccolithophore abundances decrease in the SAZ by 46 % and 31 %, respectively (Fig. A4c and d). Consequently, in the lambfe3% experiment, the EP increases by 98 % in the AZ and decreases by 17 % in the SAZ compared to the 70ka-control simulation (Figs. 5a and  6a).

The overall Southern Ocean EP south of 30 S increases by ∼9 % in the lambfe3% experiment, the NPP changes become insignificant (Table 3), and the total biological pump efficiency increases by 6.7 % (Table 3). In the northern Pacific (47–56 N, 150–220 E), EP increases by 13 %, while NPP increases by 7 % due to an 82 % increase in coccolithophores and 58 % decrease in diatoms. Both NPP and EP decrease by ∼10 % in the northern Atlantic (47–56 N, 60 W–0 E), where we see an overall decrease in coccolithophores of 22 % and an increase in diatoms of 16 %.

The increased carbon export from the surface into the deep ocean increases DIC at all depths south of 47 S. The DIC-rich bottom water is advected northward and leads to a stronger vertical DIC gradient in the ocean (Fig. 7a and b). While DIC increases by 25 mmol m−3 in the deep Southern Ocean (south of 30 S, below 3000 m), it does not increase uniformly across sectors: it increases by 22 mmol m−3 in the deep Indo-Pacific sector (below 3000 m, Fig. 7a) and by 29 mmol m−3 in the deep Atlantic sector (below 3000 m, Fig. 7b). The positive DIC anomalies north of 30 S are associated with AABW and its northward spread into the Atlantic and the Indo-Pacific basins, leading to a 10 mmol m−3 DIC increase in the deep Indo-Pacific sector and a 15 mmol m−3 increase in the deep Atlantic sector.

Figure 7Zonally averaged DIC anomalies (mmol m−3) over (a) the Indo-Pacific and (b) Atlantic sectors and zonally averaged alkalinity anomalies (mmol m−3) over (c) the Indo-Pacific and (d) the Atlantic sectors for lambfe3% compared to 70ka-control.


Figure 8Zonally averaged (a, b) remineralised carbon (mmol m−3), (c, d) carbon due to calcite dissolution (mmol m−3), and (e, f) preformed carbon (mmol m−3) anomalies for lambfe3% compared to 70ka-control in the Indo-Pacific (a, d, g, j) and Atlantic c, f, i, l) sectors. Please note that the colour scale is different for panels (c) and (d) compared to other panels.


In contrast, lower DIC is simulated north of 47 S in the upper Pacific ocean (above 2000 m) and within the North Atlantic Deep Water (NADW) and Antarctic Intermediate Water (AAIW) pathways in the Atlantic above 3000 m. Simulated changes in deep-ocean oxygen concentrations (not shown) are opposite to the DIC changes indicating an increase in re-mineralised carbon in the deep ocean (Fig. 8a and b). An increase in mixed-layer alkalinity (Fig. 7c and d) is simulated, partially due to the decrease in coccolithophores between 47 S and 50 N (Figs. 6d and A4d). As a result of these changes in surface CaCO3 formation and subsequent changes in the subduction and sinking of CaCO3 into deeper layers, reduced CaCO3 dissolution at depth leads to lower alkalinity in the northern Pacific, while higher CaCO3 dissolution in the Southern Ocean and in the deep Atlantic leads to an alkalinity increase there (Figs. 7c, d and 8c, d).

To quantify the processes leading to changes in oceanic DIC, we decompose the changes in DIC into their remineralised, carbonate, and preformed contributions (see Sect. 2). In our 70ka-control simulation, the remineralisation process leads to the generation of 131 mmol m−3 Creg in the deep (≥3000 m depth) Southern Ocean, which increases by 22.5 mmol m−3 in lambfe3% due to higher iron influx. Changes in the carbonate pump contribute marginally to the deep Southern Ocean DIC showing an increase of 3.16 mmol m−3 CCaCO3, whereas the preformed DIC contribution (Cpref) decreases by 0.66 mmol m−3 compared to 70ka-control. North of 30 S, Creg in lambfe3% increases by 14.7 mmol m−3 in the deep Indo-Pacific basin, while CCaCO3 and Cpref decrease by 2 and 2.6 mmol m−3, respectively (Fig. 8a, c and e). The deep Atlantic Ocean shows an increase in Creg of 16.3 mmol m−3 and in CCaCO3 of 2 mmol m−3 in the lambfe3% experiment, while Cpref decreases by 3.2 mmol m−3 (Fig. 8b, d and f).

4 Discussion

The colder climate at 70 ka leads to an overall decrease in NPP and EP compared to the PI period in our simulations (9.7 % and 3 %, respectively). The changes in NPP and EP are spatially heterogeneous and sensitive to changes in the length of the growing season, sea ice cover, and nutrient availability and are thus associated with shifts in plankton distribution. In the Southern Ocean, the ∼7.5 % reduction in NPP and 2 % reduction in EP might be due to a shorter growing season for diatoms and coccolithophores (Saini et al.2021), associated with a larger annual sea ice extent and lower SSTs at 70 ka. Both NPP and EP also decrease at low latitudes and mid-latitudes.

Alkenone fluxes (Lamy et al.2014; Martínez-García et al.2014), as well as opal and organic carbon fluxes (Thöle et al.2019; Amsler et al.2022), reconstructed from sub-Antarctic sediment cores suggest that EP was higher in the SAZ during MIS4 than present-day levels. In addition, bottom water oxygenation records indicate a deep-ocean oxygenation decrease during MIS4, which might suggest increased respired carbon storage in the deep ocean (Amsler et al.2022) but could also reflect a change in ocean dynamics and water residence time. It has been suggested that the sub-Antarctic EP increase during MIS4 was due to higher iron fluxes into the Southern Ocean. Higher dust deposition in the southern Pacific and the southern Atlantic at 70 ka is consistent with available palaeo-proxy records (Martínez-García et al.2011; Lambert et al.2012; Lamy et al.2014; Thöle et al.2019). At the same time, palaeoceanographic records from the AZ suggest a decrease in EP during MIS4 (Anderson et al.2009; Jaccard et al.2013; Studer et al.2015; Thöle et al.2019; Amsler et al.2022; Weber et al.2022).

In our 70 ka simulations with enhanced iron input, EP increases in the AZ and polar frontal zone (between SAF and APF) and decreases in the SAZ, in contrast with most palaeo-proxy records (Fig. A5b). Existing δ15N records suggest higher nutrient consumption at the MIS4–MIS5 transition in the AZ (Studer et al.2015; Ai et al.2020), a region where the supply of macro-nutrients to the surface is determined by ocean upwelling and mixing (Lefèvre and Watson1999; Tagliabue et al.2017). Studer et al. (2015) argue that the relative increase in nutrient utilisation in the AZ during MIS4 is due to the general decrease in the nitrate supply, resulting from greater isolation of the deep ocean during glacial periods (Francois et al.1997; Jaccard et al.2013; Sigman et al.2021), rather than iron fertilisation.

This contradicts our results, where we find a deepening of the winter mixed layer in the 70 ka simulations compared to the PI-control, which leads to a slight increase in macro-nutrients in the AZ. Higher iron supply within this diatom-dominated seasonal sea ice zone then leads to a greater utilisation of available nutrients, in agreement with previous studies (Abelmann et al.2006, 2015), resulting in greater EP within the AZ. With higher nutrient utilisation in the AZ, the nutrient advection into the SAZ is reduced, leading to a decrease in EP in the SAZ. The disagreement between the modelling results and palaeo-proxy records thus suggests that the AZ is not stratified enough in our 70 ka simulations. This could be due to a weakening and/or more equatorward position of the SH westerly winds during glacial times (e.g. Toggweiler et al.2006; Gray et al.2021), which are not included in our simulations as present-day surface winds are prescribed.

However, our simulated EP increase in the AZ during MIS4 compared to the PI period, and thus the increase in regenerated carbon storage in the deep ocean, aligns with some proxy records suggesting a decrease in deep-ocean oxygenation in the AZ (Jaccard et al.2016; Amsler et al.2022). Our EP increase is also consistent with the increase in opal flux north of the APF (Amsler et al.2022) at the MIS4 onset (Fig. A5b). In summary, while our EP responses in the SAZ and AZ are inconsistent with most of the palaeo records, the strengthening of the biological pump leading to greater storage of deep-ocean organic carbon is consistent with palaeoceanographic evidence. Since proxy data for dust flux and EP is limited to only a few marine sites and ice cores, additional palaeo-proxy records covering a larger area in the Southern Ocean during MIS4 are needed to better quantify the impact of iron fertilisation during the glaciation.

Our results suggest that enhanced iron input south of the Antarctic polar front, i.e. south of ∼47 S (Dong et al.2006; Sokolov and Rintoul2009; Giglio and Johnson2016), is responsible for 100 % of the simulated atmospheric CO2 drawdown. This is in agreement with the hypothesis put forward by Marinov et al. (2006), which suggests that changes in the efficiency of the biological pump in the Antarctic zone play a dominant role in CO2 drawdown. However, our results differ from Lambert et al. (2021), who find that ∼30 % of the simulated CO2 decrease during the last glacial termination is due to enhanced aeolian iron input into the northern Pacific in their simulations. While they use the same LGM dust flux, their iron input in the ocean differs from our study. Lambert et al. (2021) assume 3.5 wt % iron fraction in the aeolian dust flux, whereas we use a global map of percentage distribution of iron in dust (see Sect. 2). More importantly, their solubility factor is not constant but varies non-linearly with the magnitude of the dust flux.

While the two dust masks used here suggest an increase in aeolian iron deposition in the Atlantic sector of the Southern Ocean during glacial times compared to the PI period, the iron input in glacfe is lower in the southern Pacific compared to pife, whereas lambfe is higher in the Pacific sector. Despite these differences in spatial patterns, the two iron masks (glacfe and lambfe) with the same iron solubilities lead to similar decreases in atmospheric CO2 in our simulations. This indicates that changes in atmospheric CO2 are more dependent on changes in solubility, than on regional differences in aeolian iron fluxes in the Southern Ocean. The experiments using the glacfe mask are slightly more efficient in drawing down CO2. One of the major Antarctic Bottom Water formation sites is located in the Weddell Sea in our simulations, making the southern Atlantic sector a more efficient region for carbon sequestration. Higher iron dust input in the Weddell Sea in the glacfe experiments thus leads to greater CO2 drawdown compared to the lambfe experiments.

With no significant ocean circulation changes, changes in atmospheric CO2 are a linear function of the overall efficiency of the biological pump (P*, Fig. A7). In agreement with earlier studies (Matsumoto et al.2002; Tagliabue et al.2014; Lambert et al.2015; Yamamoto et al.2019), we find that the efficiency of the soft tissue pump increases south of 47 S, while it decreases in the lower latitudes due to reduced northward advection of surface nutrients. The increase in the vertical nutrient gradient leads to the eventual saturation of the ocean carbon uptake capacity.

We simulate a maximum of 16 %–18 % increase in global P* in our experiments, leading to a maximum drop of 19 to 21 ppm in CO2. The timescale of this CO2 drawdown in our idealised simulations, i.e. without transient changes in dust or transient changes in ice sheet volume, is ∼5000 years, which is consistent with the observed timescale of CO2 transitions during MIS4 (Bereiter et al.2015). If we define the plausible range of the large-scale average iron solubility in the ocean to be 1 % to 10 % at 70 ka, the corresponding range of atmospheric CO2 drawdown is ∼4 to 16 ppm. It is more likely that large-scale solubility in the Southern Ocean at 70 ka ranged between 3 % and 5 %, which leads in our simulations to a CO2 decline of 9 to 11 ppm (Fig. 4b).

We find that the biological response to changes in iron fertilisation depends not only on the iron solubility during glacial periods but also on the iron solubility during warm periods. Our results are based on the assumption that the global average iron solubility during warm periods equals ∼1 %. At higher initial values, the total potential drawdown of CO2 would be smaller. For example, for an assumed solubility during warm periods closer to ∼3 %, we simulate a range of CO2 changes between 6.4 ppm (no change in solubility and glacial fluxes based on Lambert et al. (2015)) and 16.4 ppm (change to 20 % solubility and glacial fluxes based on Ohgaito et al. (2018)). This range reduces to 6.9 to 14.4 ppm if the initial solubility was 5 % and to 6.7 to 6.9 ppm if the initial solubility was 20 %.

Previous modelling studies obtained a 2 ppm drop in CO2 using 2 % iron solubility in the Southern Ocean under PI boundary conditions (Tagliabue et al.2014), while under LGM boundary conditions, using solubility factors between 1 % and 2 %, a 2–28 ppm drop in CO2 was simulated by a hierarchy of models (Bopp et al.2003; Tagliabue et al.2009b; Lambert et al.2015; Muglia et al.2017; Yamamoto et al.2019). Using iron solubility factors of 3 % and 10 % under LGM conditions, Yamamoto et al. (2019) simulated a 15.6 ppm to 20 ppm CO2 reduction, while a solubility of 7 % with varying dust masks led to a decrease in CO2 by 12–29 ppm in Lambert et al. (2021). A 12 ppm decrease at 70 ka was simulated in a transient simulation performed with a model of intermediate complexity using 1 % solubility (Menviel et al.2012). Our results are consistent but in the lower range of previous modelling estimates. Although our lowest estimate of a 4 ppm CO2 reduction is in agreement with the studies mentioned above, the upper limit of 16 ppm is lower than some previous estimates (e.g. Yamamoto et al.2019; Lambert et al.2021). This might be due to the comparatively larger simulated decrease in EP north of 47 S and in the equatorial Pacific in our study compared to previous studies (Yamamoto et al.2019).

None of the previous modelling studies on iron fertilisation simulate coccolithophore and diatom abundances prognostically. By including four distinct classes of plankton in our model, we highlight the competitive dynamics between different major phytoplankton functional types for light and nutrient availability. Coccolithophores contribute to the total carbon export mainly in the polar frontal zone, while the diatom contribution is in the Antarctic zone. As previously mentioned, carbon export close to convection sites in the Southern Ocean can be more efficient in reducing atmospheric CO2. Furthermore, while both diatom and coccolithophores contribute to CO2 uptake in the ocean through photosynthesis, coccolithophores produce CaCO3-rich platelets, which reduce surface ocean alkalinity, thus reducing the CO2 uptake efficiency. The inclusion of this competition should be taken into account when investigating the impact of ecosystem changes on the global carbon cycle.

Our estimate of a 16 ppm CO2 decrease with 10 % iron solubility corresponds to half of the ∼32 ppm CO2 drop estimated for the MIS4 transition (Bereiter et al.2015). This suggests that enhanced aeolian iron input at the MIS4 transition could have played a significant role in the CO2 drawdown (Martínez-García et al.2014; Kohfeld and Chase2017), but other processes must also have contributed to the CO2 decline. For example, an AMOC weakening (Piotrowski et al.2005; Wilson et al.2015; O'Neill et al.2021) could have enhanced ocean stratification, contributing 15–30 ppm to the CO2 decrease during MIS4 (Thornalley et al.2013; Yu et al.2016; Menviel et al.2017). An increase in ocean stratification at the MIS4 transition due to changes in AABW properties (Adkins2013) could also have increased the deep-ocean carbon content (Menviel et al.2017), thus contributing to the atmospheric CO2 decline. Finally, both changes in oceanic circulation and a more efficient biological pump could have enhanced sediment calcite dissolution, thus increasing ocean alkalinity and enhancing the CO2 drawdown (Boyle1988; Archer and Maier-Reimer1994; Yu et al.2016; Kobayashi and Oka2018).

5 Conclusions

We use an Earth system model of intermediate complexity incorporating a newly developed ecosystem model to better constrain the oceanic carbon uptake due to enhanced aeolian iron fluxes at the MIS4 transition. Based on a series of 70 ka sensitivity experiments forced with two different aeolian iron flux masks and different iron solubility factors, we find that enhanced iron input south of 47 S could have led to a 4 to 16 ppm CO2 decline at 70 ka, with a most likely range of 9 to 11 ppm. Our results suggest that the ocean's capacity to take up carbon decreases with increasing iron input because enhanced nutrient utilisation in the Antarctic zone is compensated for by decreased nutrient availability at mid-latitudes and low latitudes. Enhanced aeolian iron input at the MIS4 transition could thus explain up to 50 % of the observed CO2 drop at that time.

Appendix A

Table A1Globally integrated NPP and EP values. Percentage changes from 70ka-control experiment are provided in brackets.

Download Print Version | Download XLSX

Figure A1Aeolian iron dust flux (µmol m−2 yr−1) for (a) pife1%, (b) lambfe10%-30S (iron flux is equal lambfe10% south of 30 S and pife1% elsewhere), (c) lambfe10%-40S (same as lambfe10%-30S but for 40 S), (d) lambfe10%-50S, (e) lambfe10%-60S, and (f) lambfe50%-47S (iron flux is equal to lambfe50% south of 47 S and pife1% elsewhere).

Figure A2Aeolian iron dust fluxes (µmol m−2 yr−1) for pife (a, d), lambfe (b, e), and glacfe (c, f) masks with 3 % (a–c) and 10% (d–f) solubility factors.

Figure A3Atlantic meridional streamfunction (Sv) (a) at 70 ka and (b) for the PI period.


Figure A4(a) Diatom and (b) coccolithophore abundance anomalies (g C m−2) at 70 ka compared to the PI period. (c) Diatom and (d) coccolithophore abundance anomalies (g C m−2) in lambfe3% compared to 70ka-control.

Figure A5(a) Surface nitrate concentrations (mmol m−3) for the 70ka-control experiment. The overlaid contours show the zero isoline of annual wind stress curl (black) and the zero contour of EP anomalies (red) between lambfe3% and 70ka-control (pife1%). Blue arrows represent surface currents. (b) Export production anomalies (g C m−2 yr−1) at 177.5 m depth for lambfe3% compared to PI-control (pife1%) with opal flux (triangles) proxy records from Amsler et al. (2022) (left of 50 E) and from Thöle et al. (2019) (right of 50 E) to show qualitative comparison between model and data. Dark (light) orange represents significantly higher (slightly higher) and dark (light) blue represents significantly lower (slightly lower) values at 70 ka compared to the PI period. Overlaid dashed contours represent modern SAF in black and APF in green based on the definition of Sokolov and Rintoul (2009).

Figure A6The lambfe3% to 70ka-control (pife1%) anomalies of zonally averaged (a–c) total phosphate and (d–f) regenerated phosphate concentration anomalies (mmol m−3) over (a, d) the Pacific, (b, e) Indian, and (c, f) Atlantic basins.


Figure A7Atmospheric CO2 (ppm) anomalies as a function of global P*, with a slope of −207 ppm (compared to −312 ppm as per Ito and Follows2005).


Figure A8(a) Ratio of EP to NPP averaged over the Southern Ocean (south of 30 S), (b) surface nitrate concentration averaged over 40–47 S (mmol m−3), and (c) surface nitrate concentration south of 47 S (mmol m−3) as a function of the aeolian iron flux (Gmol yr−1) into the Southern Ocean south of 47 S.


Figure A9Hovmöller diagrams of the proportion of ocean grid cells per latitude band for the 70ka-control experiment for which (a) light is limiting coccolithophore growth, (b) nitrate is limiting coccolithophore growth, (c) light is limiting diatom growth, (d) nitrate is limiting diatom growth, and (e) silicate is limiting diatom growth. In all panels, if shade is equal to 1, then all ocean grid cells at that latitude band are limited by the respective limiting factor; however, if shade is equal to 0.5, then half of the grid cells are limited in this way.


Figure A10Hovmöller diagrams of the lambfe20% minus 70ka-control anomaly in the proportion of ocean grid cells per latitude band for which (a) nitrate is limiting coccolithophore growth, (b) nitrate is limiting diatom growth, (c) silicate is limiting diatom growth, (d) either of the macro-nutrients is limiting coccolithophore growth, and (e) either of the macro-nutrients is limiting diatom growth.


Figure A11Ecosystem anomalies between lambfe3%-70ka-control and lambfe20%-70ka-control including EP (g C m−2 yr−1) at 177.5 m, NPP (g C m−2 yr−1), diatoms (g C m−2), coccolithophores (g C m−2), general phytoplankton types (g C m−2), and diazotrophs (g C m−2).

Data availability

All the final data from modelling simulations are published on the UNSW ResData repository (, Saini et al.2022).

Author contributions

HS performed all the simulations and analyses. KJM and LM supervised the research and assisted in the interpretation of the results. HS, LM, and KJM drafted the manuscript. KK provided scientific support for the KMBM3 model and comments on an advanced draft of the manuscript.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Climate of the Past. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Himadri Saini acknowledges funding from the University International Postgraduate Award (UIPA) scheme from UNSW. Karin Kvale is supported by the Global Change Through Time Programme administered by the New Zealand Ministry of Business, Innovation and Employment. All experiments were performed on the computational facility of National Computational Infrastructure (NCI) owned by the Australian National University through awards under the Merit Allocation Scheme and the UNSW HPC at NCI Scheme.

Financial support

Katrin J. Meissner and Laurie Menviel acknowledge support from the Australian Research Council (grant nos. DP180100048, DP180102357, SR200100008, and FT180100606).

Review statement

This paper was edited by Marisa Montoya and reviewed by Juan Muglia and two anonymous referees.


Abelmann, A., Gersonde, R., Cortese, G., Kuhn, G., and Smetacek, V.: Extensive phytoplankton blooms in the Atlantic sector of the glacial Southern Ocean, Paleoceanography, 21, PA1013,, 2006. a

Abelmann, A., Gersonde, R., Knorr, G., Zhang, X., Chapligin, B., Maier, E., Esper, O., Friedrichsen, H., Lohmann, G., Meyer, H., and Tiedemann, R.: The seasonal sea-ice zone in the glacial Southern Ocean as a carbon sink, Nat. Commun., 6, 1–13, 2015. a

Abe-Ouchi, A., Saito, F., Kawamura, K., Raymo, M. E., Okuno, J., Takahashi, K., and Blatter, H.: Insolation-driven 100,000-year glacial cycles and hysteresis of ice-sheet volume, Nature, 500, 190–193, 2013. a, b

Adkins, J. F.: The role of deep ocean circulation in setting glacial climates, Paleoceanography, 28, 539–561, 2013. a

Ahn, J. and Brook, E. J.: Atmospheric CO2 and climate on millennial time scales during the last glacial period, Science, 322, 83–85, 2008. a

Ai, X. E., Studer, A. S., Sigman, D. M., Martínez-García, A., Fripiat, F., Thöle, L. M., Michel, E., Gottschalk, J., Arnold, L., Moretti, S., Schmitt, M., Oleynik, S., Jaccard, S. L., and Haug, G. H.: Southern Ocean upwelling, Earth's obliquity, and glacial-interglacial atmospheric CO2 change, Science, 370, 1348–1352, 2020. a

Amsler, H. E., Thöle, L. M., Stimac, I., Geibert, W., Ikehara, M., Kuhn, G., Esper, O., and Jaccard, S. L.: Bottom water oxygenation changes in the southwestern Indian Ocean as an indicator for enhanced respired carbon storage since the last glacial inception, Clim. Past, 18, 1797–1813,, 2022. a, b, c, d, e, f, g, h

Anderson, R., Ali, S., Bradtmiller, L., Nielsen, S., Fleisher, M., Anderson, B., and Burckle, L.: Wind-driven upwelling in the Southern Ocean and the deglacial rise in atmospheric CO2, Science, 323, 1443–1448, 2009. a, b

Anderson, R. F., Barker, S., Fleisher, M., Gersonde, R., Goldstein, S. L., Kuhn, G., Mortyn, P. G., Pahnke, K., and Sachs, J. P.: Biological response to millennial variability of dust and nutrient supply in the Subantarctic South Atlantic Ocean, Philos. T. Roy. Soc. A, 372, 20130054,, 2014. a

Archer, D. and Maier-Reimer, E.: Effect of deep-sea sedimentary calcite preservation on atmospheric CO2 concentration, Nature, 367, 260–263, 1994. a, b

Aumont, O. and Bopp, L.: Globalizing results from ocean in situ iron fertilization studies, Global Biogeochem. Cy., 20, GB2017,, 2006. a

Aumont, O., Maier-Reimer, E., Blain, S., and Monfray, P.: An ecosystem model of the global ocean including Fe, Si, P colimitations, Global Biogeochem. Cy., 17, 1060,, 2003. a

Baker, A., Jickells, T., Witt, M., and Linge, K.: Trends in the solubility of iron, aluminium, manganese and phosphorus in aerosol collected over the Atlantic Ocean, Mar. Chem., 98, 43–58, 2006. a

Bassinot, F. C., Labeyrie, L. D., Vincent, E., Quidelleur, X., Shackleton, N. J., and Lancelot, Y.: The astronomical theory of climate and the age of the Brunhes-Matuyama magnetic reversal, Earth Planet. Sc. Lett., 126, 91–108, 1994.  a

Batchelor, C. L., Margold, M., Krapp, M., Murton, D. K., Dalton, A. S., Gibbard, P. L., Stokes, C. R., Murton, J. B., and Manica, A.: The configuration of Northern Hemisphere ice sheets through the Quaternary, Nat. Commun., 10, 1–10, 2019. a, b

Bereiter, B., Lüthi, D., Siegrist, M., Schüpbach, S., Stocker, T. F., and Fischer, H.: Mode change of millennial CO2 variability during the last glacial cycle associated with a bipolar marine carbon seesaw, P. Natl. Acad. Sci. USA, 109, 9755–9760, 2012. a

Bereiter, B., Eggleston, S., Schmitt, J., Nehrbass-Ahles, C., Stocker, T. F., Fischer, H., Kipfstuhl, S., and Chappellaz, J.: Revision of the EPICA Dome C CO2 record from 800 to 600 kyr before present, Geophys. Res. Lett., 42, 542–549, 2015. a, b, c, d, e

Berger, A.: Long-term variations of daily insolation and Quaternary climatic changes, J. Atmos. Sci., 35, 2362–2367, 1978. a, b

Bopp, L., Kohfeld, K. E., Le Quéré, C., and Aumont, O.: Dust impact on marine biota and atmospheric CO2 during glacial periods, Paleoceanography, 18, 1046,, 2003. a, b

Boyle, E. A.: The role of vertical chemical fractionation in controlling late Quaternary atmospheric carbon dioxide, J. Geophys. Res.-Oceans, 93, 15701–15714, 1988. a

Brovkin, V., Ganopolski, A., Archer, D., and Munhoven, G.: Glacial CO2 cycle as a succession of key physical and biogeochemical processes, Clim. Past, 8, 251–264,, 2012. a

Conway, T. M., Wolff, E. W., Röthlisberger, R., Mulvaney, R., and Elderfield, H.: Constraints on soluble aerosol iron flux to the Southern Ocean at the Last Glacial Maximum, Nat. Commun., 6, 1–9, 2015. a

De Deckker, P., Arnold, L. J., van der Kaars, S., Bayon, G., Stuut, J.-B. W., Perner, K., dos Santos, R. L., Uemura, R., and Demuro, M.: Marine Isotope Stage 4 in Australasia: a full glacial culminating 65,000 years ago–global connections and implications for human dispersal, Quaternary Sci. Rev., 204, 187–207, 2019. a

Dong, S., Sprintall, J., and Gille, S. T.: Location of the Antarctic polar front from AMSR-E satellite sea surface temperature measurements, J. Phys. Oceanogr., 36, 2075–2089, 2006. a

Eby, M., Zickfeld, K., Montenegro, A., Archer, D., Meissner, K., and Weaver, A. J.: Lifetime of anthropogenic climate change: Millennial time scales of potential CO2 and surface temperature perturbations, J. Climate, 22, 2501–2511,, 2009. a

Eggleston, S., Schmitt, J., Bereiter, B., Schneider, R., and Fischer, H.: Evolution of the stable carbon isotope composition of atmospheric CO2 over the last glacial cycle, Paleoceanography, 31, 434–452, 2016. a, b

EPICA Community Members, Augustin, L., Barbante, C., Barnes, P. R., Barnola, J. M., Bigler, M., Castellano, E., Cattani, O., Chappellaz, J., Dahl-Jensen, D., Delmonte, B., Dreyfus, G., Durand, G., Falourd, S., Fischer, H., Flückiger, J., Hansson, M. E., Huybrechts, P., Jugie, G., Johnsen, S. J., Jouzel, J., Kaufmann, P., Kipfstuhl, J., Lambert, F., Lipenkov, V. Y., Littot, G. C., Longinelli, A., Lorrain, R., Maggi, V., Masson-Delmotte, V., Miller, H., Mulvaney, R., Oerlemans, J., Oerter, H., Orombelli, G., Parrenin, F., Peel, D. A., Petit, J.-R., Raynaud, D., Ritz, C., Ruth, U., Schwander, J., Siegenthaler, U., Souchez, R., Stauffer, B., Steffensen, J. P., Stenni, B., Stocker, T. F., Tabacco, I. E., Udisti, R., van de Wal, R. S. W., van den Broeke, M., Weiss, J., Wilhelms, F., Winther, J.-G., Wolff, E. W., and Zucchelli, M.: Eight glacial cycles from an Antarctic ice core, Nature, 429, 623–628, 2004. a

Fanning, A. F. and Weaver, A. J.: An atmospheric energy-moisture balance model: Climatology, interpentadal climate change, and coupling to an ocean general circulation model, J. Geophys. Res.-Atmos., 101, 15111–15128,, 1996. a

Ferrari, R., Jansen, M. F., Adkins, J. F., Burke, A., Stewart, A. L., and Thompson, A. F.: Antarctic sea ice control on ocean circulation in present and glacial climates, P. Natl. Acad. Sci. USA, 111, 8753–8758, 2014. a

Francois, R., Altabet, M. A., Yu, E.-F., Sigman, D. M., Bacon, M. P., Frank, M., Bohrmann, G., Bareille, G., and Labeyrie, L. D.: Contribution of Southern Ocean surface-water stratification to low atmospheric CO2 concentrations during the last glacial period, Nature, 389, 929–935, 1997. a, b

Gebhardt, H., Sarnthein, M., Grootes, P. M., Kiefer, T., Kuehn, H., Schmieder, F., and Röhl, U.: Paleonutrient and productivity records from the subarctic North Pacific for Pleistocene glacial terminations I to V, Paleoceanography, 23, PA4212,, 2008. a

Giglio, D. and Johnson, G. C.: Subantarctic and polar fronts of the Antarctic Circumpolar Current and Southern Ocean heat and freshwater content variability: A view from Argo, J. Phys. Oceanogr., 46, 749–768, 2016. a

Grant, K., Rohling, E., Bar-Matthews, M., Ayalon, A., Medina-Elizalde, M., Ramsey, C. B., Satow, C., and Roberts, A.: Rapid coupling between ice volume and polar temperature over the past 150,000 years, Nature, 491, 744–747, 2012. a, b

Gray, W. R., Lavergne, C. d., Wills, R. C. J., Menviel, L., Spence, P., Holzer, M., Kageyama, M., and Michel, E.: Poleward shift in the Southern Hemisphere westerly winds synchronous with the deglacial rise in CO2, Paleoceanogr. Paleoclim., 38, e2023PA004666,, 2021. a

Heinze, C., Maier-Reimer, E., and Winn, K.: Glacial pCO2 reduction by the world ocean: Experiments with the Hamburg carbon cycle model, Paleoceanography, 6, 395–430, 1991. a

Hibler, W. D.: A dynamic thermodynamic sea ice model, J. Phys. Oceanogr., 9, 815–846, 1979. a

Hunke, E. C. and Dukowicz, J. K.: An elastic–viscous–plastic model for sea ice dynamics, J. Phys. Oceanogr., 27, 1849–1867, 1997. a

Ito, A., Myriokefalitakis, S., Kanakidou, M., Mahowald, N. M., Scanza, R. A., Hamilton, D. S., Baker, A. R., Jickells, T., Sarin, M., Bikkina, S., Gao, Y., Shelley, R. U., Buck, C. S., Landing, W. M., Bowie, A. R., Perron, M. M. G., Guieu, C., Meskhidze, N., Johnson, M. S., Feng, Y., Kok, J. F., Nenes, A., and Duce, R. A.: Pyrogenic iron: The missing link to high iron solubility in aerosols, Sci. Adv., 5, eaau7671,, 2019. a, b

Ito, T. and Follows, M. J.: Preformed phosphate, soft tissue pump and atmospheric CO2, J. Mar. Res., 63, 813–839, 2005. a, b, c

Jaccard, S. L., Haug, G. H., Sigman, D. M., Pedersen, T. F., Thierstein, H. R., and Rohl, U.: Glacial/interglacial changes in subarctic North Pacific stratification, Science, 308, 1003–1006, 2005. a

Jaccard, S. L., Hayes, C. T., Martinez-Garcia, A., Hodell, D. A., Anderson, R. F., Sigman, D. M., and Haug, G.: Two modes of change in Southern Ocean productivity over the past million years, Science, 339, 1419–1423, 2013. a, b, c

Jaccard, S. L., Galbraith, E. D., Martínez-García, A., and Anderson, R. F.: Covariation of deep Southern Ocean oxygenation and atmospheric CO2 through the last ice age, Nature, 530, 207–210, 2016. a, b

Jouzel, J., Masson-Delmotte, V., Cattani, O., Dreyfus, G., Falourd, S., Hoffmann, G., Minster, B., Nouet, J., Barnola, J.-M., Chappellaz, J., Fischer, H., Gallet, J. C., Johnsen, S., Leuenberger, M., Loulergue, L., Luethi, D., Oerter, H., Parrenin, F., Raisbeck, G., Raynaud, D., Schilt, A., Schwander, J., Selmo, E., Souchez, R., Spahni, R., Stauffer, B., Steffensen, J. P., Stenni, B., Stocker, T. F., Tison, J. L., Werner, M., and Wolff, E. W.: Orbital and millennial Antarctic climate variability over the past 800,000 years, Science, 317, 793–796, 2007. a

Kalnay, E., Kanamitsu, M., Kistler, R., Collins, W., Deaven, D., Gandin, L., Iredell, M., Saha, S., White, G., Woollen, J., and Zhu, Y.: The NCEP/NCAR 40-year reanalysis project, B. Am. Meteorol. Soc., 77, 437–472, 1996. a

Keller, D. P., Oschlies, A., and Eby, M.: A new marine ecosystem model for the University of Victoria earth system climate model, Geosci. Model Dev., 5, 1195–1220,, 2012. a

Khatiwala, S., Schmittner, A., and Muglia, J.: Air-sea disequilibrium enhances ocean carbon storage during glacial periods, Sci. Adv., 5, eaaw4981,, 2019. a, b

Kleman, J., Fastook, J., Ebert, K., Nilsson, J., and Caballero, R.: Pre-LGM Northern Hemisphere ice sheet topography, Clim. Past, 9, 2365–2378,, 2013. a

Kobayashi, H. and Oka, A.: Response of atmospheric pCO2 to glacial Changes in the Southern Ocean amplified by carbonate compensation, Paleoceanogr. Paleoclim., 33, 1206–1229, 2018. a

Kohfeld, K. E. and Chase, Z.: Temporal evolution of mechanisms controlling ocean carbon uptake during the last glacial cycle, Earth Planet. Sc. Lett., 472, 206–215, 2017. a, b, c, d

Kohfeld, K. E. and Harrison, S. P.: DIRTMAP: the geological record of dust, Earth-Sci. Rev., 54, 81–114, 2001. a

Kohfeld, K. E., Le Quéré, C., Harrison, S. P., and Anderson, R. F.: Role of marine biology in glacial-interglacial CO2 cycles, Science, 308, 74–78, 2005. a, b

Kohfeld, K. E., Graham, R., De Boer, A., Sime, L., Wolff, E., Le Quéré, C., and Bopp, L.: Southern Hemisphere westerly wind changes during the Last Glacial Maximum: paleo-data synthesis, Quaternary Sci. Rev., 68, 76–95, 2013. a, b

Kucera, M., Weinelt, M., Kiefer, T., Pflaumann, U., Hayes, A., Weinelt, M., Chen, M.-T., Mix, A. C., Barrows, T. T., Cortijo, E., Duprat, J., Juggins, S., and Waelbroeck, C.: Reconstruction of sea-surface temperatures from assemblages of planktonic foraminifera: multi-technique approach based on geographically constrained calibration data sets and its application to glacial Atlantic and Pacific Oceans, Quaternary Sci. Rev., 24, 951–998, 2005. a

Kvale, K., Keller, D. P., Koeve, W., Meissner, K. J., Somes, C. J., Yao, W., and Oschlies, A.: Explicit silicate cycling in the Kiel Marine Biogeochemistry Model version 3 (KMBM3) embedded in the UVic ESCM version 2.9, Geosci. Model Dev., 14, 7255–7285,, 2021. a, b, c, d

Kvale, K. F., Meissner, K., Keller, D., Eby, M., and Schmittner, A.: Explicit planktic calcifiers in the university of victoria earth system climate model, version 2.9, Atmos.-Ocean, 53, 332–350, 2015a. a, b

Kvale, K. F., Meissner, K., and Keller, D.: Potential increasing dominance of heterotrophy in the global ocean, Environ. Res. Lett., 10, 074009,, 2015b. a

Lambeck, K., Purcell, A., Zhao, J., and SVENSSON, N.-O.: The Scandinavian ice sheet: from MIS 4 to the end of the last glacial maximum, Boreas, 39, 410–435, 2010. a

Lambert, F., Delmonte, B., Petit, J.-R., Bigler, M., Kaufmann, P. R., Hutterli, M. A., Stocker, T. F., Ruth, U., Steffensen, J. P., and Maggi, V.: Dust-climate couplings over the past 800,000 years from the EPICA Dome C ice core, Nature, 452, 616–619, 2008. a

Lambert, F., Bigler, M., Steffensen, J. P., Hutterli, M., and Fischer, H.: Centennial mineral dust variability in high-resolution ice core data from Dome C, Antarctica, Clim. Past, 8, 609–623,, 2012. a, b, c

Lambert, F., Tagliabue, A., Shaffer, G., Lamy, F., Winckler, G., Farias, L., Gallardo, L., and De Pol-Holz, R.: Dust fluxes and iron fertilization in Holocene and Last Glacial Maximum climates, Geophys. Res. Lett., 42, 6014–6023, 2015. a, b, c, d, e, f

Lambert, F., Opazo, N., Ridgwell, A., Winckler, G., Lamy, F., Shaffer, G., Kohfeld, K., Ohgaito, R., Albani, S., and Abe-Ouchi, A.: Regional patterns and temporal evolution of ocean iron fertilization and CO2 drawdown during the last glacial termination, Earth Planet. Sc. Lett., 554, 116675,, 2021. a, b, c, d, e

Lamy, F., Gersonde, R., Winckler, G., Esper, O., Jaeschke, A., Kuhn, G., Ullermann, J., Martínez-Garcia, A., Lambert, F., and Kilian, R.: Increased dust deposition in the Pacific Southern Ocean during glacial periods, Science, 343, 403–407, 2014. a, b, c, d, e, f

Lefèvre, N. and Watson, A. J.: Modeling the geochemical cycle of iron in the oceans and its impact on atmospheric CO2 concentrations, Global Biogeochem. Cy., 13, 727–736, 1999. a

Lüthi, D., Le Floch, M., Bereiter, B., Blunier, T., Barnola, J.-M., Siegenthaler, U., Raynaud, D., Jouzel, J., Fischer, H., Kawamura, K., and Stocker, T. F.: High-resolution carbon dioxide concentration record 650,000–800,000 years before present, Nature, 453, 379–382, 2008. a

Maher, B., Prospero, J., Mackie, D., Gaiero, D., Hesse, P., and Balkanski, Y.: Global connections between aeolian dust, climate and ocean biogeochemistry at the present day and at the last glacial maximum, Earth-Sci. Rev., 99, 61–97, 2010. a

Mahowald, N. M., Muhs, D. R., Levis, S., Rasch, P. J., Yoshioka, M., Zender, C. S., and Luo, C.: Change in atmospheric mineral aerosols in response to climate: Last glacial period, preindustrial, modern, and doubled carbon dioxide climates, J. Geophys. Res.-Atmos., 111, D10202,, 2006. a

Marinov, I., Gnanadesikan, A., Toggweiler, J., and Sarmiento, J. L.: The southern ocean biogeochemical divide, Nature, 441, 964–967, 2006. a, b

Martin, J. H.: Iron Hypothesis of CO2 Change, Paleoceanography, 5, 1–13, 1990. a

Martínez-García, A., Rosell-Melé, A., Jaccard, S. L., Geibert, W., Sigman, D. M., and Haug, G. H.: Southern Ocean dust–climate coupling over the past four million years, Nature, 476, 312–315, 2011. a, b, c

Martínez-García, A., Sigman, D. M., Ren, H., Anderson, R. F., Straub, M., Hodell, D. A., Jaccard, S. L., Eglinton, T. I., and Haug, G. H.: Iron fertilization of the Subantarctic Ocean during the last ice age, Science, 343, 1347–1350, 2014. a, b, c, d, e, f

Matsumoto, K., Sarmiento, J. L., and Brzezinski, M. A.: Silicic acid leakage from the Southern Ocean: A possible explanation for glacial atmospheric pCO2, Global Biogeochem. Cy., 16, 5–1–5–23,, 2002. a

Meissner, K., Weaver, A. J., Matthews, H. D., and Cox, P. M.: The role of land surface dynamics in glacial inception: A study with the UVic Earth System Model, Clim. Dynam., 21, 515–537,, 2003. a, b, c

Meissner, K., McNeil, B. I., Eby, M., and Wiebe, E. C.: The importance of the terrestrial weathering feedback for multimillennial coral reef habitat recovery, Global Biogeochem. Cy., 26, 1–20,, 2012. a

Mengis, N., Keller, D. P., MacDougall, A. H., Eby, M., Wright, N., Meissner, K. J., Oschlies, A., Schmittner, A., MacIsaac, A. J., Matthews, H. D., and Zickfeld, K.: Evaluation of the University of Victoria Earth System Climate Model version 2.10 (UVic ESCM 2.10), Geosci. Model Dev., 13, 4183–4204,, 2020. a

Menking, J. A., Shackleton, S. A., Bauska, T. K., Buffen, A. M., Brook, E. J., Barker, S., Severinghaus, J. P., Dyonisius, M. N., and Petrenko, V. V.: Multiple carbon cycle mechanisms associated with the glaciation of Marine Isotope Stage 4, Nat. Commun., 13, 1–10, 2022. a, b, c

Menviel, L., Joos, F., and Ritz, S.: Simulating atmospheric CO2, 13C and the marine carbon cycle during the Last Glacial–Interglacial cycle: possible role for a deepening of the mean remineralization depth and an increase in the oceanic nutrient inventory, Quaternary Sci. Rev., 56, 46–68, 2012. a, b, c

Menviel, L., Mouchet, A., Meissner, K. J., Joos, F., and England, M. H.: Impact of oceanic circulation changes on atmospheric δ13CO2, Global Biogeochem. Cy., 29, 1944–1961, 2015. a

Menviel, L., Yu, J., Joos, F., Mouchet, A., Meissner, K., and England, M. H.: Poorly ventilated deep ocean at the Last Glacial Maximum inferred from carbon isotopes: A data-model comparison study, Paleoceanography, 32, 2–17, 2017. a, b, c

Muglia, J., Somes, C. J., Nickelsen, L., and Schmittner, A.: Combined effects of atmospheric and seafloor iron fluxes to the glacial ocean, Paleoceanography, 32, 1204–1218,, 2017. a, b, c

Nickelsen, L., Keller, D. P., and Oschlies, A.: A dynamic marine iron cycle module coupled to the University of Victoria Earth System Model: The Kiel Marine Biogeochemical Model 2 for UVic 2.9, Geosci. Model Dev., 8, 1357–1381,, 2015. a

Ohgaito, R., Abe-Ouchi, A., O'Ishi, R., Takemura, T., Ito, A., Hajima, T., Watanabe, S., and Kawamiya, M.: Effect of high dust amount on surface temperature during the Last Glacial Maximum: A modelling study using MIROC-ESM, Clim. Past, 14, 1565–1581,, 2018. a, b

Oka, A., Abe-Ouchi, A., Chikamoto, M. O., and Ide, T.: Mechanisms controlling export production at the LGM: Effects of changes in oceanic physical fields and atmospheric dust deposition, Global Biogeochem. Cy., 25,, 2011. a

O'Neill, C. M., Hogg, A. McC., Ellwood, M. J., Opdyke, B. N., and Eggins, S. M.: Sequential changes in ocean circulation and biological export productivity during the last glacial–interglacial cycle: a model–data study, Clim. Past, 17, 171–201,, 2021. a

Pacanowski, R.: MOM2 documentation user's guide and reference manual: GFDL ocean group technical report, Geophysical Fluid Dynamics Laboratory (GFDL), National Oceanic and Atmospheric Administration, Princeton, NJ, (last access: 27 July 2023), 1995. a

Petit, J.-R., Jouzel, J., Raynaud, D., Barkov, N. I., Barnola, J.-M., Basile, I., Bender, M., Chappellaz, J., Davis, M., Delaygue, G., Delmotte, M., Kotlyakov, V. M., Legrand, M., Lipenkov, V. Y., Lorius, C., Pépin, L., Ritz, C., Saltzman, E., and Stievenard, M.: Climate and atmospheric history of the past 420,000 years from the Vostok ice core, Antarctica, Nature, 399, 429,, 1999. a

Piotrowski, A. M., Goldstein, S. L., Hemming, S. R., and Fairbanks, R. G.: Temporal relationships of carbon cycling and ocean circulation at glacial boundaries, Science, 307, 1933–1938, 2005. a, b

Saini, H., Kvale, K., Chase, Z., Kohfeld, K. E., Meissner, K. J., and Menviel, L.: Southern Ocean ecosystem response to Last Glacial Maximum boundary conditions, Paleoceanogr. Paleoclim., 36, e2020PA004075,, 2021. a, b, c, d, e, f, g, h

Saini, H., Meissner, K. J., Menviel, L., and Kvale, K.: Impact of iron fertilisation on atmospheric CO2 during the glaciation, UNSW ResData [data set],, 2022. a

Schmittner, A., Oschlies, A., Giraud, X., Eby, M., and Simmons, H. L.: A global model of the marine ecosystem for long-term simulations: Sensitivity to ocean mixing, buoyancy forcing, particle sinking, and dissolved organic matter cycling, Global Biogeochem. Cy., 19, 1–17,, 2005. a

Schroth, A. W., Crusius, J., Sholkovitz, E. R., and Bostick, B. C.: Iron solubility driven by speciation in dust sources to the ocean, Nat. Geosci., 2, 337–340, 2009. a

Semtner Jr., A. J.: A model for the thermodynamic growth of sea ice in numerical investigations of climate, J. Phys. Oceanogr., 6, 379–389, 1976. a

Shoenfelt, E. M., Sun, J., Winckler, G., Kaplan, M. R., Borunda, A. L., Farrell, K. R., Moreno, P. I., Gaiero, D. M., Recasens, C., Sambrotto, R. N., and Bostick, B. C.: High particulate iron (II) content in glacially sourced dusts enhances productivity of a model diatom, Sci. Adv., 3, e1700314,, 2017. a

Shoenfelt, E. M., Winckler, G., Lamy, F., Anderson, R. F., and Bostick, B. C.: Highly bioavailable dust-borne iron delivered to the Southern Ocean during glacial periods, P. Natl. Acad. Sci. USA, 115, 11180–11185, 2018. a

Sigman, D. M. and Boyle, E. A.: Glacial/interglacial variations in atmospheric carbon dioxide, Nature, 407, 859–869, 2000. a

Sigman, D. M., Hain, M. P., and Haug, G. H.: The polar ocean and glacial cycles in atmospheric CO2 concentration, Nature, 466, 47–55, 2010. a

Sigman, D. M., Fripiat, F., Studer, A. S., Kemeny, P. C., Martínez-García, A., Hain, M. P., Ai, X., Wang, X., Ren, H., and Haug, G. H.: The Southern Ocean during the ice ages: A review of the Antarctic surface isolation hypothesis, with comparison to the North Pacific, Quaternary Sci. Rev., 254, 106732,, 2021. a

Sokolov, S. and Rintoul, S. R.: Circumpolar structure and distribution of the Antarctic Circumpolar Current fronts: 1. Mean circumpolar paths, J. Geophys. Res.-Oceans, 114, C11019,, 2009. a, b, c, d

Stephens, B. B. and Keeling, R. F.: The influence of Antarctic sea ice on glacial–interglacial CO2 variations, Nature, 404, 171–174, 2000. a

Studer, A. S., Sigman, D. M., Martínez-García, A., Benz, V., Winckler, G., Kuhn, G., Esper, O., Lamy, F., Jaccard, S. L., Wacker, L., Oleynik, S., Gersonde, R., and Haug, G. H.: Antarctic Zone nutrient conditions during the last two glacial cycles, Paleoceanography, 30, 845–862, 2015. a, b, c, d

Tagliabue, A., Bopp, L., and Aumont, O.: Evaluating the importance of atmospheric and sedimentary iron sources to Southern Ocean biogeochemistry, Geophys. Res. Lett., 36, L13601,, 2009a. a

Tagliabue, A., Bopp, L., Roche, D. M., Bouttes, N., Dutay, J.-C., Alkama, R., Kageyama, M., Michel, E., and Paillard, D.: Quantifying the roles of ocean circulation and biogeochemistry in governing ocean carbon-13 and atmospheric carbon dioxide at the last glacial maximum, Clim. Past, 5, 695–706,, 2009b. a

Tagliabue, A., Aumont, O., and Bopp, L.: The impact of different external sources of iron on the global carbon cycle, Geophys. Res. Lett., 41, 920–926, 2014. a, b, c, d

Tagliabue, A., Bowie, A. R., Boyd, P. W., Buck, K. N., Johnson, K. S., and Saito, M. A.: The integral role of iron in ocean biogeochemistry, Nature, 543, 51–59, 2017. a

Thöle, L. M., Amsler, H. E., Moretti, S., Auderset, A., Gilgannon, J., Lippold, J., Vogel, H., Crosta, X., Mazaud, A., Michel, E., Martínez-García, A., and Jaccard, S. L.: Glacial-interglacial dust and export production records from the Southern Indian Ocean, Earth Planet. Sc. Lett., 525, 115716,, 2019. a, b, c, d, e, f

Thornalley, D. J., Barker, S., Becker, J., Hall, I. R., and Knorr, G.: Abrupt changes in deep Atlantic circulation during the transition to full glacial conditions, Paleoceanography, 28, 253–262, 2013. a, b

Toggweiler, J.: Origin of the 100,000-year timescale in Antarctic temperatures and atmospheric CO2, Paleoceanography, 23, PA2211,, 2008. a

Toggweiler, J. R., Russell, J. L., and Carson, S. R.: Midlatitude westerlies, atmospheric CO2, and climate change during the ice ages, Paleoceanography, 21, PA2005,, 2006. a

Waelbroeck, C., Labeyrie, L., Michel, E., Duplessy, J. C., McManus, J., Lambeck, K., Balbon, E., and Labracherie, M.: Sea-level and deep water temperature changes derived from benthic foraminifera isotopic records, Quaternary Sci. Rev. 21, 295–305, 2002. a

Watson, A. J., Vallis, G. K., and Nikurashin, M.: Southern Ocean buoyancy forcing of ocean ventilation and glacial atmospheric CO2, Nat. Geosci., 8, 861–864, 2015. a

Weaver, A. J., Eby, M., Wiebe, E. C., Bitz, C. M., Duffy, P. B., Ewen, T. L., Fanning, A. F., Holland, M. M., MacFadyen, A., Matthews, H. D., Meissner, K., Saenko, O., Schmittner, A., Wang, H., and Yoshimori, M.: The UVic earth system climate model: Model description, climatology, and applications to past, present and future climates, Atmos.-Ocean, 39, 361–428,, 2001. a

Weber, M. E., Bailey, I., Hemming, S. R., Martos, Y. M., Reilly, B. T., Ronge, T. A., Brachfeld, S., Williams, T., Raymo, M., Belt, S. T., et al.: Antiphased dust deposition and productivity in the Antarctic Zone over 1.5 million years, Nat. Commun., 13, 2044,, 2022. a, b

Williams, R. G. and Follows, M. J.: Ocean dynamics and the carbon cycle: Principles and mechanisms, Cambridge University Press, ISBN 9780511977817,, 2011. a

Wilson, D. J., Piotrowski, A. M., Galy, A., and Banakar, V. K.: Interhemispheric controls on deep ocean circulation and carbon chemistry during the last two glacial cycles, Paleoceanography, 30, 621–641, 2015. a

Wolff, E., Barbante, C., Becagli, S., Bigler, M., Boutron, C., Castellano, E., De Angelis, M., Federer, U., Fischer, H., Fundel, F., Hansson, M., Hutterli, M., Jonsell, U., Karlin, T., Kaufmann, P., Lambert, F., Littot, G. C., Mulvaney, R., Röthlisberger, R., Ruth, U., and Wegner, A.: Changes in environment over the last 800,000 years from chemical analysis of the EPICA Dome C ice core, Quaternary Sci. Rev., 29, 285–295, 2010.  a

Wolff, E. W., Fischer, H., Fundel, F., Ruth, U., Twarloh, B., Littot, G. C., Mulvaney, R., Röthlisberger, R., de Angelis, M., Boutron, C. F., Hansson, M., Jonsell, U., Hutterli, M. A., Lambert, F., Kaufmann, P., Stauffer, B., Stocker, T. F., Steffensen, J. P., Bigler, M., Siggaard-Andersen, M. L., Udisti, R., Becagli, S., Castellano, E., Severi, M., Wagenbach, D., Barbante, C., Gabrielli, P., and Gaspari, V.: Southern Ocean sea-ice extent, productivity and iron flux over the past eight glacial cycles, Nature, 440, 491–496, 2006. a

Yamamoto, A., Abe-Ouchi, A., Ohgaito, R., Ito, A., and Oka, A.: Glacial CO2 decrease and deep-water deoxygenation by iron fertilization from glaciogenic dust, Clim. Past, 15, 981–996,, 2019. a, b, c, d, e, f, g

Yu, J., Menviel, L., Jin, Z., Thornalley, D., Barker, S., Marino, G., Rohling, E., Cai, Y., Zhang, F., Wang, X., Dai, Y., Chen, P., and Broecker, W. S.: Sequestration of carbon in the deep Atlantic during the last glaciation, Nat. Geosci., 9, 319–324, 2016. a, b, c, d

Zhang, Y., Mahowald, N., Scanza, R. A., Journet, E., Desboeufs, K., Albani, S., Kok, J. F., Zhuang, G., Chen, Y., Cohen, D. D., Paytan, A., Patey, M. D., Achterberg, E. P., Engelbrecht, J. P., and Fomba, K. W.: Modeling the global emission, transport and deposition of trace elements associated with mineral dust, Biogeosciences, 12, 5771–5792,, 2015. a, b

Short summary
Understanding the changes in atmospheric CO2 during the last glacial cycle is crucial to comprehend the impact of climate change in the future. Previous research has hypothesised a key role of greater aeolian iron input into the Southern Ocean in influencing the global atmospheric CO2 levels by impacting the changes in the marine phytoplankton response. In our study, we test this iron hypothesis using climate modelling and constrain the impact of ocean iron supply on global CO2 decrease.