Acceleration of Northern Ice Sheet Melt Induces AMOC Slowdown and Northern Cooling in Simulations of the Early Last Deglaciation

The cause of a rapid change in Atlantic Ocean circulation and northern cooling at the onset of Heinrich Stadial 1 ~18.5 ka is unclear. Previous studies have simulated the event using ice sheet and/or iceberg meltwater forcing, but these idealized freshwater fluxes have been unrealistically large. Here we use a different approach, driving a high‐resolution drainage network model with a recent time‐resolved global paleo‐ice sheet reconstruction to generate a realistic meltwater forcing. We input this flux to the Hadley Centre Coupled Model version 3 (HadCM3) climate model without adjusting the timing or amplitude and find that an acceleration in northern ice sheet melting (up to ~7.5 m/kyr global mean sea level rise equivalent) triggers a 20% reduction in the Atlantic Meridional Overturning Circulation. The simulated pattern of ocean circulation and climate change matches an array of paleoclimate and ocean circulation reconstructions for the onset of Heinrich Stadial 1, in terms of both rates and magnitude of change. This is achieved with a meltwater flux that matches constraints on sea level changes and ice sheet evolution around 19–18 ka. Since the rates of melting are similar to those projected for Greenland by 2200, constraining the melt rates and magnitude of climate change during Heinrich Stadial 1 would provide an important test of climate model sensitivity to future ice sheet melt.


Introduction
The early stage of the last transition from glacial to interglacial climate was characterized by rapid weakening of the Atlantic Meridional Overturning Circulation (AMOC; McManus et al., 2004), Northern Hemisphere cooling (Bard et al., 2000;Buizert et al., 2014), and warming of the Southern Hemisphere (Parrenin et al., 2007). This period of climate change, known as "Heinrich Stadial 1" (or even the "mystery interval"; Denton et al., 2006), remains enigmatic. Historically, the Heinrich Stadial 1 AMOC slowdown and associated surface climate change (northern cooling) have been attributed to centennial-scale meltwater pulses from Heinrich Event 1 at~16 ka (Broecker, 1994). During this event, catastrophic iceberg calving led to huge volumes of ice rafting out across the North Atlantic, where it melted in the open ocean in much more southerly locations than usual (Hemming, 2004;Ruddiman, 1977). The resulting freshwater perturbation would have increased surface buoyancy and thus may have inhibited North Atlantic Deep Water (NADW) formation and weakened the AMOC (Rahmstorf, 1994).
However, there is evidence that the onset of the Heinrich Stadial 1 AMOC slowdown started~18.5 ka (Ng et al., 2018;Stern & Lisiecki, 2013), much earlier than Heinrich Event 1 Hodell et al., 2017;Stern & Lisiecki, 2013). This chain of events, whereby the onset of the Heinrich Stadial (a multimillennial period of low Northern Hemisphere temperatures, which may be linked to a relatively weak or collapsed AMOC) preceded the Heinrich Event (a shorter period of enhanced iceberg calving, captured by sedimentary ice rafted debris records in the North Atlantic), is consistent with earlier time periods (Barker et al., 2015;Bond & Lotti, 1995); note the distinction between Heinrich Stadials and Heinrich Events. Furthermore, climate models suggest that Heinrich Stadial-like conditions can only be simulated during or close after glacial climates if unrealistically large iceberg-derived freshwater fluxes are adopted (e.g., compared to estimates by Roberts et al., 2014). What, then, could have forced the AMOC slowdown and northern cooling that began~18.5 ka?
One possible solution lies in the vast ice sheets that covered much of North America and Eurasia at the time. These ice sheets had begun to melt due to increasing Northern Hemisphere summer insolation (Dyke, 2004;Gregoire et al., 2015;Hughes et al., 2016;Roche et al., 2011) and together they accounted for~90 m global mean sea level rise of freshwater or more (Lambeck et al., 2014;Peltier et al., 2015). While most of this ice melted later in the deglaciation (Dyke, 2004;Gomez et al., 2015;Gregoire et al., 2012;Hughes et al., 2016;Ivanovic et al., 2016;Lambeck et al., 2014;Peltier et al., 2015;Tarasov et al., 2012), it is plausible that some melting in the earliest stages of the deglaciation (e.g., Toucanne et al., 2010Toucanne et al., , 2015Zaragosi et al., 2001) may have been sufficient to drive the Heinrich Stadial 1 ocean circulation and climate changes, as suggested by Hall et al. (2006). This possibility is enhanced if the melting was routed directly to sites of NADW formation, where it would be more effective in stabilizing the water column and reducing the AMOC before mixing dissipated the freshwater anomaly (Condron & Winsor, 2012;Roche et al., 2009).
Other model studies have simulated Heinrich Stadial 1 by adjusting the timing and amplitude of ice meltwater fluxes to reproduce proxy records of the AMOC and climate change (Z. Liu et al., 2009;Menviel et al., 2011), thereby fingerprinting the possible sources of freshwater forcing. However, in these simulations, the volume of meltwater needed to satisfy climate and AMOC proxies is approximately two to three times the estimates of global mean sea level rise (e.g., Peltier et al., 2015;Figure 1). This discrepancy between the previously used model forcing and sea level records may have two causes. First, the magnitude of the target AMOC reduction may have been overestimated (Bradtmiller et al., 2014) since it is difficult to calculate the quantitative change in the AMOC from existing ocean circulation proxies. Second, numerical limitations led to simplifications in defining the domains of meltwater entry to the ocean, which are unknown and thought to have changed over these timescales (e.g., Duplessy et al., 1988;Meland et al., 2008). This is important because the location of freshwater input affects the efficiency of the forcing (Condron & Winsor, 2012;Roche et al., 2009;Smith & Gregory, 2009).
In this paper we investigate the possibility that Heinrich Stadial 1 was forced by melting of the Northern Hemisphere ice sheets early in the deglacial phase, rather than by Heinrich Event 1 icebergs. Following Ivanovic et al. (2017), we adopt meltwater fluxes from a recent time-resolved global ice sheet reconstruction spanning the period (ICE-6G_C VM5a; Peltier et al., 2015) without adjusting the timing or amplitude of the fluxes. Next, we account for the spatial distribution of melting and subsequent freshwater delivery to the ocean by combining a model of past sea level change (e.g., Kendall et al., 2005) with a high-resolution drainage network routing model (Wickert, 2016;Wickert et al., 2013). This methodology yields a realistic meltwater scenario for the period 19-15 ka that, unlike previous studies, is physically consistent with the history of ice sheet evolution ( Figure 1). Finally, we input this global, transient freshwater forcing into the The meltwater fluxes used in this study were calculated according to the ICE-6G_C (VM5a) reconstruction, consistent with global mean sea level data (Peltier et al., 2015). For comparison, the total freshwater fluxes used by Menviel et al. (2011;all input to the North Atlantic) and Liu et al. (2009;all input to the North Atlantic and Gulf of Mexico) are also shown, alongside Greenland melting rates measured in modern times (0.74 mm ± 0.14 mm/year GMSLRE; McMillan et al., 2016) and projected for 2200 (6.7 mm/year GMSLRE; Lenaerts et al., 2015). coupled ocean-atmosphere-vegetation general circulation model (GCM) Hadley Centre Coupled Model version 3 (HadCM3; Gordon et al., 2000;Pope et al., 2000;Valdes et al., 2017) to simulate the effect on ocean circulation and climate.
The climate model setup and experiment design is described in section 2, including details of the meltwater histories used to drive the simulations, as well as our methods for comparing model output with climate and ocean circulation proxy data. A presentation of the main results follows in section 3, and a critical evaluation of the timing of the meltwater forcing and simulated climate event is in section 4. Section 5 compares the model results to climate and ocean circulation proxy data from Heinrich Stadial 1, while section 6 gives an overview of the potential impact of changes in the orbital forcing and rising atmospheric CO 2 . Section 7 discusses issues of model sensitivity and stability. Section 8 presents a possible link between Heinrich Stadial 1 and Heinrich Event 1, as suggested by the climate model results. Finally, section 9 summarizes our conclusions.

Climate Model Description
The climate simulations were carried out using the BRIDGE version of the HadCM3 coupled atmosphereocean-vegetation GCM, described in detail by Valdes et al. (2017). Briefly, the atmosphere model comprises 19 hybrid sigma-pressure coordinate vertical layers with a horizontal resolution of 2.5°× 3.75° (Pope et al., 2000) and is coupled to the Met Office Surface Exchange Scheme (MOSES) version 2.1, including the Topdown Representation of Interactive Foliage and Flora including Dynamics (TRIFFID) vegetation model (Cox, 2001). Once per model-day, the atmosphere GCM is coupled to the ocean GCM, which comprises 20-vertical layers on a depth (z) coordinate system designed to give maximum spatial resolution at the surface (Johns et al., 1997, Table 2) and a horizontal resolution of 1.25°× 1.25° (Gordon et al., 2000). The model is computationally efficient enough to feasibly run multimillennial simulations, while also having the complexity and physical-robustness of a full GCM.

Meltwater Fluxes and Climate Simulation Design
To generate realistic ice sheet meltwater fluxes to input to the climate model, we follow the procedure employed by Ivanovic et al. (2017Ivanovic et al. ( , 2018, with the following revisions: (1) The global drainage network model (Wickert, 2016;Wickert et al., 2013) was driven by sea level changes computed using the glacial isostatic adjustment (GIA) model ICE-6G_C (VM5a; Peltier et al., 2015; henceforth, "ICE-6G_C") for the period spanning 19.0-15.5 ka (see section 2.3). The GIA model is composed of both an ice history and a coupled (Maxwell) viscoelastic Earth model (VM5a) that yields reasonable fits to a suite of geological markers of sea level change, present-day crustal motions measured by Global Positioning System (GPS), and anomalies in the rotational state and gravity field. GIA-induced sea level changes were based on a gravitationally self-consistent sea level theory that accounts for migration of shorelines, the evolution of grounding lines in marine-based settings, and perturbation in sea level associated with load-induced changes in Earth rotation (Kendall et al., 2005). ICE-6G_C is one of two recently published global ice sheet reconstructions; the other is GLAC-1D. The latter is also constrained by surface observations (see Ivanovic et al., 2016, and citations therein for more details, including information on how the two reconstructions differ), but we adopted ICE-6G because at the time of the study the history of ice thickness changes and the associated meltwater routing had not been published for GLAC-1D.
The distribution of ice thickness through time is required both for calculating the drainage network routing and for computing the volume of meltwater discharged through the rivers via the ice mass balance (e.g., Figures S1 and S2). We chose to investigate the time window 19.0-15.5 ka to include the onset of Heinrich Stadial 1 (~18.5 ka) and a subsequent period of accelerated ice sheet melt at 16.5-15.5 ka in ICE-6G_C. Inputting the time-varying GIA-induced sea level (and topography) changes into the drainage network model yielded a fully distributed, global freshwater discharge, with 10-arc second spatial resolution and a time step of 500 years (which is the time step of the ICE-6G_C ice history). Thus, we produced a spatially detailed meltwater flux that is physically consistent with ICE-6G_C and therefore also subject to the same strengths, uncertainties, and possible biases as the parent ice sheet reconstruction. Due to technical difficulties with the implementation of this transient model boundary condition in HadCM3, the time step of the meltwater forcing was stretched to 545.36 years (precisely: 196,331 days) when input to the climate model. Since the magnitude of this stretching is <10% of the time-step of the ICE-6G_C ice sheet history, we consider it to be within the temporal uncertainty of the forcing. There are no problems with misalignment to other transient forcings because all other model boundary conditions (insolation, trace gases, ice sheet geometry, land-sea mask, topography, and bathymetry) are kept constant. During analysis (including comparison to ocean circulation and climate proxy data), we align the center of our 3,450-year climate simulations (i.e., model year 1750) with the center of the utilized ICE-6G_C routing (17.395 ka). Thus, the model results are presented as 19.145-15.645 ka to reduce the small temporal distortion of the results. This distortion has no effect on the overall interpretation of the results.
(2) In the climate simulations we adopted a previous 1,750-year long HadCM3 integration that was initialized for 21 ka (Singarayer et al., 2011) and updated the model's boundary conditions to be in agreement with current knowledge (results from the same family of simulations are presented by Morris et al., 2018, andSwindles et al., 2018). Atmospheric trace gas concentrations were set at 188 ppmv for CO 2 (Lüthi et al., 2008), 375 ppbv for CH 4 (Loulergue et al., 2008), and 200 ppbv for N 2 O (Schilt et al., 2010). The ice sheets and associated parameters (ice mask, surface topography, land sea mask, and ocean bathymetry) were updated according to ICE-6G_C at 21 ka using an anomaly method (as per Singarayer & Valdes, 2010) to maintain consistency with the standard preindustrial setup of the model . The simulation was then run with this "equilibrium" setup (i.e., without any time-evolving boundary conditions or forcings) for a further 200 years, by which time the model had reached a near-steady climate state. The 3,450-year long, transient-meltwater simulation presented in this study was then configured and initialized from the end of this 1,950-year long spin-up phase. To verify that the simulated changes (presented and analyzed below) were not the result of incomplete model spin-up (Marzocchi & Jansen, 2017), we also continued the 21 ka equilibrium simulation unchanged for a further 3,450 years. This extended equilibrium simulation confirmed that by the start of the transient-meltwater run, the model had reached a near-steady climate state and no drifts exist in the climate signals discussed herein. Furthermore, the extended equilibrium simulation demonstrated that the model response is forced by the transient meltwater scenario and not by internal variability. (3) HadCM3 has a rigid lid ocean, which means that hydrological fluxes are represented by virtual salinity fluxes (section 2.3 in Ivanovic et al., 2017). Consequently, point-source freshwater fluxes (as produced by the high-resolution drainage network model) can result in individual grid cells becoming so fresh that they exceed the valid range for the ocean model's equation of state (Bryan & Cox, 1972;Fofonoff, 1962;Fofonoff & Millard, 1983). We therefore adapted the method used by Ivanovic et al. (2017Ivanovic et al. ( , 2018 to regrid time-varying meltwater fluxes onto the HadCM3 ocean using the 21 ka land-sea mask and bathymetry, ensuring that all meltwater reached the coast. We then spread the freshwater over adjacent ocean grid boxes with depths >500 m along the coastline and off the continental shelf. In this exercise we used previous high-resolution model results (e.g., Condron & Winsor, 2011) and our experience with freshwater perturbations in the HadCM3 model to choose an appropriate mask for this redistribution ( Figure S3).
Note that meteoric water catchments are not adjusted in these simulations and the routing of ice sheet meltwater described in this section is varied independently of the model's rivers. Furthermore, although all forms of ice mass loss are captured in the ice sheet history used to generate the meltwater fluxes (e.g., surface ablation, calving at the marine margin, etc.), iceberg calving is not directly represented in our model. The location of iceberg melting is partly accounted for by the redistribution mask used in the simulation ( Figure S3) and a standard-HadCM3 iceberg melting mask around Greenland and Antarctica, but rafting further into the open ocean is not included in our forcing. Additionally, there is no thermal heat flux associated with the freshwater forcing. These factors may result in the model being over-sensitive to the iceberg component of the freshwater discharge (Levine & Bigg, 2008), although this is difficult to know due to the competing effects of associated freshwater delivery and localized cooling. Moreover, the iceberg component is small compared to ablation in the continental interior of the ice sheets (e.g., Figure 2) and much of our forcing enters the same regions in which icebergs are thought to have melted during the Last Glacial Maximum and early stages of Heinrich Stadial 1 (e.g., Stern & Lisiecki, 2013;Watkins et al., 2007). Thus, while the absence of explicit icebergs is a limitation of the methodology, it is likely of second-order importance to this study.

Model-Geological Data Comparison
We compare our model results to a range of climate and ocean circulation proxies to evaluate the extent to which the model reproduces the paleoclimate record. The following provides a description of the relevant data sets, and the method we use to compare the climate and ocean circulation proxy data to model results.
We calculate the simulated maximum AMOC strength between depths of 204 and 5,500 m, and between latitudes of 19.375 and 80.625°N. This spatial range captures the maximum AMOC stream function and is not skewed by any possible shifts in the AMOC position during the simulation. We compare the results to proxy observations of AMOC strength consisting of a recent compilation of North Atlantic sedimentary 231 Pa/ 230 Th records (Ng et al., 2018). The use of 231 Pa/ 230 Th as an AMOC rate proxy is based on the deviation of the sedimentary 231 Pa/ 230 Th ratio from the oceanic production value (0.093) due to changes in the rate of oceanic transport, and given the difference in residence time between 231 Pa and 230 Th (Yu et al., 1996). The 231 Pa/ 230 Th data used in the analysis are derived from an integrated data set of North Atlantic sediment records, which exhibit coherent trends indicative of variations in AMOC strength (Ng et al., 2018). Despite the evident link, quantifying the AMOC rate from 231 Pa/ 230 Th observations remains a work in progress given the many factors beyond oceanic advection that influence oceanic 231 Pa and 230 Th cycling and the incorporation of these isotopes into sediments (e.g., Burke et al., 2011;Deng et al., 2014;Hayes et al., 2015).
For sea surface temperature (SST) records, we compare our simulation to three North Atlantic reconstructions (Bard et al., 2000;Calvo et al., 2001;Naafs et al., 2013) that are derived from measurements of alkenone unsaturation index (U K 0 37 ; Prahl & Wakeham, 1987) and calibrated to mean annual temperature with the global coretop calibration of Müller et al. (1998). The records are presented with their initial published age models, with the exception of the data from core SU8118 (Iberian Margin; Bard et al., 2000), which we present on the updated age model of Bard et al. (2004). The reconstructed SST data are compared to simulated SST in a 3 × 3 grid of ocean model surface cells centered on the location of the measurement (Bard et al., 2000;Calvo et al., 2001;Naafs et al., 2013); Table 1 provides the exact coordinates of the records and model cells. These 3 × 3 grids cover an area large enough to account for limitations of the model's spatial resolution.
The Greenland surface air temperatures used in the analysis are taken from the most recent reconstruction based on gas phase δ 15 N-N 2 , δ 18 O in ice, and a firn densification model . The reconstruction includes measurements from three ice cores: NEEM, GISP2, and NGRIP, which span central Greenland  Figure S2). Model land is masked in gray with modern coastlines outlined in black.

Paleoceanography and Paleoclimatology
IVANOVIC ET AL.
( Table 1). Similar to the SSTs, these data were then compared to a 3 × 3 atmosphere grid at the Earth's surface (Table 1), which was aligned to the range of coordinates covered by the compilation, thus spanning central Greenland.

Simulated Changes in North Atlantic Ocean Circulation and Climate
According to the ICE-6G_C reconstruction (Peltier et al., 2015), both the North American and Eurasian ice sheets lost mass between 19 and 15.5 ka, particularly along the southern margin of the Laurentide Ice Sheet, around the perimeter of the Cordilleran Ice Sheet and over the Barents-Kara Sea (e.g., Figure 2, which shows ice mass loss for the period 17-15.5 ka, when the rate of Northern Hemisphere deglaciation was greatest in this time period, Figure 1). This ice sheet loss was likely initiated by an increase in summer insolation with positive feedback from changes in surface albedo, ice sheet elevation and atmospheric CO 2 (Gregoire et al., 2012Roche et al., 2011), and potentially from internal ice sheet dynamics in marine sectors of the Eurasian ice sheet (Patton et al., 2017).
At the start of the simulation, AMOC is active (max. strength~16.5 Sv), but high-latitude North Atlantic convection is relatively shallow (<1,800 m), and Antarctic Bottom Water fills the deep-abyssal basin ( Figure S4). Although there is a small acceleration in forcing from Eurasian ice sheet melt~19 ka (Figure 1), it is insufficient to drive the AMOC slow-down captured in sedimentary 230 Pa/ 231 Th records at this time ( Figure 3b). The rise in melting from the Barents-Kara Sea regions during this period is concurrent with a comparable decrease in rates of North American ice sheet melt, resulting in little change to North Atlantic sea surface salinity and hence negligible impact on water column stability and NADW formation.
On the other hand, the approximate doubling of melt rates for both the North American and Eurasian ice sheets after 17 ka in the ICE-6G_C reconstruction produces enough freshwater to markedly reduce high latitude sea surface salinity (Figure 2), increasing vertical stratification by 20-80% in the upper 200 m (Figure 4). This stabilizes the North Atlantic water column, acting to inhibit vertical convection and thus reduce NADW formation, which weakens the AMOC by 20% (3.5 Sv; Figure 3b).
The simulated AMOC slowdown reduces the transport of tropically sourced water to the eastern North Atlantic, cooling the ocean from the surface (Figure 3c) down to~1,600 m depth. The cold anomaly is spread by the subpolar and subtropical gyres, propagating to the overlying atmosphere ( Figure 5). One consequence of the widespread high-latitude Atlantic surface cooling is the expansion of northern sea ice cover, especially in the Boreal winter-spring, when sea ice skirts Greenland and Iceland and spreads south from Baffin Bay into the Labrador Sea. This growth in sea ice cover causes strong winter surface air cooling where the sea ice thermally isolates the atmosphere from the ocean (Figure 5c). Surface winds recirculate cool air from Baffin Bay and the Labrador Sea and drive the cold anomaly northward over central Greenland. This Greenland cooling is not as strong in the summer as in the winter-spring (Figure 5d) due to less expansive sea ice cover. Thus, Greenland's seasonality is amplified in the model, consistent with ice core and glaciological evidence of enhanced seasonality and extreme winter temperatures during Northern Hemisphere stadials Denton et al., 2005), although the onset of this climatic change is late compared to the start of Heinrich Stadial 1 (~18.5 ka). In response to the northern cooling and expansion of winter sea ice, the Intertropical Convergence Zone (ITCZ) and associated tropical rain belt are driven southward (Figures 5e and 5f), weakening the Indian summer monsoon and causing seasonal drying over sub-Saharan Africa, Central America, and parts of the Amazon (consistent with Kageyama et al., 2013).

Paleoceanography and Paleoclimatology
It is clear in the simulation that sea ice growth during the northern cooling has a strong seasonal impact on North Atlantic winter surface air temperature, which has important implications. Specifically, if winter sea ice is too extensive in the climate model, its presence may skew the local surface air temperature to be too low, seasonally. Furthermore, expanded winter sea ice reduces local evaporation from the ocean surface, decreasing local moisture availability for rain/snowfall. Coupled with the influence of surface temperature on moisture availability (whereby lower temperatures reduce evaporation through the Clausius-Clapeyron relation), high-latitude North Atlantic precipitation and evaporation anomalies would be affected by any bias in the simulated sea ice behavior.
It is difficult to constrain sea ice extent and thickness from geological records, but we can consider the possible impact of such changes on the results. We have compared precipitation-weighted surface air temperature to the Greenland ice core record of Buizert et al. (2014), see Figure 3f. This reconstruction in the early deglaciation is based largely on the δ 18 O of ice at that time , which is sourced from precipitation (e.g., snowfall), and as such may reflect precipitation-weighted temperatures (Persson et al., 2011;Sime et al., 2008). This temperature signal from ice cores would then be biased toward times of higher precipitation, usually the summer season, especially during times of amplified seasonality, such as at 17.5 ka in our simulations. During the very cold winter, temperatures were too low (and sea ice cover was too extensive) for there to be sufficient snowfall to record the strong cooling (Figure 5c), and much less cooling occurred in the summer (Figure 5d) when there was more precipitation. Therefore, it is more appropriate to compare the ice core-derived temperature reconstruction to the precipitation-weighted temperature change simulated over Greenland (which is summer biased). Using this method essentially removes any seasonal bias in simulated Greenland temperatures that arise from a potentially overly cold winter and annual mean ( Figure S5), which may be related to the overestimation of sea ice extent. For example, if winter sea ice growth was less than the climate model simulates, winter temperatures would not be so low, decreasing the seasonality of the anomaly and reducing the annual mean surface air temperature anomaly (darker blue line in Figure 3f) toward the smaller summer anomaly (Figure 5c).
Another consequence of the strong winter cooling in the mid-high latitude North Atlantic is a shift in European winds (Figure 6), which follows an increase in the pressure gradient. This brings relatively warm, southerly air northwards over Fennoscandia and causes weak surface air warming of <1.7°C (Figure 5c). Because much of Fennoscandia was covered in an ice sheet during this period, it is difficult to verify from the geological record whether such a warming took place. If it did, it likely would not have affected local ice sheet mass balance, because the small temperature increase is restricted to the Boreal winter, when surface air was already very cold. However, the warming signal is likely to be highly influenced by the strong sea ice feedback on North Atlantic surface air temperatures, which means that it may be specific to the sea ice distribution simulated by the climate model. On the other hand, Fennoscandian warming during Heinrich Stadial 1 is also predicted by another model (Z. Liu et al., 2009), possibly also due to

10.1029/2017PA003308
Paleoceanography and Paleoclimatology extended winter sea ice cover, and so this may be a robust climatic feature. The multimodel experiment for the last deglaciation (Ivanovic et al., 2016) will shed light on this potential model-dependency.
The increase in winter temperatures over Fennoscandia is not the only seasonal warming signal in our results; the 20% reduction in AMOC strength causes a weak bipolar climate anomaly, with <1°C sea surface warming in the Southern Hemisphere (Figures 5a and 5b) and <0.5°C warming in most of the overlying surface air temperature (Figures 5c and 5d). While this is a relatively small change compared to the North Atlantic surface anomalies, it is sufficient to push the Southern Ocean winter sea ice edge a few degrees southward, reducing thermal insulation around 60°S and causing weak localized warming of the overlying air (up to 1.7°C; Figure 5d). Similar to the sea ice effects discussed above, the exact magnitude and location of this warming may be model-specific. However, warming of the Southern Hemisphere in response to an AMOC slowdown has been widely observed using other climate models (Kageyama et al., 2013;Z. Liu et al., 2009;Menviel et al., 2011;Rahmstorf, 1994), and thus broadly, this effect is likely robust.
Less seasonally, the reduction in the AMOC and mid-high latitude vertical convection allows heat to build up in the subsurface North Atlantic south of 45°N, which extends to the surface between 35 and 40°N (Figures 5a  and 5b). The surface warming enhances evaporation, thus increasing atmospheric moisture availability and local precipitation, whereas cooling south and east of the warming has the opposite effect.

Timing of Heinrich Stadial 1 in the Model
Global reconstructions of ice sheet evolution from 21 ka to present-such as the ICE-6G_C data set and the GLAC-1D reconstruction discussed in section 2.2-incorporate constraints from a wide variety of sources (e.g., see Ivanovic et al., 2016, for a summary), some of which contain substantial uncertainties. A comparison of the history of global mean sea level in ICE-6G_C and other available reconstructions provides a measure of this uncertainty. Figure 7 shows the global mean sea level rise equivalent (GMSLRE) of ice mass loss in ICE-6G_C and GLAC-1D, and a sea level reconstruction by the Australian National University sea level group (Lambeck et al., 2014). Notably, both GLAC-1D and the Australian National University reconstruction show higher rates of global mean sea level change between 19 and 18 ka than occurs in ICE-6G_C (although the melting in all reconstructions is still significantly lower than the rates employed by Liu et al., 2009 andMenviel et al., 2011). This highlights the possibility that ICE-6G_C may underestimate Northern Hemisphere ice sheet melt during the period when AMOC slowdown (and northern cooling) is thought to have begun. In the construction of the ICE-6G_C model, the GMSLRE of ice volume is tuned to fit sea level records in the far field of the major ice sheets, in particular at Barbados (Peltier et al., 2015). However, the relatively shallow slope of the coral record at Barbados in the period prior to~16 ka and the large indicative range of the corals indicate that the timing of the global mean sea level rise in the ICE-6G_C model can easily be in error by several thousand years without misfitting the coral record (see Figure 2 of Peltier et al., 2015).
With regard to the North American ice cover, ICE-6G_C is constrained to match the margin chronology of Dyke (2004). However, at present, glacial geological data cannot significantly distinguish North American  Figure 3. Cross-hatching marks significance below 95% confidence in the simulation using a student t test. See Figure S5 for corresponding annual mean anomalies.

Paleoceanography and Paleoclimatology
ice sheet changes between 18.5 and 17.0 ka. The areal extent of the North American ice sheet changed little over these 1,500 years (Dyke, 2004), and this, combined with fundamental issues with relating continental ice sheet extent to volume when it lies on a deformable bed (Carlson & Clark, 2012;Clark & Walder, 1994;Iverson & Petersen, 2011), make its volume at either time, and its volumetric change, difficult to estimate. In addition, ice-marginal dates are often sparse and geographically concentrated (Dyke, 2004), and many record the millennial-scale response of the ice margin position to internal dynamics (Patterson, 1997) rather than changes in whole ice sheet mass balance. All of this makes it difficult to establish whether or not a rapid change in North American ice volume could have taken place 18.5-17 ka.
The most recent reconstruction of the evolution of Eurasian ice sheet extent, DATED-1 (Hughes et al., 2016), suggests an earlier onset of retreat in the Barents and Kara Sea sectors compared to ICE-6G_C. More precisely, the DATED-1 reconstruction shows that during the early last deglaciation, rates of Eurasian melt were likely twice as large as those in ICE-6G_C, with the largest rates of ice loss occurring between 19 and 18 ka (averaging~4 m/kyr; Figure 7). When added to contemporaneous North American ice sheet melt (which may be a lower limit, considering current limitations in the constraints on North American ice volume change, see above), the resulting freshwater forcing is 6 m/kyr. This is close to the maximum rate of freshwater forcing that drives the AMOC slowdown and northern cooling in our simulation (7.5 m/kyr). Geochemical evidence of elevated rates of continental erosion supports enhanced Northern Hemisphere melt early in the deglaciation (Chen et al., 2016). Furthermore, a peak in meltwater delivery to the ocean between 19 and 18 ka is consistent with sedimentological, geochemical and micropaleontological records of Eurasian meltwater drainage to the Nordic Seas and northeast North Atlantic (Mojtahid et al., 2017;Toucanne et al., 2010Toucanne et al., , 2015. These studies also record a second stronger and extended Eurasian meltwater pulse 18.2-16.7 ka, which could have contributed toward triggering and maintaining a weakened AMOC state.

Link Between Model Results and Heinrich Stadial 1
The climate model simulations (Figure 3) produce a good match to the proxy climate and ocean circulation records from the onset Heinrich  Comparison of ice sheet meltwater fluxes from different reconstructions and freshwater forcing scenarios (expressed as global mean sea level rise equivalent since 21 ka). The thick solid black, red, and green curves are global mean sea level change from the ice sheet and sea level reconstructions, respectively: ICE-6G_C (Peltier et al., 2015; as used in this study), GLAC-1D (Ivanovic et al., 2016;Tarasov et al., 2012) and the sea level reconstruction produced by the Australian National University group (Lambeck et al., 2014), which includes eustatic sea level data (±1σ; vertical green bars). The blue curves are the Eurasian contributions to sea level change in ICE-6G_C (solid), GLAC-1D (dot-dashed), and as inferred from the DATED-1 ice extent chronology (dashed; Hughes et al., 2016). The Northern Hemisphere freshwater fluxes used in the climate modeling studies of Menviel et al. (2011) and Liu et al. (2009) are shown in black dashed lines and dot-dashed lines, respectively (same as Figure 1).

Paleoceanography and Paleoclimatology
Stadial 1, in terms of both amplitude and rate of change, but the modeled change happens~1,500 years later (17 ka) than all the changes recorded in the climate and ocean circulation proxies (~18.5 ka). It is possible that imprecision in the exact location of freshwater input to the ocean and/or in the simulated location of NADW formation could account for this discrepancy by reducing the efficiency of the forcing. On the timescale of the first few years of applying the meltwater flux, its site of entry to the ocean has been shown to be important both for where the water ends up and its potential to affect AMOC (e.g., Condron & Winsor, 2011. However, after a few decades, continued dispersal of the anomaly means that the centennial-scale impact on AMOC and climate is not so influenced by such biases. There may be some differences in the detailed climate response (Roche et al., 2009), but these are of the order of model dependency (e.g., Kageyama et al., 2013) and the broad-scale result of AMOC weakening, northern surface cooling, and subsurface warming (reaching the surface in the western subtropical Atlantic) is more robust.
Considering uncertainties in ice sheet histories and the latest constraints on Eurasian ice sheet evolution (see section 4), a more likely explanation for the late simulation of AMOC weakening is that total melting during the Heinrich Stadial 1 period is biased too young in the ICE-6G_C reconstruction. In other words, the freshwater forcing for the climate model is mistimed and should ramp-up~1,500 years earlier. Thus, in this section, we compare model and climate proxy data on the basis that the event simulated after 17 ka in the model is the onset of Heinrich Event 1~18.5 ka.
The simulated reduction in AMOC strength is compatible with qualitative information inferred from a recent compilation of North Atlantic 230 Pa/ 231 Th measurements~18.5 ka (Ng et al., 2018; Figure 3b), although the rate of AMOC weakening cannot [yet] be reliably quantified from this ocean circulation proxy. This is a smaller change in the AMOC than that targeted by previous climate model studies (Z. Liu et al., 2009;Menviel et al., 2011), yet the effect on surface climate is strong and consistent with a range of geological records (Figures 3b-3f), demonstrating that a weakening of the AMOC and not necessarily a collapse or near-collapse is sufficient to explain surface climate signals during Heinrich Stadial 1.
The model yields absolute SSTs that are consistently a few degrees cooler than reconstructed values, but that are within uncertainty of the absolute temperature based on scatter in the global core top calibration (Müller et al., 1998). Moreover, the pattern of SST change at the onset of Heinrich Stadial 1, 18.5 ka, is very well reproduced in the model from 17 ka (Figures 3c-3e and 5a and 5b). Broadly, temperatures in the mid-high latitude North Atlantic drop (e.g., as recorded near the Iberian margin by Bard et al., 2000), while the subtropics show a more complex spatial pattern of surface warming in the west-central basin (shown by Naafs et al., 2013) surrounded by cooling anomalies and no change in between (Calvo et al., 2001). The subtropical surface warming (Figure 3e and 5a and 5b) is approximately 3°too far south and 11°too far west compared to SST proxy records (Naafs et al., 2013). This misfit may be due to imprecision in the simulated Gulf Stream and North Atlantic Drift positions caused by the relatively course horizontal resolution of the model (see section 2.1).
Similar to the surface warming, simulated intermediate-depth warming in the central North Atlantic (Figure 8) may not extend as far north as suggested by benthic temperature proxies (Mg/Ca and δ 18 O; Marcott et al., 2011) at 48°N near Newfoundland. Subsurface warming has been seen in previous modeling of Heinrich events (Z. Liu et al., 2009), although the temperature increase in that study extended through all latitudes, possibly due to the broader meltwater forcing that was adopted. In our simulation, subsurface cooling occurs in the northernmost North Atlantic (Figure 8), due to a combination of reduced northward heat transport (associated with the weakened AMOC) and downward mixing of the surface cooling~60°N. The Figure 8. Change in ocean potential temperature. Annual mean ocean potential temperature anomalies at (a) 30°W (Atlantic basin and adjacent high latitude sectors) and (b) 1 km depth, calculated from the 15.75 to 15.65 ka mean minus the 19.145 to 18.145 ka mean. Cross-hatching marks significance below 95% confidence using a student t test.

10.1029/2017PA003308
Paleoceanography and Paleoclimatology positional shifts required to match our simulated North Atlantic warming with the reconstructions (Marcott et al., 2011;Naafs et al., 2013) are small (on the order of a few degrees north) and are in line with the mismatch between the location of the simulated and reconstructed sea surface warming in the subtropics discussed above. The comparison between modeled and reconstructed subsurface temperature change is otherwise very good, in terms of both the rate of change and peak magnitude.
The slight cooling simulated over Greenland from 17 ka matches well with the reconstructed trend from 18.5 ka  when temperature is weighted by precipitation to account for enhanced seasonality (Figure 3f; see section 3). The simulated southward shift of the ITCZ is also consistent with a reduction in precipitation during Heinrich Stadial 1 as recorded in sediments from the Cariaco basin (Deplazes et al., 2013).
These results indicate that beginning at 17 ka, the model produces a coherent set of climate signals driven by feasible rates of ice sheet melt that are consistent with the recorded patterns of change that began~18.5 ka. There is little evidence in the geological record that any of the climate and ocean circulation changes discussed above took place between 17 and 15.5 ka, and this further suggests that while being realistic in terms of rates and amplitude, the timing of the meltwater forcing in ICE-6G_C is biased young within this time period.

Other Forcings
In our simulations, the only transient forcing is meltwater; the geometry of the ice sheets (ice extent and elevation), atmospheric greenhouse gases, and insolation parameters associated with Earth's varying orbit have all been kept constant. This approach isolates the impact of meltwater effects on ocean circulation and climate from the other forcings, testing the hypothesis that an acceleration in northern ice melt triggered the onset of Heinrich Stadial 1. While a full investigation into the potential impact of these other factors is beyond the scope of this study, it is important to consider their possible influence when evaluating the results.
North American ice sheet geometry has been shown to play an important role in influencing North Atlantic ocean circulation and surface climate Gregoire et al., 2018;Löfverström & Lora, 2017;Zhang et al., 2014). For example, it has been suggested that the separation of the Laurentide and Cordilleran ice sheets led to major shifts in atmospheric circulation patterns (Löfverström & Lora, 2017). The lowering of the Laurentide ice sheet could also have led to a weakening of North Atlantic gyre circulation and the AMOC (À4.5 Sv between 21 ka and present day; Gong et al., 2015). However, ice sheet surface extent (Dyke, 2004) and maximum elevation (as reconstructed in the global ice sheet models; Peltier et al., 2015;Tarasov et al., 2012) do not undergo major transitions over the period we have focused on, suggesting that they did not have a major effect on the AMOC and surface climate at the start of Heinrich Stadial 1.
Atmospheric CO 2 began to rise from~18 ka (Bereiter et al., 2015) and, due to its role as a greenhouse gas, this provided an important feedback to orbitally triggered warming and the disintegration of the North American and Eurasian ice sheets during the last deglaciation Roche et al., 2011). Therefore, in addition to having a direct effect on surface temperatures, both CO 2 and Northern Hemisphere summer insolation impacted climate and ocean circulation by driving ice sheet retreat and thus generating meltwater. In terms of their direct effect, neither are likely candidates for driving the Heinrich Stadial 1 AMOC slowdown and northern cooling~18.5 ka because (i) there is no known mechanism by which the increase in Boreal summer insolation would cause such changes, (ii) the rise in atmospheric CO 2 began~500 years after the onset of the event, and (iii) the timing and magnitude of CO 2 and Northern Hemisphere summer insolation changes are very well constrained (e.g., Bereiter et al., 2015;Berger & Loutre, 1991).
Nevertheless, CO 2 and orbital changes may be important for the detailed pattern of climate change observed during this time period, especially outside of the Atlantic sector where the direct effect of meltwater is weaker. For example, the Southern Hemisphere ocean warming simulated in this study is not strong enough to cause all of the~1.5°C surface air warming recorded by 17.5 ka in Antarctic ice cores (Parrenin et al., 2007). Assuming that the meltwater-driven climate signal does represent the earlier Heinrich Stadial 1 (section 4), the implication is that a stronger meltwater-induced bipolar see-saw is required (possibly from additional forcing by icebergs, not included in this study) and/or that dynamical atmospheric and oceanic responses to rising greenhouse gases (Bereiter et al., 2015;Marcott et al., 2014) play a key role in Southern Hemisphere temperature change at this time.
Equilibrium-type HadCM3 snapshot simulations carried out at 500-year intervals from 21 to 0 ka (utilized by Morris et al., 2018; shed some light on the combined effects of changes in ice sheet geometry, atmospheric greenhouse gases, and orbital parameters. Considering the independent effects of these changes (briefly discussed above), it is not surprising that the climate change produced by these snapshot simulations does not match the spatial or temporal patterns of climate and ocean circulation change recorded for Heinrich Stadial 1. For example, the snapshot simulations show AMOC strengthening and widespread surface warming at the onset of Heinrich Stadial 1, which is inconsistent with most records ( Figure S6). This supports our argument that ice sheet meltwater flux (which is not included in the snapshots) is required to simulate Heinrich Stadial 1. However, it is important not to overinterpret these results since the snapshot simulations are not transient and the addition of ice sheet meltwater is unlikely to produce a linear response Singarayer et al., 2011;Singarayer & Valdes, 2010). Coordinated simulations of the last deglaciation that include a full set of transient forcings and explore uncertainty in model boundary conditions will provide further insight on this issue (Ivanovic et al., 2016).

Ability of Models to Respond Abruptly to Freshwater Forcing
It has previously been suggested that the present model, as well as others used for similar studies (Kageyama et al., 2013;Z. Liu et al., 2009;Menviel et al., 2011), are intrinsically too insensitive to freshwater forcing to accurately simulate an abrupt change in ocean circulation (Valdes, 2011). Furthermore, there are open questions about AMOC stability, and whether moderate freshwater forcing pushed it into a stable "off" mode during Heinrich Stadial 1 (Rahmstorf, 1995(Rahmstorf, , 2002, whereby continual freshwater forcing was not needed to maintain a collapsed AMOC. A metric used to assess a simulation's intrinsic AMOC [bi-]stability (i.e., the ability of AMOC to persist in both a stable "on" and "off" mode) is the net import of freshwater to the Atlantic Ocean, termed F ov (as defined by Hawkins et al., 2011), where a positive F ov indicates that the AMOC is mono-stable and a negative Fov characterizes a bi-stable AMOC. We do not know what the early deglacial F ov was, since there are no proxies for this metric, but in our simulation, F ov is~0.4 Sv, suggesting that the AMOC regime is mono-stable and therefore cannot exist in an unforced collapsed state. It is worth noting that while F ov may be a reliable metric for the stability of present-day AMOC (W. Liu et al., 2013Liu et al., , 2017, it is not clear whether this is the case for glacial climate (Zhang et al., 2017;Figure S2 k). Moreover, to consider AMOC as existing either in a mono-stable or bi-stable state is likely a misleading oversimplification (as discussed by Ivanovic et al., 2018), due to there being multiple sites of deep water formation that can switch on/off. Since HadCM3 is able to simulate a Heinrich Stadial 1-like AMOC slowdown and climate change with realistic meltwater fluxes (given uncertainties in timing of the forcing, see sections 4 and 5), this implies that the model has adequate sensitivity to freshwater forcing. The simulated change in AMOC produced by the imposed meltwater may be different under alternative boundary conditions (e.g., ice sheet geometry) or if performed with another climate model, but the key conclusion from the presented results is that AMOC bi-stability is not necessary to simulate the onset of Heinrich Stadial 1. A longer integration and improved knowledge of ice sheet geometry and meltwater histories is required to establish whether the Heinrich Stadial 1 ocean circulation and climate would persist in the model under the simulated AMOC regime.

Connection Between Heinrich Stadial 1 and Heinrich Event 1
In our simulation, the Arctic and Nordic Seas have a weak subsurface warming of <1°C, which is less than was previously simulated using a different, larger, and earlier meltwater forcing (Marcott et al., 2011). However, according to recent research by Bassis et al. (2017), the simulated multicentury weak subsurface warming may still be sufficient to trigger a Heinrich Event (after a few centuries) through melting at the ice sheet's submerged grounding line. Evidence for such warming exists in the Nordic Seas (Rasmussen & Thomsen, 2004) and south of the Labrador Sea mouth (Marcott et al., 2011). That the warming does not extend much further north into the Labrador Sea in our simulation suggests a different trigger for heightened North American iceberg calving~16 ka, possibly lending weight to the "binge-purge" hypothesis (MacAyeal, 1993) for Heinrich Event 1, whereby internal glacial dynamics are responsible (Hodell et al., 2017). On the other hand, continued climate forcing beyond the acceleration in ice sheet melt examined here-for example, the input of more meltwater from the elevated level of Eurasian iceberg discharge (Benway et al., 2010;McManus et al., 1999;Peck et al., 2006)  weaken the AMOC to produce heightened and more widespread subsurface warming. This warming would then provide a positive feedback to enhanced iceberg calving during Heinrich Event 1 by basal melting of buttressing ice shelves (Alvarez-Solas et al., 2010). In this instance, the acceleration of northern ice sheet melt plays an important role in preconditioning the deglacial ocean for enhanced North American iceberg calving to begin (Peck et al., 2007).
Further to this, the reduction in SST simulated from~17 ka in this study, but which we hypothesize occurred earlier (see section 4), may have helped to create sustainable conditions for North American and Eurasian derived icebergs to reach further south before melting, as inferred for the subsequent Heinrich Event 1 (Barker et al., 2015;Bond & Lotti, 1995;Grousset et al., 1993;Hodell et al., 2017). As well as amplifying the subsurface warming, any reduction in the AMOC from the iceberg meltwater could also enhance the simulated surface climate changes, including the southward shift of the ITCZ, which would explain the larger anomalies in precipitation immediately prior to the Bølling warming (Deplazes et al., 2013;Muller et al., 2012). Thus, Heinrich Event 1 (and earlier) iceberg calving may have played an important positive feedback role in driving or at least sustaining the Heinrich Stadial 1 AMOC and climate change and could explain the peak in North Atlantic 230 Pa/ 231 Th 16-15 ka (Ng et al., 2018).

Conclusions
The precise cause of the AMOC slowdown and northern cooling that characterizes the onset of Heinrich Stadial 1 at~18.5 ka remains enigmatic. While previous work has attributed the recorded ocean circulation and climate changes to ice meltwater (associated with enhanced iceberg calving and/or terrestrial ice sheets), the amount of meltwater used in these studies has consistently been unrealistically large. In this study, we used the recent ICE-6G_C reconstruction of ice sheet evolution to derive a spatially detailed meltwater history and then used this history to drive the climate model. The small acceleration (from 1 to 1.8 m/ kyr GMSLRE) in Eurasian ice sheet melting that begins~19.0 ka in our simulation (~18.5 ka in ICE-6G_C) has negligible impact on Atlantic overturning or surface climate. On the other hand, a larger acceleration in northern ice sheet melting (from 4 to 8 m/kyr À1 GMSLRE) that occurs later in ICE-6G_C (after 17 ka in our simulation) sufficiently stabilizes the north Atlantic water column to inhibit NADW formation. This reduces AMOC strength by~20%, inducing a weak bipolar climate signal with widespread and relatively strong Northern Hemisphere cooling (enhanced in the winter due to expanded sea ice) and less pronounced Southern Hemisphere warming. The result is a complex and coherent pattern of climate change that matches paleoclimate and ocean circulation reconstructions for the onset of Heinrich Stadial 1, in terms of both rates and amplitude of change.
We argue that based on a range of paleo-ice sheet, paleo-drainage, and paleoclimate data (Bard et al., 2000;Buizert et al., 2014;Calvo et al., 2001;Chen et al., 2016;Hughes et al., 2016;Mojtahid et al., 2017;Naafs et al., 2013;Ng et al., 2018;Toucanne et al., 2010Toucanne et al., , 2015, the ICE-6G_C timing of melting in the earliest stages of the last deglaciation is biased toward being too young. Therefore, we propose that an increase in northern ice sheet melt of the order of 4 m/kyr GMSLRE was needed to trigger the ocean circulation and surface climate changes that took place at the onset of Heinrich Stadial 1~18.5 ka. This acceleration from~3.5 to~7.5 m/kyr GMSLRE is less than half of the 15-22 m/kyr GMSLRE required in previous studies to drive Heinrich Stadial 1 (Z. Liu et al., 2009;Menviel et al., 2011) and fits existing constraints on sea level and ice sheet evolution.
Moreover, accelerated Eurasian ice sheet melting has also been recorded for earlier Heinrich Stadials (e.g., 2 and 3,~31-29 ka and~26-23 ka, respectively), caused by retreat of the southern Fennoscandian ice sheet, which led to amplified freshwater discharge by the Channel River to the eastern North Atlantic (Toucanne et al., 2015). While it is not possible to reconstruct the ice sheet history in as much detail prior to the Last Glacial Maximum as after it, this does provide some indication that Eurasian meltwater could be a consistent driver of millennial-scale AMOC weakening (Böhm et al., 2015;Henry et al., 2016), North Atlantic surface cooling (Barker et al., 2015), and subsurface warming throughout the late Quaternary (Barker et al., 2015;Marcott et al., 2011), with particularly strong potential for forcing these changes during deglacial phases. This chain of events may even have acted as a precursor or trigger for enhanced iceberg calving marked by the Heinrich Events (Bond & Lotti, 1995;Grousset et al., 2000;Scourse et al., 2000;Shaffer et al., 2004), as presented in section 8.

10.1029/2017PA003308
Of wider and more pressing relevance, the rates of Northern Hemisphere deglaciation used to drive our simulation are similar to projected rates of Greenland melt by 2200 (6.7 m/kyr GMSLRE by 2200; Lenaerts et al., 2015). Since Greenland meltwater enters sites of present-day NADW formation relatively quickly and directly, understanding the trigger for Heinrich Stadial 1 will provide important constraints for climate sensitivity to future ice sheet melt and test our models' ability to correctly simulate these changes.