# 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

Eric D Galbraith, Casimir de Lavergne

► To cite this version:

Eric D Galbraith, Casimir de Lavergne. 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. Climate Dynamics, 2019, 52 (1-2), pp.653-679. 10.1007/s00382-018-4157-8. hal-02437667

HAL Id: hal-02437667

https://hal.science/hal-02437667

Submitted on 1 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# 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

Eric Galbraith¹,²,³ Ⓒ · Casimir de Lavergne⁴

Received: 29 May 2017 / Accepted: 25 February 2018 / Published online: 9 March 2018

© The Author(s) 2018. This article is an open access publication

## 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₂.

**Keywords** Ocean circulation · Climate · Earth system model · Carbon dioxide · Ice ages · Ocean biogeochemistry

## 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

<figure><img src="image_3.png" /></figure>

<footnote>Eric Galbraith
eric.d.galbraith@gmail.com</footnote>
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 (Togweiler 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

<figure><img src="image_4.png" /></figure>
the relative importance—and interactions—of the vari-
ous mechanisms over the full range of climatic states that
have held sway over the past few million years. Addition-
ally, though simplified models have illustrated how the
interaction of ocean circulation with the marine ecosys-
tem 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. Sim-
mons 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 equilib-
rium 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 CO2
forcings. The model we use is CM2Mc, a coarser-reso-
lution 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
CO2, 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 for-
ing variables. Thus, in contrast with experiments that arbi-
trarily 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 focus-
ing on the circulation of the global deep ocean. We also
compare the results with some key paleoceanographic
reconstructions from the LGM, showing encouraging con-
sistency 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 (Grif-
fies 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 (Sim-
mons 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 Labora-
tory 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 unre-
solved flows between islands. These changes altered the
pre-industrial control state slightly, but the overall simula-
tion 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 atmos-
phere, 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 there-
fore 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-indus-
trial 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 aver-
age 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; Grif-
fies 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 topo-
graphic slopes. As a result, convective chimneys, instead
of downslope currents, are the dominant mode of dense
water injection to depth, biasing the properties and sensi-
tivity of deep water masses to some degree.

Water masses were ‘tagged’ in the surface layer in order
to unambiguously identify volume contributions from dif-
ferent 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 AAPW
and NADW volumes as the global volume integral of the
Southern and North Atlantic tracer concentrations, respec-
tively, 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 applica-
ble 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 AAPW and NADW volumes. Thus, AAPW 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 AAPW 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 Dis-
solved 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 atmos-
pheric Δ¹⁴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 bio-
geochemical 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 sig-
nificantly if forced by other greenhouse gases (CH₄, N₂O)
at equivalent levels of radiative forcing.

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 lati-
tudinal 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 perihe-
lion (Held 1982). The precessional phase was set to the posi-
tions 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 sur-
face 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.

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

<figure><img src="image_5.png" /></figure>
<table><caption>**Table 1** Overview of model simulations</caption><thead><tr><th rowspan="3">CO<sub>2</sub>-ice sheets</th><th rowspan="3">Integration time<br>(prior spinup)</th><th colspan="4">Global average Δ<sup>14</sup>C drift, final century (last 5 centuries)</th></tr><tr><th colspan="2">o22, p270</th><th colspan="2">o24.5, p270</th></tr><tr><th>%c century<sup>-1</sup></th><th>%c century<sup>-1</sup></th><th>%c century<sup>-1</sup></th><th>%c century<sup>-1</sup></th></tr></thead><tbody><tr><td>147-PI</td><td>3100</td><td>-1.3 (-1.5)</td><td>-1.3 (-1.4)</td><td>-1.3 (-1.4)</td><td>-1 (-1.1)</td></tr><tr><td>180-PI</td><td>5000</td><td>-0.6 (-0.6)</td><td>-0.6 (-0.7)</td><td>-0.6 (-0.6)</td><td>-0.3 (-0.4)</td></tr><tr><td>180-G/PI</td><td>4600</td><td>-0.6 (-0.7)</td><td>-0.8 (-0.8)</td><td>-0.5 (-0.6)</td><td>-0.4 (-0.4)</td></tr><tr><td>180-G</td><td>6000 (3200)</td><td>-0.5 (-0.5)</td><td><strong>-0.5 (-0.6) g</strong></td><td>-0.5 (-0.6)</td><td>-0.4 (-0.4)</td></tr><tr><td>220-PI</td><td>5000</td><td>-0.7 (-0.8)</td><td>-0.6 (-0.7)</td><td>-0.4 (-0.5)</td><td>-0.3 (-0.3)</td></tr><tr><td>220-G</td><td>6000</td><td>-0.7 (-0.8)</td><td>-0.6 (-0.7)</td><td>-0.5 (-0.6)</td><td>-0.3 (-0.4)</td></tr><tr><td>270-PI</td><td>3200</td><td>-0.8 (-0.8)</td><td>-0.6 (-0.5)</td><td>-0.6 (-0.6)</td><td><strong>-0.4 (-0.5) ig</strong></td></tr><tr><td>405-PI</td><td>3200</td><td>-0.6 (-0.6)</td><td>-0.4 (-0.5)</td><td>-0.2 (-0.3)</td><td>-0.2 (-0.3)</td></tr><tr><td>607-PI</td><td>3200 (2000)</td><td>-0.3 (-0.3)</td><td>-0.5 (-0.5)</td><td>-0.3 (-0.4)</td><td>-0.5 (-0.6)</td></tr><tr><td>911-PI</td><td>3400 (2000)</td><td>-0.3 (-0.4)</td><td>-0.4 (-0.5)</td><td>-0.2 (-0.3)</td><td>-0.3 (-0.5)</td></tr></tbody></table>

The first column lists the different configurations of CO₂ (in ppm) and land-based ice sheets (either 'PI' for Pre-Industrial or 'G' for LGM). The second column lists the length of the simulation; in the case of an initial spinup with only one orbital configuration, the length of the initial portion is shown in parentheses. The final four columns give the changes of the whole ocean radiocarbon content over the last 100 and 500 years (the latter in parentheses) as an indication of model drift. Each column corresponds to one orbital configuration, defined by the obliquity (o22 or o24.5) and precessional phase (p270 or p90). The generic 'glacial' and 'interglacial' states are indicated by boldface

<figure><img src="image_6.png" /><figcaption>**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)</figcaption></figure>

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 $\Delta^{14}\text{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 $\text{CO}_2$ 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

<figure><img src="image_7.png" /><figcaption>**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</figcaption></figure>
of change driven by the two types of forcing are very differ-
ent, a factor which is not evident in the common simplifica-
tion 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 rota-
tional axis), intensifying the equator-to-pole gradient by
more than 3 °C. Precession produces a more complex spa-
tial 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.

Figure 3 shows the globally-averaged surface air tempera-
ture 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 inso-
lation, 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).

<figure><img src="image_8.png" /><figcaption>Fig. 3 Globally-averaged surface air temperature vs. CO₂. Each tri-
angle 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</figcaption></figure>

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 cool-
ing 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 mod-
els generally show a larger LGM-PI change than the corre-
sponding 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 con-
tribution of AABW to the deep (> 1 km) ocean (Fig. 4a),
varying from 30 to 80%. The AABW volume shows a
strong dependence on CO2, with a striking nonlinearity:
the AABW volume reaches a minimum at intermediate
CO2 (270–405 ppm) and increases toward both higher and
lower CO2. The effect of orbital forcing is also significant,
and displays a further nonlinearity with CO2. For CO2 of
270 ppm and below, low obliquity (blue) increases the Ant-
arctic ventilation volume, whereas only precession plays a
role at higher CO2, producing a greater fraction of AABW
when boreal seasonality is weak (downward-pointing trian-
gles). 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 CO2, low obliq-
uity 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
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**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</figcaption></figure>

Pacific, and d the maximum AMOC transport at 40°N. All symbols are as defined in Fig. 3

decreases both with increasing CO₂ and with larger ice sheets (Fig. 4c).

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

<figure><img src="image_13.png" /></figure>
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).

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₂

<figure><img src="image_14.png" /><figcaption>**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</figcaption></figure>
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.

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.

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

<figure><img src="image_15.png" /><figcaption>**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</figcaption></figure>

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

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

<figure><img src="image_16.png" /></figure>
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>**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⁻¹) of the southern westerlies does not show a significant relationship with the AABW fraction</figcaption></figure>

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

<figure><img src="image_19.png" /><figcaption>**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</figcaption></figure>

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.

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,
<figure><img src="image_20.png" /><figcaption>(c) Southern Ocean sea ice transport</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(d) Sea surface freshwater flux south of 60S</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>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</figcaption></figure>

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

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₂, 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).

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.

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

<figure><img src="image_24.png" /><figcaption>**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</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**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</figcaption></figure>

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

<figure><img src="image_26.png" /></figure>
<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

**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

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
<figure><img src="image_31.png" /><figcaption>**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</figcaption></figure>

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 \pm 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 \pm 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.

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 mineralized 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).

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

<figure><img src="image_32.png" /></figure>
<figure><img src="image_33.png" /><figcaption>**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)</figcaption></figure>

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.

### 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
<figure><img src="image_34.png" /></figure>

**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 ($\sigma_1$) 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 $\sigma_1$ 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

(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 abysmal 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 $^{14}$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 $\Delta^{14}$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_2$ 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 $^{14}$C age caused by a reduction of the pre-formed $^{14}$C activity of AABW, its volumetric dominance of the glacial ocean, and the $pCO_2$ 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_2$ 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_3$ 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_3$ 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_2$ (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

<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /><figcaption>**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</figcaption></figure>

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.

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

<figure><img src="image_38.png" /></figure>
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.

**Acknowledgements** We are grateful to Compute Canada and Scinet (University of Toronto) for providing the computational resources to carry out the simulations, and to the Canadian Foundation for Innovation for local computing infrastructure. We thank Rick Slater, Simon Yang and Nicolas Brown for assistance with configuring and running the model, and Tim Merlis, Luke Skinner, Sam Jaccard, Sarah Eggleston, Robbie Toggweiler, Daniel Sigman, Robert Anderson, Jorge Sarmiento and Jaime Palter for discussions. Pearse Buchanan and three anonymous reviewers provided constructive comments on the manuscript. EDG acknowledges financial support from the Spanish Ministry of Economy and Competitiveness, through the María de Maeztu Programme for Centres/Units of Excellence in R&D (MDM-2015-0552). CdL acknowledges Australian Research Council support through grant FL150100090. All model runscripts, code and simulation output are freely available from the authors, or can be obtained by download from https://earthsystemdynamics.org/cm2mc-simulation-library/.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
## References

Adkins JF (2013) The role of deep ocean circulation in setting glacial climates. Paleooceanography. https://doi.org/10.1002/palo.20046

Adkins JF, Schrag DP (2003) Reconstructing Last Glacial Maximum bottom water salinities from deep-sea sediment pore fluid profiles. Earth Planet Sci Lett 216:109–123

Adkins JF, McIntyre K, Schrag DP (2002) The salinity, temperature, and delta O-18 of the glacial deep ocean. Science 298:1769–1773

Anderson RF, Ali S, Bradtmiller LI, Nielsen SHH, Fleisher MQ, Anderson BE, Burckle LH (2009) Wind-driven upwelling in the southern ocean and the deglacial rise in atmospheric CO₂. Science 323:1443–1448. https://doi.org/10.1126/Science.1167441

Annan J, Hargreaves J (2013) A new global reconstruction of temperature changes at the Last Glacial Maximum. Clim Past 9:367–376

Bard E (1988) Correction of accelerator mass spectrometry ¹⁴C ages measured in planktonic foraminifera: Paleooceanographic implications. Paleooceanography 3:635–695

Bell DB, Jung SJA, Kroon D (2015) The Plio-Pleistocene development of Atlantic deep-water circulation and its influence on climate trends. Quatern Sci Rev 123:265–282. https://doi.org/10.1016/j.quascirev.2015.06.026

Bereiter B, Shackleton S, Baggenstos D, Kawamura K, Severinghaus J (2018) Mean global ocean temperatures during the last glacial transition. Nature 553:39–44

Böhm E et al (2015) Strong and deep Atlantic meridional overturning circulation during the last glacial cycle. Nature 517:73–76

Bopp L, Resplandy L, Untersee A, Le Mezo P, Kageyama M (2017) Ocean (de)oxygenation from the Last Glacial Maximum to the twenty-first century: insights from Earth System models. Philos Trans R Soc A Math Phys Eng Sci. https://doi.org/10.1098/rsta.2016.0323

Bouttes N, Paillard D, Roche D (2010) Impact of brine-induced stratification on the glacial carbon cycle. Clim Past 6:575–589

Boyle EA, Keigwin L (1982) Deep circulation of the North Atlantic over the Last 200,000 years: geochemical evidence. Science 218:784–787

Braconnot P et al (2007) Results of PMIP2 coupled simulations of the Mid-Holocene and Last Glacial Maximum—Part 1: experiments and large-scale features. Clim Past 3:261–277

Bradtmiller LI, Anderson RF, Sachs JP, Fleisher MQ (2010) A deeper respired carbon pool in the glacial equatorial Pacific Ocean. Earth Planet Sci Lett 299:417–425. https://doi.org/10.1016/j.epsl.2010.09.022

Broccoli A, Manabe S (1987) The influence of continental ice, atmospheric CO₂, and land albedo on the climate of the last glacial maximum. Clim Dyn 1:87–99

Brovkin V, Ganopolski A, Archer D, Munhoven G (2012) Glacial CO₂ cycle as a succession of key physical and biogeochemical processes. Clim Past 8:251–264. https://doi.org/10.5194/Cp-8-251-2012

Brown N, Galbraith E (2015) Hosed vs. unhosed: global response to interruptions of the Atlantic Meridional Overturning, with and without freshwater forcing. Clim Past Discuss 11:4669–4700

Buchanan PJ, Matear RJ, Lenton A, Phipps SJ, Chase Z, Etheridge DM (2016) The simulated climate of the Last Glacial Maximum and insights into the global marine carbon cycle. Clim Past 12:2271

Buizert C, Schmittner A (2015) Southern Ocean control of glacial AMOC stability and Dansgaard–Oeschger interstadial duration. Paleooceanography 30:1595–1612. https://doi.org/10.1002/2015P A002795

Burke A, Robinson LF (2012) The Southern Ocean's role in carbon exchange during the last deglaciation. Science 335:557–561

Butzin M, Prange M, Lohmann G (2005) Radiocarbon simulations for the glacial ocean: the effects of wind stress, Southern Ocean sea ice and Heinrich events. Earth Planet Sci Lett 235:45–61

Cessi P, Jones CS (2017) Warm-route versus cold-route interbasin exchange in the meridional overturning circulation. J Phys Oceanogr 47:1981–1997

Cheng H et al (2009) Ice age terminations. Science 326:248–252. https://doi.org/10.1126/Science.1177840

Crowley TJ (1992) North Atlantic deep water cools the southern hemisphere. Paleooceanography 7:489–497

Curry WB, Oppo D (2005) Glacial water mass geometry and the distribution of δ¹³C of sigma-CO₂ in the western. Atlantic Ocean Paleooceanography. https://doi.org/10.1029/2004PA001021

De Boer AM, Sigman DM, Toggweiler JR, Russell JL (2007) The effect of global ocean temperature change on deep ocean ventilation. Paleooceanography 22

de Lavergne C, Palter JB, Galbraith ED, Bernardello R, Marinov I (2014) Cessation of deep convection in the open Southern Ocean under anthropogenic climate change. Nat Clim Change 4:278–282. https://doi.org/10.1038/nclimate2132

de Lavergne C, Madec G, Roquet F, Holmes RM, McDougall TJ (2017) Abyssal ocean overturning shaped by seafloor distribution. Nature 551:181–186

DeVries T, Primeau F (2010) An improved method for estimating water-mass ventilation age from radiocarbon data. Earth Planet Sci Lett 295:367–378

Dunne JP et al (2012) GFDL's ESM2 global coupled climate—carbon earth system models. Part I: Physical formulation and baseline simulation characteristics. J Clim 25:6646–6665. https://doi.org/10.1175/jcli-d-11-00560.1

Duplessy JC, Shackleton NJ, Matthews RK, Prell WL, Ruddiman WF, Caralp MH, Hendy CH (1984) C-13 record of benthic foraminifera in the last interglacial ocean—implications for the carbon cycle and the global deep-water circulation. Quatern Res 21:225–243

Durack PJ, Wijffels SE, Matear RJ (2012) Ocean salinities reveal strong global water cycle intensification during 1950–2000. Science 336:455–458

Eggleston S, Galbraith ED (2018) The devil's in the disequilibrium: sensitivity of ocean carbon storage to climate state and iron fertilization in a general circulation model. Biogeosci Discuss. https://doi.org/10.5194/bg-2017-328

Emiliani C (1955) Pleistocene temperatures. J Geol 63:537–578

Ferrari R, Jansen MF, Adkins JF, Burke A, Stewart AL, Thompson AF (2014) Antarctic sea ice control on ocean circulation in present and glacial climates. Proc Natl Acad Sci 201323922

Fichefet T, Hovine S, Duplessy JC (1994) A model study of the Atlantic thermohaline circulation during the last glacial maximum. Nature 372:252–255

Fischer H et al (2010) The role of Southern Ocean processes in orbital and millennial CO₂ variations—a synthesis. Quatern Sci Rev 29:193–205

Fox-Kemper B, Ferrari R, Hallberg R (2008) Parameterization of mixed layer eddies. Part I: Theory and diagnosis. J Phys Oceanogr 38:1145–1165

François R et al (1997) Contribution of Southern Ocean surface-water stratification to low atmospheric CO₂ concentrations during the last glacial period. Nature 389:929–935

Frölicher TL, Winton M, Sarmiento JL (2014) Continued global warming after CO₂ emissions stoppage. Nat Clim Change 4:40–44

Galbraith ED, Jaccard SL (2015) Deglacial weakening of the oceanic soft tissue pump: global constraints from sedimentary nitrogen isotopes and oxygenation proxies. Quatern Sci Rev 109:38–48

Galbraith ED et al (2007) Carbon dioxide release from the North Pacific abyss during the last deglaciation. Nature 449:890–893. https://doi.org/10.1038/Nature06227
Galbraith ED, Gnanadesikan A, Dunne JP, Hiscock MR (2010) Regional impacts of iron-light colimitation in a global biogeochemical model. *Biogeosciences* 7:1043–1064

Galbraith ED et al (2011) Climate variability and radiocarbon in the CM2Mc earth system model. *J Clim* 24:4230–4254. https://doi.org/10.1175/2011jcli3919.1

Galbraith ED, Kwon EY, Bianchi D, Hain MP, Sarmiento JL (2015) The impact of atmospheric pCO₂ on carbon isotope ratios of the atmosphere and ocean. *Global Biogeochem Cycles* 29:307–324

Galbraith ED, Merlis TM, Palter JB (2016) Destabilization of glacial climate by the radiative impact of Atlantic Meridional Overturning Circulation disruptions. *Geophys Res Lett* 43:8214–8221

Ganopolski A, Winkelmann R, Schellnhuber H (2016) Critical insolation—CO₂ relation for diagnosing past and future glacial inception. *Nature* 529:200–203

Gebbie G (2014) How much did Glacial North Atlantic Water shoal? *Paleoceanography* 29:190–209

Gent PR, McWilliams JC (1990) Isopycnal mixing in ocean circulation models. *J Phys Oceanogr* 20:150–155

Gildor H, Tziperman E, Toggweiler JR (2002) Sea ice switch mechanism and glacial-interglacial CO₂ variations. *Global Biogeochem Cycles* 16:art.no.-1032

Gnanadesikan A (1999) A simple predictive model for the structure of the oceanic Pycnocline. *Science* 283:2077–2079

Gong X, Zhang X, Lohmann G, Wei W, Zhang X, Pfeiffer M (2015) Higher Laurentide and Greenland ice sheets strengthen the North Atlantic ocean circulation. *Clim Dyn* 45:139–150

Goosse H, Fichefet T (1999) Importance of ice-ocean interactions for the global ocean circulation: a model study. *J Geophys Res Oceans* 104:23337–23355

Griffies SM (1998) The Gent-McWilliams skew flux. *J Phys Oceanogr* 28:831–841

Hain MP, Sigman DM, Haug GH (2010) Carbon dioxide effects of Antarctic stratification, North Atlantic Intermediate Water formation, and subantarctic nutrient drawdown during the last ice age: diagnosis and synthesis in a geochemical box model. *Global Biogeochem Cycles* 24 doi:https://doi.org/10.1029/2010gb003790

Hall IR, McCave IN, Shackleton NJ, Weedon GP, Harris SE (2001) Intensified deep Pacific inflow and ventilation in Pleistocene glacial times. *Nature* 412:809–812

Hayes CT et al (2014) A stagnation event in the deep South Atlantic during the last interglacial period. *Science* 346:1514–1517

Held I (1982) Climate models and the astronomical theory of the ice Ages. *Icarus* 50:449–461

Held IM, Soden BJ (2006) Robust responses of the hydrological cycle to global warming. *J Clim* 19:5686–5699

Heuzé C, Heywood KJ, Stevens DP, Ridley JK (2013) Southern Ocean bottom water characteristics in CMIP5 models. *Geo-phys Res Lett* 40:1409–1414

Hewitt CD, Stouffer RJ, Broccoli AJ, Mitchell JFB, Valdes PJ (2003) The effect of ocean dynamics in a coupled GCM simulation of the Last Glacial Maximum. *Clim Dyn* 20:203–218

Ho SL, Laepple T (2015) Glacial cooling as inferred from marine temperature proxies TEX86 and UK' 37. *Earth Planet Sci Lett* 409:15–22

Hodell DA, Venz-Curtis KA (2006) Late Neogene history of deepwater ventilation in the Southern Ocean. *Geochem Geophys Geosyst* 7

Hoogakker BA, Chapman MR, McCave IN, Hillaire-Marcel C, Ellison CR, Hall IR, Telford RJ (2011) Dynamics of North Atlantic deep water masses during the Holocene. *Paleoceanography* 26

Hoogakker BA, Elderfield H, Schmiedl G, McCave IN, Rickaby RE (2015) Glacial-interglacial changes in bottom-water oxygen content on the Portuguese margin. *Nat Geosci* 8:40–43

Howe JN, Piotrowski AM, Noble TL, Mulitza S, Chiessi CM, Bayon G (2016) North Atlantic deep water production during the last glacial maximum. *Nat Commun* 7

Huybers P, Wunsch C (2005) Obliquity pacing of the late Pleistocene glacial terminations. *Nature* 434:491–494

Insua TL, Spivack AJ, Graham D, D'Hondt S, Moran K (2014) Reconstruction of Pacific Ocean bottom water salinity during the Last Glacial Maximum. *Geophys Res Lett* 41:2914–2920

Jaccard SL, Hayes CT, Martinez-Garcia A, Hodell DA, Anderson RF, Sigman DM, Haug GH (2013) Two modes of change in southern ocean productivity over the past million years. *Science* 339:1419–1423. https://doi.org/10.1126/Science.1227545

Jaccard SL, Galbraith ED, Martínez-García A, Anderson RF (2016) Covariation of deep Southern Ocean oxygenation and atmospheric CO₂ through the last ice age. *Nature* 530:207–210

Jansen MF (2017) Glacial ocean circulation and stratification explained by reduced atmospheric temperature. *Proc Natl Acad Sci* 114:45–50

Jansen MF, Nadeau L-P (2016) The effect of southern ocean surface Buoyancy loss on the deep-ocean circulation and stratification. *J Phys Oceanogr* 46:3455–3470

Kawamura K et al (2017) State dependence of climatic instability over the past 720,000 years from Antarctic ice cores and climate modeling. *Sci Adv* 3:e1600446

Keeling RF, Stephens BB (2001) Antarctic sea ice and the control of Pleistocene climate instability. *Paleoceanography* 16:112–131

Keigwin LD, Swift SA (2017) Carbon isotope evidence for a northern source of deep water in the glacial western North Atlantic. *Proc Natl Acad Sci* 114:2831–2835

Kienast SS, Kienast M, Mix AC, Calvert SE, François R (2007) Thoria-230 normalized particle flux and sediment focusing in the Panama Basin region during the last 30,000 years. *Paleoceanography*. https://doi.org/10.1029/2006PA001357

Klockmann M, Mikolajewicz U, Marotzke J (2016) The effect of greenhouse gas concentrations and ice sheets on the glacial AMOC in a coupled climate model. *Clim Past* 12:1829–1846

Köhler P, Fischer H, Munhoven G, Zeebe RE (2005) Quantitative interpretation of atmospheric carbon records over the last glacial termination. *Global Biogeochem Cycles* 19 https://doi.org/10.1029/2004gb002345

Korff L, Dobeneck T, Frederichs T, Kasten S, Kuhn G, Gersonde R, Diekmann B (2016) Cyclic magnetite dissolution in Pleistocene sediments of the abyssal northwest Pacific Ocean: evidence for glacial oxygen depletion and carbon trapping. *Paleoceanography*

Krueger S, Leuschner D, Ehrmann W, Schmiedl G, Mackensen A (2012) North Atlantic Deep Water and Antarctic Bottom Water variability during the last 200 ka recorded in an abyssal sediment core off South Africa. *Global Planet Change* 80:180–189

Kutzbach J, Liu X, Liu Z, Chen G (2008) Simulation of the evolutionary response of global summer monsoons to orbital forcing over the past 280,000 years. *Clim Dyn* 30:567–579

Kutzbach JE, He F, Vavrus SJ, Ruddiman WF (2013) The dependence of equilibrium climate sensitivity on climate state: Applications to studies of climates colder than present. *Geophys Res Lett* 40:3721–3726

Kwon EY, Hain MP, Sigman DM, Galbraith ED, Sarmiento JL, Togweiler JR (2012) North Atlantic ventilation of “southern-sourced” deep water in the glacial ocean. *Paleoceanography* 27

Large WG, McWilliams JC, Doney SC (1994) Oceanic vertical mixing—a review and a model with a nonlocal boundary-layer parameterization. *Rev Geophys* 32:363–403

Laskar J, Robutel P, Joutel F, Gastineau M, Correia A, Levrard B (2004) A long-term numerical solution for the insolation quantities of the Earth. *Astron Astrophys* 428:261–285

Lee JE, Shen A, Fox-Kemper B, Ming Y (2017) Hemispheric sea ice distribution sets the glacial tempo. *Geophys Res Lett* 44:1008–1014

Lefebvre V, Donnadieu Y, Sepulchre P, Swingedouw D, Zhang ZS (2012) Deciphering the role of southern gateways and carbon
dioxide on the onset of the Antarctic Circumpolar Current. *Paleoceanography* 27

Li C, Battisti DS (2008) Reduced Atlantic storminess during Last Glacial Maximum: evidence from a coupled climate model. *J Clim* 21:3561–3579

Lippold J et al (2016) Deep water provenance and dynamics of the (de) glacial Atlantic meridional overturning circulation. *Earth Planet Sci Lett* 445:68–78

Lisiecki LE, Raymo ME (2007) Plio-Pleistocene climate evolution: trends and transitions in glacial cycle dynamics. *Quatern Sci Rev* 26:56–69. https://doi.org/10.1016/j.quascirev.2006.09.005

Lu Z et al. (2016) Oxygen depletion recorded in upper waters of the glacial Southern Ocean. *Nat Commun* 7

MacKinnon JA et al (2017) Climate process team on internal-wave driven ocean mixing. *Bull Am Meteorol Soc*

Manabe S, Bryan K (1985) CO₂-induced change in a coupled ocean–atmosphere model and its paleoclimatic implications. *J Geophys Res Oceans* 90:11689–11707

Marcott SA, Shakun JD, Clark PU, Mix AC (2013) A reconstruction of regional and global temperature for the past 11,300 years. *Science* 339:1198–1201

Martinez-Garcia A, Rosell-Mele A, Jaccard S, Geibert W, Sigman DM, Haug GH (2011) Southern Ocean dust-climate coupling over the past four million years. *Nature* 476

Martínez-García A et al (2014) Iron fertilization of the Subantarctic Ocean during the last ice age. *Science* 343:1347–1350

Martinson DG (1990) Evolution of the Southern Ocean winter mixed layer and sea ice: open-ocean deep water formation and ventilation. *J Geophys Res* 95:114614–11654

Marzocchi A, Jansen MF (2017) Connecting Antarctic sea ice to deep-ocean circulation in modern and glacial climate simulations. *Geophys Res Lett* 44:6286–6295

Maslin MA, Li XS, Loutre M-F, A. B (1998) The contribution of orbital forcing to the progressive intensification of Northern Hemisphere Glaciation. *Quatern Sci Rev* 17:411–426

Matsumoto K, Hashioka T, Yamanaka Y (2007) Effect of temperature-dependent organic carbon decay on atmospheric *pCO₂*. *J Geophys Res Biogeosci* 112

McManus JF, Oppo DW, Cullen JL (1999) A 0.5-million-year record of millennial-scale climate variability in the North Atlantic. *Science* 283:971–975

McManus JF, Francois R, Gherardi JM, Keigwin LD, Brown-Leger S (2004) Collapse and rapid resumption of Atlantic meridional circulation linked to deglacial climate changes. *Nature* 428:834–837

Menviel L, Yu J, Joos F, Mouchet A, Meissner K, England M (2017) Poorly ventilated deep ocean at the Last Glacial Maximum inferred from carbon isotopes: a data-model comparison study. *Paleoceanography* 32:2–17

Mix AC, Bard E, Schneider R (2001) Environmental processes of the ice age: land, oceans, glaciers (EPILOG). *Quatern Sci Rev* 20:627–657

Muglia J, Schmittner A (2015) Glacial Atlantic overturning increased by wind stress in climate models. *Geophys Res Lett* 42:9862–9868

Munday DR, Johnson HL, Marshall DP (2013) Eddy saturation of equilibrated circumpolar currents. *J Phys Oceanogr* 43:507–532

Nikurashin M, Vallis G (2012) A theory of the interhemispheric meridional overturning circulation and associated stratification. *J Phys Oceanogr* 42:1652–1667

Nilsson J, Langen PL, Ferreira D, Marshall J (2013) Ocean basin geometry and the salinification of the Atlantic Ocean. *J Clim* 26:6163–6184

Orsi AH, Johnson GC, Bullister JL (1999) Circulation, mixing, and production of Antarctic Bottom Water. *Prog Oceanogr* 43:55–109

Otto-Bliesner B et al. (2007) Last Glacial Maximum ocean thermohaline circulation: PMIP2 model intercomparisons and data constraints. *Geophys Res Lett* 34

Paillard D, Parrenin F (2004) The Antarctic ice sheet and the triggering of deglaciations. *Earth Planet Sci Lett* 227:263–271

Pausata FS, Li C, Wettstein J, Kageyama M, Nisancioglu KH (2011) The key role of topography in altering North Atlantic atmospheric circulation during the last glacial period. *Clim Past* 7:1089–1101

Pena LD, Goldstein SL (2014) Thermohaline circulation crisis and impacts during the mid-Pleistocene transition. *Science* 345:318–322

Piotrowski AM, Goldstein SL, Hemming SR, Fairbanks RG (2005) Temporal relationships of carbon cycling and ocean circulation at glacial boundaries. *Science* 307:1933–1938

Rahmstorf S (2002) Ocean circulation and climate during the past 120,000 years. *Nature* 419:207–214

Raymo M, Lisiecki L, Nisancioglu KH (2006) Plio-Pleistocene ice volume, Antarctic climate, and the global $\delta^{18}$O record. *Science* 313:492–495

Ritz SP, Stocker TF, Grimalt JO, Menviel L, Timmermann A (2013) Estimated strength of the Atlantic overturning circulation during the last deglaciation. *Nat Geosci* 6:208–212. https://doi.org/10.1038/Ngeo1723

Roe GH, Lindzen RS (2001) The mutual interaction between continental-scale ice sheets and atmospheric stationary waves. *J Clim* 14:1450–1465

Russell JL, Stouffer RJ, Dixon KW (2006) Intercomparison of the Southern Ocean circulations in IPCC coupled model control simulations. *J Clim* 19:4560–4575

Sarmiento JL, Toggweiler JR (1984) A new model for the role of the oceans in determining atmospheric *pCO₂*. *Nature* 308:621–624

Sarnthein M, Winn K, Jung SJ, Duplessy JC, Labeyrie L, Erlenkeuser H, Ganssen G (1994) Changes in east Atlantic deepwater circulation over the last 30,000 years: eight time slice reconstructions. *Paleoceanography* 9:209–267

Sarnthein M, Schneider B, Grootes PM (2013) Peak glacial $^{14}$C ventilation ages suggest major draw-down of carbon into the abyssal ocean. *Clim Past*

Schmittner A (2003) Southern Ocean sea ice and radiocarbon ages of glacial bottom waters. *Earth Planet Sci Lett* 213:53–62

Schmittner A, Clement AC (2002) Sensitivity of the thermohaline circulation to tropical and high latitude freshwater forcing during the last glacial-interglacial cycle. *Paleoceanography*

Schmittner A, Green J, Wilmes SB (2015) Glacial ocean overturning intensified by tidal mixing in a global circulation model. *Geophys Res Lett* 42:4014–4022

Schrag D, De Paolo D (1993) Determination of $\delta^{18}$O of seawater in the deep ocean during the Last Glacial Maximum. *Paleoceanography* 8:1–6

Sherriff-Tadano S, Abe-Ouchi A, Yoshimori M, Oka A, Chan W-L (2017) Influence of glacial ice sheets on the Atlantic meridional overturning circulation through surface wind change. *Clim Dyn* 1–23

Shin SI, Liu Z, Otto-Bliesner BL, Kutzbach JE, Vavrus SJ (2003) Southern Ocean sea-ice control of the glacial North Atlantic thermohaline circulation. *Geophys Res Lett* 30

Sigman DM, Boyle EA (2000) Glacial/interglacial variations in atmospheric carbon dioxide. *Nature* 407:859–869

Sigman DM, Jaccard SL, Haug GH (2004) Polar ocean stratification in a cold climate. *Nature* 428:59–63

Sigman DM, De Boer AM, Haug GH (2007) Antarctic stratification, atmospheric water vapor, and Heinrich events: a hypothesis for late Pleistocene deglaciations. *Mechanisms and impacts-past and future changes of meridional overturning, ocean circulation*, pp 335–349
Sigman DM, Hain MP, Haug GH (2010) The polar ocean and glacial cycles in atmospheric CO₂ concentration. Nature 466:47–55. https://doi.org/10.1038/Nature09149

Simmons HL, Jayne SR, St Laurent LC, Weaver AJ (2004) Tidally driven mixing in a numerical model of the ocean general circulation. Ocean Model 6:245–263

Simmons C, Matthews H, Mysak L (2016) Deglacial climate, carbon cycle and ocean chemistry changes in response to a terrestrial carbon release. Clim Dyn 46:1287–1299

Skinner L (2009) Glacial-interglacial atmospheric CO₂ change: a possible standing volume effect on deep-ocean carbon sequestration. Clim Past 5:537–550

Skinner LC, Fallon S, Waelbroeck C, Michel E, Barker S (2010) Ventilation of the Deep Southern Ocean and Deglacial CO₂ Rise. Science 328:1147–1151. https://doi.org/10.1126/Science.1183627

Skinner L, McCave I, Carter L, Fallon S, Scrivner A, Primeau F (2015) Reduced ventilation and enhanced magnitude of the deep Pacific carbon pool during the last glacial period. Earth Planet Sci Lett 411:45–52

Stephens BB, Keeling RF (2000) The influence of Antarctic sea ice on glacial-interglacial CO₂ variations. Nature 404:171–174

Stocker TF, Wright DG, Broecker WS (1992) The influence of high-latitude surface forcing on the global thermohaline circulation. Paleooceanography 7:529–541

Stommel H (1961) Thermohaline convection with two stable regimes of flow. Tellus 13:224–230

Stouffer RJ, Manabe S (2003) Equilibrium response of thermohaline circulation to large changes in atmospheric CO₂ concentration. Clim Dyn 20:759–773

Swingedouw D, Braconnot P, Delécluse P, Guilyardi E, Marti O (2007) Quantifying the AMOC feedbacks during a 2 × CO₂ stabilization experiment with land-ice melting. Clim Dyn 29:521–534

Telford R, Li C, Kucera M (2013) Mismatch between the depth habitat of planktonic foraminifera and the calibration depth of SST transfer functions may bias reconstructions. Clim Past 9:859–870

Thornalley DJ, Elderfield H, McCave IN (2009) Holocene oscillations in temperature and salinity of the surface subpolar North Atlantic. Nature 457:711

Toggweiler JR (1999) Variation of atmospheric CO₂ by ventilation of the ocean’s deepest water. Paleooceanography 14:571–588

Toggweiler JR, Russell JL, Carson SR (2006) Midlatitude westerlies, atmospheric CO₂, and climate change during the ice ages. Paleooceanography 21

Wang XT et al (2017) Deep-sea coral evidence for lower Southern Ocean surface nitrate concentrations during the last ice age. Proc Natl Acad Sci 201615718

Watson AJ, Vallis GK, Nikurashin M (2015) Southern Ocean buoyancy forcing of ocean ventilation and glacial atmospheric CO₂. Nat Geosci

Winton M (1997) The effect of cold climate upon North Atlantic deep water formation in a simple ocean–atmosphere model. J Clim:37–51

Wolfe CL, Cessi P (2009) Overturning circulation in an eddy-resolving model: the effect of the pole-to-pole temperature gradient. J Phys Oceanogr 39:125–142

Wolfe CL, Cessi P (2010) What sets the strength of the middepth stratification and overturning circulation in eddying ocean models? J Phys Oceanogr 40:1520–1538

Wolfe CL, Cessi P (2011) The adiabatic pole-to-pole overturning circulation. J Phys Oceanogr 41:1795–1810

Wunsch C (2003) Determining paleoceanographic circulations, with emphasis on the Last Glacial Maximum. Quatern Sci Rev 22:371–385

Wunsch C (2016) Pore fluids and the LGM ocean salinity—Reconsidered. Quatern Sci Rev 135:154–170

Yamamoto A, Abe-Ouchi A, Shigemitsu M, Oka A, Takahashi K, Ohgaito R, Yamanaka Y (2015) Global deep ocean oxygenation by enhanced ventilation in the Southern Ocean under long-term global warming. Global Biogeochem Cycles 29:1801–1815

Zhang X, Lohmann G, Knorr G, Purcell C (2014) Abrupt glacial climate shifts controlled by ice sheet changes. Nature 512:290–294

Zhu J, Liu Z, Zhang X, Eisenman I, Liu W (2014) Linear weakening of the AMOC in response to receding glacial ice sheets in CCSM3. Geophys Res Lett 41:6252–6258. https://doi.org/10.1002/2014GL060891

## Affiliations

Eric Galbraith¹,²,³ • Casimir de Lavergne⁴

¹ ICREA, Pg. Lluís Companys 23, 08010 Barcelona, Spain

² Department of Mathematics, Institut de Ciència i Tecnologia Ambientals (ICTA), Universitat Autonoma de Barcelona, 08193 Barcelona, Spain

³ Department of Earth and Planetary Sciences, McGill University, Montreal, Canada

⁴ School of Mathematics and Statistics, University of New South Wales, Sydney, NSW 2052, Australia