Climate Dynamics (2019) 52:653–679 https://doi.org/10.1007/s00382-018-4157-8 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 Galbraith1,2,3 · Casimir de Lavergne4 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 C ­ O2 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 ­CO2, 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 differ- ence and the AABW volume occur near interglacial ­CO2 (270–400 ppm). At low ­CO2, 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 C ­ O2, 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, enhanc- ing Antarctic sea ice export and leading to very salty and expanded AABW. Together, these results suggest that a ‘sweet spot’ of low ­CO2, 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 pos- sibility 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 ­CO2. 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 ­CO2. Keywords Ocean circulation · Climate · Earth system model · Carbon dioxide · Ice ages · Ocean biogeochemistry 1 Introduction Since the first quantitative reconstructions of ocean tem- * Eric Galbraith peratures were made (Emiliani 1955), investigators have eric.d.galbraith@gmail.com pondered how changes in climate state may have altered Extended author information available on the last page of the article deep ocean circulation during the ice age cycles of the past 13 Vol.:(0123456789) 654 E. Galbraith, C. de Lavergne few million years (Adkins 2013). The greatest focus has indication that the penetration depth of North Atlantic Deep traditionally been placed on understanding the most recent Water (NADW) was reduced during the LGM (Lippold et al. extreme state, the Last Glacial Maximum (LGM), typically 2016). This coherent dual-tracer signal of AABW expansion defined as the period from 19,000 to 23,000 years ago (Mix and contraction over glacial cycles extends back to the mid- et al. 2001). Numerous mechanisms have been championed dle Pleistocene, and appears to have developed alongside the as key drivers behind glacial-interglacial changes in ocean intensification of ice ages about 1 million years ago (Pena state, including changes in southern westerly winds (Tog- and Goldstein 2014). But despite these apparent changes gweiler et al. 2006; Anderson et al. 2009), sea ice dynamics in water-mass geometry, the ratio of another pair of iso- (Stephens and Keeling 2000; Schmittner 2003; Shin et al. topes in sediments, 230Pa/231Th, suggests that the volume 2003; Bouttes et al. 2010; Ferrari et al. 2014), ice sheet transport of the Atlantic Meridional Overturning Circulation topography (Pausata et al. 2011; Gong et al. 2015), tidal (AMOC) during the LGM was not necessarily very different mixing (Wunsch 2003; Schmittner et al. 2015), and the than today (McManus et al. 2004; Ritz et al. 2013; Böhm nonlinear dependence of density on temperature (Winton et al. 2015). 1997; Sigman et al. 2004; De Boer et al. 2007). In turn, Looking further back in time to the Pliocene and pre- such physical oceanographic changes may have contributed ceding eras, before the northern hemisphere ice age cycles to changes in atmospheric C ­ O2 concentrations—hereafter began, ­CO2 was as high, or higher than today. It would referred to simply as “CO2”—measured in ice core records. appear that the production of dense waters from polar If so, the physical oceanographic changes would have been oceans was also different under these warmer states, and integral to the dynamics of ice age cycles (Sigman and changed during the onset of northern hemisphere glaciation Boyle 2000; Ganopolski et al. 2016). The close correlation 2.7 million years ago (Sigman et al. 2004, 2007), when δ13C of atmospheric C ­ O2 with ice core records of Antarctic tem- records imply that Antarctic deep water formation weak- perature over the past 800,000 years, combined with proxy ened alongside a ­CO2 drop (Hodell and Venz-Curtis 2006). evidence of accompanying changes in the Southern Ocean, Although tectonics may have played a role in these more has inspired many authors to link changes in Southern Ocean ancient changes, it is possible that a decrease in Antarctic circulation to oceanic ­CO2 sequestration (Paillard and Par- water formation was caused by the climatic effect of the C ­ O2 renin 2004; Fischer et al. 2010; Sigman et al. 2010; Skinner decrease, via its impact on surface ocean fluxes. Indeed, the et al. 2010; Jaccard et al. 2013; Ferrari et al. 2014; Watson sense of change is consistent with model simulations forced et al. 2015). by anthropogenic ­CO2 rise that have shown, after the rapid A number of popular hypotheses call on saline brine, transient warming, increased production of Antarctic deep rejected from a massive apron of sea ice around Antarc- water in the equilibrated warm state (Yamamoto et al. 2015). tica, to have produced an abundant source of dense water Yet this is in apparent contradiction with the evidence from during the last ice age that filled the global abyss (Adkins ice age cycles, in which Antarctic waters dominate the deep et al. 2002; Bouttes et al. 2010; Adkins 2013). The idea that ocean in the colder glacial state. Antarctic Bottom Water (AABW) was much saltier during Many prior studies have explored how these and other ice ages arose from reconstructions made by extracting fos- changes in steady-state ocean circulation might have sil pore-water from deep sea sediment cores (Schrag and resulted from C ­ O 2, orbital and ice sheet changes, using de Paolo 1993; Adkins and Schrag 2003). After account- General Circulation Models (GCMs). Most of these stud- ing for changes in global ice volume, these reconstructions ies either used coupled ocean and atmosphere GCMs, suggest that the salinity of the deep ocean was significantly forced by ‘time-slice’ combinations of radiation and elevated relative to the upper ocean during the LGM (Adkins topography (most commonly those of the mid-Holocene et al. 2002; Insua et al. 2014). An Antarctic origin of these and LGM, e.g. Hewitt et al. 2003; Braconnot et al. 2007), salty bottom waters would be consistent with reconstruc- or used ocean GCMs to which experimental changes of tions of the chemical composition of deep waters, primar- wind, heat, and freshwater fluxes were applied (e.g. But- ily reflected by isotopic ratios in the shells of fossil benthic zin et al. 2005; Toggweiler et al. 2006; Jansen 2017). foraminifera. Both the stable carbon isotope ratio, δ 13C These studies have identified and illustrated many poten- (Boyle and Keigwin 1982; Duplessy et al. 1984; Sarnthein tially important mechanisms. However, when multiple et al. 1994; Curry and Oppo 2005) and neodymium isotope models have been subjected to the same boundary condi- ratio, εNd (Piotrowski et al. 2005; Howe et al. 2016) indicate tions, they have not always produced the same results, that dramatic changes occurred in the Atlantic basin, with because of differences in model physics and numerics a volumetric expansion of the Antarctic chemical finger- (Otto-Bliesner et al. 2007; Muglia and Schmittner 2015), print under glacial conditions. Although the δ13C evidence insufficient model equilibration (Marzocchi and Jansen is not unambiguous (Gebbie 2014), when viewed together 2017), or both. Furthermore, the heterogeneity of prior with the spatial distribution of εNd changes, there is strong experiments makes it difficult to gain an appreciation for 13 Response of a comprehensive climate model to a broad range of external forcings: relevance for… 655 the relative importance—and interactions—of the vari- 2 Methods ous mechanisms over the full range of climatic states that have held sway over the past few million years. Addition- 2.1 Model ally, though simplified models have illustrated how the interaction of ocean circulation with the marine ecosys- The simulations use the coupled ocean–atmosphere–ice- tem may have driven the ice age cycles of atmospheric biogeochemistry model CM2Mc.v2. In brief, this includes: ­CO2 (Köhler et al. 2005; Bouttes et al. 2010; Hain et al. a finite volume atmospheric model, similar to that used 2010; Brovkin et al. 2012), comprehensive models have in the GFDL ESM2 models (Dunne et al. 2012); MOM5, thus far failed to do so without modifications (e.g. Sim- a non-Boussinesq ocean model with a fully-nonlinear mons et al. 2016; Buchanan et al. 2016). Presumably, equation of state, parameterizations for mesoscale (Grif- the mechanisms that drove the past changes will also be fies 1998) and submesoscale (Fox-Kemper et al. 2008) relevant for the long-term future of the deep ocean as it turbulence, and parameterizations for vertical mixing due adjusts to the high ­CO2 of the Anthropocene. to boundary layer turbulence (Large et al. 1994) and to the Here we aim to remedy the conceptual fragmentation breaking of internal tides above rough topography (Sim- of prior works by presenting a large number of equilib- mons et al. 2004); a sea-ice module; a static land module; rium ocean states in a single coupled ocean–atmosphere and a coupler to exchange fluxes between the components. GCM. This approach is inspired by the pioneering work It differs from the model described by Galbraith et al. of Manabe and Bryan (1985), who explored the response (2011) due to the following small changes: the codebase of an early ocean–atmosphere GCM to a range of six C ­ O2 was updated to the Geophysical Fluid Dynamics Labora- forcings. The model we use is CM2Mc, a coarser-reso- tory Siena version; the albedo of snow-covered sea ice lution configuration of the Geophysical Fluid Dynamics was increased to a more realistic value of 0.85 from 0.80; Laboratory’s ESM2M model (Dunne et al. 2012), which the bottom drag was increased from 0.001 to 0.002 (closer captures critical aspects of the Southern Ocean circulation to its default value); the background vertical mixing was (Russell et al. 2006). We simultaneously vary atmospheric reduced from 0.1 to 0.05 cm 2 s−1; additional cross-land ­CO2, orbital configuration and terrestrial ice sheet size in mixing in Indonesia was included to account for unre- a matrix of 40 experiments, mapping out the large-scale solved flows between islands. These changes altered the behaviour of the model over an unusually-wide parameter pre-industrial control state slightly, but the overall simula- space. The model does not use flux corrections, and freely tion characteristics remain similar to the GFDL ESM2M evolves to an equilibrium state determined only by the in spite of the 1.5- to threefold reduction of horizontal response of the internal model physics to the external forc- resolution (nominally 3° in the ocean and 3° in the atmos- ing variables. Thus, in contrast with experiments that arbi- phere, vs. 1° and 2°, respectively, in ESM2M). The ocean trarily alter surface boundary fluxes in order to illustrate and atmosphere comprise 28 and 24 layers, respectively. hypotheses, we show the model’s response to the external The model does not use flux corrections, and conserves forcings without preconceptions. Although the model is freshwater by routing precipitation over land to coastal subject to important biases, the simulation suite provides a ‘river’ discharge cells, including as ‘iceberg’ discharge comprehensive overview of how these three drivers could where that precipitation falls as snow. Sea ice has a bulk feasibly influence climate given Earth’s present tectonic salinity of 5 PSU, and its formation from seawater rejects configuration. brine in the top ocean layer. The high precision of the The simulations show a rich diversity of behaviour in freshwater balance is evidenced by the small drift of whole the ocean, atmosphere and land surface, among which we ocean salinity, which changes by < 0.001 PSU over the last briefly discuss general atmospheric aspects before focus- 5 centuries in all simulations. ing on the circulation of the global deep ocean. We also The model does not apply flux corrections, and there- compare the results with some key paleoceanographic fore has significant biases. Primary among these are a reconstructions from the LGM, showing encouraging con- global cold bias in the pre-industrial control surface air sistency with physical observations as well as important temperature (12 °C, rather than the estimated pre-indus- biogeochemical differences. We conclude by proposing trial average of roughly 14 °C), which is partly due to a simplifying principle, based on the relative density of the increased sea ice albedo mentioned above, but which polar waters, that makes testable predictions regarding the was deemed acceptable given that the global ocean aver- volume filling of the deep ocean under different climate age temperature bias is very small (+ 0.04 °C). The North states, and offer speculation regarding how this may have Pacific has a particularly cold bias, which tends to weaken contributed to the late Cenozoic intensification of the ice the stratification there and encourage more sea ice and ages. stronger North Pacific Intermediate Water formation than 13 656 E. Galbraith, C. de Lavergne observed. The Southern Ocean has a warm bias, though radioactive decay. Air–sea fluxes are calculated assuming a its hydrographic structure agrees well with observations ­CO2 partial pressure of 270 µatm and a pre-industrial atmos- (Galbraith et al. 2011). It should also be borne in mind that pheric Δ14C of 0‰. Because the biogeochemistry module the ocean model does not resolve geostrophic turbulence ­ O2 and 14CO2, it does uses fixed pre-industrial atmospheric C and therefore parameterizes the effect of mesoscale eddies. not include any air–sea exchange effect of variable pCO2 This parameterization (Gent and McWilliams 1990; Grif- (Bard 1988; Galbraith et al. 2015). Note also that the bio- fies 1998) is an imperfect surrogate for eddy effects whose geochemical model simulates changes in dissolved nutrients response to changes in the mean climate state is certainly and gases that are caused by the changes in ocean circulation biased (e.g. Munday et al. 2013). In addition, the model— only, thus excluding changes in dust fertilization or any other like all global climate models—is incapable of correctly ecosystem change not driven directly by the orbital, ­CO2 and simulating the flow of dense waters along steep topo- ice sheet forcing impact on ocean circulation. graphic slopes. As a result, convective chimneys, instead of downslope currents, are the dominant mode of dense 2.2 Simulations water injection to depth, biasing the properties and sensi- tivity of deep water masses to some degree. The simulations populate a matrix that is defined by three Water masses were ‘tagged’ in the surface layer in order external forcing variables (Table 1). First, atmospheric ­CO2 to unambiguously identify volume contributions from dif- was prescribed for each simulation, at levels ranging from ferent source regions. Three tags are defined: a Southern 147 to 911 ppm. The levels from 270 (pre-industrial) to 911 tracer south of 30°S, a North Atlantic tracer north of 30°N are each 1.5-fold greater than the prior level, thus presenting in the Atlantic, and a North Pacific tracer north of 30°N in a nearly-regular spacing of radiative forcing, while the four the Pacific. In the surface layer, the tracers are set to 1 within lowest (spanning the ­CO2 range over recent glacial cycles) their tagging region and to zero outside. In all other ocean are spaced at half that distance. Although we prescribe only layers, they are carried as passive tracers. We define AABW changes in ­CO2, the model response would not differ sig- and NADW volumes as the global volume integral of the nificantly if forced by other greenhouse gases (­ CH4, ­N2O) Southern and North Atlantic tracer concentrations, respec- at equivalent levels of radiative forcing. tively, below 1 km depth. Because the model’s analogue of Second, the two astronomical parameters were each set to Antarctic Intermediate Water does not penetrate deeper than two values, resulting in four different configurations (Fig. 1). 1 km, restricting the analysis to depths > 1 km effectively Obliquity, i.e. the tilt of the rotational axis, alters the lati- excludes the volume contribution of intermediate and lighter tudinal distribution of radiation across the surface of the waters. These definitions were chosen so as to be applica- Earth, increasing annual incident radiation at high latitudes ble under a wide range of climate states, and differ from under high obliquity. Obliquity was set to either 22.0° or classical water mass definitions based on density and basin 24.5°, spanning the calculated range of the last 5 million criteria: all basins and all ocean depths > 1 km contribute to years (Laskar et al. 2004). The precessional phase changes computed AABW and NADW volumes. Thus, AABW and the seasonal distribution of insolation at any given location, NADW are broadly defined as the deep waters last ventilated without altering the total annual incident radiation, and is in the Southern Ocean and North Atlantic, and their volumes defined here as the angle between the Earth’s position during essentially quantify the total volumetric influence of each the northern hemisphere autumnal equinox and the perihe- surface high latitude region, irrespective of ventilation rates lion (Held 1982). The precessional phase was set to the posi- and pathways. The temperature, salinity and in situ density tions at which the insolation seasonalities are most and least of the AABW and NADW end-members are calculated at severe: the most severe northern seasonality occurs at 270°, 100 m-depth, averaged over the Southern Ocean or North while the most severe southern seasonality occurs at 90°. Atlantic locations and years where the annual mean mixing It should be noted that these precessional phases will not layer depth exceeds 1 km. necessarily produce the most extreme responses in climate An ‘ideal age’ tracer is defined as zero in the global metrics, given that the seasonal cycles of insolation and sur- surface layer, and increases in all other layers at a rate of face temperatures are not perfectly in phase (Kutzbach et al. 1 year per year. The biogeochemical module analyzed here is 2008); nonetheless they provide strong and easily interpreted BLINGv0 (Galbraith et al. 2010). Radiocarbon is simulated forcings. The eccentricity was maintained at a value of 0.03, by two ocean tracers, Dissolved Inorganic Carbon and Dis- which is close to the mean for the range of 0.000–0.058 solved Organic Carbon (­ DI14C and D ­ O14C), and compared to spanned over the last 5 million years (Laskar et al. 2004) and the non-radiogenic tracers DIC and DOC. These tracers are nearly double the modern value of 0.01671. cycled implicitly through organic matter, assuming a fixed Third, simulations were run using non-interactive land C:P of uptake, but the 14C tracers are not cycled through cover type and topography corresponding to either pre- ­CaCO3 for simplicity. Both D ­ I14C and D ­ O14C undergo industrial or ‘glacial’ conditions. For the latter, the surface 13 Response of a comprehensive climate model to a broad range of external forcings: relevance for… 657 Table 1  Overview of model CO2-ice sheets Integration time Global average ∆14C drift, final century (last 5 centuries) simulations (prior spinup) o22, p270 o22, p90 o24.5, p270 o24.5, p90 −1 −1 −1 Years: ‰ ­century ‰ ­century ‰ ­century ‰ ­century−1 147-PI 3100 − 1.3 (− 1.5) − 1.3 (− 1.4) − 1.3 (− 1.4) − 1 (− 1.1) 180-PI 5000 − 0.6 (− 0.6) − 0.6 (− 0.7) − 0.6 (− 0.6) − 0.3 (− 0.4) 180-G/PI 4600 − 0.6 (− 0.7) − 0.8 (− 0.8) − 0.5 (− 0.6) − 0.4 (− 0.4) 180-G 6000 (3200) − 0.5 (− 0.5) − 0.5 (− 0.6) g − 0.5 (− 0.6) − 0.4 (− 0.4) 220-PI 5000 − 0.7 (− 0.8) − 0.6 (− 0.7) − 0.4 (− 0.5) − 0.3 (− 0.3) 220-G 6000 − 0.7 (− 0.8) − 0.6 (− 0.7) − 0.5 (− 0.6) − 0.3 (− 0.4) 270-PI 3200 − 0.8 (− 0.8) − 0.6 (− 0.5) − 0.6 (− 0.6) − 0.4 (− 0.5) ig 405-PI 3200 − 0.6 (− 0.6) − 0.4 (− 0.5) − 0.2 (− 0.3) − 0.2 (− 0.3) 607-PI 3200 (2000) − 0.3 (− 0.3) − 0.5 (− 0.5) − 0.3 (− 0.4) − 0.5 (− 0.6) 911-PI 3400 (2000) − 0.3 (− 0.4) − 0.4 (− 0.5) − 0.2 (− 0.3) − 0.3 (− 0.5) The first column lists the different configurations of ­CO2 (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 ini- tial 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 configu- ration, defined by the obliquity (o22 or o24.5) and precessional phase (p270 or p90). The generic ‘glacial’ and ‘interglacial’ states are indicated by boldface of 180 ppm, but for which the land surface type was kept at 90° pre-industrial, thereby allowing a separation of the effects 22° of ice sheet topography vs. albedo. Note that ice sheet sizes 24.5° are held constant throughout any given simulation here, perihelion assuming all mass additions (by snow) are balanced by abla- tion (iceberg calving at river discharge points) and thereby ignoring the possible impacts of transient ice sheet growth aphelion 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’) 270° experiments by Brown and Galbraith (2015) and Galbraith et al. (2016). We do not reduce ice sheet size or mean ocean Fig. 1  Orbital forcing configurations. Earth is shown at the closest salinity in the high C ­ O2 runs: only the roles of radiative and and furthest approaches to the sun (perihelion and aphelion, respec- orbital forcings are explored within the 270–911 ppm C ­ O2 tively, not to scale). The tilt of the Earth’s rotational axis (obliquity) range. 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 All simulations were initialized from either an equili- north pole points towards/away from the sun at perihelion/aphelion, brated pre-industrial control state or modern observa- boreal seasons are more intense than austral seasons (upward-point- tions, and integrated until deep ocean salinity gradients ing arrows) approached a stable quasi-equilibrium. Some of the simu- lations were branched from 2000-year integrations with topography was altered to reflect the geometry of terrestrial the same C ­ O2 and modern orbital configuration; because ice sheets as reconstructed for the LGM and land surface the ocean temperature was already somewhat equilibrated types were changed according to the Paleoclimate Mod- to the C­ O2, these required a shorter subsequent integration elling Intercomparison Project 3 (PMIP3) protocol, the with the four different orbital forcings (1200–1400 years). Bering Strait was closed, coastal bathymetry and land-sea Similarly, the ‘180-glacial’ simulations were integrated for masks were altered to prevent ocean cells occurring under- 3200 years with pre-industrial albedo, before swapping to neath terrestrial ice sheets, and average ocean salinity was the glacial land surface albedo. Total spinup times (start- increased by 1 PSU [that which would be produced by a ing from the pre-industrial control) were between 3100 110 m sea level fall, typical of glacial states, and ~ 0.1 PSU and 6000 years (Table 1). The simulations were run using less than estimated for the LGM (Adkins et al. 2002)]. We the Scinet general purpose cluster at the University of also include four simulations at glacial conditions and C ­ O2 Toronto between 2013 and 2015, at a computational cost 13 658 E. Galbraith, C. de Lavergne of 23 CPU hours per simulation year, requiring a total of 3 Results approximately 3.6 million CPU hours. The analyses pre- sented here are all averages of the final 100 years of the 3.1 Global atmospheric response simulations. Because of the very long timescale of radio- carbon equilibration, the simulations all have substantial Although this paper is focused on the ocean, we begin with downward drift in the global average Δ 14C of 0.2–1.3‰ a brief description of some general features of the global per century (Table 1). The fact that the downward drifts atmospheric response, since all oceanic changes are medi- are very similar for the ‘glacial’ and ‘interglacial’ states ated by the atmospheric response to the external forcings. we analyse closely (boldface in Table 1) suggests that the Figure 2 shows the global pattern of surface air tem- difference between them is likely to be representative of perature change driven by the two orbital forcing variables the steady state. However, this assumption cannot be tested of obliquity and precession, after full equilibration of the here, placing a caveat on the interpretation of the radiocar- climate system, including the slow adjustment of the deep bon contrasts (Sect. 4.4.3). ocean and high latitude clouds and sea ice. The two patterns 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 (a) Obliquity (low - high) Jansen (2017). The length of our simulations may be a major 80°N factor in explaining both their consistent behaviour and the relatively close agreement with observations. 40°N 2.3 Definition of ‘glacial’ and ‘interglacial’ states 0° Our simulations are designed to test model behaviour under external forcings that approximately bracket their variabil- 40°S ity 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 80°S of comparing interglacials and glacials in general, we take two representative sets of boundary conditions. The onset of (b) Precession (strong boreal - weak boreal seasons) full glacial conditions generally occurs under low obliquity, whereas interglacials generally occur under high obliquity 80°N (Huybers and Wunsch 2005). The chosen ‘glacial’ state is therefore a low obliquity simulation and the ‘interglacial’ 40°N is a high obliquity simulation, extrema that occurred most recently at 28 and 9 ka, respectively. Although deglacia- tions generally occur during the precessional phase with 0° intense boreal seasonality (Cheng et al. 2009), precession is less strongly tied to glacial vs. interglacial states, given 40°S that the duration of the states frequently exceeds half a pre- cessional cycle. We therefore choose the precessional state of 90° between boreal autumnal equinox and perihelion for 80°S both states, corresponding closely with both the present-day 50 °E 150 °E 110 °W 10 °W and the LGM. The ‘glacial’ state has an atmospheric ­CO2 concentration of 180 ppm and full LGM ice sheets, altered bathymetry and a 1 PSU increase in mean ocean salinity. -3 -2 -1 0 1 2 3 Where appropriate, we also compare observations with a Temperature anomaly (°C) pre-industrial control simulation, identical to the ‘intergla- cial’ simulation but using present-day orbital forcing (obliq- Fig. 2  Effect of orbital forcings on surface air temperatures. Each panel shows the surface air temperature anomaly produced by one uity of 23.439°, 102.932° longitude between perihelion and of the two types of orbital forcing, calculated from the average of boreal autumnal equinox). all simulations with pre-industrial ice sheets. a Low obliquity minus high obliquity. b Strong boreal seasons minus weak boreal seasons 13 Response of a comprehensive climate model to a broad range of external forcings: relevance for… 659 of change driven by the two types of forcing are very differ- The global temperature response to the forcings is subject ent, a factor which is not evident in the common simplifica- to model uncertainties, most notably the response of clouds. tion of orbital forcing to metrics such as peak summertime For reference, the equilibrium climate sensitivity of CM2Mc insolation at 65°N. Obliquity has a very strong impact on (the response to a C ­ O2 doubling from 270 to 540 ppm in the the zonal mean temperature profile: high latitudes are much modern orbital configuration) is 3.6 K (not shown). This cooler under low obliquity (i.e. a small tilt of Earth’s rota- value is somewhat higher than the 3.1 K/doubling sensitivity tional axis), intensifying the equator-to-pole gradient by of the one-degree ESM2M model (Frölicher et al. 2014). As more than 3 °C. Precession produces a more complex spa- shown by Fig. 3, the climate sensitivity increases with cool- tial pattern of temperature change. Note that the magnitude ing of the background state (i.e. the slope is steeper at low of the precessional response depends on eccentricity (set ­CO2), agreeing with the findings of Kutzbach et al. (2013), here at an approximate Quaternary average of 0.03) and has though this tendency may not be uniform among models therefore changed over time. and is sensitive to uncertain aspects of sea ice and clouds. Figure 3 shows the globally-averaged surface air tempera- The global mean surface air temperature difference ture of each final state in the suite of simulations. ­CO2 is the between the ‘glacial’ simulation and a pre-industrial (PI) leading control of the global mean temperature but orbital control simulation is 5.4 K. This can be compared with forcing also has a significant impact. This orbital impact is PMIP2 models, which simulated a range of 3.6–5.7 K for due only to the spatial and temporal distribution of inso- the LGM to PI change (Braconnot et al. 2007). PMIP2 mod- lation, since the total integrated insolation is unaffected els generally show a larger LGM-PI change than the corre- by either obliquity or precession. In general, the world is sponding proxy-based estimates of 4.0 ± 0.8 K (Annan and warmer when obliquity is high (orange) and when boreal Hargreaves 2013), although the latter might underestimate seasonality is weak (downward-pointing triangles). Ice the true change (Telford et al. 2013; Ho and Laepple 2015). sheets exert a significant cooling at the global scale, most of The relatively large glacial-interglacial contrast simulated which is attributed to the direct albedo impact rather than by CM2Mc, compared to the PMIP2 LGM-PI average, can dynamical topographic effects (compare filled, outlined, and be partly explained by the fact that the envelope of orbital unfilled triangles in the 180 ppm experiments). 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. 20 3.2 Changes in water mass volumes and the AMOC 18 The experiments show a large range in the volumetric con- Surface air temperature (C) 16 tribution of AABW to the deep (> 1 km) ocean (Fig. 4a), varying from 30 to 80%. The AABW volume shows a 14 strong dependence on ­CO2, with a striking nonlinearity: ig the AABW volume reaches a minimum at intermediate 12 low obliquity, strong boreal seasons ­CO2 (270–405 ppm) and increases toward both higher and low obliquity, strong austral seasons lower ­CO2. The effect of orbital forcing is also significant, 10 high obliquity, strong boreal seasons and displays a further nonlinearity with C ­ O2. For C ­ O2 of high obliquity, strong austral seasons 8 LGM topography, PI land albedo 270 ppm and below, low obliquity (blue) increases the Ant- LGM topography and albedo arctic ventilation volume, whereas only precession plays a g 6 role at higher C­ O2, producing a greater fraction of AABW 147 180 220 270 405 607 911 when boreal seasonality is weak (downward-pointing trian- CO 2 (ppmv) gles). The inclusion of large LGM ice sheets also increases the AABW fraction, due to both the topographic effect and Fig. 3  Globally-averaged surface air temperature vs. ­CO2. Each tri- the ice sheet albedo (compare filled, outlined, and unfilled angle shows the final state (100-year average) of one equilibrium triangles in Fig. 4a). North Atlantic ventilation volumes are simulation. Blue and orange symbols indicate low and high obliquity, respectively. Upward- and downward-pointing arrows indicate strong a mirror image of the Antarctic volumes (Fig. 4b), because boreal and austral seasons, respectively, according to the precessional North Pacific waters make up a small fraction of the global phase. Filled solid outlines indicate simulations using pre-industrial deep ocean in all simulations. Under low ­CO2, low obliq- land surface cover with LGM ice sheet topography, bathymetry, and uity and pre-industrial ice sheets the North Pacific fraction ocean salinity, while unfilled solid outlines replace the land surface type with an LGM estimate. Small letters indicate the ‘glacial’ (g) reaches a maximum, overlapping with the North Atlantic and ‘interglacial’ (ig) simulations discussed in the text fraction which is at its minimum. The North Pacific fraction 13 660 E. Galbraith, C. de Lavergne (a) Antarctic (b) North Atlantic 80 60 g 70 50 ig Volume fraction (%) Volume fraction (%) 60 40 50 30 40 ig 20 g 30 10 20 0 147 180 220 270 405 607 911 147 180 220 270 405 607 911 CO 2 (ppmv) CO 2 (ppmv) (c) North Pacific (d) AMOC 12 22 11 20 18 10 ig 16 g Volume fraction (%) 9 transport (Sv) 14 8 12 ig 7 10 6 8 5 g 6 4 4 3 2 147 180 220 270 405 607 911 147 180 220 270 405 607 911 CO 2 (ppmv) CO 2 (ppmv) Fig. 4  Global deep water volumes and AMOC. Panels show the sim- Pacific, and d the maximum AMOC transport at 40°N. All symbols ulated percentage of the deep (> 1 km) ocean volume occupied by are as defined in Fig. 3 waters originating in the a Antarctic, b North Atlantic and c North decreases both with increasing C ­ O2 and with larger ice oscillates between strong and weak states, reminiscent of sheets (Fig. 4c). Dansgaard–Oeschger variability, as discussed in detail by The flow of water transported by the AMOC is relatively Brown and Galbraith (2015). Note that, in all cases shown stable between 16 and 22 Sv for all simulations in which here, the freshwater input to ice sheets from precipitation is ­CO2 is 270 ppm or above, and varies with both precession balanced by calving in the ocean, so that we do not include any and obliquity (Fig. 4d). However, below 270 ppm the AMOC possible transient effects of ice sheet mass balance on ocean transports span a range of 3–22 Sv. Reduced AMOC transport circulation (e.g. Swingedouw et al. 2007). under low C­ O2 (150 ppm) was also simulated by Stouffer and Manabe (2003) using an earlier version of the GFDL coupled 3.3 Deep ocean ventilation rate model. The weakest AMOC transports occur under low obliq- and biogeochemistry uity with pre-industrial ice sheets, whereas strong transports occur under high obliquity with LGM ice sheets. In two of the Like the AABW volume fraction, the simulated ventila- intermediate-strength simulations the AMOC spontaneously tion rate of the global deep ocean (quantified as the mean 13 Response of a comprehensive climate model to a broad range of external forcings: relevance for… 661 14 ideal age below 1 km) has a nonlinear dependence on C ­ O2 C—that is, an increase in the radiocarbon age. Hence, (Fig. 5a). The mean age is lowest (i.e. the ventilation is most radiocarbon is not a straightforward monitor of ideal age rapid) under low ­CO2, passes through a maximum at inter- (Schmittner 2003; DeVries and Primeau 2010), although mediate ­CO2, and decreases somewhat at high ­CO2. Orbital it still contains valuable information, as discussed below. forcing has a large impact at low ­CO2, and very little impact Deep ocean oxygen is much more consistently linked at high C ­ O2. LGM ice sheets increase the ventilation rate of to ideal age than 14C: elevated age co-occurs with low the deep ocean significantly, mostly due to the topography oxygen content, and vice-versa (Fig. 5a, c). Southern effect (i.e. there is little difference between filled-outlined Ocean nitrate, which has been linked to the efficiency of and unfilled triangles at 180 ppm). the global biological carbon pump through its control of Perhaps surprisingly, whole ocean radiocarbon does preformed nutrients (Sarmiento and Toggweiler 1984; Sig- not follow the expectation given by ideal age (Fig. 5a, man and Boyle 2000), varies most strongly with obliquity b). Ideal age decreases with increasing C ­ O2 from 405 to at low C ­ O2, and generally decreases with increasing C ­ O2 911 ppm, but radiocarbon is relatively invariant across this (Fig. 5d). This is opposite to what would be expected from range. Further, whereas ideal age decreases with declining the argument of a Southern Ocean control of ice age ­CO2 ­CO2 below 405 ppm, the simulations show a decrease of Fig. 5  Deep ocean ventilation rates and biogeochemistry. a Mean ideal age deeper than 1 km. b Mean radiocarbon content (Δ14C, 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 13 662 E. Galbraith, C. de Lavergne tied to a reduction of preformed nutrients via reduced ver- tical exchange (Sigman et al. 2010). 4 Discussion Our results show large and nonlinear responses of AABW volume, AMOC and ocean ventilation to ­CO2, 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 cir- culation. 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 increas- ing sensitivity of the climate system to orbital forcing under decreasing ­CO2 that may have contributed to the develop- ment of ice age cycles. 4.1 CO2 and orbital influence The simulated changes in AABW vs. NADW volume frac- tions can be predicted, to first order, by the density differ- ence between the two formation regions (Fig. 6). At interme- diate ­CO2 (270–400 ppm), for which the AABW volume is Fig. 6  Antarctic and North Atlantic volume contributions to the deep minimal, the density at deep water formation regions of the ocean as a function of the Antarctic–North Atlantic density differ- North Atlantic tends to match or slightly exceed its Antarctic ence. Simulations with LGM (pre-industrial) topography are shown as unfilled (filled) circles. Symbol size reflects atmospheric C ­ O2 (in counterpart. At higher and lower ­CO2, for which AABW ppm) while symbol colours reflect the Antarctic-North Atlantic a dominates, the density of surface waters in AABW forma- salinity and b temperature difference. Densities, temperatures and tion regions substantially exceeds that of NADW formation salinities are averaged at 100 m depth over the Southern Ocean or regions. The fact that similar volume fractions co-occur with North Atlantic locations and years where deep (> 1 km) convection occurs similar density differences despite very different combina- tions of temperature and salinity (colours in Fig. 6) under- lines the predictive strength of the density difference. of deep convection sites in the North Atlantic causes them to In contrast, the volume changes do not correlate with the decrease in density more rapidly than their Antarctic coun- strength of southern hemisphere westerlies, though they do terparts, despite an increase of salinity. Meanwhile, with correlate with the volume transport of the Antarctic Circum- decreasing ­CO2 below 270 ppm, Antarctic waters grow rela- polar Current (ACC, Fig. 7). This indicates that changing tively denser despite North Atlantic cooling because of large westerlies are not a main driving force behind the simulated ­ O2 (< 180 ppm), low salinity increases. Note that three low C changes in volume fractions. Rather, it is changes in surface obliquity experiments with pre-industrial ice sheets make buoyancy forcing, largely independent from the winds, that an exception to the general North Atlantic trend: unusually control the variations in water mass volumes. The large sim- warm and salty NADW (Fig. 8) results from an equatorward ulated changes in density reflect changes in both temperature migration of the North Atlantic convection centre, following and salinity, which we discuss in turn. a sea ice expansion (not shown). Despite this large change in The temperature ranges of waters in convective regions convective regions, the density difference still explains the of the North Pacific and North Atlantic span 12 °C among relative volume fractions (Fig. 6). the simulations, whereas the Antarctic convection sites Salinity is influenced at the sea surface by the transport remain close to freezing in most simulations, and only rise of freshwater in the atmosphere and by sea ice cycling (e.g., 5 °C above the freezing point in the highest ­CO2 simula- Durack et al. 2012; Goosse and Fichefet 1999). At steady tions (Fig. 8). Thus, at high ­CO2, the more rapid warming state, these surface freshwater fluxes are balanced by the 13 Response of a comprehensive climate model to a broad range of external forcings: relevance for… 663 (a) Antarctic Circumpolar Current Polar ocean deep convection regions 80 12 g 10 26.0 * 70 26.5 27.0 * * 8 60 AABW fraction (%) temperature 28.0 6 50 4 40 ig 2 27.5 28.5 30 0 -2 20 33.5 34 34.5 35 35.5 36 120 140 160 180 200 220 salinity ACC transport (Sv) (b) Southern westerly winds Fig. 8  Surface properties in convective regions of the three polar oceans. For each simulation there is one value for each of the North 80 Pacific, Antarctic and North Atlantic, calculated from 100-year time g series. Warm/cool shades indicate high/low obliquity, triangle ori- entation indicates precessional phase as in Fig. 3, larger symbol 70 size indicates higher atmospheric C ­ O2, solid fills indicate preindus- trial topography and outlines indicate LGM topography and albedo. 60 Temperatures and salinities, at 100 m depth, were averaged within AABW fraction (%) each basin at the locations and years where the deepest convection occurred (mixing layer depth > 1 km or > 90% of its overall maxi- 50 mum). Under pre-industrial conditions, Antarctic and North Atlantic waters have similar densities, whereas they diverge under both high and low ­CO2. The North Pacific has lower density than the Antarc- 40 ig tic and North Atlantic under all simulated conditions. Small asterisks indicate the three simulations with low-latitude North Atlantic deep 30 convection, discussed in the text 20 7.5 8 8.5 9 9.5 10 Ocean and North Pacific but increases the Atlantic salin- Westerly windspeed ity (Fig. 9). In the model experiments, the net impact of lowering ­CO2 on precipitation and evaporation is therefore Fig. 7  AABW fraction vs. transport by the Antarctic Circumpolar to decrease the surface freshening of the Southern Ocean Current (ACC) and vs. the strength of Southern Hemisphere wester- lies. Symbols are defined as in Fig. 3. In general, the ACC transport is relative to the North Atlantic. increased under low C ­ O2 due to the strong meridional density gradi- Even more important here are changes in the extent of ent induced by Antarctic sea ice export, and at high ­CO2 due to the southern hemisphere sea ice. As discussed in earlier studies thermal density gradient. The mean surface velocity (in m s­ − 1) of the (Keeling and Stephens 2001; Shin et al. 2003; Ferrari et al. southern westerlies does not show a significant relationship with the AABW fraction 2014; Watson et al. 2015; Jansen and Nadeau 2016; Klock- mann et al. 2016), low atmospheric temperatures encourage sea ice production and export, extracting freshwater around net freshwater transport by ocean currents (Nilsson et al. Antarctica (where freezing dominates) and releasing it fur- 2013; Cessi and Jones 2017). In the model, as observed ther north (where melting dominates). Sea ice also acts as a (Durack et al. 2012), the high latitude Southern Ocean and barrier to precipitation, collecting snowfall and preventing North Pacific are net ‘freshening’, receiving more freshwater it from entering the ocean until the ice beneath it melts. from precipitation than they lose to evaporation, while the Because of the open geometry and strong zonal winds of the Atlantic basin as a whole is net evaporative, due to atmos- Southern Ocean, Ekman transport can potentially move large pheric freshwater export. Thus, a stronger atmospheric volumes of sea ice across long distances, unlike the North hydrological cycle, as occurs under warmer planetary tem- Atlantic where the basin shape limits the ability to export sea peratures given the Clausius–Clapeyron relationship (Held ice, causing the meridional transport to be smaller (Goosse and Soden 2006), freshens surface waters in the Southern and Fichefet 1999; Shin et al. 2003; Lee et al. 2017). Thus, 13 664 E. Galbraith, C. de Lavergne (a) Net atmospheric freshwater fluxes, Atlantic basin (b) Net atmospheric freshwater fluxes, south of 40 °S -0.35 1.15 1.1 Atmospheric freshwater input + runoff (Sv) -0.4 Net freshwater input + runoff (Sv) 1.05 -0.45 ig g 1 -0.5 0.95 -0.55 0.9 -0.6 0.85 -0.65 ig 0.8 -0.7 0.75 g -0.75 0.7 -0.8 0.65 147 180 220 270 405 607 911 147 180 220 270 405 607 911 CO 2 (ppmv) CO 2 (ppmv) (c) Southern Ocean sea ice transport (d) Sea surface freshwater flux south of 60°S 0.8 0.4 0.3 0.7 g Freshwater transport by sea ice (Sv) Net freshwater input + runoff (Sv) 0.2 ig 0.6 0.1 freshening 0.5 0 salinifying 0.4 -0.1 -0.2 0.3 -0.3 0.2 ig -0.4 g 0.1 -0.5 0 -0.6 180 220 270 405 607 911 147 180 220 270 405 607 911 CO 2 (ppmv) CO 2 (ppmv) Fig. 9  Surface freshwater balance of the Atlantic and Southern Ocean, south of 40°S. c Total freshwater transport by Antarctic sea Oceans. a Net atmospheric freshwater input to the Atlantic Ocean, ice. d Total freshwater input to the ocean south of 60°S, including including the Mediterranean and Arctic (note negative values indicate all air–sea exchange and sea ice cycling. Symbols are defined as in net evaporation). b Net atmospheric freshwater input to the Southern Fig. 3. All quantities are annually-averaged when low ­CO2 cools the global atmosphere, sea ice forms waters close to Antarctica but freshens the lower-latitude more rapidly at high latitudes of the Southern Ocean and is intermediate and mode waters—whose northward transport transported further north, causing the high latitude Southern ultimately feeds North Atlantic sinking. Thus, when Ant- Ocean to transition from a net-freshening to a net-salinifying arctic waters become saltier and denser, and consequently domain (Fig. 9c, d). Low obliquity (blue) also increases the expand their domain in the deep ocean, they force a reduc- growth of sea ice, since it leads to colder temperatures at tion in the salinity and therefore density of North Atlantic high latitudes. waters, amplifying the impact of the initial density gain. The Because the model does not use freshwater flux correc- volume of Antarctic sea ice shows a strong, approximately tions, we can observe an additional impact of the Antarctic linear relationship with the Antarctic-North Atlantic salinity salinity increase: for a given amount of salt in the ocean, if difference (Fig. 10). deep waters formed in the Antarctic become saltier, deep The salt balance can also explain part of the observed waters formed in other regions (primarily the North Atlantic) changes in the strength of the AMOC at low C ­ O2. The must become fresher. Amplified sea ice export salinifies the freshening of the North Atlantic surface brings the Atlantic 13 Response of a comprehensive climate model to a broad range of external forcings: relevance for… 665 circulation closer to a bifurcation point where the AMOC is Battisti 2008). However, over the Labrador Sea, the ice sheet bi-stable: the poleward advection of salty subtropical waters topography generates a strong northwesterly flow, accompa- becomes prone to interruption, allowing the subpolar gyre nied by a very strong cooling (Broccoli and Manabe 1987) to develop a strong halocline and impeding the formation that is amplified by the ice sheet albedo effect (Fig. 11c). of dense waters (Stommel 1961). This behaviour permits An anomalous northward wind stress occurs in the north- unforced oscillations of the AMOC to develop in two of the east Atlantic. The wind shift and cooling appear to drive the simulations at 180 ppm C ­ O2 as discussed in detail by Brown AMOC strengthening shown in Fig. 4d, as has previously and Galbraith (2015), and provides a potential explanation been observed in other models (Pausata et al. 2011; Zhang for the fact that the Dansgaard–Oeschger “bi-polar seesaw” et al. 2014; Zhu et al. 2014; Gong et al. 2015; Muglia and variability does not occur under the high ­CO2 of intergla- Schmittner 2015; Klockmann et al. 2016; Sherriff-Tadano cial periods (McManus et al. 1999; Rahmstorf 2002; Buiz- et al. 2017). When compared with the influences of ­CO2 and ert and Schmittner 2015; Kawamura et al. 2017). However, orbital forcing, this suggests that Dansgaard–Oeschger-like Dansgaard–Oeschger variability is not observed under full AMOC disruptions would be most likely to have occurred glacial maxima either (McManus et al. 1999; Kawamura in a ‘sweet spot’ of low ­CO2, low obliquity and a relatively et al. 2017), a fact that can be explained by the simulated small (low-elevation) Laurentide ice sheet. We therefore ice sheet effects. hypothesize that, under these conditions, the AMOC was more sensitive to variations in freshwater forcing arising 4.2 Effects of Laurentide ice sheet topography from tropical climate variability or freshwater discharge and albedo from ice sheets (Schmittner and Clement 2002). Counter intuitively, the ice sheet-driven AMOC strength- Although the AMOC is weakened by the inter-basin salt ening does not increase the volume contribution of NADW balance under low ­CO2 and low obliquity, this tendency is to the global ocean, but instead increases the contribution counteracted in all simulations that include the full LGM ice of AABW. Whereas low ­CO2 has mixed effects on the tem- sheets. Figure 11a confirms that low ­CO2 causes the mixed perature of NADW (Fig. 12a), the intense convection under layer depths to shoal in the northeast Atlantic, as a strength- extremely cold atmospheric conditions induced by the Lau- ened halocline impedes convection. As previously noted by rentide ice sheet cools NADW dramatically (Fig. 12b, c). other authors, steering of the atmospheric circulation by a The cold NADW flows at < 2.5 km depths through the Atlan- large Laurentide ice sheet (Fig. 11b) strengthens the zonal- tic and feeds into waters upwelling in the ACC, similar to the ity of the westerly jet along the southern margin of the ice mechanism proposed by Crowley (1992). A further contribu- sheet and central Atlantic (Roe and Lindzen 2001; Li and tion 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 1.2 of the near-surface Southern Ocean increases the produc- 1 tion of sea ice by reducing the heat loss required for winter g ice growth (Martinson 1990) and increases the northward 0.8 export of Antarctic sea ice by translating the melting regions Salinity difference 0.6 northward. Thus, like low ­CO2 and low obliquity, a large Laurentide ice sheet increases Antarctic salinity and ulti- 0.4 mately the density and volume of AABW (Figs. 4, 8, 9). This 0.2 mechanism is similar to previous suggestions that colder 0 NADW may cool the Southern Ocean remotely (Gildor et al. 2002) and that it may increase Antarctic salinity by chang- -0.2 ing the melting rates of land-based glaciers (Adkins 2013). -0.4 The simulations here confirm the importance of the remote ig cooling effect, but its impact on salinity occurs via Antarctic -0.6 0 5 10 15 20 25 30 35 40 sea ice formation and export. Annual average ice volume, Southern Ocean (103 km3 ) The overall result is a decoupling between meridional volume transport and water mass geometry, showing that Fig. 10  The inter-basin salinity gradient and Southern Ocean sea ice. a strong glacial AMOC can feasibly coexist with, and Difference in salinity between Antarctic and North Atlantic deep con- even contribute to, an expanded AABW. Importantly, this vection sites vs. the total annually-averaged Antarctic sea ice volume. decoupling shows that variations in the AABW and NADW A larger Antarctic sea ice volume is accompanied by larger northward freshwater transport, increasing the salinity of AABW and indirectly steady-state circulations cannot be universally reduced to decreasing the salinity of NADW an overturning see-saw: as we discuss further below, a large 13 666 E. Galbraith, C. de Lavergne (a) Glacial CO2 (b) Glacial topography 70° 70° 60° 60° 50° 50° 400 40° 40° 300 30° 30° 200 100 20° 20° 80°W 40°W 0° 40°E 80°W 40°W 0° 40°E 0 -100 (c) Glacial albedo (d) All glacial -200 70° 70° -300 60° 60° -400 50° 50° 40° 40° 5.0 m s-1 30° 30° 20° 20° 80°W 40°W 0° 40°E 80°W 40°W 0° 40°E Fig. 11  CO2 and ice sheet impacts on convection and surface winds ence of ice sheet topography (LGM topography minus PI topogra- in the North Atlantic. Panels show the glacial-interglacial change in phy). d Influence of ­CO2, ice sheet topography and albedo. Low ­CO2 mixed layer depth (m, shaded) and winds (arrows) averaged over the impedes deep convection in the northeast Atlantic, which we attribute 4 orbital configurations, due to individual forcing aspects (a–c) and to the inter-basin salt balance, while the ice sheet topography causes due to all (d). a Influence of low ­CO2 (180 ppm minus 270 ppm). b deep convection to greatly increase in the Labrador Sea in response Influence of ice sheet albedo (LGM albedo minus PI albedo). c Influ- to the wind change, an effect that is accentuated by the albedo change Laurentide ice sheet size leads to a simultaneous increase in difference in the mechanistic pathways by which density var- the formation rates and volume transports of both NADW iations occur. In general, the greater the density of AABW and AABW. relative to NADW, the more AABW dominates the ocean’s deeper layers. Our results also show that trends in volume 4.3 The unifying role of relative polar ocean contributions do not simply reflect a see-saw between north- densities ern and southern overturning, since the AMOC does not weaken at high C ­ O2 despite the shrinking NADW fraction, As we have shown, the general trends in deep ocean circula- and since the AMOC and NADW fraction respond differ- tion can be rationalized by the influence of surface buoyancy ently to the inclusion of a large Laurentide ice sheet. fluxes on the surface density difference between the Ant- These results are subject to important caveats: climate mod- arctic and North Atlantic dense water formation regions. els are incapable of simulating many aspects of dense water In contrast, the model does not show significant changes in formation and sinking, and poorly represent the mechanisms deep ocean ventilation in response to changes in southern that consume dense waters at depth. However, we expect these westerly winds, in disagreement with some prior sugges- biases to affect ventilation rates and pathways more strongly tions (e.g. Toggweiler et al. 2006; Anderson et al. 2009). The than the global volumetric contributions diagnosed here. We unifying role of relative density variations encapsulates the also expect the link between volumetric contributions and changes brought about by all three climate forcings, despite the inter-basin density competition to be relatively robust to 13 Response of a comprehensive climate model to a broad range of external forcings: relevance for… 667 Fig. 12  CO2 and ice sheet impacts on global ocean temperature are zonal averages from 50°W to the prime meridian (Atlantic) and and salinity. Panels show the change in temperature (top) and salin- from 180°W to 140°W (Pacific). The increases in AABW salinity, in ity (bottom) averaged over the 4 orbital configurations, due to indi- response to all three forcings, are primarily attributable to increases vidual forcing aspects (a–c) and due to all (d), as defined in Fig. 11, of sea ice formation and transport in the Southern Ocean along sections through the Atlantic and Pacific oceans. The sections biases in the location and mechanism of dense water forma- 4.4 ‘Glacial’ vs. ‘interglacial’ states tion. Indeed, the obtained relationship between volumes and formation-site densities holds across the vast parameter space 4.4.1 Glacial water mass characteristics explored (Fig. 6), despite major variations in the locations of deep water formation, in the global temperature and salinity We first compare the ‘glacial’ ocean temperature and salin- distributions and in the water-mass transport rates. Addition- ity to porewater and ice core-derived reconstructions of the ally, the relationship echoes earlier results obtained with sim- LGM ocean (Fig. 13). The glacial simulation reproduces the plified ocean models (Stocker et al. 1992; Fichefet et al. 1994) remarkable switch from a thermally-stratified to a salinity- and idealized eddy-resolving numerical experiments (Wolfe stratified deep ocean identified by Adkins et al. (2002). The and Cessi 2009, 2010, 2011). In particular, the volumetric simulated glacial salinity between 3500 and 5500 m in the expansion of NADW as its formation density approaches that Pacific Ocean is 36.2, indistinguishable from the LGM salin- of AABW appears to be consistent with deeper NADW pene- ity of 36.1 ± 0.1 reconstructed from porewater extractions at tration for a larger window of shared surface densities between five sites spanning the same depth range (Insua et al. 2014). the North Atlantic and the Southern Ocean, as proposed by We are not aware of any other comprehensive climate model Wolfe and Cessi (2009, 2010, 2011). that has previously reproduced the reconstructed salinities of 13 668 E. Galbraith, C. de Lavergne (a) Adkins reconstructions Atlantic vertical salinity gradient (Jansen 2017). There is also a large contrast in the depth to which NADW pen- 5 etrates (Fig. 14c). In the interglacial case, the full water 4 Interglacial model column of the North Atlantic is filled by NADW. This observations remains true in the glacial state for the Nordic and Lab- 3 rador seas, but further south the contribution of NADW Potential temperature (C) dwindles rapidly below 2 km depth. This simulated con- .0 37 2 traction of NADW co-occurs with the volumetric expan- sion of underlying Antarctic waters, and is consistent with .5 .0 .5 38 38 37 1 the many paleoceanographic observations for shoaled .0 39 0 NADW in the glacial North Atlantic (Curry and Oppo 2005; Adkins 2013; Lippold et al. 2016), although it does -1 Glacial not preclude the possibility that part of the apparent shoal- ( ) ing reflects a change in the remineralized carbon content of -2 AABW (Gebbie 2014; Howe et al. 2016). We also note the 34.5 35 35.5 36 36.5 37 37.5 interesting possibility that, under glacial conditions, the Salinity Labrador Sea produced dense bottom waters that trickled downslope analogous to the modern Antarctic (Keigwin (b) Insua reconstructions 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 con- 34.5 35 35.5 36 36.5 37 37.5 tribution of North Atlantic waters to the Indo-Pacific in Salinity the glacial state, arguing against the idea of North-Atlan- tic ventilation of the southern source waters (Kwon et al. Fig. 13  Simulated and reconstructed deep ocean temperature and salinity. Outlined circles show glacial reconstructions in blue and 2012). modern values at the same locations in red. Filled circles show the The glacial vs. interglacial contrast is consistent with simulated glacial and interglacial temperatures and salinities at the theoretical suggestions (Ferrari et al. 2014; Watson et al. same locations. a Temperature and salinity reconstructions of Adkins 2015; Jansen 2017; Marzocchi and Jansen 2017) that et al. (2002). b Salinity reconstructions of Insua et al. (2014). The model captures the basic transition from the modern, thermally- Southern Ocean surface buoyancy fluxes explain the stratified deep ocean to the glacial, salinity-stratified deep ocean. The expansion of Antarctic deep waters in the global ocean, Shona Rise site, for which the salinity reconstruction has been ques- as detailed above. However there is also a key difference tioned as discussed in the text, is identified by parentheses 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 the glacial ocean with comparable accuracy (e.g. see Fig. S7 under the glacial state, following a simple geometric con- of Bopp et al. 2017). The simulated salinity at Shona Rise straint. In our simulations, the corresponding isopycnal is only 36.4, rather than the reconstructed 37.1 ± 0.2, but surface does not move appreciably (Fig. 15), which we the robustness of the Shona Rise reconstruction has been attribute to the impact of southern subpolar gyres on the questioned (Wunsch 2016) so we are unsure whether or not Southern Ocean density structure. The gyre circulations this discrepancy is significant. The model also simulates a act to homogenize density south of the ACC. In the inter- change in Mean Ocean Temperature of 2.47 °C between glacial state, sea ice cycling and surface buoyancy forcing the pre-industrial control and the ‘glacial’ state, which is are moderate and the gyres maintain a roughly uniform indistinguishable from the recent estimate of 2.57 ± 0.24 zonal mean upper-ocean density south of 60°S (Fig. 15a, ºC derived from noble gas ratios in Antarctic ice (Bereiter c, e). As a result, the isopycnal connecting with the zero- et al. 2018). Thus, we are unable to identify a significant bias buoyancy-flux line does not dive to the deep ocean until in the model simulation of the glacial-interglacial change in it reaches about 60°S, near the northern edge of the gyres deep ocean temperature and salinity. (Fig. 15c, e). In the glacial, the zero-buoyancy-flux line Figure 14 shows a global slice through the ‘glacial’ and shifts to a more northern location (~ 60°S; Fig. 15b), as ‘interglacial’ oceans. The deep ocean of the interglacial predicted by Ferrari et al. (2014), but since this simply state is stably stratified in temperature and, in the Atlan- moves to the northern edge of the gyre domains, the lati- tic, unstably stratified in salinity. By contrast, the glacial tude at which the dividing isopycnal plunges to depth does ocean is essentially isothermal beneath 1.5 km, and the not change appreciably (Fig. 15d, f). The absence of a salinification of glacial Antarctic waters reverses deep strong upward displacement of the dividing isopycnal in 13 Response of a comprehensive climate model to a broad range of external forcings: relevance for… 669 Interglacial Glacial (a) Temperature (C) 14 1 km 2 km 10 3 km 6 4 km 2 5 km -2 60N 30N eq 30S 60S 60S 30S eq 30N 60N 60N 30N eq 30S 60S 60S 30S eq 30N 60N (b) Salinity 1 km 36.7 2 km 35.9 3 km 35.1 4 km 34.3 5 km 60N 30N eq 30S 60S 60S 30S eq 30N 60N 60N 30N eq 30S 60S 60S 30S eq 30N 60N (c) North Atlantic water fraction 0.9 1 km 0.7 2 km 0.5 3 km 0.3 4 km 5 km 0.1 60N 30N eq 30S 60S 60S 30S eq 30N 60N 60N 30N eq 30S 60S 60S 30S eq 30N 60N (d) Pre-formed radiocarbon age (y) 1650 1 km 2 km 1250 3 km 850 4 km 450 5 km 50 60N 30N eq 30S 60S 60S 30S eq 30N 60N 60N 30N eq 30S 60S 60S 30S eq 30N 60N (e) Ideal age (y) 1400 1 km 1000 2 km 3 km 600 4 km 200 5 km 60N 30N eq 30S 60S 60S 30S eq 30N 60N 60N 30N eq 30S 60S 60S 30S eq 30N 60N N Atlantic Southern Ocean N Pacific N Atlantic Southern Ocean N Pacific Fig. 14  Sections through the Atlantic and Pacific oceans of (left) ‘interglacial’ and (right) ‘glacial’ simulations. The sections are zonal averages from 50°W to the prime meridian (Atlantic) and from 180°W to 140°W (Pacific) the model is robust using different measures of density and Toggweiler 1984; Toggweiler 1999; Sigman et al. 2010; (not shown) and considering seasonal effects (Fig. 15c–f, Skinner et al. 2010). If deep waters communicate more black and grey). Thus, we do not find a straightforward slowly with atmosphere-equilibrated surface waters, they relationship between the latitude of the zero-buoyancy- can accumulate more carbon from sinking organic detritus, flux line and the depth of its matching isopycnal in the thereby lowering atmospheric C ­ O2. Indeed, as shown by deep ocean. Eggleston and Galbraith (2018), the soft tissue pump stor- age in CM2Mc is well-approximated by a linear product 4.4.2 Glacial deep ocean ventilation and diapycnal mixing of the global ventilation age and the global export produc- tion. Many have argued that the low C ­ O2 of glacial periods The ventilation rate of the deep ocean can significantly alter was due, at least in part, to sluggish deep ocean ventilation carbon storage by the biological soft tissue pump (Sarmiento 13 670 E. Galbraith, C. de Lavergne Fig. 15  The zero-buoyancy-flux line and its connection with the deep c–f. Key features include: (1) upper-ocean density does not increase ocean under (a, c, e) interglacial and (b, d, f) glacial conditions. The monotonically with latitude; (2) the zero-buoyancy-flux line crosses annual mean potential density referenced to 1000 dbars (σ1) is shown a significant range of surface densities, particularly in the interglacial (a,b) at the surface of the Southern Ocean and as a zonal mean within state; (3) the depth of the mean dividing isopycnal at low latitudes is the c, d Indo-Pacific and e, f Atlantic basins. The black curves in a, comparable in the interglacial and glacial. These features are robust to b mark the southernmost latitude of zero annual mean surface buoy- different measures of density and result mainly from the presence of ancy flux. The (black) zonal annual and (grey) zonal austral-winter subpolar gyres means of σ1 along this zero-buoyancy-flux line are contoured in (Toggweiler 1999; Galbraith et al. 2007; Sigman et al. 2010; allows the interior buoyancy supply to bottom waters, and Skinner et al. 2010; Menviel et al. 2017). their resultant upwelling rate, to match the strong surface The model simulations do not accord with this picture of buoyancy loss and convective downwelling occurring around a sluggish deep ocean. In the glacial state, although most of Antarctica. the global ocean is ventilated from the Antarctic, it actually The maintenance of a rapid circulation across large gla- does so at a rate that exceeds the combined input of North cial density contrasts thus hinges on a roughly threefold Atlantic plus Antarctic waters under the interglacial state, increase of the deep mixing energy supplied implicitly by leading to very low ideal ages in the deep Southern Ocean tides. On the one hand, this does happen to be compatible (Fig. 14). with stronger tidal energy input to the deep sea under lower This remarkably rapid ventilation arises as a consequence sea level, as might be expected from the fact that dissipation of the strong Antarctic density gradients and the model’s on shelves was less due to sea level fall during the glacial diapycnal mixing parameterization (Simmons et al. 2004). (Wunsch 2003), and as simulated with tidal models of the The large density excess of Antarctic-sourced bottom waters LGM (Schmittner et al. 2015). Sedimentological features relative to other deep water masses increases the deep strati- have suggested that fast bottom flows existed in the abys- fication, and the mixing parameterization assumes that dif- sal ocean during the glacial (Hall et al. 2001; Kienast et al. fusive buoyancy fluxes, driven by tides, are proportional to 2007; Krueger et al. 2012) which could be consistent with the local bottom stratification. This stratification-dependence a strong AABW source. On the other hand, the magnitude, 13 Response of a comprehensive climate model to a broad range of external forcings: relevance for… 671 distribution and stratification-dependence of the parameter- of atmospherically-produced 14C to waters forming AABW, ized glacial tidal mixing in the model are certainly incorrect further increasing the pre-formed age of AABW. In addition, to a significant degree (MacKinnon et al. 2017). the volumetric increase of “old” AABW, at the expense of In summary, a strong interior buoyancy source is required “young” NADW, drives a whole-ocean shift to higher pre- to balance the large surface buoyancy loss in the glacial Ant- formed age (Schmittner 2003). arctic freezing zone at steady state, but the representation of The overall result of the changes in ocean circulation and subsurface ocean mixing in this and other climate models sea ice is an increase in the pre-formed ventilation age of the remains a major source of uncertainty for how and where ocean by 440 years between the glacial and pre-industrial this occurs, and for the consequent relationship between sur- states, as well as can be determined given the remaining drift face forcing and the deep ocean ventilation. Further work is in Δ14C (cf Sect. 2.2). In other words, if the ideal age of the needed to evaluate the assumptions behind the ocean mixing ocean had been exactly the same between the two states, the parameterizations, and the simulated ventilation rates should change in pre-formed reservoir ages would have made the be viewed with caution. Nonetheless, we can compare the ocean more depleted in radiocarbon by an amount equiva- simulated radiocarbon with observations to evaluate the lent to roughly 440 years of decay. When this is combined modelled ventilation rate. with the pCO2 effect on air–sea equilibration of radiocar- bon, expected to be approximately 250 years for the LGM 4.4.3 Glacial radiocarbon simulation (Galbraith et al. 2015), it implies an apparent aging of the glacial ocean by 690 years relative to the pre-industrial. This One of the main pillars of observational support for sluggish is very similar to the whole-ocean average radiocarbon age ventilation of the deep glacial ocean has been the apparently increase of 600 years reported by Sarnthein et al. (2013) and low radiocarbon content of benthic foraminifera in glacial- 700 years reported by Skinner et al. (2015). aged sediments, when compared to the contemporary atmos- Thus, comparison of the simulated radiocarbon with the phere (Galbraith et al. 2007; Skinner et al. 2010; Burke and global compilations of radiocarbon reconstructions suggests Robinson 2012; Sarnthein et al. 2013). Radiocarbon is pro- that the true ideal age of the LGM ocean could have been duced exclusively in the atmosphere, and is supplied to the similar to the interglacial, with the apparent difference in 14C ocean through air–sea exchange, so the benthic-atmosphere age caused by a reduction of the pre-formed 14C activity of radiocarbon gradient is typically interpreted as representing AABW, its volumetric dominance of the glacial ocean, and the ‘ventilation age’ of bottom waters, often linked to the the pCO2 effect. The fact that the simulated ideal age under rate of ocean overturning. However, as previously shown, the glacial state is roughly 200 years younger than that of the large changes in radiocarbon could have also occurred due interglacial state (Figs. 5 and 14), rather than the same (as to changes in air–sea flux and ocean circulation patterns, suggested by the radiocarbon), would be consistent with a rather than overturning rate (e.g. Schmittner 2003; Butzin significant overestimate of tidal mixing in the glacial simula- et al. 2005; Devries and Primeau 2010). Our model simula- tion, as suggested above in Sect. 4.4.2. An additional, non- tions show that, under a broad range of forcings relevant exclusive possibility is that a significant decrease in ventila- to the late Quaternary, changes in air–sea flux and ocean tion occurred following the obliquity minimum at 28 ka, so circulation patterns place the dominant controls on deep sea that the LGM was not directly equivalent to the simulated radiocarbon, with only a secondary imprint of changes in the ‘glacial’ state, as discussed further below. circulation rate, quantified here by the ideal age. The air–sea exchange timescale for radiocarbon is very 4.4.4 Glacial nitrate, deep ocean ­O2 and carbon storage slow, requiring multiple years of exposure to approach satu- ration. Because deep water upwelling and subduction can As shown above, the glacial simulation agrees well with the occur on time scales shorter than this, the exposure history reconstructed temperature, salinity and water mass proxies of waters at the ocean surface plays a role in determining the for the LGM. However, it simulates high summertime N ­ O3 initial, pre-formed radiocarbon activity of sinking waters. In concentrations at the Southern Ocean surface, and an oxy- the pre-industrial ocean, NADW had high radiocarbon activ- gen-rich deep ocean (Fig. 5). In contrast, nitrogen isotopic ity due to the long exposure of its source waters in the sub- evidence has been interpreted as showing that the LGM tropical Atlantic, whereas the Southern Ocean waters were Southern Ocean surface was almost completely devoid of relatively poorly equilibrated. In the glacial circulation, the ­NO3 in summer (François et al. 1997; Martínez-García et al. input of radiocarbon-rich NADW to the AABW formation 2014; Wang et al. 2017), and the LGM deep ocean appears regions is diminished, because it is insufficiently dense to to have had significantly lower O ­ 2 (Bradtmiller et al. 2010; upwell at high latitudes of the Southern Ocean. In addition, Galbraith and Jaccard 2015; Hoogakker et al. 2015; Korff the Southern Ocean is persistently sea-ice covered, blocking et al. 2016) as expected given greater storage of respired local air–sea exchange. Both these factors restrict the supply carbon by the biological soft tissue pump. Recent evidence 13 672 E. Galbraith, C. de Lavergne even suggests that ­O2 was quite low at 3 km depth in the explain why the model does not simulate the low deep ocean LGM Southern Ocean (Jaccard et al. 2016), and that this ­O2 indicated by multiple proxy records. depletion may have extended to relatively shallow depths Second, in addition to the ventilation age bias, it is pos- (Lu et al. 2016). It seems possible that the Southern Ocean sible that the absence of slope convection in the model con- nitrogen isotopic records could have been misinterpreted if, tributes to excessive oxygen uptake in the glacial state rela- for example, there had been isotopic enrichment due to deni- tive to the interglacial state. CM2Mc can only form deep trification in the ­O2-depleted glacial deep ocean. However, waters through vertical convection, in grid cells that are the diversity of oxygen proxies makes the glacial oxygen more than 100 km across, and does not simulate ice-shelf depletion likely very robust: the model is clearly missing interactions. This mode of ventilation ‘aerates’ the whole something. Buchanan et al. (2016) found a similar increase water column, rapidly pumping oxygen into the deep sea. In in deep ocean O ­ 2 in a well-equilibrated simulation with the contrast, most of the production of AABW in the modern CSIRO coupled model under glacial boundary conditions, ocean occurs via cascading plumes of dense shelf waters suggesting one or more missing elements may be common down the continental slope (Orsi et al. 1999; Heuzé et al. in the response of existing models to LGM radiative and ice 2013), which climate models are incapable of simulating. sheet forcing. We propose three non-exclusive hypotheses This bias likely introduces an error in the simulation of the for what might cause these models to grossly overestimate modern Southern Ocean, but could greatly exaggerate the the glacial deep ocean oxygen content. oxygen input to the deep sea under glacial conditions, under First, as argued in the prior section, it appears that the which the rapid production of dense brines produces very simulated glacial deep ocean ideal age is roughly 200 years large, persistent polynyas, mostly in the Weddell Sea. The too young, when compared to the LGM. A 200-year bias in interpretation that the glacial Southern Ocean was charac- ventilation age would suggest an underestimate of ­O2 utiliza- terized by low rates of near-surface vertical mixing, based tion in the deep sea of 60–80 µM, according to a regression on sedimentary observations (Francois et al. 1997; Sigman of simulated age vs. O ­ 2 utilization (Eggleston and Galbraith et al. 2010), would certainly suggest that the model’s rapid 2018), which would bring the simulated glacial ­O2 back vertical aeration of the water column is incorrect, and may down to interglacial levels. Others have argued that ventila- indicate additional biases in the simulation of mixed layer tion of the LGM deep ocean was significantly slower than dynamics. PI, based on model simulations forced with varying fresh- Third, it is likely that significant changes in oxygen con- water inputs (Bopp et al. 2017; Menviel et al. 2017), which sumption occurred via ecosystem processes not included would clearly help to achieve low glacial ­O2 concentrations; here. One likely candidate is iron fertilization, which cer- however it is not obvious that a glacial steady-state deep tainly contributed to Southern Ocean N ­ O3 drawdown to ocean could have been substantially older than PI without some degree (Martínez-García et al. 2014), though proxy violating the observed constraints on deep ocean radiocar- evidence for low Antarctic export production precludes bon and salinity. this being the sole explanation of widespread N ­ O3 deple- It is conceivable that the LGM deep ocean diverged sig- tion (Sigman et al. 2010; Jaccard et al. 2013; Wang et al. nificantly from a steady state: the obliquity minimum, more 2017). Deepening of the remineralization profile, as would directly comparable to our ‘glacial’ simulation, occurred at have resulted from slow ecosystem metabolism under low 28 ka, and it is possible that the gradually increasing axial temperatures (Matsumoto et al. 2007), could have contrib- tilt over the ensuing millennia reduced sea ice freshwater uted to low deep ocean oxygen concentrations, though again transport, driving a freshening trend in the Southern Ocean it seems unlikely to explain the surface nitrate depletion. that truncated AABW formation during the LGM. A related Simulations previously made by Buchanan et al. (2016) and possibility is that the advance of the Antarctic ice sheet to Galbraith and Jaccard (2015) suggest that the combination the shelf margins eliminated important coastal AABW for- of increased export production and deeper remineralization mation regions (Paillard and Parrenin 2004), a mechanism (together with the 200-year age underestimate) could pos- which climate models cannot yet capture. An LGM decrease sibly close the gap between the CM2Mc glacial state and the of AABW formation could perhaps explain the otherwise- existing LGM oxygen observations. Testing these hypoth- puzzling drop of deep ocean radiocarbon that appears to eses would benefit from improved quantitative estimates of occur between 25 and 20 ka (e.g. Burke and Robinson 2012, dissolved oxygen and export production, and from future their Fig. 2b). In this case, the 200-year lower-than-PI age model simulations that include additional tracers (e.g. δ13C, of our ‘glacial’ simulation may not actually be a model bias, εNd) that can be compared directly with observations. 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 13 Response of a comprehensive climate model to a broad range of external forcings: relevance for… 673 4.5 Deep ocean circulation in a warmer world (a) Orbital sensitivity of global surface air temperature 3 Surface air temperature range (C) The discussion thus far has focused on the glacial-intergla- 2.5 cial contrast, given the relatively strong observational con- LGM ice sheets straints available. We now point out a few salient features of 2 the high ­CO2 simulations. As discussed above, minima of the AABW volume con- 1.5 PI ice sheets tribution and ventilation rate occur at approximately pre- sent-day pCO2 (Figs. 4, 5). The lowest AABW and highest 1 NADW fractions of any simulation occur at 405 ppm, with high obliquity and strong boreal seasonality. This configu- 0.5 147 180 220 270 405 607 911 ration—interglacial ­CO2, high obliquity and strong boreal CO 2 (ppmv) seasons—occurred most recently in the early Holocene (ca. 9 ka), when NADW formation may have been stronger based (b) Orbital sensitivity of planetary albedo 1.2 on both ocean circulation proxies (Hoogakker et al. 2011) and surface water characteristics (Thornalley et al. 2009; 1 Albedo range (%) Marcott et al. 2013). The NADW-strengthening orbital forc- ing configuration would have been particularly strong during 0.8 MIS 5e due to the greater eccentricity at that time (> 0.04), with a maximum effect near peak obliquity (131 ka) and 0.6 strongest boreal seasons (127 ka), implying that NADW for- 0.4 mation and volume would have been even larger then. This model-based prediction coincides remarkably well with an 0.2 apparent stagnation event recorded in Southern Ocean sedi- 147 180 220 270 405 607 911 CO 2 (ppmv) ments dated at 127 ka and interpreted as a brief cessation of AABW formation (Hayes et al. 2014). We therefore suggest (c) Orbital sensitivity of AABW fraction the possibility that this stagnation event was enabled by a 25 ‘perfect storm’ of astronomical configuration and relatively AABW fraction range (%) high ­CO2, though it may also have been pushed over the 20 edge by glacial melting as proposed by Hayes et al. (2014). As pCO2 rises above present-day levels, the AABW vol- 15 ume fraction once again increases. Although this nonlinear- ity is surprising, it would appear to agree with the observa- 10 tional inferences of dominant AABW under the high ­CO2 of the Pliocene (Hodell and Venz-Curtis 2006; Bell et al. 2015). We stress that this is an equilibrium response to warming, 5 147 180 220 270 405 607 911 and is very different from the transient response to rapid CO 2 (ppmv) (e.g. anthropogenic) warming, under which AABW forma- tion ceases due to freshwater stratification (de Lavergne et al. Fig. 16  Enhanced sensitivity to orbital forcing at low atmospheric 2014). ­CO2. a Each symbol shows the range of globally-averaged surface We attribute the return of AABW dominance at high air temperatures among the 4 orbital configurations under each ­CO2 concentration. Unfilled circles indicate glacial ice sheet topography temperature to the relatively rapid warming of North Atlan- and albedo, filled circles indicate pre-industrial ice sheet topography tic waters, which leads to AABW once again gaining the and albedo. b As in panel a, but for the range of planetary albedo, i.e. upper hand of density (Figs. 4, 6, 8), aided by the nonlinear the ratio of outgoing to incident shortwave radiation at the top of the dependence of density on temperature (De Boer et al. 2007). atmosphere. c As in panel a, but for the AABW volume fraction in the global deep ocean In the thermally-stratified deep ocean, precession becomes more important than obliquity, generating a greater con- tribution of NADW when boreal winters are particularly stratification-dependent buoyancy fluxes in the interior, the cold. As described for the low pCO2 trends, the increased larger density difference enhances the formation and con- ventilation follows from the increase in the surface den- sumption rate of abyssal waters, lowering ideal age once sity difference between Antarctic and North Atlantic deep again (Fig. 5). We note that in a warm world with no Antarc- water formation sites. By promoting convective instabilities tic ice sheet (e.g. the Eocene), much higher Southern Ocean around Antarctica (Yamamoto et al. 2015) and increasing the 13 674 E. Galbraith, C. de Lavergne temperatures may have led to a different relationship with (Martinez-Garcia et al. 2011). Thus a gradual drawdown of ­CO2 than that shown here (Lefebvre et al. 2012). ­CO2, resulting from an excess of carbon burial relative to ­CO2 supply from outgassing and weathering, would have increased the leverage of orbital forcing on global climate, 4.6 CO2‑dependent sensitivity to astronomical a key ingredient in the development of the high-amplitude forcing ice age cycles of the Pleistocene. One of the mysteries of the late Cenozoic is why the response of the climate system to astronomical forcing inten- 5 Summary sified from the Pliocene into the Pleistocene (Lisiecki and Raymo 2007). A number of explanations focus on the direct Our simulations provide a very broad view on the equi- influence of orbital forcings on the Laurentide ice sheet librium response of a comprehensive ocean–atmosphere (Maslin et al. 1998) and Antarctic ice sheet (Raymo et al. model to the three most important external forcings of the 2006). Our model simulations provide a different perspec- late Cenozoic: greenhouse gas concentrations, astronomi- tive, illustrating the degree to which an inherent increase in cal changes, and changes in ice sheet size. The simula- Earth system sensitivity to orbital forcing could have con- tions do not account for changes in vegetation or aero- tributed to this shift. sols, and hold some important parameters (e.g. the orbit Figure 16 shows that the response of the climate system eccentricity) fixed. Moreover, the coarse resolution of the to a given orbital forcing change is strongly amplified at low model, the absence of important processes (e.g. interac- ­CO2. For the global surface air temperature, this involves tive ice-shelves, downslope flows of dense waters) and a fourfold increase in the combined impact of obliquity the limitations of the parameterizations (e.g. clouds in the and precession, increasing from 0.6 to 2.5 °C as the ­CO2 atmosphere, tidal mixing and mesoscale eddy effects in the decreases from 911 to 147 ppm (Fig. 16a). The reason ocean), introduce important biases in the results. None- underlying the ­CO2-dependence of orbital sensitivity is that theless, the simulations provide a bird’s-eye view of the albedo feedbacks become stronger in a colder world, as sea expectations of a modern climate model for equilibrium ice and snow cover expand, increasing the leverage of orbital changes under a broad range of boundary conditions. As changes on global surface air temperature (Fig. 16b). The such, they can be seen as networks of internally-consistent sensitivity to obliquity is particularly responsive to C ­ O2 hypotheses that can be tested against observations. (Fig. 3, compare orange vs. blue symbols), which we attrib- We summarize the main points as follows: ute primarily to its important role as a modulator of sea ice extent. The impact of obliquity on sea ice is only significant 1. The globally-averaged surface air temperature is rela- when high latitude waters are cold enough for extensive tively high when obliquity is high and when the preces- sea ice to develop. Replacing pre-industrial with glacial ice sional phase results in weak boreal seasonality. sheets slightly increases this sensitivity, by further increas- 2. The volume contribution of AABW to the global deep ing the potential response of albedo. ocean varies nonlinearly with global temperature. The volume fraction of AABW is also more sensitive to AABW maxima occur at both low and high C ­ O2, with orbital forcing at low C­ O2, although it shows a pronounced a minimum at interglacial C ­ O2. The deep ocean vol- nonlinearity, peaking at the average late Pleistocene ­CO2 ume ventilated from the North Atlantic varies as a mir- of 220 ppm (Fig. 16c). Given that the volume fraction of ror image of this, showing a maximum at interglacial AABW can play a role in modulating atmospheric ­CO2 ­CO2. through disequilibrium carbon storage (Skinner 2009; 3. The simulated changes in deep water-mass volumes Eggleston and Galbraith 2018), it seems possible that this can be predicted from the density difference between enhanced sensitivity of AABW fraction to orbital forcing surface waters of the AABW and NADW formation at Pleistocene temperatures increased the carbon cycle regions. This appears to be consistent with the den- response to orbital forcing relative to earlier in the Cenozoic. sity-overlap mechanism discussed by Wolfe and Cessi We suggest that increased orbital sensitivity of albedo- (2009, 2010, 2011). driven changes in surface air temperature, combined with 4. Under low ­CO2, relatively dense, salty AABW is increased sensitivity of C­ O2 to orbital forcing (particularly produced. Cold atmospheric temperatures decrease obliquity) via variations in the volume fraction of AABW, precipitation in the Southern Ocean and enhance the contributed to the intensification of the ice ages through the formation and equatorward export of Antarctic sea ice, Plio–Pleistocene. This increased sensitivity may have been turning the high latitude Southern Ocean into a net further amplified through additional Earth system feed- salinifying environment, instead of the net freshening backs, such as the production of dust by Patagonian glaciers environment it is today. In a cold world, low obliquity 13 Response of a comprehensive climate model to a broad range of external forcings: relevance for… 675 further amplifies sea ice cycling and Southern Ocean the marine ecosystem in addition to the changes in salinification, while precession has a smaller direct deep water-mass geometry. Unresolved changes in the impact on the ocean. Southern Ocean near-surface environment (e.g. sea- 5. Full LGM ice sheets also increase the salinity and vol- sonal stratification and/or iron fertilization) or more ume contribution of AABW, through the downstream widely-distributed changes (e.g. changes in the vertical impact of cold NADW on Southern Ocean sea ice profile of remineralization or low latitude export) could cycling. However, the large ice sheets simultaneously play important roles. accelerate the AMOC by increasing convection in the 11. Under lower ­CO2, the model shows a greater cli- North Atlantic. mate sensitivity to orbital forcing, due to an increased 6. The combination of low C ­ O2, low obliquity and a low- response of albedo, and an enhanced sensitivity of elevation Laurentide ice sheet leads to a ‘sweet spot’ deep ocean ventilation via the control of sea ice on at which the AMOC is most unstable. We suggest that AABW salinity. We suggest that an effect of AABW the AMOC was most susceptible to freshwater per- volume on ocean carbon storage, together with the turbations under this state, contributing to enhanced enhanced albedo sensitivity, gave greater leverage to Dansgaard–Oeschger style abrupt climate change when orbital forcing as atmospheric C ­ O2 declined over the the Southern Ocean was very cold due to low C ­ O2 Cenozoic. This tendency would have contributed to and/or low obliquity but ice sheets were not at their the Plio–Pleistocene development and intensification maxima. of the ice age cycles. 7. At higher C­ O2, the influence of sea ice is diminished as its areal extent contracts. Above 405 ppm the AABW The analysis provided here has only scratched the sur- volume expands once again, which we attribute to the face of the behaviour exhibited in this large suite of model fact that the Antarctic convective regions warm more simulations. It may be valuable to undertake further analysis, slowly than their North Atlantic counterparts. such as for comparison with specific proxy sites, comparison 8. The simulated deep ocean temperature and salinity of with other models, exploration of the atmospheric dynamics, the ‘glacial’ state agree very closely with reconstructed or consideration of the changes in ocean circulation accord- deep ocean temperature and salinity for the LGM. The ing to existing theoretical frameworks (e.g. Gnanadesikan simulated expansion of AABW is also qualitatively 1999; Wolfe and Cessi 2011; Nikurashin and Vallis 2012; de consistent with common interpretations of δ13C and Lavergne et al. 2017). To this end we will make all simula- εNd. We take this consistency as an indication that tion output freely available. prior disagreements between model LGM simulations and reconstructions could have arisen from model-spe- Acknowledgements We are grateful to Compute Canada and Scinet (University of Toronto) for providing the computational resources to cific biases or insufficient equilibration, rather than a carry out the simulations, and to the Canadian Foundation for Innova- general inability for modern GCMs to reproduce the tion for local computing infrastructure. We thank Rick Slater, Simon glacial climate state—as far as deep water temperature, Yang and Nicolas Brown for assistance with configuring and running salinity and surface origin are concerned. the model, and Tim Merlis, Luke Skinner, Sam Jaccard, Sarah Egg- leston, Robbie Toggweiler, Daniel Sigman, Robert Anderson, Jorge 9. Despite the close agreement with LGM temperature Sarmiento and Jaime Palter for discussions. Pearse Buchanan and three and salinity, the radiocarbon simulation suggests that anonymous reviewers provided constructive comments on the manu- the glacial deep ocean ideal age is on the order of script. EDG acknowledges financial support from the Spanish Ministry 200 years too low, implying that it should have been of Economy and Competitiveness, through the María de Maeztu Pro- gramme for Centres/Units of Excellence in R&D (MDM-2015-0552). quite similar to the pre-industrial, despite the strong CdL acknowledges Australian Research Council support through grant sea-ice driven buoyancy fluxes in the Southern Ocean. FL150100090. All model runscripts, code and simulation output are This discrepancy may arise from an inadequate repre- freely available from the authors, or can be obtained by download from sentation of deep ocean mixing as parameterized in https​://earth​syste​mdyna​mics.org/cm2mc​-simul​ation​-libra​ry/. the model. It is also possible that the LGM was not in Open Access This article is distributed under the terms of the Crea- steady state. We highlight the unreliability of simulated tive Commons Attribution 4.0 International License (http://creat​iveco​ ventilation rates in existing climate models, given the mmons​.org/licen​ses/by/4.0/), which permits unrestricted use, distribu- poor representation of deep mixing processes. tion, and reproduction in any medium, provided you give appropriate 10. The simulated glacial deep ocean ­O2 concentrations credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 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 pro- cesses of dense waters, and/or important changes in 13 676 E. Galbraith, C. de Lavergne References 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 Adkins JF (2013) The role of deep ocean circulation in setting glacial Cessi P, Jones CS (2017) Warm-route versus cold-route interbasin climates. Paleoceanography. https​://doi.org/10.1002/palo.20046​ exchange in the meridional overturning circulation. J Phys Adkins JF, Schrag DP (2003) Reconstructing Last Glacial Maximum Oceanogr 47:1981–1997 bottom water salinities from deep-sea sediment pore fluid pro- Cheng H et al (2009) Ice age terminations. Science 326:248–252. https​ files. Earth Planet Sci Lett 216:109–123 ://doi.org/10.1126/Scien​ce.11778​40 Adkins JF, McIntyre K, Schrag DP (2002) The salinity, temperature, Crowley TJ (1992) North Atlantic deep water cools the southern hemi- and delta O-18 of the glacial deep ocean. Science 298:1769–1773 sphere. Paleoceanography 7:489–497 Anderson RF, Ali S, Bradtmiller LI, Nielsen SHH, Fleisher MQ, Curry WB, Oppo D (2005) Glacial water mass geometry and the dis- Anderson BE, Burckle LH (2009) Wind-driven upwelling in the tribution of δ13C of sigma-CO2 in the western. Atlantic Ocean southern ocean and the deglacial rise in atmospheric C ­ O2. Sci- Paleoceanography. https​://doi.org/10.1029/2004P​A0010​21 ence 323:1443–1448. https​://doi.org/10.1126/Scien​ce.11674​41 De Boer AM, Sigman DM, Toggweiler JR, Russell JL (2007) The Annan J, Hargreaves J (2013) A new global reconstruction of tempera- effect of global ocean temperature change on deep ocean ventila- ture changes at the Last Glacial Maximum. Clim Past 9:367–376 tion. Paleoceanography 22 Bard E (1988) Correction of accelerator mass spectrometry 14C ages de Lavergne C, Palter JB, Galbraith ED, Bernardello R, Marinov I measured in planktonic foraminifera: Paleoceanographic implica- (2014) Cessation of deep convection in the open Southern Ocean tions. Paleoceanography 3:635–695 under anthropogenic climate change. Nat Clim Change 4:278– Bell DB, Jung SJA, Kroon D (2015) The Plio–Pleistocene development 282. https​://doi.org/10.1038/nclim​ate21​32 of Atlantic deep-water circulation and its influence on climate de Lavergne C, Madec G, Roquet F, Holmes RM, McDougall TJ trends. Quatern Sci Rev 123:265–282. https​://doi.org/10.1016/j. (2017) Abyssal ocean overturning shaped by seafloor distribu- quasc​irev.2015.06.026 tion. Nature 551:181–186 Bereiter B, Shackleton S, Baggenstos D, Kawamura K, Severinghaus DeVries T, Primeau F (2010) An improved method for estimating J (2018) Mean global ocean temperatures duing the last glacial water-mass ventilation age from radiocarbon data. Earth Planet transition. Nature 553:39–44 Sci Lett 295:367–378 Böhm E et al (2015) Strong and deep Atlantic meridional overturning Dunne JP et al (2012) GFDL’s ESM2 global coupled climate—carbon circulation during the last glacial cycle. Nature 517:73–76 earth system models. Part I: Physical formulation and baseline Bopp L, Resplandy L, Untersee A, Le Mezo P, Kageyama M (2017) simulation characteristics. J Clim 25:6646–6665. https​://doi. Ocean (de)oxygenation from the Last Glacial Maximum to the org/10.1175/jcli-d-11-00560​.1 twenty-first century: insights from Earth System models. Philos Duplessy JC, Shackleton NJ, Matthews RK, Prell WL, Ruddiman Trans R Soc A Math Phys Eng Sci. https​://doi.org/10.1098/ WF, Caralp MH, Hendy CH (1984) C-13 record of benthic rsta.2016.0323 foraminifera in the last interglacial ocean—implications for the Bouttes N, Paillard D, Roche D (2010) Impact of brine-induced strati- carbon cycle and the global deep-water circulation. Quatern Res fication on the glacial carbon cycle. Clim Past 6:575–589 21:225–243 Boyle EA, Keigwin L (1982) Deep circulation of the North Atlan- Durack PJ, Wijffels SE, Matear RJ (2012) Ocean salinities reveal strong tic over the Last 200,000 years: geochemical evidence. Science global water cycle intensification during 1950–2000. Science 218:784–787 336:455–458 Braconnot P et al (2007) Results of PMIP2 coupled simulations of the Eggleston S, Galbraith ED (2018) The devil’s in the disequilibrium: Mid-Holocene and Last Glacial Maximum—Part 1: experiments sensitivity of ocean carbon storage to climate state and iron fer- and large-scale features. Clim Past 3:261–277 tilization in a general circulation model. Biogeosci Discuss. https​ Bradtmiller LI, Anderson RF, Sachs JP, Fleisher MQ (2010) A deeper ://doi.org/10.5194/bg-2017-328 respired carbon pool in the glacial equatorial Pacific Ocean. Emiliani C (1955) Pleistocene temperatures. J Geol 63:537–578 Earth Planet Sci Lett 299:417–425. https​://doi.org/10.1016/J. Ferrari R, Jansen MF, Adkins JF, Burke A, Stewart AL, Thompson AF Epsl.2010.09.022 (2014) Antarctic sea ice control on ocean circulation in present Broccoli A, Manabe S (1987) The influence of continental ice, atmos- and glacial climates. Proc Natl Acad Sci 201323922 pheric ­CO2, and land albedo on the climate of the last glacial Fichefet T, Hovine S, Duplessy JC (1994) A model study of the Atlan- maximum. Clim Dyn 1:87–99 tic thermohaline circulation during the last glacial maximum. Brovkin V, Ganopolski A, Archer D, Munhoven G (2012) Glacial Nature 372:252–255 ­CO2 cycle as a succession of key physical and biogeochemi- Fischer H et al (2010) The role of Southern Ocean processes in orbital cal processes. Clim Past 8:251–264. https​://doi.org/10.5194/ and millennial ­CO2 variations—a synthesis. Quatern Sci Rev Cp-8-251-2012 29:193–205 Brown N, Galbraith E (2015) Hosed vs. unhosed: global response to Fox-Kemper B, Ferrari R, Hallberg R (2008) Parameterization of interruptions of the Atlantic Meridional Overturning, with and mixed layer eddies. Part I: Theory and diagnosis. J Phys Ocean- without freshwater forcing. Clim Past Discuss 11:4669–4700 ogr 38:1145–1165 Buchanan PJ, Matear RJ, Lenton A, Phipps SJ, Chase Z, Etheridge DM François R et al (1997) Contribution of Southern Ocean surface-water (2016) The simulated climate of the Last Glacial Maximum and stratification to low atmospheric ­CO2 concentrations during the insights into the global marine carbon cycle. Clim Past 12:2271 last glacial period. Nature 389:929–935 Buizert C, Schmittner A (2015) Southern Ocean control of glacial Frölicher TL, Winton M, Sarmiento JL (2014) Continued global warm- AMOC stability and Dansgaard–Oeschger interstadial duration. ing after ­CO2 emissions stoppage. Nat Clim Change 4:40–44 Paleoceanography 30:1595–1612. https:​ //doi.org/10.1002/2015P​ Galbraith ED, Jaccard SL (2015) Deglacial weakening of the oceanic A0027​95 soft tissue pump: global constraints from sedimentary nitrogen Burke A, Robinson LF (2012) The Southern Ocean’s role in carbon isotopes and oxygenation proxies. Quatern Sci Rev 109:38–48 exchange during the last deglaciation. Science 335:557–561 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/Natur​e0622​7 13 Response of a comprehensive climate model to a broad range of external forcings: relevance for… 677 Galbraith ED, Gnanadesikan A, Dunne JP, Hiscock MR (2010) Huybers P, Wunsch C (2005) Obliquity pacing of the late Pleistocene Regional impacts of iron-light colimitation in a global biogeo- glacial terminations. Nature 434:491–494 chemical model. Biogeosciences 7:1043–1064 Insua TL, Spivack AJ, Graham D, D’Hondt S, Moran K (2014) Recon- Galbraith ED et al (2011) Climate variability and radiocarbon in the struction of Pacific Ocean bottom water salinity during the Last CM2Mc earth system model. J Clim 24:4230–4254. https​://doi. Glacial Maximum. Geophys Res Lett 41:2914–2920 org/10.1175/2011j​cli39​19.1 Jaccard SL, Hayes CT, Martinez-Garcia A, Hodell DA, Anderson Galbraith ED, Kwon EY, Bianchi D, Hain MP, Sarmiento JL (2015) RF, Sigman DM, Haug GH (2013) Two modes of change in The impact of atmospheric pCO2 on carbon isotope ratios of the southern ocean productivity over the past million years. Science atmosphere and ocean. Global Biogeochem Cycles 29:307–324 339:1419–1423. https​://doi.org/10.1126/Scien​ce.12275​45 Galbraith ED, Merlis TM, Palter JB (2016) Destabilization of glacial Jaccard SL, Galbraith ED, Martínez-García A, Anderson RF (2016) climate by the radiative impact of Atlantic Meridional Overturn- Covariation of deep Southern Ocean oxygenation and atmos- ing Circulation disruptions. Geophys Res Lett 43:8214–8221 pheric ­CO2 through the last ice age. Nature 530:207–210 Ganopolski A, Winkelmann R, Schellnhuber H (2016) Critical insola- Jansen MF (2017) Glacial ocean circulation and stratification explained tion—CO2 relation for diagnosing past and future glacial incep- by reduced atmospheric temperature. Proc Natl Acad Sci tion. Nature 529:200–203 114:45–50 Gebbie G (2014) How much did Glacial North Atlantic Water shoal? Jansen MF, Nadeau L-P (2016) The effect of southern ocean surface Paleoceanography 29:190–209 Buoyancy loss on the deep-ocean circulation and stratification. J Gent PR, McWilliams JC (1990) Isopycnal mixing in ocean circulation Phys Oceanogr 46:3455–3470 models. J Phys Oceanogr 20:150–155 Kawamura K et al (2017) State dependence of climatic instability over Gildor H, Tziperman E, Toggweiler JR (2002) Sea ice switch mecha- the past 720,000 years from Antarctic ice cores and climate mod- nism and glacial-interglacial ­CO2 variations. Global Biogeochem eling. Sci Adv 3:e1600446 Cycles 16:art.no.-1032 Keeling RF, Stephens BB (2001) Antarctic sea ice and the control of Gnanadesikan A (1999) A simple predictive model for the structure of Pleistocene climate instability. Paleoceanography 16:112–131 the oceanic Pycnocline. Science 283:2077–2079 Keigwin LD, Swift SA (2017) Carbon isotope evidence for a northern Gong X, Zhang X, Lohmann G, Wei W, Zhang X, Pfeiffer M (2015) source of deep water in the glacial western North Atlantic. Proc Higher Laurentide and Greenland ice sheets strengthen the North Natl Acad Sci 114:2831–2835 Atlantic ocean circulation. Clim Dyn 45:139–150 Kienast SS, Kienast M, Mix AC, Calvert SE, François R (2007) Tho- Goosse H, Fichefet T (1999) Importance of ice–ocean interactions rium-230 normalized particle flux and sediment focusing in the for the global ocean circulation: a model study. J Geophys Res Panama Basin region during the last 30,000 years. Paleoceanog- Oceans 104:23337–23355 raphy. https​://doi.org/10.1029/2006P​A0013​57 Griffies SM (1998) The Gent–McWilliams skew flux. J Phys Oceanogr Klockmann M, Mikolajewicz U, Marotzke J (2016) The effect of green- 28:831–841 house gas concentrations and ice sheets on the glacial AMOC in Hain MP, Sigman DM, Haug GH (2010) Carbon dioxide effects of a coupled climate model. Clim Past 12:1829–1846 Antarctic stratification, North Atlantic Intermediate Water forma- Köhler P, Fischer H, Munhoven G, Zeebe RE (2005) Quantitative tion, and subantarctic nutrient drawdown during the last ice age: interpretation of atmospheric carbon records over the last gla- diagnosis and synthesis in a geochemical box model. Global Bio- cial termination. Global Biogeochem Cycles 19 https​://doi. geochem Cycles 24 doi:https​://doi.org/10.1029/2010g​b0037​90 org/10.1029/2004g​b0023​45 Hall IR, McCave IN, Shackleton NJ, Weedon GP, Harris SE (2001) Korff L, Dobeneck T, Frederichs T, Kasten S, Kuhn G, Gersonde R, Intensified deep Pacific inflow and ventilation in Pleistocene Diekmann B (2016) Cyclic magnetite dissolution in Pleistocene glacial times. Nature 412:809–812 sediments of the abyssal northwest Pacific Ocean: evidence for Hayes CT et al (2014) A stagnation event in the deep South Atlantic glacial oxygen depletion and carbon trapping. Paleoceanography during the last interglacial period. Science 346:1514–1517 Krueger S, Leuschner D, Ehrmann W, Schmiedl G, Mackensen A Held I (1982) Climate models and the astronomical theory of the ice (2012) North Atlantic Deep Water and Antarctic Bottom Water Ages. Icarus 50:449–461 variability during the last 200 ka recorded in an abyssal sediment Held IM, Soden BJ (2006) Robust responses of the hydrological core off South Africa. Global Planet Change 80:180–189 cycle to global warming. J Clim 19:5686–5699 Kutzbach J, Liu X, Liu Z, Chen G (2008) Simulation of the evolution- Heuzé C, Heywood KJ, Stevens DP, Ridley JK (2013) Southern ary response of global summer monsoons to orbital forcing over Ocean bottom water characteristics in CMIP5 models. Geo- the past 280,000 years. Clim Dyn 30:567–579 phys Res Lett 40:1409–1414 Kutzbach JE, He F, Vavrus SJ, Ruddiman WF (2013) The dependence Hewitt CD, Stouffer RJ, Broccoli AJ, Mitchell JFB, Valdes PJ (2003) of equilibrium climate sensitivity on climate state: Applications The effect of ocean dynamics in a coupled GCM simulation of to studies of climates colder than present. Geophys Res Lett the Last Glacial Maximum. Clim Dyn 20:203–218 40:3721–3726 Ho SL, Laepple T (2015) Glacial cooling as inferred from marine Kwon EY, Hain MP, Sigman DM, Galbraith ED, Sarmiento JL, Tog- temperature proxies TEX86 and UK′ 37. Earth Planet Sci Lett gweiler JR (2012) North Atlantic ventilation of “southern- 409:15–22 sourced” deep water in the glacial ocean. Paleoceanography 27 Hodell DA, Venz-Curtis KA (2006) Late Neogene history of deepwater Large WG, McWilliams JC, Doney SC (1994) Oceanic vertical ventilation in the Southern Ocean. Geochem Geophys Geosyst 7 mixing—a review and a model with a nonlocal boundary-layer Hoogakker BA, Chapman MR, McCave IN, Hillaire-Marcel C, Ellison parameterization. Rev Geophys 32:363–403 CR, Hall IR, Telford RJ (2011) Dynamics of North Atlantic deep Laskar J, Robutel P, Joutel F, Gastineau M, Correia A, Levrard B water masses during the Holocene. Paleoceanography 26 (2004) A long-term numerical solution for the insolation quan- Hoogakker BA, Elderfield H, Schmiedl G, McCave IN, Rickaby RE tities of the Earth. Astron Astrophys 428:261–285 (2015) Glacial-interglacial changes in bottom-water oxygen con- Lee JE, Shen A, Fox-Kemper B, Ming Y (2017) Hemispheric sea tent on the Portuguese margin. Nat Geosci 8:40–43 ice distribution sets the glacial tempo. Geophys Res Lett Howe JN, Piotrowski AM, Noble TL, Mulitza S, Chiessi CM, Bayon 44:1008–1014 G (2016) North Atlantic deep water production during the last Lefebvre V, Donnadieu Y, Sepulchre P, Swingedouw D, Zhang ZS glacial maximum. Nat Commun 7 (2012) Deciphering the role of southern gateways and carbon 13 678 E. Galbraith, C. de Lavergne dioxide on the onset of the Antarctic Circumpolar Current. Pale- Otto-Bliesner B et al. (2007) Last Glacial Maximum ocean thermoha- oceanography 27 line circulation: PMIP2 model intercomparisons and data con- Li C, Battisti DS (2008) Reduced Atlantic storminess during Last Gla- straints. Geophys Res Lett 34 cial Maximum: evidence from a coupled climate model. J Clim Paillard D, Parrenin F (2004) The Antarctic ice sheet and the triggering 21:3561–3579 of deglaciations. Earth Planet Sci Lett 227:263–271 Lippold J et al (2016) Deep water provenance and dynamics of the (de) Pausata FS, Li C, Wettstein J, Kageyama M, Nisancioglu KH (2011) glacial Atlantic meridional overturning circulation. Earth Planet The key role of topography in altering North Atlantic atmos- Sci Lett 445:68–78 pheric circulation during the last glacial period. Clim Past Lisiecki LE, Raymo ME (2007) Plio–Pleistocene climate evolution: 7:1089–1101 trends and transitions in glacial cycle dynamics. Quatern Sci Pena LD, Goldstein SL (2014) Thermohaline circulation crisis Rev 26:56–69. https​://doi.org/10.1016/j.quasc​irev.2006.09.005 and impacts during the mid-Pleistocene transition. Science Lu Z et al. (2016) Oxygen depletion recorded in upper waters of the 345:318–322 glacial Southern Ocean. Nat Commun 7 Piotrowski AM, Goldstein SL, Hemming SR, Fairbanks RG (2005) MacKinnon JA et al (2017) Climate process team on internal-wave Temporal relationships of carbon cycling and ocean circulation driven ocean mixing. Bull Am Meteorol Soc at glacial boundaries. Science 307:1933–1938 Manabe S, Bryan K (1985) C ­ O2-induced change in a coupled ocean– Rahmstorf S (2002) Ocean circulation and climate during the past atmosphere model and its paleoclimatic implications. J Geophys 120,000 years. Nature 419:207–214 Res Oceans 90:11689–11707 Raymo M, Lisiecki L, Nisancioglu KH (2006) Plio–Pleistocene ice Marcott SA, Shakun JD, Clark PU, Mix AC (2013) A reconstruction volume, Antarctic climate, and the global δ18O record. Science of regional and global temperature for the past 11,300 years. 313:492–495 Science 339:1198–1201 Ritz SP, Stocker TF, Grimalt JO, Menviel L, Timmermann A (2013) Martinez-Garcia A, Rosell-Mele A, Jaccard S, Geibert W, Sigman DM, Estimated strength of the Atlantic overturning circulation dur- Haug GH (2011) Southern Ocean dust-climate coupling over the ing the last deglaciation. Nat Geosci 6:208–212. https​://doi. past four million years. Nature 476 org/10.1038/Ngeo1​723 Martínez-García A et al (2014) Iron fertilization of the Subantarctic Roe GH, Lindzen RS (2001) The mutual interaction between conti- Ocean during the last ice age. Science 343:1347–1350 nental-scale ice sheets and atmospheric stationary waves. J Clim Martinson DG (1990) Evolution of the Southern Ocean winter mixed 14:1450–1465 layer and sea ice: open-ocean deep water formation and ventila- Russell JL, Stouffer RJ, Dixon KW (2006) Intercomparison of the tion. J Geophys Res 95:114614–11654 Southern Ocean circulations in IPCC coupled model control Marzocchi A, Jansen MF (2017) Connecting Antarctic sea ice to deep- simulations. J Clim 19:4560–4575 ocean circulation in modern and glacial climate simulations. Sarmiento JL, Toggweiler JR (1984) A new model for the role of the Geophys Res Lett 44:6286–6295 oceans in determining atmopheric pCO2. Nature 308:621–624 Maslin MA, Li XS, Loutre M-F, A. B (1998) The contribution of Sarnthein M, Winn K, Jung SJ, Duplessy JC, Labeyrie L, Erlenkeuser orbital forcing to the progressive intensification of Northern H, Ganssen G (1994) Changes in east Atlantic deepwater circula- Hemisphere Glaciation. Quatern Sci Rev 17:411–426 tion over the last 30,000 years: eight time slice reconstructions. Matsumoto K, Hashioka T, Yamanaka Y (2007) Effect of temperature- Paleoceanography 9:209–267 dependent organic carbon decay on atmospheric pCO2. J Geo- Sarnthein M, Schneider B, Grootes PM (2013) Peak glacial 14C venti- phys Res Biogeosci 112 lation ages suggest major draw-down of carbon into the abyssal McManus JF, Oppo DW, Cullen JL (1999) A 0.5-million-year record ocean. Clim Past of millennial-scale climate variability in the North Atlantic. Sci- Schmittner A (2003) Southern Ocean sea ice and radiocarbon ages of ence 283:971–975 glacial bottom waters. Earth Planet Sci Lett 213:53–62 McManus JF, Francois R, Gherardi JM, Keigwin LD, Brown-Leger S Schmittner A, Clement AC (2002) Sensitivity of the thermohaline (2004) Collapse and rapid resumption of Atlantic meridional cir- circulation to tropical and high latitude freshwater forcing dur- culation linked to deglacial climate changes. Nature 428:834–837 ing the last glacial–interglacial cycle. Paleoceanography Menviel L, Yu J, Joos F, Mouchet A, Meissner K, England M (2017) Schmittner A, Green J, Wilmes SB (2015) Glacial ocean overturn- Poorly ventilated deep ocean at the Last Glacial Maximum ing intensified by tidal mixing in a global circulation model. inferred from carbon isotopes: a data-model comparison study. Geophys Res Lett 42:4014–4022 Paleoceanography 32:2–17 Schrag D, De Paolo D (1993) Determination of δ18O of seawater in Mix AC, Bard E, Schneider R (2001) Environmental processes of the deep ocean during the Last Glacial Maximum. Paleocean- the ice age: land, oceans, glaciers (EPILOG). Quatern Sci Rev ography 8:1–6 20:627–657 Sherriff-Tadano S, Abe-Ouchi A, Yoshimori M, Oka A, Chan W-L Muglia J, Schmittner A (2015) Glacial Atlantic overturning (2017) Influence of glacial ice sheets on the Atlantic meridi- increased by wind stress in climate models. Geophys Res Lett onal overturning circulation through surface wind change. 42:9862–9868 Clim Dyn 1–23 Munday DR, Johnson HL, Marshall DP (2013) Eddy saturation of Shin SI, Liu Z, Otto-Bliesner BL, Kutzbach JE, Vavrus SJ (2003) equilibrated circumpolar currents. J Phys Oceanogr 43:507–532 Southern Ocean sea-ice control of the glacial North Atlantic Nikurashin M, Vallis G (2012) A theory of the interhemispheric merid- thermohaline circulation. Geophys Res Lett 30 ional overturning circulation and associated stratification. J Phys Sigman DM, Boyle EA (2000) Glacial/interglacial variations in Oceanogr 42:1652–1667 atmospheric carbon dioxide. Nature 407:859–869 Nilsson J, Langen PL, Ferreira D, Marshall J (2013) Ocean basin Sigman DM, Jaccard SL, Haug GH (2004) Polar ocean stratification geometry and the salinification of the Atlantic Ocean. J Clim in a cold climate. Nature 428:59–63 26:6163–6184 Sigman DM, De Boer AM, Haug GH (2007) Antarctic stratification, Orsi AH, Johnson GC, Bullister JL (1999) Circulation, mixing, and atmospheric water vapor, and Heinrich events: a hypothesis for production of Antarctic Bottom Water. Prog Oceanogr 43:55–109 late Pleistocene deglaciations. Mechanisms and impacts-past and future changes of meridional overturning, ocean circula- tion, pp 335–349 13 Response of a comprehensive climate model to a broad range of external forcings: relevance for… 679 Sigman DM, Hain MP, Haug GH (2010) The polar ocean and glacial Toggweiler JR (1999) Variation of atmospheric ­CO2 by ventilation of cycles in atmospheric C ­ O2 concentration. Nature 466:47–55. the ocean’s deepest water. Paleoceanography 14:571–588 https​://doi.org/10.1038/Natur​e0914​9 Toggweiler JR, Russell JL, Carson SR (2006) Midlatitude westerlies, Simmons HL, Jayne SR, St Laurent LC, Weaver AJ (2004) Tidally atmospheric ­CO2, and climate change during the ice ages. Pale- driven mixing in a numerical model of the ocean general cir- oceanography 21 culation. Ocean Model 6:245–263 Wang XT et al (2017) Deep-sea coral evidence for lower Southern Simmons C, Matthews H, Mysak L (2016) Deglacial climate, carbon Ocean surface nitrate concentrations during the last ice age. Proc cycle and ocean chemistry changes in response to a terrestrial Natl Acad Sci 201615718 carbon release. Clim Dyn 46:1287–1299 Watson AJ, Vallis GK, Nikurashin M (2015) Southern Ocean buoy- Skinner L (2009) Glacial-interglacial atmospheric C­ O2 change: a possi- ancy forcing of ocean ventilation and glacial atmospheric C ­ O 2. ble” standing volume” effect on deep-ocean carbon sequestration. Nat Geosci Clim Past 5:537–550 Winton M (1997) The effect of cold climate upon North Atlantic Skinner LC, Fallon S, Waelbroeck C, Michel E, Barker S (2010) Venti- deep water formation in a simple ocean–atmosphere model. J lation of the Deep Southern Ocean and Deglacial ­CO2 Rise. Sci- Clim:37–51 ence 328:1147–1151. https​://doi.org/10.1126/Scien​ce.11836​27 Wolfe CL, Cessi P (2009) Overturning circulation in an eddy-resolving Skinner L, McCave I, Carter L, Fallon S, Scrivner A, Primeau F (2015) model: the effect of the pole-to-pole temperature gradient. J Phys Reduced ventilation and enhanced magnitude of the deep Pacific Oceanogr 39:125–142 carbon pool during the last glacial period. Earth Planet Sci Lett Wolfe CL, Cessi P (2010) What sets the strength of the middepth strati- 411:45–52 fication and overturning circulation in eddying ocean models? J Stephens BB, Keeling RF (2000) The influence of Antarctic sea ice on Phys Oceanogr 40:1520–1538 glacial-interglacial ­CO2 variations. Nature 404:171–174 Wolfe CL, Cessi P (2011) The adiabatic pole-to-pole overturning cir- Stocker TF, Wright DG, Broecker WS (1992) The influence of high- culation. J Phys Oceanogr 41:1795–1810 latitude surface forcing on the global thermohaline circulation. Wunsch C (2003) Determining paleoceanographic circulations, with Paleoceanography 7:529–541 emphasis on the Last Glacial Maximum. Quatern Sci Rev Stommel H (1961) Thermohaline convection with two stable regimes 22:371–385 of flow. Tellus 13:224–230 Wunsch C (2016) Pore fluids and the LGM ocean salinity—Reconsid- Stouffer RJ, Manabe S (2003) Equilibrium response of thermohaline ered. Quatern Sci Rev 135:154–170 circulation to large changes in atmospheric C ­ O2 concentration. Yamamoto A, Abe-Ouchi A, Shigemitsu M, Oka A, Takahashi K, Clim Dyn 20:759–773 Ohgaito R, Yamanaka Y (2015) Global deep ocean oxygenation Swingedouw D, Braconnot P, Delécluse P, Guilyardi E, Marti O (2007) by enhanced ventilation in the Southern Ocean under long-term Quantifying the AMOC feedbacks during a 2 × ­CO2 stabilization global warming. Global Biogeochem Cycles 29:1801–1815 experiment with land-ice melting. Clim Dyn 29:521–534 Zhang X, Lohmann G, Knorr G, Purcell C (2014) Abrupt glacial cli- Telford R, Li C, Kucera M (2013) Mismatch between the depth habi- mate shifts controlled by ice sheet changes. Nature 512:290–294 tat of planktonic foraminifera and the calibration depth of SST Zhu J, Liu Z, Zhang X, Eisenman I, Liu W (2014) Linear weakening of transfer functions may bias reconstructions. Clim Past 9:859–870 the AMOC in response to receding glacial ice sheets in CCSM3. Thornalley DJ, Elderfield H, McCave IN (2009) Holocene oscillations Geophys Res Lett 41:6252–6258. https:​ //doi.org/10.1002/2014G​ in temperature and salinity of the surface subpolar North Atlan- L0608​91 tic. Nature 457:711 Affiliations Eric Galbraith1,2,3 · Casimir de Lavergne4 1 3 ICREA, Pg. Lluís Companys 23, 08010 Barcelona, Spain Department of Earth and Planetary Sciences, McGill 2 University, Montreal, Canada Department of Mathematics, Institut de Ciència i Tecnologia 4 Ambientals (ICTA), Universitat Autonoma de Barcelona, School of Mathematics and Statistics, University of New 08193 Barcelona, Spain South Wales, Sydney, NSW 2052, Australia 13