Transient Response of Southern Ocean Ecosystems During Heinrich Stadials

Antarctic ice core records suggest that atmospheric CO2 increased by 15–20 ppm during Heinrich stadials (HS). These periods of abrupt CO2 increase are associated with a significant weakening of the Atlantic meridional overturning circulation (AMOC), and a warming at high southern latitudes. As such, modeling studies have explored the link between changes in AMOC, high southern latitude climate and atmospheric CO2. While proxy records suggest that the aeolian iron input to the Southern Ocean decreased significantly during HS, the potential impact on CO2 of reduced iron input combined with oceanic circulation changes has not been studied in detail. Here, we quantify the respective and combined impacts of reduced iron fertilization and AMOC weakening on CO2 by performing numerical experiments with an Earth system model under boundary conditions representing 40,000 years before present (ka). Our study indicates that reduced iron input can contribute up to 6 ppm increase in CO2 during an idealized Heinrich stadial. This is caused by a 5% reduction in nutrient utilization in the Southern Ocean, leading to reduced export production and increased carbon outgassing from the Southern Ocean. An AMOC weakening under 40ka conditions and without changes in surface winds leads to a ∼0.5 ppm CO2 increase. The combined impact of AMOC shutdown and weakened iron fertilization is almost linear, leading to a total CO2 increase of 7 ppm. Therefore, this study highlights the need of including changes in aeolian iron input when studying the processes leading to changes in atmospheric CO2 concentration during HS.


Introduction
Paleoproxy records from Greenland ice cores have shown the occurrence of abrupt climate variability during the last glacial period (Dansgaard et al., 1982;Johnsen et al., 1972;Kindler et al., 2014).Such variability, referred to as Dansgaard-Oeschger (D/O) variability, is characterized by an abrupt rise in Greenland air temperatures on decadal timescales (D/O event), followed by a gradual cooling over hundreds or thousands of years to cold (stadial) stages (Bond et al., 1993;Dansgaard et al., 1993).About 25 such climatic oscillations have been recorded in Greenland ice cores (Kindler et al., 2014) out of which 14 D/O events occurred during Marine Isotope Stage 3 (MIS 3, 65-27 thousand years before present (ka BP hereafter)).D/O stadials are associated with layers of ice rafted debris in North Atlantic marine sediment cores, indicating iceberg discharges into the ocean from the surrounding ice sheets (Van Kreveld et al., 2000).Some stadials are characterized by especially thick layers of ice rafted debris originating from the Laurentide ice-sheet (Broecker et al., 1992;Elliot et al., 1998;Grousset et al., 2000;Heinrich, 1988;Snoeckx et al., 1999), and are referred to as Heinrich Events (HEs).
Observations suggest that concurrent to each HS in the North Atlantic during MIS 3, there was an Antarctic warm event (also known as Antarctic Isotope Maxima, AIM) (Bauska et al., 2021;Buizert et al., 2015).Antarctic ice core records also reveal that there was a rise in atmospheric CO 2 concentration during these events (Ahn & • Simulated marine ecosystems and CO 2 respond linearly to the combined impact of freshwater addition and reduced iron flux into the ocean • A maximum of 7 ppm increase in CO 2 is simulated due to the AMOC shutdown and weakened iron fertilization under MIS 3 boundary conditions Brook, 2008Brook, , 2014;;Bauska et al., 2021;Bereiter et al., 2015;Parrenin et al., 2013).For example, during Heinrich stadial 4 (H4), which is associated with an intense ice rafted debris deposition in the North Atlantic (Hemming, 2004), the temperature in Greenland decreased by ∼5°C (Huber et al., 2006;Kindler et al., 2014), while Antarctic temperatures rose by 2-3°C (Jouzel et al., 2007), and global CO 2 concentration increased by 15-20 ppm (Ahn & Brook, 2014;Bauska et al., 2018;Bereiter et al., 2015) (Figures 1a-1c).
A compilation of numerical studies by Gottschalk et al. (2019) suggests that simulated atmospheric CO 2 levels can vary between 7 and +14 ppm in response to an AMOC shutdown under interglacial boundary conditions, whereas under glacial boundary conditions, this range can vary between 5 and +6 ppm.Modeling studies have further suggested that the rise in CO 2 during HS was caused by a decrease in deep Pacific and Southern Ocean carbon resulting from stronger North Pacific Deep Water (NPDW) formation and enhanced Southern Ocean outgassing (Huiskamp & Meissner, 2012;Menviel et al., 2014Menviel et al., , 2018)).Enhanced Southern Ocean outgassing  (Kindler et al., 2014); (b) Antarctic deuterium (δD) concentration (‰) as recorded in European Project for Ice coring in Antarctica (EPICA) DOME C ice core (Jouzel et al., 2007), representing temperature variability in Antarctica; (c) atmospheric CO 2 concentrations (ppm) from EPICA DOME C ice core (Bereiter et al., 2015); (d) dust flux (mg/m 2 /yr) from EPICA DOME C ice core (Lambert et al., 2012); Gray shade represents Heinrich stadials during MIS 3 (65-27 ka).could result from stronger convection or enhanced wind-driven upwelling as also suggested from paleo-proxy ventilation records (Anderson et al., 2009;Skinner et al., 2020).
Marine sediment cores from the South Atlantic (Martínez-García et al., 2011) and the South Pacific (Lamy et al., 2014) as well as Antarctic ice cores (Fischer et al., 2010;Lambert et al., 2012) have highlighted millennialscale fluctuations in the atmospheric dust flux deposition at high southern latitudes that are opposite to changes in atmospheric CO 2 levels (Figure 1d).The associated changes in aeolian iron input to the ocean have been suggested to impact CO 2 concentration by modulating Southern Ocean nutrient utilization (Martin, 1990;Martínez-García et al., 2014).Foraminifera-bound δ 15 N (FB-δ 15 N) records (Martínez-García et al., 2014;Studer et al., 2015) suggest that the greater iron supply can stimulate higher nutrient utilization in the ocean, by reducing the iron limitation effect, and causing higher export production (EP).
Many modeling studies have assessed the impact of changes in iron input to the ocean on Southern Ocean EP and atmospheric CO 2 by performing sensitivity experiments under Last Glacial Maximum (Aumont et al., 2003;Bopp et al., 2003;Khatiwala et al., 2019;Lambert et al., 2015Lambert et al., , 2021;;Muglia et al., 2017;Oka et al., 2011;Saini et al., 2021;Yamamoto et al., 2019) and pre-industrial (PI) (Aumont & Bopp, 2006;Tagliabue et al., 2009Tagliabue et al., , 2014) ) boundary conditions.However, even though the link between aeolian iron input and Southern Ocean EP should also hold for millennial scale climate changes (Anderson et al., 2014), the response of marine ecosystems and ocean biogeochemistry to changes in iron dust flux during millennial-scale events of the last glacial period has not been simulated in detail.So far, only one experiment studied the impact of millennial-scale changes in dust flux on atmospheric CO 2 by forcing the Bern3D model with reconstructed dust deposition from DOME C ice core under PI boundary conditions (Parekh et al., 2008).They simulated the largest increase in atmospheric CO 2 (+9 ppm) during H4 as a response to the lowest reconstructed aeolian dust flux.
Alkenone based EP evidence from a marine sediment core east of New Zealand, indicates higher productivity during HS over the past 70,000 years in this region (Sachs & Anderson, 2005).This is in agreement with the higher amount of opal flux during HS obtained from four cores in the Antarctic zone (south of ∼60°S in the Pacific sector and south of ∼50°S in the Atlantic sector) (Anderson et al., 2009).However, recent studies (Anderson et al., 2014;Martínez-García et al., 2014) based on alkenone concentrations, opal and organic carbon fluxes measured in sub-Antarctic Atlantic sediment cores (close to ∼40°S), point toward a reduction in biological productivity during HS when the dust supply declined (Lambert et al., 2012).This is also in agreement with the records of higher deep ocean oxygenation during HS in the sub-Antarctic Atlantic (Gottschalk et al., 2016;Jaccard et al., 2016), which is an indicator of less remineralized organic carbon in the deep ocean.However, lower carbon sequestration in the ocean can also be due to increased ventilation as recorded from radiocarbon changes in the deep Southern Ocean during H4 (Gottschalk et al., 2016).Therefore, the Southern Ocean biological efficiency is impacted by more than just iron limitation.As there are only a few cores at sparse locations, it is thus challenging to draw robust conclusions about EP changes in the Southern Ocean during HS.
In addition, while many studies assessed the impact of an AMOC shutdown or an increase in aeolian iron input on the carbon cycle, no modeling study has yet looked into the combined impact of an AMOC shutdown and a decrease in aeolian iron flux on the carbon cycle under appropriate glacial boundary conditions.A study by Jochum et al. (2022) analyzed the impact of longer stadials and reduced iron flux on global atmospheric CO 2 concentration.However, they investigated this response with an experiment resembling LGM climate conditions and iron flux reduced by 50% from PI dust flux.Here, we investigate the response of marine ecosystems and atmospheric CO 2 to the individual and combined impacts of AMOC shutdown and reduced iron flux input in the Southern Ocean under MIS 3 boundary conditions.We focus on the changes in ecosystems and atmosphere-ocean CO 2 exchange in the Southern Ocean as the iron fertilization effect is crucial in this region (Saini et al., 2023).

Model
The experiments are performed with the version 2.9 University of Victoria Earth System Climate Model of intermediate complexity.It includes the MOM2 ocean general circulation model (Pacanowski, 1995) coupled to a sea ice model (Hibler, 1979;Hunke & Dukowicz, 1997;Semtner, 1976), a land surface scheme and a dynamic vegetation model (Meissner et al., 2003), a sediment model (Archer & Maier-Reimer, 1994;Meissner et al., 2012) with fixed weathering flux, and a vertically integrated 2-D energy and moisture balance atmospheric model (Fanning & Weaver, 1996).Precipitation occurs when relative humidity exceeds 85%.All model components are on a grid with a spatial resolution of 3.6 by 1.8°and there are 19 vertical depth levels in the ocean.Seasonally varying solar insolation (Berger, 1978) at the top of the atmosphere and seasonal variations in present day wind stress and wind fields (Kalnay et al., 1996) are used as external forcing to drive the model.An ideal age tracer is added to track water mass age (Koeve et al., 2015).A detailed description of the model structure and physics can be obtained from (Eby et al., 2009;Meissner et al., 2003;Mengis et al., 2020;Weaver et al., 2001).
The climate model incorporates the recently developed ecosystem model called Kiel Marine Biogeochemistry Model, v3 (Kvale et al., 2021;Saini et al., 2021), which represents the ocean carbon cycle.The ecosystem model was first developed by Ewen et al. (2004) and Schmittner et al. (2005) as a relatively simple marine inorganic and organic carbon cycle module, coupled to the UVic ESCM.The module has then been refined in Schmittner and Galbraith (2008); Keller et al. (2012); Kvale et al. (2015aKvale et al. ( , 2021)).
This model now includes prognostic equations for diatoms (silicifying plankton) and coccolithophores (calcifying plankton), in addition to general phytoplankton and diazotrophs.It also includes a prognostic iron cycle (Nickelsen et al., 2015) featuring hydrothermal sources (Yao et al., 2019), and fully incorporated silica and CaCO 3 cycles.The concentrations of dissolved iron and particulate iron in each grid box are determined by the physical transport and source minus sink terms (Equations 1-3 in Nickelsen et al., 2015).The physical transport includes advection, diapycnal and isopycnal diffusion and convection.Particulate iron sinks vertically with a velocity that progressively increases with depth.Iron that settles on the ocean floor is either promptly remineralized or gets incorporated into the sediments, depending on the oxygen state of the surrounding water masses with a cut off value of 5 mmol m 3 .The mean residence time of iron is 38 years in this model under PI conditions.
We employ a subgrid bathymetric parameterization corresponding to LGM sea-level, as described by Somes et al. (2017), to depict processes associated with iron cycling at the ocean's seabed.It is integrated into the biogeochemical component of the ocean model, ensuring that it is consistent with the model's grid for oceanic physics.

Experimental Design
Using this latest version, described in Kvale et al. (2021), we performed an equilibrated PI control simulation (Saini et al., 2021) with orbital parameters corresponding to year 1800 and a CO 2 forcing of 283.86 ppm.Here, we integrate a 40ka-control simulation (Figure 2), with orbital parameters (Berger, 1978), ice sheet topography (obtained from an offline ice sheet simulation (Abe-Ouchi et al., 2013)) and global mean alkalinity corresponding to 40ka BP, and CO 2 levels set at 200 ppm (Bereiter et al., 2015).The global mean ocean alkalinity is 0.1 mol m 3 higher in the 40ka-control simulation than in the PI simulation; this increase was determined based on estimates of sea level difference.
We use the LGM dust deposition map from Lambert et al. (2015) (Figure S1a in Supporting Information S1), which was derived by running a 2-day interpolation technique to globally distributed observational records of dust flux, the majority of which are included in the Dust Indicators and Records of Terrestrial and Marine Paleoenvironments database (Kohfeld & Harrison, 2001;Maher et al., 2010).
From this global dust map, we create the global maps of iron (lambfe) and silica (lambsi) fluxes by multiplying the LGM dust flux (Figure S1a in Supporting Information S1) with 2-D fields of iron/silica percentages in dust based on Zhang et al. (2015) (Figure S1 in Supporting Information S1).The iron and silica content in dust are modern estimates and they could have been different during different climates.Furthermore, we also note that using the dust deposition map based on the LGM (and not 40ka BP) records could be another source of uncertainty and highlight the need for additional proxy records with more spatial and temporal coverage during MIS3.For our 40ka-control experiment, we use the lambfe fluxes with an assumed 3% seawater solubility (Figure S1c in Supporting Information S1).
After 10,000 years of spin up, the model is run with prognostic CO 2 for another 1,000 years (Figure 3, black line).The PI simulation is performed with a PI iron mask with 1% iron solubility factor (pife1%), created from the BASE-PI dust deposition simulation from Mahowald et al. (2006) and the percentage map from Zhang et al. (2015) as described above.
To simulate an idealized HS, we add a meltwater pulse of 0.05 Sv (equivalent to a 2 m global sea level rise) into the North Atlantic (52°W:5.4°E ; 78°N:86°N) for 500 years (FW, Figure 3a, orange line, Phase 1).This meltwater causes an AMOC shutdown.The model is run for another 500 years without any meltwater input (Phase 2).To simulate an AMOC recovery, a salt flux addition to the North Atlantic is required in the UVic ESM.A negative meltwater input ( 0.15 Sv) is thus added in the North Atlantic for 200 years.The model is then run for another 300 years (Phase 3).
To understand the impact of reduced iron flux, we run another experiment called FE (Figure 3, red line), during which there is no freshwater forcing, but the global iron flux into the ocean is abruptly decreased by 90% to 2.744 Gmol yr 1 (Figure 3g).This experimental design was chosen based on European Project for Ice coring in Antarctica DOME C ice core records, where the optical particle counts as well as soluble calcium (Ca 2+ ) and nonsea salt soluble calcium (nssCa 2+ ) records indicate a reduction of dust flux by about 84% during H4, between 40ka and 38ka (Lambert et al., 2012).Furthermore, as per the records from the Berkner ice core (Conway et al., 2015) and two marine sediment cores (South Atlantic (Martínez-García et al., 2011); South Pacific (Lamy et al., 2014)) the iron flux decreased by ∼67%, ∼30% and ∼40%, respectively, during H4.Therefore, our approach provides an upper estimate of CO 2 change due to iron flux reduction.We run the model for 1,000 years and then switch the iron flux back to the original lambfe flux value (27.44 Gmol yr 1 ) in Phase 3.
Finally, in our third experiment, FWFE (Figure 3, green line), we analyze the impact of the combination of a freshwater discharge into the North Atlantic and a simultaneous reduction in the iron input.Therefore, the iron flux is reduced by 90% and 0.05 Sv of meltwater is added in the North Atlantic during the first 500 years (Phase 1).For the following 500 years, the iron flux stays reduced while no additional meltwater is added in the North Atlantic (Phase 2).In Phase 3, the salt is applied, the iron flux is changed back to the original value and the model is run for 200 years.After that the model continues to run for another 300 years with the same iron flux and no freshwater perturbation.
In each of these simulations, the North Atlantic Deep Water (NADW) transport is calculated as the maximum meridional Atlantic overturning streamfunction at 30°N between 1,500 and 3,500 m depth.Antarctic Bottom Water transport (AABW) is calculated as the minimum meridional global overturning streamfunction at 35°S below 3,000 m depth.North Pacific Deep Water transport is calculated as the maximum meridional Indo-Pacific overturning streamfunction at 30°N below 700 m depth.
In addition, simulated dissolved inorganic carbon (DIC) is decomposed into its three main contributors: respired organic carbon (C reg ), DIC generated by dissolution of calcium carbonate (C CaCO 3 ) and preformed carbon (C pref ) following the methods described below.

CO 2 Decomposition
Changes in sea surface pressure of CO 2 (pCO 2 ), which influence the magnitude and direction of air-sea flux exchange, are controlled by changes in sea surface temperature (SST), sea surface salinity (SSS), alkalinity (ALK) and DIC.
The sea surface pCO 2 anomalies can be decomposed into the individual contributions using the equations given below (Sarmiento & Gruber, 2006): The SST contribution is derived from: where γ SST = 0.0423 (Sarmiento & Gruber, 2006), and pCO ref 2 represents the sea surface pCO 2 value of the 40kacontrol simulation.

Simulated Oceanic Conditions at 40ka
We first briefly describe the simulated oceanic conditions for the 40ka-control experiment compared to the PI experiment.
The globally averaged ocean temperature is ∼1.2°C lower in the 40ka-control simulation compared to PI, and the globally averaged SSTs are ∼1.8°Clower.The SSTs are particularly lower at mid to high latitudes, being 1.5°C lower south of 30°S and 3°C lower north of 42°N (Figure 2a).The high latitude cooling is associated with an increase in sea ice cover: the annual sea ice edge in the Atlantic and Indian sectors of the Southern Ocean is located 4°northward compared to PI, and 2°northward in the Pacific sector of the Southern Ocean.In the North Atlantic, the sea ice cover is ∼6°southward compared to PI, and ∼10°southward in the sea of Okhotsk.The simulated NADW transport in the 40ka-control simulation is 8 Sv, compared to 12 Sv in the PI simulation.
The globally colder and drier conditions at 40ka combined with the lower CO 2 concentration result in a total loss of terrestrial carbon by 424 gigatons (GtC) (Table 1), out of which 153 GtC is from vegetation and 271 GtC from soil carbon changes.
Due to higher iron availability, the diatom concentration south of the annual sea ice edge and coccolithophores north of the annual sea ice edge increase at 40ka compared to PI (Figures 2b and 2c), leading to a 6% increase in the total Southern Ocean EP (Table 1).Overall, the ocean carbon budget is increased by approximately 750 GtC at 40ka compared to PI (Table 1).This is due to a more efficient biological pump, higher ocean alkalinity (+0.1 mol m 3 ) and colder seawater temperatures.

Carbon Cycle Response to an AMOC Shutdown
As a response to the freshwater discharge in the North Atlantic in experiment FW (Figure 3a) during Phase 1, the NADW transport decreases gradually from a maximum of ∼8 Sv and shuts off completely after 400 years of meltwater addition (Figure 3c) and until the end of Phase 2. NPDW transport strengthens by ∼3 Sv over 1,000 years, while AABW transport weakens initially by ∼1.As described elsewhere (Okazaki et al., 2010;Saenko et al., 2004), a weakening in NADW transport is accompanied by an enhanced NPDW formation in the North Pacific and AAIW formation in the Southern Ocean due to an increase in surface salinity arising from oceanic and atmospheric teleconnections with the North Atlantic (Figure 4b).
In agreement with earlier studies (Menviel et al., 2020), an inter-hemispheric thermal seesaw pattern, with a 3.3°d ecrease in North Atlantic (35°W:10°W, 47°N:55°N) sea surface temperatures (SSTs) and a 0.4°increase in Southern Ocean SST, is simulated during the AMOC shutdown (Figure 4a).There is no significant change in the Southern Hemisphere annual mean sea ice edge, which is situated at 53°S in the Indian and Atlantic sectors and at 62°S in the Pacific sector (Figure 4a).
The combination of greater ocean stratification (attributable to higher SSTs and SSS), moderate increase in AAIW formation and the reduction in AABW transport lead to less surface nutrient availability during Phase 1.Therefore, there is an overall decrease in both diatoms (by 3%) and coccolithophores (by 6%) abundance in the Southern Ocean during Phase 1 (Figures 3h-3j).
However, during Phase 2, a further increase in AAIW formation, combined with a slow recovery of AABW leads to an increase in ventilation and deepening of the mixed layer depth (MLD) in the Southern Ocean (Figure 4c), which causes greater nutrient availability particularly in the Pacific sector.The locally increased nutrients are quickly taken up by the diatoms (being prescribed higher growth rates and half saturation constants for nutrients than coccolithophores (Kvale et al., 2021;Saini et al., 2021)), increasing diatom's abundance by ∼16% in the Pacific sector of the Southern Ocean, while coccolithophores decline (Figures 5a and 5b).This increase in diatoms is responsible for the total Southern Ocean diatom increase (+6.7%) shown in Figure 3h during Phase 2. On the contrary, coccolithophores decrease by ∼9.6% in the Southern Ocean as they are outcompeted by diatoms (Figures 5b and 3j).Overall, EP is 8.6% lower (Figures 3i and 5c) and net primary productivity (NPP) 5% lower in the Southern Ocean (not shown) during the AMOC shutdown.
In contrast, significantly weaker NADW transport and higher stratification, thus higher residence time, leads to a greater amount of regenerated phosphate (P reg ) at intermediate depths in the North Atlantic and in the deep Atlantic ocean (Figure S2 in Supporting Information S1).Correspondingly, DIC, in the form of regenerated carbon (C reg ) (Figure S3 in Supporting Information S1), accumulates in the Atlantic Ocean.Positive DIC anomalies are simulated below 500 m depth in the North Atlantic and below 1500 m depth in the South Atlantic (Figure 6a).At the end of Phase 2, the Atlantic Ocean carbon reservoir has thus increased by 140 GtC (Table 2).On the contrary, the Pacific Ocean loses 120 GtC, because of increased ventilation associated with stronger NPDW and AAIW transports (Figures 3e, 3f, and 4d-4f, Table 2).
Enhanced intermediate depth ventilation in the Southern Ocean (Figure 3f), and to some extent lower Southern Ocean EP, leads to a decrease in DIC, and particularly in C reg , in the upper 1000 m of the Southern Ocean (Figures 6a-6c, Figures S3a-S3c in Supporting Information S1).This anomaly spreads at intermediate depths within AAIW (Figures 6a-6c).The positive DIC anomalies in the deep Pacific and Indian Ocean (Figures 6b and  6c) are mainly due to increased preformed carbon (C pref ) resulting from the lower nutrient utilization at high southern latitudes (Figures S2 and S3d-S3f in Supporting Information S1).In the Indian Ocean, an increase in C reg is also simulated below 2500 m depth (Figure S3c in Supporting Information S1).Overall, DIC decreases by 19 GtC in the Indian Ocean and by 9 GtC in the Southern Ocean (Table 2, Figure 6c).
At the surface of the Southern Ocean, between ∼47°S and 60°S, higher DIC concentration due to reduction in EP and higher ventilation, increases the surface ocean partial pressure of CO 2 (pCO 2 , Figure 7b).Conversely, an increase in surface alkalinity in the Southern Ocean (Figure S4 in Supporting Information S1), associated with the decrease in coccolithophore population (Figure 5b), partly compensates the effect of higher DIC and ventilation on surface ocean pCO 2 (Figure 7b).In addition, higher SST and SSS lead to a reduction of CO 2 solubility in sea water (Figure 7b).Therefore, higher surface DIC and a weakened solubility pump cause a CO 2 outgassing in the Southern Ocean (red regions south of 30°S in Figure 8b), and particularly in the southeast Pacific between ∼47°F S:60°S.In contrast, the increase in surface alkalinity leads to carbon drawdown from the atmosphere into the Southern Ocean (blue regions south of 40°S in Figure 8b).Overall, the Southern Ocean carbon outgassing decreases by ∼8%.We also see an anomalous outgassing from the northeast Atlantic (Figure 8b), which is due to a decrease in surface alkalinity (Figure 7a) in response to the freshwater flux addition in this region.
The AMOC shutdown in our study leads to colder and drier conditions over the North Atlantic (Figures S5a and S5b in Supporting Information S1), extending into northern Europe and northwest Asia (40°N to 80°N), and to warmer and wetter conditions in the Southern Hemisphere (20°S to 40°S; Figures S5a and S5b in Supporting Information S1).The changes in precipitation in the tropics are small as the model cannot simulate a shift in the inter-tropical convergence zone (ITCZ).Due to the colder and drier conditions, vegetation carbon decreases by 4.5 GtC north of 18°N (Figure S5d in Supporting Information S1).The fraction of needle-leaf trees in northern  (Sachs & Anderson, 2003, 2005) and organic carbon flux (square) (Anderson et al., 2014;Thöle et al., 2019).Dark (light) orange represents significantly (slightly) higher values, while dark (light) blue represents significantly (slightly) lower values.
Europe, C4 grasses in Europe and mid-latitude Asia, and shrub in the high northern latitudes is reduced (resulting in an approximately 6.3%, 8.8%, and 9% loss in total vegetation carbon, respectively) while C3 grasses expand and contribute +4.7% to the total vegetation carbon change.C3 grasses and broadleaf trees increase between 10°S and 18°N to the detriment of C4 grasses.Vegetation carbon in that region increases by ∼4% for broadleaf trees, ∼1.7% for C3 grasses, and decreases by ∼3% for C4 grasses.This leads to a net increase of ∼3.8GtC in total vegetation carbon in the tropical regions (Figure S5d in Supporting Information S1).South of 10°S, the slightly warmer conditions lead to slightly wetter conditions and a 2.2 GtC vegetation carbon increase (Figures S5a, S5b, and S5d in Supporting Information S1).Globally, vegetation carbon increases by about 1.5 GtC (0.3%) in the FW simulation compared to the 40ka-control.The sign of the soil carbon anomaly mostly follows the sign of the vegetation carbon anomaly.For example, the increase in litter fall in the tropical sector leads to a ∼1.7% increase in soil carbon in this region (Figure S5e in Supporting Information S1) and the soil carbon decrease in northern Europe is due to a decrease in the fraction and litterfall of needle-leaf trees (not shown).One exception is eastern Asia, where there is a 3.6% soil carbon increase while vegetation carbon decreases.This is due to a reduction in soil respiration, resulting from colder conditions (Figure S5b in Supporting Information S1) and from an increase in soil moisture above optimal conditions (Figure S5c in Supporting Information S1).Globally, soil carbon increases by 6 GtC in the FW simulation (0.6% compared to 40ka-control).
The changes in carbon exchange with the atmosphere causes the global ocean carbon to decrease by 8 GtC by the end of Phase 2 (Table 2, Figure 3l).This is accompanied by a 7.5 GtC increase in terrestrial carbon reservoir (Figure 3k), and about 1 GtC carbon uptake by the atmosphere leading to an insignificant rise in CO 2 (Figure 3b) at the end of Phase 2.
During the AMOC recovery in Phase 3, the NADW transport increases to ∼14 Sv before equilibrating at 7 Sv once the salt flux is turned off (Figure 3c).NPDW transport also reduces back to near initial conditions, while AABW transport stabilizes to a slightly higher value, at 11 Sv (Figures 3d and 3e).The addition of the salt flux causes a ∼1 ppm rise in CO 2 , after which it reduces back to the initial value.
Diatoms increase by 21% compared to Phase 2 in the Indian Ocean sector south of the annual sea ice edge, and coccolithophores increase by 6% mostly everywhere north of the annual sea ice edge (Figures S6a and S6b in Supporting Information S1).This increases the total Southern Ocean EP by 7.6% at the end of Phase 3 compared to Phase 2 and brings EP back to a similar value as in 40ka-control.

Carbon Cycle Response to a Reduction in Aeolian Iron Flux
When the aeolian iron input to the ocean is abruptly reduced (Figure 3g), the CO 2 increases by 6 ppm over 1,000 years in experiment FE (Figure 3b).There are no changes in ocean circulation in this experiment, and therefore the changes in CO 2 are due to the reduction in the biological pump efficiency and the land carbon changes responding to higher CO 2 and temperatures.Plankton abundances south of ∼47°S are reduced due to decreased iron flux in the Southern Ocean (Figures 5d and 5e).Diatoms decrease by 19% in the Atlantic, 5.7% in the Pacific and ∼9.5% in the Indian sector south of ∼47°S, while coccolithophores decrease by 27% south of ∼47°S .This leads to reductions in EP (Figure 5f) and NPP south of ∼47°S.
Reduced nutrient utilization within the Antarctic zone causes an advection of excess nutrients northward (Figure S7 in Supporting Information S1), increasing diatoms (by 88% in the Indian region) and coccolithophores (by 23% in total) abundances and thus EP and NPP north of 47°S in the Southern Ocean (Figures 5d-5f).Thereby, EP decreases by 17% south of 47°S and increases by 29% between 47°S and 30°S as a response to nutrient reorganization, with an overall 5.5% decrease in diatoms and 1.6% decrease in coccolithophores in the Southern Ocean.The total Southern Ocean EP decreases by 7% at the end of Phase 2 (Figure 3i) due to a 5% decrease in nutrient utilization in this sector.
The lower EP south of 47°S causes an increase in surface DIC.However, this also leads to reduced DIC, in the form of regenerated DIC (C reg ) within AABW (Figure 6d-6f and Figures S8a-S8c in Supporting Information S1).These negative DIC anomalies are further advected into the deep Atlantic and Indo-Pacific basins (Figures 6d-6f).In contrast, greater DIC content is simulated north of 47°S in the upper 2500 m depth in all the basins.There is also a small increase in surface alkalinity in the Southern Ocean (Figure S4 in Supporting Information S1) due to the overall decrease in Southern Ocean coccolithophores (Figure 3j).
No significant changes in surface alkalinity, SST and SSS are simulated between 47°S and 60°S (Figure 7c).In contrast, higher surface DIC leads to ∼48% increase in CO 2 outgassing between 47°S and 60°S (Figure 7c and red regions south of 40°S in Figure 8c).This leads to a 50 GtC carbon loss from the Southern Ocean in experiment FE (Table 2).
In this simulation, there are no significant changes in global precipitation (Figure S5f in Supporting Information S1).Changes in terrestrial carbon are primarily driven by higher temperatures (Figure S5g in Supporting Information S1) and the CO 2 fertilization effect.Elevated CO 2 levels and higher temperatures increase net canopy photosynthesis, leading to an increase in both broadleaf and needleleaf trees, each contributing to a global vegetation carbon increase of ∼3%, with some increase (∼0.7%) in the northern high latitudes due to C4 grasses.
There is also a reduction in C3 grasses in the northern high latitudes, which is compensated by an increase in shrubs in that region.This results in a 10 GtC rise in total vegetation carbon (+2% compared to 40ka-control; Figure S5i in Supporting Information S1).These changes also contribute to a 2% increase in litter fall globally, leading to a subsequent 13 GtC increase in soil carbon in our simulation (1.3% compared to 40ka-control, Figure S5j in Supporting Information S1).
Overall, the ocean loses 35 GtC (Figure 3l), the atmosphere gains 12 GtC, and 23 GtC are absorbed by the terrestrial biosphere (Figure 3k).During Phase 3, CO 2 decreases gradually as the ocean biogeochemistry recovers to the initial conditions of 40ka-control.

Carbon Cycle Response to an AMOC Shutdown and Iron Flux Reduction
In this section we analyze the combined effect of an AMOC shutdown due to freshwater input in the North Atlantic (Figure 3a, experiment FWFE), and a reduced global iron flux (Figure 3g, experiment FWFE).A 7 ppm CO 2 increase is simulated in this experiment (1 ppm higher than FE and 6.5 ppm higher than FW).
As can be seen in Figures 3c-3f, the changes in ocean circulation in FWFE are similar to those of the FW experiment, with slightly larger anomalies in AABW transport and slightly smaller anomalies in NPDW transport.Figures 3 and 5 clearly show that the diatom, coccolithophore and EP responses in FWFE experiment are a linear combination of the responses simulated in FW and FE.
Diatoms increase south of 47°S in the Pacific sector due to circulation changes as in experiment FW.However, because of the parallel reduction in iron supply (which causes 6% diatom reduction in this sector), the net increase is limited to 9.5% (instead of 16% as simulated in FW) compared to 40ka-control.Similarly, the diatoms decrease south of 47°S in the Atlantic ( 16%, +3% from circulation and 19% from the iron effect) and Indian sectors ( 14%, 3.2% from circulation and 9.5% from the iron effect) and the increase north of 47°S in the Indian sector (+63%) is also a net outcome from changes in both circulation and iron supply.Consequently, due to the weak positive anomalies in the Pacific, and the stronger negative anomalies in the Atlantic and Indian sectors, the Southern Ocean diatom abundance is 6% lower compared to FW but 5.8% higher compared to FE (Figure 3h) at the end of Phase 2. Coccolithophores on the other hand decrease by 11.4% over the Southern Ocean (Figure 3j) compared to the 40ka-control (where 9.6% decrease is simulated due to the circulation changes as in FW and 1.6% is simulated due to changes in iron as in FE).With these changes in the plankton population, the total Southern Ocean EP reduces by 14% in FWFE at the end of Phase 2 (Figure 5i).
The larger simulated decrease in C reg in the FWFE simulation (Figures S9a-S9c in Supporting Information S1) compared to that in FW and FE, due to the weaker ( 13%) Southern Ocean nutrient utilization (Figure S10 in Supporting Information S1), causes a larger decrease in DIC below 200 m depth in the Southern Ocean (Figures 6g-6i).Simulated air-sea CO 2 fluxes in the FWFE experiment (Figures 7 and 8d) as well as the changes in ocean carbon (Table 2) in the Southern Ocean are a combination of FW and FE carbon changes.
The cumulative impact of the processes described for the FW and FE simulations above leads to a 12 GtC increase in vegetation carbon (+3% compared to 40ka-control) and a 19 GtC increase in soil carbon (+3% compared to 40ka-control) in the FWFE simulation globally (Figures S5n and S5o in Supporting Information S1).In addition to the increase in broadleaf trees and C3 grasses in the tropics, there is also an increase in C3 grasses in eastern and north-eastern Asia.Furthermore, in western North America, needleleaf trees expand and replace shrubs, contributing to 20% of the total vegetation carbon increase north of 18°N in FWFE.
Therefore, a total of ∼45 GtC are lost from the ocean in experiment FWFE (Table 2, Figure 3l).Terrestrial carbon increases by 31 GtC (Figure 3k), while atmospheric carbon increases by 14 GtC by the end of Phase 2.

Discussion
Our study indicates that to understand the mechanisms behind atmospheric CO 2 increase during HS, an investigation of dust related impacts is crucial in addition to previously explored feedbacks from oceanic circulation changes.
In our idealized experiment simulating a HS under MIS 3 boundary conditions, the simulated SST anomalies in the North Atlantic (∼ 3°C) are consistent with the alkenone based SST proxies from the Iberian Margin (Martrat et al., 2007) and the North Atlantic (Cacho et al., 1999).On the contrary, the simulated Southern Ocean SST anomalies (+0.4°C) are significantly lower than what is inferred from marine records (+∼2°C) (Barker et al., 2009;Caniupán et al., 2011;Kaiser et al., 2005;Pahnke et al., 2003) but are still within the range of previously simulated responses to an AMOC shutdown (Kageyama et al., 2013;Menviel et al., 2014Menviel et al., , 2015)).
The changes in global circulation (Figure S11 in Supporting Information S1), and therefore the changes in Southern Ocean physics, impact marine ecosystems in this region.As pointed out earlier, high resolution paleo evidence of changes in EP during HS are limited.The simulated decrease in coccolithophores and EP at ∼45°S east of the South Island of New Zealand in our FW simulation (Figures 5b and 5c) is in contradiction with the increase in alkenone (proxy for EP) and C 37 methyl ketones (a marker for coccolithophores) concentrations recorded in a marine sediment core, MD97-2120 (Sachs & Anderson, 2005), at this location near Chatham Rise during H4 (Figure 5c).Our experiment FWFE, which includes the combined impact from the iron flux reduction (Figures 5e and 5f) and freshwater forcing (Figures 5b and 5c), improves this model data agreement by simulating an increase in coccolithophores and EP at this location (Figures 5h and 5i).The simulated increase in diatom abundance east of New Zealand, due to deeper MLD in experiments FW and FWFE, is in agreement with brassicasterol concentrations (a marker for diatoms) (Figures 5a-5g) observed in the Chatham rise core (Sachs & Anderson, 2005).Increased diatom abundance south of ∼50°S in the Atlantic in FW and FWFE and a decrease in diatoms between ∼50°S and ∼55°S in the Indian sector in FE and FWFE, are also in agreement with the marine sediment opal flux records from the Southern Ocean (Amsler et al., 2022;Anderson et al., 2009;Thöle et al., 2019).In the sub-Antarctic, simulated diatom changes north of ∼50°S in FW are consistent with the proxy records from Atlantic (Anderson et al., 2014;Sachs & Anderson, 2003) and Indian sectors (Amsler et al., 2022;Thöle et al., 2019), while in the FE and FWFE simulations, the diatom increase is consistent with only one record (Amsler et al., 2022) in the Indian sector.
The simulated EP increase at ∼45°S in simulations FE and FWFE is in contradiction with the reduced alkenone and organic carbon fluxes (Anderson et al., 2014;Sachs & Anderson, 2003) recorded in marine sediment cores in the subantarctic Atlantic.However, the simulated EP decrease south of 40°S in experiment FW is in agreement with these paleo EP records.The simulated increase in EP in the Indian sector is consistent with the organic carbon flux proxy (Thöle et al., 2019).In experiment FWFE, there is a 60% agreement between the changes in simulated diatom abundance and the available proxy records, while there is a 40% agreement between changes in EP and the proxy records.
The simulated EP decrease in FW, along with the simulated decline south of ∼47°S in FE and FWFE experiments, is in agreement with the EP response to an AMOC shutdown obtained with another version of the UVic model (Menviel et al., 2014).Based on the above-mentioned proxy evidence and our simulated results, we suggest that an EP decrease due to reduced iron supply during HS could most likely have occurred south of ∼45°S in the Southern Ocean.Clearly, better constraints on changes in EP and plankton abundances during HS, along with additional proxy records of the last glacial millennial-scale climate variability are required.
Enhanced AAIW transport and reduction in Southern Ocean EP lead to an increase in intermediate and deep ocean oxygen concentrations south of ∼40°S in all of our simulations (Figure S12 in Supporting Information S1).This is in agreement with the reconstructed bottom water oxygen concentrations (Gottschalk et al., 2016), and authigenic Uranium based deep oxygenation records (Jaccard et al., 2016) in sub-Antarctic Atlantic marine sediment cores at ∼45°S and ∼55°S respectively.
A total of 7 ppm CO 2 increase is simulated in experiment FWFE due to the compounding effects of changes in oceanic circulation and a reduced iron dust flux in the Southern Ocean.Most of this CO 2 increase is due to the reduced iron availability (experiment FE).This is caused by a larger carbon loss (50 GtC compared to 9 GtC) from the Southern Ocean in experiment FE compared to FW.In simulation FW, outgassing in the North Atlantic is compensated by the higher carbon sequestration (140 GtC in FW compared to 8 GtC in FE, Table 2) in the deep North Atlantic Ocean due to the AMOC shutdown, and by the decrease in outgassing in the Southern Ocean.
Earlier model studies did not simulate a consistent terrestrial carbon reservoir response to a shutdown of the AMOC.While some studies simulate an overall terrestrial carbon decrease (Jochum et al., 2022;Menviel et al., 2008;Nielsen et al., 2019), some simulate an overall terrestrial carbon increase (Huiskamp & Meissner, 2012).This is not surprising given that the processes included in the models differ. in each grid box and do not compete with each other.These models cannot represent regional shifts in vegetation types.Our study and Huiskamp and Meissner (2012) do not include a dynamic atmosphere, and therefore cannot represent shifts of the ITCZ and associated changes in precipitation.All of these experiments started with different initial climate conditions (preindustrial, 12ka BP, Last Glacial Maximum, 40ka BP) and used different forcings for the transient simulations.While our simulated precipitation patterns broadly align with earlier modeling studies (Bozbiyik et al., 2011;Jochum et al., 2022;Menviel et al., 2008;Nielsen et al., 2019), the precipitation anomalies are much smaller than those simulated in dynamic atmosphere models (Kageyama et al., 2013).
There is a 23 GtC increase in terrestrial carbon in our FE simulation associated with a 6 ppm increase in CO 2 and a 0.15°C global and annual mean warming.A recent study by Keenan et al. (2023) indicates that a 80 ppm increase in CO 2 during the modern climate period (1981-2020) might have resulted in a 16 GtC increase in terrestrial carbon.Compared to this study, our simulations show a stronger CO 2 fertilization effect.This is due to the nonlinear dependency of net canopy photosynthesis with CO 2 (Cox et al., 1998) with a steeper slope for lower CO 2 values.The dependency of photosynthesis on temperature is also non-linear (Cox et al., 1998).The CO 2 fertilization effect is therefore highly dependent on the background climate state and the impact of CO 2 fertilization during glacial times should be further investigated.
In agreement with Jochum et al. (2022), the rate of change in atmospheric CO 2 in simulations FE and FWFE (Phase 1) decreases over time.While the decrease in iron flux leads to a decrease in C reg in the ocean, this reduction in C reg is partly compensated by an increase in oceanic C pref (Huiskamp et al., 2016) leading to a total CO 2 increase of 6 ppm in FE.
The 7 ppm rise in FWFE is significantly lower than the maximum observed CO 2 increase of 20 ppm during HS within MIS 3 (Ahn & Brook, 2008;Bauska et al., 2021;Bereiter et al., 2015).It is also lower than the simulated increase of 9 ppm simulated by Parekh et al. (2008) as a response to reduced dust flux, although using a much simpler biogeochemical model than our study and under PI conditions.Jochum et al. ( 2022) simulated an increase of ∼8 ppm in CO 2 in their experiment with a longer stadial and a reduction in PI iron flux by 50%.
Previous modeling studies have suggested that enhanced Southern Ocean ventilation could have contributed to the atmospheric CO 2 increase (Menviel et al., 2014(Menviel et al., , 2015) ) during HS.Changes in Southern Ocean ventilation could result from intensified Southern Hemisphere westerlies (Anderson et al., 2009;Huiskamp et al., 2016;Menviel et al., 2018;Sigman et al., 2021;Toggweiler et al., 2006;Tschumi et al., 2008).In addition, recent studies have also highlighted the role of changes in AAIW in driving atmospheric CO 2 increases (Menviel et al., 2018;Yu et al., 2022).Due to the coarse resolution of the ocean model and lack of dynamic atmosphere (and thus constant Southern Hemisphere westerlies) in the experiments presented above, it is possible that changes in Southern Ocean circulation in our study are not appropriately captured.The major Southern Ocean processes impacting atmospheric CO 2 are therefore related to iron fertilization in our study with no significant contribution from circulation changes.
There is also uncertainty involving the total amount of dissolved iron available for uptake by phytoplankton under different climates.In the experiments presented here, we use 3% iron solubility, which is based on the most plausible estimate of glacial iron solubility (3%-5%) (Conway et al., 2015;Saini et al., 2023;Shoenfelt et al., 2018).Based on the relationship between iron flux and CO 2 change presented in Saini et al. (2023), we estimate that if the 40ka-control was run with a solubility of 5%-10%, a reduction in iron flux such as the one performed in FE would lead to a 10-13 ppm CO 2 increase.Furthermore, the use of a uniform mean solubility is a simplification and does not reflect regional variations in iron solubility due to different regional dust characteristics and sources.We acknowledge this as a limitation and a caveat in our research.
We also do not consider the impact of potential changes in ligand concentrations in the ocean which can also impact iron solubility (Parekh et al., 2006(Parekh et al., , 2008)).Additionally, better constraints on the pattern and amount of dust fluxes under different climates are required to constrain the contribution of iron fertilization during millennial events.

Conclusions
Our study quantifies the individual and combined effects of reduced aeolian iron input and weaker AMOC on atmospheric CO 2 concentration and ecosystems during the last glacial period to understand the processes leading to the observed atmospheric CO 2 increase during HS.A total increase of 7 ppm in atmospheric CO 2 is simulated

Figure 3 .
Figure 3. Transient evolution of simulated biological and physical variables.Time series of (a) freshwater flux (Sv) applied to the North Atlantic and (g) global iron flux (Gmol yr 1 ), used as forcings in experiments (orange) FW, (red) FE, and (green) FWFE.Please note that the green line overlaps the orange line in (a) and the red line overlaps the green line in (g).The black line represents 200 years of 40ka-control experiment with prognostic CO 2 .Time series of simulated (b) global atmospheric CO 2 concentration (ppm) anomaly compared to 40ka-control, (c) North Atlantic Deep Water (Sv), (d) AABW (Sv), and (e) North Pacific Deep Water transport (Sv), (f) Southern Ocean ideal age anomaly averaged over intermediate depths (600-1700 m, years) compared to 40ka-control; time series anomalies of simulated (h) Southern Ocean depth integrated diatoms (PgC), (i) Southern Ocean export production at 177.5 m depth (PgC), (j) Southern Ocean depth integrated coccolithophores (PgC), (k) global terrestrial carbon (GtC) and (l) global ocean carbon (GtC) compared to 40ka-control.The Southern Ocean (SO in the figure) is defined as the region south of 30°S, and all the biological tracers are integrated over the Southern Ocean.

igure 4 .
FW experiment at year 1,000 (end of Phase 2) minus 40ka-control anomalies of (a) sea surface temperature (°C), (b) sea surface salinity (psu), (c) ventilation depth (m), and zonally averaged ideal age (years) over the (d) Atlantic, (e) the Pacific, and (f) Indian sector; where red in panels (d-f) indicates reduction in ventilation.Overlaid contours in (a) represent the 15% annual mean sea ice concentration in (red) 40ka-control and in (green) FW experiment at the end of Phase 2.

Figure 6 .
Figure 6.Zonally averaged dissolved inorganic carbon anomalies (mmol m 3 ) over (left) the Atlantic, (center) the Pacific, and (right) Indian ocean for experiments (top) FW, (center) FE, and (bottom) FWFE at year 1,000 (end of Phase 2) compared to the 40ka-control simulation.Please note that the color scale is non linear.

Figure 7 .
Figure 7. Time series of sea surface pressure of CO 2 (pCO 2 ) (ppm) anomalies resulting from changes in (blue) ALK, (orange) dissolved inorganic carbon (DIC), (green) sea surface salinity (SSS), and (red) SST for (a) FW averaged over the North Atlantic (north of 35°N), (b) FW averaged over the Southern Ocean (47°S:60°S), (c) FE over the Southern Ocean (47°S:60°S), and (d) FWFE averaged over the Southern Ocean (47°S:60°S), compared to the 40ka-control simulation.The dashed black line represents simulated sea surface pCO 2 (ppm), while the solid black line shows the sum of the ALK, DIC, SSS and SST contributions.Please note the scale differs between the sub-panels.
For example, Nielsen et al. (2019) and Jochum et al. (2022) do not include a dynamic vegetation model.Their vegetation types are fixed

Table 1
Simulated Export Production and Carbon Reservoirs in 40ka-Control and Pre-Industrial Simulations SAINI ET AL.

Table 2
Changes in Atmospheric CO 2 , Atmospheric, Terrestrial and Ocean Carbon at the End of Phase 2 (Yr 1000) Compared to the 40ka-Control Simulation an Earth System model of intermediate complexity when the AMOC is shutdown and the iron flux to the Southern reduced.About 86% of this CO 2 increase is due to the reduced iron input, which leads to reduced nutrient utilization in the Southern Ocean, and an increased carbon outgassing from this region.Himadri Saini would like to acknowledge funding from the University International Postgraduate Award scheme at University of New South Wales.Model runs were undertaken with the assistance of resources from the National Computational Infrastructure (NCI), which is supported by the Australian Government.Laurie Menviel acknowledges funding from the Australian Research Council (ARC) grants FT180100606, DP220102134 and SR200100008.Karin Kvale acknowledges support from the New Zealand Ministry of Business, Innovation and Employment (MBIE) within the Antarctic Science Platform, grant ANTA1801.Karin Kvale also acknowledges support from the New Zealand MBIE through the Global Change through Time programme (Strategic Science Investment Fund, contract C05X1702).Katrin J. Meissner acknowledges funding from the ARC grants DP180100048 and DP180102357 and high-performance computing resources through the Merit Allocation Scheme of the National Computational Infrastructure (NCI) and the UNSW HPC at NCI Scheme.
SAINI ET AL. in