A salty deep ocean as a prerequisite for glacial termination

Abstract

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.

Main

A multitude of ice-core and marine records show that the past 800,000 years of Earth’s history were dominated by global-scale alternations between glacial and interglacial conditions with a high degree of nonlinearity in the climatic response to gradual external forcing^1,2. This is particularly relevant for glacial terminations, which are thought to be driven by relatively weak insolation forcing due to changes in Earth’s orbital configuration^3,4. In this context, it has been proposed that disruptions to circulation patterns in the coupled atmosphere–ocean system might have acted to accelerate Antarctic warming and the outgassing of CO₂ from the Southern Ocean (SO) to cross a critical level necessary to maintain interglacial conditions at a global scale^5,6,7. A possible mechanism to reach such a point is an abrupt weakening of the Atlantic meridional overturning circulation (AMOC) and a corresponding reduction of northward oceanic heat transport associated with SO and Antarctic warming, according to the bipolar-seesaw concept^8,9,10, which has been detected in a wide spectrum of models^11,12. However, since abrupt oscillations between relatively weak and strong AMOC states and corresponding SO changes are thought to be a ubiquitous feature of millennial-scale climate variability throughout most glacial periods^13,14, it is open to question what facilitates such a strong temperature response during deglaciation (Fig. 1).

Fig. 1: Records of Antarctic temperature and the rate of temperature change.

a, Time series of Antarctic temperature anomaly⁵⁷ (AAT ΔT). b, Rate of change (dT/dt) calculated after a temporal smooth of 500 years (red curve) and 200 years (blue curve). Deglacial events display higher rates of change on average (see Methods for discussion).

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²⁰. 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.²¹ 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²² 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²³. 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.

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²⁶. Such salinity changes could be brought about by enhanced equatorward sea-ice export in the glacial SO²⁸ in combination with reduced melting of Antarctic ice as a consequence of colder NADW²⁰. Analyses of different boundary conditions in atmosphere–ocean general circulation model (AOGCM) equilibrium simulations suggest that especially the interplay of low CO₂, low obliquity²⁹ and large Northern Hemisphere ice sheets³⁰ provides favourable conditions for the formation of salty deep waters around Antarctica.

Using the fully coupled Earth system model (COSMOS)^29,31 (Methods), Zhang et al. have previously shown²⁹ that the application of glacial boundary conditions can lead to particularly salty deep waters in the SO region compared with pre-industrial conditions. In this study, we augment the glacial (LGM_CTRL) and pre-industrial (PID_CTRL) control experiments of Zhang et al.²⁹ to include an intermediate climate state (40ka_CTRL, Methods, Figs. 3 and 4 and Extended Data Figs. 2–9) to investigate any differences in the response to ‘terminal’ and non-terminal AMOC perturbations (which we describe later).

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 δ¹⁸O values includes the possibility that NADW formation might actually have been weaker and even shallower^32,33,34 than in our simulated LGM state³⁵).

Fig. 2: Interglacial and glacial temperature and salinity distribution and changes between them as zonal averages in the Atlantic sector.

a,b, Interglacial (pre-industrial, PID_CTRL) (a) and glacial (Last Glacial Maximum, LGM_CTRL) (b) temperature and salinity (practical salinity units (PSU), white contours) distribution. c,d, Glacial–interglacial difference (LGM_CTRL – PID_CTRL) for temperature (c) and salinity (d). All fields represent 100 yr averages. Plotted using Ocean Data View (http://odv.awi.de).

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²⁰ (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. ²⁰, 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:

$$\frac{{{\mathrm{d}}T}}{{{\mathrm{d}}z}} = - \frac{{{\mathrm{d}}\rho/{\mathrm{d}}z}}{{\alpha_{\mathrm{T}}\rho _0}}\left( {1 + \frac{{{\mathrm{d}}\rho_{\mathrm{S}}/{\mathrm{d}}z}}{{{\mathrm{d}}\rho_{\mathrm{T}}/{\mathrm{d}}z}}} \right)^{ - 1}$$

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 ρ₀ 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 ($\frac{{{\mathrm{d}}\rho_{\mathrm{S}}/{\mathrm{d}}z}}{{{\mathrm{d}}\rho_{\mathrm{T}}/{\mathrm{d}}z}}$) for a variety of background states (Fig. 3).

Fig. 3: Temperature changes with depth in dependence of haline to thermal density ratio and background climate.

Temperature gradients (dT/dz) along lines of constant density stratification (dρ/dz) for the individual climate states with reference conditions from the zonally averaged Atlantic at 30° S and 3,000 m depth. The pressure reference level is 3,000 dbar to account for the impact of pressure conditions on temperature and stability properties at the chosen reference level of 3,000 m depth. The circles show the position of the simulated interglacial (PID_CTRL, black), 40 ka (40ka_CTRL, green) and Last Glacial Maximum (LGM_CTRL, blue) states (compare Extended Data 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).

For any of the three different background states with a given density structure (dρ/dz, ρ₀ 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⁶ m³ 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 high-latitude Southern Hemisphere response to an AMOC weakening only under LGM conditions (Fig. 4d,e).

Fig. 4: Changes in circulation strength, southern high-latitude temperature and sea-ice area fraction in response to North Atlantic freshwater perturbation.

a, NADW export (Sv) at 30° S. b–d, Surface temperature (°C) and sea-ice area fraction changes (iso-lines) for interglacial (PID_PERT – PID_CTRL) (b), 40 ka (40ka_PERT – 40ka_CTRL) (c) and Last Glacial Maximum (LGM_PERT – LGM_CTRL) (d) conditions, averaged between model years 700 and 800. e, The averaged Antarctic surface temperature changes as a 100 yr running mean with the same colour code as in a. Panels b–d plotted using Ocean Data View (http://odv.awi.de).

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⁴⁵ 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⁴⁶. 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).

Fig. 5: Temperature anomalies and the difference between them.

a,b, Temperature anomalies by AMOC weakening for interglacial (PID_PERT − PID_CTRL) (a) and Last Glacial Maximum (LGM_PERT − LGM_CTRL) (b) conditions between model years 100 and 200 as a 100 yr mean, zonally averaged in the Atlantic sector. c, The difference between the glacial and interglacial anomalies (LGM − PID). Plotted using Ocean Data View (http://odv.awi.de).

In all experiments, the Atlantic basin temperature response below ~2 km is one of cooling (Fig. 5 and Extended Data Fig. 8) as the NADW cell weakens (Fig. 4a) and shoals (Extended Data Fig. 7) during the first ~200 years of hosing. Compared with PID_PERT, the more pronounced LGM_PERT cooling between 2.0 and 2.5 km (Fig. 5c) reflects the initially shallower NADW cell (Extended Data Fig. 3) as 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 and PID_PERT. Only in the glacial state does the AMOC perturbation result in atmospheric warming in the southernmost ocean region. When combined with atmospheric 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¹², 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⁴⁷. 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²⁹ as employed by Kageyama et al.⁴⁷.

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₂ variability^{48,49,50,51,52}. Furthermore, it has also been suggested that AMOC-induced CO₂ changes⁵³ might provide an internal feedback with the potential to promote transitions between different climate states⁵⁴. Abrupt changes in the mode of AMOC occurred throughout the last glacial period¹³ and were linked to temperature changes across Antarctica and (on occasion) variations in CO₂⁵⁵, 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² 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.

Methods

Rate of change in Antarctic temperature records

In 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.

Model description

We use the comprehensive fully coupled Earth system model COSMOS (ECHAM5-JSBACH-MPIOM) in this study. The atmospheric model ECHAM5⁵⁸, complemented by a land surface component with dynamic vegetation JSBACH⁵⁹, is used at T31 resolution (~3.75°), with 19 vertical hybrid sigma-pressure levels. The ocean model MPIOM⁶⁰, including sea-ice dynamics that are formulated using viscous-plastic rheology⁶¹, 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³¹, the Miocene warm climate^62,63, the Pliocene⁶⁴, the internal variability of the climate system⁶⁵, the Last Glacial Maximum climate^29,47 and glacial millennial-scale variability at intermediate climate conditions^54,66. The glacial (LGM_CTRL) and pre-industrial (PID_CTRL) states used in our study are experiments LGMW and PI from Zhang et al.²⁹. 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²⁰ 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.²⁹, 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)⁶⁷. The greenhouse gas concentrations are 195 ppm (CO₂), 413 ppb (CH₄) and 231 ppb (N₂O)^68,69. The ice-sheet configuration is a combination of ice-sheet reconstructions from ICE-5G⁷⁰ 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²⁹. 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 and 40ka_PERT. A freshwater flux of 0.15 Sv is added to the Ruddiman Belt (40°–55° N, 45°–20° W) for 800 years to ensure a quasi-equilibrated response in the SO. A freshwater perturbation of 0.15 Sv agrees with estimates of deglacial meltwater release during Heinrich Stadial 1 as assessed in different studies^71,72,73. Global mean sea-level height is corrected for the input of additional freshwater. Hence, coastlines remain the same during a freshwater perturbation experiment. Averages of 100 years are shown in our results.

Vertical temperature and salinity changes

On the basis of the linear approximation of the density equation and neglected pressure variations, we can calculate density changes with depth by:

$$\frac{{{\mathrm{d}}\rho }}{{{\mathrm{d}}z}} = \rho _0 \left( { - \alpha_{\mathrm{T}}\frac{{{\mathrm{d}}T}}{{{\mathrm{d}}z}}} \right. + \beta_{\mathrm{S}}\left. {\frac{{{\mathrm{d}}S}}{{{\mathrm{d}}z}}} \right)$$

(1)

where ρ₀ 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:

$$\frac{{{\mathrm{d}}\rho_{\mathrm{T}}}}{{{\mathrm{d}}z}} = - \rho _0 \alpha_{\mathrm{T}}\frac{{{\mathrm{d}}T}}{{{\mathrm{d}}z}}$$

(2)

and a salinity contribution equates to:

$$\frac{{{\mathrm{d}}\rho _{\mathrm{S}}}}{{{\mathrm{d}}z}} = \rho _0 \beta _{\mathrm{S}}\frac{{{\mathrm{d}}S}}{{{\mathrm{d}}z}}$$

(3)

Rearrangement of equation (2) leads to:

$$\frac{{{\mathrm{d}}T}}{{{\mathrm{d}}z}} = - \frac{{{\mathrm{d}}\rho_{\mathrm{T}}/{\mathrm{d}}z}}{{\alpha_{\mathrm{T}} \rho _0}}$$

(4)

Further substituting dρ_T/dz in equation (4) by:

$$\frac{{{\mathrm{d}}\rho _{\mathrm{T}}}}{{{\mathrm{d}}z}} = \frac{{{\mathrm{d}}\rho }}{{{\mathrm{d}}z}} - \frac{{{\mathrm{d}}\rho_{\mathrm{S}}}}{{{\mathrm{d}}z}}$$

(5)

temperature changes with depth can be expressed as a function of the ratio between the haline and thermal density changes with depth:

$$\frac{{{\mathrm{d}}T}}{{{\mathrm{d}}z}} = - \frac{{{\mathrm{d}}\rho/{\mathrm{d}}z}}{{\alpha_{\mathrm{T}}\rho _0}}\left( {1 + \frac{{{\mathrm{d}}\rho_{\mathrm{S}}/{\mathrm{d}}z}}{{{\mathrm{d}}\rho_{\mathrm{T}}/{\mathrm{d}}z}}} \right)^{ - 1}$$

(6)

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⁷⁴, which limits the ability to investigate a fuller spectrum of relevant processes and interactions such as the ‘pre-freshening’²⁰ 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⁷⁴, which represents a fundamental source of abyssal ocean circulation in the modern ocean. However, there is another mode associated with open-ocean convection⁷⁵ 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⁷⁵.

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⁷⁸ 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²⁰, 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)⁷⁹ to quantify changes in the subshelf melt rates. PICO is developed from the ocean box model of Olbers and Hellmer⁸⁰ and is implemented as a module of the Parallel Ice Sheet Model (PISM)⁸¹ 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.⁷⁹ to obtain a realistic pre-industrial ice sheet. The PISM settings were taken from Sutter et al.⁸², and the atmospheric forcing is from regridded yearly RACMO2.3 2 m air temperature and precipitation fields⁸³.

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⁻¹ (Extended Data Fig. 4a), which is not far from the observed value 1,500 ± 237 Gt yr⁻¹ for present-day conditions⁸⁴. 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 δ¹⁸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.²⁹. From the six other models in that study, two models (CCSM3⁸⁵ and HadCM3M2⁸⁶) 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²⁹. 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⁴⁰ generations, the north–south deep-ocean salinity contrast is dominated by a relatively salty North Atlantic, except within CCSM4⁸⁷. 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⁸⁸.

• 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²⁹. In particular, orbital configurations incorporating low obliquity and low atmospheric CO₂ 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⁴¹. 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²⁰ (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⁴². 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²⁶ 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⁸⁸. Similarly to the ‘glacial’ state in Galbraith and Lavergene³⁰, we simulate a lower salinity at this core location (36.04 PSU). Aside from this extreme salinity value, our simulated salinity in the deep Pacific between 3,500 m and 5,500 m is 36.03. This value is within the uncertainty of the reconstructed salinity (36.1 ± 0.1) based on chloride concentrations of deep Pacific bottom waters at five different core locations containing the same depth range²⁷.

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:

$$\frac{{{\mathrm{d}}\rho _{{\mathrm{T}}\_{\mathrm{LGM}}}}}{{{\mathrm{d}}z}} + \frac{{{\mathrm{d}}\rho _{{\mathrm{S}}\_{\mathrm{PID}}}}}{{{\mathrm{d}}z}} < 0$$

(7)

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:

$$\frac{{{\mathrm{d}}T_{\_{\mathrm{LGM}}}}}{{{\mathrm{d}}z}} > - \frac{{{\mathrm{d}}\rho /{\mathrm{d}}z}}{{\alpha _{{\mathrm{T}}\_{\mathrm{PID}}}\ \rho _{0\_{\mathrm{PID}}}}}$$

(8)

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 ρ₀__PID still requires a cooling with depth, which is stronger than dT__LGM/dz.