Collapse of the North American ice saddle 14,500 years ago caused widespread cooling and reduced ocean overturning circulation

Collapse of ice sheets can cause significant sea level rise and widespread climate change. We examine the climatic response to meltwater generated by the collapse of the Cordilleran‐Laurentide ice saddle (North America) ~14.5 thousand years ago (ka) using a high‐resolution drainage model coupled to an ocean‐atmosphere‐vegetation general circulation model. Equivalent to 7.26 m global mean sea level rise in 340 years, the meltwater caused a 6 sverdrup weakening of Atlantic Meridional Overturning Circulation (AMOC) and widespread Northern Hemisphere cooling of 1–5°C. The greatest cooling is in the Atlantic sector high latitudes during Boreal winter (by 5–10°C), but there is also strong summer warming of 1–3°C over eastern North America. Following recent suggestions that the saddle collapse was triggered by the Bølling warming event at ~14.7–14.5 ka, we conclude that this robust submillennial mechanism may have initiated the end of the warming and/or the Older Dryas cooling through a forced AMOC weakening.


Introduction
During the Last Glacial Maximum 26-19 thousand years ago (ka), a vast ice sheet stretched over North America [Clark et al., 2009]. In subsequent millennia, as climate warmed and this ice sheet decayed, large volumes of meltwater flooded to the oceans [Tarasov and Peltier, 2006;Wickert, 2016]. This period, known as the "last deglaciation," included episodes of abrupt climate change, such as the Bølling warming, when Northern Hemisphere temperatures increased by 4-5°C in just a few decades [Lea et al., 2003;Buizert et al., 2014], coinciding with a 12-22 m sea level rise in less than 340 years (Meltwater Pulse 1a (MWP1a)) [Deschamps et al., 2012].
Currently, we do not know how (or even if) MWP1a, the Bølling warming, and a possible concurrent rapid strengthening of the Atlantic Meridional Overturning Circulation (AMOC) [e.g., McManus et al., 2004;Roberts et al., 2010] were related. Knowing the source(s) of MWP1a is an important step for resolving these open questions because the influence of freshwater on ocean circulation and surface climate has been shown to depend on where it enters the oceans [e.g., Weaver et al., 2003;Roche et al., 2009;Smith and Gregory, 2009;Menviel et al., 2011;Condron and Winsor, 2012]. For example, it has been suggested that an Antarctic source of MWP1a could have caused the Bølling warming [Weaver et al., 2003], though possibly only if the total flux is <0.2 sverdrup (Sv) [Swingedouw et al., 2008], while a Northern Hemisphere source would lead to Northern Hemisphere cooling, potentially associated with the Older Dryas cold event at~14 ka [e.g., Stanford et al., 2006;Menviel et al., 2011].
From dynamic ice sheet model results forced with a transient simulation of surface climate, Gregoire et al. [2012;henceforth G12] suggested that the North American ice sheet contributed~7 m of sea level rise in 350 years at the time of MWP1a, mostly due to the Cordilleran and Laurentide ice sheets' "saddle collapse" during their separation (Figure 1a). Recent sea level fingerprinting also supports this mechanism as a major, possibly dominant source of MWP1a [Gomez et al., 2015;Liu et al., 2016], although exclusively Antarctic provenance is not ruled out. Thus, the Cordilleran-Laurentide ice saddle collapse is a plausible source of meltwater to the deglacial ocean at~14.5 ka.
Previous modeling of ice sheet meltwater's influence on ocean circulation and climate has often employed highly idealized meltwater fluxes (in terms of timing, location, and their evolution through time) that may have been unrealistically large [e.g., Ganopolski Roberts et al. [2014]) or that attempted to constrain meltwater source, magnitude, and duration by their climatic fingerprint [e.g., Weaver et al., 2003;Liu et al., 2009;Menviel et al., 2011]. Here, for the first time, we evaluate the impact of an increased ice meltwater flux to the ocean from the Cordilleran-Laurentide saddle collapse. Also, for the first time, we connect global ice sheet mass balance to the ocean by using a drainage network model that produces fully distributed meltwater inputs. To do this, we combine results from the ice sheet modeling of G12 with the high-resolution drainage model of Wickert [2016] to determine how the ice melt is routed and where it reaches the ocean. For the remaining ice sheet meltwater fluxes (including Eurasia and Antarctica), the same drainage model is driven by the ICE-6G_C (VM5a) reconstruction [Argus et al., 2014;Peltier et al., 2015].    [Gregoire et al., 2016;henceforth G16]. For the purpose of this study, we accept the G12 hypothesis and shift the saddle collapse 2900 years backward in time such that the peak in North American ice melt corresponds to the timing of MWP1a at~14.5 ka (see section 4 for discussion on timing). As well as the event itself, the simulation includes the more gradual, ongoing "background" melt of the North American and Greenland ice sheets.
There remains considerable uncertainty on the source of MWP1a. A major contribution from the North American ice sheet is supported by a large body of evidence [Keigwin et al., 1991;Marshall and Clarke, 1999;Peltier, 2005;Tarasov and Peltier, 2005;G12;Tarasov et al., 2012]. While this scenario is compatible with fingerprinting of sea level change, additional contributions from other ice sheets are needed to explain the observations [Gomez et al., 2015;Liu et al., 2016]. Although Antarctica is not thought to have made a major contribution to sea level rise until much later (around 12 ka onward [Whitehouse et al., 2012;Gomez et al., 2013;Argus et al., 2014;Briggs et al., 2014;Mackintosh et al., 2014]), several lines of reasoning and evidence suggest that it was at least a partial source of MWP1a and could have contributed 2 m of sea level rise in 340 years [Golledge et al., 2014;Weber et al., 2014].
Here, longer-term melt histories from Antarctica and Eurasia are incorporated in the meltwater flux routed to the ocean GCM, producing 0.41 and 1.36 m global eustatic sea level rise (respectively) in 340 years (centered at 14.5 ka; the peak in global meltwater; Figure 1b). These fluxes are calculated from the recent ICE-6G_C (VM5a) reconstruction [Argus et al., 2014;Peltier et al., 2015] with interpolation between the 500 year timesteps. Between 15.5 and 13.5 ka, Eurasia and Antarctica provide a maximum discharge of less than 0.05 Sv and 0.02 Sv, respectively, remaining well below North America's discharge (always >0.1 Sv) and contributing to the background flux. Thus, it is possible that additional melt from Antarctica and Eurasia, not included in this study, could have had further influence [e.g., Golledge et al., 2014]. This work focusses on the climatic impacts of melt from the North American ice sheet saddle collapse as a necessary first step toward disentangling the events that took place at~14.5 ka.

Calculating Meltwater Drainage
Time-varying ice mass balance calculated from G12 (À2900 years) for North America and Greenland (100 year time step), and ICE-6G_C (500 year time step) for Eurasia [Peltier et al., 2015] and Antarctica [Argus et al., 2014] was dynamically routed to the ocean on a 30 arcsec (<1 km) grid as a function of surface topography, ice sheet thickness, and Glacial-Isostatic Adjustment (GIA). The computed transient meltwater flux was output as a distributed grid of coastal flow to the oceans (m 3 s À1 ) at 100 year time steps. In this way, the saddle collapse event and the background melt of ice into the oceans were combined into a global ice meltwater flux that retains its geographic distribution. Major drainage basins and outlets at 14.5 ka are indicated in Figure 1a.
Pairing the G12 ice sheet with the VM2 solid Earth model in the drainage network routing, Wickert [2016] produced excessive GIA that formed a marine embayment in the vicinity of the Laurentian Great Lakes, resulting in midcontinental river mouths for ice melting around and to the west of Lake Superior. Much of the Great Lakes corridor was ice-free in G12 because its southern Laurentide margin closely tracked the geologic record of deglaciation, meaning that a À2900 year shift in G12's ice evolution produced ice retreat past the St. Lawrence lowlands before 14.5 ka. In reality, the southern Laurentide ice melt was drained by the Mississippi and Hudson Rivers at this time ( Figure 1a), with most of the flow routed toward the Mississippi River until~12.9 ka when St. Lawrence drainage began [Williams et al., 2012;Wickert et al., 2013;Wickert, 2016]. To evaluate the effect of this uncertainty in North American drainage during MWP1a, we created two meltwater scenarios. In SC_south, we route all meltwater from Lake Superior down the Mississippi River. In SC_east, 40% of this meltwater is diverted to the Hudson River. The shape of other North American drainage basins are more robust as they are less sensitive to ice sheet geometry at this time.

The Climate Model
The  [Cox, 2001]. The ocean GCM [Gordon et al., 2000] has a horizontal resolution of 1.25°× 1.25°and 20 vertical layers. It has a fixed lid; the volume of ocean grid boxes cannot vary. Therefore, hydrological fluxes (such as evaporation, precipitation, river runoff, and ice melt) are represented as virtual salinity fluxes. It is possible that steric behavior would otherwise influence the propagation pathways of freshwater input to the ocean. However, such effects are likely to be small for pulses 0.1 Sv or less [Yin et al., 2010], and on the time scale studied here, it is likely that the large-scale enddestination of such "real" surface fluxes would be similar to those of the virtual fluxes.
The land-surface has grid-defined river catchments that instantaneously deliver meteoric water to the coast. These catchments remain constant throughout our experiments, while the routing of ice sheet meltwater is independently varied (section 2.2) as a separate input to the model.

Climate Model Experiment Design
The control simulation is based on the 15 ka simulation from Singarayer et al. [2011]. Orbital parameters follow Berger and Loutre [1991]. For the atmospheric trace gases, CO 2 is 225 ppmv, CH 4 is 473 ppbv, and N 2 O is 241 ppbv [Petit et al., 1999;Spahni et al., 2005;Parrenin et al., 2007;Loulergue et al., 2008]. The ice sheets, bathymetry, coastlines, and topography match the ICE-5G (VM2) reconstruction [Peltier, 2004], using an anomaly method for consistency with existing preindustrial boundary conditions [Singarayer and Valdes, 2010]. This is an earlier ice sheet reconstruction from the same family as ICE-6G_C (VM5a) and was the most up-to-date version when the climate simulation was initially carried out [Singarayer et al., 2011]. All of these boundary conditions were kept constant throughout the runs, so that the simulations were not complicated by their effect on climate. No ice meltwater was included in this control simulation.
With this setup, the control simulation was run for 1000 years with prescribed (preindustrial) vegetation, followed by a further 750 years with the TRIFFID dynamic vegetation model, providing a total spin-up of 1750 years. The result is a climate, land-surface, and ocean in near steady state, from which the meltwater simulations (described below) were initialized. The control simulation was continued for a further 2000 years in parallel with the meltwater simulations, confirming that the model was in an equilibrium state for the duration of the experiment.
We ran two 2000 year experiments (15.5-13.5 ka) to capture the saddle collapse event, driving the climate model with the two transient meltwater scenarios described in section 2.2: SC_south and SC_east. To isolate the effect of the century-scale saddle collapse meltwater pulse from the multimillennial scale background melt of the ice sheets, we ran two further simulations: NoSC_south and NoSC_east, each with no saddle collapse meltwater pulse. In these simulations, meltwater discharge from the SC_south and SC_east scenarios (respectively) was held constant at their 14.8 ka values for the remainder of the simulation.
For each scenario (Figures 1b-1e), the time-varying meltwater inputs were transformed onto the climate model grid and spread over adjacent ocean grid boxes with depths >500 m, along the coastline and off the continental shelf. Distributing the water prevents cells from reaching negative salinity, which would be unphysical, and keeps them within the valid range of the equation of state [Fofonoff, 1962;Bryan and Cox, 1972;Fofonoff and Millard, 1983]. Hyperpycnal meltwater flows were not simulated because these occur rarely, and typically in small, steep catchments that contribute an insignificant amount of water to the global ocean (see discussion in Wickert et al. [2013]).

Sea Surface Salinity and Ocean Circulation
In SC_south,~50% of the saddle collapse meltwater is routed down the Mackenzie River into the Arctic Ocean, with just under 50% routed down the Mississippi River into the Gulf of Mexico and Atlantic (Figure 1b). In SC_east, the Mississippi-bound meltwater is reduced to~30% of the total saddle collapse pulse, with~20% of the flux now heading to the North American East Coast (mostly down the St. Lawrence and Hudson Rivers; Figure 1d).

10.1002/2016GL071849
On the centennial time scale investigated here, there is no discernible effect on sea surface salinity from diverting 40% of the Mississippi discharge in SC_south eastward in SC_east. This is because freshwater that enters the western North Atlantic is rapidly mixed and dispersed by the subtropical (and to a lesser extent, the subpolar) gyre(s). Consequently, there is also little difference in large-scale ocean circulation patterns, such as the AMOC (e.g., Figure 2b), suggesting that there is low sensitivity to precisely where meltwater is delivered to the Atlantic subtropical gyre region. In deference to the geologic record, the remaining presentation of results and discussion will focus on SC_south (with NoSC_south used for reference), bearing in mind that the results for SC_east (and NoSC_east) are equivalent.
Meltwater runoff from the saddle collapse caps the Arctic Ocean, making it 1-6 practical salinity units (psu) fresher in the upper 200 m than in NoSC_south and less than 0.5 psu fresher below that (100 year mean, 14.4 ka). Minimum mean Arctic and Greenland-Iceland-Norwegian (GIN) sea surface salinity (Figure 2a) occurs around the same time as maximum Mackenzie River meltwater discharge is reached (14.4 ka; Figure 1b), 100 years after the peak saddle collapse event.
The relatively fresh water propagates across the Arctic, through Fram Strait and into the North Atlantic ( Figure 1a). Here it is joined by~0.07 Sv of melt from Hudson Strait, the North American East Coast, and the Mississippi River (Figure 1b). Due to vigorous mixing in the gyres, the North Atlantic has a weaker maximum freshwater anomaly and is less stratified than in the Arctic, especially in the subtropical gyre, where freshening of 1-3 psu (with respect to NoSC_south) extends to 400-500 m deep.
Broadly, this North Atlantic-Arctic freshwater cap increases the vertical density gradient, thus acting to stabilize the water column in mid-high latitudes. The change in ocean buoyancy reduces North Atlantic Deep Water (NADW) formation. By 14.4 ka, AMOC is 6 Sv weaker (À40%) and 200-400 m shallower than NoSC_south, coincident with minimum Arctic salinity and maximum meltwater discharge to the Arctic, after which it begins to recover (Figure 2b). While this significant weakening is a robust response to such freshwater forcing, the exact magnitude of the change may be model-specific, depending on the initial AMOC strength and depth, ocean structure, and sites of deepwater formation. After 14.1 ka, AMOC becomes stronger than when it started at 15.0 ka. There are signs that this strong AMOC recovery leads to an "overshoot" effect [Liu et al., 2009] at 14.1-14.0 ka, followed by a more gradual increase caused by the continued reduction in Arctic-, Atlantic-, and Pacific-bound meltwater (Figure 1b). That a more stratified high-latitude freshwater cap is produced by Arctic-draining meltwater than by Atlantic/Gulf of Mexico sources is consistent with previous findings [e.g., Condron and Winsor, 2012, who ran shorter, higherresolution simulations].
It has previously been demonstrated that a strong reduction in AMOC can lead to the development of North Pacific overturning [e.g., Saenko et al., 2004;Okazaki et al., 2010;Chikamoto et al., 2012]. Here the effect on North Pacific overturning is negligible, likely due to a smaller freshwater forcing.

Surface Climate
All surface climate anomalies discussed in this section are calculated from the 100 year means for SC_south (with respect to NoSC_south), centered at 14.4 ka.
The strong 14.4 ka AMOC reduction reduces northward flow of relatively warm shallow-intermediate ocean water, leading to northern hemisphere cooling and slight Southern Hemisphere warming (up to +1°C; Figures 2d and 2e). This bipolar anomaly is particularly strong in Boreal winter (December-January-February), when surface air temperatures widely cool by 1-5°C, locally dropping by up to 20°C and 12°C in the GIN and Labrador Seas, respectively.
During Boreal summer (June-July-August), sea ice increases in the North Atlantic-GIN Seas (above 40°N; Figure S1b in the supporting information), amplifying the initial high-latitude cooling through an ice-albedo feedback (+0.1-0.3 albedo change). In Boreal winter, sea ice expansion is limited to the Labrador and GIN Seas ( Figure S1a), corresponding to regions with greatest surface air cooling because the enhanced sea ice thermally insulates the relatively warm underlying ocean water. The elevated albedo (+0.2-0.5 below 70°N) may also produce a modest positive feedback.
Widespread winter cooling reduces precipitation across the Northern Hemisphere, particularly over Greenland and the Labrador and GIN Seas (Figures S2a and S2b). Expanded sea ice over the Labrador and In the tropics, the bipolar surface climate anomaly causes the Intertropical Convergence Zone to shift southward year-round ( Figure S2) in order to scavenge heat [Chiang and Bitz, 2005;Broccoli et al., 2006]. These changes are very similar to the effect of 100 years of 0.4 Sv freshwater hosing across the HadCM3 North Atlantic under Last Glacial Maximum conditions [Kageyama et al., 2013].
Contrasting with overall northern winter temperature reductions is 1-3°C warming over eastern North America (Figure 2d), whereby widespread North Atlantic surface cooling develops higher pressures over the eastern North Atlantic, drawing relatively warm air northward from the Gulf of Mexico. This dipole anomaly of cooling in the high-latitude North Atlantic and warming over North America is broadly reproduced by different freshwater-forced experiments performed with multiple GCMs when similar surface air temperature conditions are met, specifically, high-latitude cooling of more than 5°C and a tongue of cooling reaching westward It is unlikely that such warming would provide a significant positive feedback to Southern Laurentide melt because it is most powerful in the winter, when local temperatures are À40 to À50°C, and is nonexistent in summer. The warming results in greater precipitation locally (Figures S2a and S2b) and could increase southern Laurentide ice sheet mass balance, but coupled climate-ice sheet model simulations are needed to test this.
The temperature increase observed from Panama to Brazil (Figures 2d and 2e) is from regional drying ( Figure  S2), which causes dieback of the broadleaf rainforest to be replaced by C4 grasslands and bare soil. Locally, this acts as a positive feedback, amplifying the warming and drying. The signal arises from large seasonal latitudinal swings in Amazonian precipitation belts, resulting in some regions struggling to achieve sufficient year-round moisture to maintain a tropical rain forest. The same pattern is simulated by Community Climate System Model 1.4-carbon (CSM1.4-carbon) when AMOC is weakened by northern-sourced freshwater [Bozbiyik et al., 2011]. However, this signal may not be robust; a revised representation of leaf respiration (based on the Joint UK Land Environment Simulator (JULES) [Clark et al., 2011]) reduces the sensitivity of the rainforest but shows that it does not significantly impact the main results (away from the Amazon) discussed above.

Discussion and Conclusions
To assess the effect of the North American ice saddle collapse, we developed a physically consistent numerical approach whereby a transient GCM climate simulation was used to drive a dynamical ice sheet model, which in turn drove a high-resolution drainage model, providing freshwater forcing scenarios to plug back into a GCM.
We find that half of the melt from the Cordilleran-Laurentide saddle collapse is routed to the Arctic Ocean. Biases in the timing and distribution of meltwater routing from the ice sheet model (possibly passed through from the initial GCM used to drive the dynamical ice sheet) introduce some temporal error in the routing of southern Laurentide ice sheet melt. However, our results show that the climatic response is insensitive to this error, which affects whether meltwater is routed to the Mississippi River or the North American East Coast. Furthermore, in agreement with previous research [Condron and Winsor, 2012], the climate model displays greater sensitivity to Arctic-bound meltwater than outlets draining directly into the North Atlantic or via the Gulf of Mexico.
These results have been quantitatively derived from a feasible, dynamically simulated mechanism of acceleration in North American ice sheet melt that is supported by geological data [Dyke, 2004;Tarasov and Peltier, 2004;G12;Tarasov et al., 2012;Gomez et al., 2015;Wickert, 2016]. They provide compelling evidence that through the production of meltwater, the separation of the Cordilleran and Laurentide ice sheets was capable of weakening the AMOC (À6 Sv in HadCM3) and cooling the Northern Hemisphere (up to À6°C in the winter over central Greenland).
In our simulations, both the meltwater pulse and the climatic response (Figures 1 and 2) last approximately 600 years, but the timing, duration, and amplitude of the cooling event modeled here are affected by uncertainties in the timing, duration, and amplitude of the Cordilleran-Laurentide ice saddle collapse meltwater pulse. Quantification of this uncertainty by G16 demonstrated that the Bølling warming can trigger the saddle collapse, with peak melting occurring at the maximum warming (14.5 ka) or a few centuries later. The resulting meltwater pulse is 200-600 years long and~70-100% of the magnitude used in SC_South (G16). In this scenario (Scenario1; Figure 3b), a shorter meltwater pulse would likely result in a shorter climatic signal. However, the magnitude of cooling would not necessarily decrease, since changes in ocean circulation (and Greenland temperature) may depend more on discharge rate than total meltwater volume. Another possibility arising from G16 is that the saddle collapse happened earlier than the Bølling warming (also compatible with geological records [Dyke, 2004;Munyikwa et al., 2011]), producing a 400-600 yearlong pulse at~16 ka that was 15-70% of the magnitude used in SC_South. We could assume that this would produce a climatic signal 15-70% of the SC_South magnitude (Scenario2; Figure 3b), though again, this is likely affected by meltwater discharge rate. Simulations run by G16 show that shorter saddle collapse events can release a similar volume of meltwater as longer ones. Furthermore, shorter events have faster meltwater Geophysical Research Letters 10.1002/2016GL071849 discharges, which could result in a strong[er] AMOC reduction and Greenland cooling, even if the total volume of meltwater is less (and vice versa). Hence, it is possible that the climatic impact of a shorter saddle collapse would be similar to the impact of a longer one.
Scenario1 would explain the cooling period after the Bølling warming and/or the Older Dryas. Scenario2 could correspond to cooling recorded in Greenland~16 ka. Both scenarios are feasible from a surface climate perspective (Figure 3), but Scenario1 could also explain 40-60% of the MWP1a rapid sea level rise at 14.5 ka (12-22 m in 340 years [Deschamps et al., 2012]). This is further supported by high-resolution sea surface temperature records from North Atlantic faunal assemblages [e.g., Thornalley et al., 2010], which show pronounced cooling for a few hundred years following the abrupt Bølling warming.
Since the Cordilleran-Laurentide ice saddle collapse only accounts for 40-60% of MWP1a, it is possible that Antarctica also made a significant contribution [Golledge et al., 2014], compatible with sea level fingerprinting scenarios [Gomez et al., 2015;Liu et al., 2016]. Swingedouw et al. [2008] noted that a Southern Ocean bound meltwater flux at <0.2 Sv invigorates NADW formation by reducing Antarctic Bottom Water formation (deepening the Atlantic pycnocline and reducing Southern Ocean density) and increasing Southern Hemisphere wind stress, with both processes acting to enhance NADW export. The precise influence of such additional Antarctic melt on the scenario simulated here has yet to be tested; it may reduce the duration and possibly the magnitude of the effect, pushing the surface climate signal closer to observed events (Figure 3b). In conclusion, we propose a chain of events whereby the abrupt Bølling warming triggered the North American ice saddle collapse as demonstrated by Gregoire et al. [2016]. Delivery of the resultant meltwater to the oceans subsequently ended the Bølling warming and/or caused the Older Dryas cooling via a multicentennial reduction in AMOC and surface temperature (Scenario1; Figure 3b) that would have lasted 200-400 years.