Dynamics of Spontaneous (Multi) Centennial‐Scale Variations of the Atlantic Meridional Overturning Circulation Strength During the Last Interglacial

Recent reconstructions of bottom water δ13C during the last interglacial (LIG) suggest short‐lived variability in the Atlantic meridional overturning circulation (AMOC). Spontaneous (multi) centennial‐scale variability of the AMOC simulated in the Earth system model of intermediate complexity iLOVECLIM are investigated for that period. The model simulates abrupt AMOC transitions occurring at 300 years frequency and correspond to a switch of the AMOC vigor between a strong (∼17 Sv) and a weak (∼11 Sv) state. The onset of these abrupt transitions is associated with changes in orbital forcings resulting in the decline of summer insolation in the high latitudes of the North Atlantic and affecting the sea ice cover in two key deep convection regions: (1) the northern Nordic Seas (NNS) and (2) the northwest North Atlantic (NWNA). Northward inflow of Atlantic surface water increases the convection depth in (1) and strengthens the Greenland Iceland Norway (GIN) Seas overturning circulation. Subsequent ocean‐atmosphere interactions involving sea ice, ocean heat release, anomalies of evaporation‐precipitation, and wind stress over the Nordic Seas lead also to an increase in deep convection in (2), followed by increase in the AMOC strength.


Introduction
Large and abrupt changes in the Atlantic meridional overturning circulation (AMOC) associated with dramatic climate shifts may have repeatedly occurred during the past Glacial cycle (Rahmstorf, 2002) and have been associated with freshwater input from ice-rafted detritus into the North Atlantic (Bond & Lotti, 1995;Broecker, 1994;Heinrich, 1988) and disturbances in the North Atlantic Deep Water (NADW) formation (Henry et al., 2016). However, the specific role of icebergs and freshwater forcing in driving the deep circulation anomalies is still not resolved (Barker et al., 2015). In contrast, abrupt AMOC variability is less evident during warmer periods such as the late Pleistocene interglacial periods with both model and proxy studies pointing toward a generally vigorous NADW production on millennial and longer time scales (Adkins et al., 1997;Kuhlbrodt et al., 2007;Oppo et al., 1997). The Intergovernmental Panel on Climate Change (IPCC) Assessment Report 2013 and an increasing number of studies (Cheng et al., 2013;Rugenstein et al., 2013;Schmittner et al., 2005;Weaver et al., 2012) conclude that rapid changes in the AMOC are very unlikely to occur by the end of the 21st century, suggesting rather a slow decline from the current state (11% for RCP2.6 and 32% for RCP8.5, in the ensemble means). However, Castellana et al. (2019) evaluate the risk of sudden decline and collapse of the AMOC at 15% and recent studies have revealed bias in the current climate models concerning their AMOC stability. Models have been classified as having either monostable or multistable regimes depending on the sign of their overturning-related freshwater transport at 30-32°S in the Atlantic (Weaver et al., 2012). Models with different regimes would likely simulate different AMOC responses to a given climate forcing. It is found that a majority of climate models feature an AMOC that imports fresh water into the North Atlantic, that is, monostable (Drijfhout et al., 2011;Weaver et al., 2012), in contrast to observations-based analysis, which suggest that the current AMOC is in bistable regime (Bryden et al., 2011;Garzoli et al., 2013;McDonagh & King, 2005;Weijer et al., 1999). This means that the current thermohaline (THC) component of the AMOC, driven by the salinity and temperature fields, would have multiple modes as it may have had in the past (Rahmstorf, 2002).
In this context, an increasing number of model studies, from the pioneer box model of Stommel (1961) to more recent and complex climate model (Liu et al., 2017;Yin & Stouffer, 2007), have been performed to analyze the capability of climate models to switch from monostable to multistable regimes; highlighting the salt advection feedback as the responsible process for multiple equilibria (Weijer et al., 2019). Other model studies showed that internal processes in the ocean-atmosphere system could generate spontaneous abrupt variations of the AMOC through their influence on the vertical density gradient in the North Atlantic (Arzel et al., 2012;Brown & Galbraith, 2016;Drijfhout et al., 2013;Friedrich et al., 2010;Ganopolski & Rahmstorf, 2001;Peltier & Vettoretti, 2014;Sakai & Peltier, 1997;Schulz et al., 2002;Vettoretti & Peltier, 2015;Winton, 1993). Thus, AMOC oscillations associated with the weak and strong modes of the THC have been analyzed for various climate conditions. However, to the authors knowledge, there have been few efforts to explore the potential for spontaneous changes of the AMOC during warm interglacial boundary conditions. Yet reconstructions reveal that large short-lived changes in NADW may have been plausible during this period (Galaasen et al., 2014;Hodell et al., 2009) and potentially the AMOC  may have occurred during the last interglacial period.
Here we use a set of experiments under the last interglacial (LIG) forcing conditions to provide insights into the possible reasons for (multi) centennial oscillations of the AMOC using the Earth system model of intermediate complexity (EMIC) iLOVECLIM (Bouttes et al., 2015;Goosse et al., 2010). Specifically, we aim to delineate the thresholds and internal mechanisms within the model responsible for simulating spontaneous, large variations in AMOC strength under the LIG boundary conditions. While Kessler et al. (2020) discuss AMOC-related changes in ocean carbon cycling, biogeochemistry, and compare to proxy (benthic carbon isotope reconstructions, e.g., Galaasen et al., 2014;Hodell et al., 2009), the internal mechanisms behind these oscillations were not thoroughly evaluated due to the limitation of this previous study design. Here we provide additional experiments in order to delineate the importance of different forcing mechanisms and achieve a mechanistic understanding of the abrupt AMOC variations that could have occurred under the LIG climate conditions. Thus, we use the iLOVECLIM model to provide a dynamically consistent interpretation of the mechanisms underlying rapid carbon cycle and climate variability found in recent reconstructions.

Model Description
This study uses the iLOVECLIM model, a carbon isotope-enabled version of LOVECLIM composed with similar atmosphere, ocean, and vegetation components to that in LOVECLIM version 1.2 (Goosse et al., 2010;Roche et al., 2007). The atmospheric component ECBilt was developed at the Dutch Royal Meteorological Institute (Opsteegh et al., 1998) and has a resolution of 5.6°with three vertical layers at 800, 500, and 200 hPa. The atmosphere is represented by quasi-geostrophic approximation to which ageostrophic terms are added in order to better represent the circulation at low latitude, that is, the Hadley cell. The oceanic component CLIO (Coupled Large-scale Ice Ocean) adopts 20 levels on the vertical coordinate and has a horizontal resolution of ∼3°. The oceanic general circulation model (Goosse & Fichefet, 1999) is based on the Navier-Stokes equations written in a rotating frame of reference including a parameterization for downsloping currents (Campin & Goosse, 1999). This means that shelf water grid box with higher density than the neighboring will flow along the slope until it reaches a depth of equal density while a vertical and horizontal return flow from the deep ocean to the shelf compensate this volume transport. In addition, the usual approximations are employed such as the hydrostatic and Boussinesq approximations, while an updated version of the thermodynamical sea ice component of Fichefet & Maqueda (1997 is also included. The sea ice growth/decay rates are governed by bottom and surface sea ice prognostic energy budgets in leads, where sea ice thickness increases when the snow cover is heavy enough to push the entity under the water level and freeze with the infiltration of seawater. The parameterization of McPhee (1992) is applied for the heat flux from the ocean to the ice and mass conservation controls salt and freshwater surface exchanges. The terrestrial biosphere model (VECODE) has three submodels that exchange heat, stress, water, and carbon (Brovkin et al., 1997). Trees and grass are simulated in addition to deserts and subdivided into four components (leaves, wood, litter, and soil) that exchange carbon, while the fraction of plant functional type (PFT) are computed at equilibrium with the simulated climate. The photosynthesis depends on precipitation, temperature and atmospheric CO 2 . Finally, the ocean carbon cycle model applies a nutrient-phytoplankton-zooplankton-detritus (NPZD) model described in Six and Maier-Reimer (1996) and in addition includes organic and inorganic components as described in Bouttes et al. (2015).

Experimental Setup
A set of experiments was conducted to study the mechanisms that generate internal (multi) centennial-scale variability in the AMOC observed in the transient simulation of Kessler et al. (2020) (hereafter referred as CTRL). We have performed three additional simulations over the 125-115 ka period of the LIG (summarized in Table 1) using prescribed yearly interpolated values of greenhouse gases (CO 2 , CH 4 , and N 2 O) and orbital forcings from the third phase of the Paleoclimate Modelling Intercomparison Project (PMIP3; https://pmip3. lsce.ipsl.fr/) and compare them to CTRL. The CTRL experiment consists of a transient simulation under natural greenhouse gas and orbital forcings. In the three additional experiments, we repeated the CTRL simulation but with (1) fixed greenhouse gases concentration to the value at 125 ka (CO 2 = 276 ppmv, CH 4 = 640 ppb, N 2 0 = 263 ppb; this experiment is referred as GHG fixed ); (2) prescribed land vegetation to the value at 125 ka (hereafter referred as VEG fixed ); and (3) fixed orbital forcings to their values at 125 ka (hereafter referred as ORB fixed . While ORB and GHG experiments allow us to differentiate the influence of orbital and greenhouse gas forcings, respectively, the VEGF experiment allows us to determine the role of the terrestrial vegetation-induced climate feedback on the AMOC. The prescribed land vegetation value is taken from a 5,000 years equilibrium simulation under the 125 ka boundary conditions. As in the CTRL experiment, these three additional experiments are branched after two consecutive runs in order to set the 125 ka initial conditions; a preindustrial spin-up followed by an equilibrium experiment using 125 ka boundary conditions. Both experiments are simulated over 5,000 model years and are described in more detail in section 2.2 of Kessler et al. (2020).

AMOC Response to Varying Forcings
Deep water formation in our model occurs south and north of the Greenland-Scotland Ridge. The northern formation region is primarily located in the eastern Greenland Iceland Norway (GIN) Sea, while the southern formation region is near south Greenland and south of Iceland. We define the AMOC index as the maximum of the overturning stream function at 27°N. Furthermore, we define the northern branch of the AMOC as the GIN Sea overturning circulation (GSOC) index, which we characterize by the maximum of the overturning stream function north of 70°N. Figure 1 depicts the simulated variations in AMOC (pink) and GSOC (cyan) for the four experiments (CTRL, GHG fixed , VEG fixed , and ORB fixed ). The obliquity and June insolation at 60°N forcings are shown in panels (a) and (b) and are decreasing from their maximum values (23.8°and 528.6 W m −2 , respectively) toward 115 ka and the Glacial inception.
During the first 3,000 years the overturning indices are relatively stable in all experiments after which the model simulates an increase in variability and numerous occurrence of spontaneous (multi) centennial-scale transitions for the experiments CTRL, GHG fixed and VEG fixed (Figures 1c-1e). In these three experiments, the second half of the LIG is characterized by an AMOC that preferentially operates at a weak    Figure 1c). The dashed vertical lines denote averaging intervals for panels P1-P6 of this figure. Panels P1-P6 represent the SSS (shade), the convection depth (black thick lines), the sea ice limit (light blue lines), and the averaged values of the AMOC (pink) and GSOC (cyan) indices for the corresponding intervals.

10.1029/2020PA003913
Paleoceanography and Paleoclimatology and) preindustrial boundary conditions. They obtained a similar AMOC behavior in their model experiment when reducing the obliquity from 23.446°to 22.1°. Similar behavior has also been found by Jongma et al. (2007) when applying external freshwater forcings using the ECBilt-CLIO coupled model. In our simulation, the shift in AMOC variability occurs at a higher obliquity of 23°, which occurs at around 120 ka in CTRL ( Figure 1c) and GHG fixed (Figure 1d) experiments and corresponds to a solar insolation at 60°N similar to present day near 460 W m −2 ( Figure 1b). In VEG fixed the system is unstable around 120 ka and seems to switch toward a preferentially weak AMOC state later, at around 118.5 ka ( Figure 1e) corresponding to an obliquity of about ∼22.8°. In our simulations the initiation and ending of the AMOC/GSOC oscillations show some distinct differences. While the oscillation amplitudes seem to gradually increase in CTRL and VEG fixed (Figures 1c and 1e), they abruptly begin around 120 ka in GHG fixed (Figure 1d). At the end of the simulations the abrupt variations remain relatively strong in VEG fixed , moderate in GHG fixed while they stabilize at weak state in CTRL. No abrupt transitions are simulated for the run only forced by varying greenhouse gases (ORB fixed , Figure 1f), which also show a significant reduction in standard deviation (±1.1 Sv against ±2.5/2.6/2.7 Sv in CTRL/VEG fixed /GHG fixed , respectively). This suggests that the climate generated by the changes in orbital forcings is the main driver for (multi) centennial-scale variations to occur in our experiments and that these unforced transitions can occur spontaneously (in the model) under LIG climate conditions. Also similar to Friedrich et al. (2010), we note that the changes in GSOC indices are simulated several decades ahead of the AMOC indices (this can be seen in more detail in Figure 2 (top panel) later in the manuscript), suggesting that changes in the GIN Seas are antecedents of the changes south of Greenland (and the global ocean circulation), highlighting the importance of this region as a key area for altering AMOC in the model. Figure 3 shows the main mode for these abrupt transitions by depicting the power spectral density of the AMOC for each experiment. The LIG is divided into two window periods taking the 120 ka as a pivotal period, which represents a transition for an ocean circulation with a preference for a strong AMOC state to another with a preference for a weak AMOC state. Figure 3 shows that, for experiments CTRL,GHG fixed , and VEG fixed , the second part of LIG (120-115 ka, blue lines) is dominated by these abrupt transitions with a main oscillation frequency of around 300 years per cycle, depicting a power spectrum density up to 10 times stronger as compared to the early LIG window (125-120 ka, red lines). This centennial time scale is similar to that depicted by Friedrich et al. (2010). The experiment with fixed orbital forcing (ORB fixed ) shows similar variability in both time windows.

Regional Deep Convection
Here we only show the analysis over the CTRL experiment, but the same analysis has been performed for the other experiments and revealed similar results. We analyze four domains (NWNA, NENA, SNS, and NNS; see  contrast, the variations of deep convection in NENA and SNS (Figure 4b, red and light orange curves) are in antiphase with the overturning circulation indices; decreasing when the AMOC strengthens. In these regions, the convection remains active during the entire simulation and even increases by twofold in SNS toward 115 ka (Figure 4b, red curve). As a result, the convection occurring in NENA and SNS is mainly suggested to sustain the mean overturning circulation state, as they oscillate between a strong and a weak state (Figure 1c). On the other hand, the NWNA and NNS regions are suggested to mainly contribute to the amplitude of the overturning circulation through enhancement (reduction) of strength during strong (weak) state.

Physical Ocean Surface Changes
Figures 5 and 6 depict the changes in physical properties at the ocean surface such as the vertical gradient of density over the first 600 m (gradρ), the sea surface salinity (SSS), the percentage of sea ice cover, and the sea surface temperature (SST) in each region.
In the regions that show an evolution of convection depth in phase with the overturning indices (NWNA and NNS, Figure 4b, top panel), the early LIG is characterized by a relative low stratification (i.e., Figure 5a) while it is increased by twofold in the late LIG in NWNA and NNS with values around 3 and 2 g m −4 , respectively. This is consistent with the previously mentioned changes in convection depth occurring in these regions. The changes in SSS (Figure 5b) covary (antiphase) with the vertical stability (r = 0.98) suggesting that the surface salinity plays a key role in triggering the convection in these regions. In the early LIG relatively high SSS are simulated in NNS and NWNA with values around 34 and 35 psu, respectively. This period corresponds to low sea ice cover (Figure 5c) around 50% in NNS and 10% in NWNA, in the model. High sea ice cover is associated with fresher surface water; hence, when sea ice cover extends, SSS drops, and vice versa. For example, when NNS is fully covered by sea ice (Figure 5c) the SSS is reduced to 33 psu. Here the simulated SST (Figure 5d) is in phase with convection depth and in antiphase with vertical stratification. When the stratification is weak (strong), warmer (colder) SST are simulated. From 125 to 115 ka, the SST declines by about 4°C in both regions and is consistent with the simulated increase of sea ice cover (Figure 5c). However, this cooling does not seem to affect directly (i.e., increase) the convection in these regions, where the influence of surface salinity and fresh water appear dominant.
In regions where the evolution of the convection depth does not correspond to the changes in overturning indices (NENA and SNS, Figure 4b, bottom panel), the vertical stratification ( Figure 6a) remains weak during the entire simulation depicting values below 1 g m −4 . This is consistent with the persistent simulated convection in these regions. Similarly, the simulated SSS (Figure 6b) varies very little and remains above 35 psu, setting favorable conditions for convection to persist. The influence of the sea ice in these regions remains minor as no sea ice is simulated in NENA while it grows only up to 25% of the total area in SNS. As SSS varies little, the changes in convection depth are more affected by the changes in SST (Figure 6d). Consequently, the simulated long-term surface cooling toward 115 ka (∼2.5°C) in each of these two regions is in good agreement with the slight decreasing trend in vertical stratification (∼0.3 g m −4 in 10,000 years, Figure 6a). The increase of the stratification during the abrupt transitions depicted in NENA (Figure 6a)  (Figures 7c and 7g).
In our model, the simulated abrupt increases in SSS in the NNS region, are induced by northward intrusion of Atlantic water (more saline) into the GIN Seas also referred to as "Atlantification" of the polar regions, a phenomenon which has been observed in the past decades (Årthun et al., 2012) and is important for setting the boundaries of sea ice in the Barents Sea. As this region becomes fresher toward 115 ka, the intrusion of saline Atlantic surface water generates even larger modifications. This sudden increase in SSS generates deep convection in NNS and in the regions south of Greenland (NWNA) few decades later. The complete set of processes linking NNS to NWNA will be discussed in section 3.5. On the other hand, the NENA and SNS regions depict, in general, higher mean surface salinity and appear to be more resilient to the climate evolution over the simulation period, with a higher influence from changes in surface temperature than changes in surface salinity.

Role of Summer Air Temperature
The high-latitude cooling toward 115 ka appears to be a key factor for initiating the high-frequency AMOC transitions at the mid-LIG. This cooling is particularly important during the summer as it modulates the melting season of the sea ice and consequently the remaining total sea ice cover during the deep convection season (winter). Figure 8 shows the time series of summer air temperature over the latitude bands 69-80°N and depicts a general and persistent cooling in all of the experiments that exhibit abrupt AMOC variations (CTRL, GHG fixed , and VEG fixed ). The CTRL and GHG fixed experiments (Figure 8, black and light blue lines) simulate relatively similar summer temperature where GHG fixed is only 0.3°C warmer during weak AMOC states at the end of the LIG. This may be due to the small difference in atmospheric CO 2 concentration used in the radiative budget calculation. Indeed, at the end of the LIG the constant value of atmospheric CO 2 used in GHG fixed is slightly higher (up to 10 ppm) than the one used in CTRL (prescribed values from PMIP3). The VEG fixed experiment is, in general, warmer than CTRL and GHG fixed . Consequently, the oscillations are delayed until around 118.5 ka. The temperature in ORB fixed remains relatively stable and warm during the entire simulation period, allowing most of the sea ice produced in the winter to melt during the

10.1029/2020PA003913
Paleoceanography and Paleoclimatology summer. Notably, in each experiment with high AMOC variability (CTRL, GHG fixed , and VEG fixed ) the pronounced AMOC oscillatory state begins in each case around the time period when summer temperature decreases to roughly −4°C at 75°N. This occurs slightly before 120 ka in CTRL and GHG fixed and later at 118.5 ka in VEG fixed due to the higher temperatures in this simulation. This threshold marks the limit below which the summer sea ice does not melt in response to the atmospheric forcing and is therefore only affected by water mass changes. This confirms that the NNS, which is the most impacted by sea ice cover, is the key region for initiating the high-frequency AMOC variations in the model.
Our model simulations demonstrate the key role of air temperature in establishing the seasonal sea ice coverage and consequently the timing of AMOC phase transition. Figure 9 highlights for each experiment the mean seasonal air temperature at 75°N at both the beginning and end of the large AMOC oscillation phase. The CTRL experiment simulates the coldest seasonal cycle at both the onset of the oscillations and at the end of the LIG (Figure 9, black solid and dashed lines). The GHG fixed and VEG fixed experiments are on averaged warmer than CTRL by 1°C and 3.6°C, respectively. As previously mentioned, the warmest period (July) remains below −4°C. The gray shading (Figure 9) highlights the range of temperatures for which abrupt transitions are simulated in our experiments. Above this range, the atmospheric temperature is warm enough during the summer to considerably reduce the summer sea ice extent and therefore the averaged sea ice extent over the year. With a retreated sea ice, more Atlantic water is able to penetrate into the northern Nordic Seas and south of Greenland increasing the surface salinity ( Figure 5b) and setting the AMOC into its strong state. While the initiation of the AMOC oscillation seems to be linked to a cold enough climate for pushing the system into its weak state, the termination period of the large oscillations phases remains less clear. Do the oscillations stop as suggested by the CTRL run (Figure 1c) or do they just become less frequent as the climate continues to cool down? One hypothesis could be that the sea ice becomes too thick to be melted only by ocean heat flux. The sea ice is indeed approximately 10 to 15 cm thicker in CTRL than in GHG fixed . Additional experiment on sea ice and ocean heat flux sensitivity would be needed to give more insight into the long-term AMOC stability beyond our experimental period.

Centennial-Scale AMOC Variability and Mechanisms
In the previous subsections, we highlight that changes in GSOC are antecedent to changes (increase) in total AMOC which are associated to greater convection depth south of Greenland (also defined here as the NWNA region). In this subsection, we connect the processes occurring in the NNS region, that are responsible for the simulated increase in GSOC indices, to the changes occurring in NWNA and in AMOC indices. To do so we show the example of the last oscillation of the CTRL run (Figure 1c, gray shade), which we divide into several phases. Figures , 2,10, and 11 illustrate atmosphere-ocean

Paleoceanography and Paleoclimatology
interactions during each of these phases and Figure 11 shows the evolution of the temperature (T) and salinity (S) vertical profiles in each of the four regions presented in Figure 4a.
The phase "P0" (Figure 2, top panel) corresponds to a weak AMOC state characterized by an extended sea ice cover over NNS and NWNA comparable to what is simulated during P6 (Figure 2, P6, light blue line). During this phase, heat accumulates in NNS below 200 m depth as the sea ice covers up to 100% of the region reducing the heat flux toward the atmosphere. This process weakens the stratification of the water column, Figure 10. Anomalies evaporation-precipitation (E-P) in cm yr −1 (shade) and air temperature in°C (contour lines) for each of the phases described in Figure 2 (top panel). These anomalies take the interval P0 for reference. P6 does not depict any significant changes in air temperature as compared to P0.

Paleoceanography and Paleoclimatology
preconditioning the NNS region for deep convection. The first phase (P1) is characterized by the activation of deep convection in NNS (Figure 2, P1). It is initiated by the inflow of saline (>34.75 psu) and warm Atlantic surface waters, which abruptly increase the SSS in this region by 1 psu (Figure 7e) and the SST by 2°C (Figure 7a). The SSS increase overwhelms the SST increase, therefore the convection starts and GSOC index increases (7 Sv). Phase 2 (P2) is characterized by a reduced sea ice cover induced by both the advection of warm Atlantic surface waters and the mixing of the interior ocean heat. As a result, more heat escapes from the ocean to the atmosphere generating a positive anomaly in air temperature up to +10°C over the NNS (Figure 10, P2, green contours). Associated to the increase in surface air Figure 11. Anomalies wind stress in N m −2 (shade) and its direction (arrows) for each of the phases described in Figure 2 (top panel). These anomalies take the interval P0 for reference.

10.1029/2020PA003913
Paleoceanography and Paleoclimatology temperature, the model simulates an increase in evaporation-precipitation (E-P, Figure 10, shade), which gradually extends southward toward the NWNA region. This surface flux response acts as a positive feedback further increasing the surface salinity in both regions, although the largest positive SSS anomalies occur in regions where the fraction of high-salinity Atlantic water increases. The warm anomaly generates a modification in the wind stress ( Figure 11, P2-P4) increasing the easterlies flow south of Iceland and advecting more Atlantic waters toward the NWNA. Consequently, both the E-P and the Atlantic water influx increase the SSS in NWNA by up to 1.5 psu (Figure 7f) in P3 relative to P0 and the convection depth deepens. This gradually strengthens the total AMOC, which reaches its maximum in P3 (17 Sv, Figure 2, P3). Hence, P3 is characterized by both maximum AMOC and GSOC strengths, which lasts for about 50 years. The exact reason for this time scale remains difficult to determine, however, a persistent AMOC-induced northward freshwater transport at 33°(∼0.3 Sv) sets the model simulations in a monostable regime. This freshwater flow increases by about 0.05 Sv when the AMOC strengthens, slowly acting as a restoring force (e.g., Liu et al., 2013Liu et al., ,2015. As a result, GSOC strength begins to decline in P4 (Figure 2, P4) and the associated heat loss in the NNS reduces. We note that the temperature below 200 m depth in NNS decreases by about 2°C between the beginning of P1 to the end of P3 (Figure 7a). This cooling could also negatively affect the magnitude of the heat loss toward the atmosphere and impact the time scale of AMOC recovery. The air temperature anomaly diminishes over the NNS from 10 to 4°C (Figure 10, P4), and the sea ice formation gradually returns contributing to the surface freshening ( Figure 7e). The convection in the NWNA remains active for the proceeding 60 years as the anomaly of salinity remains high (Figure 7f). This is the result of both persistent strong anomalies of wind stress (Figure 11, P4) and E-P ( Figure 10, P4) near Greenland, responding to the remaining +6°C air temperature anomaly (Figure 10, P4). In phase 5 (P5), the salinity anomaly continues to reduce in the NNS area (/sim0.4 psu) with a complete return of the sea ice and GSOC declines to 3 Sv. At this point, the anomalously warm air temperature as well as the induced wind stress anomalies over the entire North Atlantic disappear. Finally, the sea ice starts to regrow in the NWNA region. The AMOC strength declines toward its weak state (13 Sv). In phase 6 (P6), the total AMOC returns to its weak state (11 Sv) and the GSOC is almost shut down (1 Sv) as the Atlantic surface waters retreated further south, leaving the NNS and NWNA regions covered by sea ice and fresher surface water.
Figure 7 allows us to further analyze the effect of the temperature and salinity regarding each region. The NNS and NWNA regions depict near-surface salinity changes of about 1.5 psu (Figures 7e and 7f) which overwhelm the associated temperature increase (Figures 7a and 7b) and control the activation of the deep convection. On the other hand, SNS and NENA depict stronger changes in temperature (Figures 7c and  7d, up to 3.8°C) relative to salinity (Figures 7g and 7h, only 0.5 psu maximum), which sets the temperature as main driver for deep convection in these regions. This is slightly less pronounced in NENA where the anomalies in both temperatures and salinity are smaller.

Discussion and Conclusion
We investigated the mechanisms responsible for plausible spontaneous (multi) centennial-scale AMOC shifts in the iLOVECLIM model under the LIG boundary conditions. Our findings reveal that abrupt transitions of the AMOC from strong to weak states may also have been plausible during the LIG under natural greenhouse gas and orbital forcings. These AMOC shifts are comparable in time scale to that suggested by proxy marine records from that period (Galaasen et al., 2014Hodell et al., 2009), which have been interpreted as modifications of deep water mass geometry and potentially AMOC strength .
The two simulated AMOC states are characterized by the activation and deactivation of deep convection areas in (i) the northern Nordic Seas (NNS) followed several decades later by changes in (ii) the Northwest North Atlantic (NWNA). These (de)activations modulate the AMOC strength and are shown to be associated with the alteration of the vertical stratification due to salinity/fresh water influence. In our model, the changes in water mass stability is initially driven by the growing presence of sea ice, which hinders the "Atlantification" process (advection of warm and saline surface water) into the NNS and NWNA regions. The major role of the salinity/fresh water depicted by our study in affecting the vertical stability has also been shown in previous studies applying ECBilt-CLIO coupled model (Jongma et al., 2007;Schulz et al., 2007) and LOVECLIM (Friedrich et al., 2010). Although the AMOC oscillations in Jongma et al. (2007) arise from external fresh water forcing into the Labrador Sea, the general processes are still similar to those found in our study with a change in surface salinity/fresh water into this region (here approximately corresponding to NWNA) strongly modifies the AMOC strength and can occur spontaneously in the model under LIG boundary conditions.
In this study, the decrease of high-latitude insolation, primarily due to the declining obliquity, is critical for initiating the (multi) centennial-scale oscillations of the AMOC strength. These variations occur effectively as jumps between weak (∼11 Sv) and strong (∼17 Sv) AMOC states. As insolation declines there is a gradually increasing preference toward the AMOC weak state with less frequent episodes of strong AMOC. We found a common threshold to all our experiments depicting spontaneous AMOC strength variations (CTRL,GHG fixed , and VEG fixed ): ∼ −4°C in summer's (JJA) air temperature over the latitude bands 69-80°N. We suggest that below that temperature, the sea ice only responds to changes in surface water temperature (i.e., advection of warm Atlantic surface water) as the atmosphere is cold enough to maintain the sea ice cover. Consequently, the annual sea ice extends abruptly in NNS and NWNA, lowering the SSS and reducing the deep convection in these regions, and thereby altering GSOC and AMOC. Yet it is not clear if with further cooling, the oscillations of the AMOC would stop (as suggested by CTRL experiment) or become more rare as it would demand more energy from the ocean to melt a thicker sea ice cover. We note that Bakker et al. (2013) also reported an abrupt decrease of the AMOC strength around 120 ka in their LIG experiment with the LOVECLIM model. However, they linked it to a sudden decrease of January's air temperature over the latitude bands 60-90°N. In our simulations, a similar cooling is simulated during the winter months (DJF, 7°C, not shown here) when the AMOC starts to oscillate, but it warms back up to its previous value until the next abrupt transition. On the contrary, the summer temperature (JJA) shows less variability which is the reason why we suggest it as a crucial threshold. We note that the July air temperature in Bakker et al. (2013) seems also to be slightly above −4°C when the AMOC abruptly weakens, which is consistent with our finding.
We show that internal processes in the atmosphere-ocean system, strongly supported by the presence of sea ice, trigger the simulated abrupt transitions. First, the sea ice presence in key convection regions initiates the AMOC swings by turning off the convection. When the NNS is fully covered by sea ice, heat accumulates below 200 m depth for up to 150 years ( Figure 12). This heat corresponds to an average temperature increase of about 2°C, which helps destabilizing the stratification and creating a precondition for deep convection in the NNS region. Friedrich et al. (2010) suggested the deep decoupling as the trigger for deep convection, while it is not the case in our study. In their study, this arose from a subsurface temperature increase of about 4°C which is stronger than what is simulated in our study. Our results support the idea that the main trigger affecting the convection regions remains the changes in salinity/fresh water occurring in the upper ocean layers. Here, we refer to this generically as the relative influence of Atlantic surface waters in NNS and NWNA, although E-P feedbacks enhance the initial changes. We note that this advection of saline Atlantic water seems to dominate over the potential contribution of Arctic sea ice melt as shown in Figure 2, P1 and P2 where the SSS still increases in NNS while the sea ice melts. Second, the ∼25 years delay between the variations in GSOC and AMOC emanates from a series of ocean-atmosphere processes (namely, ocean heat release-anomalously warmer air temperature-southerly increase in the Denmark Straitadvection of warmer Atlantic surface water) that link the NNS and NWNA regions and ultimately result in reduced sea ice extent around Greenland. Finally, we suggest that the monostable regime of our model simulation, and its negative feedback on the AMOC increase, provides a potential restoring force, assisted by the reduction of oceanic heat release from the NNS region as the temperature below 200 m decreases by about 2°C. These processes are summarized in Figure 12 and correspond relatively well to that described for the LOVECLIM model in Friedrich et al. (2010) even though the AMOC strength they simulated is more than 50% stronger than in our study and that the abrupt AMOC transition they described is on the opposite sign as the one described here (i.e., they showed reduction, whereas we show increase in AMOC strength). This suggests that for a certain range of AMOC state, the model is able to generate spontaneous and large variations in AMOC involving changes in both strength and structure of the overturning. This higher-instability regime arises suddenly under the influence of gradual forcing due to key nonlinearities in the response to this gradual forcing.
This study shows the importance of the sea ice cover in the high latitudes and how the ocean-atmosphere system can trigger short-lived climate transitions and link two regions of high importance for the overturning circulation. While the warmer (early) part of the LIG remains stable, it seems that there is a threshold when the sea ice has extended further southward to feedback on the atmospheric temperature and ocean heat content. We suggest that such mechanisms could have been involved during the late LIG, for instance around the 119 ka event reported in Galaasen et al. (2014), where no ice rafted detritus could be linked to the large variations in δ 13 C, which suggest the possibility that a large water mass reorganization occurred. These mechanisms could also be potential candidates for the shifts recorded during previous interglacials such as the MIS 11c and MIS 9e, where both occur while the IRD record is low . Further model analysis under previous interglacial boundary conditions are required to test the range of background conditions, especially in sea ice extent, capable of triggering such variability/instability.
In addition, setting these spontaneous AMOC variations into a broader context, we note that no significant changes in SSTs are simulated in the Atlantic section of the Southern Ocean (not shown here) when the AMOC switches from one operating mode to another. This potentially distinguishes the centennial AMOC oscillations simulated in this study from the multimillennial AMOC oscillations of glacial periods (Broecker, 1997;Toggweiler & Lea, 2010). Furthermore, the key regions in this model (NNS and NWNA) where sea ice, surface ocean fluxes, and water mass advection changes are crucial for setting the vertical stratification and the overturning strength are projected to be ice free as early as the second half of this century. This suggests that such mechanisms responsible for the increase in AMOC strength are unlikely to occur toward the end of this century.