In this study, we evaluate the simulated properties of the NADW and AABW in climate model simulations against glacial reconstructions through direct model-data comparison. We also explore the influences of an improved AABW in the models on the AMOC geometry through a set of dedicated numerical experiments. As changes in the density of deep waters potentially affect the pathways via which deep water returns to the surface, this set of experiments permits to assess the dependence of AMOC strength on various overturning pathways across climate states, eventually allowing to constrain AMOC strength with observation-based estimates of the pathways in both the LGM and warming climate.

Changes in the properties of the NADW and AABW influence the AMOC depth and strength and therefore became a central topic of paleoclimate modelling. The potential density of seawater, which is a function of temperature, salinity and pressure, is a direct indicator for the properties of the NADW and AABW and hence could provide information about their influences on the AMOC. We first discuss the mean temperature and salinity, and the climate-model biases relative to observations at Feni Drift, Bermuda Rise and Shona Rise (Fig. 1a) in pre-industrial simulations with the Kiel Climate Model (KCM, "Methods") and climate models participating in the Paleoclimate Modelling Intercomparison Project (PMIP). Temperature and salinity at these locations reveal the approximate properties of the NADW and AABW. Two pre-industrial simulations were conducted with the KCM: one (PI) is performed with the standard version and the other (PI) with salinity in the deep Southern Ocean (south of 60°S and below 3250 m) relaxed to the Levitus' salinity climatology ("Methods", Supplementary Fig. 2). All simulations (KCM and PMIP) presented in this study are integrated for at least a few hundred years to reduce the effects of spin-up.

There are noticeable warm biases in the pre-industrial control integrations of a number of climate models amounting up to about 3 °C in the temperatures at Feni Drift, Bermuda Rise and Shona Rise relative to the EN4.2.2 and WOA18 datasets (Fig. 1b-d). Yet the climatological mean salinity at the three sites is generally in good agreement with the observations (Fig. 1b-d). Overall, most PMIP models and the KCM reproduce reasonably well the potential density σ (contours in Fig. 1b-d), which is referenced to 2000 m depth, in the deep North Atlantic and Southern Ocean in the pre-industrial control integrations.

Following the PMIP protocol, the LGM simulations with the KCM and the PMIP models are performed by implementing glacial boundary conditions including among others continental ice sheets, greenhouse gas concentrations and orbital parameters (Methods). The potential density σ in the LGM simulations is significantly underestimated compared to the proxy data at the three sites (Fig. 1e-g). At Feni Drift and Bermuda Rise, the bias in σ can be mainly attributed to a warm temperature and to a lesser extent fresh salinity bias. The biased potential density at Shona Rise is largely related to fresh salinity biases in the deep Southern Ocean. Fresh biases there exceed 1 practical salinity unit (psu) in the LGM simulations performed with the standard KCM (Methods) and the PMIP models. We note that the salinity bias at Shona Rise in FGOLAS-g2 and MIROC-ESM exceeds 2 psu. This is likely because these models do not fully implement the PMIP protocol and do not add 1 psu to the oceanic initial conditions of the LGM simulations. The two CESM2 models show relatively dense waters at Shona Rise, which may be due to the parameterisation of density-driven overflow that moves dense water from the shelves to the deep ocean without cascading.

The σ of the NADW and AABW at their source regions is suggested to be a major factor controlling the AMOC in the LGM and modern climate. We examine the relationship of the AMOC in the region 20°S-10°N to the σ of the NADW and AABW. Here the AMOC depth is defined as the depth at which the overturning streamfunction is zero. The density of the NADW and AABW is represented by the vertical and area-weighted mean from 800 m to 2800 m over 50°N-55°N and 50°W-10°W in the North Atlantic, and below 3000 m over 50°S-55°S and 50°W-10°E in the South Atlantic, respectively (yellow rectangles in Fig. 1a). The AMOC depth exhibits a statistically significant correlation of 0.70 at the 95% confidence level with the σ contrast between NADW and AABW in the LGM simulations (Fig. 2a). The significance of correlation coefficient is tested using the p value approach. It is primarily the σ of the AABW that exhibits a strong relationship with the AMOC depth, with a correlation of -0.71 (Fig. 2b). The NADW hardly influences the AMOC depth, as illustrated by its small and insignificant correlation with σ (now shown). Clearly, the shoaling of the AMOC in the LGM relative to the pre-industrial simulations is linked to the enhanced stratification due to the rather dense AABW in the deep ocean.

The AMOC is projected to shoal under enhanced greenhouse-gas forcing in the models participating in the Coupled Model Intercomparison Project phase 6 (CMIP6). We examine the influence of the AABW density on the AMOC geometry in 90 global warming simulations under the Shared Socioeconomic Pathway 5-8.5 (SSP5-8.5) from 25 CMIP6 ensembles (Table 1, Supplementary Fig. 3). The individual members in each ensemble start from different atmospheric and oceanic initial conditions covering a wide range of climate states. We compute from all available ensemble members the mean AMOC depth between 20°S-10°N and the σ of the NADW and AABW averaged over the years 2071-2100. Similar to the LGM, the density contrast between the NADW and AABW exhibits a statistically significant correlation with the AMOC depth under global warming, which, however, is smaller amounting to 0.46 (Fig. 2c). Again, the σ of the AABW is the dominating factor determining the AMOC depth with a correlation of -0.72 (Fig. 2d), while the correlation with the σ of the NADW is again small and insignificant (not shown).

We perform 3 sensitivity experiments with the KCM in order to further understand the role of the AABW density in the LGM. As in many other climate models, the AABW in the KCM is formed by open-ocean convection and not over the continental shelf. Previous studies suggest that reducing warm temperature biases in the deep Southern Ocean improves the simulated ocean stratification and the AMOC. We correct the σ bias in the deep Southern Ocean by relaxing the KCM's salinity below 3250 m south of 60°S towards the palaeoceanographic salinity reconstructions (Methods). This approach has been proven useful to reproduce the vertical stratification of the glacial ocean. The initial conditions for the sensitivity experiments are the temperature and salinity climatology calculated from the LGM simulation with the standard KCM, to which either 1.0 psu, 1.5 psu or 2.0 psu is added (Methods). Each salinity-restoring experiment is integrated for 1600 years, with the monthly output from the last 300 years analysed here (Fig. 3a, b).

The σ biases at Shona Rise are relatively small in the salinity-restoring experiment, and the AMOC shoals considerably (Fig. 3b). In particular, the temperature and salinity in LGM lie within the uncertainty of the reconstructions (Fig. 1c). The AMOC shoals by about 800 m relative to PI to a depth of 2050 m in LGM (Fig. 3c, d), which agrees well with various proxy data. The AMOC strength reduces from 16.8 Sv in the standard LGM simulation to 13.8 Sv in the LGM experiment, which is still slightly larger with respect to the pre-industrial state (12.7 Sv). Moreover, the dominant role of the AABW on the AMOC depth is also found in the salinity-restoring experiments (Fig. 2b), whereas the density contrast between the NADW and AABW appears to be slightly less important (Fig. 2a). Overall, the restoring experiments support the argument that strong ocean stratification drives the shoaling of the AMOC during the LGM by limiting the deep convection in the North Atlantic.

The KCM suffers from systematic model biases in the North Atlantic as other PMIP models (Supplementary Fig. 4). To alleviate the impact of systematic biases, we compute the simulated SST responses in the LGM simulations relative to the pre-industrial experiment. The KCM's response of sea surface temperature (SST) over the North Atlantic benefits from the reduced σ bias in the Southern Ocean. The late-winter mean (January-March, JFM) SST changes during the LGM from MARGO, which is a proxy-based reconstruction for the glacial ocean conditions, depict a cooling pattern relative to the modern climate over the eastern North Atlantic, which extends from 60°N to 40°N (Fig. 4a). In comparison to MARGO, the standard KCM simulates a weaker response of cooling in the eastern North Atlantic that is restricted to north of 50°N from 35°W to 5°W (Fig. 4b). Correcting the salinity bias in the deep Southern Ocean in the experiment LGM, not only improves the AMOC but also broadens the meridional extent of the cooling in the eastern North Atlantic (Fig. 4c). Relative to MARGO, the root-mean-square-error averaged over 40°N-60°N and 35°W-5°W reduces from 5.32 °C to 3.87 °C in the restoring experiment. Consistently, the restoring experiment exhibits improvements in the simulated JFM mean SST compared to MARGO (Supplementary Fig. 5). We note that there are still some significant discrepancies in the SST-cooling response between MARGO and the salinity-restoring experiment, which may partly arise from the mean-state biases over the North Atlantic associated with deficiencies in the dynamical processes in the pre-industrial simulation with the KCM. There is an eastward expansion of the sea ice associated with the cooling over the North Atlantic in LGM (Fig. 4d-f), which reduces the local surface-heat flux and wind-driven oceanic salt transport. As a result, the deep oceanic convection, as measured by the JFM mean mixed layer depth, is weakened south of Iceland in LGM relative to the standard LGM simulation (Fig. 4g-i). We note that LGM and LGM show similar but smaller changes as LGM.

One remaining issue is the strength of the glacial AMOC. As proxy data disagree on the changes in the AMOC strength during the LGM, we seek to constrain the AMOC changes using present-day observations, which is similar to the concept of emergent constraints. This approach assumes that despite major differences among the climate models, relationships between the present-day and glacial climate are implicit within the model solutions of the partial differential equations governing the physical processes and parameterisations. Thus, the relationship between the present-day constraint variable () and glacial change in AMOC () can be expressed as where represents a negligible deviation from function . The function f can be obtained through analysis of a diverse model ensemble with varying resolution and parameterisations, such as the PMIP/CMIP ensembles. This allows us to establish a relationship between the glacial AMOC changes and the pre-industrial overturning pathways, where the pathways are quantified by the overturning streamfunction, which stands for the zonally integrated volume transport that returns the NADW to the surface (Methods).

In latitude-depth space, the AMOC can be regarded as a superposition of the shallower adiabatic cell and the deeper diabatic cell separated in the Atlantic and the Indo-Pacific basin. An acceleration (deceleration) of either cell strengthens (weakens) the AMOC. As suggested by Cessi (2019), about half of the volume of the modern NADW returns adiabatically to the surface through the wind-driven upwelling pathway in the Southern Ocean (SO_Wind). The remaining half of the NADW completes its circulation diabetically via diffusive upwelling in the Atlantic (At) and Indo-Pacific (IP). There is compelling evidence from modern hydrography that a geostrophic inter-basin transport connects the Atlantic with the Indo-Pacific, resulting from differences in the isopycnal depth on the order of 100 m between the two basins. Specifically, analysis of water mass properties shows that the NADW rising to the surface of the Southern Ocean is transformed into AABW and then transported to the bottom of the Indo-Pacific. These dense waters are upwelled through diapycnal mixing and eventually return to the Atlantic.

Each overturning pathway exhibits a large spread among the PMIP and CMIP models. Motivated by the influence of deep waters on the AMOC, we investigate relationships of the overturning pathways with the potential density of the NADW and AABW (Fig. 5). The Southern Ocean wind-driven pathway is positively correlated to the potential density contrast between the NADW and AABW in the LGM simulations (Fig. 5a). The correlation coefficient amounts to 0.62, which is significant at the 95% confidence level. Consistent with the LGM simulations, the potential density contrast between the NADW and AABW is associated with the magnitude of the Southern Ocean wind-driven pathway under global warming in CMIP simulations, with a correlation of 0.42 that is significant at the 95% confidence level (Fig. 5c).

The Indo-Pacific diffusive pathway exhibits no significant dependence on the potential density of the NADW and AABW in the LGM simulations with the PMIP models (blue lines Fig. 5b). Salinity-restoring experiments with the KCM suggest that enhanced potential density in the Atlantic basin tends to weaken the Indo-Pacific pathway (green lines in Fig. 5b). In contrast, global warming simulations with the CMIP6 models exhibits a positive correlation of the Indo-Pacific diffusive pathway with the potential density of the AABW, amounting to 0.53 which is significant at the 99% confidence level (Fig. 5d). Differences in the influence of the AABW on the Indo-Pacific diffusive pathway across the glacial and warming climate states are likely associated with the boundary conditions, particularly in the Northern Hemisphere. In summary, the two overturning pathways exhibit a significant relationship to the potential density of the AABW and NADW in the models across climate states.

Then each of the two pathways is examined with respect to the response of the AMOC strength, defined as the change in the LGM simulation relative to the pre-industrial simulation. Models with strong pre-industrial Southern Ocean wind-driven pathway tend to exhibit strengthened glacial AMOC driven by intensified westerly winds (Supplementary Fig 6), as indicated by the positive correlation (r = 0.58), explaining about one third of the total variance among the models (PMIP and KCM, Fig. 6a). PMIP models with a stronger pre-industrial Indo-Pacific diffusive pathway are sensitive to the reduced inter-basin transport, which is affected by the glacial ocean stratification and depict an opposite relationship with the glacial AMOC response (r = -0.54; PMIP and KCM, Fig. 6b). Furthermore, the Indo-Pacific diffusive pathway is highly correlated with the glacial AMOC response (r = 0.8, PMIP and KCM; Supplementary Fig 7c), consistent with previous studies.

The above results suggest that the Southern Ocean wind-driven and the Indo-Pacific diffusive pathway can be used as a constraint on the glacial AMOC response. This hypothesis is tested using 110 global warming simulations from 12 CMIP5 (Table 2, Supplementary Fig 8) and 25 CMIP6 ensembles (Table 1, Supplementary Fig 3) employing the SSP5-8.5 scenario. The ensemble mean is computed using all available members and the model data are averaged over the years 2071-2100. Consistent with the LGM simulations, the weakening of the AMOC under global warming is dependent on the magnitude of the Southern Ocean wind-driven and Indo-Pacific diffusive pathways in the pre-industrial simulations (Fig. 6c, d). In particular, the Indo-Pacific diffusive pathway plays a dominant role in the AMOC response to global warming, explaining about 50% of the variance among the CMIP models (Fig. 6d). The weakening of the AMOC under global warming is largely attributed to the restricted Indian Ocean inflow induced by the shoaling of the AMOC and compensated by the Southern Ocean pathway driven by the intensified westerly winds (Supplementary Fig 6). We note that the CMIP6 models exhibit a stronger dependence of the AMOC response on the pre-industrial pathways than the CMIP5 models. This could be due to improvements in ocean stratification and abyssal circulation in the CMIP6 models.

Since the responses of glacial AMOC in the models exhibits a dependence on the pre-industrial Southern Ocean wind-driven as well as Indo-Pacific diffusive pathway, we use this dependence, as represented by linear regression, to estimate the glacial AMOC strength (Fig. 6a, b). Two global ocean-reanalysis systems employing state-of-the-art data assimilation are used to compute the present-day constraining variables: GloRanV14 and ECCOv4. As shown in Fig. 6a, the dependence of the glacial AMOC response () on the pre-industrial Southern Ocean wind-driven pathway () can be simplified as using linear regression (green lines in Fig. 6a). The transport by the Southern Ocean wind-driven pathway amounts to 7.6 Sv (GloRanV14) and 14.0 Sv (ECCOv4). Substituting estimates for the Southern Ocean wind-driven pathway in the above constraint relationship, we obtain positive changes of the glacial AMOC amounting to 3.2 Sv (GloRanV14) and 8.6 Sv (ECCOv4). Both reanalysis datasets suggest a stronger glacial AMOC relative to the pre-industrial state (Fig. 6a). Similarly, regression analysis in Fig. 6b shows the constraint relationship between changes in the glacial AMOC and the present-day Indo-Pacific diffusive pathway as (green lines in Fig. 6b). The transport by the Indo-Pacific diffusive pathway in GloRanV14 amounts to 8.1 Sv averaged over the years 2000-2021. In ECCOv4, the transport is 0.0 Sv averaged over the years 1992-2015. Both datasets reveal a stronger glacial AMOC relative to its pre-industrial state, with a large uncertainty ranging from 2.5 Sv (GloRanV14) to 9.1 Sv (ECCOv4) (Fig. 6b). A stronger glacial AMOC is consistent with enhanced surface-buoyancy forcing over the oceanic deep convection regions and also with previous studies employing isotope modelling.