Antarctic Icebergs Reorganize Ocean Circulation During Pleistocene Glacials

in

The dominant feature of large-scale mass transfer in the modern ocean is the Atlantic meridional overturning circulation (AMOC).The geometry and vigour of this circulation influences global climate on various timescales.Paleoceanographic evidence suggests that during glacial periods of the past 1.5 million years the AMOC was markedly different from modern 1 : in the Atlantic basin, deep waters of Southern Ocean origin increased in volume whilst above them the core of the North Atlantic Deep Water (NADW) shoaled 2 .An absence of evidence to elucidate the origin of this phenomenon means that the sequence of events during the descent into global glacial conditions remains unclear.Here, we present multi-proxy evidence showing that northward shifts in Antarctic iceberg melt in the Indian-Atlantic Southern Ocean systematically preceded deep-water mass reorganisations by 1-2 thousand-years (kyr) during Pleistocene glaciations.With the aid of new iceberg trajectory model experiments, we demonstrate such a shift in iceberg trajectories during glacial periods can result in a considerable redistribution of freshwater in the Southern Ocean.We suggest that this, in concert with increased sea-ice cover, allowed positive buoyancy anomalies to 'escape' into the upper limb of the AMOC, providing a teleconnection between surface Southern Ocean conditions and the formation of NADW.The magnitude and pacing of this mechanism evolved substantially across the Mid-Pleistocene Transition (MPT), and the coeval increase in magnitude of the 'southern escape' and deep circulation perturbations implicate this mechanism as a key feedback in the transition to the 100-kyr world.
In the modern ocean, the AMOC is characterised by the deep, southward spread of NADW towards the Southern Ocean (SO), balanced by the northward return of surface, mode, intermediate, and bottom waters 3,4 .Whilst the widely documented shoaling of NADW during Pleistocene glacial periods has previously been explained through changes in North Atlantic processes 5 , this paradigm has been challenged by studies invoking 'upstream' disruptions (e.g.variable Agulhas Leakage 6,7 ), modifying the shallow return of waters to the North Atlantic required to complete this 'Upper Cell' of the overturning circulation.Furthermore, the southward and northward components of the Upper Cell are connected via wind-and buoyancy-related processes in the SO, with much of the outcropping NADW first taking an indirect route, via the so-called 'Lower Cell' of overturning circulation, before joining the AMOC return limb via the upper cell 4 (Extended Data Figure 1).SO conditions have therefore been increasingly invoked by modelling studies [7][8][9] as key in setting NADW dynamics, however paleoceanographic evidence supporting such a causal link remains scarce 10 .
Here we show a tight surface-deep coupling between the presence of far-travelled Antarctic icebergs and deep-water mass structure at the northern limit of the modern Subantarctic Zone (SAZ) during the past 1.65 Million Years (Ma).We examine new and previously published records of ice-rafted debris mass accumulation rate 11,12 (IRD MAR ), and carbon and oxygen isotope records 13 from benthic foraminifera (δ 13 C benthic , δ 18 O benthic ) from a new composite record of sediment cores recovered from the Agulhas Plateau in the south-west Indian Ocean.Spanning the past 1.65

Iceberg Trajectories and Deep Hydrography
The Agulhas Plateau is situated at the southern boundary of the Indian-Atlantic Ocean Gateway (Figure 1a); the upper water column is dominated by eastward flowing Indian Ocean surface waters retroflected as the Agulhas Return Current (i.e.not 'leaked' into the South Atlantic).To the south, the relatively cold and fresh waters of the SAZ meet their northern limit, forming steep hydrographic gradients at the subtropical frontal zone (STFZ; Figure 1).Iceberg presence (and subsequently IRD deposition) at this location can therefore be influenced by regional hydrography (i.e shifts in the STFZ) as well as the export, survivability, and transport of icebergs into and across the SO.To aid in our interpretation of the IRD MAR proxy, we perform iceberg trajectory modelling experiments using Pyberg, an offline implementation of the Finite Element Sea Ice-Ocean Model (FESOM-IB) 14 iceberg drift and decay module.We force Pyberg with the Community Earth System Models (COSMOS) 15 run to quasi-equilibrium under Pre-Industrial (PI) and Last Glacial Maximum (LGM) climate conditions (see online methods).In each experiment, 111 icebergs were seeded from observationally-derived locations in the Weddell Sea 'iceberg alley' (online methods; shaded box in Figure 1a) and are treated as passive particles, moving and decaying according to oceanic (including sea-ice) and atmospheric conditions.In reality the mass balance and discharge of icebergs from the Antarctic Ice Sheet was likely variable in the past, however by initializing these experiments with consistent iceberg seeding, and neglecting a possible Subantarctic source of icebergs 16 (on the petrological grounds outlined in the online methods), we are able to assess the relative redistribution of iceberg melt resulting only from changes in SO conditions.
Much as the surface hydrography at the Agulhas Plateau represents a junction between subantarctic and subtropical regimes, the deep circulation is characterised by competing Northern-and Southernsourced water masses (NSW and SSW, respectively).NSW here refers to NADW, and glacial analogues, which today extends to 3000m water depth as it exits the Atlantic Ocean into the SO 3 .A weakening and shoaling of NSW, as in some glacial reconstructions 2 , would result in increased presence of SSW at the Agulhas Plateau, consisting of some combination of Lower Circumpolar Deep Water and Antarctic Bottom Water, flowing clockwise around the SO, and spreading northwards into the neighbouring subtropical basins 4 .We examine the relative contribution of NSW vs SSW at the Agulhas Plateau through

Starr et al., ACCEPTED MANUSCRIPT in Nature
Post-Print the δ 13 C of Cibicidoides wuellerstorfi, an epifaunal benthic foraminifer shown previously to be closely related to the δ 13 C DIC of ambient seawater 17 , and therefore a widely applied proxy to trace past changes in deep ocean ventilation by ocean circulation 1,2 .
As with the IRD proxy, the δ 13 C benthic record may be influenced by numerous environmental processes.In other words, it is not a simple water mass tracer, and may be affected additionally by productivity signals 18 and non-stationary water mass end-members 19 .After accounting for time-varying end-member signatures using a binary mixing model (Extended Data Figure 4) we find that the %NSW is well correlated to the raw δ 13 C benthic record at the AP comp (r 2 =0.39), supporting our interpretation that relative NSW versus SSW presence at the site is the dominant control on δ 13 C benthic .This inference is in agreement with complementary water mass tracers 20 and basin-wide data compilations and modelling studies 2 .The timing of δ 13 C benthic changes recorded in the AP comp can therefore be instructive in identifying the origin of deep circulation changes with respect to surface SO conditions recorded in the co-registered IRD MAR data.

Surface-Deep Phasing
The glacial AP comp is characterised by generally higher IRD MAR and lower δ 13 C benthic (shaded bars in Figure 2) compared to interglacials, reflecting the increased presence of Antarctic icebergs and the reduced presence of NSW, respectively.A significant negative correlation is present between IRD MAR and δ 13 C benthic (r 2 = -0.39),becoming stronger still when IRD MAR is log-transformed (r 2 = -0.54).The reason for log-transforming the IRD MAR data is twofold: firstly, it has a strongly skewed distribution with a small number of very high values.Secondly, we consider log-IRD MAR a more sensitive indicator of paleoiceberg drift, for example the probability of iceberg presence in the SO decreases by orders of magnitude with distance north of the main iceberg belt in the Weddell Sea sector 21 .To explore the phasing of the close relationship between the co-registered surface and deep AP comp conditions we employ three independent approaches.The strongest correlation between the two records (r 2 = -0.59)occurs when IRD MAR leads δ 13 C benthic by 2±1 kyrs (Extended Data Figure 5d), and spectral coherency analysis reveals a lead for IRD MAR in all three major Milankovitch frequency bands (Figure 3a).Additionally, a numerical algorithm which determines the offsets between peaks in rate of change reveals that maxima in the IRD MAR occur on average 1.65±0.59kyrs before maxima in the rate of δ 13 C benthic decrease (Figure 3b).Together, these approaches independently demonstrate that IRD MAR systematically leads δ 13 C benthic at the AP comp ; the unidirectional result from the 'peak lag algorithm' best allows us to quantify this lead as 1-2 kyrs, providing a valuable constraint on the link between the SO surface variability and deep water-mass geometry during transitions from interglacial to glacial conditions over the past 1.65 Ma.
Understanding this relationship further requires a robust interpretation of both the IRD MAR and δ 13 C benthic proxies.For example, we might interpet higher glacial IRD MAR as indicating a glacial shift in iceberg trajectories.This interpretation is coherent with published IRD records which show that deposition close to Antarctica is generally highest during interglacial periods 22,23 and periods of ice sheet retreat 23,24 (e.g.glacial terminations), whilst at SAZ locations IRD maxima typically occur during glacial periods (this study; 22,25 ).From this we can deduce that glacial conditions in the SO favour the increased transport of icebergs away from Antarctica (Figure 1), as the early decay of icebergs is reduced and the overall survivability and transport enhanced, due to cooler surface temperatures 26 , increased seaice 26,27 , and changes in surface circulation, possibly relating to the westerly-wind belt.Iceberg drift in the ocean is associated with, and can facilitate sea-ice formation 21 , and sea-ice presence in turn can inhibit the wave-driven erosion of icebergs 28 and govern iceberg movement under severe ice conditions 29 .Enhanced sea-ice concentration and extent is a well-documented feature of the glacial SO 26 , and the close correspondence between the AP comp IRD MAR record and an Antarctic ice core proxy for sea ice extent 27 (Figure 2) over the past 8 glacial cycles suggests that such a feedback may be important.

Starr et al., ACCEPTED MANUSCRIPT in Nature
Post-Print Furthermore, our modelling results indicate that LGM conditions are associated with an equatorward shift and lengthening of trajectories relative to PI (Figure 1) and modern observations 30 .Whilst long-term IRD MAR variability in the SAZ is likely influenced by varied calving rates, a dramatically altered distribution of icebergs does not require a substantial change to ice sheet mass balance.We acknowledge, however, that iceberg presence and therefore IRD deposition at the AP comp does require some calving of icebergs to occur, and therefore necessitates a sufficiently large ice sheet.Considering evidence that East Antarctica maintained a marine-terminating ice sheet across the warm interglacial Marine Isotope stage (MIS) 31 31 , we reason that the calving-ice-sheet-condition would have been consistently met during the glacial intervals relevant to our proposed mechanism 32 .

The Southern Escape
In the modern ocean the conveyor of icebergs traversing the SO constitutes a major transport of meteoric freshwater away from the Antarctic Ice Sheet, particularly in the Weddell Sea sector 33 : a sizeable fraction of modern iceberg mass (up to 35% for giant icebergs 14,33 ) can be exported north of 63 • S. By estimating iceberg meltwater inputs from our trajectory experiments (online methods), we find a substantial redistribution of freshwater in the SO, as the latitude at which this meltwater is greatest shifts several degrees to the north under glacial conditions, with meltwater input south of 50 • S diminishing (Figure 1).Furthermore, we find nearly twice as much meltwater in the region between 0-50 • E under LGM conditions compared to PI conditions (22.7% and 40.7%, respectively, of the total initial iceberg mass (online methods)).Whilst this is clearly an idealized representation, with no change in calving-rate imposed between experiments, it demonstrates that changes in SO conditions can dramatically alter the buoyancy budget of the SO.For example in the modern SO, the meltwater associated with iceberg maxima can exceed the local precipitation-evaporation balance 14 .
A connection between SO surface forcing and the deep ocean has been explored through a variety of theoretical lenses.A common theme occurs wherein enhanced sea-ice (formation and extent) results in an expansion of the region of net surface buoyancy loss 8,9 , due to increased brine rejection close to, and subsequent melt occurring away from, the continent.However, the ability of sea-ice buoyancy forcing alone to explain the full NSW shoaling implied by paleoceanographic evidence has been questioned 34 and such a direct link is challenged by the failure of climate models without additional cooling in the SO 15 to simulate a shoaled glacial NSW.We suggest that the equatorward shift in Antarctic iceberg melt, a feature absent from the aforementioned models, may be a key, previously unconsidered component in connecting the surface SO to overturning circulation.This would allow positive buoyancy anomalies to effectively 'escape' the SO into the return AMOC limb; ultimately influencing deepwater formation in the Upper (NSW) instead of the Lower (SSW) overturning cell.The return limb of AMOC consists of surface to intermediate waters entering the South Atlantic gyre 4 through northward Ekman transport across the SO and entrainment in South Atlantic Subtropical Gyre, or via the shedding of warm and saline mesoscale eddies by the so-called Agulhas Leakage 35 (Figure 1a).The importance of this return limb is illustrated by the role of Agulhas Leakage in the resumption of NADW formation at glacial terminations, when an increase in the inter-ocean exchange of salt transfers positive density anomalies into the Atlantic basin, triggering the resumption of a strong, deep AMOC mode 7,36 .The corollary would suggest that reducing the density export into the South Atlantic, for example through the 'southern escape' of iceberg meltwater, would have the inverse effect: suppressing NADW formation.
Such a mechanism finds support from SO 'hosing' experiments which find that freshwater addition in the SO (in these experiments, south of 60 • S) can affect NADW formation 37,38 due to Ekman spread of surface waters away from the SO.However, these experiments also generate a weaker Lower Cell: the paleoceanographic evidence for which remains equivocal.Whilst sea-ice formation may have stimulated the Lower Cell 39 , grounded ice cover in coastal regions of AABW formation may have an opposing

Starr et al., ACCEPTED MANUSCRIPT in Nature
Post-Print effect 32,40 .At the Agulhas Plateau, negative δ 13 C benthic excursions during glacial periods of the past 200 ka are associated with increased near-bottom flow speeds 41 , indicating a vigorous yet poorly-ventilated glacial SSW (i.e.lower overturning cell) consistent with increased AABW convection 10 fed by surface waters isolated from the atmosphere by enhanced sea-ice 39 .That being said, the 1-2 kyr lag time observed here is difficult to reconcile with a direct SO mechanism and the temporal offset between IRD MAR and δ 13 C benthic may instead favour a teleconnection operating with a threshold behaviour.For example, while the northward export of icebergs to the SAZ would directly decrease the buoyancy flux in AABW formation regions, NADW formation may initially remain strong enough to prevent the northward/upward expansion of SSW.This would continue until the associated 'escape' of buoyancy anomalies into the Upper Cell has gradually weakened NADW formation across some threshold whereby the expansion of SSW into the deep Atlantic is possible.

Orbital Pacing and Sequence of Events
With this mechanism in mind, we next examine the timing of events at the AP comp in the context of global climate.In agreement with evidence for an early climate signal in the SO during interglacialglacial transitions 10,42,43 we find that IRD MAR leads δ 18 O benthic at all three significant orbital frequency bands (Figure 3a).Under the assumption that δ 18 O benthic changes synchronously with global climate, this would imply that IRD MAR peaks lead relative to global climate (or at least Northern Hemisphere ice volume).Moreover, a distinct inter-hemispheric asynchrony is suggested by the phasing of IRD between the SO and the North Atlantic: IRD peaks at the AP comp typically precede IRD peaks at ODP Site 982 44 (Figure 4d).A sequence of events emerges in which at the onset of glacial periods a northward expansion of Antarctic iceberg melt occurs, facilitated by increased sea ice extent and surface cooling in the SO 10 , preceding the descent into full global glacial conditions as indicated by maxima in the global δ 18 O benthic record.The external forcing of this sequence of events may therefore be identified by the external forcing of SST and sea-ice in the SO; a role for which obliquity has been invoked 45 .Indeed, we observe cooler surface conditions in our low obliquity model experiment (LGM27ka), compared to LGM (Extended Data 6) and find that icebergs travel further and decay more slowly between 'iceberg alley' and the Atlantic-Indian SAZ.This is consistent with existing SST reconstructions from the Agulhas Plateau 46 showing colder SST during lower obliquity (Figure 2b).
The pacing of IRD MAR peaks by obliquity is apparent during cycles before ∼1.2 Ma (MIS 36) and after ∼0.4 Ma (MIS 11; Figure 4b), however is somewhat obscured inbetween during the Mid-Pleistocene Transition (MPT).During the MPT, the dominant cyclicity in IRD MAR shifts from ∼ 40-kyr −1 to ∼ 100kyr −1 , possibly reflecting a transition to a nonlinear response to orbital forcing associated with larger ice-sheets 47 .We note that the surface-deep lead (IRD MAR versus δ 13 C benthic ) is relatively constant at 1-2 kyr for the entire 1.65 Ma record, however the lead of SO processes over maxima in δ 18 O benthic (and hence glacial conditions) is variable, increasing between ∼1.2 Ma and ∼0.4 Ma.For example, during MIS 24-26, >80% of the total IRD MAR accumulated within the first 50% of the glacial cycle (Figure 4a).This is compared to ∼60% accumulation by halfway through the pre-and post-MPT cycles.An increase in this lead of SO processes over global ice volume may have resulted from growing larger ice sheets across the MPT (which would take longer to develop), with the return to a shorter lead time in the Late Pleistocene implying that the rate of ice sheet growth had increased to accommodate this larger size.Interestingly, recent evidence from the North Atlantic suggests that the rate of northern ice sheet growth may have accelerated across the MPT, reaching a maximum after ∼0.6 Ma, as a consequence of increasing Atlantic Inflow into the Nordic Seas since ∼1.2 Ma (SB, GK, XZ, et al., submitted).Furthermore, glacial IRD MAR maxima at the AP comp increase across the MPT (Figure 4d), beginning at MIS 36 (∼1.2 Ma) and likely reflecting enhanced export of icebergs across the SAZ, coincident with surface cooling and freshening in the Atlantic SO 48 .

Starr et al., ACCEPTED MANUSCRIPT in Nature
Post-Print While this long-term trend may be driven by increasingly iceberg-favourable (cooler) conditions in the region 48 , we cannot rule out a role for some change in the discharge of Antarctic icebergs amid a background of global cooling and possible transitions in the character of Antarctic ice sheets 49 .Regardless, the southern escape of freshwater (IRD MAR ) into the SAZ appears to scale with the magnitude of glacial AMOC perturbations as indicated by δ 13 C benthic minima; notably, the IRD MAR maximum (and δ 13 C benthic minimum) at ∼0.9 Ma coincides with the implied collapse in AMOC as recorded in the nearby deep Cape Basin 50 .15% concentration contour for September) is shown for PI and LGM.Surface ocean circulation is modified from 3 ; AR = Agulhas Retroflection, AL = Agulhas Leakage, ARC = Agulhas Return Current.The shaded grey band represents the Subtropical Frontal Zone (STFZ) 51 .b. 0-50 • E zonal mean of mean SST under PI (pink) and LGM (blue) conditions.SST proxy data are from MARGO 26 (blue symbols; triangle = diatom, circle = radiolarian, square = alkenone).Model SST data are from COSMOS experiments (see online methods).The estimated meltwater distribution is given as total meltwater input in the 0-50 • E zone as simulated by Pyberg (online methods), with values for the PI (pink) and LGM (blue) being the percentage of the total initial iceberg mass (added as icebergs in the yellow 'seeding region') which melts in the region.   O benthic time series are divided into glacial intervals, then the cummulative integral is taken within each interval and normalized to 100% in time and value.Glacial Progress % refers to the relative time between the onset and termination of each glacial interval.Accumulation refers to the normalized cumulative integral, for example 50% indicates that half of the total IRD MAR deposited in that cycle has been deposited; if this falls to the right of the δ 18 O benthic curve, it suggests that IRD MAR is accumulating relatively early within the glacial (see online methods).The three panels show the average and 95% bootstrap confidence intervals (s.e.m) for three periods: Marine Isotope Stage 1 (MIS 1) to MIS 10 (left), MIS 11 to MIS 36 (middle), and MIS 37 to MIS 59 (right).b.IRD MAR on a log-scale (black) with orbital eccentricity (red, dashed line) and orbital obliquity (blue, dashed line) from 52 .Shaded grey bars (top) highlight glacial MISs and the blue gradient bar (second from top) indicate periods with dominant obliquity pacing.c.Wavelet analysis (using a Mortlet wave with zero-padding) of IRD MAR , revealing the temporal evolution of relative power in the eccentricity (red-orange lines) and obliquity (blue lines) frequency bands.d.Inter-hemispheric comparison between the AP comp IRD MAR record (green; top) and a 1 Ma record of IRD (% IRD relative to planktonic foraminifera) from the North Atlantic site ODP 982 44 (bottom; grey; note reversed y-axis).

AP comp
The Agulhas Plateau composite (AP comp ) is a stratigraphic framework consisting of proximal sediment core sites MD02-2588 (41 We utilise existing records from the upper 10.27m of MD02-2588 (δ 18 O benthic and δ 13 C benthic from 13 , IRD counts from 11 ) and correlated the lowest deglaciation in MD02-2588 to the respective deglaciation in Site U1475 using δ 18 O benthic and δ 13 C benthic from both cores.This tie point is at 10.27m in MD02-2588 and 8.16m on the U1475 Splice Extended Data Figure 2a).The chronology is based on graphically aligning δ 18 O benthic to the 'Prob-stack' (a probabilistic stack of δ 18 O benthic using a profile hidden Markov model; 53 ).The Prob-stack includes 180 globally distributed δ 18 O benthic records and provides an update of the LR04 stack and timescale 54

Ice-Rafted Debris
IRD MAR was determined by counting detrital mineral grains in the >150µm sediment fraction (subsampled using a micropaleontology splitter to yield 500-1000 entities) before normalization to sample weight (IRD concentration in #g) and then multiplication by apparent bulk mass accumulation rate (Extended Data Figure 2c), derived from estimates of dry bulk density from 57 and linear sedimentation rate.The temporal evolution, spectral characteristics, and timing of peaks in the IRD record are largely unaffected by conversion from concentration to Mass Accumulation Rate (MAR).This methodology is equivalent to that employed for MD02-2588 (the measurements for which were previously published by 11,12 n = 522), however as no dry bulk density (DBD) data is available for MD02-2588, we estimate it by using a 2nd order polynomial fitted to DBD in the overlapping interval of Site U1475 (r 2 = 0.57): No quantitative distinction is presented between mineralogies of phenocrysts present, however grains that were clearly of volcanic (e.g.tephra) or authigenic (e.g.pyrite) were excluded from the IRD counts.Energy Dispersive X-Ray Spectrometry (EDS) point analysis was used to identify the mineralogy of several grains from a selection of samples (total grains = 31).Of the 31 measured, 23 were quartz (mostly 'clean' quartz with some Fe-K and Fe members), 7 were orthoclase (K-feldspar), and 1 was garnet (Almandine member) (Extended Data Figure 3).The origin of ice-rafted minerals in the Atlantic SO has been the topic of some discussion, as the presence of sea-ice rafted volcanic and mafic material in sediment cores close to volcanic Subantarctic islands can complicate the interpretation of IRD records 58,59 .However, the observed AP comp assemblage of predominantly quartz with some K-feldspar and garnet is distinct from the plagioclase-dominated clear mineral assemblages found close to Bouvet Island 25,59 .This mineralogy implies a continental Antarctic origin 60 : the high quartz proportion may suggest a Ronne Ice Shelf,

Starr et al., ACCEPTED MANUSCRIPT in Nature
Post-Print Filchner Ice Shelf, or Antarctic Peninsula origin 60 whilst the presence of garnet indicates an East Antarctic contribution 25,60 .Micro-textures such as striations and step-like fractures identified on selected quartz grains further support a glacial origin 61 .

Climate Model Experiments
We use a comprehensive fully coupled atmosphere-ocean general circulation model (AOGCM), COS-MOS (ECHAM5-JSBACH-MPI-OM) in this study.The atmospheric model ECHAM5 62 complemented by the land surface component JSBACH 63 , is used at T31 resolution (3.75 • ), with 19 vertical layers.The ocean model MPI-OM 64 , including sea-ice dynamics that is formulated using viscous-plastic rheology 65 , has a resolution of GR30 (3 • ×1.8 • ) in the horizontal, with 40 uneven vertical layers.The climate model has already been used to investigate a range of paleoclimate phenomena 15,66,67 .This indicates that it can capture key features of warm and cold climate and is thus a very suitable climate model for this study, especially glacial SO conditions and AABW formation 15 .In this study, to investigate responses of Antarctic iceberg transport to different climate backgrounds of the Southern Ocean, we utilized published equilibrium experiments under Pre-industrial (PI) and Last Glacial Maximum (LGM) boundary conditions 15 .
To further assess roles of a lower-than-LGM obliquity (i.e.additional cooling on the glacial Southern Ocean), we further conducted one LGM sensitivity experiment with 27ka obliquity (22.25 • ) that is lower than the LGM value (22.949 • ).The three experiments are integrated to quasi-equilibrium and the average of the last 100 years is calculated to represent the corresponding climatology and used for modelling iceberg trajectory.The simulated climatological SO conditions of these runs can be found in Extended Data Figure 6.

Iceberg Trajectory Model
We use Pyberg (https://github.com/trackow/pyberg),an offline Python implementation of the FESOM-IB iceberg drift and decay module 14 .The model simulates iceberg trajectories and alongtrack rates of melting using established iceberg physics [68][69][70] .Pyberg reads monthly forcing data from different climates as simulated by the COSMOS climate model (10m winds and ocean currents for atmospheric and oceanic drags, sea surface height for the surface slope term, sea surface temperature for melting parameterizations).Monthly sea ice fields are also read and allow to lock icebergs into the sea ice when ice concentrations exceed 86% 29,70 .In that case, icebergs are advected along with the sea ice velocity field.Erosion of icebergs by surface waves is damped at high ice concentrations 69 .Pyberg uses the constant density roll-over criterion (Eq. 3 in 71 ), which can be viewed as a more physical version of the widely adopted original ref. 72formulation.Iceberg fracture is not included as parametrization of this process for climate models is still an active area of research 73,74 .
Icebergs are initialized between 63 • W-50 • W in the iceberg alley of the Weddell Sea from an observationallyderived modern distribution of near-coastal iceberg positions and horizontal sizes (from 75 ).The initial iceberg height is set to 250 m and the density to 850 kg/m3 33,75,76 .The total initial iceberg mass in the iceberg alley is thus 49.59 Gt.For consistency, and to start the icebergs within the ocean model domain, the initial positions are shifted 10 • to the East for all experiments to account for the larger spatial extent of the Antarctic Ice Sheet during the last glacial maximum.Along-track melt-rates are gridded and summed on a 1 • x1 • regular grid to produce the meltwater distribution.The modelled PI trajectories are consistent with modern observational datasets 29,30 (Extended Data 7).The LGM results are broadly consistent with the previously published results from the Fine ResolUtion Greenland And Labrador (FRUGAL) intermediate complexity climate model iceberg module 16 , with an equatorward shift and lengthening of trajectories in the Atlantic-Indian Southern Ocean Sector, however discrepancies are apparent such as the more extreme equatorward spread of LGM trajectories in the FRUGAL model compared to Pyberg.

Starr et al., ACCEPTED MANUSCRIPT in Nature
Post-Print This is the result of different iceberg seeding configurations (i.e.we introduce icebergs only in the 'iceberg alley' region and neglect possible Subantarctic sources) and variations in the forcing provided by the underlying climate models.We note that our interpretations based on the Pyberg results are entirely compatible with the results from FRUGAL, and the results of the latter suggest that our estimated meltwater redistribution is likely conservative.Furthermore, the Pyberg LGM results are supported by available ice-rafted debris records from the SO [22][23][24][25]77,78 , which typically show a latitudinal divide between higher (lower) IRD accumulation at close to (away from) Antarctica during glacial compared to interglacial intervals.

Water-Mass Mixing Model
To identify the water mass signal on the AP comp δ 13 C benthic , we apply a simple binary mixing model (following 79 ) to estimate the % of NSW from stacked δ 13 C benthic records representing NSW and SSW over the past 1.5 Ma: Where δ 13 C N SW is a shallow North Atlantic stack (consisting of ODP Sites 980, 982, 983; 80 and references therein) and δ 13 C SSW is a deep South Atlantic stack (consisting of ODP Sites 1089 and 1090; 81 ).Stacks were created by converting all data onto the LR04 54 timescale, smoothing with a 5-kyr gaussian filter, and averaging.We apply this technique using both a SSW stack and a PDW stack as the second end-member, as in reality the "Southern" deep water-mass present in the glacial Atlantic is likely some combination of the two.For example, whilst several studies employ a PDW end-member 79,82,83 , Raymo et al. 82 note that this is due to the absence of suitable Southern Ocean end-member records.Indeed, Venz and Hodell 84 later apply the same approach with newly available deep Southern Ocean δ 13 C benthic records.However, recent evidence indicates that PDW was present in the glacial Atlantic Ocean in larger proportions than previously thought 85 and so a PDW end-member may be preferable after all.

Times Series Analysis
The lagged cross-correlation between IRD MAR against δ 18 O benthic and δ 13 C benthic is determined using the Gaussian-kernel-based cross-correlation algorithm of 86 (Extended Data figure 5b).This allows crosscorrelation analysis to be performed on the original, irregularly-spaced data, in order to avoid the possibility of spurious phasing introduced by linear interpolation.
The 'peak-lag' algorithm we design to iteratively measure the offset between peaks in the rate of change of IRD MAR and δ 13 C benthic is modified from the approach of 87 and analogous to 'Event Synchronisation' techniques.It smooths the time series, finds the first difference, and then measures the offset in local maxima of both signals.We tested the sensitivity of the algorithm to the smoothing filter selection, finding that the average lag computed for a range of filter designs (moving-average and Savitsky-Golay filters; orders between 1 and 21) give similar results within 1σ of each other.We additionally test the confidence that the result from our final chosen parameters (7-point moving-average) is significantly non-zero by performing 1000 'Monte Carlo' simulations with zero-lag red noise surrogates (generated from an autoregressive model of the real data), finding that the lag identified in our data (IRD MAR vs δ 13 C benthic ) is significantly different from the surrogate series (p<10 -10 ); this additionally yields a 99% red noise confidence level of ±0.4 kyrs for the lag calculated for our data.We test the utility of the algorithm to detect lead-lag relationships in a series of surrogate time series with known lags applied by constructing synthetic time series from 3 sine waves (100, 41, 23-kyr periods), with higher-frequency oscillation (7-kyr period) and white-noise components added to emulate sub-orbital features.Known lags are imposed between pairs of the surrogate series.The results show that the algorithm performs well and is capable of detecting a lag of down to 0.4 kyrs to a confidence of p<10 -10 .The algorithm shows an

Starr et al., ACCEPTED MANUSCRIPT in Nature
Post-Print inherent drift at higher lags (Extended Data figure 5a and b), meaning it underestimates lag times greater than 4 kyrs, however the 1.5 kyrs lag observed in our data should be short enough to be unaffected by this drift.
To divide the time series into glacial cycles (as in figure 4) we define intervals between peak interglacial conditions (identified as local minima in δ 18 O benthic separated by a minimum distance of 25-kyrs) and the next time δ 18 O benthic crosses an 'interglacial threshold' of 3.32 (C.wuellerstorfi scale, equivalent to 3.75 on the LR04-scale), which is lower than adopted by 88 , but ensures that we do not define the interstadial MIS 5c as an interglacial.We manually add a peak interglacial marker at 506 ka (based on Northern Hemisphere perihelion) to include the 'cool' interglacial MIS 13a.Within each interval we calculate the cumulative integral of IRD MAR and δ 18 O benthic , before normalizing the values to percentages of total 'accumulation' for that interval (Extended Data figure 5d).
Finally, for the phase wheels in figure 3, we use the Blackman-Tukey method with a Bartlett window (using the Analyseries software 89 ) after linearly interpolating (dt=1.5 kyrs) and standardizing (mean = 0, standard deviation = 1) the data.

Orbital Forcing and the Mid-Pleistocene Transition
We estimate the time-frequency spectral power in the AP comp IRD MAR time series by applying a continuous wavelet transform (using zero-padding and a Mortlet mother wave) 90 .To visualize the results, we extract the power (in units of normalized variance) at the 'Fourier periods' between 39 to 43 kyr −1 to represent frequencies associated with obliquity, and 92 to 108 kyr −1 to represent frequencies associated with eccentricity (figure 4c).We emphasise that this gives an estimate of relative not absolute spectral power in these bands.Between 1600 and ∼900 ka, obliquity-band power is dominant, decreasing after 900 ka in favour of the ∼100 kyr −1 power.This transition then reverses after ∼400 ka, when obliquity-band power again dominates.It is important to note that spectral power at frequencies of ∼40 kyr −1 and ∼100 kyr −1 is not evidence for a direct forcing by the obliquity and eccentricity of Earth's orbit (with dominant periodicities at these values, respectively).Furthermore, it is instructive to also evaluate the timing of cycles in the IRD MAR data relative to cycles in orbital parameters.In figure 4b-c, we show that while the spectral power is dominated by ∼40 kyr −1 frequencies (pre-and post-MPT) there is a close coupling between high IRD MAR and low obliquity.A mechanistic relationship here is supported by our climate model results (Extended Data 6) which indicate that the low obliquity 27ka experiment yields cooler SO SST and increased sea-ice extent relative to LGM.Indeed, a relationship between obliquity and SO seaice 15,91 , as well as stronger westerlies (related to the increased latitudinal insolation gradient 45,92 ) has been previously demonstrated; both would act to increase iceberg survivability and transport into and across the SAZ (Extended Data 7).However, the transient emergence of ∼100 kyr −1 power between ∼1.2 Ma and ∼0.4 Ma coincides with the concealment of this apparent obliquity pacing.Recent modelling results have shown that SO sea ice is more sensitive to obliquity than its Northern Hemisphere counterpart 92 , and one explanation for the obscuration could be a shift to more complex, non-linear pacing of SO climate as NH ice sheets expand and glacial-interglacial variations in atmospheric CO 2 increase 47 .The post-MBE return to obliquity-pacing might then indicate a shift in the nature of NH ice sheet growth, for example higher rate of growth relating to an increase in the Atlantic Inflow(SB, GK, XZ, et al., submitted).

Extended Data Figure Captions
Extended Data Figure 1: Schematic representation of the 'Southern Escape' mechanism described in this study.a. Shows an idealized representation of the Southern Ocean during Interglacial conditions.The front panel gives meridional-averaged overturning circulation (following 93 ), with circulation divided into an Upper Cell and Lower Cell.Arrow width represents the relative strength of each circulation cell (wider being stronger).The top panel represents the ocean surface, with icebergs calving from the Antarctic Ice Sheet (AIS) and following a northward then eastward trajectory.The blue shading represents iceberg meltwater being entrained into the Lower Cell as icebergs melt south of the main Antarctic Circumpolar Current (ACC) belt.Orange wavy arrows represent brine rejection from sea-ice formation.The peach-coloured band represents the Subtropical Frontal Zone (STFZ) delineating the Subtropical regime to the north and the (sub)Antarctic regime to the South.b.Same as a., but for conditions during Glacial Inception; i.e. the transition from interglacial to glacial conditions.The major change from a. is the northward displacement of iceberg trajectories, resulting in much of the meltwater spreading northwards and being entrained in the Upper Cell, instead of Lower Cell.Secondly, colder surface conditions facilitate extended sea-ice cover, and subsequent increase in brine rejection.The Lower Cell thus experiences less positive buoyancy forcing from the combination of less iceberg meltwater and more brine rejection, becomes stronger.The Southern Escape of meltwater into the Upper Cell has not yet occurred in high enough amounts to perturb the Upper Cell and hence the geometries of the overturning cells are unchanged from a. c.Same as a. and b., but for full glacial conditions (occurring after b.).Here, the Southern Escape of meltwater has successfully perturbed the Upper Cell, which is now weaker and has contracted upwards.The now-dominant Lower Cell has thus increased in volume.
Extended Data Figure 2: AP comp Age model and composite record construction.a.The appending of U1475 to the bottom of MD02-2588 using δ 18 O benthic with the new AP comp composite depth scale on the right-hand axis.b.Age-depth tie points (circles) with the median and 95% confidence bounds for the deterministic age model.Depth is given on the AP comp scale.c.The δ 18 O benthic from this study (purple) and the tuning target, the global delo stack 53 (light grey), and the implied linear sedimentation rates (black; dashed line) and calculated mass accumulation rate (green; solid line) from the final age model in the lower subplot.Post-Print and a PDW (grey) or SSW (purple) end-members.b.Scatter plots of AP comp δ 13 C benthic vs. %NSW calculated with an SSW (left) and PDW (middle) end-member, and vs. the ǫNd isotope record of 20 (right).Pearson's correlation coefficient r 2 is given.c.Shows Pre-Industrial δ 13 C of Dissolved Organic Carbon (from 94 along a North-South (blue) and then East-West (red) transect, with the AP comp position shown as a white circle.

Extended Data
Extended Data Figure 5: Time Series Analysis Algorithm Tests.a. and b.Testing the peak-lag algorithm to detect relationships in a series of surrogate time-series with known lags applied.The 'actual lag' axes represent the known lag-time imposed between pairs of the surrogate series and the 'calculated lag' axes show the lag estimated by the algorithm.A perfect performance would manifest as a 1:1 straight line through the scatter points.The violin plots in b. show the mean, median and kernel probability density estimates of calculated lags from 10 iterations of the test.c. Results of the 'glacial accumulation' algorithm with each shaded curve representing IRD MAR (green) and δ 18 O benthic (purple) integrated and normalized to 100% within each glacial cycle.Above, the δ 18 O benthic record is shown (purple, solid) with green triangles denoting the peak interglacials identified and the dashed black line showing the δ 18 O benthic threshold above which the transitions from glacial to interglacial conditions are defined.d.Gaussian-kernel-based cross-correlation (gXCF, following 86 ) function for AP comp IRD MAR vs. δ 18 O benthic and δ 13 C benthic .The horizontal lines show the 95% Monte Carlo confidence levels for significant crosscorrelation.LGM, c. 27ka).Solid/dashed contours represent 90%/15% sea ice concentration.d-f.As for a-c.but for climatological annual mean Sea Surface Temperature (SST).g.Anomaly of mean annual SSH between 27ka and LGM experiments.Solid/dashed contours give 90%/15% sea ice concentration (black and white lines represent LGM and 27ka, respectively).h.Same as g. but for SST.Extended Data 7: Pyberg model experiment results.Pyberg results.a., c., and e. Spatial distribution of meltwater input estimated for 1 x 1 • cells by Pyberg when forced by COSMOS outputs for Pre-Industrial (a; PI), Last Glacial Maximum (c; LGM), and Last Glacial Maximum 27ka (e; LGM27ka) conditions.In purple are observed modern iceberg trajectories from the BYU Giant Iceberg database 95 (QSCAT); the modelled PI trajectories appear substantially longer than the observed likely because icebergs are tracked in Pyberg even after then become too small to be identified and hence tracked by modern observational techniques.b., d., and f.Zonally averaged meltwater estimates for PI (b), LGM (d), and LGM27ka (f) experiments.The average is taken for each latitude between 0 and 50 • E; in other words, this shows the latitudinal distribution of meltwater across the Indian-Atlantic Ocean Gateway.Extended Data Figure 3:

Figure 1 :
Figure 1: Modelled iceberg trajectories and meridionally-averaged SST and meltwater distribution for the PI and LGM Atlantic Southern Ocean.a. Iceberg trajectories simulated by Pyberg for Pre-Industrial (PI) and Last Glacial Maximum (LGM) conditions (see online methods).Winter sea ice extent from 26 (15% concentration contour for September) is shown for PI and LGM.Surface ocean circulation is modified from 3 ; AR = Agulhas Retroflection, AL = Agulhas Leakage, ARC = Agulhas Return Current.The shaded

Figure 2 :
Figure 2: Paleoceanographic proxy records from Mid-to Late-Pleistocene AP comp .a. Time-series of δ 13 C benthic (top; pink, note reversed scale), IRD MAR (middle; teal), and δ 18 O benthic (lower; grey, note reversed scale) from the AP comp for the past 1.65 Ma, with glacial Marine Isotope Stages shaded grey.The interval 0-800 ka is expanded in panel b. which additionally includes the Na sea-salt accumulation from the EPICA Dome C Antarctic ice core27 , a proxy for SO sea-ice extent (top, purple).The AP comp U k' 37 SST record from46 is also shown (second from top, orange).Note the log-scale for IRD MAR in b.

Figure 3 :
Figure 3: IRD MAR vs δ 13 C benthic lead-lag relationship.a. Blackman-Tukey coherency phase of IRD MAR against δ13 C benthic (red; triangular-head arrows) and δ 18 O benthic (blue; diamond-head arrows) in the 100 kyr −1 , 40 kyr −1 , and 23 kyr −1 frequency bands.These phase wheels are cropped at 0-90 • , with clockwise rotation indicating a larger lead for IRD MAR (given in kyrs on the curved axes).The arrow length represents coherency and the shaded areas give the 95% uncertainty in phase.IRD MAR leads δ 13 C benthic and δ18 O benthic in phase for all 3 frequencies.b.Summary of lead times for IRD MAR over δ13 C benthic determined by three independent approaches, with estimates of 95% uncertainty given by the back bars (see online methods).

Figure 4 :
Figure 4: The evolution of the AP MAR IRD MAR record in the time and frequency domain orbital forcing, global climate, and interhemispheric phasing.a. Results from the 'glacial accumulation' analysis: IRD MAR and δ18 O benthic time series are divided into glacial intervals, then the cummulative integral is taken within each interval and normalized to 100% in time and value.Glacial Progress % refers to the relative time between the onset and termination of each glacial interval.Accumulation refers to the normalized cumulative integral, for example 50% indicates that half of the total IRD MAR deposited in that cycle has been deposited; if this falls to the right of the δ18 O benthic curve, it suggests that IRD MAR is accumulating relatively early within the glacial (see online methods).The three panels show the average and 95% bootstrap confidence intervals (s.e.m) for three periods: Marine Isotope Stage 1 (MIS 1) to MIS 10 (left), MIS 11 to MIS 36 (middle), and MIS 37 to MIS 59 (right).b.IRD MAR on a log-scale (black) with orbital eccentricity (red, dashed line) and orbital obliquity (blue, dashed line) from52 .Shaded grey bars (top) highlight glacial MISs and the blue gradient bar (second from top) indicate periods with dominant obliquity pacing.c.Wavelet analysis (using a Mortlet wave with zero-padding) of IRD MAR , revealing the temporal evolution of relative power in the eccentricity (red-orange lines) and obliquity (blue lines) frequency bands.d.Inter-hemispheric comparison between the AP comp IRD MAR record (green; top) and a 1 Ma record of IRD (% IRD relative to planktonic foraminifera) from the North Atlantic site ODP 98244 (bottom; grey; note reversed y-axis).

Figure 3 :Extended Data Figure 4 :
SEM Imaging and Mineralogy of IRD Grains from the AP comp .a.An SEM image of a quartz grain from Site U1475 with an enlarged view of surface microtextures in the inset panel.b-d.SEM images of IRD grains with examples of relative peak intensities for elements derived from EDS point analyses.b shows a spectrum typical of quartz, c of garnet (almandine member), and d of Kfeldspar (orthoclase).Scale bars are given in white on all SEM images.Comparison of δ 13 C benthic and water mass 'end-members'.a. Authigenic ǫNd isotope record from the deep equatorial Atlantic (ODP Site 154-929 from 20 ; top, green), Smoothed δ 13 C benthic stacks for Northern-Sourced Water (NSW), Pacific Deep Water (PDW) and Southern-Sourced Water (SSW) end-members (see Methods for stack construction and constituent core sites) and the AP comp δ 13 C benthic record.Selected Marine Isotope Stages (MIS) are shown as grey vertical shading.The bottom time series in a. shows % NSW at the AP comp calculated using a binary mixing model with NSW Starr et al., ACCEPTED MANUSCRIPT in Nature
56The Prob-stack and LR04 are similar over the MPT interval studied here.From selected age-depth markers (n = 45; consisting 12 radiocarbon dates and 33 δ 18 O benthic ties points; Extended Data Figure2b), the final age-depth model was constructed using the deterministic age modelling routine Undatable55with 105 bootstrap simulations.The resultant sedimentation rates range from 1.2 cm/kyr to 7.6 cm/kyr, averaging 3.1 cm/kyr.The Brunhes/Matuyama magnetic reversal (0.781 Ma) and the Jaramillo Subchron (1.072 -0.988 Ma) are identified in the U1475 down-core magnetic inclination56, confirming the robustness of the final age-depth model.The records presented have a temporal resolution ranging from ∼0.3 kyrs in the Late Pleistocene to ∼1.5 kyrs in the early-mid Pleis- tocene.Stable isotope measurements were performed on 2-3 well-preserved Cibicidoides wuellerstorfi specimens (250-315µm size fraction) and made at Cardiff University using a Thermo Finnigan MAT 253 mass spectrometer linked online to a Carbo Kiel carbonate preparation device with long-term precision of ±0.05 for δ 18 O and±0.021 for δ 13 C (±1σ).Results are calibrated to an internal laboratory standard (BCT63) and presented relative to the Vienna Pee Dee Belemnite scale.