Research ArticleCLIMATOLOGY

The role of Northeast Pacific meltwater events in deglacial climate change

See allHide authors and affiliations

Science Advances  26 Feb 2020:
Vol. 6, no. 9, eaay2915
DOI: 10.1126/sciadv.aay2915


Columbia River megafloods occurred repeatedly during the last deglaciation, but the impacts of this fresh water on Pacific hydrography are largely unknown. To reconstruct changes in ocean circulation during this period, we used a numerical model to simulate the flow trajectory of Columbia River megafloods and compiled records of sea surface temperature, paleo-salinity, and deep-water radiocarbon from marine sediment cores in the Northeast Pacific. The North Pacific sea surface cooled and freshened during the early deglacial (19.0-16.5 ka) and Younger Dryas (12.9-11.7 ka) intervals, coincident with the appearance of subsurface water masses depleted in radiocarbon relative to the sea surface. We infer that Pacific meltwater fluxes contributed to net Northern Hemisphere cooling prior to North Atlantic Heinrich Events, and again during the Younger Dryas stadial. Abrupt warming in the Northeast Pacific similarly contributed to hemispheric warming during the Bølling and Holocene transitions. These findings underscore the importance of changes in North Pacific freshwater fluxes and circulation in deglacial climate events.


Abrupt Northern Hemisphere climate fluctuations during the last glacial and deglacial periods are widely assumed to track changes in freshwater fluxes and sea-ice dynamics of the high-latitude North Atlantic (15); buoyancy forcing provides a mechanism for cooling via suppression of deepwater formation and reduction of net northward oceanic heat transport (1, 2). Deglacial warming is often attributed to rising CO2 (6), but this cannot account for the abrupt onset of Bølling warming at 14.7 thousand years (ka) (7). An “overshoot” of the Atlantic meridional overturning circulation (AMOC), forced by rapid reduction of local freshwater input, has been hypothesized (7), but evidence for such an AMOC overshoot is equivocal (5, 8), and sea surface temperature (SST) reconstructions near the inferred site of deepwater formation show muted warming during the Bølling event (9). Moreover, the warming overshoot forced by abrupt removal of fresh water in the model (7) is not consistent with reconstructed freshwater fluxes (10). Thus, the case for forcing of North Atlantic changes based on local meltwater remains uncertain.

So far, changes in the North Pacific have received less attention as a forcing agent of abrupt climate fluctuations in the past, in part due to the lack of deepwater formation (11) and a dearth of high-resolution records from the region. Nevertheless, models show that perturbations to the surface ocean in the North Pacific can have substantial climate impacts on Northern Hemisphere climate through changes in sea-ice extent and poleward moisture transport (2, 12, 13).

A key controversy is whether a deep or intermediate water mass in the Northeast Pacific developed during times when the strength of the AMOC was reduced, such as during Heinrich Stadial 1 (HS1; ~18 to 15 ka) (14, 15) or the Younger Dryas (YD; 12.9 to 11.7 ka), and how changes in Northeast Pacific ventilation relate to surface ocean changes. The development of an oceanic thermal seesaw between the North Pacific and North Atlantic (16) would tend to buffer Northern Hemisphere temperature changes; reduced heat transport in one basin would be partly compensated by the other basin (14, 17). In some models, enhanced Pacific meridional overturning circulation warms the Northeast Pacific sea surface and increases surface salinity during reductions in the strength of AMOC (14, 15). These models generally neglect deglacial meltwater inputs from the retreating Cordilleran Ice Sheet (18).

Cordilleran meltwater inputs came from glacial calving and runoff locally around the Gulf of Alaska (GOA), into the Bering Sea via the Yukon River drainage, and from glaciers on Vancouver Island, as well as from the western North American interior through the Fraser and Columbia River basins. There were also a series of late glacial and deglacial megafloods drained via the Columbia, funneling 0.5 to 17 sverdrups (106 m3 s−1) of fresh water over a period of days to weeks (1924). Flood repeat times are uncertain, but the largest events from ice-dammed Lake Missoula occurred at least 89 times in a few thousand years (25), implying decadal-to-multidecadal frequency. The primary flooding interval is between ~19 and ~14 ka (24), but other events of unknown origin exist, including one near ~23 ka (21), as well as later events that extend into the YD (1926). Evidence for cyclical ~80-year ice-rafting and meltwater plume events between 19 and 17 ka also occurs in marine sediments near Vancouver Island, possibly associated with outburst floods (27).

Sustained freshwater inputs into the Northeast Pacific from a melting Cordilleran Ice Sheet are estimated at 0.02 to 0.06 sverdrups from 20 to 12 ka (18), although short-term variations in flux could be larger. Closure of the Bering Strait when sea level was lower would have impounded an additional 0.03 sverdrups of fresh water in the Northeast Pacific based on estimates of modern freshwater loss through the Bering Strait (11). Transient outburst flood events into the Pacific superimposed on this baseline freshening could have plausibly exceeded 0.1 sverdrups on decadal-to-centennial time scales, large enough to elicit regional SST decreases (28) and widespread Northern Hemisphere cooling and sea-ice expansion (2). For comparison, sustained long-term flows of Laurentide meltwater to the North Atlantic are thought to have been 0.05 to 0.2 sverdrups (10), and glacial outburst floods to the North Atlantic may have included 5 sverdrups of freshwater flows over a half year (29); i.e., although the long-term flows to the North Pacific are likely smaller than those to the Atlantic, peak Pacific flooding events may have been as large or larger.

Here, we model the flow trajectory of fresh water exiting the Columbia River under both glacial and deglacial climate boundary conditions, with the Bering Strait closed and open, respectively, to assess the likely pathways that meltwater takes once injected into the Northeast Pacific (Fig. 1 and Materials and Methods). We compile reconstructions of the alkenone-based Uk′37 index for SST, records of near-surface paleosalinity based on the δ18Oseawater−ice volume corrected18Osw-ivc) derived from planktic foraminifera δ18O data, and records of water column reservoir ages based on paired planktic and benthic foraminifera radiocarbon ages in marine sediment cores from 40° to 60°N in the Northeast Pacific (Fig. 1), which span the last deglaciation. We then compare our regional Northeast Pacific SST compilation with a global SST compilation to reveal SST patterns during major deglacial climate intervals that help to identify plausible forcing mechanisms of climate variability.

Fig. 1 Modeled flood from the Columbia River in the glacial simulation, along with the locations of paleoclimate datasets shown/referenced in this study.

SSS anomaly relative to the control is shown after 1 year. Marine sediment cores include EW0408-85JC, EW0408-87JC, EW0408-66JC, EW0408-26JC, ODP 887, and SO2020-27-6 from the GOA; JT96-09PC and MD02-2496 from the Vancouver margin; and ODP 1019 and W8709-13PC from the California margin, along with locations of the speleothem records from the Oregon Cave National Monument and the Cave of Bells.

We have four aims: (i) to develop a modeling framework for understanding the fate and far-field impacts of increased freshwater fluxes in the North Pacific, (ii) to reconstruct deglacial changes in sea surface salinity (SSS) in the Northeast Pacific, (iii) to determine water column radiocarbon age changes in the Northeast Pacific to evaluate whether changes in upper ocean density are related to changes in subsurface ventilation, and (iv) to synthesize Northeast Pacific deglacial SST changes and their relationship to global climate patterns.


Floods and meltwater routing

In the glacial simulation, with a closed Bering Strait, fresh water released from the Columbia River (see Materials and Methods for details) is primarily routed northward into the cyclonic Alaskan Coastal Current, where it reaches all the marine core locations in our study within the first year (Fig. 1). Some southward transport is also observed, consistent with inferences of large salinity anomalies [4 practical salinity units (psu)] recorded in marine cores at 42°N as evidence of Columbia River flood events (21). From the GOA, the fresh water is transported into the Northwest Pacific along boundary currents and subsequently entrained into the North Pacific current, which recirculates it back to the eastern margin, resulting in dispersal of a lower magnitude SSS anomaly across the entire subpolar gyre (Fig. 2). In the model, this generates an SST cooling in the Kuroshio Current region (2°C), as well as modest cooling (0.5°C) along the North Pacific margins.

Fig. 2 Modeled trajectory of a flood from the Columbia River in the glacial and deglacial simulations, with a closed and open Bering Strait, respectively.

SSS anomalies are shown for model years 5 and 10 for each case, along with the SST anomalies at year 5 for the closed Bering Strait simulation and year 15 for the open Bering Strait simulation. The Kuroshio Current cools in the closed Bering Strait simulation, whereas the Gulf Stream and subpolar North Atlantic cool in response to freshwater migration to this region in the open Bering Strait simulation.

In the simulation with an open Bering Strait (see Materials and Methods for details), a substantial fraction of the injected fresh water is routed into the Arctic and then travels southward along the western boundary currents of the North Atlantic (Fig. 2), where it reaches the Gulf Stream and Atlantic deepwater formation regions of the Labrador Sea within a decade. This results in an export of Arctic sea ice southward into the North Atlantic (fig. S1) and SST cooling in the Gulf Stream and Labrador Sea regions (Fig. 2).

These results raise the question as to whether the Bering Strait was open during any of the late deglacial flood or meltwater events from the Cordilleran region. While most of the large Columbia River megaflood events are thought to have occurred early in the deglaciation (19 to 14 ka) (24), evidence for flood events through the Columbia extends into the YD period (19, 22, 26). Our δ18Osw reconstructions indicate freshening along the Northeast Pacific margin during the YD (Fig. 3). Dates for the opening of the Bering Strait have been controversial, with the earliest estimates near the start of the YD (~13 ka) (30) and later estimates clustering around the early Holocene (~11 ka) (31). Pico et al. (32) show that the Bering Strait was partially open but much shallower than today between 13 and 11.5 ka because of the interplay of global sea level rise and glacial isostatic adjustment. This overlap creates the possibility that YD-age meltwater events into the Pacific may have been partially routed into the Arctic and North Atlantic deepwater formation regions.

Fig. 3 Records of Northeast Pacific alkenone-derived Uk′37 SST reconstructions compared with Greenland ice core records and speleothem records from the western United States spanning the deglaciation and Holocene (left) along with a synthesis of Northeast Pacific paleosalinity records (right).

(Left) Records from top to bottom: δ18O records from NEEM [dark green; (46)] and NGRIP [light green; (53)] ice cores; Uk′37 records from EW0408-85JC [dark blue; (33)], EW0408-87JC (teal; this paper), EW0408-66JC [light purple; (36)] and EW0408-26JC [purple; (36)], JT96-09PC [magenta; (34)], MD02-2496 [pink; (38)], and ODP 1019 [blue; (35, 37)]; an average of the Northeast (NE) Pacific SST records for the margin sites with a 400-year running average (black); error bars are SEM; speleothem δ18O records from the Oregon Caves National Monument [bright pink; (69, 70)] and Cave of Bells [maroon; (71)]. (Right) δ18Osw-ivc records from EW0408-85JC [dark blue; (33)], EW0408-87JC (teal; this paper), EW0408-66JC [light purple; (36)], EW0408-26JC [purple; (36)], JT96-09PC (magenta; this paper), MD02-2496 [pink; (38)], and ODP 1019 (blue; this paper); an average (standardized, 400-year running average) of the δ18Osw-ivc records for the margin sites (errors bars are SEM), with blue shading denoting negative anomalies (low relative salinity) and maroon shading denoting positive anomalies (higher relative salinity); a reconstruction of surface salinity based on freshwater diatom assemblages from core W8709-13PC [black; (21)] and modeled freshwater discharge into the Pacific [100-year mean, 13P ensemble, dark green; (18)]. Reconstructions of the δ18Osw-ivc are based on paired SST and planktic G. bulloides δ18O records, with the exception of ODP site 1019, which is calculated on the basis of the thermocline-dwelling Nps δ18O record. Direct sample pairs were rare between datasets (plotted as data points), so continuous time series were calculated on the basis of 100-year interpolated records of SST and δ18O and then corrected for the δ18O of changes in global ice volume.

Paleotemperature and paleosalinity reconstructions

The Uk′37 SST reconstructions for the Northeast Pacific (Fig. 3) indicate that SSTs were, on average, 4° to 5°C colder during the late glacial (18 to 16.5 ka) relative to early Holocene values (11.5 to 9 ka); minima in SST occur near 17.5 ± 0.5 ka. Gradual warming of ~1°C occurred from 16.5 to 14.8 ka and was followed by an abrupt warming of 3° to 5°C at the Bølling transition from ~14.8 to 14.6 ka in the Northeast Pacific margin sites.

Most of the Northeast Pacific SST records show two warm intervals during the Bølling-Allerød (BA) punctuated by cooling from 14.2 to 14.0 ka. A 3° to 4°C cooling into the YD occurs in all cores except the GOA gyre site EW0408-87JC, which has relatively low SSTs during the BA (see Materials and Methods). An abrupt ~5°C warming into the Holocene is seen in cores EW0408-85JC (33), JT96-09PC (34), and ODP (Ocean Drilling Program) 1019 (35). The magnitude of the abrupt warming into the Holocene is slightly reduced in EW0408-66JC (36) relative to the other margin sites, and core EW0408-87JC in the Pacific subpolar gyre shows a more gradual warming into the Holocene that peaks at ~10 ka. Most cores show peak early Holocene warmth until ~9 ka, when a ~1°C cooling trend is initiated, reaching a Holocene minimum near ~6 ka, followed by warming near ~4 ka. The SST records from the Northeast Pacific show no evidence for a long-term cooling through the entire Holocene, as has been inferred from the North Atlantic region (37).

The δ18Osw-ivc values are calculated by removing the Uk′37 temperature effect and global δ18O changes related to ice volume from the planktic foraminiferal δ18O measurements (fig. S2); it represents a δ18O anomaly relative to a modern surface water and is roughly proportional to salinity anomalies if the δ18O of the freshwater source is constant. The δ18Osw-ivc values show greater spatial variability between sites than the SST records, especially during the early deglaciation (Fig. 3). For example, ODP site 1019 (~400 km south of the Columbia River off southern Oregon in the California Current) has relatively low δ18Osw-ivc values throughout the Last Glacial Maximum (LGM) and early deglaciation (Fig. 3), consistent with our modeled flow trajectories and previous interpretations of some southward transport of Columbia River outflow at the time of megaflood events (19, 21). A reconstruction of surface salinity based on freshwater diatom assemblages along the Oregon margin shows a prominent episode of surface freshening at ~18 ± 1 ka, inferred to reflect freshening from Columbia River floods (21), and these dates have been confirmed by cosmogenic exposure dates of flood-transported boulders on land (24). A δ18Osw-ivc record from the Vancouver margin (MD02-2496) (38) shows similar trends to the ODP 1019 δ18Osw-ivc record, with relatively fresh near-surface salinities throughout the early deglacial period. In contrast, core sites in the GOA show more distinct “pulses” of freshening that punctuate this interval (18 to 16 ka), likely reflecting their proximal position relative to local meltwater sources along the Southeast Alaska margin. For example, a prominent low δ18Osw-ivc event occurs in a GOA margin site (EW0408-26JC) at 17.5 ka (36), and two sites in the GOA gyre (EW0408-87JC and SO202-27-6 (28) indicate freshening events at ~16 ka.

To estimate regionally representative salinity variations, an average of all these records provides an approximation of deglacial changes in surface salinity (Fig. 3). This record exhibits a series of low δ18Osw-ivc (lower salinity) intervals throughout the deglaciation, with the lowest δ18Osw-ivc values of the last 20 ka occurring during the early deglaciation (19.0 to 16.5 ka) and YD periods (13 to 11.5 ka). The timing of generally low δ18Osw-ivc in our compilation during the early deglaciation is consistent with the timing of the main phase of the Columbia River megafloods sourced from Lake Missoula (24).

The BA period exhibits a large degree of variability among the paleosalinity records; the δ18Osw-ivc records of the Oregon (ODP 1019) and Vancouver [MD02-2496 (38) and JT96-09PC] margin sites and a diatom-based δ18Osw record from the GOA (SO202-27-6) (28) all show a general progression toward higher salinity during the BA, whereas the GOA δ18Osw-ivc records [EW0408-26JC (36), EW0408-85JC (33), and EW0408-87JC] show a brief increase in surface salinity at the start of the Bølling (~14.8 ka), followed by surface freshening of varying degrees among the sites through the Allerød. The relatively low SSTs in the EW0408-87JC record during the BA (relative to other sites) may produce anomalously low δ18Osw-ivc values for this interval, especially if low SST values in this core were related to a shifting seasonality of the alkenone producers to cooler months relative to the planktic foraminifera, which were most likely blooming during the spring and summer (39). In the average of the sites, the BA period has generally high δ18Osw-ivc (higher salinity), within the range of Holocene values.

Nearly all δ18Osw-ivc records and the average of all the sites show a strong and sustained interval of lower δ18O (surface freshening) associated with regional cooling during the YD, with some sites exhibiting the lowest values at this time (EW0408-85JC, JT96-09PC, and EW040-66JC). Records of the δ18O offset between the near surface–dwelling Globigerina bulloides (Gb) and thermocline-dwelling Neogloboquadrina pachyderma sinistral (Gb-Nps) indicate a greater offset during the early deglacial period and YD, consistent with greater near-surface stratification at these times (fig. S2). Diatom-based salinity reconstructions show no clear evidence for a freshening event during the YD (21, 28), but both of these records are of much lower resolution than the foraminiferal δ18O-based records presented here, with only one data point for the SO202-27-6 record (28) and two data points for the W8709-13PC record (21) within the YD interval, making it plausible either that these records are missing the event or that freshening during the YD was more restricted to the margins.

The δ18Osw-ivc records for the Holocene are generally not as complete as the SST records due to poor carbonate preservation in some sites, such as JT96-09PC and EW0408-87JC. This limits assessment of Holocene surface salinity changes in comparison with the reconstructed deglacial variability. In general, however, the available records indicate higher δ18Osw-ivc (elevated surface salinity) during the Holocene relative to the deglaciation (Fig. 3), as would be expected from modeled freshwater discharge to the region (18) in addition to the transitory impact of flood events.

Deepwater ventilation

The Northeast Pacific radiocarbon compilation is shown as a vertical slice of the water column, illustrating changes in benthic-planktic (B-P) age differences in both depth and time through the deglaciation (Fig. 4). Larger B-P ages reflect the presence of subsurface water masses that have been isolated from the atmosphere for longer durations (i.e., more poorly ventilated). B-P age differences were comparable during the BA and Holocene at intermediate depth sites [EW0408-85JC (40), JT96-09PC (41), ODP 1019 (42, 43), and EW0408-26JC; Fig. 4]. Data during the LGM are sparse, but B-P ages at sites ODP1019 and W8709A-13PC (42, 43) are similar to B-P ages during the BA and Holocene (fig. S3). The most pronounced increases in B-P age differences occur from 19.0 to 16.5 ka and 12.7 to 11.5 ka, roughly aligned with the intervals of peak freshening recorded in the Northeast Pacific paleosalinity records (Fig. 4). A transition to lower B-P ages occurred at ~16.5 ka and persisted through ~12.8 ka. During the YD, large B-P ages of up to 1500 years appeared in waters as shallow as 680-m depth (40). Low B-P ages (∼500 years) occur during both the Bølling and early Holocene transitions in the upper 1000 m [EW0408-85JC (40)]. The 14C differences relative to the contemporaneous atmosphere were also calculated (∆14C ‰) and result in similar patterns to those inferred based on B-P age differences, with the most pronounced increases in subsurface watermass radiocarbon ages occurring at the same time, but most prominently near 18-ka (fig. S4).

Fig. 4 Compilation of Northeastern Pacific surface hydrography and subsurface ventilation records from the LGM to early Holocene in the context of major deglacial climate fluctuations.

Records from top to bottom: NEEM ice core δ18O record (46) (A), the ODP 1019 SST record (35, 68) (B), the average Northeast Pacific SST record shown in Fig. 3 (C), the NGRIP ice core δ18O record (53) (D), the average Northeast Pacific δ18Osw-ivc record shown in Fig. 3 (E), and the x-ray fluorescence (XRF) Ca/Sr ratio in core U1308 in the North Atlantic inferred to reflect Heinrich events (45) (F). The shaded depth-time plot shows records of B-P radiocarbon age differences from marine sediment cores spanning 680- to 3680-m water depth in the Northeast Pacific (G). Data from shallow to deep: EW0408-85JC (680 m) (40), JT96-09PC (920 m) (41), ODP 1019 (980 m) (42, 43), EW0408-26JC/TC (1620 m), W8709-13PC (2,710 m) (42, 43), ODP 887 (3650 m) (65), and EW0408-87JC (3680 m). Blue bars denote periods of cooling, freshening, and increased B-P age in the North Pacific; the younger event is approximately coeval with the YD.

While our radiocarbon compilation shows a general progression toward reduced abyssal 14C ages at 16 ka, in conjunction with SST warming, we do not find any evidence for a pronounced Northeast Pacific deepwater formation event before the Bølling, as suggested by Rae et al. (15). In contrast, we see the most 14C-rich intermediate and deep waters during the BA and early Holocene, implying improved ventilation of North Pacific watermasses roughly in tandem with the reinvigoration of North Atlantic circulation (5), and at the same time as the appearance of hypoxia at intermediate depths along the Northeast Pacific margin (40, 41), which supports a scenario in which changes in export productivity and ocean warming are the primary drivers of regional hypoxia rather than reduced ventilation (33).

On the basis of the reconstructed εNd of intermediate and deep waters in the GOA, Du et al. (44) infer that the depleted radiocarbon water masses during the early deglaciation and YD in the North Pacific were sourced from end-member changes in the ventilation of Antarctic Bottom Water. In this scenario, we cannot directly attribute the increase in radiocarbon age of deep water to regional surface freshening. However, warming in the Southern Ocean is seen in response to freshwater forcing in the North Pacific in one model (2), so it remains a possibility that ventilation anomalies generated in the Southern Ocean may have been mechanistically linked to North Pacific freshening.

Global SST synthesis

To place our regional Northeast Pacific data compilation in a global framework, we compiled high-resolution SST reconstructions (n = 135) and generated maps of SST anomalies for major deglacial climate events to reveal likely sources of climate forcing (Fig. 5). Climate intervals shown include the early and late parts of HS1, the BA, YD, and the early Holocene (see Materials and Methods for time interval designations). The early HS1 (16.5 to 18.0 ka) anomaly shows modest cooling (~0.5° to 1°C) across the North Pacific, whereas cooling in the North Atlantic is most pronounced in a band between ~30° and 45°N, especially in the eastern basin, near the Iberian margin. Core sites northward of 45°N in the Atlantic show either a slight warming or no discernible SST change associated with the early HS1 interval. This may imply that the large SST response near the Iberian margin is partly a function of migrating fronts related to the position of the Gulf Stream current.

Fig. 5 Global SST anomalies for various deglacial climate intervals.

Climate intervals shown are early HS1 relative to the LGM (A and B), late HS1 relative to early HS1 (C and D), the BA relative to early HS1 (E and F), the YD relative to the BA (G and H), and the early Holocene relative to the YD (I and J) (see Materials and Methods for data references and dates of climate intervals). SST anomalies for each core site are plotted on the left panels, and an interpolated version using a weighted average grid scheme is shown in the right panels. Maps were generated using Ocean Data View (82).

The SST anomalies associated with the late HS1 interval (15.0 to 16.4 ka) are more spatially heterogeneous across the Northern Hemisphere than the early HS1 period, with evidence for modest warming along the Northeast Pacific, in contrast to an enhanced or sustained cooling in the Labrador Sea and midlatitude North Atlantic, which is most likely related to ice rafting from the Hudson Strait associated with Heinrich event 1 (45).

The BA (13.0 to 14.6 ka) and early Holocene (10.0 to 11.5 ka) periods show greater SST warming in both the Northeast Pacific and midlatitude North Atlantic. SST anomalies north of 60°N in the North Atlantic show either no change or reverse SST anomalies to the dominant trends in the midlatitude North Atlantic for each of the major deglacial climate transitions (HS1, BA, and YD), whereas the Holocene transition is the only one that shows a clear penetration of the warm North Atlantic current into the Nordic Seas (9). These observations call into question whether deglacial changes observed in the strength of the AMOC at lower latitudes (5) are related to changes in deepwater formation in the Nordic Seas, such as predicted in many models driven by freshwater forcing to the high-latitude North Atlantic (7, 46). Although the underlying mechanisms driving the SST change along the Northeast Pacific margin remain somewhat enigmatic, warming of this region would have contributed to Arctic warming (13) during both of these rapid Northern Hemisphere climate events.

Cooling during the YD is most pronounced along the Northeast Pacific and midlatitude North Atlantic. Both the Northeast Pacific (Fig. 4) and North Atlantic (5) show evidence for changes in ocean circulation during the YD, implying some degree of local forcing in each region. Freshwater routing has been suggested through the St. Lawrence River (47) and via the Arctic (18, 48) during the YD, which is inferred to have contributed to a decline in AMOC strength. Here, we also show evidence for freshening and cooling along the Northeast Pacific during the YD, indicating that multiple drainage pathways were activated at this time or that closely sequenced meltwater events helped to prolong the YD cooling event. While the timing of the submergence of the Bering Strait remains controversial, the presence of a Pacific bivalve dated to 13 ka in the Arctic suggests a possible YD-age breach of the Bering gateway (30, 32). The regions of greatest SST cooling during the YD in the global compilation largely correspond to the modeled regions of maximum surface freshening in our flood simulation with an open Bering Strait, raising the possibility that Pacific meltwater could have also been partially routed into the Arctic. However, estimates for the early phase of the Bering Strait submergence (13 to 11.5 ka) suggest that it may have only been a few meters below sea level (32), leaving it unclear whether significant throughflow could have occurred. Future work will be required to better constrain the timing of the Bering Strait opening and identify the specific sources of fresh water observed in our compilation (i.e., Columbia River floods or direct input from marine-terminating Cordilleran glaciers).


While many modeling studies have examined the global climate impacts of freshwater input into the North Atlantic, to date, only two models have examined the global effects of freshwater input to the North Pacific (2, 28). In one model, “hosing” (0.1 to 0.3 sverdrups over 100 to 200 years) in the Northeast Pacific elicits broad-scale Northern Hemisphere cooling and pronounced Northern Hemisphere sea-ice expansion, as well as a bipolar seesaw warming effect in the Southern Ocean (2). Modeling studies that impose a cooling in the North Pacific also yield a downwind cooling effect on the North Atlantic region and significant cooling and sea-ice expansion in the Arctic (12, 13). Some modeling studies that induce a freshwater anomaly in the North Atlantic result in a simulated warming in the Northeast Pacific (2, 14, 49), while others produce diffuse downwind cooling across the North Pacific (49, 50). However, the models that drive cooling in the North Pacific as a response to Atlantic cooling exhibit a stronger SST response in the Northwestern Pacific (where cold westerly winds exit Asia) (49, 50), which is at odds with our proxy reconstructions that show that SST anomalies are often amplified in the Northeastern Pacific during abrupt deglacial climate fluctuations (Fig. 5).

Synchronous cooling in the North Atlantic and Northeast Pacific from ~18.0 to 16.5 ka (i.e., early HS1) and again during the YD is thus consistent with a contribution to climate forcing from Northeast Pacific freshwater input from the Cordilleran region (Fig. 5). Reconstructions of sea ice across the North Pacific show the most extensive sea-ice cover during both HS1 and the YD (51), corroborating regional cooling during these periods. The SST anomaly patterns during early HS1 do not appear consistent with freshwater forcing from the Fennoscandian Ice Sheet (2), which is often suggested as a possible driver of circulation changes before Heinrich ice-rafting events (5). Both the early HS1 (18.0 to 16.5 ka) and YD intervals show an overall north-south seesaw in SST patterns, with warming in the Southern Ocean more pronounced in the Pacific basin than in the Atlantic sector (Fig. 5). Given that a bipolar oceanic seesaw may emerge in response to both North Atlantic (49) and North Pacific freshwater hosing (2), this SST pattern is not necessarily diagnostic of a single freshwater forcing region.

Cooling of the surface ocean in the North Atlantic preceded Heinrich events, implying that massive ice-rafting discharge events may be a response to Northern Hemisphere cooling, not the initial trigger of cooling (3, 45, 52). The δ18O records from various Greenland ice cores show contrasting histories during the early deglaciation (46, 53), which are not explained by time scale uncertainties. For example, the NEEM (North Greenland Eemian Ice Drilling) δ18O record (46) shows the coldest temperatures during the early part of HS1, followed by a gradual warming leading up to the abrupt BA transition, whereas the NGRIP (North Greenland Ice Core Project) δ18O record (53) shows the coldest period in the second half of HS1—the time period associated with Heinrich events and the most significant decline in AMOC strength (Fig. 4) (5). The timing and pattern of change in Northeast Pacific SST records are most consistent with δ18O trends in the NEEM ice core (77°N) (Fig. 4), which receives a greater fraction of Pacific moisture than the more southern Greenland ice cores, such as NGRIP (75°N) and GRIP (72°N) (46). Therefore, differences in North Pacific and North Atlantic hydrography may, in part, be reflected in differing δ18O histories in the various Greenland ice cores during HS1.

Evidence for changes in Northeast Pacific deepwater ventilation during the earliest part of the deglaciation (Fig. 4) occurs with surface cooling, major meltwater floods, and ice-rafting events, pointing toward a role for changes in the Northeast Pacific region as an instigator of early deglacial circulation changes. These events were coeval with both the onset of deglacial warming in Antarctica (54) and North Atlantic cooling that predates Heinrich ice-rafting events (3). Pulses of ice-rafted debris in the GOA between 17.5 to 16.5 ka (17, 27, 55) occurred before those in the North Atlantic associated with Heinrich event 1 (16.2 ± 1.6 ka and 15.1 ± 1.6 ka) (45). Amplifying responses in the North Atlantic region associated with AMOC reductions (5) may have helped to prolong Northern Hemisphere cooling that was initially triggered by Pacific freshening/cooling events. Similarly, we find that a gradual sea surface warming before the abrupt Bølling warming occurred in the Northeast Pacific at ~16.0 ka (17, 33, 38, 55, 56), which was 500 to 1000 years earlier than warming in the central North Atlantic (6). We hypothesize that North Pacific warming may have been driven, in part, by Heinrich event 1 in the North Atlantic at 16.2 ka (45), consistent with climate models that show that Northeast Pacific warming would have occurred in response to North Atlantic freshwater forcing (5, 14, 49), in addition to deglacial warming driven by atmospheric CO2 rise (6).


Rather than a single freshening event at ~16 ka (28), we find evidence for regional heterogeneities and multiple intervals of reduced surface salinity along the Northeast Pacific during the last deglaciation. The most sustained anomalies occurred during the early deglaciation (19.0 to 16.5 ka) and YD periods, with accompanying decreases in SST and increases in deepwater radiocarbon age (i.e., older deepwater mass indicating a reduction in ventilation). The co-occurrence of surface cooling with surface freshening and larger B-P age differences in the North Pacific points toward changes in North Pacific Ocean circulation as a contributor to regional climate changes, rather than a simple “downwind” atmospheric transfer of surface air temperature anomalies from the North Atlantic across Europe and Asia. Abrupt increases in SST and reduction of the apparent radiocarbon age of abyssal waters were paralleled between the North Pacific and North Atlantic during the transitions into the BA and Holocene, supporting a dual role for North Pacific and North Atlantic heat transport in abrupt Northern Hemisphere warming events (17).

While the timing of major meltwater events from the Cordilleran and Laurentide ice sheets appears to have been somewhat asynchronous during the early deglaciation, cooling during the YD was generally synchronous and of similar magnitude between sites in the Northeast Pacific and midlatitude North Atlantic, pointing to essentially synchronous meltwater influences in both ocean basins in this abrupt Northern Hemisphere cooling event. While dates for the submergence of the Bering Strait remain controversial, our results suggest that if the strait was breached at the beginning of the YD (30, 32), then some of the meltwater funneled into the Northeast Pacific may have been transported through the Arctic and into North Atlantic deepwater formation regions. Thus, major meltwater pulses to the Northeast Pacific and attendant impacts on ocean circulation and sea-ice formation may have played a hitherto underrepresented role in abrupt deglacial climate variability.


Modeled freshwater routing

All numerical model simulations were performed using the Massachusetts Institute of Technology General Circulation Model (MITgcm) (57). The model has an eddy-permitting 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 were 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 was coupled to a dynamic-thermodynamic sea-ice model that assumes a viscous plastic ice rheology and computes ice thickness, ice concentration, and snow cover (58, 59). The simulations performed to study the propagation of meltwater from the Columbia River in the geological past were integrated under both glacial and deglacial boundary conditions. The glacial configuration follows that of Hill and Condron (60): Sea level is 120 m lower than modern day, and the atmospheric boundary conditions (10-m wind, 2-m air temperature, surface humidity, downward longwave and shortwave radiation, precipitation, and runoff) were provided from output from the fully coupled Community Climate System Model version 3 (CCSM3) LGM integration (61); as the Bering Strait is above sea level, there is no oceanic exchange between the Arctic and the North Pacific. The deglacial configuration follows that of Condron and Winsor (62) with flow through the Canadian Archipelago blocked by the Laurentide Ice Sheet joining the adjacent Greenland Ice Sheet. For simplicity, sea level was set to modern day, and atmospheric boundary conditions were prescribed from European Centre for Medium-Range Weather Forecasts ERA-40 reanalysis data. Both model configurations were spun up for 50 years before control and perturbation experiments were performed. In each glacial and deglacial perturbation simulation, 3 × 106 m3 s−1 of meltwater was released for 1 year with a temperature and salinity of 0°C and 0 psu, respectively, into the four model grid cells closest to the mouth of the Columbia River.

The size of individual Missoula Flood events has been estimated between 0.5 and 17 sverdrups (1924, 63), with up to 89 events within a period of a few thousand years (25). Alho et al. (63) estimate peak flows of 17 sverdrups for a 70-hour duration for the largest Missoula Flood events. These floods—among others sourced from glacial lakes Bonneville and Columbia (26)—would have been superimposed on elevated deglacial freshwater discharge ranging from 0.02 to 0.06 sverdrups, with the higher estimates (0.06 sverdrups) expected during the early deglaciation (18), during times of peak volume and reoccurrence of Missoula Flood events (Fig. 3). A large flood event with a volumetric equivalent of 105 km3 has also been hypothesized to have included subglacial reservoirs from British Columbia in addition to glacial Lake Missoula based on geomorphic and sedimentary features of the Channeled Scablands (64), although this estimate is controversial (21). Modeling multiple flood events over a multicentennial time scale was not computationally feasible with existing resources; therefore, a single large flood event scaled to the estimates of Shaw et al. (64) was used to trace the freshwater trajectory and gauge the decadal persistence of a single flood event. While this event exceeds the volume of water estimated for a single flood sourced from glacial Lake Missoula, the modeled salinity anomaly is on par with paleosalinity reconstructions that indicate a 4-psu anomaly in marine cores at 42°N early in the deglaciation (~18 ka), inferred to reflect the largest of the Missoula Floods (21). If multiple flood events occurred within the time frame of freshwater dissipation, their impact on Pacific hydrography would likely reinforce each other. Future work will explore the persistence and aggregation effects of multiple smaller flood events but is beyond the scope of this paper. Here, we present our modeling results to illustrate the trajectory of fresh water injected into the Northeast Pacific, rather than as a quantitative assessment of regional salinity variations in response to a single Missoula Flood event.

Marine sediment cores

New and previously published paleoceanographic proxy data were synthesized from marine sediment cores along the Northeast Pacific margin, including core sites EW0408-85JC (59.56°N, 144.15°W, 682 m), EW0408-87JC (58.77°N, 144.50°W, 3,680 m), EW0408-66JC (58.45°N, 137.17°W, 426 m), and EW0408-26JC (57.60°N, 136.72°W, 1,623 m) in the GOA; JT96-09 (48.90°N, 126.88°W, 920 m) off Vancouver Island; and ODP site 1019 (41.68°N, 124.93°W, 980 m). These cores were chosen because each site had an alkenone-based SST reconstruction, a planktic δ18O record, and paired benthic and planktic radiocarbon data for the last deglacial period. A Mg/Ca-based SST record and planktic δ18O records from a marine core near Vancouver Island, MD02-2496 (48.98°N, 127.04°W, 1243 m) (38), were also used for the SST and δ18Osw-ivc compilations to provide better data coverage for the early part of the deglaciation. B-P radiocarbon data from core sites W8709A-13PC (42.12°N, 125.75°W, 2,712 m) (42, 43) along the California margin and ODP 887 (54.37°N, 148.45°W, 3,647 m) (65) in the GOA were also included in the radiocarbon depth transect (in addition to the cores listed above) to provide additional ventilation estimates for the LGM and greater depth coverage for intermediate and deepwater ventilation changes.

Age models

Radiocarbon ages for the individual core sites in the alkenone compilation were first calibrated with the Marine13 calibration curve (66) using the Calib7.1 software (67) with a constant marine reservoir correction. The age model for EW0408-85JC is based on 36 mixed planktic foraminiferal radiocarbon dates assuming a marine reservoir correction of 880 ± 80 years (ΔR = 480 ± 80 years) (40). The age models for cores EW0408-26JC/TC consist of 10 radiocarbon dates on mixed planktic foraminifera using a marine reservoir correction of 735 ± 50 years (ΔR = 335 ± 50 years) (17). The age model for EW0408-66JC consists of 19 radiocarbon dates on mixed planktic foraminifera using a marine reservoir correction of 595 ± 50 years (ΔR = 195 ± 50 years) (17). The age model for core EW0408-87JC is based on 18 mixed planktic radiocarbon dates using a marine reservoir correction of 850 ± 100 years (ΔR = 450 ± 100 years) (33). The age model for JT96-09PC (34) was updated with the Marine13 calibration curve (66) using the Calib 7.1 program (67) and using 10 radiocarbon dates on mixed planktic foraminifera with a marine reservoir of 900 ± 100 years (ΔR = 500 ± 100 years). The age model for ODP site 1019A is based on 11 radiocarbon dates, 10 of which are on mixed planktic foraminifera and 1 of which is a piece of bark (42), recalibrated here with the Marine13 calibration curve using a marine reservoir age of 900 ± 50 years (ΔR = 500 ± 50 years). A modified mean core depth splice was created between ODP 1019A and 1019C based on the overlapping alkenone SST records (table S1 and fig. S5) to refine the age model transfer from core 1019A to 1019C in the composite alkenone SST record, which splices data from core 1019C (35) with data from core 1019A (68). The age model for ODP site 887 was constructed from seven calibrated planktic radiocarbon dates (65) with a total marine reservoir correction of 900 ± 100 years (ΔR = 500 ± 100 years). The age model for core W8709-13PC is based on 22 planktic radiocarbon dates (42) with a marine reservoir correction of 720 ± 100 years (ΔR = 320 ± 100 years) (43). The age model for core MD02-2496 was adopted from Taylor et al. (38).

Most of the Northeast Pacific alkenone SST records show similar patterns and magnitude of abrupt warming and cooling transitions during the deglaciation. However, there are apparent differences in the timing of these abrupt transitions between core sites when applying constant marine reservoir corrections in the construction of age models (fig. S6). Therefore, it is likely that apparent age differences between sites (exceeding 200 years) reflect changes in local marine reservoir age rather than true leads and lags of several hundreds of years between sites that are within the same geographic region. We therefore devised alternative age models for these cores by allowing relatively minor changes to the original marine reservoir correction (±100 to 300) to produce age models that better aligned with the abrupt δ18O shifts in the absolutely dated (U/Th) speleothem records of the western United States during the BA, YD, and early Holocene transitions. This includes the Oregon Caves National Monument δ18O record that spans 13 to 0 ka (69, 70) and the Cave of Bells speleothem record that spans the entire deglaciation (71). This is justified by observations that hydrologic changes in the western United States were strongly correlated to SST variability in the North Pacific (70, 71). While these two speleothem records are not located in close proximity to one another, the similar timing for the transition into the YD and Holocene in overlapping time intervals suggests generally cohesive hydrological changes affecting the entire western United States in conjunction with abrupt climate transitions.

The timing of the δ18O transitions in the western U.S. speleothem records is also consistent (within error) to the timing of abrupt δ18O transitions observed in the Greenland ice cores for the Inter-Allerød cooling, YD, and Holocene transitions. However, the midpoint of the “Bølling” transition in the Cave of Bells speleothem occurred earlier (15 ± 0.2 ka) (71) than the abrupt Bølling transition in Greenland (14.7 ± 0.2 ka) (53). It is unclear whether this age offset reflects a true lead in the abrupt hydrological changes that occurred in the mid-deglaciation in the southwestern United States or there are unresolved chronological errors. We used a “conservative” approach and adopt marine reservoir corrections that better align with the later timing of the Bølling transition in the Greenland chronology; however, this apparent discrepancy in timing of the BA transition between these records underscores the need for future studies to resolve uncertainties in the relative timing of abrupt climate fluctuations between the Northeast Pacific and Greenland/North Atlantic regions.

Age models were constructed using the midpoint age of the 2-sigma calibrated radiocarbon age range, with total marine reservoir ages ranging from 500 to 1100 years (∆R = 100 to 700 years in Marine13). Inferred changes in the surface marine reservoir age were similar among the Northeast Pacific margin sites from 48° to 60°N; marine reservoir ages were 100 to 300 years younger during the YD and early Holocene (between 9.5 and 11.5 ka in 14C age) relative to the Bølling transition. In contrast, site 1019 shows evidence for anomalously old planktic radiocarbon ages near the Bølling and Holocene transitions, which may indicate brief events in which older deep waters are upwelled to the surface, possibly associated with the reinvigoration of North Pacific Ocean circulation. An additional tie point for the Holocene transition (11.58 ka) was incorporated into the age model for core JT96-09PC due to sparse radiocarbon dates through this interval. The constant marine reservoir age models for core ODP 887 and W8709-13PC were retained because there were no alkenone-based Uk′37 SST records available to refine estimates for possible marine reservoir variations through time. We note that in no cases did the revised “variable marine reservoir” age model exceed ±380 years from the original applied (i.e., constant) marine reservoir correction. All radiocarbon dates, marine reservoir corrections, and relevant age model data are provided in a supplementary Excel file.

Oxygen isotope records

New planktic oxygen isotope records were generated for core EW0408-87JC (fig. S2). Stable isotope analyses were analyzed at the Oregon State University Stable Isotope Laboratory using a Kiel III carbonate preparation device connected to a Thermo Finnigan MAT 252 mass spectrometer. Planktic foraminifera Nps and Gb species were picked from the 150- to 255-μm size fraction. Approximately 25 specimens of Nps and 15 specimens of Gb were used for each sample measurement. New planktic Gb and Nps δ18O records from core JT96-09PC were measured at Geotop, University of Quebec at Montreal, by dual inlet mass spectrometry using a MultiCarb acidification system connected to an Isoprime isotope ratio mass spectrometer, as well as at the Oregon State University. Data were reported in δ notation, δ18O = [(18O/16O)sample/(18O/16O)standard] − 1), relative to the Vienna Pee Dee Belemnite standard, with an average external precision of ±0.04‰ based on replicate analyses of standards.

A previously published planktic Nps δ18O record was used for calculating the δ18Osw-ivc in ODP 1019 (42). Published records of Nps and Gb δ18O from cores EW0408-85JC (55), EW0408-26JC (17), and EW0408-66JC (17) were used for the δ18Osw-ivc records in these cores (33, 36). The δ18Osw was calculated using the 11-chambered Gb temperature equation of Bemis et al. (72) [δ18Osw = 0.27 + δ18Oc + (T − 13.2)/4.89)] and corrected for global ice volume using the δ18O-equivalent sea-level curve (denoted δ18Osw-ivc) (73). The global ice correction is relatively low resolution and may underestimate some transient effects (73). There is also the potential for some nonstationarity in the δ18Osw-salinity relationship through time (74), but the magnitude of these uncertainties appears to be in the range of the estimated uncertainties for the δ18Osw-ivc calculation.

Because many datasets did not have precise sample pairs for the SST and δ18O measurements, records were linearly interpolated on 100-year time steps, and an average record of all the Northeast Pacific margin sites (excluding core EW0408-87JC for reasons discussed below) was produced by averaging overlapping time intervals of the 100-year interpolated records, with the requirement that the average contained at least two overlapping records for any given time interval. Uncertainties for the reconstructed δ18Osw-ivc values were calculated as the root mean square propagation of the δ18O equivalent of the SST uncertainty (average = 0.15 ‰), the reconstructed δ18O of sea level (0.10 ‰), and the analytical uncertainty of the foraminiferal δ18O measurements (0.04 ‰) for each sample. For the interpolated records, the average uncertainty for each parameter (defined above) was used in constructing the total estimated δ18Osw-ivc uncertainty of ±0.19‰.

Alkenone SST estimates

Measurements of Uk′37 from EW0408 cores in the GOA were conducted at the Oregon State University in Fredrick Prahl’s laboratory. Data from cores EW0408-66JC and EW0408-26JC/TC were published by Praetorius et al. (36). Data from EW0408-85JC were published by Praetorius et al. (33). New data from core EW0408-87JC are presented here, following the same protocols as outlined by Praetorius et al. (33, 36), including standard solvent extraction methods (75) and urea adduction techniques (76). Uncertainties in the quantification of di-unsaturated (K37:2) and tri-unsaturated (K37:3) alkenones were calculated as the sum of an average baseline variability component and a 5% uncertainty in the integrated area and then propagated into calculation of the Uk′37 temperature index (77).

Previously published Uk′37 records from ODP 1019C (16 to 0 ka) are published by Barron et al. (35) and are combined and extended with data from 1019A and 1019C (68); average replicate error among samples was ±0.1°C. The alkenone Uk′37 record from core JT96-09PC was published by Kienast and McKay (34); average replicate error among samples was ±0.3°C (34). All Uk′37 data were calibrated to SST (°C) based on the equation of Prahl et al. (77). A total Uk′37 SST uncertainty of ±0.6°C was applied to the JT96-09PC and ODP 1019 SST records, in line with the average uncertainty (calibration and analytical) of the GOA Uk′37 SST records.

While most of the Uk′37 SST records from the Northeast Pacific margin show overall similar deglacial and Holocene trends, there are a few notable differences among records. For example, the magnitude of warming at the Bølling onset, and for the subsequent duration of the BA interval, was reduced in core EW0408-87JC relative to all other sites. The concentrations of alkenones were generally >0.5 μg/g in the BA interval in this core, and there were no other obvious issues with the chromatography in these samples that may lead to a spurious cold bias in the Uk′37 measurements. As this is the only core that is situated in an open-ocean upwelling gyre environment, contrasting with net downwelling along the highly stratified GOA margin, we hypothesized that the lower SSTs at this site relative to the margin sites during the BA reflect greater mixing with subsurface waters and/or a shift of the alkenone producers dominant bloom season to cooler months when wind-driven mixing enhances nutrient delivery. This inference was supported by sediment trap data that show a greater fraction of coccolithophore production during winter months in open-water environments in the North Pacific (78) relative to shallower marginal sites that are dominated by summer blooms (79). Additional geochemical or faunal-based SST reconstructions from this and other nearby sites will hopefully provide more insight into these regional differences.

Most of the Northeast Pacific SST records show maximum Holocene warmth occurring from 11.5 to 9.5 ka, followed by cooling into the mid-Holocene. The exception to this trend is core EW0408-66JC, which does not show a clear cooling trend after 9.5 ka; however, this core has low resolution and poor age controls through the Holocene, along with a missing core top that truncates the record ~4 ka. An inflection to gradual warming throughout the late Holocene is seen in ODP 1019 and EW0408-85JC at ~6 ka and is mirrored by the decrease in δ18O in the Oregon Caves speleothem record (70). Cores EW0408-87JC and JT96-09 show a more flatline SST trend in the late Holocene.

Given the limited alkenone-based SST data coverage for the LGM and earliest part of the deglaciation (20 to 18 ka), we also included a high-resolution Mg/Ca-based SST reconstruction from core MD02-2496 (38), offshore of Vancouver Island, to extend our Northeast Pacific SST and δ18Osw-ivc compilations. This record spanned ~21 to 13 ka and was derived from the Mg/Ca ratio of the near-surface dwelling planktic species Gb using the calibration of Elderfield and Ganssen (80). The average Northeast Pacific Uk′37 SST record (Fig. 3) was produced by linear interpolation of each record on a 100-year time step, followed by an average of the records (with the exception of EW0408-87JC) at each time step, and then a five-point (400-year) running average was applied to the record. Errors shown are the SEM. Core EW0408-87JC was excluded from the compilation due to the fact that this core had lower resolution (~300 years) than the other records (~100 years) and also because the BA interval in this record appears to have anomalously low SST values relative to the other sites; it is unclear at this time whether this reflects true SST gradients between the gyre and the margin sites or it is a function of shifting seasonality in the alkenone producers at site EW0408-87JC. Both the Northeast Pacific SST and δ18Osw-ivc average records are limited to data from sites ODP 1019 and MD02-2496 in the 18- to 20-ka age range because of the fact that the other core sites do not extend beyond 18 ka.

Deepwater radiocarbon

Paired planktic and benthic foraminiferal radiocarbon measurements were used for estimates of changes in ventilation and stratification. New records were produced from cores EW0408-87JC (3680 m) and EW0408-26JC (1623 m) in the GOA and compiled with previously published records from different depth core sites, including EW0408-85JC (682 m) (40) and ODP 887 (3647 m) (65) in the GOA, JT96-09PC off the Vancouver margin (920 m) (41), and ODP 1019 (980 m) (42, 43) and W8709A-13PC (2712 m) (42, 43) along the California margin. Ventilation ages were also calculated for planktic and benthic radiocarbon data based on the tuned SST chronology [∆14C (‰) relative to the contemporaneous atmosphere], which showed similar trends in deepwater ventilation throughout the deglaciation, with the most pronounced decreases in ventilation occurring during 23 to 22, 19 to 17, and 12.7 to 12 ka (fig. S4).

Surface salinity estimates

Given that the δ18Osw calculation is based on SST and δ18O records produced by different biota, there is potential for seasonal offsets between the alkenone-bearing coccolithophores and the planktic foraminifera to affect the inferred δ18Osw surface salinity reconstruction. Records of δ18Osw-ivc based on paired foraminiferal Mg/Ca-derived SST and δ18O from the same species (Gb) in a core from the Vancouver margin show similar SST and δ18Osw-ivc trends for the early deglacial period (here recalculated following the same procedure as the other δ18Osw-ivc records to maintain consistency; Fig. 3) (38). This provides greater confidence that the δ18Osw-ivc records based on alkenone-derived SST and planktic foraminiferal δ18O are monitoring regional trends in surface freshening and are not overly affected by seasonality differences in the biota.

We also calculated the difference in the δ18O between the near surface–dwelling symbiont Gb and the deeper thermocline-dwelling Nps species for the GOA core sites to provide an alternative method for evaluating differences in upper water column stratification and salinity (fig. S2). Differences in the dominant bloom season between Gb and Nps may still influence the Gb-Nps δ18O record; however, convergence of trends between the Gb-Nps δ18O record and the δ18Osw-ivc reconstruction provides greater confidence in broad-scale surface salinity trends throughout the deglaciation and Holocene. The Gb-Nps δ18O record shows a tendency for Gb to get lighter during the early deglaciation and the YD, whereas Gb gets heavier relative to Nps during warm climate periods, such as the Bølling and Holocene. These trends have previously been inferred to reflect surface freshening when Gb is depleted relative to Nps, especially when considering that the light Gb δ18O excursions in the early deglaciation are accompanied by pulses of ice-rafted debris (17). The averaged GOA Gb-Nps δ18O record reveals enhanced surface stratification during 18 to 16.0 ka and 13 to 12 ka, largely following trends in most of the Northeast Pacific δ18Osw-ivc records (fig. S2).

Changes in precipitation may also have contributed to regional variability in the δ18Osw. For example, increases in precipitation occur in the Northeast Pacific in response to modeled freshwater forcing in both the North Atlantic and North Pacific (28). However, the net change in surface δ18Osw is predicted to be positive in the North Pacific in response to North Atlantic hosing, whereas it is predicted to be negative in the North Pacific in response to North Pacific hosing and/or a combination of North Atlantic and North Pacific hosing (28). This implies that episodes of lowered δ18Osw in the Northeast Pacific are most consistent with regional freshwater forcing. Furthermore, the δ18O recorded in the Mt. Logan ice core in Alaska shows relatively depleted values during the HS1 and YD intervals, which are interpreted to reflect more zonal atmospheric flow and lower regional precipitation (81). These trends in precipitation would therefore appear to be inverse of the major trends in surface ocean δ18Osw observed in our compilation.

While we have attempted to assess and integrate multiple proxy reconstructions of surface salinity variations throughout the Northeast Pacific, greater spatial and temporal coverage of paleosalinity reconstructions will be necessary to refine our understanding of earlier (glacial) freshening events relative to the deglacial variations we present here, as well as to distinguish among specific sources of meltwater and the geographic extent of their impact.

Global SST compilation

A compilation of SST proxy reconstructions from marine sediment cores around the globe (n = 135) was assembled to examine the spatial patterns of SST anomalies during major deglacial climate events. These records include SST estimates based on the alkenone Uk′37 unsaturation index, planktic foraminiferal Mg/Ca, foraminiferal assemblage transfer functions, and TEX86. This database included most of the records incorporated into the SST compilation (n = 65) of Shakun et al. (6) but significantly expanded on the regional coverage of proxy reconstructions from the North Pacific Ocean (n = 27) and incorporated new records from different regions throughout the globe. Records were excluded if the average resolution for the deglacial interval was <500 years.

The SST anomalies were calculated on the basis of the following time intervals: LGM SST estimates are based on the interval 23 to 19 ka, the early HS1 SST anomaly is 18.0 to 16.5 ka relative to the LGM, the late HS1 SST anomaly is 16.4 to 15.0 ka relative to the early HS1, the BA SST anomaly is 14.6 to 13.0 ka relative to the early HS1, the YD SST anomaly is 12.7 to 12.0 ka relative to the BA, and the early Holocene SST anomaly is 11.5 to 10.0 ka relative to the YD. The climate intervals defined for the SST anomalies were typically shorter than the actual climate intervals (based on the Greenland ice core chronology) (53) to avoid climate transition boundaries given typical multicentennial age model uncertainties for marine cores. Metadata and calculated SST anomalies for each core site are included in a supplementary Excel file.

Data were plotted using the Ocean Data View software (82) with both the individual SST estimates at each site, as well as with a weighted average gridding scheme with a 90‰ cutoff in latitude and longitude (Fig. 5).


Supplementary material for this article is available at

Fig. S1. Modeled SSS, sea-ice, and SST anomalies in the deglacial simulation with an open Bering Strait for model years 5, 10, and 15.

Fig. S2. Planktic δ18O records from the Northeast Pacific.

Fig. S3. Compilation of Northeastern Pacific B-P radiocarbon age differences from the LGM through the early Holocene.

Fig. S4. Compilation of Northeastern Pacific subsurface ventilation records from the LGM through the early Holocene [with ∆14C (‰)] instead of B-P age differences.

Fig. S5. SST measurements for ODP cores 1019C and 1019A on the original mean core depth splice (red, 1019C; green, 1019A) show misalignment near the Holocene transition by ~30 cm and the BA and YD transitions by ~5 cm.

Fig. S6. Compilation of Northeast Pacific Uk′37 SST reconstructions on age models with constant marine reservoir correction, in comparison with annual layer counted chronologies of the Greenland ice cores (46, 53) and the absolute U-Th chronologies of the western U.S speleothem records (6971).

Table S1. A revised mean composite depth for the upper 8 m of ODP core 1019A based on a revised correlation of GRA and MS data to the other cores within the splice for the upper 7 m.

Table S2. Proxy data used for the global SST anomaly maps in Fig. 5 [(38, 45, 83165) listed in alphabetic order].

Data file S1. Northeast Pacific SST data.

Data file S2. Northeast Pacific oxygen isotope data.

Data file S3. Northeast Pacific radiocarbon data.

Data file S4. Metadata for global SST compilation.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank K. Brewster, J. Padman, F. Prahl, A. Ross, and M. Wolhowe for assistance with laboratory procedures; J. Rodysill, C. Davis, J. Addison, J. Barron, A. Schmittner, J. Alder, and the anonymous reviewers for suggestions that improved the manuscript; and all the researchers who made their data available, either upon request or through public databases. Funding: The research was partly supported by the NSF through grants ARC-257 1204045 and PLR-1417667. The numerical model simulations used resources from the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under contract no. DE-AC02-05CH11231. Author contributions: S.K.P. designed the study with input from A.C.M. A.C. designed and ran the freshwater routing simulations. S.K.P. produced new alkenone Uk′37 SST data from core EW0408-87JC and new radiocarbon data from cores EW0408-87JC and EW0408-26JC and compiled all data. J.L.M. produced new δ18O measurements from core JT96-96. S.K.P. wrote the paper with input from all coauthors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article