A salty deep ocean as a prerequisite for glacial termination

Deglacial transitions of the middle to late Pleistocene (terminations) are linked to gradual changes in insolation accompanied by abrupt shifts in ocean circulation. However, the reason these deglacial abrupt events are so special compared with their sub-glacial-maximum analogues, in particular with respect to the exaggerated warming observed across Antarctica, remains unclear. Here we show that an increase in the relative importance of salt versus temperature stratification in the glacial deep South Atlantic decreases the potential cooling effect of waters that may be upwelled in response to abrupt perturbations in ocean circulation, as compared with sub-glacial-maximum conditions. Using a comprehensive coupled atmosphere–ocean general circulation model, we then demonstrate that an increase in deep-ocean salinity stratification stabilizes relatively warm waters in the glacial deep ocean, which amplifies the high southern latitude surface ocean temperature response to an abrupt weakening of the Atlantic meridional overturning circulation during deglaciation. The mechanism can produce a doubling in the net rate of warming across Antarctica on a multicentennial timescale when starting from full glacial conditions (as compared with interglacial or subglacial conditions) and therefore helps to explain the large magnitude and rapidity of glacial terminations during the late Quaternary. Heat stored in the deep ocean due to salinity stratification contributed to rapid Antarctic warming during middle and late Pleistocene glacial terminations, according to coupled atmosphere–ocean general circulation model simulations.


Interglacial and glacial water-mass configurations
The large-scale meridional water-mass distribution in the modern ocean is characterized by deep waters originating in the North Atlantic (North Atlantic Deep Water; NADW) and the SO (Antarctic Bottom Water; AABW). These water masses are connected by mixing processes and wind-driven upwelling within the SO [15][16][17][18] . AABW represents the densest water mass in the modern ocean due to its low temperature (AABW is relatively fresh and cold compared with the overlying NADW) 18,19 . Only at great depth does the pressure dependence of the thermal expansion coefficient in the equation of state (the thermobaric effect) mean that the warmth of NADW makes it less dense than AABW 20 . Thus, a reduction in the temperature contrast between the two water masses will tend to destabilize the water column, ultimately resulting in NADW replacing AABW as the deepest water mass.
With respect to ocean temperature, the most accurate/robust evidence for glacial ocean cooling on a global scale stems from noble gases trapped in ice cores. Bereiter et al. 21 show that mean ocean cooling (2.5 ± 0.24 °C) was greater than the maximum possible cooling of AABW, since modern AABW has a temperature of −0.88 °C 22 and seawater has a freezing point of about −2 °C. Of particular relevance to our study is the implication that 'warmer' water masses such as NADW and intermediate water masses within the Southern Atlantic must have cooled to a greater extent than AABW 23 . This implies that thermal deep-ocean stratification in the glacial ocean was weaker than it is today. In a hypothetical scenario with no change to the salinity distribution, we can see that such heterogeneous cooling of the glacial ocean could lead to an unstable water column (Extended Data Fig. 1), and yet proxy reconstructions suggest that the glacial version of AABW maintained its position as the deepest water mass at least in the South Atlantic region 24,25 . Therefore, something probably acted to compensate for the loss of stability from thermal stratification in the glacial ocean.

A salty deep ocean as a prerequisite for glacial termination
Gregor Knorr 1 ✉ , Stephen Barker 2 , Xu Zhang 1,3,4 , Gerrit Lohmann 1,5 , Xun Gong 1,6 , Paul Gierz 1 , Christian Stepanek 1 and Lennert B. Stap 1,7 Deglacial transitions of the middle to late Pleistocene (terminations) are linked to gradual changes in insolation accompanied by abrupt shifts in ocean circulation. However, the reason these deglacial abrupt events are so special compared with their sub-glacial-maximum analogues, in particular with respect to the exaggerated warming observed across Antarctica, remains unclear. Here we show that an increase in the relative importance of salt versus temperature stratification in the glacial deep South Atlantic decreases the potential cooling effect of waters that may be upwelled in response to abrupt perturbations in ocean circulation, as compared with sub-glacial-maximum conditions. Using a comprehensive coupled atmosphere-ocean general circulation model, we then demonstrate that an increase in deep-ocean salinity stratification stabilizes relatively warm waters in the glacial deep ocean, which amplifies the high southern latitude surface ocean temperature response to an abrupt weakening of the Atlantic meridional overturning circulation during deglaciation. The mechanism can produce a doubling in the net rate of warming across Antarctica on a multicentennial timescale when starting from full glacial conditions (as compared with interglacial or subglacial conditions) and therefore helps to explain the large magnitude and rapidity of glacial terminations during the late Quaternary.
Indeed, other proxy evidence 26,27 suggests that the salinity of glacial AABW could have been much greater than today, even to the point that deep-ocean salinity stratification may have been reversed with respect to the modern 26 . Such salinity changes could be brought about by enhanced equatorward sea-ice export in the glacial SO 28 in combination with reduced melting of Antarctic ice as a consequence of colder NADW 20 . Analyses of different boundary conditions in atmosphere-ocean general circulation model (AOGCM) equilibrium simulations suggest that especially the interplay of low CO 2 , low obliquity 29 and large Northern Hemisphere ice sheets 30 provides favourable conditions for the formation of salty deep waters around Antarctica.
Using the fully coupled Earth system model (COSMOS) 29 In the glacial ocean simulation LGM CTRL , the strongest temperature decrease (of more than 4 °C with respect to both PID CTRL and 40ka CTRL ) occurs within levels of the southward NADW return flow ( Fig. 2 and Extended Data Figs. 2 and 3). The local LGM CTRL cooling maximum at ~2.5 km depth (Fig. 2c) can be linked to a combination of a shoaling (~200 m) of the glacial NADW cell and production of colder NADW (although we note that a model-data comparison for our glacial state with marine calcite δ 18 O values includes the possibility that NADW formation might actually have been weaker and even shallower 32-34 than in our simulated LGM state 35 ).
In stark contrast to the changes at intermediate depths, the weakest temperature and strongest salinity differences between the pre-industrial and glacial states are detected in the deep SO (Fig.  2c,d). Here we observe a shift from a purely temperature-controlled   stratification to one that is also characterized by a stable haline density gradient (at least within the SO). A similar contrast is observed between glacial and intermediate (40 thousand years ago (ka)) conditions (Extended Data Fig. 2). We note that the mechanism of deep-ocean salinification proposed by Adkins 20 (which calls on the cooling of NADW as a means for reducing the melting of Antarctic land ice and related input of freshwater to the SO) could work in tandem with increased sea-ice export (as observed in our model simulation) to produce highly saline bottom waters.
An overview of simulated salinity structure in LGM states (unperturbed) encompassing the last three Paleoclimate Model Intercomparison Project (PMIP)/Coupled Model Intercomparison Project (CMIP) model generations [36][37][38][39][40] shows that a simulation of a salty glacial deep ocean (Extended Data Fig. 10) by AABW on a global scale is a persistent challenge, which partly reflects limits of current models to simulate key processes 41,42 that influence AABW characteristics. This also includes the mechanism described by ref. 20 , and while the current generation of coupled AOGCMs cannot resolve that mechanism explicitly, we can at least demonstrate its viability using a coupled ice-sheet/ice-shelf cavity model (Methods and Extended Data Fig. 4). However, the precise mechanism by which a salty deep ocean is generated is not important for the following discussion, in which we investigate the theoretical implications of a relative increase in the importance of salinity versus temperature stratification in the glacial ocean. For further discussion, see Methods.

Salinity gradient and temperature changes with depth
Deep-ocean stratification depends on the interplay between salinity and temperature changes with depth. We therefore derive an analytical expression (which is independent of the model simulations) for the thermal gradient (dT/dz) that accounts for the interplay between thermal and haline contributions to the net density change with depth: Here, dρ S /dz and dρ T /dz represent the haline and thermally driven density changes with depth and dρ/dz is the combined (net) change in density with depth (which must remain positive to maintain stability). Density and the thermal expansion coefficient are indicated by ρ 0 and α T . We now explore the implications of a glacial ocean state in which an increase in salinity stratification accommodates a decrease in the thermal contrast between northern-and southern-sourced deep-water masses. In particular, we are interested in the surface SO temperature response that would result from the upwelling of deeper waters (for example, due to a perturbation of the AMOC) that are relatively warm (not as cold) compared with a non-glacial state.
With the aid of the different terms in the equation, we can evaluate the importance of the stable glacial salinity stratification for the reduced thermal stratification in LGM CTRL . We focus our attention on the deep South Atlantic (zonally averaged at 30° S) to assess the export characteristics out of the Atlantic basin and to avoid effects of intra-basin scale upwelling on Southern Hemisphere communication. In particular, the choice of the depth level influences the thermohaline characteristics of the reference water masses. Therefore, we choose a level at 3,000 m depth, which is within a broader core of waters with similar stratification characteristics that represent an unstable (PID CTRL , 40ka CTRL ) or stable (LGM CTRL ) salinity stratification. The first term on the right-hand side of the equation can then be diagnosed from our modelled control runs, enabling us to plot dT/dz versus the ratio of haline to thermally driven density changes ( dρ S /dz dρ T /dz ) for a variety of background states (Fig. 3).
For any of the three different background states with a given density structure (dρ/dz, ρ 0 and α T ), it can be seen that for constant dρ/dz the ratio of haline to thermally driven density stratification must increase to accommodate any decrease in the vertical temperature contrast between northern-and southern-sourced deep-water components in the South Atlantic (coloured curves in Fig. 3). The key difference between our glacial simulation and a pre-industrial or intermediate (40 ka) state is the large increase in salinity stratification (which is actually reversed with respect to the non-LGM states) relative to the thermal stratification. Together with the palaeo-evidence described earlier, this suggests that the cooling potential of upwelled water to the surface SO in a full glacial is relatively much less than for either the pre-industrial or intermediate state. Naturally, this has implications for the oceanic and atmospheric response to a perturbation of the AMOC that might be expected as a necessary ingredient of glacial termination.

Southern high-latitude response to AMoC weakening
We next perform perturbation simulations to investigate the implications of different glacial-interglacial water-mass configurations with respect to SO and Antarctic temperature responses to an abrupt weakening of the AMOC ( Fig. 4 and Extended Data Fig. 5). We use the same model set-up as in the control simulations and apply a freshwater perturbation of 0.15 Sv (1 Sv = 10 6 m 3 s -1 ) to the North Atlantic ice-rafted debris belt for 800 years (further details are provided in Methods). In summary, we find that warming sea surface conditions and declining sea-ice cover occur as the dominant LGM   Fig. 6). The curves show solutions for stable density stratification (dρ/dz > 0; that is, density increases with depth) with a stable thermal stratification (dρ T /dz > 0; that is, temperature decreases with depth). Hence, positive values on the x axis indicate a stable haline stratification (dρ S /dz > 0; that is, salinity increases with depth, pink shaded phase space) and negative values indicate an unstable haline stratification (dρ S /dz < 0; that is, salinity decreases with depth, grey shaded phase space).
high-latitude Southern Hemisphere response to an AMOC weakening only under LGM conditions (Fig. 4d,e).
In the following discussion, we focus on the comparison of the perturbation simulations LGM PERT and PID PERT (the 40ka PERT responses are similar to those of PID PERT ) (Fig. 4). On weakening of AMOC, surface warming at the southern boundary of the Atlantic basin (~30° S) penetrates to a depth of ~2 km irrespective of the background conditions ( Fig. 5a,b). This warming response is a common characteristic among NADW perturbation simulations 43,44 , which can be linked to heat storage 45 in the low latitudes and South Atlantic and reduced ventilation of intermediate and deep waters, combined with the downward mixing of heat from the thermocline 46 . This 'advection-diffusion balance' concept is also consistent with the relatively strong LGM PERT warming (Fig. 5) along the perturbed glacial AMOC-return pathway (Extended Data Fig. 7), reflecting the more pronounced vertical temperature contrast across the NADW cell in the unperturbed LGM CTRL state ( Fig. 2 and Extended Data Fig. 3).
In all experiments, the Atlantic basin temperature response below ~2 km is one of cooling ( the AMOC is perturbed (Extended Data Fig. 7). Below this depth, however, cooling in LGM PERT is weaker than in the PID PERT (Fig. 5c). At first glance, this is perplexing since shoaling of the overturning cell (Extended Data Fig. 7) and isopycnals at these depth levels is more pronounced in LGM PERT (the density levels at 3,000 m depth in PID CTRL and LGM CTRL shoal by ~160 m and ~230 m, not shown). However, the phenomenon is a direct consequence of the build-up of relatively warmer waters in the glacial deep ocean (with respect to intermediate depths) permitted by the increase in salt stratification ( Fig. 3 and Extended Data Fig. 6). As predicted, this fundamental change in deep-ocean stratification (from temperature controlled in PID CTRL to one that is dependent on a stable haline density gradient in LGM CTRL ) leads to a very different SO surface temperature response to North Atlantic freshwater forcing and AMOC perturbation. In PID PERT , cooling is observed across the surface SO due to the upwelling of cold deep waters when the AMOC is in a perturbed state. By contrast, for LGM PERT between ~100 and 200 years we observe a continuous surface warming across Antarctic circumpolar latitudes south of ~45° S and extending into the Drake Passage (south of ~55° S), reflecting the upwelling of relatively warmer waters that characterize the glacial deep ocean in our simulations ( Fig. 5 and Extended Data Figs. 8 and 9). It is important to emphasize that cooling does occur in the deep ocean during LGM PERT (Fig. 5b), but this cooling is reduced relative to PID PERT (Fig. 5c) and, in combination with warming at intermediate depths, produces a net warming at the SO surface.
The contrasting SO surface temperature responses lead to very different temperature responses over Antarctica (Fig. 4e). In fact, the integrated surface heat flux anomaly south of 55° S is of opposite signs in LGM PERT  heat transport, the positive ocean-atmosphere heat flux observed in LGM PERT produces by far the largest and fastest temperature rise across Antarctica within the first ~200 years when compared with PID PERT and 40ka PERT (Fig. 4e). After ~200-300 years, NADW strength and geometry have largely adjusted to the persistent freshwater perturbation, and Antarctic warming proceeds at a reduced rate.

Salty glacial deep ocean and deglaciation
We have demonstrated that a fundamental change in deep-ocean stratification can lead to a greatly enhanced climatic response to AMOC weakening while in a glacial state. The simulated Antarctic seesaw response is amplified by a factor of ~2 on a multicentennial timescale compared with interglacial or subglacial conditions. The surface cooling (as opposed to warming) that we observe in the SO high latitudes under pre-industrial and subglacial conditions might be a model-specific characteristic 12 , but it is important to note that the key to an amplified seesaw warming is a more pronounced haline to thermal density stratification ratio (which is predicted from proxy reconstructions, as described earlier). Increasing this ratio has the same effect irrespective of the climate background conditions, as shown in the phase space analysis in Fig. 3. The underlying mechanism might also therefore represent a key to understanding the spread of SO high-latitude responses in previous glacial AMOC perturbation experiments using a range of models 47 . For the specific model used in our study, a damped SO surface temperature response is obtained for LGM experiments using a modern-like salinity stratification 29 as employed by Kageyama et al. 47 . Our simulations do not include a representation of the global carbon cycle, but several previous studies have called on AMOC-induced stratification changes within the SO as a key player in glacial-interglacial CO 2 variability 48-52 . Furthermore, it has also been suggested that AMOC-induced CO 2 changes 53 might provide an internal feedback with the potential to promote transitions between different climate states 54 . Abrupt changes in the mode of AMOC occurred throughout the last glacial period 13 and were linked to temperature changes across Antarctica and (on occasion) variations in CO 2

55
, and yet only those changes associated with the last termination were sufficient to mark the transition back to an interglacial state. Previous studies have suggested various combinations of ice-sheet size 2 and specific orbital configurations 1,56 in an attempt to explain why glacial termination occurs only when it does. Our findings provide an alternative explanation: that global-scale water-mass salinity characteristics represent a key player to promote deglaciation once a salty deep-ocean reservoir has been formed at full glacial conditions.

online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41561-021-00857-3. Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/. © The Author(s) 2021 Fig. 1, the relatively large rates of change observed in the upper sections of the record (particularly the Holocene, which has not experienced high-amplitude millennial-scale variability) using a 200 yr pre-smooth reflects the much higher temporal resolution of the upper sections of ice cores relative to deeper sections. Using a pre-smoothing of 500 yr alleviates this problem while (in our opinion) preserving an accurate (not overly smoothed) representation of the relative rates of change associated with millennial-scale events.

Rate of change in Antarctic temperature records. In
Model description. We use the comprehensive fully coupled Earth system model COSMOS (ECHAM5-JSBACH-MPIOM) in this study. The atmospheric model ECHAM5 58 , complemented by a land surface component with dynamic vegetation JSBACH 59 , is used at T31 resolution (~3.75°), with 19 vertical hybrid sigma-pressure levels. The ocean model MPIOM 60 , including sea-ice dynamics that are formulated using viscous-plastic rheology 61 , runs on a bipolar curvilinear Arakawa-C grid. The grid size (120 × 101) corresponds to a formal longitudinal resolution of ~3° (360/120) and a latitudinal resolution of ~1.8° (180/101). The vertical domain is represented by 40 unevenly spaced depth layers. The climate model has already been used to simulate, for example, the past millennium 31 29 . The glacial state is characterized by a salty deep-ocean reservoir in reasonable agreement with available salinity reconstructions. Note that the generation of an interglacial-to-glacial salinity reversal is expected to be dependent on, for example, the model itself, the model configuration and the experimental set-up (Extended Data Fig. 10). However, for the mechanism described here, a salinity reversal per se is not necessary. The data-based demand that the temperature gradient between northern and southern sources of deep water within the Atlantic was reduced during the LGM compared with today 20 implies that a reduction in the thermal deep-ocean stratification was accompanied by a stronger glacial salinity increase of AABW relative to NADW. Were this not the case, then the density stratification of the glacial South Atlantic/SO would necessarily have been weaker than today. A simple evaluation for the importance of a glacial salinity reversal for glacial-interglacial density changes is provided in the corresponding section.
Experimental design. Besides the pre-industrial (PID CTRL ) and the Last Glacial Maximum (LGM CTRL ) control climates from Zhang et al. 29 , we conduct an equilibrated MIS3 control simulation (40ka CTRL ) by imposing fixed boundary conditions for 40 ka. In detail, the three orbital parameters in experiment 40ka CTRL are 0.013146 (eccentricity), 358.17° (precession of the equinoxes, the angle between the Earth's position during the Northern Hemisphere vernal equinox and the orbit perihelion) and 23.61° (obliquity) 67 . The greenhouse gas concentrations are 195 ppm (CO 2 ), 413 ppb (CH 4 ) and 231 ppb (N 2 O) 68,69 . The ice-sheet configuration is a combination of ice-sheet reconstructions from ICE-5G 70 and PMIP3. That is, the topography anomaly between 40 ka and the LGM in ICE-5G is first calculated and then is added to the PMIP3 LGM topography 29 . Global mean sea level is ~80 m lower at 40 ka than the pre-industrial level. Experiment 40ka CTRL is integrated for 5,000 years to an equilibrated state that serves as the starting point for the perturbation experiment.
Identical hosing experiments are then conducted under the three different climate backgrounds in experiments PID PERT , LGM PERT where ρ 0 is density, T is temperature, S is salinity, α T is the thermal expansion coefficient and β S is the haline contraction coefficient. So that a temperature contribution to density changes with depth equates to: and a salinity contribution equates to: Rearrangement of equation (2) leads to: Further substituting dρ T /dz in equation (4) by: temperature changes with depth can be expressed as a function of the ratio between the haline and thermal density changes with depth: Formation of a salty glacial deep ocean: model limitations. State-of-the-art AOGCMs lack fully interactive ice-sheet components and do not allow to simulate increased glacial melt in response to ocean warming. Furthermore, essential processes such as ice shelf-ice cavity ocean interactions are not included in the current generation of fully coupled AOGCMs 74 , which limits the ability to investigate a fuller spectrum of relevant processes and interactions such as the 'pre-freshening' 20 of AABW via NADW. Furthermore, the ocean model component in our study and other global climate models is too coarse to represent the northward transport of dense shelf waters across the Antarctic shelf edge and their subsequent sinking to depth 74 , which represents a fundamental source of abyssal ocean circulation in the modern ocean. However, there is another mode associated with open-ocean convection 75 that has been identified as an additional surface to deep-ocean exchange process based on satellite data 76,77 . This modern mode with notable deep (>2,000 m) open-ocean convection between 90° S and 55° S under pre-industrial conditions has been diagnosed in more than two-thirds of the CMIP5 models. While this polynya-induced mode is not dominant today, it could have been the main source of AABW formation during full glacial conditions with a widely ice-covered continental shelf and a lower sea level than today. This situation might have favoured a glacial open-ocean type convection mode that is free to reach the seafloor 75 .
In that sense, a dominant glacial production regime of AABW by open-ocean convection, predominantly in the centres of the Weddell and Ross Sea gyres, as modelled in our model set-up 78 might capture the actual glacial water-mass conversion much better than under modern conditions.
To test the viability of the mechanism for the generation of a salty glacial deep ocean as proposed by Adkins 20 , we have evaluated how far our simulated control states are consistent with less melting of Antarctic ice by colder NADW during glacial conditions. In this context, we have applied the Potsdam Ice-Shelf Cavity Model (PICO) 79 to quantify changes in the subshelf melt rates. PICO is developed from the ocean box model of Olbers and Hellmer 80 and is implemented as a module of the Parallel Ice Sheet Model (PISM) 81 to simulate the vertical overturning circulation in ice-shelf cavities and thus enables the computation of subshelf melt rates (Extended Data Fig. 4) consistent with this circulation. We run the simulations at a spatial resolution of 16 × 16 km. The PICO settings are the standard settings that were used by Reese et al. 79 to obtain a realistic pre-industrial ice sheet. The PISM settings were taken from Sutter et al. 82 , and the atmospheric forcing is from regridded yearly RACMO2.3 2 m air temperature and precipitation fields 83 .
Driving this model set-up with the thermohaline ocean fields of experiment PID CTRL in equilibrium simulations of 200 kyr we end up with an aggregated melt flux of ~1,400 Gt yr −1 (Extended Data Fig. 4a), which is not far from the observed value 1,500 ± 237 Gt yr −1 for present-day conditions 84 . Keeping the same experimental set-up but replacing the thermohaline fields by the 40ka CTRL and LGM CTRL fields, results in a reduction of the aggregated melt flux by ~26% and 50%, respectively. This suggests that a reduction of Antarctic ice melt by colder NADW could indeed contribute to the generation and maintenance of a salty deep-ocean reservoir, especially since surface melt rates are negligible in the contribution to Antarctic ice-mass loss for modern climate conditions and colder climates. Hence, attempts should be made to incorporate this mechanism in future Earth system model approaches to explore essential interactions in particular at the ice-ocean interface, and it is open to question how far a fuller representation of the AABW formation mechanism would lead to a different model response. These investigations should be combined with the analysis of, for example, benthic foraminiferal δ 18 O and Mg/Ca records to discriminate sea-ice-derived versus Antarctic-ice-derived salinity stratification.
Glacial deep-ocean salinity structure in models of different PMIP/CMIP generations. Currently, there is no collective framework of perturbation experiments for glacial and interglacial conditions that would allow us to test our perturbation experiments within a model intercomparison. Nevertheless, we have generated an overview of simulated glacial deep-ocean salinity structure in LGM states (unperturbed) encompassing the last three PMIP/CMIP model generations (Extended Data Fig. 10).
Using the same model as employed in our study, the glacial deep-ocean salinity structure (including a PMIP2/CMIP4 36,37 model intercomparison) was investigated in the study of Zhang et al. 29 . From the six other models in that study, two models (CCSM3 85 and HadCM3M2 86 ) also generated a stable glacial salt stratification (Extended Data Fig. 10). Note that PMIP utilized no specific protocol concerning the initial ocean condition for LGM simulations 36,37 , and notably the three models (including COSMOS) that were initialized from a previous glacial ocean state were the only ones to simulate a stable glacial deep-ocean salt stratification 29 . We suggest that this difference is a key to the different outcomes.
In model simulations at hand for the PMIP3/CMIP5 38,39 and PMIP4/CMIP6 40 generations, the north-south deep-ocean salinity contrast is dominated by a relatively salty North Atlantic, except within CCSM4 87 . In that model, saltier deep SO conditions are simulated, but these salinity characteristics do not extend into the Atlantic basin. Hence, the simulation of a salty glacial deep ocean via AABW at a global scale is not a coherent characteristic in LGM simulations with current coupled AOGCMs.
Where do we stand now? We think at least three possibilities should be taken into account for further research: • The inference of a salinity-stratified glacial ocean via AABW has been questioned, and the scarce pore-water chlorinity data availability necessitates an increase in data constraints 88 . • The simulation of a salty deep ocean seems to depend on conditions before the LGM, which therefore need to be incorporated into transient simulations rather than initiating models with LGM boundary conditions 29 . In particular, orbital configurations incorporating low obliquity and low atmospheric CO 2 concentrations in combination with Northern Hemisphere ice sheets in a later stage of a glacial cycle (but before the glacial maxima) should be considered as an alternative framework for simulating a salty glacial deep ocean 29,30 . • Models need to be improved; accurate formation of AABW via shelf processes is an ongoing challenge 41 . Furthermore, important processes such as interactive ice sheets and shelf-ice cavity ocean interactions need to be included for detailed studies, and we note that the mechanism of deep-ocean salinification proposed by Adkins 20 (Extended Data Fig. 4) could work in tandem with increased sea-ice export to produce highly saline bottom waters. Recently, it has also been suggested that the export of icebergs from the SO to the Atlantic basin might modify interglacial-to-glacial salt balance between NADW and AABW 42 . As yet, these key processes and mechanisms are missing in coupled climate models, and therefore the conjecture that models that are simulating the salty reservoir might do it for the wrong reason cannot be excluded.
Glacial deep-ocean salinity constraints. The observation/reconstruction of extremely high glacial SO salinities and abyssal waters is very scarce, and a salinity value of 37.1 ± 0.2 PSU at Shona Rise 26 represents a singularity that has been questioned by the analyses of pore fluid chlorinity/salinity data from deep-sea cores using estimation methods deriving from linear control theory 88 . Similarly to the 'glacial' state in Galbraith and  Estimates for the importance of salinity changes. To estimate the importance of interglacial-to-glacial salinity changes for the maintenance of deep-ocean stratification, we suggest two simple criteria. One criterion can help to evaluate the necessity of salinity changes in the evolution of the simulated glacial temperature structure. For example, in the PID CTRL and LGM CTRL experiment pair (Fig. 2), an unstable glacial situation applies if: That is, the glacial temperature decrease with depth is too weak to balance the effect of pre-industrial salinity decrease with depth on density (compare Extended Data Fig. 1).
While the latter criterion can help to assess the necessity of salinity changes, it does not answer the importance of a salinity reversal and the transition from an unstable to a stable salinity structure. Therefore, we additionally suggest to compare the LGM CTRL with a PID CTRL temperature gradient that would result from the shift of the unstable to a neutral salt stratification. This can be depicted in Fig.  3. Here the glacial cooling with depth (blue dot) is weaker than the intercept of the PID CTRL solution with the y axis. The corresponding criterion can be expressed as: That is, even if the unstable haline stratification contribution to dρ_ PID /dz would be replaced by a neutral haline density stratification (dρ S _ PID /dz = 0), the reduced thermal stratification along solutions with constant dρ_ PID /dz, α T_PID and ρ 0 _ PID still requires a cooling with depth, which is stronger than dT_ LGM /dz.

Data availability
The ice-core data used in this paper are available via the referenced sources.  conditions, averaged between model years 700-800. The north Atlantic cooling is most pronounced in the LGM perturbation experiment, which is linked to relatively pronounced sea ice cover increase and associated feedbacks. Together with the southern high latitude responses this is indicating a more pronounced interhemispheric seesaw response in surface temperature in the LGM perturbation experiment than under 40 ka or pre-industrial conditions. Extended Data Fig. 9 | Meridional overturning circulation in the Southern ocean. Shown are the Eulerian-mean meridional overturning stream function (Sv, 10 6 m 3 s −1 ) in depth coordinates in the perturbation experiments representing the average between model years 100-200 for a) interglacial (PID PERT ), b) intermediate (40ka PERT ) and c) glacial maximum (LGM PERT ) conditions (cf. Extended Data Fig. 6). Positive (negative) values indicate clockwise (anticlockwise) circulation. The latitudinal variation in zonal wind-stress induces upwelling between the easterly and westerly wind maxima, and downwelling outside this corridor. Upwelled water can be returned as part of the wind-driven cells south or north of the Antarctic divergence (Ant. Div.) that is linked to a southbound Ekman transport under the influence of Polar easterlies and a northbound Ekman transport under the influence of the westerlies. Along these transport pathways heat and freshwater fluxes at the surface are part of the generalized circulation pattern that feed into the lower or upper limb of the recirculation, respectively. However, the actual upwelling transport largely occurs along surfaces of equal density, with wind and eddy processes playing a central role 16 . To illustrate such an isopycnal upwelling pathway, we additionally show the σ 3 density surface (black line) that characterizes conditions at 3000 m depth and 30°S and upwells towards the latitude belt of the Antarctic Divergence that is subject to strong mixing during austral fall and winter seasons 78 .