This study analyzes the response of the Atlantic meridional overturning circulation (AMOC) to different CO2 concentrations and two ice sheet configurations in simulations with the coupled climate model MPI-ESM. With preindustrial (PI) ice sheets, there are two different AMOC states within the studied CO2 range: one state with a strong and deep upper overturning cell at high CO2 concentrations and one state with a weak and shallow upper cell at low CO2 concentrations. Changes in AMOC variability with decreasing CO2 indicate two stability thresholds. The strong state is stable above the first threshold near 217 ppm, and the weak state is stable below the second threshold near 190 ppm. Between the two thresholds, both states are marginally unstable, and the AMOC oscillates between them on millennial time scales. The weak AMOC state is stable when Antarctic Bottom Water becomes dense and salty enough to replace North Atlantic Deep Water (NADW) in the deep North Atlantic and when the density gain over the North Atlantic becomes too weak to sustain continuous NADW formation. With Last Glacial Maximum (LGM) ice sheets, the density gain over the North Atlantic and the northward salt transport are enhanced with respect to the PI ice sheet case. This enables active NADW formation and a strong AMOC for the entire range of studied CO2 concentrations. The AMOC variability indicates that the simulated AMOC is far away from a stability threshold with LGM ice sheets. The nonlinear relationship among AMOC, CO2, and prescribed ice sheets provides an explanation for the large intermodel spread of AMOC states found in previous coupled LGM simulations.
The Atlantic meridional overturning circulation (AMOC) plays an important role in the climate system and is a fundamental component of the ocean circulation. It contributes significantly to the heat transport from the South Atlantic to the North Atlantic and thus influences the climate of the Northern Hemisphere. On the other hand, the climate influences the AMOC because surface fluxes of heat and freshwater affect the formation of North Atlantic Deep Water (NADW) and Antarctic Bottom Water (AABW). These two water masses are tightly coupled to the two overturning cells of the AMOC, and their properties and formation rates can thus alter the strength and geometry of the AMOC. Wind-driven processes, such as upwelling in the Southern Ocean, can also have a strong impact on the AMOC [see, e.g., Kuhlbrodt et al. (2007) for a review]. Because of these feedback processes, the interaction between the AMOC and climate is strongly nonlinear. Reconstructions of the past climate based on ice cores and marine sediment cores have shown that both the climate and the AMOC have varied strongly in the past. While the climate and AMOC were relatively stable during the Holocene and during periods of maximum ice sheet volume, they were very variable in periods of intermediate ice sheet volume (e.g., Grootes et al. 1993; McManus et al. 2004; Böhm et al. 2015). Changes in the stability of the AMOC could play an important role in controlling the variability of the climate, and it is still not fully understood how changes in background climate affect the AMOC stability. We therefore investigate the effect of different background climates on the mean state and variability of the AMOC. To do so, we use a coupled climate model in which we prescribe different combinations of climate forcings ranging from glacial to modern conditions.
Prominent features in the Greenland ice core records were the so-called Dansgaard–Oeschger (DO) events, which describe abrupt warming events over Greenland followed by a gradual cooling period during the more moderate glacial periods of the last glacial cycle (e.g., Dansgaard et al. 1984; Grootes et al. 1993). Heinrich events are another important feature of the glacial period—abrupt cooling events associated with large iceberg discharges from the Laurentide ice sheet during the cold phases of a DO cycle (Heinrich 1988; Hemming 2004). It is a widely accepted hypothesis that these abrupt temperature changes were linked to abrupt state transitions of the AMOC and variations in deep-water formation in the North Atlantic (e.g., Broecker et al. 1985; Kageyama et al. 2010). Records of the ratio of sedimentary 231Pa/230Th over the last glacial cycle indeed indicated that the AMOC was strong during DO warm periods and weakened or even collapsed during some Heinrich events (McManus et al. 2004; Böhm et al. 2015). It is, however, still a topic of debate whether the AMOC is driving these abrupt transitions or whether it is responding to them [see, e.g., Clement and Peterson (2008); Kageyama et al. (2010) for reviews)].
In periods of maximum ice sheet volume, such as the Last Glacial Maximum (LGM; 21 kyr before present), proxy data suggested that the AMOC was shallower and that AABW reached much farther into the North Atlantic than today (Duplessy et al. 1988; Lynch-Stieglitz et al. 2007; Böhm et al. 2015; Lippold et al. 2016). The strength of the AMOC during this time is more difficult to constrain. Reconstructions of the overturning rate based on 231Pa/230Th could be consistent with a weaker glacial AMOC (e.g., McManus et al. 2004) as well as with a glacial AMOC that was at least as strong as today (e.g., Yu et al. 1996; Lippold et al. 2012, 2016). Higher radiocarbon ages in the deep North Atlantic, on the other hand, would support a weaker glacial AMOC (e.g., Skinner et al. 2017). The state of the AMOC in coupled simulations of the LGM differed strongly between different models that participated in the second and third phases of the Paleoclimate Modeling Intercomparison Project (PMIP). Most models also failed to simulate the shallow upper overturning cell (Weber et al. 2007; Muglia and Schmittner 2015; Marzocchi and Jansen 2017) even though the mean surface climate was in reasonable agreement with reconstructions (Braconnot et al. 2007; Braconnot and Kageyama 2015). Different reasons for the intermodel spread have been proposed, such as differences in sea ice formation and brine release in the Southern Ocean (Marzocchi and Jansen 2017) or compensating effects of the LGM greenhouse gas (GHG) concentrations and LGM ice sheets on the AMOC (Klockmann et al. 2016, hereafter K16). Simulations with the ocean model COCO suggested that the state of the glacial AMOC was controlled by a threshold in atmospheric cooling (Oka et al. 2012). It is, however, not clear whether a similar threshold would also exist in coupled climate models (Marotzke 2012). If the threshold existed, and if the LGM climate in the models was near the threshold, it could explain the very different AMOC states in the different models.
The existence of such a threshold would not only have implications for the AMOC during the LGM, but also for the abrupt AMOC transitions documented in the proxy records. A threshold is closely linked to the existence of multiple AMOC states and their stability, and the existence of multiple AMOC states is one possible explanation for abrupt AMOC transitions. In a coupled climate model, such a cooling threshold would be determined by the atmospheric GHG concentrations. Previous studies that looked at the AMOC response to changing GHG concentrations have mostly focused on warm climates with simulations in which the atmospheric CO2 concentration was doubled or quadrupled (e.g., Manabe and Stouffer 1994; Voss and Mikolajewicz 2001a; Stouffer and Manabe 2003; Li et al. 2013; Rugenstein et al. 2016). A few studies have also addressed the AMOC response to reduced GHG concentrations. The AMOC strengthened slightly in response to a reduction of atmospheric CO2 concentrations from 345 to 280 ppm (Voss and Mikolajewicz 2001b), weakened in response to a CO2 reduction from 354 to 200 ppm (Liu et al. 2005), and almost shut down in response to a CO2 reduction from 300 to 150 ppm (Stouffer and Manabe 2003). K16, on the other hand, found an almost-linear AMOC weakening in response to CO2 reductions between 284 and 149 ppm in simulations with a coupled climate model and prescribed LGM ice sheets.
The simulated AMOC response to GHG changes likely also depends on the prescribed ice sheets, including changes in bathymetry and land–sea mask. The amount of cooling necessary to cross the threshold in the simulations with COCO depended on whether preindustrial (PI) or LGM wind stress forcing was applied to the model (Oka et al. 2012), and the wind stress is strongly influenced by the prescribed ice sheets. The previously published studies with changed GHG concentrations used a model setup with either modern ice sheets (Voss and Mikolajewicz 2001b; Stouffer and Manabe 2003) or LGM ice sheets (K16). But none of these studies have assessed a possible nonlinear dependence of the AMOC response to changing CO2 concentrations on the ice sheet configuration. In simulations with CCSM, the AMOC decreased in response to a CO2 doubling with modern ice sheets and increased in response to a CO2 doubling with LGM ice sheets (Zhu et al. 2015). In the simulation with modern ice sheets, the CO2 concentration was doubled from 355 to 710 ppm, while it was doubled from 185 to 370 ppm in the simulations with LGM ice sheets. It is, therefore, difficult to separate whether the different AMOC response was caused by the different ice sheets or by the different absolute GHG concentrations in the two simulations. Based on simulations with transient CO2, Zhang et al. (2017) suggested that the ice sheet volume may affect the CO2 range of a possible AMOC bistability. Brown and Galbraith (2016) found that LGM ice sheets increased the AMOC strength and suppressed millennial-scale AMOC variability that occurred at low CO2 concentrations with preindustrial ice sheets. However, they did not analyze this nonlinear response further because their study focused on the AMOC response to hosing.
We present two sets of simulations with the Max Planck Institute Earth System Model (MPI-ESM) that allow us to analyze the full nonlinearity of the AMOC response to CO2 and ice sheet forcing. The first set of simulations has prescribed PI ice sheets with different CO2 concentrations; the second set has prescribed LGM ice sheets with different CO2 concentrations. The simulations with LGM ice sheets were already introduced by K16, but our present analysis goes substantially beyond that of K16. With the two sets of simulations, we cover a much larger part of the forcing space between the glacial and the modern climate and can now compare how different combinations of ice sheets and CO2 concentrations induce changes in the mean state and variability of the AMOC. We show that the AMOC response to decreasing CO2 concentrations is highly sensitive to the prescribed ice sheets and that a similar threshold, as described by Oka et al. (2012), controls the AMOC state when PI ice sheets are prescribed. We identify the mechanisms that control the threshold behavior and discuss its dependence on the prescribed ice sheets.
We use the MPI-ESM in the coarse-resolution setup without ocean biogeochemistry. The coarse resolution provides a useful compromise between computational cost and model accuracy. The model consists of the atmospheric component ECHAM6.1 (Stevens et al. 2013), the land surface component JSBACH with dynamic vegetation (Reick et al. 2013), and the ocean component MPI-OM (Marsland et al. 2003; Jungclaus et al. 2006). A sea ice model is incorporated into MPI-OM [see Marsland et al. (2003); Notz et al. (2013) for details]. Except for the dynamical vegetation and the coarse resolution, the model configuration applied here corresponds to MPI-ESM-P, which can be found in the CMIP5/PMIP3 database. The model configuration is the same as described in greater detail in K16 and Klockmann (2017). The spectral resolution of ECHAM in the coarse-resolution setup is T31, which corresponds to a horizontal resolution of approximately 3.75° × 3.75° in grid space. In the vertical, there are 31 σ-hybrid layers. MPI-OM uses a curvilinear grid configuration with a nominal resolution of 3° × 3° in the coarse-resolution setup. The resolution increases toward the grid poles, which are located over southeast Greenland and central Antarctica. The grid spacing increases from approximately 350 km in the tropics to 31 km around Greenland and 86 km around Antarctica. The water column is resolved by 40 levels in the vertical. The first model level has a thickness of 15 m. Below that, the level thickness increases from 10 m in the near-surface levels to approximately 550 m in the deepest level.
We perform a set of five simulations with prescribed PI ice sheets and different GHG concentrations. The radiative forcing between the respective simulations decreases approximately linearly. The atmospheric CO2 concentrations range from 149 to 353 ppm. The simulations are named Pxxx, where xxx denotes the atmospheric CO2 concentration in ppm (see Table 1). We compare the Pxxx simulations with their equivalent simulations with prescribed PMIP3 LGM ice sheets (Abe-Ouchi et al. 2015), which were presented in K16. The simulations with LGM ice sheets are named LGMxxx (see Table 1). The simulations P284 and LGM185 are identical to the simulations piORB and LGMref in K16. The orbital parameters are set to LGM values in all simulations to exclude the orbital effect in the comparison with the simulations with LGM ice sheets (see K16 for details on the orbital effect). A detailed description and evaluation of the simulated preindustrial and LGM climate was also given by K16. When we refer to the PI or LGM ice sheets in this paper, we always implicitly refer to the corresponding orography, land–sea mask, and ocean bathymetry. When we refer to CO2 or CO2 concentrations, we always refer to the atmospheric CO2 concentration.
We perform three additional simulations with PI ice sheets and a smaller difference in radiative forcing between 230 and 185 ppm in order to investigate a nonlinear AMOC transition from a strong to a weak state that occurs in this CO2 range (see P217, P206, and P195 in Table 1). In addition, we perform a sensitivity experiment in which we reduce the brine release in the Southern Ocean (see P185_lessbrine in Table 1) in order to investigate the role of brine release in the AMOC transition.
The ocean is initialized from fields that were adjusted to be reasonably close to the expected equilibrium climate to shorten the spinup period. The initial states are obtained either by globally adjusting the ocean temperature by a best-guess estimate or by linearly combining states from previous simulations (see “Initial state” in Table 1). All simulations are integrated until the AMOC reached a quasi-steady state (i.e., until there was no discernible AMOC drift over the last 500 simulation years). To analyze the steady-state response, we use averages of the last 300 years of the respective simulations. To analyze the temporal variability, we use time series of yearly means, which are smoothed with a 31-yr running mean.
3. Two AMOC states
The AMOC weakens and shoals in response to decreasing CO2 concentrations, independent of the prescribed ice sheets (Fig. 1). But the magnitude and the degree of nonlinearity of the response depend strongly on the ice sheets. With LGM ice sheets, the weakening is gradual and quasi-linear for CO2 concentrations below 284 ppm (Fig. 1b). A gradual shoaling of the interface between the two overturning cells (defined as the depth of zero transport) sets in for CO2 concentrations below 230 ppm. With PI ice sheets, on the other hand, the response is strongly nonlinear; there are two distinct AMOC states. At high CO2 concentrations (in P353, P284, and P230), the AMOC is in a state with a strong and deep upper overturning cell (Fig. 1a). The overturning strength and geometry in P353 and P284 are nearly indistinguishable; in both simulations, the upper overturning cell has a strength of approximately 16 Sv (1 Sv ≡ 106 m3 s−1) and extends down to 2900 m. In P230, the upper overturning cell has a strength of 13 Sv. At low CO2 concentrations (in P185 and P149), the AMOC is in a state with a weak and shallow upper overturning cell. When the CO2 concentration is reduced from 230 to 185 ppm, the AMOC strength decreases abruptly from 13 Sv in P230 to 7 Sv in P185, and the extent of the upper overturning cell shoals from approximately 2800 to 2550 m. In P149, the strength of the upper overturning cell is 6 Sv, and the upper overturning cell reaches down to approximately 2500 m.
In summary, the AMOC is more sensitive to a CO2 reduction with prescribed PI ice sheets than with prescribed LGM ice sheets. It should be noted that the AMOC in the LGM185 simulation is too deep and likely too strong in comparison with the AMOC reconstructions for the LGM (e.g., Lynch-Stieglitz et al. 2007) and would probably not be reconcilable with the increased radiocarbon ages in the glacial Atlantic (Skinner et al. 2017). This is due to the ice sheet effect, which strengthens and deepens the AMOC in MPI-ESM in combination with insufficient brine release in the Southern Ocean under LGM boundary conditions (K16). The overall stronger AMOC may explain why the AMOC is less sensitive to the CO2 reduction in the simulations with LGM ice sheets. A stronger AMOC may be less likely to become unstable than a weaker AMOC (e.g., Tziperman 2000). We will return to this issue in section 5. In the following sections, we will investigate the transition from the strong to the weak AMOC state in detail and focus on the processes that control the AMOC state at low CO2 concentrations with PI and LGM ice sheets.
4. Two stability thresholds
The abrupt transition from the strong to the weak AMOC state in the simulations with PI ice sheets indicates the crossing of a stability threshold. If a system is in a stable state, a perturbation will always be dampened by negative feedbacks, and the perturbed system will quickly return to its unperturbed state. Close to a bifurcation point beyond which the system becomes unstable, the negative feedbacks weaken, and small perturbations can result in large changes in the system (e.g., Held and Kleinen 2004). A good indicator for the presence of a stability threshold is, therefore, the increase of noise or temporal variability as the threshold is approached (e.g., Tziperman 2000; Knutti and Stocker 2002; Kleinen et al. 2003; Armstrong et al. 2017). As a simple measure for the temporal variability, we compare the temporal standard deviation of the AMOC strength in the simulations with PI ice sheets. We include the three additional simulations (P217, P206, and P195) into the analysis to have more data points in the CO2 range in which the transition occurs. We estimate the standard deviation from time series smoothed with a 31-yr running mean because we are mainly interested in the low-frequency variability.
In the range of 353–217 ppm, the temporal standard deviation of the overturning strength increases almost exponentially with decreasing CO2 concentrations (black line in Fig. 2, read from right to left). At 353 ppm, the standard deviation is 0.5 Sv. At 284 ppm, it increases by 0.1 Sv to 0.6 Sv. At 230 ppm, it increases further by 0.7 Sv to 1.3 Sv, and at 217 ppm, it increases by 1.4 Sv to 2.7 Sv. The standard deviation peaks at 206 ppm with about 3.7 Sv. Below 206 ppm, the standard deviation decreases with decreasing CO2 concentrations. At 195 ppm, the standard deviation is only slightly smaller, with 3.4 Sv. Below 195 ppm, the standard deviation then decreases sharply to about 1.1 Sv at 185 ppm and 0.8 Sv at 149 ppm.
By estimating the time between two AMOC maxima in the respective AMOC time series (Fig. 3), we can give an estimate of the time scale of the low-frequency variability. The larger the amplitude of the AMOC variability, the longer the associated time scale becomes. At 353 and 284 ppm, the dominant time scale is 200–250 years (Fig. 3a). At 230 ppm, the time scale increases to 500–700 years. There is even one large AMOC excursion after approximately 3900 simulation years, which takes about 1000 years to complete. In the range of 217–195 ppm, the AMOC starts to oscillate between a strong and a weak AMOC state (Fig. 3b). The high standard deviation corresponds roughly to the amplitude of the oscillations. At 217 ppm, the AMOC oscillates around a mean state of approximately 12 Sv; one oscillation takes about 750 years to complete. At 206 ppm, the AMOC switches quite abruptly from a weak state to a strong state and back again. The transition period is about 200–300 years long, and the AMOC spends 500–700 years in the strong state and 1000–1200 years in the weak state. The strong state has a mean strength of approximately 16 Sv, and the weak state has a mean state of about 7.5 Sv with strong variability on centennial time scales. At 195 ppm, similar state transitions occur, but the AMOC spends less time in the strong state and more time in the weak state. At 185 and 149 ppm, the AMOC switches into the weak state after a spinup period, during which the AMOC is still strong, and then remains in the weak state for the rest of the respective simulations (Fig. 3a). When the AMOC is in the weak state, the time scale of variability is 300–500 years.
We conclude from the changes in the temporal variability that there are two thresholds of CO2-induced cooling if PI ice sheets are prescribed. The first threshold is located near 217 ppm, and the strong AMOC state is stable above this threshold and unstable below. The second threshold is located between 185 and 195 ppm, and the weak AMOC state is unstable above this threshold and stable below. In between the two thresholds, both the strong and the weak states are marginally unstable, and the AMOC oscillations in P206 and P195 are a result of oscillations between two unstable states (e.g., Colin de Verdière 2007). With LGM ice sheets, on the other hand, the AMOC remains in the regime of the stable strong state; the temporal variability is almost constant at about 0.8 Sv in the CO2 range of 353–230 ppm and then decreases slightly to 0.6 and 0.4 Sv at 185 and 149 ppm, respectively (Fig. 2, gray line). In the remainder of the study, we focus on the conditions that generate the stable weak AMOC state in P185 and P149 and discuss why the weak AMOC state does not occur in the simulations with LGM ice sheets.
5. Conditions for the stability of the weak AMOC state
A first indication of why the AMOC does not switch into a weak state with LGM ice sheets is the AMOC strength itself. A weak AMOC with reduced heat transport and a weaker thermal forcing was found to be inherently less stable than a strong and thermally dominated AMOC (e.g., Tziperman 2000; Arzel et al. 2012). In the simulations with LGM ice sheets, the AMOC is stronger than in all simulations with PI ice sheets. Even the weakest AMOC with LGM ice sheets (18 Sv in LGM149) is still stronger than the strongest AMOC with PI ice sheets (16 Sv in P353 and P284). The bifurcation could thus be linked to a critical AMOC strength. According to the simulations with PI ice sheets, the critical AMOC strength would be slightly weaker than 13 Sv. With a strength of approximately 18 Sv, the AMOC in LGM149 is still very far away from this critical point.
The divergence of the freshwater transport associated with the AMOC can also be used as an indicator of the AMOC stability [e.g., Liu et al. (2017) and references therein]. The freshwater-transport divergence dMOV is given by the difference of the AMOC-induced freshwater transport at 80°N and 34°S. A positive dMOV indicates a monostable AMOC, while a negative dMOV indicates a bistable AMOC. For our simulations, the stability indicator suggests that the AMOC is in a strongly monostable regime with LGM ice sheets and in a neutral to weakly bistable regime with PI ice sheets (Fig. 4). This gives additional support to the hypothesis that the strong AMOC in the simulations with LGM ice sheets is far away from the bifurcation point beyond which the weak AMOC state occurs. However, based on the AMOC strength and the stability indicator alone, we cannot say anything about the factors that determine the AMOC strength and geometry in our simulations. In the following, we therefore analyze the properties and formation of NADW and AABW.
a. The north–south salinity difference
The overturning strength and geometry is to a large extent set by the density difference and properties of NADW and AABW (e.g., Kuhlbrodt et al. 2007). We therefore analyze the changes in NADW and AABW with decreasing CO2 concentrations. The analysis is analogous to K16. We determine the temperature and salinity of the two water masses at the in situ density maximum at approximately 1500 m in the North Atlantic and the Weddell Sea, respectively (see small map in Fig. 5b). Though 1500 m may seem rather shallow for the core properties of present-day NADW, we choose this depth because the core of NADW shifts upward in the simulations with a weak AMOC state, and a depth of, for example, 2000 m would be no longer representative of NADW properties. In the simulations with a strong AMOC state, the water column in the North Atlantic is quite homogeneous between 1200 and 2000 m, and thus 1500 m is also representative of the NADW properties of the strong AMOC state. We show the same TS diagram for the simulations with LGM ice sheets for comparison. We have also included a comparison of the water masses in P284 and LGM185 with observations (PHC3 climatology; Steele et al. 2001) and reconstructions (Adkins et al. 2002).
In the simulations with PI ice sheets, NADW becomes colder and fresher with decreasing CO2 (circles in Fig. 5a). In the range of 353–230 ppm, the temperature effect on density dominates, and the potential density increases. In the range of 230–149 ppm, the cooling weakens, and the salinity effect on density dominates so that the potential density decreases. AABW becomes colder and saltier with decreasing CO2 (triangles in Fig. 5a), and the potential density increases. The cooling takes place until the temperature reaches the freezing point at 185 ppm. Afterward, the density increase is due to the salinity increase alone. At 185 ppm, AABW and NADW have approximately the same salinity; at 149 ppm, the north–south salinity difference reverses, and AABW becomes saltier than NADW. In the simulations with LGM ice sheets, the changes in the water mass properties are qualitatively very similar (Fig. 5b).
In comparison with observations, the density difference between NADW and AABW is represented quite well in P284, even though the simulated NADW is somewhat warmer and saltier than in the PHC3 climatology, and AABW is slightly colder and fresher (Fig. 5c, orange symbols). Compared to reconstructions, the increase in the density difference between AABW and NADW at the LGM is strongly underestimated in LGM185 (Fig. 5c, blue symbols). This underestimation is mainly due to the too-fresh AABW, which is a consequence of insufficient Southern Ocean brine release at glacial CO2 concentrations (K16). In K16, we have further shown that a shoaling of the glacial upper overturning cell with respect to the preindustrial state could only be simulated at lower CO2 concentrations, when the strengthened brine release induced a reversal of the salinity difference between NADW and AABW. In the following, we show that the north–south salinity difference is also a key factor in the transition from the strong to the weak AMOC state. We focus on salinity rather than density here because the salinity changes control the changes in the density difference between NADW and AABW.
We hypothesize that AABW needs to be at least as salty as NADW to have a stable weak and shallow AMOC state. This hypothesis is supported by the temporal evolution of the salinity during the spinup of P185 (Fig. 6a). The salinity of AABW (defined as Weddell Sea salinity at 1500 m) increases by 0.1 g kg−1 during the first 1400 years of the simulation and then remains constant for the remaining 2600 years. The salinity increase is caused by brine release and shelf convection in the Southern Ocean in combination with a weak freshwater input from precipitation and runoff (not shown). At the same time, the salinity of NADW decreases by about 0.1 g kg−1 during the first 1200 years of the simulation. The AMOC decreases slowly by about 3 Sv during the first 1300 years and then rapidly by 8 Sv in the following 200 years. As long as the AMOC is strong, NADW is formed both in the North Atlantic and in the southern Nordic Seas. The rapid AMOC decrease sets in when the Nordic Seas stop contributing to the NADW formation (see also next section). The rapid AMOC decrease is accompanied by a freshening of 0.2 g kg−1 of NADW. Once the AMOC has reached a stable weak state, NADW and AABW have approximately the same salinity.
To test whether the salinity in the Southern Ocean controls the stability of the AMOC states, we perform a sensitivity experiment with reduced brine release in the Southern Ocean (P185_lessbrine; see Table 1). To reduce the brine release, we increased the salinity of sea ice in the Southern Ocean from 5 to 20 g kg−1. This way, the brine release reduces by approximately 50%, and the oceanic salt budget is conserved. Otherwise, the setup is identical to P185. The P185_lessbrine simulation is branched off P185 during the spinup phase after 900 years when the AMOC is still strong. The reduced brine release causes a freshening of 0.06 g kg−1 in the Weddell Sea. The associated density decrease of AABW allows for the NADW to penetrate deeper levels, and the percentage of NADW in the deep North Atlantic increases with respect to the P185 simulation, in which the deep North Atlantic is dominated by AABW (not shown). The resulting changes in stratification then affect the AMOC state (e.g., K16; Jansen and Nadeau 2016). The AMOC in P185_lessbrine remains in the strong state for the entire length of the simulation (Fig. 6b, black line). The P185_lessbrine simulation is integrated for 4000 years to ensure that no transition to the weak AMOC state occurs at a later point. The upper overturning cell deepens by approximately 400 m (not shown).
The strong and deep AMOC in P185_lessbrine supports the hypothesis that the stability of the strong and weak AMOC states is controlled by the north–south salinity difference and that AABW needs to be at least as salty as NADW in order to have a stable weak and shallow AMOC state. It is, however, a necessary but not a sufficient condition. In the simulations with the LGM ice sheets, the north–south salinity difference also reverses between 185 and 149 ppm, but there is no transition to a weak AMOC state. In the next section, we show that this is due to processes in the North Atlantic.
b. NADW formation
The density budget over the convective regions in the North Atlantic can give insight into the processes that drive the formation of NADW and help to characterize the two AMOC states further. We define the net surface density flux as the sum of the density fluxes due to heat fluxes, sea ice formation and melt, and atmospheric freshwater input. The flux components are then integrated separately over the Nordic Seas and over the North Atlantic, including the Labrador Sea, in order to estimate the relative importance of the NADW formation sites north and south of the Greenland–Scotland Ridge. To keep the area the same for all simulations, we define the integration area as all grid points in which the mixed layer depth in any of the Pxxx simulations exceeds 1000 m.
The Nordic Seas contribute to NADW formation via the Nordic Sea overflows. A stronger overflow is known to increase the AMOC strength and NADW penetration depth (e.g., Danabasoglu et al. 2010; Zhang et al. 2011). Recent observation-based estimates indicate a volume transport of dense water (σΘ > 27.8 kg m−3) through the Denmark Strait of approximately 3.2 ± 0.5 Sv (Jochumsen et al. 2017). In P284, the water that is transported through the Denmark Strait is slightly lighter (σΘ > 27.7 kg m−3), but the volume transport is comparable with about 2.8 Sv. Given that also the NADW in P284 is about 0.1 kg m−3 lighter than in the PHC3 climatology (Fig. 5c), we are confident that the Nordic Seas contribute significantly to NADW formation in P284 and that changes in the surface density budget over the Nordic Seas will affect the total NADW formation.
Over the Nordic Seas, the net density flux is positive (i.e., a density gain) when the AMOC is in the strong state (Fig. 7a, right panel). At 353 ppm, the dominant contribution to the net density flux is a density gain due to heat loss. Sea ice melt and atmospheric freshwater input induce a density loss, but they are outweighed by the heat flux contribution. With decreasing CO2, the heat flux contribution weakens. The water that enters the Nordic Seas becomes colder (not shown), and its cooling potential decreases while the water temperature approaches the freezing point. As the heat flux contribution decreases, the net density gain also decreases. At 284 ppm, the heat flux contribution is still strong enough to outweigh the density loss due to sea ice melt and freshwater input. At 230 ppm, the three contributions add up to a very small net density gain. When the AMOC is in the weak state, the heat flux contribution is zero. The Nordic Seas become completely ice covered, and the sea ice formation induces a small density gain, which is compensated for by the atmospheric freshwater input. The resulting net density flux is almost zero at 185 and 149 ppm, and the Nordic Seas do not contribute to NADW formation when the AMOC is in the weak state.
Over the North Atlantic and Labrador Sea, the net density gain is also mostly due to heat loss when the AMOC is in the strong state (Fig. 7a, left panel). At 353 and 284 ppm, the density gain due to heat loss is counteracted only by a comparatively small density loss due to freshwater input. The net sea ice contribution is close to zero. At 230 ppm, the net density gain decreases because of decreased heat loss and increased sea ice melt. When the AMOC is in the weak state, the net density gain is negative (i.e., a density loss) in the long-term mean. The heat flux contribution is still positive, which is due to sporadic deep convection events. At 185 and 149 ppm, deep convection events in the eastern North Atlantic occur only every 10–15 and 20–50 years, respectively. During these events, there is a net density gain, but in the long-term mean, the heat flux contribution is outweighed by the sum of contributions from sea ice melt and atmospheric freshwater input (not shown).
In the simulations with LGM ice sheets, the density gain over the North Atlantic and Labrador Sea is about 3 times stronger than in the simulations with PI ice sheets (cf. Figs. 7a and 7b). The heat loss is strong enough to outweigh the freshwater input and sea ice melt also at low CO2 concentrations, and there is a net density gain in all simulations with LGM ice sheets. In LGM149, the net density gain over the North Atlantic is still of a similar magnitude as that in P284. Over the Nordic Seas, the net density flux becomes zero at 230 ppm, just like in the simulations with PI ice sheets. But it appears that the AMOC is less sensitive to deep convection in the Nordic Seas when LGM ice sheets are prescribed. With LGM ice sheets, the density gain over the North Atlantic is strong enough to sustain a strong AMOC also when the Nordic Seas do not contribute to NADW formation.
The difference in density fluxes between the simulations with PI ice sheets and with LGM ice sheets shows that the density gain in the NADW formation areas plays a crucial role in regulating the AMOC strength. However, the coupling between the AMOC and the surface density budget involves a number of self-amplifying feedbacks that make it difficult to directly attribute causality. A strong AMOC transports warm and salty water into the subpolar North Atlantic, thus weakening the stratification and keeping the sea ice cover reduced. As a consequence, there is a strong heat loss and active deep convection that sustains the strong AMOC. A weak AMOC, on the other hand, is associated with a weaker northward transport of salt and heat. This increases the stratification in the subpolar North Atlantic and thus favors the expansion of sea ice. The weak northward heat transport and the insulating effect of the sea ice keep the density gain due to heat loss small and the AMOC in a weak state.
There are, however, other factors that impact the surface density gain over the North Atlantic and that are not directly related to the AMOC strength. The LGM ice sheets modify the wind stress over the North Atlantic and thus affect also the wind-driven horizontal ocean circulation (i.e., the subtropical–subpolar gyre system). At latitudes north of about 45°N, the subpolar gyre plays a more important role in the northward transport of heat and salt than the AMOC (e.g., Jungclaus et al. 2013); thus, changes in the subpolar gyre will have a strong impact on the density gain in the North Atlantic. The mean salinity and velocity of the upper 500 m illustrate this very well (Figs. 8c,d; shown are P185 and LGM185). In the simulations with PI ice sheets, the subpolar gyre (SPG) extends far eastward, and the North Atlantic and Nordic Seas are dominated by relatively fresh subpolar water. Only very little saline water of subtropical origin reaches the deep convection areas in the Iceland Basin and Irminger Sea. The inflow is visible as a narrow band of saline water off the Bay of Biscay and the Irish coast. In the simulations with LGM ice sheets, on the other hand, the subtropical gyre (STG) shifts northward and the SPG contracts. The inflow of salty subtropical water extends over the entire Atlantic basin east of the Mid-Atlantic Ridge. This inflow supplies salty water to the deep-convection sites in the Iceland Basin and Irminger Sea and thus helps to maintain continuous deep convection and a strong AMOC also at low CO2 concentrations. A similar relationship among the northward salt transport, North Atlantic wind stress, and the glacial overturning strength was described by Muglia and Schmittner (2015) and Sherriff-Tadano et al. (2018). The surface wind stress may also affect the surface density budget through modification of the sea ice distribution [see also Zhu et al. (2014); Figs. 8a,b]. In addition, the advection of very cold air from the Laurentide ice sheet to the open ocean in the North Atlantic can further enhance the heat loss and deep convection.
Previous studies have found multiple AMOC states and AMOC transitions for a variety of different climate states and forcing combinations. The AMOC in simulations with the coupled model CM2Mc behaved very similarly to the AMOC in MPI-ESM (Brown and Galbraith 2016). The AMOC oscillated between a weak and a strong state in a simulation with 180 ppm and PI ice sheets. At higher CO2 concentrations, the AMOC was strong, and at lower CO2 concentrations, the AMOC was weak. In simulations with LGM ice sheets, the AMOC strength increased, and no oscillations occurred. Small CO2 changes from 210 to 225 ppm in transient simulations with the coupled model ECHAM5/MPI-OM induced abrupt AMOC transitions from a weak to a strong state, and vice versa (Zhang et al. 2017). The ice sheets were prescribed at intermediate volume. Zhang et al. (2017) suggested that lower CO2 concentrations would be needed to force abrupt AMOC transitions with increasing ice sheet volume. The CO2 range of 210–225 ppm is in good agreement with the CO2 range in which both AMOC states are unstable in our simulations with PI ice sheets. But while the AMOC transitions in ECHAM5/MPI-OM were forced by transient CO2 changes, the AMOC transitions in our simulations occur unforced because there are two marginally unstable AMOC states. While there seems to be a consensus that intermediate to low CO2 concentrations are needed for abrupt AMOC transitions, AMOC oscillations and abrupt state transitions have been found in simulations with both LGM ice sheets (Wang and Mysak 2006; Arzel et al. 2012; Peltier and Vettoretti 2014) and PI ice sheets (e.g., Friedrich et al. 2010; Brown and Galbraith 2016). In addition to changing the surface wind stress, different ice sheets can also open or close oceanic gateways such as Bering Strait or the Canadian Archipelago and change the throughflow depth of passages such as the Denmark Strait. These changes can affect the freshwater and density budget of the North Atlantic and thus also the AMOC stability (e.g., Hu et al. 2012). In the different models, the total effect of the ice sheets on the AMOC stability then depends on the exact ice sheet reconstruction, the model-specific implementation of the ice sheets, and the resulting circulation pattern in the North Atlantic.
The two stability thresholds in our simulations with PI ice sheets are consistent with the thermal threshold proposed by Oka et al. (2012). In their simulations with PI wind stress forcing, the AMOC switched into a weak state when the glacial share of the heat flux forcing was approximately 60%. In their simulations with LGM wind stress forcing, the threshold lay at a glacial share of 140% (i.e., a much stronger cooling was needed to force the AMOC into a weak state when the LGM wind stress was used to force the model). In our simulations with LGM ice sheets, we do not find any stability threshold within the studied CO2 range (Figs. 1b, 2). A rough extrapolation of the density fluxes in Fig. 7b suggests that a further CO2 reduction might reduce the heat loss over the North Atlantic enough to reach the first stability threshold at approximately 100 ppm. On the other hand, the very low temporal variability in all simulations with LGM ice sheets indicates that the AMOC is still far from becoming unstable (Fig. 2). Oka et al. (2012) proposed that the existence of the thermal threshold could be an explanation for the different AMOC states in the LGM simulations of the PMIP2 ensemble. In the following, we argue that the nonlinear AMOC response to decreasing GHG concentrations and the opposing effects of the ice sheets and CO2 concentrations may indeed explain the large AMOC spread in the PMIP simulations.
K16 showed that the partial compensation of ice sheet effect and GHG effect was the reason for the lack of shoaling of the glacial upper overturning cell in MPI-ESM. K16 have estimated the GHG effect as the difference between LGM185 and LGM284 and the ice sheet effect as the difference between LGM284 and P284 (see label “i” in Fig. 9). We can now add a second estimate of the GHG effect by comparing P185 and P284 and a second estimate of the ice sheet effect by comparing LGM185 and P185 (see label “ii” in Fig. 9). Comparing the two sets of estimates gives a measure of how much the respective effects depend on the background climate. The GHG effect on the AMOC strength is stronger with PI ice sheets (−9 Sv) than with LGM ice sheets (−4 Sv). In the same way, the ice sheet effect on the AMOC strength is stronger at low CO2 (+13 Sv) than at high CO2 (+8 Sv). On the other hand, the change in the depth of the interface between the two overturning cells is quite similar in both sets of estimates (±300 and ±350 m, respectively).
Regardless of the nonlinear response of the AMOC strength, the glacial AMOC response in MPI-ESM is always the combination of two larger opposing effects. A similar compensation might happen in other models, too. But even if all models agreed on the sign of the individual responses, small differences in the magnitudes of either effects can change the sign of the total response. Just as the magnitude of the respective effects depends on the background climate in MPI-ESM, the magnitude of the respective effects in different models very likely also depends on model specifics. Brine release in the Southern Ocean is crucial for the AMOC response to low GHG concentrations (e.g., Shin et al. 2003; Ferrari et al. 2014; Jansen and Nadeau 2016; K16), and it seems to be insufficient in many PMIP3 models (Marzocchi and Jansen 2017). Because the brine release is a function of sea ice formation, the GHG effect in different models depends crucially on the implementation of the sea ice and the representation of the wind field over the Southern Ocean. Insufficient spinup duration may also play a role (Marzocchi and Jansen 2017). The slow equilibration of the deep ocean is also illustrated by the spinup of P185. It takes about 1300 years of spinup until the transition to the weak state occurs. For the ice sheet effect in different models, both the specific reconstruction of the ice sheets as well as the implementation of the ice sheets may play a large role. Even the same ice sheet reconstruction resulted in a very different representation of the ice sheets in individual PMIP3 models because of differences in the implementation (Chavaillaz et al. 2013). This likely affects the wind stress over the North Atlantic and the northward salt transport and thus the formation of NADW and the magnitude of the ice sheet effect. In the upcoming fourth PMIP phase, the modeling groups will have a choice between different ice sheet reconstructions (Kageyama et al. 2017). Comparing the northward salt transport for the different ice sheet reconstructions may increase the understanding of the simulated glacial AMOC and may further help to explain why AMOC oscillations were found in simulations with both LGM and PI ice sheets.
Based on the two sets of simulations with PI and LGM ice sheets, we conclude the following:
The effect of low GHG concentrations on the AMOC depends on the prescribed ice sheets. With LGM ice sheets, the AMOC strength decreases gradually when the CO2 concentration is lower than 284 ppm. With PI ice sheets, there are two distinct states: a state with a strong and deep upper overturning cell at high CO2 concentrations and a state with a weak and shallow upper overturning cell at low CO2 concentrations.
We find two stability thresholds that determine the stability of the two AMOC states with PI ice sheets. The location of the thresholds will likely depend on model specifics. In our simulations, the first threshold is located near 217 ppm. Above this threshold, the strong AMOC state is stable within the studied CO2 range. The second threshold is located between 195 and 185 ppm. Below this threshold, the weak AMOC state is stable within the studied CO2 range. Between the thresholds, both states are marginally unstable, and the AMOC oscillates between them on millennial time scales.
A stable weak AMOC state occurs in response to changes in both the Southern Ocean and the North Atlantic. First, AABW needs to be salty and dense enough to replace NADW in the North Atlantic. Second, the density gain over the North Atlantic is weak, and the Nordic Seas do not contribute to the formation of NADW.
The LGM ice sheets increase the upper-ocean salt transport into the NADW formation sites and enhance the density gain over the North Atlantic. With LGM ice sheets, the AMOC is thus far away from the first stability threshold, and the AMOC remains in the strong state also at low CO2 concentrations.
The glacial AMOC response in MPI-ESM is the sum of two large opposing effects. The sign of the GHG and ice sheet effects is robust, but their respective magnitudes depend on the background climate and likely on model specifics. Small changes in either of the two effects can change the sign of the total response. If a similar compensation of effects takes place in other models, this provides an explanation for the large spread of glacial AMOC states in the PMIP2/PMIP3 ensemble.
The authors thank Chao Li and three anonymous reviewers for their helpful comments. This study was part of the Ph.D. thesis of M.K. (Klockmann 2017) and was funded by the International Max Planck Research School on Earth System Modelling. This study is a contribution to the BMBF project PalMod. All simulations were performed on the supercomputers of the Deutsches Klimarechenzentrum (DKRZ). The code of MPI-ESM is available at http://www.mpimet.mpg.de/en/science/models/mpi-esm.html. The scripts for analysis and processing can be obtained at http://hdl.handle.net/11858/00-001M-0000-002E-9BE0-9. Selected model output can be obtained from the CERA database at https://cera-www.dkrz.de/WDCC/ui/cerasearch/entry?acronym=DKRZ_LTA_033_ds00002. Additional model output can be requested directly from M.K.
Current affiliation: Institute of Coastal Research, Helmholtz-Zentrum Geesthacht, Centre for Materials and Coastal Research, Geesthacht, Germany.