Response of a comprehensive climate model to a broad range of external forcings: relevance for deep ocean ventilation and the development of late Cenozoic ice ages

Abstract

Over the past few million years, the Earth descended from the relatively warm and stable climate of the Pliocene into the increasingly dramatic ice age cycles of the Pleistocene. The influences of orbital forcing and atmospheric CO₂ on land-based ice sheets have long been considered as the key drivers of the ice ages, but less attention has been paid to their direct influences on the circulation of the deep ocean. Here we provide a broad view on the influences of CO₂, orbital forcing and ice sheet size according to a comprehensive Earth system model, by integrating the model to equilibrium under 40 different combinations of the three external forcings. We find that the volume contribution of Antarctic (AABW) vs. North Atlantic (NADW) waters to the deep ocean varies widely among the simulations, and can be predicted from the difference between the surface densities at AABW and NADW deep water formation sites. Minima of both the AABW-NADW density difference and the AABW volume occur near interglacial CO₂ (270–400 ppm). At low CO₂, abundant formation and northward export of sea ice in the Southern Ocean contributes to very salty and dense Antarctic waters that dominate the global deep ocean. Furthermore, when the Earth is cold, low obliquity (i.e. a reduced tilt of Earth’s rotational axis) enhances the Antarctic water volume by expanding sea ice further. At high CO₂, AABW dominance is favoured due to relatively warm subpolar North Atlantic waters, with more dependence on precession. Meanwhile, a large Laurentide ice sheet steers atmospheric circulation as to strengthen the Atlantic Meridional Overturning Circulation, but cools the Southern Ocean remotely, enhancing Antarctic sea ice export and leading to very salty and expanded AABW. Together, these results suggest that a ‘sweet spot’ of low CO₂, low obliquity and relatively small ice sheets would have poised the AMOC for interruption, promoting Dansgaard–Oeschger-type abrupt change. The deep ocean temperature and salinity simulated under the most representative ‘glacial’ state agree very well with reconstructions from the Last Glacial Maximum (LGM), which lends confidence in the ability of the model to estimate large-scale changes in water-mass geometry. The model also simulates a circulation-driven increase of preformed radiocarbon reservoir age, which could explain most of the reconstructed LGM-preindustrial ocean radiocarbon change. However, the radiocarbon content of the simulated glacial ocean is still higher than reconstructed for the LGM, and the model does not reproduce reconstructed LGM deep ocean oxygen depletions. These ventilation-related disagreements probably reflect unresolved physical aspects of ventilation and ecosystem processes, but also raise the possibility that the LGM ocean circulation was not in equilibrium. Finally, the simulations display an increased sensitivity of both surface air temperature and AABW volume to orbital forcing under low CO₂. We suggest that this enhanced orbital sensitivity contributed to the development of the ice age cycles by amplifying the responses of climate and the carbon cycle to orbital forcing, following a gradual downward trend of CO₂.

1 Introduction

Since the first quantitative reconstructions of ocean temperatures were made (Emiliani 1955), investigators have pondered how changes in climate state may have altered deep ocean circulation during the ice age cycles of the past few million years (Adkins 2013). The greatest focus has traditionally been placed on understanding the most recent extreme state, the Last Glacial Maximum (LGM), typically defined as the period from 19,000 to 23,000 years ago (Mix et al. 2001). Numerous mechanisms have been championed as key drivers behind glacial-interglacial changes in ocean state, including changes in southern westerly winds (Toggweiler et al. 2006; Anderson et al. 2009), sea ice dynamics (Stephens and Keeling 2000; Schmittner 2003; Shin et al. 2003; Bouttes et al. 2010; Ferrari et al. 2014), ice sheet topography (Pausata et al. 2011; Gong et al. 2015), tidal mixing (Wunsch 2003; Schmittner et al. 2015), and the nonlinear dependence of density on temperature (Winton 1997; Sigman et al. 2004; De Boer et al. 2007). In turn, such physical oceanographic changes may have contributed to changes in atmospheric CO₂ concentrations—hereafter referred to simply as “CO₂”—measured in ice core records. If so, the physical oceanographic changes would have been integral to the dynamics of ice age cycles (Sigman and Boyle 2000; Ganopolski et al. 2016). The close correlation of atmospheric CO₂ with ice core records of Antarctic temperature over the past 800,000 years, combined with proxy evidence of accompanying changes in the Southern Ocean, has inspired many authors to link changes in Southern Ocean circulation to oceanic CO₂ sequestration (Paillard and Parrenin 2004; Fischer et al. 2010; Sigman et al. 2010; Skinner et al. 2010; Jaccard et al. 2013; Ferrari et al. 2014; Watson et al. 2015).

A number of popular hypotheses call on saline brine, rejected from a massive apron of sea ice around Antarctica, to have produced an abundant source of dense water during the last ice age that filled the global abyss (Adkins et al. 2002; Bouttes et al. 2010; Adkins 2013). The idea that Antarctic Bottom Water (AABW) was much saltier during ice ages arose from reconstructions made by extracting fossil pore-water from deep sea sediment cores (Schrag and de Paolo 1993; Adkins and Schrag 2003). After accounting for changes in global ice volume, these reconstructions suggest that the salinity of the deep ocean was significantly elevated relative to the upper ocean during the LGM (Adkins et al. 2002; Insua et al. 2014). An Antarctic origin of these salty bottom waters would be consistent with reconstructions of the chemical composition of deep waters, primarily reflected by isotopic ratios in the shells of fossil benthic foraminifera. Both the stable carbon isotope ratio, δ¹³C (Boyle and Keigwin 1982; Duplessy et al. 1984; Sarnthein et al. 1994; Curry and Oppo 2005) and neodymium isotope ratio, εNd (Piotrowski et al. 2005; Howe et al. 2016) indicate that dramatic changes occurred in the Atlantic basin, with a volumetric expansion of the Antarctic chemical fingerprint under glacial conditions. Although the δ¹³C evidence is not unambiguous (Gebbie 2014), when viewed together with the spatial distribution of εNd changes, there is strong indication that the penetration depth of North Atlantic Deep Water (NADW) was reduced during the LGM (Lippold et al. 2016). This coherent dual-tracer signal of AABW expansion and contraction over glacial cycles extends back to the middle Pleistocene, and appears to have developed alongside the intensification of ice ages about 1 million years ago (Pena and Goldstein 2014). But despite these apparent changes in water-mass geometry, the ratio of another pair of isotopes in sediments, ²³⁰Pa/²³¹Th, suggests that the volume transport of the Atlantic Meridional Overturning Circulation (AMOC) during the LGM was not necessarily very different than today (McManus et al. 2004; Ritz et al. 2013; Böhm et al. 2015).

Looking further back in time to the Pliocene and preceding eras, before the northern hemisphere ice age cycles began, CO₂ was as high, or higher than today. It would appear that the production of dense waters from polar oceans was also different under these warmer states, and changed during the onset of northern hemisphere glaciation 2.7 million years ago (Sigman et al. 2004, 2007), when δ¹³C records imply that Antarctic deep water formation weakened alongside a CO₂ drop (Hodell and Venz-Curtis 2006). Although tectonics may have played a role in these more ancient changes, it is possible that a decrease in Antarctic water formation was caused by the climatic effect of the CO₂ decrease, via its impact on surface ocean fluxes. Indeed, the sense of change is consistent with model simulations forced by anthropogenic CO₂ rise that have shown, after the rapid transient warming, increased production of Antarctic deep water in the equilibrated warm state (Yamamoto et al. 2015). Yet this is in apparent contradiction with the evidence from ice age cycles, in which Antarctic waters dominate the deep ocean in the colder glacial state.

Many prior studies have explored how these and other changes in steady-state ocean circulation might have resulted from CO₂, orbital and ice sheet changes, using General Circulation Models (GCMs). Most of these studies either used coupled ocean and atmosphere GCMs, forced by ‘time-slice’ combinations of radiation and topography (most commonly those of the mid-Holocene and LGM, e.g. Hewitt et al. 2003; Braconnot et al. 2007), or used ocean GCMs to which experimental changes of wind, heat, and freshwater fluxes were applied (e.g. Butzin et al. 2005; Toggweiler et al. 2006; Jansen 2017). These studies have identified and illustrated many potentially important mechanisms. However, when multiple models have been subjected to the same boundary conditions, they have not always produced the same results, because of differences in model physics and numerics (Otto-Bliesner et al. 2007; Muglia and Schmittner 2015), insufficient model equilibration (Marzocchi and Jansen 2017), or both. Furthermore, the heterogeneity of prior experiments makes it difficult to gain an appreciation for the relative importance—and interactions—of the various mechanisms over the full range of climatic states that have held sway over the past few million years. Additionally, though simplified models have illustrated how the interaction of ocean circulation with the marine ecosystem may have driven the ice age cycles of atmospheric CO₂ (Köhler et al. 2005; Bouttes et al. 2010; Hain et al. 2010; Brovkin et al. 2012), comprehensive models have thus far failed to do so without modifications (e.g. Simmons et al. 2016; Buchanan et al. 2016). Presumably, the mechanisms that drove the past changes will also be relevant for the long-term future of the deep ocean as it adjusts to the high CO₂ of the Anthropocene.

Here we aim to remedy the conceptual fragmentation of prior works by presenting a large number of equilibrium ocean states in a single coupled ocean–atmosphere GCM. This approach is inspired by the pioneering work of Manabe and Bryan (1985), who explored the response of an early ocean–atmosphere GCM to a range of six CO₂ forcings. The model we use is CM2Mc, a coarser-resolution configuration of the Geophysical Fluid Dynamics Laboratory’s ESM2M model (Dunne et al. 2012), which captures critical aspects of the Southern Ocean circulation (Russell et al. 2006). We simultaneously vary atmospheric CO₂, orbital configuration and terrestrial ice sheet size in a matrix of 40 experiments, mapping out the large-scale behaviour of the model over an unusually-wide parameter space. The model does not use flux corrections, and freely evolves to an equilibrium state determined only by the response of the internal model physics to the external forcing variables. Thus, in contrast with experiments that arbitrarily alter surface boundary fluxes in order to illustrate hypotheses, we show the model’s response to the external forcings without preconceptions. Although the model is subject to important biases, the simulation suite provides a comprehensive overview of how these three drivers could feasibly influence climate given Earth’s present tectonic configuration.

The simulations show a rich diversity of behaviour in the ocean, atmosphere and land surface, among which we briefly discuss general atmospheric aspects before focusing on the circulation of the global deep ocean. We also compare the results with some key paleoceanographic reconstructions from the LGM, showing encouraging consistency with physical observations as well as important biogeochemical differences. We conclude by proposing a simplifying principle, based on the relative density of polar waters, that makes testable predictions regarding the volume filling of the deep ocean under different climate states, and offer speculation regarding how this may have contributed to the late Cenozoic intensification of the ice ages.

2 Methods

2.1 Model

The simulations use the coupled ocean–atmosphere–ice-biogeochemistry model CM2Mc.v2. In brief, this includes: a finite volume atmospheric model, similar to that used in the GFDL ESM2 models (Dunne et al. 2012); MOM5, a non-Boussinesq ocean model with a fully-nonlinear equation of state, parameterizations for mesoscale (Griffies 1998) and submesoscale (Fox-Kemper et al. 2008) turbulence, and parameterizations for vertical mixing due to boundary layer turbulence (Large et al. 1994) and to the breaking of internal tides above rough topography (Simmons et al. 2004); a sea-ice module; a static land module; and a coupler to exchange fluxes between the components. It differs from the model described by Galbraith et al. (2011) due to the following small changes: the codebase was updated to the Geophysical Fluid Dynamics Laboratory Siena version; the albedo of snow-covered sea ice was increased to a more realistic value of 0.85 from 0.80; the bottom drag was increased from 0.001 to 0.002 (closer to its default value); the background vertical mixing was reduced from 0.1 to 0.05 cm² s⁻¹; additional cross-land mixing in Indonesia was included to account for unresolved flows between islands. These changes altered the pre-industrial control state slightly, but the overall simulation characteristics remain similar to the GFDL ESM2M in spite of the 1.5- to threefold reduction of horizontal resolution (nominally 3° in the ocean and 3° in the atmosphere, vs. 1° and 2°, respectively, in ESM2M). The ocean and atmosphere comprise 28 and 24 layers, respectively. The model does not use flux corrections, and conserves freshwater by routing precipitation over land to coastal ‘river’ discharge cells, including as ‘iceberg’ discharge where that precipitation falls as snow. Sea ice has a bulk salinity of 5 PSU, and its formation from seawater rejects brine in the top ocean layer. The high precision of the freshwater balance is evidenced by the small drift of whole ocean salinity, which changes by < 0.001 PSU over the last 5 centuries in all simulations.

The model does not apply flux corrections, and therefore has significant biases. Primary among these are a global cold bias in the pre-industrial control surface air temperature (12 °C, rather than the estimated pre-industrial average of roughly 14 °C), which is partly due to the increased sea ice albedo mentioned above, but which was deemed acceptable given that the global ocean average temperature bias is very small (+ 0.04 °C). The North Pacific has a particularly cold bias, which tends to weaken the stratification there and encourage more sea ice and stronger North Pacific Intermediate Water formation than observed. The Southern Ocean has a warm bias, though its hydrographic structure agrees well with observations (Galbraith et al. 2011). It should also be borne in mind that the ocean model does not resolve geostrophic turbulence and therefore parameterizes the effect of mesoscale eddies. This parameterization (Gent and McWilliams 1990; Griffies 1998) is an imperfect surrogate for eddy effects whose response to changes in the mean climate state is certainly biased (e.g. Munday et al. 2013). In addition, the model—like all global climate models—is incapable of correctly simulating the flow of dense waters along steep topographic slopes. As a result, convective chimneys, instead of downslope currents, are the dominant mode of dense water injection to depth, biasing the properties and sensitivity of deep water masses to some degree.

Water masses were ‘tagged’ in the surface layer in order to unambiguously identify volume contributions from different source regions. Three tags are defined: a Southern tracer south of 30°S, a North Atlantic tracer north of 30°N in the Atlantic, and a North Pacific tracer north of 30°N in the Pacific. In the surface layer, the tracers are set to 1 within their tagging region and to zero outside. In all other ocean layers, they are carried as passive tracers. We define AABW and NADW volumes as the global volume integral of the Southern and North Atlantic tracer concentrations, respectively, below 1 km depth. Because the model’s analogue of Antarctic Intermediate Water does not penetrate deeper than 1 km, restricting the analysis to depths > 1 km effectively excludes the volume contribution of intermediate and lighter waters. These definitions were chosen so as to be applicable under a wide range of climate states, and differ from classical water mass definitions based on density and basin criteria: all basins and all ocean depths > 1 km contribute to computed AABW and NADW volumes. Thus, AABW and NADW are broadly defined as the deep waters last ventilated in the Southern Ocean and North Atlantic, and their volumes essentially quantify the total volumetric influence of each surface high latitude region, irrespective of ventilation rates and pathways. The temperature, salinity and in situ density of the AABW and NADW end-members are calculated at 100 m-depth, averaged over the Southern Ocean or North Atlantic locations and years where the annual mean mixing layer depth exceeds 1 km.

An ‘ideal age’ tracer is defined as zero in the global surface layer, and increases in all other layers at a rate of 1 year per year. The biogeochemical module analyzed here is BLINGv0 (Galbraith et al. 2010). Radiocarbon is simulated by two ocean tracers, Dissolved Inorganic Carbon and Dissolved Organic Carbon (DI¹⁴C and DO¹⁴C), and compared to the non-radiogenic tracers DIC and DOC. These tracers are cycled implicitly through organic matter, assuming a fixed C:P of uptake, but the ¹⁴C tracers are not cycled through CaCO₃ for simplicity. Both DI¹⁴C and DO¹⁴C undergo radioactive decay. Air–sea fluxes are calculated assuming a CO₂ partial pressure of 270 µatm and a pre-industrial atmospheric Δ¹⁴C of 0‰. Because the biogeochemistry module uses fixed pre-industrial atmospheric CO₂ and ¹⁴CO₂, it does not include any air–sea exchange effect of variable pCO₂ (Bard 1988; Galbraith et al. 2015). Note also that the biogeochemical model simulates changes in dissolved nutrients and gases that are caused by the changes in ocean circulation only, thus excluding changes in dust fertilization or any other ecosystem change not driven directly by the orbital, CO₂ and ice sheet forcing impact on ocean circulation.

2.2 Simulations

The simulations populate a matrix that is defined by three external forcing variables (Table 1). First, atmospheric CO₂ was prescribed for each simulation, at levels ranging from 147 to 911 ppm. The levels from 270 (pre-industrial) to 911 are each 1.5-fold greater than the prior level, thus presenting a nearly-regular spacing of radiative forcing, while the four lowest (spanning the CO₂ range over recent glacial cycles) are spaced at half that distance. Although we prescribe only changes in CO₂, the model response would not differ significantly if forced by other greenhouse gases (CH₄, N₂O) at equivalent levels of radiative forcing.

Table 1 Overview of model simulations

Second, the two astronomical parameters were each set to two values, resulting in four different configurations (Fig. 1). Obliquity, i.e. the tilt of the rotational axis, alters the latitudinal distribution of radiation across the surface of the Earth, increasing annual incident radiation at high latitudes under high obliquity. Obliquity was set to either 22.0° or 24.5°, spanning the calculated range of the last 5 million years (Laskar et al. 2004). The precessional phase changes the seasonal distribution of insolation at any given location, without altering the total annual incident radiation, and is defined here as the angle between the Earth’s position during the northern hemisphere autumnal equinox and the perihelion (Held 1982). The precessional phase was set to the positions at which the insolation seasonalities are most and least severe: the most severe northern seasonality occurs at 270°, while the most severe southern seasonality occurs at 90°. It should be noted that these precessional phases will not necessarily produce the most extreme responses in climate metrics, given that the seasonal cycles of insolation and surface temperatures are not perfectly in phase (Kutzbach et al. 2008); nonetheless they provide strong and easily interpreted forcings. The eccentricity was maintained at a value of 0.03, which is close to the mean for the range of 0.000–0.058 spanned over the last 5 million years (Laskar et al. 2004) and nearly double the modern value of 0.01671.

Fig. 1

Orbital forcing configurations. Earth is shown at the closest and furthest approaches to the sun (perihelion and aphelion, respectively, not to scale). The tilt of the Earth’s rotational axis (obliquity) is set to either 22° or 24.5°. The precessional phase is defined by the orientation of the rotational axis relative to perihelion/aphelion; if the north pole points towards/away from the sun at perihelion/aphelion, boreal seasons are more intense than austral seasons (upward-pointing arrows)

Third, simulations were run using non-interactive land cover type and topography corresponding to either pre-industrial or ‘glacial’ conditions. For the latter, the surface topography was altered to reflect the geometry of terrestrial ice sheets as reconstructed for the LGM and land surface types were changed according to the Paleoclimate Modelling Intercomparison Project 3 (PMIP3) protocol, the Bering Strait was closed, coastal bathymetry and land-sea masks were altered to prevent ocean cells occurring underneath terrestrial ice sheets, and average ocean salinity was increased by 1 PSU [that which would be produced by a 110 m sea level fall, typical of glacial states, and ~ 0.1 PSU less than estimated for the LGM (Adkins et al. 2002)]. We also include four simulations at glacial conditions and CO₂ of 180 ppm, but for which the land surface type was kept at pre-industrial, thereby allowing a separation of the effects of ice sheet topography vs. albedo. Note that ice sheet sizes are held constant throughout any given simulation here, assuming all mass additions (by snow) are balanced by ablation (iceberg calving at river discharge points) and thereby ignoring the possible impacts of transient ice sheet growth and decay on freshwater inputs to the ocean. The latter impacts were previously explored using a small subset of these simulations together with freshwater forced (‘hosed’) experiments by Brown and Galbraith (2015) and Galbraith et al. (2016). We do not reduce ice sheet size or mean ocean salinity in the high CO₂ runs: only the roles of radiative and orbital forcings are explored within the 270–911 ppm CO₂ range.

All simulations were initialized from either an equilibrated pre-industrial control state or modern observations, and integrated until deep ocean salinity gradients approached a stable quasi-equilibrium. Some of the simulations were branched from 2000-year integrations with the same CO₂ and modern orbital configuration; because the ocean temperature was already somewhat equilibrated to the CO₂, these required a shorter subsequent integration with the four different orbital forcings (1200–1400 years). Similarly, the ‘180-glacial’ simulations were integrated for 3200 years with pre-industrial albedo, before swapping to the glacial land surface albedo. Total spinup times (starting from the pre-industrial control) were between 3100 and 6000 years (Table 1). The simulations were run using the Scinet general purpose cluster at the University of Toronto between 2013 and 2015, at a computational cost of 23 CPU hours per simulation year, requiring a total of approximately 3.6 million CPU hours. The analyses presented here are all averages of the final 100 years of the simulations. Because of the very long timescale of radiocarbon equilibration, the simulations all have substantial downward drift in the global average Δ¹⁴C of 0.2–1.3‰ per century (Table 1). The fact that the downward drifts are very similar for the ‘glacial’ and ‘interglacial’ states we analyse closely (boldface in Table 1) suggests that the difference between them is likely to be representative of the steady state. However, this assumption cannot be tested here, placing a caveat on the interpretation of the radiocarbon contrasts (Sect. 4.4.3).

It should be noted that these simulations were integrated for much longer than typical experiments published with full-complexity climate models, many of which have been run for less than 1000 years, as discussed by Marzocchi and Jansen (2017). The length of our simulations may be a major factor in explaining both their consistent behaviour and the relatively close agreement with observations.

2.3 Definition of ‘glacial’ and ‘interglacial’ states

Our simulations are designed to test model behaviour under external forcings that approximately bracket their variability over the past few million years. As such, none of the simulations is specifically designed to simulate the present interglacial (the Holocene) or the LGM. For the purpose of comparing interglacials and glacials in general, we take two representative sets of boundary conditions. The onset of full glacial conditions generally occurs under low obliquity, whereas interglacials generally occur under high obliquity (Huybers and Wunsch 2005). The chosen ‘glacial’ state is therefore a low obliquity simulation and the ‘interglacial’ is a high obliquity simulation, extrema that occurred most recently at 28 and 9 ka, respectively. Although deglaciations generally occur during the precessional phase with intense boreal seasonality (Cheng et al. 2009), precession is less strongly tied to glacial vs. interglacial states, given that the duration of the states frequently exceeds half a precessional cycle. We therefore choose the precessional state of 90° between boreal autumnal equinox and perihelion for both states, corresponding closely with both the present-day and the LGM. The ‘glacial’ state has an atmospheric CO₂ concentration of 180 ppm and full LGM ice sheets, altered bathymetry and a 1 PSU increase in mean ocean salinity. Where appropriate, we also compare observations with a pre-industrial control simulation, identical to the ‘interglacial’ simulation but using present-day orbital forcing (obliquity of 23.439°, 102.932° longitude between perihelion and boreal autumnal equinox).

3 Results

3.1 Global atmospheric response

Although this paper is focused on the ocean, we begin with a brief description of some general features of the global atmospheric response, since all oceanic changes are mediated by the atmospheric response to the external forcings.

Figure 2 shows the global pattern of surface air temperature change driven by the two orbital forcing variables of obliquity and precession, after full equilibration of the climate system, including the slow adjustment of the deep ocean and high latitude clouds and sea ice. The two patterns of change driven by the two types of forcing are very different, a factor which is not evident in the common simplification of orbital forcing to metrics such as peak summertime insolation at 65°N. Obliquity has a very strong impact on the zonal mean temperature profile: high latitudes are much cooler under low obliquity (i.e. a small tilt of Earth’s rotational axis), intensifying the equator-to-pole gradient by more than 3 °C. Precession produces a more complex spatial pattern of temperature change. Note that the magnitude of the precessional response depends on eccentricity (set here at an approximate Quaternary average of 0.03) and has therefore changed over time.

Fig. 2

Effect of orbital forcings on surface air temperatures. Each panel shows the surface air temperature anomaly produced by one of the two types of orbital forcing, calculated from the average of all simulations with pre-industrial ice sheets. a Low obliquity minus high obliquity. b Strong boreal seasons minus weak boreal seasons

Figure 3 shows the globally-averaged surface air temperature of each final state in the suite of simulations. CO₂ is the leading control of the global mean temperature but orbital forcing also has a significant impact. This orbital impact is due only to the spatial and temporal distribution of insolation, since the total integrated insolation is unaffected by either obliquity or precession. In general, the world is warmer when obliquity is high (orange) and when boreal seasonality is weak (downward-pointing triangles). Ice sheets exert a significant cooling at the global scale, most of which is attributed to the direct albedo impact rather than dynamical topographic effects (compare filled, outlined, and unfilled triangles in the 180 ppm experiments).

Fig. 3

Globally-averaged surface air temperature vs. CO₂. Each triangle shows the final state (100-year average) of one equilibrium simulation. Blue and orange symbols indicate low and high obliquity, respectively. Upward- and downward-pointing arrows indicate strong boreal and austral seasons, respectively, according to the precessional phase. Filled solid outlines indicate simulations using pre-industrial land surface cover with LGM ice sheet topography, bathymetry, and ocean salinity, while unfilled solid outlines replace the land surface type with an LGM estimate. Small letters indicate the ‘glacial’ (g) and ‘interglacial’ (ig) simulations discussed in the text

The global temperature response to the forcings is subject to model uncertainties, most notably the response of clouds. For reference, the equilibrium climate sensitivity of CM2Mc (the response to a CO₂ doubling from 270 to 540 ppm in the modern orbital configuration) is 3.6 K (not shown). This value is somewhat higher than the 3.1 K/doubling sensitivity of the one-degree ESM2M model (Frölicher et al. 2014). As shown by Fig. 3, the climate sensitivity increases with cooling of the background state (i.e. the slope is steeper at low CO₂), agreeing with the findings of Kutzbach et al. (2013), though this tendency may not be uniform among models and is sensitive to uncertain aspects of sea ice and clouds.

The global mean surface air temperature difference between the ‘glacial’ simulation and a pre-industrial (PI) control simulation is 5.4 K. This can be compared with PMIP2 models, which simulated a range of 3.6–5.7 K for the LGM to PI change (Braconnot et al. 2007). PMIP2 models generally show a larger LGM-PI change than the corresponding proxy-based estimates of 4.0 ± 0.8 K (Annan and Hargreaves 2013), although the latter might underestimate the true change (Telford et al. 2013; Ho and Laepple 2015). The relatively large glacial-interglacial contrast simulated by CM2Mc, compared to the PMIP2 LGM-PI average, can be partly explained by the fact that the envelope of orbital changes applied in our simulations is larger than the LGM-PI change (the orbit is more eccentric, and the obliquity change is larger), though the change in ocean temperatures appears to be well-simulated by the model, as discussed below.

3.2 Changes in water mass volumes and the AMOC

The experiments show a large range in the volumetric contribution of AABW to the deep (> 1 km) ocean (Fig. 4a), varying from 30 to 80%. The AABW volume shows a strong dependence on CO₂, with a striking nonlinearity: the AABW volume reaches a minimum at intermediate CO₂ (270–405 ppm) and increases toward both higher and lower CO₂. The effect of orbital forcing is also significant, and displays a further nonlinearity with CO₂. For CO₂ of 270 ppm and below, low obliquity (blue) increases the Antarctic ventilation volume, whereas only precession plays a role at higher CO₂, producing a greater fraction of AABW when boreal seasonality is weak (downward-pointing triangles). The inclusion of large LGM ice sheets also increases the AABW fraction, due to both the topographic effect and the ice sheet albedo (compare filled, outlined, and unfilled triangles in Fig. 4a). North Atlantic ventilation volumes are a mirror image of the Antarctic volumes (Fig. 4b), because North Pacific waters make up a small fraction of the global deep ocean in all simulations. Under low CO₂, low obliquity and pre-industrial ice sheets the North Pacific fraction reaches a maximum, overlapping with the North Atlantic fraction which is at its minimum. The North Pacific fraction decreases both with increasing CO₂ and with larger ice sheets (Fig. 4c).

Fig. 4

Global deep water volumes and AMOC. Panels show the simulated percentage of the deep (> 1 km) ocean volume occupied by waters originating in the a Antarctic, b North Atlantic and c North Pacific, and d the maximum AMOC transport at 40°N. All symbols are as defined in Fig. 3

The flow of water transported by the AMOC is relatively stable between 16 and 22 Sv for all simulations in which CO₂ is 270 ppm or above, and varies with both precession and obliquity (Fig. 4d). However, below 270 ppm the AMOC transports span a range of 3–22 Sv. Reduced AMOC transport under low CO₂ (150 ppm) was also simulated by Stouffer and Manabe (2003) using an earlier version of the GFDL coupled model. The weakest AMOC transports occur under low obliquity with pre-industrial ice sheets, whereas strong transports occur under high obliquity with LGM ice sheets. In two of the intermediate-strength simulations the AMOC spontaneously oscillates between strong and weak states, reminiscent of Dansgaard–Oeschger variability, as discussed in detail by Brown and Galbraith (2015). Note that, in all cases shown here, the freshwater input to ice sheets from precipitation is balanced by calving in the ocean, so that we do not include any possible transient effects of ice sheet mass balance on ocean circulation (e.g. Swingedouw et al. 2007).

3.3 Deep ocean ventilation rate and biogeochemistry

Like the AABW volume fraction, the simulated ventilation rate of the global deep ocean (quantified as the mean ideal age below 1 km) has a nonlinear dependence on CO₂ (Fig. 5a). The mean age is lowest (i.e. the ventilation is most rapid) under low CO₂, passes through a maximum at intermediate CO₂, and decreases somewhat at high CO₂. Orbital forcing has a large impact at low CO₂, and very little impact at high CO₂. LGM ice sheets increase the ventilation rate of the deep ocean significantly, mostly due to the topography effect (i.e. there is little difference between filled-outlined and unfilled triangles at 180 ppm).

Fig. 5

Deep ocean ventilation rates and biogeochemistry. a Mean ideal age deeper than 1 km. b Mean radiocarbon content (Δ¹⁴C, in per mil) deeper than 1 km. c Mean oxygen concentration deeper than 1 km. d Mean surface nitrate south of 40°S. Symbols are defined as in Fig. 3

Perhaps surprisingly, whole ocean radiocarbon does not follow the expectation given by ideal age (Fig. 5a, b). Ideal age decreases with increasing CO₂ from 405 to 911 ppm, but radiocarbon is relatively invariant across this range. Further, whereas ideal age decreases with declining CO₂ below 405 ppm, the simulations show a decrease of ¹⁴C—that is, an increase in the radiocarbon age. Hence, radiocarbon is not a straightforward monitor of ideal age (Schmittner 2003; DeVries and Primeau 2010), although it still contains valuable information, as discussed below.

Deep ocean oxygen is much more consistently linked to ideal age than ¹⁴C: elevated age co-occurs with low oxygen content, and vice-versa (Fig. 5a, c). Southern Ocean nitrate, which has been linked to the efficiency of the global biological carbon pump through its control of preformed nutrients (Sarmiento and Toggweiler 1984; Sigman and Boyle 2000), varies most strongly with obliquity at low CO₂, and generally decreases with increasing CO₂ (Fig. 5d). This is opposite to what would be expected from the argument of a Southern Ocean control of ice age CO₂ tied to a reduction of preformed nutrients via reduced vertical exchange (Sigman et al. 2010).

4 Discussion

Our results show large and nonlinear responses of AABW volume, AMOC and ocean ventilation to CO₂, overprinted with additional influences of orbital forcing and ice sheets. We first discuss the mechanisms by which the three types of forcing achieve these nonlinear responses in ocean circulation. Next we discuss two of the simulations in greater detail, most representative of ‘glacial’ and ‘interglacial’ states, in order to test the model behaviour against the well-documented LGM to Holocene change. We then briefly discuss the simulations of warmer climates, of relevance both to pre-Quaternary warm periods and to the long-term future of the Anthropocene. Finally, we identify an increasing sensitivity of the climate system to orbital forcing under decreasing CO₂ that may have contributed to the development of ice age cycles.

4.1 CO₂ and orbital influence

The simulated changes in AABW vs. NADW volume fractions can be predicted, to first order, by the density difference between the two formation regions (Fig. 6). At intermediate CO₂ (270–400 ppm), for which the AABW volume is minimal, the density at deep water formation regions of the North Atlantic tends to match or slightly exceed its Antarctic counterpart. At higher and lower CO₂, for which AABW dominates, the density of surface waters in AABW formation regions substantially exceeds that of NADW formation regions. The fact that similar volume fractions co-occur with similar density differences despite very different combinations of temperature and salinity (colours in Fig. 6) underlines the predictive strength of the density difference.

Fig. 6

Antarctic and North Atlantic volume contributions to the deep ocean as a function of the Antarctic–North Atlantic density difference. Simulations with LGM (pre-industrial) topography are shown as unfilled (filled) circles. Symbol size reflects atmospheric CO₂ (in ppm) while symbol colours reflect the Antarctic-North Atlantic a salinity and b temperature difference. Densities, temperatures and salinities are averaged at 100 m depth over the Southern Ocean or North Atlantic locations and years where deep (> 1 km) convection occurs

In contrast, the volume changes do not correlate with the strength of southern hemisphere westerlies, though they do correlate with the volume transport of the Antarctic Circumpolar Current (ACC, Fig. 7). This indicates that changing westerlies are not a main driving force behind the simulated changes in volume fractions. Rather, it is changes in surface buoyancy forcing, largely independent from the winds, that control the variations in water mass volumes. The large simulated changes in density reflect changes in both temperature and salinity, which we discuss in turn.

Fig. 7

AABW fraction vs. transport by the Antarctic Circumpolar Current (ACC) and vs. the strength of Southern Hemisphere westerlies. Symbols are defined as in Fig. 3. In general, the ACC transport is increased under low CO₂ due to the strong meridional density gradient induced by Antarctic sea ice export, and at high CO₂ due to the thermal density gradient. The mean surface velocity (in m s^− 1) of the southern westerlies does not show a significant relationship with the AABW fraction

The temperature ranges of waters in convective regions of the North Pacific and North Atlantic span 12 °C among the simulations, whereas the Antarctic convection sites remain close to freezing in most simulations, and only rise 5 °C above the freezing point in the highest CO₂ simulations (Fig. 8). Thus, at high CO₂, the more rapid warming of deep convection sites in the North Atlantic causes them to decrease in density more rapidly than their Antarctic counterparts, despite an increase of salinity. Meanwhile, with decreasing CO₂ below 270 ppm, Antarctic waters grow relatively denser despite North Atlantic cooling because of large salinity increases. Note that three low CO₂ (< 180 ppm), low obliquity experiments with pre-industrial ice sheets make an exception to the general North Atlantic trend: unusually warm and salty NADW (Fig. 8) results from an equatorward migration of the North Atlantic convection centre, following a sea ice expansion (not shown). Despite this large change in convective regions, the density difference still explains the relative volume fractions (Fig. 6).

Fig. 8

Surface properties in convective regions of the three polar oceans. For each simulation there is one value for each of the North Pacific, Antarctic and North Atlantic, calculated from 100-year time series. Warm/cool shades indicate high/low obliquity, triangle orientation indicates precessional phase as in Fig. 3, larger symbol size indicates higher atmospheric CO₂, solid fills indicate preindustrial topography and outlines indicate LGM topography and albedo. Temperatures and salinities, at 100 m depth, were averaged within each basin at the locations and years where the deepest convection occurred (mixing layer depth > 1 km or > 90% of its overall maximum). Under pre-industrial conditions, Antarctic and North Atlantic waters have similar densities, whereas they diverge under both high and low CO₂. The North Pacific has lower density than the Antarctic and North Atlantic under all simulated conditions. Small asterisks indicate the three simulations with low-latitude North Atlantic deep convection, discussed in the text

Salinity is influenced at the sea surface by the transport of freshwater in the atmosphere and by sea ice cycling (e.g., Durack et al. 2012; Goosse and Fichefet 1999). At steady state, these surface freshwater fluxes are balanced by the net freshwater transport by ocean currents (Nilsson et al. 2013; Cessi and Jones 2017). In the model, as observed (Durack et al. 2012), the high latitude Southern Ocean and North Pacific are net ‘freshening’, receiving more freshwater from precipitation than they lose to evaporation, while the Atlantic basin as a whole is net evaporative, due to atmospheric freshwater export. Thus, a stronger atmospheric hydrological cycle, as occurs under warmer planetary temperatures given the Clausius–Clapeyron relationship (Held and Soden 2006), freshens surface waters in the Southern Ocean and North Pacific but increases the Atlantic salinity (Fig. 9). In the model experiments, the net impact of lowering CO₂ on precipitation and evaporation is therefore to decrease the surface freshening of the Southern Ocean relative to the North Atlantic.

Fig. 9

Surface freshwater balance of the Atlantic and Southern Oceans. a Net atmospheric freshwater input to the Atlantic Ocean, including the Mediterranean and Arctic (note negative values indicate net evaporation). b Net atmospheric freshwater input to the Southern Ocean, south of 40°S. c Total freshwater transport by Antarctic sea ice. d Total freshwater input to the ocean south of 60°S, including all air–sea exchange and sea ice cycling. Symbols are defined as in Fig. 3. All quantities are annually-averaged

Even more important here are changes in the extent of southern hemisphere sea ice. As discussed in earlier studies (Keeling and Stephens 2001; Shin et al. 2003; Ferrari et al. 2014; Watson et al. 2015; Jansen and Nadeau 2016; Klockmann et al. 2016), low atmospheric temperatures encourage sea ice production and export, extracting freshwater around Antarctica (where freezing dominates) and releasing it further north (where melting dominates). Sea ice also acts as a barrier to precipitation, collecting snowfall and preventing it from entering the ocean until the ice beneath it melts. Because of the open geometry and strong zonal winds of the Southern Ocean, Ekman transport can potentially move large volumes of sea ice across long distances, unlike the North Atlantic where the basin shape limits the ability to export sea ice, causing the meridional transport to be smaller (Goosse and Fichefet 1999; Shin et al. 2003; Lee et al. 2017). Thus, when low CO₂ cools the global atmosphere, sea ice forms more rapidly at high latitudes of the Southern Ocean and is transported further north, causing the high latitude Southern Ocean to transition from a net-freshening to a net-salinifying domain (Fig. 9c, d). Low obliquity (blue) also increases the growth of sea ice, since it leads to colder temperatures at high latitudes.

Because the model does not use freshwater flux corrections, we can observe an additional impact of the Antarctic salinity increase: for a given amount of salt in the ocean, if deep waters formed in the Antarctic become saltier, deep waters formed in other regions (primarily the North Atlantic) must become fresher. Amplified sea ice export salinifies the waters close to Antarctica but freshens the lower-latitude intermediate and mode waters—whose northward transport ultimately feeds North Atlantic sinking. Thus, when Antarctic waters become saltier and denser, and consequently expand their domain in the deep ocean, they force a reduction in the salinity and therefore density of North Atlantic waters, amplifying the impact of the initial density gain. The volume of Antarctic sea ice shows a strong, approximately linear relationship with the Antarctic-North Atlantic salinity difference (Fig. 10).

Fig. 10

The inter-basin salinity gradient and Southern Ocean sea ice. Difference in salinity between Antarctic and North Atlantic deep convection sites vs. the total annually-averaged Antarctic sea ice volume. A larger Antarctic sea ice volume is accompanied by larger northward freshwater transport, increasing the salinity of AABW and indirectly decreasing the salinity of NADW

The salt balance can also explain part of the observed changes in the strength of the AMOC at low CO₂. The freshening of the North Atlantic surface brings the Atlantic circulation closer to a bifurcation point where the AMOC is bi-stable: the poleward advection of salty subtropical waters becomes prone to interruption, allowing the subpolar gyre to develop a strong halocline and impeding the formation of dense waters (Stommel 1961). This behaviour permits unforced oscillations of the AMOC to develop in two of the simulations at 180 ppm CO₂ as discussed in detail by Brown and Galbraith (2015), and provides a potential explanation for the fact that the Dansgaard–Oeschger “bi-polar seesaw” variability does not occur under the high CO₂ of interglacial periods (McManus et al. 1999; Rahmstorf 2002; Buizert and Schmittner 2015; Kawamura et al. 2017). However, Dansgaard–Oeschger variability is not observed under full glacial maxima either (McManus et al. 1999; Kawamura et al. 2017), a fact that can be explained by the simulated ice sheet effects.

4.2 Effects of Laurentide ice sheet topography and albedo

Although the AMOC is weakened by the inter-basin salt balance under low CO₂ and low obliquity, this tendency is counteracted in all simulations that include the full LGM ice sheets. Figure 11a confirms that low CO₂ causes the mixed layer depths to shoal in the northeast Atlantic, as a strengthened halocline impedes convection. As previously noted by other authors, steering of the atmospheric circulation by a large Laurentide ice sheet (Fig. 11b) strengthens the zonality of the westerly jet along the southern margin of the ice sheet and central Atlantic (Roe and Lindzen 2001; Li and Battisti 2008). However, over the Labrador Sea, the ice sheet topography generates a strong northwesterly flow, accompanied by a very strong cooling (Broccoli and Manabe 1987) that is amplified by the ice sheet albedo effect (Fig. 11c). An anomalous northward wind stress occurs in the northeast Atlantic. The wind shift and cooling appear to drive the AMOC strengthening shown in Fig. 4d, as has previously been observed in other models (Pausata et al. 2011; Zhang et al. 2014; Zhu et al. 2014; Gong et al. 2015; Muglia and Schmittner 2015; Klockmann et al. 2016; Sherriff-Tadano et al. 2017). When compared with the influences of CO₂ and orbital forcing, this suggests that Dansgaard–Oeschger-like AMOC disruptions would be most likely to have occurred in a ‘sweet spot’ of low CO_2, low obliquity and a relatively small (low-elevation) Laurentide ice sheet. We therefore hypothesize that, under these conditions, the AMOC was more sensitive to variations in freshwater forcing arising from tropical climate variability or freshwater discharge from ice sheets (Schmittner and Clement 2002).

Fig. 11

CO₂ and ice sheet impacts on convection and surface winds in the North Atlantic. Panels show the glacial-interglacial change in mixed layer depth (m, shaded) and winds (arrows) averaged over the 4 orbital configurations, due to individual forcing aspects (a–c) and due to all (d). a Influence of low CO₂ (180 ppm minus 270 ppm). b Influence of ice sheet albedo (LGM albedo minus PI albedo). c Influence of ice sheet topography (LGM topography minus PI topography). d Influence of CO₂, ice sheet topography and albedo. Low CO₂ impedes deep convection in the northeast Atlantic, which we attribute to the inter-basin salt balance, while the ice sheet topography causes deep convection to greatly increase in the Labrador Sea in response to the wind change, an effect that is accentuated by the albedo change

Counter intuitively, the ice sheet-driven AMOC strengthening does not increase the volume contribution of NADW to the global ocean, but instead increases the contribution of AABW. Whereas low CO₂ has mixed effects on the temperature of NADW (Fig. 12a), the intense convection under extremely cold atmospheric conditions induced by the Laurentide ice sheet cools NADW dramatically (Fig. 12b, c). The cold NADW flows at < 2.5 km depths through the Atlantic and feeds into waters upwelling in the ACC, similar to the mechanism proposed by Crowley (1992). A further contribution to cooling, attributed to the ice sheet albedo, occurs in the upper layers of the Pacific ocean and is also transmitted to the Southern Ocean (Fig. 12c). The consequent cooling of the near-surface Southern Ocean increases the production of sea ice by reducing the heat loss required for winter ice growth (Martinson 1990) and increases the northward export of Antarctic sea ice by translating the melting regions northward. Thus, like low CO₂ and low obliquity, a large Laurentide ice sheet increases Antarctic salinity and ultimately the density and volume of AABW (Figs. 4, 8, 9). This mechanism is similar to previous suggestions that colder NADW may cool the Southern Ocean remotely (Gildor et al. 2002) and that it may increase Antarctic salinity by changing the melting rates of land-based glaciers (Adkins 2013). The simulations here confirm the importance of the remote cooling effect, but its impact on salinity occurs via Antarctic sea ice formation and export.

Fig. 12

CO₂ and ice sheet impacts on global ocean temperature and salinity. Panels show the change in temperature (top) and salinity (bottom) averaged over the 4 orbital configurations, due to individual forcing aspects (a–c) and due to all (d), as defined in Fig. 11, along sections through the Atlantic and Pacific oceans. The sections are zonal averages from 50°W to the prime meridian (Atlantic) and from 180°W to 140°W (Pacific). The increases in AABW salinity, in response to all three forcings, are primarily attributable to increases of sea ice formation and transport in the Southern Ocean

The overall result is a decoupling between meridional volume transport and water mass geometry, showing that a strong glacial AMOC can feasibly coexist with, and even contribute to, an expanded AABW. Importantly, this decoupling shows that variations in the AABW and NADW steady-state circulations cannot be universally reduced to an overturning see-saw: as we discuss further below, a large Laurentide ice sheet size leads to a simultaneous increase in the formation rates and volume transports of both NADW and AABW.

4.3 The unifying role of relative polar ocean densities

As we have shown, the general trends in deep ocean circulation can be rationalized by the influence of surface buoyancy fluxes on the surface density difference between the Antarctic and North Atlantic dense water formation regions. In contrast, the model does not show significant changes in deep ocean ventilation in response to changes in southern westerly winds, in disagreement with some prior suggestions (e.g. Toggweiler et al. 2006; Anderson et al. 2009). The unifying role of relative density variations encapsulates the changes brought about by all three climate forcings, despite difference in the mechanistic pathways by which density variations occur. In general, the greater the density of AABW relative to NADW, the more AABW dominates the ocean’s deeper layers. Our results also show that trends in volume contributions do not simply reflect a see-saw between northern and southern overturning, since the AMOC does not weaken at high CO₂ despite the shrinking NADW fraction, and since the AMOC and NADW fraction respond differently to the inclusion of a large Laurentide ice sheet.

These results are subject to important caveats: climate models are incapable of simulating many aspects of dense water formation and sinking, and poorly represent the mechanisms that consume dense waters at depth. However, we expect these biases to affect ventilation rates and pathways more strongly than the global volumetric contributions diagnosed here. We also expect the link between volumetric contributions and the inter-basin density competition to be relatively robust to biases in the location and mechanism of dense water formation. Indeed, the obtained relationship between volumes and formation-site densities holds across the vast parameter space explored (Fig. 6), despite major variations in the locations of deep water formation, in the global temperature and salinity distributions and in the water-mass transport rates. Additionally, the relationship echoes earlier results obtained with simplified ocean models (Stocker et al. 1992; Fichefet et al. 1994) and idealized eddy-resolving numerical experiments (Wolfe and Cessi 2009, 2010, 2011). In particular, the volumetric expansion of NADW as its formation density approaches that of AABW appears to be consistent with deeper NADW penetration for a larger window of shared surface densities between the North Atlantic and the Southern Ocean, as proposed by Wolfe and Cessi (2009, 2010, 2011).

4.4 ‘Glacial’ vs. ‘interglacial’ states

4.4.1 Glacial water mass characteristics

We first compare the ‘glacial’ ocean temperature and salinity to porewater and ice core-derived reconstructions of the LGM ocean (Fig. 13). The glacial simulation reproduces the remarkable switch from a thermally-stratified to a salinity-stratified deep ocean identified by Adkins et al. (2002). The simulated glacial salinity between 3500 and 5500 m in the Pacific Ocean is 36.2, indistinguishable from the LGM salinity of 36.1 ± 0.1 reconstructed from porewater extractions at five sites spanning the same depth range (Insua et al. 2014). We are not aware of any other comprehensive climate model that has previously reproduced the reconstructed salinities of the glacial ocean with comparable accuracy (e.g. see Fig. S7 of Bopp et al. 2017). The simulated salinity at Shona Rise is only 36.4, rather than the reconstructed 37.1 ± 0.2, but the robustness of the Shona Rise reconstruction has been questioned (Wunsch 2016) so we are unsure whether or not this discrepancy is significant. The model also simulates a change in Mean Ocean Temperature of 2.47 °C between the pre-industrial control and the ‘glacial’ state, which is indistinguishable from the recent estimate of 2.57 ± 0.24 ºC derived from noble gas ratios in Antarctic ice (Bereiter et al. 2018). Thus, we are unable to identify a significant bias in the model simulation of the glacial-interglacial change in deep ocean temperature and salinity.

Fig. 13

Simulated and reconstructed deep ocean temperature and salinity. Outlined circles show glacial reconstructions in blue and modern values at the same locations in red. Filled circles show the simulated glacial and interglacial temperatures and salinities at the same locations. a Temperature and salinity reconstructions of Adkins et al. (2002). b Salinity reconstructions of Insua et al. (2014). The model captures the basic transition from the modern, thermally-stratified deep ocean to the glacial, salinity-stratified deep ocean. The Shona Rise site, for which the salinity reconstruction has been questioned as discussed in the text, is identified by parentheses

Figure 14 shows a global slice through the ‘glacial’ and ‘interglacial’ oceans. The deep ocean of the interglacial state is stably stratified in temperature and, in the Atlantic, unstably stratified in salinity. By contrast, the glacial ocean is essentially isothermal beneath 1.5 km, and the salinification of glacial Antarctic waters reverses deep Atlantic vertical salinity gradient (Jansen 2017). There is also a large contrast in the depth to which NADW penetrates (Fig. 14c). In the interglacial case, the full water column of the North Atlantic is filled by NADW. This remains true in the glacial state for the Nordic and Labrador seas, but further south the contribution of NADW dwindles rapidly below 2 km depth. This simulated contraction of NADW co-occurs with the volumetric expansion of underlying Antarctic waters, and is consistent with the many paleoceanographic observations for shoaled NADW in the glacial North Atlantic (Curry and Oppo 2005; Adkins 2013; Lippold et al. 2016), although it does not preclude the possibility that part of the apparent shoaling reflects a change in the remineralized carbon content of AABW (Gebbie 2014; Howe et al. 2016). We also note the interesting possibility that, under glacial conditions, the Labrador Sea produced dense bottom waters that trickled downslope analogous to the modern Antarctic (Keigwin and Swift 2017), and which the model is incapable of simulating. The reduced input of NADW to dense waters upwelling in the Southern Ocean leads to very little contribution of North Atlantic waters to the Indo-Pacific in the glacial state, arguing against the idea of North-Atlantic ventilation of the southern source waters (Kwon et al. 2012).

Fig. 14

Sections through the Atlantic and Pacific oceans of (left) ‘interglacial’ and (right) ‘glacial’ simulations. The sections are zonal averages from 50°W to the prime meridian (Atlantic) and from 180°W to 140°W (Pacific)

The glacial vs. interglacial contrast is consistent with theoretical suggestions (Ferrari et al. 2014; Watson et al. 2015; Jansen 2017; Marzocchi and Jansen 2017) that Southern Ocean surface buoyancy fluxes explain the expansion of Antarctic deep waters in the global ocean, as detailed above. However there is also a key difference with Ferrari et al. (2014), who elegantly postulated that the isopycnal surface connecting with the zero-buoyancy-flux line at the Southern Ocean surface would migrate upwards under the glacial state, following a simple geometric constraint. In our simulations, the corresponding isopycnal surface does not move appreciably (Fig. 15), which we attribute to the impact of southern subpolar gyres on the Southern Ocean density structure. The gyre circulations act to homogenize density south of the ACC. In the interglacial state, sea ice cycling and surface buoyancy forcing are moderate and the gyres maintain a roughly uniform zonal mean upper-ocean density south of 60°S (Fig. 15a, c, e). As a result, the isopycnal connecting with the zero-buoyancy-flux line does not dive to the deep ocean until it reaches about 60°S, near the northern edge of the gyres (Fig. 15c, e). In the glacial, the zero-buoyancy-flux line shifts to a more northern location (~ 60°S; Fig. 15b), as predicted by Ferrari et al. (2014), but since this simply moves to the northern edge of the gyre domains, the latitude at which the dividing isopycnal plunges to depth does not change appreciably (Fig. 15d, f). The absence of a strong upward displacement of the dividing isopycnal in the model is robust using different measures of density (not shown) and considering seasonal effects (Fig. 15c–f, black and grey). Thus, we do not find a straightforward relationship between the latitude of the zero-buoyancy-flux line and the depth of its matching isopycnal in the deep ocean.

Fig. 15

The zero-buoyancy-flux line and its connection with the deep ocean under (a, c, e) interglacial and (b, d, f) glacial conditions. The annual mean potential density referenced to 1000 dbars (σ₁) is shown (a,b) at the surface of the Southern Ocean and as a zonal mean within the c, d Indo-Pacific and e, f Atlantic basins. The black curves in a, b mark the southernmost latitude of zero annual mean surface buoyancy flux. The (black) zonal annual and (grey) zonal austral-winter means of σ₁ along this zero-buoyancy-flux line are contoured in c–f. Key features include: (1) upper-ocean density does not increase monotonically with latitude; (2) the zero-buoyancy-flux line crosses a significant range of surface densities, particularly in the interglacial state; (3) the depth of the mean dividing isopycnal at low latitudes is comparable in the interglacial and glacial. These features are robust to different measures of density and result mainly from the presence of subpolar gyres

4.4.2 Glacial deep ocean ventilation and diapycnal mixing

The ventilation rate of the deep ocean can significantly alter carbon storage by the biological soft tissue pump (Sarmiento and Toggweiler 1984; Toggweiler 1999; Sigman et al. 2010; Skinner et al. 2010). If deep waters communicate more slowly with atmosphere-equilibrated surface waters, they can accumulate more carbon from sinking organic detritus, thereby lowering atmospheric CO₂. Indeed, as shown by Eggleston and Galbraith (2018), the soft tissue pump storage in CM2Mc is well-approximated by a linear product of the global ventilation age and the global export production. Many have argued that the low CO₂ of glacial periods was due, at least in part, to sluggish deep ocean ventilation (Toggweiler 1999; Galbraith et al. 2007; Sigman et al. 2010; Skinner et al. 2010; Menviel et al. 2017).

The model simulations do not accord with this picture of a sluggish deep ocean. In the glacial state, although most of the global ocean is ventilated from the Antarctic, it actually does so at a rate that exceeds the combined input of North Atlantic plus Antarctic waters under the interglacial state, leading to very low ideal ages in the deep Southern Ocean (Fig. 14).

This remarkably rapid ventilation arises as a consequence of the strong Antarctic density gradients and the model’s diapycnal mixing parameterization (Simmons et al. 2004). The large density excess of Antarctic-sourced bottom waters relative to other deep water masses increases the deep stratification, and the mixing parameterization assumes that diffusive buoyancy fluxes, driven by tides, are proportional to the local bottom stratification. This stratification-dependence allows the interior buoyancy supply to bottom waters, and their resultant upwelling rate, to match the strong surface buoyancy loss and convective downwelling occurring around Antarctica.

The maintenance of a rapid circulation across large glacial density contrasts thus hinges on a roughly threefold increase of the deep mixing energy supplied implicitly by tides. On the one hand, this does happen to be compatible with stronger tidal energy input to the deep sea under lower sea level, as might be expected from the fact that dissipation on shelves was less due to sea level fall during the glacial (Wunsch 2003), and as simulated with tidal models of the LGM (Schmittner et al. 2015). Sedimentological features have suggested that fast bottom flows existed in the abyssal ocean during the glacial (Hall et al. 2001; Kienast et al. 2007; Krueger et al. 2012) which could be consistent with a strong AABW source. On the other hand, the magnitude, distribution and stratification-dependence of the parameterized glacial tidal mixing in the model are certainly incorrect to a significant degree (MacKinnon et al. 2017).

In summary, a strong interior buoyancy source is required to balance the large surface buoyancy loss in the glacial Antarctic freezing zone at steady state, but the representation of subsurface ocean mixing in this and other climate models remains a major source of uncertainty for how and where this occurs, and for the consequent relationship between surface forcing and the deep ocean ventilation. Further work is needed to evaluate the assumptions behind the ocean mixing parameterizations, and the simulated ventilation rates should be viewed with caution. Nonetheless, we can compare the simulated radiocarbon with observations to evaluate the modelled ventilation rate.

4.4.3 Glacial radiocarbon simulation

One of the main pillars of observational support for sluggish ventilation of the deep glacial ocean has been the apparently low radiocarbon content of benthic foraminifera in glacial-aged sediments, when compared to the contemporary atmosphere (Galbraith et al. 2007; Skinner et al. 2010; Burke and Robinson 2012; Sarnthein et al. 2013). Radiocarbon is produced exclusively in the atmosphere, and is supplied to the ocean through air–sea exchange, so the benthic-atmosphere radiocarbon gradient is typically interpreted as representing the ‘ventilation age’ of bottom waters, often linked to the rate of ocean overturning. However, as previously shown, large changes in radiocarbon could have also occurred due to changes in air–sea flux and ocean circulation patterns, rather than overturning rate (e.g. Schmittner 2003; Butzin et al. 2005; Devries and Primeau 2010). Our model simulations show that, under a broad range of forcings relevant to the late Quaternary, changes in air–sea flux and ocean circulation patterns place the dominant controls on deep sea radiocarbon, with only a secondary imprint of changes in the circulation rate, quantified here by the ideal age.

The air–sea exchange timescale for radiocarbon is very slow, requiring multiple years of exposure to approach saturation. Because deep water upwelling and subduction can occur on time scales shorter than this, the exposure history of waters at the ocean surface plays a role in determining the initial, pre-formed radiocarbon activity of sinking waters. In the pre-industrial ocean, NADW had high radiocarbon activity due to the long exposure of its source waters in the subtropical Atlantic, whereas the Southern Ocean waters were relatively poorly equilibrated. In the glacial circulation, the input of radiocarbon-rich NADW to the AABW formation regions is diminished, because it is insufficiently dense to upwell at high latitudes of the Southern Ocean. In addition, the Southern Ocean is persistently sea-ice covered, blocking local air–sea exchange. Both these factors restrict the supply of atmospherically-produced ¹⁴C to waters forming AABW, further increasing the pre-formed age of AABW. In addition, the volumetric increase of “old” AABW, at the expense of “young” NADW, drives a whole-ocean shift to higher pre-formed age (Schmittner 2003).

The overall result of the changes in ocean circulation and sea ice is an increase in the pre-formed ventilation age of the ocean by 440 years between the glacial and pre-industrial states, as well as can be determined given the remaining drift in Δ¹⁴C (cf Sect. 2.2). In other words, if the ideal age of the ocean had been exactly the same between the two states, the change in pre-formed reservoir ages would have made the ocean more depleted in radiocarbon by an amount equivalent to roughly 440 years of decay. When this is combined with the pCO₂ effect on air–sea equilibration of radiocarbon, expected to be approximately 250 years for the LGM (Galbraith et al. 2015), it implies an apparent aging of the glacial ocean by 690 years relative to the pre-industrial. This is very similar to the whole-ocean average radiocarbon age increase of 600 years reported by Sarnthein et al. (2013) and 700 years reported by Skinner et al. (2015).

Thus, comparison of the simulated radiocarbon with the global compilations of radiocarbon reconstructions suggests that the true ideal age of the LGM ocean could have been similar to the interglacial, with the apparent difference in ¹⁴C age caused by a reduction of the pre-formed ¹⁴C activity of AABW, its volumetric dominance of the glacial ocean, and the pCO₂ effect. The fact that the simulated ideal age under the glacial state is roughly 200 years younger than that of the interglacial state (Figs. 5 and 14), rather than the same (as suggested by the radiocarbon), would be consistent with a significant overestimate of tidal mixing in the glacial simulation, as suggested above in Sect. 4.4.2. An additional, non-exclusive possibility is that a significant decrease in ventilation occurred following the obliquity minimum at 28 ka, so that the LGM was not directly equivalent to the simulated ‘glacial’ state, as discussed further below.

4.4.4 Glacial nitrate, deep ocean O₂ and carbon storage

As shown above, the glacial simulation agrees well with the reconstructed temperature, salinity and water mass proxies for the LGM. However, it simulates high summertime NO₃ concentrations at the Southern Ocean surface, and an oxygen-rich deep ocean (Fig. 5). In contrast, nitrogen isotopic evidence has been interpreted as showing that the LGM Southern Ocean surface was almost completely devoid of NO₃ in summer (François et al. 1997; Martínez-García et al. 2014; Wang et al. 2017), and the LGM deep ocean appears to have had significantly lower O₂ (Bradtmiller et al. 2010; Galbraith and Jaccard 2015; Hoogakker et al. 2015; Korff et al. 2016) as expected given greater storage of respired carbon by the biological soft tissue pump. Recent evidence even suggests that O₂ was quite low at 3 km depth in the LGM Southern Ocean (Jaccard et al. 2016), and that this depletion may have extended to relatively shallow depths (Lu et al. 2016). It seems possible that the Southern Ocean nitrogen isotopic records could have been misinterpreted if, for example, there had been isotopic enrichment due to denitrification in the O₂-depleted glacial deep ocean. However, the diversity of oxygen proxies makes the glacial oxygen depletion likely very robust: the model is clearly missing something. Buchanan et al. (2016) found a similar increase in deep ocean O₂ in a well-equilibrated simulation with the CSIRO coupled model under glacial boundary conditions, suggesting one or more missing elements may be common in the response of existing models to LGM radiative and ice sheet forcing. We propose three non-exclusive hypotheses for what might cause these models to grossly overestimate the glacial deep ocean oxygen content.

First, as argued in the prior section, it appears that the simulated glacial deep ocean ideal age is roughly 200 years too young, when compared to the LGM. A 200-year bias in ventilation age would suggest an underestimate of O₂ utilization in the deep sea of 60–80 µM, according to a regression of simulated age vs. O₂ utilization (Eggleston and Galbraith 2018), which would bring the simulated glacial O₂ back down to interglacial levels. Others have argued that ventilation of the LGM deep ocean was significantly slower than PI, based on model simulations forced with varying freshwater inputs (Bopp et al. 2017; Menviel et al. 2017), which would clearly help to achieve low glacial O₂ concentrations; however it is not obvious that a glacial steady-state deep ocean could have been substantially older than PI without violating the observed constraints on deep ocean radiocarbon and salinity.

It is conceivable that the LGM deep ocean diverged significantly from a steady state: the obliquity minimum, more directly comparable to our ‘glacial’ simulation, occurred at 28 ka, and it is possible that the gradually increasing axial tilt over the ensuing millennia reduced sea ice freshwater transport, driving a freshening trend in the Southern Ocean that truncated AABW formation during the LGM. A related possibility is that the advance of the Antarctic ice sheet to the shelf margins eliminated important coastal AABW formation regions (Paillard and Parrenin 2004), a mechanism which climate models cannot yet capture. An LGM decrease of AABW formation could perhaps explain the otherwise-puzzling drop of deep ocean radiocarbon that appears to occur between 25 and 20 ka (e.g. Burke and Robinson 2012, their Fig. 2b). In this case, the 200-year lower-than-PI age of our ‘glacial’ simulation may not actually be a model bias, but simply correspond more closely to 28 ka than to the LGM. However, regardless of whether the 200-year LGM ventilation age discrepancy is a physical model inaccuracy or reflects a non-steady-state LGM ocean, it is insufficient to explain why the model does not simulate the low deep ocean O₂ indicated by multiple proxy records.

Second, in addition to the ventilation age bias, it is possible that the absence of slope convection in the model contributes to excessive oxygen uptake in the glacial state relative to the interglacial state. CM2Mc can only form deep waters through vertical convection, in grid cells that are more than 100 km across, and does not simulate ice-shelf interactions. This mode of ventilation ‘aerates’ the whole water column, rapidly pumping oxygen into the deep sea. In contrast, most of the production of AABW in the modern ocean occurs via cascading plumes of dense shelf waters down the continental slope (Orsi et al. 1999; Heuzé et al. 2013), which climate models are incapable of simulating. This bias likely introduces an error in the simulation of the modern Southern Ocean, but could greatly exaggerate the oxygen input to the deep sea under glacial conditions, under which the rapid production of dense brines produces very large, persistent polynyas, mostly in the Weddell Sea. The interpretation that the glacial Southern Ocean was characterized by low rates of near-surface vertical mixing, based on sedimentary observations (Francois et al. 1997; Sigman et al. 2010), would certainly suggest that the model’s rapid vertical aeration of the water column is incorrect, and may indicate additional biases in the simulation of mixed layer dynamics.

Third, it is likely that significant changes in oxygen consumption occurred via ecosystem processes not included here. One likely candidate is iron fertilization, which certainly contributed to Southern Ocean NO₃ drawdown to some degree (Martínez-García et al. 2014), though proxy evidence for low Antarctic export production precludes this being the sole explanation of widespread NO₃ depletion (Sigman et al. 2010; Jaccard et al. 2013; Wang et al. 2017). Deepening of the remineralization profile, as would have resulted from slow ecosystem metabolism under low temperatures (Matsumoto et al. 2007), could have contributed to low deep ocean oxygen concentrations, though again it seems unlikely to explain the surface nitrate depletion. Simulations previously made by Buchanan et al. (2016) and Galbraith and Jaccard (2015) suggest that the combination of increased export production and deeper remineralization (together with the 200-year age underestimate) could possibly close the gap between the CM2Mc glacial state and the existing LGM oxygen observations. Testing these hypotheses would benefit from improved quantitative estimates of dissolved oxygen and export production, and from future model simulations that include additional tracers (e.g. δ¹³C, εNd) that can be compared directly with observations.

4.5 Deep ocean circulation in a warmer world

The discussion thus far has focused on the glacial-interglacial contrast, given the relatively strong observational constraints available. We now point out a few salient features of the high CO₂ simulations.

As discussed above, minima of the AABW volume contribution and ventilation rate occur at approximately present-day pCO₂ (Figs. 4, 5). The lowest AABW and highest NADW fractions of any simulation occur at 405 ppm, with high obliquity and strong boreal seasonality. This configuration—interglacial CO₂, high obliquity and strong boreal seasons—occurred most recently in the early Holocene (ca. 9 ka), when NADW formation may have been stronger based on both ocean circulation proxies (Hoogakker et al. 2011) and surface water characteristics (Thornalley et al. 2009; Marcott et al. 2013). The NADW-strengthening orbital forcing configuration would have been particularly strong during MIS 5e due to the greater eccentricity at that time (> 0.04), with a maximum effect near peak obliquity (131 ka) and strongest boreal seasons (127 ka), implying that NADW formation and volume would have been even larger then. This model-based prediction coincides remarkably well with an apparent stagnation event recorded in Southern Ocean sediments dated at 127 ka and interpreted as a brief cessation of AABW formation (Hayes et al. 2014). We therefore suggest the possibility that this stagnation event was enabled by a ‘perfect storm’ of astronomical configuration and relatively high CO₂, though it may also have been pushed over the edge by glacial melting as proposed by Hayes et al. (2014).

As pCO₂ rises above present-day levels, the AABW volume fraction once again increases. Although this nonlinearity is surprising, it would appear to agree with the observational inferences of dominant AABW under the high CO₂ of the Pliocene (Hodell and Venz-Curtis 2006; Bell et al. 2015). We stress that this is an equilibrium response to warming, and is very different from the transient response to rapid (e.g. anthropogenic) warming, under which AABW formation ceases due to freshwater stratification (de Lavergne et al. 2014).

We attribute the return of AABW dominance at high temperature to the relatively rapid warming of North Atlantic waters, which leads to AABW once again gaining the upper hand of density (Figs. 4, 6, 8), aided by the nonlinear dependence of density on temperature (De Boer et al. 2007). In the thermally-stratified deep ocean, precession becomes more important than obliquity, generating a greater contribution of NADW when boreal winters are particularly cold. As described for the low pCO₂ trends, the increased ventilation follows from the increase in the surface density difference between Antarctic and North Atlantic deep water formation sites. By promoting convective instabilities around Antarctica (Yamamoto et al. 2015) and increasing the stratification-dependent buoyancy fluxes in the interior, the larger density difference enhances the formation and consumption rate of abyssal waters, lowering ideal age once again (Fig. 5). We note that in a warm world with no Antarctic ice sheet (e.g. the Eocene), much higher Southern Ocean temperatures may have led to a different relationship with CO₂ than that shown here (Lefebvre et al. 2012).

4.6 CO₂-dependent sensitivity to astronomical forcing

One of the mysteries of the late Cenozoic is why the response of the climate system to astronomical forcing intensified from the Pliocene into the Pleistocene (Lisiecki and Raymo 2007). A number of explanations focus on the direct influence of orbital forcings on the Laurentide ice sheet (Maslin et al. 1998) and Antarctic ice sheet (Raymo et al. 2006). Our model simulations provide a different perspective, illustrating the degree to which an inherent increase in Earth system sensitivity to orbital forcing could have contributed to this shift.

Figure 16 shows that the response of the climate system to a given orbital forcing change is strongly amplified at low CO₂. For the global surface air temperature, this involves a fourfold increase in the combined impact of obliquity and precession, increasing from 0.6 to 2.5 °C as the CO₂ decreases from 911 to 147 ppm (Fig. 16a). The reason underlying the CO₂-dependence of orbital sensitivity is that albedo feedbacks become stronger in a colder world, as sea ice and snow cover expand, increasing the leverage of orbital changes on global surface air temperature (Fig. 16b). The sensitivity to obliquity is particularly responsive to CO₂ (Fig. 3, compare orange vs. blue symbols), which we attribute primarily to its important role as a modulator of sea ice extent. The impact of obliquity on sea ice is only significant when high latitude waters are cold enough for extensive sea ice to develop. Replacing pre-industrial with glacial ice sheets slightly increases this sensitivity, by further increasing the potential response of albedo.

Fig. 16

Enhanced sensitivity to orbital forcing at low atmospheric CO₂. a Each symbol shows the range of globally-averaged surface air temperatures among the 4 orbital configurations under each CO₂ concentration. Unfilled circles indicate glacial ice sheet topography and albedo, filled circles indicate pre-industrial ice sheet topography and albedo. b As in panel a, but for the range of planetary albedo, i.e. the ratio of outgoing to incident shortwave radiation at the top of the atmosphere. c As in panel a, but for the AABW volume fraction in the global deep ocean

The volume fraction of AABW is also more sensitive to orbital forcing at low CO₂, although it shows a pronounced nonlinearity, peaking at the average late Pleistocene CO₂ of 220 ppm (Fig. 16c). Given that the volume fraction of AABW can play a role in modulating atmospheric CO₂ through disequilibrium carbon storage (Skinner 2009; Eggleston and Galbraith 2018), it seems possible that this enhanced sensitivity of AABW fraction to orbital forcing at Pleistocene temperatures increased the carbon cycle response to orbital forcing relative to earlier in the Cenozoic.

We suggest that increased orbital sensitivity of albedo-driven changes in surface air temperature, combined with increased sensitivity of CO₂ to orbital forcing (particularly obliquity) via variations in the volume fraction of AABW, contributed to the intensification of the ice ages through the Plio–Pleistocene. This increased sensitivity may have been further amplified through additional Earth system feedbacks, such as the production of dust by Patagonian glaciers (Martinez-Garcia et al. 2011). Thus a gradual drawdown of CO₂, resulting from an excess of carbon burial relative to CO₂ supply from outgassing and weathering, would have increased the leverage of orbital forcing on global climate, a key ingredient in the development of the high-amplitude ice age cycles of the Pleistocene.

5 Summary

Our simulations provide a very broad view on the equilibrium response of a comprehensive ocean–atmosphere model to the three most important external forcings of the late Cenozoic: greenhouse gas concentrations, astronomical changes, and changes in ice sheet size. The simulations do not account for changes in vegetation or aerosols, and hold some important parameters (e.g. the orbit eccentricity) fixed. Moreover, the coarse resolution of the model, the absence of important processes (e.g. interactive ice-shelves, downslope flows of dense waters) and the limitations of the parameterizations (e.g. clouds in the atmosphere, tidal mixing and mesoscale eddy effects in the ocean), introduce important biases in the results. Nonetheless, the simulations provide a bird’s-eye view of the expectations of a modern climate model for equilibrium changes under a broad range of boundary conditions. As such, they can be seen as networks of internally-consistent hypotheses that can be tested against observations.

We summarize the main points as follows:

• 1. The globally-averaged surface air temperature is relatively high when obliquity is high and when the precessional phase results in weak boreal seasonality.

• 2. The volume contribution of AABW to the global deep ocean varies nonlinearly with global temperature. AABW maxima occur at both low and high CO₂, with a minimum at interglacial CO₂. The deep ocean volume ventilated from the North Atlantic varies as a mirror image of this, showing a maximum at interglacial CO₂.

• 3. The simulated changes in deep water-mass volumes can be predicted from the density difference between surface waters of the AABW and NADW formation regions. This appears to be consistent with the density-overlap mechanism discussed by Wolfe and Cessi (2009, 2010, 2011).

• 4. Under low CO₂, relatively dense, salty AABW is produced. Cold atmospheric temperatures decrease precipitation in the Southern Ocean and enhance the formation and equatorward export of Antarctic sea ice, turning the high latitude Southern Ocean into a net salinifying environment, instead of the net freshening environment it is today. In a cold world, low obliquity further amplifies sea ice cycling and Southern Ocean salinification, while precession has a smaller direct impact on the ocean.

• 5. Full LGM ice sheets also increase the salinity and volume contribution of AABW, through the downstream impact of cold NADW on Southern Ocean sea ice cycling. However, the large ice sheets simultaneously accelerate the AMOC by increasing convection in the North Atlantic.

• 6. The combination of low CO₂, low obliquity and a low-elevation Laurentide ice sheet leads to a ‘sweet spot’ at which the AMOC is most unstable. We suggest that the AMOC was most susceptible to freshwater perturbations under this state, contributing to enhanced Dansgaard–Oeschger style abrupt climate change when the Southern Ocean was very cold due to low CO₂ and/or low obliquity but ice sheets were not at their maxima.

• 7. At higher CO₂, the influence of sea ice is diminished as its areal extent contracts. Above 405 ppm the AABW volume expands once again, which we attribute to the fact that the Antarctic convective regions warm more slowly than their North Atlantic counterparts.

• 8. The simulated deep ocean temperature and salinity of the ‘glacial’ state agree very closely with reconstructed deep ocean temperature and salinity for the LGM. The simulated expansion of AABW is also qualitatively consistent with common interpretations of δ¹³C and εNd. We take this consistency as an indication that prior disagreements between model LGM simulations and reconstructions could have arisen from model-specific biases or insufficient equilibration, rather than a general inability for modern GCMs to reproduce the glacial climate state—as far as deep water temperature, salinity and surface origin are concerned.

• 9. Despite the close agreement with LGM temperature and salinity, the radiocarbon simulation suggests that the glacial deep ocean ideal age is on the order of 200 years too low, implying that it should have been quite similar to the pre-industrial, despite the strong sea-ice driven buoyancy fluxes in the Southern Ocean. This discrepancy may arise from an inadequate representation of deep ocean mixing as parameterized in the model. It is also possible that the LGM was not in steady state. We highlight the unreliability of simulated ventilation rates in existing climate models, given the poor representation of deep mixing processes.

• 10. The simulated glacial deep ocean O₂ concentrations also appear to be incorrect, with significantly higher concentrations than expected based on proxy evidence. We propose that this biogeochemical error reflects model biases in the formation and consumption processes of dense waters, and/or important changes in the marine ecosystem in addition to the changes in deep water-mass geometry. Unresolved changes in the Southern Ocean near-surface environment (e.g. seasonal stratification and/or iron fertilization) or more widely-distributed changes (e.g. changes in the vertical profile of remineralization or low latitude export) could play important roles.

• 11. Under lower CO₂, the model shows a greater climate sensitivity to orbital forcing, due to an increased response of albedo, and an enhanced sensitivity of deep ocean ventilation via the control of sea ice on AABW salinity. We suggest that an effect of AABW volume on ocean carbon storage, together with the enhanced albedo sensitivity, gave greater leverage to orbital forcing as atmospheric CO₂ declined over the Cenozoic. This tendency would have contributed to the Plio–Pleistocene development and intensification of the ice age cycles.

The analysis provided here has only scratched the surface of the behaviour exhibited in this large suite of model simulations. It may be valuable to undertake further analysis, such as for comparison with specific proxy sites, comparison with other models, exploration of the atmospheric dynamics, or consideration of the changes in ocean circulation according to existing theoretical frameworks (e.g. Gnanadesikan 1999; Wolfe and Cessi 2011; Nikurashin and Vallis 2012; de Lavergne et al. 2017). To this end we will make all simulation output freely available.