Millennial‐Scale Climate Oscillations Triggered by Deglacial Meltwater Discharge in Last Glacial Maximum Simulations

Our limited understanding of millennial‐scale variability in the context of the last glacial period can be explained by the lack of a reliable modeling framework to study abrupt climate changes under realistic glacial backgrounds. In this article, we describe a new set of long‐run Last Glacial Maximum experiments where such climate shifts were triggered by different snapshots of ice‐sheet meltwater derived from the early stages of the last deglaciation. Depending on the location and the magnitude of the forcing, we observe three distinct dynamical regimes and highlight a subtle window of opportunity where the climate can sustain oscillations between cold and warm modes. We identify the Eurasian Arctic and Nordic Seas regions as being most sensitive to meltwater discharge in the context of switching to a cold mode, compared to freshwater fluxes from the Laurentide ice sheets. These cold climates follow a consistent pattern in temperature, sea ice, and convection, and are largely independent from freshwater release as a result of effective AMOC collapse. Warm modes, on the other hand, show more complexity in their response to the regional pattern of the meltwater input, and within them, we observe significant differences linked to the reorganization of deep water formation sites and the subpolar gyre. Broadly, the main characteristics of the oscillations, obtained under full‐glacial conditions with ice‐sheet reconstruction derived meltwater patterns, share similar characteristics with δ18O records of the last glacial period, although our experiment design prevents detailed conclusions from being drawn on whether these represent actual Dansgaard‐Oeschger events.

During decades of study, numerous hypotheses have been put forward to understand the underlying mechanisms behind D-O events (a comprehensive list can be found in Li and Born [2019]), and more generally, millennial scale variability. Despite this effort, a consistent and comprehensive theory is yet to be firmly established. Nonetheless, at the crossroads of all theories lies the crucial role of the Atlantic Meridional Overturning Circulation (AMOC) (Burckel et al., 2015;Henry et al., 2016;Rahmstorf, 2002). A modification of the thermohaline circulation affects heat and salt redistribution between the tropics and the poles, and consequently has a global-scale impact on the climate (Clark et al., 2002;Rahmstorf, 2002). There is substantial evidence that the AMOC has existed in other configurations (or "modes") than the one we observe at present times (e.g., Böhm et al., 2015), and that AMOC may thus have the capacity to exist in multiple stable states, as predicted theoretically (Stommel, 1961) and supported by early observations (Broecker et al., 1985) and climate models (Manabe & Stouffer, 1988). We believe that uncovering the causes of AMOC modes switches would result in greatly improving our understanding of abrupt climate changes.
The AMOC can be disrupted by freshwater release events in the North Atlantic-Arctic region. They have the power to target vital points of the thermohaline circulation by affecting the ocean density profile at North Atlantic Deep Water (NADW) formation sites (Broecker et al., 1985;Paillard & Labeyriet, 1994;Vidal et al., 1997). In models, freshwater hosing experiments have been widely used to force abrupt climate transitions (e.g., Ganopolski & Rahmstorf, 2001;Kageyama et al., 2010;Manabe & Stouffer, 1997) and observe hysteresis cycles (e.g., Schmittner et al., 2002). They also highlighted the large sensitivity of the climate to the strength and the location of the release, especially in the Greenland, Iceland, and Nordic (GIN) Seas Smith & Gregory, 2009). Consequently, it is valuable to explore the different sources of freshwater that had the potential to lead to millennial-scale variability.
Iceberg surges during Heinrich events (H events ;Heinrich, 1988) recorded by Ice Rafted Debris (IRD) in the North Atlantic (Hemming, 2004), were first candidate to be held responsible for initiating stadial climates. It is now widely accepted that H events are not at the origin of D-O events, they are triggered within stadial states (Barker et al., 2015) and are not recorded at every D-O occurrence (Lynch-Stieglitz, 2017). Instead, we can conceive of them as a likely response to the earlier climate-ocean perturbation or even a positive feedback mechanism for perpetuating/amplifying stadial climates . Meltwater released from the long term decline of ice sheets was another significant source of freshwater during the last glacial period (Gregoire et al., 2012), although only a few studies have investigated the influence of such "background" melt (e.g., Ivanovic et al., 2018;Kapsch et al., 2022;Matero et al., 2017), probably because it requires precise constraints on the ice sheet geometry and history of melt/growth (Bethke et al., 2012). Holding the most complete records of ice sheet evolution, both in terms of spatial and temporal resolution (e.g., Bradwell et al., 2021;Briggs et al., 2014;Dyke, 2004;Hughes et al., 2016), the last deglaciation (and especially its early phase, ∼21-16 ka BP, thousand years before present) offers the perfect setting to assess the ability of the early phase of continental deglaciation (i.e., the long-term background melt from disintegrating ice sheets) to generate millennial-scale variability in glacial conditions. The last deglaciation initiated from the Last Glacial Maximum (LGM; ∼21 ka, Clark et al., 2009), which is the time when the Northern Hemisphere ice sheets shaped during the preceding glacial period reached their maximum extent (Batchelor et al., 2019). The upper cell of the AMOC was likely shallower, but it is not known whether it was stronger or weaker than present day (Gebbie, 2014;Lynch-Stieglitz, 2017;Muglia & Schmittner, 2021). A steady increase in Northern Hemisphere summer insolation (Berger, 1978) triggered the long-term demise of the Laurentide and Eurasian ice sheets, with rising concentrations of atmospheric CO 2 positively reinforcing the deglaciation (Gregoire et al., 2015). However, while Southern Hemisphere temperatures gradually rose (Parrenin et al., 2007), the climate in the North remained cold for several thousand years; a period known as Heinrich Stadial 1 (∼18-15 ka BP) (Denton et al., 2006;Ng et al., 2018;Roche et al., 2011). The most recent of Heinrich events, H1 (Hemming, 2004;Stanford et al., 2011), began some two thousand years after the onset of Henirich Stadial 1 (Hodell et al., 2017;Stern & Lisiecki, 2013). In the years of deglaciation that followed the LGM, several millennial scale events were observed (Weber et al., 2014). Two episodes are particularly relevant to our study: the sudden Bølling Warming (∼14.5-13 ka BP; Severinghaus & Brook, 1999) concurrent with an intensification of the AMOC (Du et al., 2020;Ng et al., 2018), and the ensuing Younger Dryas, when Northern Hemisphere climate abruptly returned to a stadial state with glacial readvance (∼13-12 ka; Liu et al., 2012;Murton et al., 2010). While not formally identified as D-O events, similarities in climate and ocean evolution between Simulating climate oscillations in glacial conditions has proven to be very challenging, and even more so during the LGM. This is because of the strong feedback between the large ice sheets and wind stress, deep water formation and energy balance (Beghin et al., 2015;Oka et al., 2012;Ullman et al., 2014), which act to intensify, or at least stabilize, the AMOC (Klockmann et al., 2016;Oka et al., 2012;Sherriff-Tadano et al., 2018). Most models from both the Paleoclimate Modeling Intercomparison Project Phase 3 (PMIP3; Muglia & Schmittner, 2015) and Phase 4 (PMIP4; Kageyama et al., 2021) tend to simulate a deeper and stronger NADW than inferred from paleo records, which could explain why few modeling studies have observed millennial scale variability in glacial background (e.g., Klockmann et al., 2018). In order to trigger abrupt climate transitions, freshwater hosing experiments have historically needed to overestimate fluxes as reviewed by Kageyama et al. (2010) (e.g., Liu et al., 2009Menviel et al., 2011), and have not succeeded in simulating abrupt changes when using "realistic" fluxes (Bethke et al., 2012;Gregoire et al., 2012;Snoll et al., 2022). Obase and Abe-Ouchi (2019) have arguably come the closest to overcoming this meltwater "paradox" by simulating the Bølling Warming even with some deglacial meltwater forcing. However, even they require a significantly lower than likely freshwater discharge from the deglaciating ice sheets (e.g., Peltier et al., 2015).
The dispute over what could feasibly cause abrupt climate changes not directly driven by freshwater fluxes led the community to start actively searching for oscillating behaviors in their models. At the same time, the criticism that simulations integrated for only a few hundred or a thousand years should not be considered to have a steady-state or "spun-up" ocean circulation began to gain traction (Dentith et al., 2019;Marzocchi & Jansen, 2017), prompting modelers to run long simulations with higher-order climate models-made possible by the increase of computational power-in order to examine long-term drifts. It is therefore probably not a coincidence that more and more coupled Atmosphere-Ocean General Circulation Models (AOGCMs) have reported observing AMOC mode oscillations in recent years (e.g., Brown & Galbraith, 2016;Klockmann et al., 2018;Peltier & Vettoretti, 2014;Sherriff-Tadano & Abe-Ouchi, 2020). They have been achieved under a range of different freshwater hosing scenarios (e.g., Cheng et al., 2011), atmospheric CO 2 concentrations (e.g., Zhang et al., 2017), and ice sheet geometries (e.g., Klockmann et al., 2018), although, to our knowledge, only Peltier and Vettoretti (2014) and Kuniyoshi et al. (2022) managed to obtain AMOC oscillations under glacial conditions.
To sum-up the combined results from these studies, there seems to exist a window of opportunity (Barker & Knorr, 2021;Vettoretti et al., 2022) in each model's inputs (parameter values, boundary conditions, and forcings) and background climates where oscillations can establish and sustain (e.g., Brown & Galbraith, 2016;Klockmann et al., 2018;Peltier & Vettoretti, 2014). The ice sheets' layout in particular has a strong influence on the local and global climate, including on the atmospheric circulation (Löfverström et al., 2014;Roberts et al., 2014;Sherriff-Tadano et al., 2021), the gyres , the energy balance , and freshwater fluxes (Matero et al., 2017). As a result, the new generation of better constrained and more detailed ice sheet reconstructions such as ICE-6G_C (Argus et al., 2014;Peltier et al., 2015) and GLAC-1D (Briggs et al., 2014;Ivanovic et al., 2016;Tarasov & Peltier, 2002;Tarasov et al., 2012) may prove to be decisive in whether or not the "right" conditions for triggering abrupt climate changes are obtained.
In this paper, we present our contribution to this initiative in the form of a new set of LGM simulations forced with deglacial meltwater. Inspired by an initial experiment that showed millennial-scale variability under a transient meltwater forcing, we designed our simulations with fixed meltwater inputs in order to be able to describe the oscillations in detail and evaluate the sensitivity of the oscillatory behavior to meltwater patterns. These inputs were derived from snapshots of the early deglaciation meltwater history (specifically, between 21.5 and 17.8 ka BP) calculated from GLAC-1D ice sheet reconstruction.
Depending on the freshwater pattern, we observe three different dynamical regimes, including regular and self-sustained climate oscillations. The oscillations are characterized by switches between strong, shallow glacial AMOC, and near-or completely collapsed AMOC modes, a Greenland surface cooling/warming of ∼10°C, and a periodicity of about 1,500 ka. The cold/warm modes resemble, but should not be considered too strictly to be stadial/interstadial states, due to the use of a quasi-ideal experimental set-up. The oscillating regime can be sustained for about 10,000 years (the maximum length of the experiments). These are, to our knowledge, the first general circulation model simulations to use ice sheet reconstruction-derived distributions of meltwater to produce strong AMOC oscillations under glacial maximum climate conditions, and to investigate the sensitivity of the oscillations to patterns of meltwater discharged to the ocean.
The non-oscillating regimes inform us of the pre-requisite conditions for passing through the oscillating window. Here, we describe the various simulations in detail, their oscillatory or non-oscillatory climate/ocean states, and the different Earth system components involved in the abrupt events. We conclude with a discussion on the relevance of our simulations in the context of known past abrupt climate changes. The design introduced in this study allows us to undertake a relatively systematic set of sensitivity tests of the impact of realistic freshwater distributions (albeit for unrealistic lengths of time) on ocean circulation, with the multi-millennial integrations enabling us to explore the long-term effect of each pattern.

The Model
The simulations introduced in this article were completed using the BRIDGE (Bristol Research Initiative for the Dynamic Global Environment group) version of the HadCM3 atmosphere-ocean general circulation model (GCM) . This GCM consists of a 19 layers × 2.5° × 3.75° atmosphere model more completely described by Pope et al. (2000), coupled every simulation day with a 20 layers (up to 5,500 m deep) × 1.25° × 1.25° ocean model, described by Gordon et al. (2000) (Bryan & Cox, 1972;Fofonoff, 1985;Fofonoff & Millard, 1983). This version of HadCM3 includes the MOSES 2.1 land model (P. M. Cox et al., 1999), and the TRIFFID dynamic vegetation model (P. Cox, 2001). HadCM3 has been tested in many different scenarios (I.P.C.C., 2014; Reichler & Kim, 2008), and was optimized for running multi-millennial simulations .

Experimental Design
The LGM simulation that makes up the base climate state for all simulations presented here was created following the PMIP4 protocol for 21 ka BP (Kageyama et al., 2017) using the GLAC-1D ice sheet reconstruction (Briggs et al., 2014;Ivanovic et al., 2016;Tarasov & Peltier, 2002;Tarasov et al., 2012); see Section S1 and Figure S1 in Supporting Information S1 (Bereiter et al., 2015;Loulergue et al., 2008;Schilt et al., 2010). This new HadCM3 LGM simulation was initialized from a chain of existing multi-millennial HadCM3 PMIP3 LGM simulations, which were started from multi-millennial continuations of earlier HadCM3 LGM simulations , giving a pre-PMIP4 LGM spin-up of several thousand years. The new PMIP4 GLAC-1D set-up was integrated for 3,500 years, and the end of this final spin-up phase provides the initial condition for all simulations presented here. We continued the LGM simulation for a further 4,000 years in parallel with our other experiments to provide a reference climate state (CTRL) for comparison to the other simulations. There are small, steady drifts in the ocean over the run ( Figure S2 in Supporting Information S1), but the signal of the trends are dwarfed in comparison to the changes of interest described below, and little is gained for this study by extending CTRL further. The starting year of CTRL is defined as year 0.
GLAC-1D has rarely been used for LGM simulations compared to the other ice sheet reconstructions (Kageyama et al., 2021) such as ICE-6G_C (Peltier et al., 2015) and the PMIP3 ice sheet (Abe-Ouchi et al., 2015). It was preferred for this study because compared to the alternative reconstructions, it includes more recent constraints on the Eurasian ice sheets (provided by the DATED-1 project; Hughes et al., 2016), a region that could be crucial for accurately capturing the early deglacial climate history . A transient meltwater history was derived from GLAC-1D's representation of the deglaciation (Section S2 in Supporting Information S1). We decided against using the transient meltwater flux, because the added complexity introduced by the temporal variability and possible ocean "memory" of the preceding [uncertain] meltwater history would have convoluted the physical interpretation of our results. Instead, we examined the triggering of abrupt climate changes using a simpler approach; by selecting six different, fixed-forcing scenarios (our "snapshots"), that allow us to investigate the sensitivity of the glacial ocean and surface climate to early deglacial freshwater inputs (Figure 1a).
The snapshots were identified for their ability to collectively capture a broad range of possible situations that may have led to changes in ocean circulation and surface climate. The six scenarios correspond to different modes of discharge and are named after the period they were extracted from (see Figure 1a); see Figure 1c for the spatial distribution of the fluxes. The 21.5k, 21k, and 20.7k snapshots were chosen for being close to the LGM, sharing a similar distribution, but with different rates of meltwater discharge. The 19.4k snapshot hosts a strong Labrador Sea/North Eastern American coast/Gulf of Mexico discharge (shortened to "North American" discharge hereafter), but has a relatively small meltwater flux to the Arctic. Conversely, 18.2k and 17.8k have high Arctic and low North Atlantic discharge, with 18.2k having the most freshwater entering the Arctic. These six snapshots of the deglacial meltwater history were used as forcing for six new equilibrium-type (i.e., constant-forcing) simulations, started from year 0. The meltwater was discharged into the surface layer of the ocean and kept constant throughout the runs. A justification is given in Section S2 of Supporting Information S1 (Quadfasel et al., 1990;Roche et al., 2007). Table 1 presents a summary of all experiments. The difference between any of them is the prescribed ice sheet meltwater (or absence of it, in CTRL). . This plot incorporates the 200 years smoothing described in Section S2 of Supporting Information S1. Vertical bars represent the time steps chosen for calculating each constant meltwater-forcing snapshot (see Section 2.2 and Table 1). Regions of low discharge are plotted in dotted lines. (b) Map of ice meltwater collection and spreading areas. Each individual box corresponds to a freshwater collection area, redistributed to the corresponding spreading areas (within the same box) indicated by the bold contours. Six main regions were defined for the presented analysis, as labeled on the right (colors). Note that these regions do not correspond to individual regions but rather to clusters of spreading areas. The color coding matches panel ( The idea of having a continuous fixed meltwater discharge for thousands of years is unrealistic by nature. However, in the most extreme scenario (18.2k), the total forcing corresponds to a sea level rise of 102 m in 10,000 years, which for context, is still less than what has been reconstructed for the whole of the last deglaciation (Lambeck et al., 2014). The fixed global flux, between 0.039 and 0.109 Sv, is of the same order of magnitude as the mean global discharge in Ivanovic et al. (2018): ∼0.05 Sv, derived from ICE-6G ice sheet reconstruction (Peltier et al., 2015). It is also lower than the mean estimations of previous Heinrich Stadial 1 transient simulations (e.g., Liu et al., 2009;Menviel et al., 2011). It therefore remains appropriate to consider our results in light of glacial and deglacial variability in order to understand the effect of the forcing, though we are careful to highlight that they are not transient simulations of deglacial meltwater neither do they correspond to any specific meltwater discharge event. To avoid long-term drifts in mean ocean salinity caused by the long freshwater forcing, we impose a constant global mean salinity target (following the VFLUX method of Dentith et al. [2019]) commensurate with the starting condition for each "snapshot" ( Table 1). The salinity target conserves water in relation to terrestrial ice volume (applied as relative to the present day) and thus, in the context of these simulations, counteracts global freshening by removing the excess water as a very small proportion of freshwater from every ocean grid cell at every ocean model timestep (1 hr). This approach is in keeping with the snapshot/equilibrium experimental design, while still allowing the ocean to "feel" the surface forcing. See Section S3 in Supporting Information S1 for details.
Most simulations were run for 10,000 years; long enough to characterize their climate behaviors. However, like CTRL, two simulations (21.5k and 20.7k), were terminated after 4,000 years. At this point, little was changing in those simulations, and since no further time series were required for the analysis, we opted to conserve the computing resource.
Some simulations experienced numerical instability after a few thousand years observed through a stream function runaway off the coast of the Philippines. This quirk was resolved by smoothing the bathymetry in the region of the instability and restarting the run a few years before the instability arose. More information, including the detail of the smoothing algorithm and its very minor impact on the climate response, is given in Section S4 of Supporting Information S1.

Characterizing Oscillations and Defining Warm/Cold-Mode Composites
All experiments apart from CTRL show some kind of alternation between weaker and stronger AMOC phases or "modes." Changes in the AMOC are correlated to increases and decreases in NGRIP temperatures (North Greenland Ice Core Project, 42.32°W, 75.01°N) and so we will also refer to these phases as "cold" and "warm" modes, respectively. We do not use the terms stadial and interstadial to describe the cold and warm states because of the complicated connotations associated with these terms, but they may be thought of in such a light.
In order to characterize oscillations in the simulations, we applied a filtering algorithm and Fourier analysis to derive the spectrum of the temperature time series at the location of NGRIP, see Section S5 in Supporting Information S1 for details. If there is a peak in the spectrum, then oscillations can be defined and described, and we apply a Butterworth low-pass filter to screen-out the frequencies lower than the millennial-scale variability of interest.
We also found it useful, in our analysis, to examine characteristics common to all warm and cold modes in the suite of simulations. Thus, to build a composite of the two modes from the time series of results, we defined quantitative boundaries bespoke to each simulation (it proved ineffective to adopt a consistent definition for all simulations because of their differences). Points below the weak limits (in AMOC strength/NGRIP temperature) were added to the composite cold mode and points above the strong limits were added to the warm modes. This approach is described in Section S6 of Supporting Information S1, where we demonstrate that the choice of how  to define the composite modes does not significantly impact the results, and that to manually set up the weak and strong limits was an easy and robust method to build the composite states.

A New Weak, Shallow AMOC LGM Simulation
The CTRL run replicates and continues the HadCM3-GLAC-1D LGM simulation presented in Kageyama et al. (2021). The global mean surface temperature is 6.6°C colder than Pre-Industrial (PI). Compared to other PMIP4 simulations, this simulation is in the coolest range, almost 2°C below the average mean temperature, and is colder than any PMIP3 simulations analyzed by Kageyama et al. (2021), yet close to the current estimate from global temperature reconstructions (∼6.1 ± 0.4°C cooler than PI in Tierney et al. [2020], ∼7.0 ± 1.0°C cooler than PI in Osman et al. [2021]). The global mean ocean surface temperature cools by 3.4°C, which is significantly cooler than the ∼1.7 ± 0.1°C cooler than PI inferred by Paul et al. (2021), but again is a good match to the reconstruction by Tierney et al. (2020) (∼3.1 ± 0.3°C cooler than PI). Contrary to most PMIP4 models, HadCM3 produces an AMOC that is both shallower and weaker in a full-glacial background compared to its preindustrial state, although AMOC is still vigorous with a maximum strength of 20 Sv at 30° N. We also note that the simulation gets slightly cooler when using GLAC-1D compared to ICE-6G_C , despite similar values for the maximum overturning circulation.
A summary of the CTRL equilibrium North Atlantic climate is shown by Figure 2. The thermohaline circulation is fueled by intense convection in the northeast Atlantic, with deep water formation sites located primarily south of Iceland and west of the British Isles ( Figure 2b). This creates a corridor in the eastern part of the Atlantic where warm waters can transit to high latitudes, while the western part of the ocean gets covered by winter sea ice and observes a much cooler climate (Figure 2c). The winter sea ice layer extends toward the East of the North Atlantic basin, creating a strong North-South gradient in atmospheric temperature in this region (Figure 2a). This is slightly more extensive than what is predicted in Tierney et al. (2020) but matches the general layout, including a sea-ice free Iceland Basin. The Arctic Ocean is covered in sea ice all year long. Dense waters sink in the Labrador Sea when the sea ice is less extensive in this region in late autumn. For information, Figure 2 was reproduced for the Southern Ocean in Figure S3 of Supporting Information S1.

Climate Response to the Meltwater Perturbations
We observe significant differences in the climate response of the six meltwater experiments (Figure 3 The meltwater simulations can be assigned to three different regimes, according to the AMOC index ( Figure 3b). In the first regime, simulations 21.5k and 21k returned to the reference state after a short cooling event during the initial years of the forcing, when the ocean adjusts to the introduction of the weak meltwater fluxes. This AMOC  decline lasted for approximately 500 years, with the index weakening by as much as 5 Sv for 21k and 1.5 Sv for 21.5k approximately halfway through this initial cooling. While 21.5k hardly recorded a change of polar temperatures (Figure 3c), the drop in 21k AMOC strength drives up to 5°C cooling over Greenland. Neither simulation shows a strong response over Antarctica (Figure 3d). Accordingly, 21.5k and 21k will be referred to as warm simulations (note that with an LGM baseline climate, this term is relative).
In a second regime of its own, the AMOC in 18.2k almost entirely collapsed as soon as meltwater was discharged, causing Northern Hemisphere cooling and a significant shift southward of the ITCZ (Figure 3e). Short recovery episodes occur at irregular intervals through the 10,000 years of the simulation, but they cannot be sustained for longer than a few hundred years. This simulation will be labeled a cold simulation.
Finally, through Fourier analysis (Figure 4c; Section 2.3) we identify three oscillating simulations (20.7k, 19.4k, and 17.8k). They switch to a cold state at the onset of the run and continue in a quasi-oscillating regime between cold (near collapsed AMOC, similar to 18.2k) and warm modes (equivalent or stronger AMOC with respect to CTRL). The three oscillating simulations behave in a similar way and will be described in more detail in the following sections.
The oscillating and warm simulations show a clear distinction between an initial transient state and a steady state. The transient states corresponds to the adjustment of the model to the sudden introduction of meltwater. In the warm simulations, the transient state is seen in the initial drop, lasting ∼500 years. In the oscillating simulations, the transient state is manifested by an irregular first oscillation, lasting ∼1500 years. We cannot determine such distinction in the cold simulation.
During the cold modes ( Figure 5, top), we observe strong consistency between the oscillating simulations and the cold simulation, to the point where it becomes almost impossible to distinguish between them. These experiments show a maximum in atmospheric and oceanic cooling around 60°N (Figures 5a and 5b), where sea ice cover is now present both in winter and summer after the deep water formation sites vanished (Figures 5c-5e). A second peak of winter sea ice is noticeable around 40°N, but is not as clear in summer sea ice nor mixed layer depth. This second peak corresponds to the closing of the warm water corridor off the western coast of Europe and the spread of winter sea ice in this region. It is around these latitudes that we observe a maximal reduction of the thermohaline circulation by up to 14 Sv (Figure 5f). Because of the loss of convection at high latitudes, the upper cell of the AMOC largely dwindles north and south of 20°N. During the transient state, 21k follows a similar pattern. The soft decline of the AMOC is consistent with a shift southward of the convection sites, most likely resulting from a slightly cooler North Atlantic climate (Figures 5c-5f). We observe a clear reduction of sea surface temperature by as much as 2°C and of surface air temperature by up to 7.5°C at high latitudes (Figures 5a,  5b, and 5f). It is remarkable that 21k's winter sea ice expansion and mixed layer depth shallowing are comparable to the oscillating and cold simulations, demonstrating an increased seasonality compared to the reference CTRL state (Figures 5c and 5d). From this analysis, 21.5k is omitted (see Section S6 in Supporting Information S1).
We do not observe such consistency between simulations during their warm modes ( Figure 5, bottom). They all show a significant recovery from the cold modes, but it is impossible to underline a single common behavior. For example, despite a couple of periodic recovery phases of the AMOC in 18.2k (around 3,500, 7,000, and 9,000 years into the run; Figure 3b), the warm modes remain in a relatively cold-climate state (Figures 5h and 5i). At 60°N, sea surface temperatures are down 2°C and surface air temperatures drop by 6.5°C compared to CTRL. Shallower mixed layer depths around the same latitude indicate that Iceland/Irminger Basin and Labrador Sea convection sites are still greatly disrupted (Figure 5j). This keeps the sea ice edge far south in both summer and winter (Figures 5k and 5l). Amongst the oscillating simulations, the AMOC index is stronger than in CTRL (Figure 5m), increasing by as much as 2 Sv in the subpolar region. However, for the most part, these simulations also show disparities in key climate descriptors, at least in terms of amplitude of the anomalies. Surface temperatures and sea ice extent in 18.2k and 17.8k bear the closest resemblance across this subset of simulations, whereas  Figure S4 of Supporting Information S1. temperatures and sea ice extent in 20.7k and 19.4k indicate a slightly warmer climate (Figures 5h-5l) despite there being stronger ocean convection in 17.8k (Figure 5j). The oscillating simulations all show an increase in winter sea ice around 40°N compared to CTRL, corresponding to the narrowing of the warm water corridor along the coast of western Europe (Figure 5k). In summary, none of the simulations exactly returned to the CTRL reference climate during their respective warm modes, indicating significant regional legacy induced by the imposed meltwater patterns. We did not include the warm modes of the warm simulations as they would be indistinguishable from CTRL with this method.

Influence of the Meltwater Discharge
Abrupt climate changes are triggered by constant meltwater discharge in our simulations. Yet, we observe a strong nonlinearity between the climate response and the location and the influx of freshwater (Figure 3). For instance, despite a similar total influx of about 0.1 Sv, 19.4k is tipped into an oscillating regime while 18.2k ends up in a cold regime. On the other hand, 20.7k and 19.4k display similar oscillating dynamics even though the total Figure 5. Composite cold and warm modes' mean zonal anomalies between the meltwater simulations and the reference state in the Atlantic (70°W and 10°E). For cold modes (top), panels show the zonally averaged (a) surface air temperature, (b) sea surface temperature, (c) mixed layer depth, (d) winter sea ice concentration, (e) summer sea ice, and (f) maximum overturning circulation flow over the water column. For warm modes (bottom), panels show the zonally averaged (h) surface air temperature, (i) sea surface temperature, (j) mixed layer depth, (k) winter sea ice concentration, (l) summer sea ice, and (m) maximum overturning circulation flow over the water column. For orientation, an AMOC time series highlighting the periods contributing to the composite cold and warm modes is shown by panels (g and n), respectively. The warm and cold modes of 21.5k and the warm mode of 21k were excluded. flux is around 20% weaker in 20.7k than 19.4k. All of this hints at the importance of differences in the meltwater discharge pattern. While the spatial distribution of freshwater forcing is comparable between 21k and 20.7k, only 20.7k manages to generate oscillations. It demonstrates that not only the spatial distribution of the freshwater flux is an important control on the oceanic response, but that there is also a sweet spot in the perturbation conditions where the forcing needs to be strong enough to trigger a switch to stadial states, but not so strong that it (semi-)permanently suppresses a recovery as in 18.2k.
A threshold is reached in the cold climate modes, when the climate cools so far that it becomes insensitive to further meltwater discharge ( Figure 5). Independent of the forcing, all simulations produce a similar spatial pattern in temperature, sea ice, and deep water formation layout during the weak phases, but only when the surface atmosphere cools by as much as 15°C does the climate cross a tipping point where the cold phases can be sustained for a few hundred years. This corresponds to the vanishing of all oceanic convection north of 40°N (Figure 5c), resulting in an almost collapsed AMOC (up to 12 Sv weaker) in the North Atlantic (Figure 5f). When all the deep water formation sites have vanished, the response becomes decoupled from the forcing and only smaller regional effects can be induced. This phenomenon resonates with the conclusion of Smith and Gregory (2009). The effect of regional patterns of discharge on the climate response can be tracked by creating a composite of warm modes mixed layer depth and sea surface salinity ( Figure S5 in Supporting Information S1). It shows clearly distinct results between the simulations with prevalent North American meltwater inputs (20.7k and 19.4k) and the simulations dominated by Arctic discharge (18.2k and 17.8k). The warm simulations (21.5k and 21k) show almost no signal.
The effect of Arctic/GIN Seas discharge is decisive for triggering the shift from warm to cold AMOC states and leads to the strongest modification of the warm modes ( Figure 5). This could be the result of two factors. First, meltwater released from the Northern Eurasian ice-sheet leads to large salinity modifications by up to 8 PSU in the Arctic (Figures S5k and S5l in Supporting Information S1), which can massively disrupt the meridional salinity gradient. Second, meltwater is released relatively close to the main deep water formation sites in the GIN seas and in the Irminger/Iceland Basins, which it can reach by following the Greenland currents without having been significantly mixed with saltier warmer and saltier waters (Born & Levermann, 2010). Meltwater entering the GIN Seas should play a similar role, but the relatively low discharge in this area in our simulations makes it hard to conclude with certainty.
On the other hand, meltwater released off the northeast coast of North America has a weaker impact. Simulation 19.4k has greater discharge in the region and similar fluxes to the Arctic and GIN Seas compared to 20.7k. However, this increase in North Atlantic meltwater does not drive any significant ocean or climate response. From Figures S5i and S5j in Supporting Information S1, we infer that North American meltwater discharge is quickly dispersed. It has a direct effect on closing the Labrador Sea deep water formation site, but because this site is easily shut down in our simulations this does not impact significantly the AMOC. Because the water has to transit all around the subpolar gyre to target the more crucial sites in the eastern North Atlantic, the water will be mixed with tropical waters, weakening the forcing. The increased sensitivity to Arctic discharge compared to North American discharge in our simulations ties in with the conclusions of Roche et al. (2010) and Condron and Winsor (2012).
The existence of a sweet spot in the rate and location of meltwater discharge to the ocean for triggering climate transitions also depends on the background climate. Understanding the conditions leading to the creation of such a window of opportunity (Barker & Knorr, 2021) where abrupt climate changes can arise has been widely discussed (e.g., Brown & Galbraith, 2016;Klockmann et al., 2018;Zhang et al., 2017). Among the parameters likely to influence it, the choice of ice sheet reconstruction seems key, even more so in our simulations, as we rely on it both for the background climate state and for the meltwater forcing scenarios. In previous HadCM3-family deglaciation studies (e.g., Ivanovic et al., 2018;Gregoire et al., 2012), different ice sheets reconstructions did not yield as significant climate transitions, in spite of having a comparable baseline climate state (Kageyama et al., 2021) and magnitude of forcing. Other variables are also considered. Klockmann et al. (2018) and Zhang et al. (2017) only observed abrupt climate transitions under certain CO 2 concentrations. Model parameterization is also crucial, as for instance Peltier and Vettoretti (2014) managed to obtain oscillations by changing the vertical mixing coefficient. Precisely what about this particular reference state provides such a compliant framework for simulating AMOC oscillations may form the basis of future study.

Bimodal Warm States Linked to Reorganization of Deep Water Formation and Subpolar Gyre Layout
Intriguingly, the oscillating simulations frequently undergo a double warm peak during times of strong AMOC (Figure 4b). This is particularly clear in the late cycles of 17.8k (i.e., after 3,000 years), where the first peak in each warm phase is slightly cooler than the second. It indicates two distinct climate states separated by a cooling and then warming transition that is lower in amplitude than the cold-AMOC to warm-AMOC mode shift and which do not resemble the classical two-stage recovery hypothesis described by Renold et al. (2010) and Cheng et al. (2011). Changes in sites of deep water formation and the SubPolar Gyre (SPG) have been at the center of recent studies of abrupt climate and ocean circulation changes (Klockmann et al., 2020;Li & Born, 2019) and could shed some light on this warm-mode transition stage. Hence, we next examine the relationship between the AMOC index, the longitude of the Center of Mass of the Mixed Layer Depth (MLD COM) and different physical indicators. The center of mass does not indicate the actual location of the deep-water formation sites but should be seen as a tool to visualize the changes in the convection layout.
There is a linear relationship of roughly 1°C Sv −1 between the AMOC index and the temperature at NGRIP (Figure 6a), demonstrating once again that AMOC index and NGRIP temperatures are interchangeable when it comes to identifying climate modes in our simulations. The intensity of the convection sites follow a similar trend (Figure 6b), with roughly a 10 m increase of the maximum mixed layer depth in the North Atlantic for each degree of warming over Greenland. Both these relationships remain consistent irrespective of the specific warm modes that the oscillating simulations are in. However, for equivalent maximum AMOC strengths, the location of the deep water formation sites and the geometry of the subpolar gyre are clearly distinct between the two different warm modes. This is manifested in the ESPG index and latitude of the SPG COM, which split into two branches for high AMOC indices (Figures 6c and 6d). To further describe this phenomenon, we examine the different convection layouts in Figure 7. Two warm modes can be distinguished, each showing different deep water formation activity in the GIN and the Irminger seas. We use the peaks in their respective time series (Figures 7b and 7c) to build a composite mean of the MLD over the three modes (Figures 7d and 7f). The analysis was performed on 17.8k as it shows the clearest distinction between the two warm modes.
We follow the dynamics of the warm mode shifts by tracking the arrows in Figure 6, omitting in this analysis the initial spin-up period corresponding to values falling outside of the cycles (e.g., the first 300 years in 20.7k). Starting from a cold mode, the recovery of the AMOC is first fueled by convection in the GIN Seas, with very little signal in the Irminger basin, which is covered in sea ice in winter (Figure 7f). This results in a shift eastward and northward of the center of mass of the convection sites, as indicated by a deeper mixed layer (Figure 6h). Although situated in an area of intense sea ice formation, the Labrador Sea deep water formation site is not always reactivated during this mode. As a result, the subpolar gyre, which is initially weak and contracted during the cold phase, gets slightly stronger and extends eastward (Figures 6g, 6i, and 6j), entering the warm-meridional mode (so labeled for the disposition of the deep water formation sites during this phase). After sustaining a warm-meridional state for a few hundred years, the deep water formation layout is disrupted again to return to a state that resembles CTRL (comparing Figures 2-7d), with convection occurring primarily in the Iceland/Irminger basin and the resumption of Labrador Sea deep water formation (Figure 7e). The convection is distributed over a larger region and consequently the MLD index slightly dwindles. Conversely, the subpolar gyre intensifies, and moves southward and westward (Figures 6i and 6j). We call this state the warm-zonal phase.
Interestingly, we also observe short episodes of AMOC overshoot occurring immediately after the transition from cold to warm-meridional modes. The overshoots only exist at the onset of warm-meridional modes and are associated with stronger convection in the GIN Seas (Figure 7b).
Overall, these two warm states of strong AMOC are different to the two modes described by Cheng et al. (2011), where there is a transfer of deep water formation across the Atlantic from the Labrador Sea to the GIN Seas, the latter associated with an AMOC overshot. In all three of oscillating simulations, a relatively strong convection is maintained throughout both warm phases of the cycle, with permanent activity in the Iceland Basin and a transfer of deep water formation sites from the GIN seas to the Irminger Sea.

A Good Example of Dansgaard-Oeschger Events?
At first glance, the oscillating simulations resemble some recorded D-O events. From a purely descriptive point of view, presented in Table 2, we observe a periodicity (defined as the inverse of the dominant frequency in Figure 4) of between 1,540 and 1,930 years (the 18.2k simulation has a periodicity of 1,290 years, but strictly we  Figure S6 in Supporting Information S1 for zones definition), (d) the latitude of the center of mass of the mixed layer depth (MLD COM latitude); and between the longitude of the center of mass of the mixed layer depth (MLD COM longitude), (e) the temperature at NGRIP, (f) the MLD index, (g) the ESPG index, (h) MLD COM latitude, (i) the latitude of the center of mass of the barotropic stream function in the subpolar region (SPG COM latitude), and (j) the longitude of the center of mass of the barotropic stream function in the subpolar region (SPG COM longitude). The calculation of the center of mass is detailed in Section S8 of Supporting Information S1. Color shading indicates the location of the main areas of activity of the center of mass of the mixed layer depth during the different modes, as plotted on panel (k). Arrows indicate the direction of flow (through time) over the phase space; from green (cold) to blue (warm-meridional) to red (warm-zonal) modes. All time series were filtered following the algorithm presented in Section S5 of Supporting Information S1 and a 30-year running mean line is shown here, for readability. The black stars indicate the mean annual values calculated from the last 100 years of CTRL.
do not define this as an oscillating simulation). In terms of the duration of D-O cycles, this simulated periodicity is close to the range approximated from paleo records (about 1,500 years) during times of regular occurrence (Lohmann & Ditlevsen, 2019;Thomas et al., 2009). As an example, we compared our simulated cycles with DO 9-11, which occurred between 44 and 40 ka BP (Figure 8) and look the most like our oscillating simulations. In Table 2, we also compiled the main characteristics of oscillating simulations from three other models (Klockmann et al., 2020;Kuniyoshi et al., 2022;Peltier & Vettoretti, 2014). The information were extracted from the text and/or estimated from the plots when it was possible. Note that biases in the comparison may arise from the different resolution of the models (e.g., both CESM1 and MPI-ESM have higher spatial atmospheric resolu tion over Greenland compared to HadCM3) as well as from the definition of NGRIP temperatures (e.g., whether a smaller of larger area is used) and the filtering methods. The duration of the oscillations is also very similar to the simulated events of Klockmann et al. (2020) and Kuniyoshi et al. (2022), but two to three times longer than Peltier and Vettoretti (2014). As already discussed, 18.2k does not qualify as an oscillating simulation due to its cold characteristics (Section 4). Yet, in Figure 4, we observe two smaller peaks identified by the frequency analysis algorithm, one at ∼1,300 years and one at ∼3,500 years, which also correspond to typical D-O values (e.g., Kindler et al., 2014).
The temperature ranges between warm and cold climate/AMOC modes are similar, with changes of about 10°C recorded (Huber et al., 2006;Kindler et al., 2014) and simulated (this study) at NGRIP. From the other model studies cited in this study, only Peltier and Vettoretti (2014) obtained a similar amplitude of Greenland temperature change. Both Klockmann et al. (2020) and Kuniyoshi et al. (2022) observed smaller transitions similar more to the 6°C amplitude of 18.2k oscillations, but still within the lower range of reconstructed D-O events. The sea ice extent extends down to the Iberian margin in winter, which is what is consistent with what is observed during some Heinrich stadials (Voelker & de Abreu, 2011).
Notwithstanding some similarities with recorded D-O events, because all simulations were realized in a maximum glacial background (specifically, the LGM), the analogy between these simulations and D-O events is not straight forward. The comparison between our simulation and other occurrences of D-O events than DO 9-11 is often inconclusive. This is all the more true as we can also identify significant discrepancies between the model results and observations. For example, the simulated cycles do not match the typical shapes of D-O events (Lohmann & Ditlevsen, 2019); an abrupt warming followed by a slow cooling over a few hundred to a  Figure S6 in Supporting Information S1). (c) Mixed Layer Depth in the Irminger Sea (see Figure S6 in Supporting Information S1). The meridional, zonal, and cold modes were identified from the time series (pink, blue, and green vertical bands, respectively, in panels [a-c] few thousand years within the warm phase followed by a sharper cooling to finish off the cycle. Our simulations displayed a gradual transition between the cold modes to the warm modes, and the strong-AMOC phases were maintained rather steadily for ∼500 years before undergoing a slow cooling into the weak-AMOC phase of the Note. Periodicity was identified using the analysis presented in Figure 4c. Amplitude was defined as the difference between the warmest and coldest point of the filtered NGRIP temperature time series Figure 4b. Warming and cooling rates were defined as the maximum and minimum of the derivative of the filtered NGRIP temperatures time series, respectively. Descriptors of the simulations performed by Peltier and Vettoretti (2014) were estimated from Figure 1 by Vettoretti and Peltier (2018). Descriptors of the simulations performed by Klockmann et al. (2020) were taken from their main text. Descriptors of the simulations performed by Kuniyoshi et al. (2022) were taken from experiment P270 and estimated from Figure 3 by Kuniyoshi et al. (2022).  (Figure 4) and Compared to Simulations From Peltier and Vettoretti (2014), Klockmann et al. (2020), and Kuniyoshi et al. (2022)  cycle ( Table 2). The warming and cooling rates are less sharp than paleo records suggest, with a typical rate of 3°C temperature change per 100 years, about 10 times slower than indicated in Lohmann and Ditlevsen (2019). It has to be noted that despite a potential damping of the warming and cooling rates due to the filtering, the duration of filtered signals matched the unfiltered ones well, indicating that our conclusion on the shape of the simulated cycles are not an artifact of the analytical method. Kuniyoshi et al. (2022) similarly produced an oscillating ocean/climate with warming rates not exceeding 6°C 100 years −1 , in contrast to the 10 times faster warming rates simulated by Peltier and Vettoretti (2014) and Klockmann et al. (2020) (their cooling rates overlap with ours).
We observe an asynchronicity between the North Atlantic and Antarctica signal. This asynchronicity recalls the bipolar see-saw phenomenon described in Stocker (1998) and Stocker and Johnsen (2003). The Antarctica temperatures seem to lead the Greenland temperatures by 200-500 years in both the warming and cooling phases in Figures 8a and 8b, although it is difficult to identify a clear pattern. This does not fit the classical estimation (Antarctica lags the Greenland signal by ∼150 years in Svensson et al. [2020]) and challenges once again the resemblance of our oscillations to DO events. At last, it is possible that the periodicity observed is an artifact of the model when forced with constant sustained freshwater. By using a dynamical freshwater release pattern, we would be able to conclude on the comparability of these oscillations to actual paleo records.
Presently, we cannot categorically conclude whether or not our oscillations relate to D-O events. Differences in the shape of real and simulated D-O cycles may be explained by the quasi-idealized nature of our experiment design, specifically, the glacial maximum and fixed nature of our climate model boundary conditions/forcings (set to 21 ka, PMIP4 LGM protocol with GLAC-1D ice sheet plus the respective meltwater scenarios from the early deglaciation in GLAC-1D). The framework for our simulations was designed to simplify the identification of different behaviors in response to early deglacial meltwater forcing, mainly inspired by the hypotheses formed in conclusion to earlier work by Ivanovic et al. (2018). It provides a solid and systematic foundation for further work to study the mechanisms at play in the climate simulations presented here. However, this framework may interfere with or block some of the complex dynamics of D-O events and damp the abrupt climate changes. We indeed observe the sharpest cooling events at the onset of each simulation, indicating that the initial reorganization of deep water formation sites in response to a change in freshwater forcing may lead to the strongest and fastest climate disruption. Also, our setup does not include feedbacks between ice sheet melt and temperature changes, which could influence the periodicity of amplitude of changes Ivanovic et al., 2017). We therefore speculate that implementing a transient meltwater pattern consistent with what is known about past ice sheets during times of abrupt climate change could be a good way to better account for the dynamical interactions and abrupt reorganizations of the earth system, but ultimately, transient coupled climate-ice sheet simulations are needed to fully unlock the challenge of understanding D-O cycles and abrupt deglacial climate change. Finally, the discrepancies between observed D-O cycles and our simulations may also be related to weaknesses in the climate model itself. This version of HadCM3 seems to be unable to capture the fast physics component that has been observed in Vettoretti and Peltier (2018), and this may be related to the representation of ocean vertical diffusion (Peltier & Vettoretti, 2014).
Furthermore, and specifically to address D-O cycles, we also need to be able to reproduce the phenomenon of oscillating weak-strong AMOC modes outside of a glacial maximum background. Marine Isotope Stage 3 (MIS3, 29-57 ka BP, Lisiecki & Raymo, 2005), as depicted in Figures 8c and 8d for comparison to our results, has been often considered an appropriate candidate for such studies. Stadial conditions then were warmer than at the LGM, and it contains the most regular occurrences of D-O events recorded. Currently, one major challenge for setting off such a suite of simulations is that our oscillations are triggered by meltwater discharge, with a strong dependency on nuanced differences in the rate and location of freshwater inputs to the ocean. Thus, robust and detailed constraints on ice sheet extent are necessary to design an appropriate model experiment.
Apart from Heinrich stadials, MIS3 was not a time when ice sheet melting is thought to have been particularly strong, and there is likely to have been a lower meltwater flux than in our oscillating simulations (Batchelor et al., 2019;Hughes et al., 2016). However, DO 15-17 are believed to have been associated with changes of ice sheet extent (Lambeck, 2004), and because their shape remarkably resembles the observed oscillatory cycles of our simulations (Barbante et al., 2006;Erhardt et al., 2019;Rasmussen et al., 2016), they could be good candidates for being triggered by relatively low-levels of Northern Hemisphere ice sheet melt; it is possible that a weaker glacial climate state had a more sensitive ocean to smaller meltwater fluxes than our LGM-based simulations. However, the further back in time we go, the less constrained ice sheet extent (and geometry more broadly) and our model boundary conditions become, which again poses a practical limitation for how well such a climate model experiment could be designed to explore the detail of real past events. We think, therefore, that it is essential to obtain a deeper understanding of the physical mechanisms at stake during the climate oscillation before investigating further the relationship between our simulations and observed DO cycles.

Conclusion
Using snapshots of the meltwater discharge derived from the GLAC-1D ice sheet history of the early last deglaciation, we produced a set of oscillating simulations in the HadCM3 climate model under glacial maximum (PMIP4 LGM) conditions. Switching regularly between cold and warm modes, this behavior can only be attained in a narrow range of circumstances: if the freshwater forcing is too strong and/or applied in a particularly sensitive part of the ocean, the climate cannot fully recover to a warm mode and stays cold for the majority of the simulation. On the other hand, if the forcing is too weak, the transition to a cold mode is incomplete and the system rebounds to stay permanently in a warm state.
Understanding what can be defined as a "weak" or "strong" freshwater flux is not straightforward, as the response to the forcing is nonlinear. All simulations are in the range of plausible magnitudes for meltwater discharge over the last deglaciation-although the scenarios are made quasi-idealized by sustaining a constant meltwater discharge over thousands of years rather than following the transient history-but differ significantly in terms of total amplitude, and the geographical distribution of the freshwater fluxes. We observe that meltwater released in regions close to the main convection sites, namely the Nordic Seas and the Irminger/Iceland Basins, is the most effective at disrupting the AMOC. One possible inference from this finding is that in spite of being smaller in size compared to the North American ice sheet, Eurasian ice sheet demise punches above its weight in having the potential to trigger abrupt climate changes over the last deglaciation. Conversely, Laurentide ice sheet meltwater is required to be much more substantial to produce similar disruptions to Atlantic Ocean circulation and climate. The impact of the initial climate and ocean state on our results must also be considered. The cold interstadial state obtained with GLAC-1D ice sheets, with a relatively weak AMOC, convection concentrated around Iceland and extensive sea ice in the Western North Atlantic, creates a background that may favor the triggering of abrupt climate changes from Eurasian ice sheet melting.
When reaching a cold mode in our simulations, the ocean and climate response to meltwater forcing becomes decoupled from the direct influence of freshwater input at high latitudes. Because the AMOC is collapsed in these states, and the susceptible North Atlantic deep water formation sites vanished, meltwater perturbations will not propagate any longer, producing only smaller, regional perturbations. Cold modes correspond to an extreme cooling of the North Atlantic, and especially over former convection regions, with winter sea ice stretching down to Spain. We note a strong consistency in the climate change pattern of all simulations irrespective of whether they belong to the oscillating, cold, or to a lesser extent, warm regimes. This is not true for the warm modes, where the response is very dependent on the meltwater forcing scenario. We observe a range of different AMOC responses, from overshooting to damping, and the resulting climate is very different for each of the simulations. In particular, the simulations dominated by Arctic discharge never manage to completely return to the reference (CTRL) state. We also observe a two-phase, bi-modal warm condition in the oscillating simulations, with shifts between deep water-formation sites from a more meridional distribution (with primary convection situated in the Nordic seas) to a more zonal configuration comparable to CTRL. Resumption of Labrador Sea convection is inconsistent, sometimes being short lived and intense and sometimes being wholly absent, and in these glacial simulations, its connection to the upper cell of the overturning circulation is tenuous.
In summary, while we cannot conclusively determine that the oscillations we simulate are comparable to D-O events, there are resemblances that make such a comparison inviting. D-O events were rarely observed in full-glacial climates and are not believed to be directly forced by changes in ice sheet meltwater discharge. Still, the oscillating simulations presented here offer a valuable framework to analyze further the mechanisms behind the millennial-scale variability in glacial backgrounds. Intriguingly, they provide a set of insightful simulations that are so far unique in that relatively small influxes of freshwater trigger self-sustaining AMOC and Greenland temperature oscillations from glacial maximum conditions with a very large LGM North American ice sheet. Finally, the presented simulations pave the way for further study of ice sheet meltwater as a trigger, but not a direct driver, of abrupt climate changes within the last deglaciation.

Data Availability Statement
The new model data presented are available from the NERC Centre for Environmental Data Analysis repository (Rome et al., 2022). The presented plots realized were produced using the matplotlib python library and some of the colormaps were taken from Crameri et al. (2020). The code to reproduce the Figures from this article is accessible online (Olnavy, 2022c). It uses on the packages pylaeoclim_leeds v1.0 (Olnavy, 2022b) and mw_protocol v1.0 (Olnavy, 2022a), also created by Y.M. Romé.

Acknowledgments
YMR was supported by the Leeds-York-Hull Natural Environment Research Council (NERC) Doctoral Training Partnership (DTP) Panorama under grant NE/S007458/1. LJG was funded by MR/ S016961/1 and SST was funded by NE/ T007443/1. Simulations were run on the high performance facilities ARC3 and ARC4 at the University of Leeds, UK. All climate and ocean circulation proxy data and model simulations compared to are published elsewhere. Thanks are due to CEMAC for technical assistance on the numerical resource, with special thanks to Richard Rigby for assistance with the scripts for managing and processing the model output. Thanks to Julia Tindal and Steve Hunter for helpful technical discussions and investigations. Thanks to Kerim Nisancioglu and Chuncheng Guo for insightful discussion on the paper, and to Heather Andres whose presentation at vEGU21 inspired Figure 6.