The potential of Hudson Valley glacial floods to drive abrupt climate change

The periodic input of meltwater into the ocean from a retreating Laurentide Ice Sheet is often hypothesized to have weakened the Atlantic meridional overturning circulation (AMOC) and triggered several cold periods during the last deglaciation (21,000 to 8,000 years before present). Here, we use a numerical model to investigate whether the Intra-Allerød Cold Period was triggered by the drainage of Glacial Lake Iroquois, ~13,300 years ago. Performing a large suite of experiments with various combinations of single and successive, short (1 month) and long (1 year) duration flood events, we were unable to find any significant weakening of the AMOC. This result suggests that although the Hudson Valley floods occurred close to the beginning of the Intra-Allerød Cold Period, they were unlikely the sole cause. Our results have implications for re-evaluating the relationship of meltwater flood events (past and future) to periods of climatic cooling, particularly with regards to flood input location, volume, frequency, and duration. Drainage of Glacial Lake Iroquois via the Hudson Valley about 13,300 years ago did not weaken the Atlantic Meridional Overturning Circulation substantially, suggests a suite of simulations with a global-ocean general circulation model.

I n the late 1980s periodic meltwater inputs from a retreating Laurentide Ice Sheet (LIS) were hypothesized to have hindered North Atlantic deep water (NADW) formation, leading to a weakened Atlantic meridional overturning circulation (AMOC), and ultimately triggering periods of cooling in the North Atlantic over the course of the most recent deglaciation (21-8 ka yrs BP) 1 .
In the intervening 30 years, geomorphic evidence of many meltwater discharge events from the LIS has been uncovered and their timing shows broad coincidence with centennial-to-millennial length cool periods (e.g., Younger Dryas, the 8.2-kyr-event). However, despite this temporal correspondence, uncertainty surrounding meltwater flood characteristics (e.g., volume, duration) have made it difficult to assign a specific meltwater event as the trigger for a specific cold period or that each flood event even elicits an AMOC response. Recent modeling also suggests that the location and duration of meltwater input significantly influence the extent to which AMOC is impacted 2,3 . These realizations highlight the need to better understand the impact of meltwater inputs on ocean circulation and the climate system, especially as ice loss from Greenland and Antarctica continues to accelerate, adding significant quantities of freshwater to the oceans [4][5][6] . Here, we perform a suite of numerical model simulations to test the hypothesis that meltwater floods originating from the Hudson Valley triggered the Intra-Allerød Cold Period (IACP) 7 .
Glacial meltwater floods and the IACP The IACP is a cold climate period occurring towards the end of the Bølling-Allerød interstadial, just prior to the onset of Younger Dryas (YD) cooling and is defined in synchronized Greenland ice core records as the negative isotope excursion between 13,311 and 13,099 yr b2k (before 2000 CE) 8 . At this time, Glacial Lake Iroquois was impounded by the retreating LIS in the modern St. Lawrence valley, between the US state of New York and the Canadian province of Ontario (Fig. 1). Continued northward retreat of the LIS during deglaciation is thought to have eventually allowed Lake Iroquois to catastrophically drain into Glacial Lake Coveville in the Champlain Valley, creating a combined 700 km 3 flood down the Hudson Valley 7, [9][10][11] . Current evidence supports a second larger flood (~2500 km 3 ) then drained the successor to Lake Iroquois, Lake Frontenac, along with the glacial lake in the Champlain Valley (Lake Fort Ann) down the Hudson River 9-12 ( Fig. 1). Available estimates indicate that discharge from the Ontario -St. Lawrence lowlands reached at least 0.1 Sv during the drainage of Lake Iroquois with sustained outflows of at least 0.06 Sv during the existence of Lake Frontenac 10 . Absolute and chronostratigraphic dating of geological sequences associated with the initial outburst of Glacial Lake Iroquois constrains the timing of the event to~13,300-13,400 cal BP 7 . Additional available temporal constraints from paleoshorelines and flood deposits place the two floods both between 13,050-13,310 cal BP 13 . During this period there is also evidence of a suppressed AMOC coinciding with the IACP from reconstructions of subtropical North Atlantic intermediate-water temperature and benthic foraminifera 14 C ventilation ages in the Nordic Seas 14,15 . Therefore, the location, timing, and estimated magnitude of these flood events has thus led to the hypothesis that meltwater emanating from the Hudson Valley weakened the AMOC and triggered the cooling associated with the IACP 7,16 .

Assessing meltwater impacts on AMOC
Here we assess the sensitivity of the AMOC to a range of meltwater discharge scenarios from the Hudson Valley using a highresolution (1/6°; 18-km) eddy-resolving global configuration of the coupled Massachusetts Institute of Technology general circulation model (MITgcm 17 ), with sea-level 75 m lower to approximate conditions 13,300 years ago 18 (Fig. 1). Modern atmospheric boundary conditions are used to approximate the relative warmth of the Allerød period (see methods). The experiments performed were specifically designed to cover a range of possible meltwater scenarios and uncertainties related to the reoccurrence interval (frequency), duration, and volumes of the flood events~13,300 years ago. Initially, we ran an experiment with two, 1-month duration floods of 0.31 Sv (i.e. 700 km 3 ) and 1.09 Sv (i.e. 2500 km 3 ) spaced 10 years apart to produce a total freshwater discharge of~3700 km 3 to reflect the two floods emanating from the Hudson Valley around the time of the   10,12,13 . This is referred to as the "Realistic" scenario. In all remaining simulations, the total discharge of each flood event used a volume that combines both flood events and includes an estimated maximum uncertainty of %15 [10] (i.e., [700 + 2500] x 15% =~3700 km 3 ). To evaluate the effect of flood duration on AMOC, this volume of water was released once (at the start of the model run) over a period of either 1 month ("1 mo") or 1 year ("1 yr"), to create single meltwater flood events with fluxes of 1.40 and 0.12 Sv, respectively. To assess the impact of flood frequency on AMOC, the second set of experiments were performed with the same volume of water (3700 km 3 ) released in three successive meltwater floods again with durations of either 1-month ("3x_1mo") or 1-year ("3x_1yr") with a return interval (i.e., frequency) of 2 years. Finally, to assess the role of any background runoff in preconditioning the AMOC to weaken or collapse from additional glacial meltwater input, an identical set of flood experiments were performed with 0.05 Sv of meltwater also continuously released from the Hudson River 10,19 (Table 1).

Results
In all experiments (regardless of flood magnitude, duration, reoccurrence interval, or background flux), meltwater released from the Hudson Valley initially leads to a sea surface freshening near the mouth of the paleo Hudson valley of~1.5 PSU. The entrainment of this meltwater in the near-shore shelf region into the Gulf Stream then leads to freshwater being advected in a northeast direction towards the subpolar North Atlantic where processes important for modulating the strength of the AMOC occur ( Fig. 2a, b). Despite this, and importantly from a climatic perspective, we find that the fate of the meltwater in all of our experimental scenarios is similar: regardless of flood duration or flood reoccurrence interval, turbulent mixing within the Gulf Stream results in meltwater from the Hudson River being rapidly mixed into the interior ocean so that by the time the freshwater is advected to the regions of NADW formation in the Labrador and Greenland Seas it is no longer detectable (Supplementary Figs. 1-5). Indeed, Fig. 2c confirms that there is no significant reduction in surface salinity in these regions over time, compared to the control simulation. Even under the most extreme scenario in which three, 1-month duration floods were released successively (3x_1mo), the freshwater, while able to persist for longer and expand farther north away from the mouth of the paleo Hudson River, rapidly became mixed into the ambient ocean after each flood event (Fig. 2, see Supplementary Movie 1). This rapid dissipation of the meltwater prevents the establishment of any lower-density surface layer in the subpolar gyre that might otherwise have been capable of inhibiting "i.e., capping" NADW formation ( Supplementary Fig. 6). While some small sea surface salinity (SSS) anomalies appear in the Greenland Sea, these represent local circulation variability compared to the Control. The lack of any change in salinity and reduction in NADW formation is further confirmed by calculating mixed layer depth-a model proxy for deep water formation (see Methods)-in the Labrador and Greenland Seas. In all of our experiments, regardless of flood duration, reoccurrence interval, or background flux, the mean mixed layer depth varied from the Control experiment by less than 4% as a result of the imposed flood perturbations ( Table 2, Supplementary Fig. 7, and Table 1). Furthermore, we find that in all of our experiments mean AMOC strength deviates by <1% from the Control simulation and exhibits a standard deviation similar to modern-day variability 20 (Fig. 3 and Table 2). This result also holds for the experiments run without the background flux of 0.05 Sv from the Hudson Valley (Supplementary Fig. 8 and Table 1).

Discussion
There appear to be several aspects of the Hudson Valley meltwater flood events that prevent them from influencing the strength of the AMOC, as well as experimental design considerations that warrant further discussion in light of the above results. To begin, our simulations show that, in all cases, meltwater is rapidly entrained, and diffused by, the Gulf Stream, causing any meltwater to sufficiently mix into the interior ocean before it reaches the subpolar gyre and locations of overturning circulation ( Fig. 2 and Supplementary Movie 1). The length of the transport distance between the Hudson Valley input location and the Labrador and Greenland Seas likely plays an overarching role in the dissipation of the freshwater anomaly. Additionally, the location and strength of the Gulf Stream plays an important role in the transportation of the meltwaters after exiting the Hudson Valley. In our simulations, the energetic and turbulent nature of the Gulf Stream combined with the relatively small volumes of water released from the Hudson River mean that any fresh water is mixed into the ocean extremely quickly and very far south of any potential regionals of NADW formation. This reinforces the results of previous modeling studies that showed the eventual impact of meltwater events on overturning strength is largely controlled by the location of meltwater input and the downstream flow paths to areas of deep water convection 3,21 .
Combined with the long transport distance, the rapid dissipation of meltwater in the turbulent Gulf Stream also echoes previous work showing that meltwater inputs can be more rapidly mixed and dissipated in higher-resolution eddy-resolving models than in lower-resolution, non-eddy-resolving models, limiting the transport of freshwater anomalies and leading to smaller reductions in AMOC 3 . Conversely, some eddy-resolving models have shown more persistent AMOC reductions compared with lowerresolution simulations 22,23 , however, these simulations use meltwater input locations from Hudson Bay, Canada, or coastal Greenland, which have significantly shorter and different flow paths to areas of the modern overturning circulation. In all likelihood, the input location, longer transport distance, and diffusion in the presence of eddies are realistic aspects of the floods presented here that play a role in limiting the transport of meltwater to locations of NADW formation. It should be noted that even in the "Realistic" scenario, a faint SSS anomaly can be seen migrating north along the east coast of modern Newfoundland (Fig. 2). This indicates that transport of meltwaters into the Labrador Sea is likely possible in our model conditions, but even our "3x_1mo" extreme scenarios are incapable of transporting sufficient freshwater into the Labrador Sea to weaken NADW formation (Fig. 2) and the strength of AMOC remains unaffected by the meltwater floods emanating from the paleo Hudson Valley.  Another aspect explored by the above simulations is the duration of individual flood events. Previous geomorphic studies of the Hudson Valley paleo-channel indicate that the largest discharge events would produce floods lasting approximately 3 months 16 . Our use of both 1-month and 1-year flood durations explores this range of possible flood event lengths and shows that neither impacts AMOC in any appreciable manner. Furthermore, the simulations testing the impact of flood frequency with combined floods repeating on a 2-year interval (3x_1mo and 3x_1yr) do not produce a salinity anomaly in NADW formation locations or provoke an AMOC response despite releasing three times more freshwater (11,100 km 3 ) into the North Atlantic than the "single" flood experiments. Of significance is that this volume of meltwater is comparable to the~9500 km 3 freshwater discharge from lake Agassiz often invoked as the trigger for the millennial length YD 24 , yet there is still no reduction in AMOC strength similar to that expected for this particular large-scale cooling event.
The meltwater flood experiments provide context for the 0.1 Sv meltwater flux suggested by many modeling studies as the threshold for slowing the AMOC. In these studies a constant 0.1 Sv freshwater flux into the North Atlantic is needed to both weaken the AMOC and for it to remain suppressed 25 ; once the freshwater forcing is removed (turned off) these models quickly show a resumption of AMOC 16,[25][26][27] . As glacial outburst floods were either single events or at most a few events spaced several years apart, the flux of meltwater to the ocean would have only exceeded 0.1 Sv for a year or so at most, rather than being a continual centennial-to-millennial forcing. Other studies exploring the AMOC response to a Hudson Bay meltwater input associated with the 8.2 ka event also found that despite using a greater flood volume and duration (2.5 Sv for 1 year), AMOC reduction could not be sustained for the observed duration of the 8.2 ka event (~150 years) 28,29 . The lack of AMOC response in our model simulations reasserts that the additional meltwater supplied to the North Atlantic by outburst floods emanating from glacial lakes was unlikely to have been large enough to weaken AMOC or invoke a large-scale climate cooling. A similar finding was also made by Meissner and Clark 30 .
In our freshwater forcing simulations investigating the role of background runoff in preconditioning the AMOC to collapse in response to short-lived glacial outburst floods, the AMOC strength (with and without flood events) was found to be~1-2% weaker than in simulations with no background discharge (Supplementary Fig. 8 and Table 1). Significantly, however, we were again not able to find a prolonged weakening in the overturning cell in any of our experiments releasing either single or multiple outburst floods on top of this background flux that suggest a threshold in the climate system capable of causing the prolonged IACP cooling had occurred.
It is worth noting that our choice of background flux is only half the aforementioned 0.1 Sv flux required in many models to significantly weaken AMOC but this is constrained by geomorphic evidence suggesting that higher background discharge from this region was unlikely 10,19 . We thus speculate that background drainage from other LIS outlets, such as the Gulf of St. Lawrence, Hudson Bay, Mackenzie River, might have played an additional role in further "pre-conditioning" the overturning cell so that smaller floods emanating from the Hudson Valley may have been sufficient to significantly weaken AMOC. However, as it stands, our results indicate that meltwater emanating solely from the Hudson Valley~13,300 BP (in the form of either outburst floods or persistent drainage from the LIS) is unlikely to have been the main cause of the cooling associated with the IACP. Even in the most extreme scenario (which likely far exceeds the actual flood events at this time), the meltwater events appear either too small, too infrequent, or originate too far south of regions of NADW formation to impact AMOC and trigger the IACP.
Atmospheric conditions should also be considered in light of our simulations. Recent modeling work has suggested that changes in atmospheric circulation due to the presence of the LIS have the potential to set up and trigger abrupt climate shifts 31,32 . The prescribed, uncoupled atmosphere in our simulations leaves out this potential "pre-conditioning" mechanism but allows us to test whether or not the Hudson Valley floods could have been responsible for significant AMOC weakening and instigation of the IACP. Having established the unlikely nature of the Hudson Valley floods as the sole driver of a weakened AMOC, future experiments can explore the interactions between competing mechanisms and potential AMOC and climate responses.
Future explorations should also investigate how separate meltwater flood events from multiple outlets may have preconditioned the AMOC to weaken in response to additional meltwater input to better understand the complexities and sensitivities of AMOC during deglaciation. Recent studies have pointed towards longer periods of meltwater input from accelerated ice sheet collapse or the combined effects of a meltwater input and accompany ice discharge from retreating ice sheets as a mechanism for driving the magnitude and duration of AMOC reduction during abrupt climate shifts 28,29 . Furthermore, recent modeling studies have shown that the combination of sea level rise, changes in atmospheric circulation, and sea-ice export, along with glacial flood inputs can significantly weaken AMOC 33 . Additionally, meltwater inputs (or a background flux) have the potential to lower SSTs and perhaps change the mixed layer depth 28,33 . However, our Control simulation with a constant 0.05 Sv flux from the Hudson Valley shows similar SSTs to our Control simulation without 0.05 Sv flux (Supplementary Fig. 9). SST changes due to flood events are also short-lived and localized to the input area and do not impact the position of the sea ice edge or mixed layer depth (Supplementary Fig. 9). Finally, in concert with assessing more complex AMOC forcings, explorations of the contribution or role of other mechanisms such as volcanism and sea-ice feedbacks in abrupt climate shifts should be considered 34,35 . Understanding the full impact of meltwater/freshwater inputs on the North Atlantic climate is all the more important as river discharge and melting of the Greenland and Antarctic Ice sheets is projected to increase and drive complex changes across the Northern Hemisphere 6,36,37 . Inclusive of flood duration, frequency, and volume, this study highlights the importance of the geographical location of any meltwater input on the potential to alter large-scale ocean circulation and climate. Finally, future model studies should be aimed at exploring the spatial and temporal interactions of multiple meltwater/freshwater sources on AMOC as well as the role of other associated and linked climate mechanisms.

Methods
All numerical model simulations were performed using the Massachusetts Institute of Technology General Circulation Model (MITgcm) 17 . The model has an eddypermitting horizontal global grid resolution of 1/6°(~18 km) and 50-levels in the vertical with spacing set from~10 m in the near-surface to~450 m at a depth of 6000 m. Ocean tracer transport equations are solved using a seventh-order monotonicity preserving advection scheme. There is no explicit horizontal diffusion, and vertical mixing follows the K-Profile Parameterization. The ocean model is coupled to a dynamic-thermodynamic sea ice model that assumes viscous-plastic ice rheology and computes ice thickness, ice concentration, and snow cover 38 . Atmospheric boundary conditions were taken from modern (1979-2002 monthly means) ERA-40 reanalysis data from the European Centre for Medium-range Weather Forecasts. Without access to atmospheric boundary conditions for the specific IACP interval, modern atmospheric conditions provide a reasonable approximation of the relatively warm conditions of the overarching Bølling-Allerød. The configuration follows that of Condron and Winsor 3 but with sea-level set 75 m lower which results in the Bering Strait being closed. Control simulations for meltwater simulations with and without a 0.05 Sv background meltwater flux 19 were spun up for 20 years. In all experiments, the glacial meltwater has a temperature and salinity of 0°C and 0 PSU. Geomorphic evidence indicates that floodwaters flowed southward through the breached moraine dam near the modern-day Narrows-Verrazzano Bridge and to the southeast along the modern continental shelf 16,39 . Therefore, floodwaters were released into the four model grid cells closest to the intersection of the 75 m bathymetric contour and the modern-day undersea Hudson canyon. The Labrador and Greenland Seas are defined as the regions between 52°N to 65°N, 70°W to 20°W, and 70°N to 80°N, 20°W to 20°E, respectively. The 37-43°N region encapsulating the outflow from the paleo-Hudson channel is bounded by longitudes of 7-74°W. Mixed layer depth is calculated at each grid cell where water depth is >2000 m and is defined as the depth at which the density of water 40 exceeds the surface density by 0.125 kg m −3 .

Data availability
Bathymetric data utilized in model runs are sourced from the General Bathymetric Chart of the Oceans (GEBCO) gridded dataset available at https://www.gebco.net/ data_and_products/gridded_bathymetry_data/. ERA-40 monthly atmospheric data is available through the European Centre for Medium-Range Weather Forecasts (ECMWF) at https://apps.ecmwf.int/datasets/data/era40-moda/levtype=sfc/.

Code availability
The model code (MITgcm) utilized to perform all of the experiments and generate the figures is publicly available from mitgcm.org.