Introduction

The LGM is traditionally defined as the interval in 23−18 thousand years ago (ka) between the HS2 and HS1 chronozones, when the maximum amount of ice was sequestered in global continental ice sheets, corresponding to the lowest stand of the global sea level1,2. The current database indicates that the world’s largest LIS, reached a maximum at approximately 22 ka3, while the lowest global sea level occurred at approximately 21 ka4. However, recent works show that the LIS advanced to its southernmost margin during the late phase of the HS2, defined as the peak glaciation, to its second southernmost margins during the early phase of traditionally defined LGM, and to its third southernmost margin during the early phase of HS15. Here, the LGM definition is slightly modified as the climate chronozone spanning from 22.5−17.8 ka between the warm Greenland Interstadial 2 (GI2) phase from 23.5−22.5 ka, the tipping interval when the LIS begins to shrink, and the cold early HS1 phase from 17.8−16.3 ka. The traditional LGM contained a very cold early and a warm late phase, while the warm phase caused a large-scale meltwater flood5. To ascertain the exact timing of these cooling and warming phases during the LGM and the associated LIS advance and retreat behaviors, we compiled and synthesized the most representative and accurate AMS 14C dates from our database to determine the exact timing and duration of these early LGM/HS1 advances (Fig. 1a). The timing of 2 separate ice sheet advances constrains the timing and duration of LIS retreat with a large-scale meltwater flood in the late LGM. Our compilation further reveals that the Eurasian Ice Sheet (EIS), Tibetan Mountain glaciers, world eolian and lacustrine sediments, and global stalagmite δ18O records all show cold early and warm late LGM phases followed by a less-cold early HS1 phase. Global ice sheet became thinning at 21 ka6, corroborated by the substantial decay of the LIS in the late LGM from 20.5-17.8 ka (this study), raises the question of how much LIS meltwater drained into the North Atlantic Ocean to influence NH climate changes. For this purpose, we applied the geophysically constrained ICE-6G model to reconstruct the LIS and measure the vertical motion of the crust constraining the thickness of the local ice cover and find that the release of meltwater in 19 ka contributed only 1.6 m to global sea level rise, thus having a limited impact on the North Atlantic Meridional Circulation (AMOC). However, the decay of the southern LIS margins could have altered the LIS morphology, thereby impacting global westerly propagation. We also analyzed summer soil temperature variations near the southern margins of the LIS and found mean summer temperatures of 10.4 and 11.8 °C in the early and late LGM, respectively, and 11.4 °C in the early HS1 phase. Climate stratigraphic records on the LIS retreating route have shown poleward-equatorward migrations of westerlies near the ice margins5. Here, we adopted a numerical model to simulate this shift in North American westerlies over North America and the North Atlantic in response to an LIS reduction from its full LGM boundaries. We restate that the LIS grew to its largest size during the late HS2 phase from 24.5−23.5 ka5, representing the real LGM, when the LIS southern margin was approximately 80 km south from its position in the so-called LGM defined in this study (Fig. 1b). By the late LGM, the southernmost margins of the LIS had retreated at least 400 km north of Milwaukee, Wisconsin, U.S.A. to form Glacial Lake Milwaukee, which is much larger than the Glacial Lake Chicago. Here, we propose that the poleward−equatorward displacement of North American westerlies coupled with the advance and retreat of the LIS amplified the small changes in NH seasonal insolation responsible for the cold-warm-cool climate shifts. In our intra D-O/Heinrich scale climate stratigraphy analysis, the first LIS retreat was found to have occurred in the GI2 phase, designated as the default deglaciation because the LIS never readvanced to this position again after GI2 (Fig. 1b). We reference the framework of the intra D-O/Heinrich-scale climate changes in Supplementary Fig. 15.

Fig. 1: Maps of the ice-margin positions of North America.
figure 1

The Lake Michigan lobe (red dashed square), and study sites in the Illinois River Valley (IRV) and Mississippi River Valley (MRV) during the last deglacial period. a Sites with data compiled in the paper: 1, Hennepin; 2, Henry; 3, Lacon; 4, Chillicothe; 5, Oswego Channel; 6, Bryn Mawr Lake; 7, Rice Lake; 8, Manito; 9, Kearney Cemetary Lake; 10, Jules; 11, Canteen Creek; and 12, Keller Farm. b AMS 14C dates of tundra plant leaves, tree stump rings and α-cellulose (Supplementary Table 1) indicate that the Lake Michigan lobe of the southern LIS margin advanced to the outermost margins during the HS2a phase, forming the Shelby-Bloomington Moraine systems (green) and to the 2nd- and 3rd- outermost margins during the LGMb stage, forming the Eureka-Arlington-Minonk-Marseilles Moraine system (red), and during the HS1b phase, forming the Valparaiso Moraine system (yellow), respectively. c Location of research area at the southern LIS sector. See Supplementary Fig. 1 for the definition of intra D-O/Heinrich-scale climate stratigraphies5. Note: LGMb-phased ice margin advances (3 red lines) are not differentiable by AMS 14C dates. The LGMb very-cold phase contains 2 minor ice-decay lukewarm subphases.

Full size image

Results

Accurate timing of LIS advance and retreat

Short-lived fossil materials recovered and identified from proglacial and ice-walled lakes in front of the ice sheet, e.g., kettle lakes on glacial sediments and ice-walled lakes on moraines, yield the best and most representative 14C dates in our large 14C age database7. To ascertain the accurate timing and duration of LIS margin advances in the central U.S.A., AMS 14C dates derived from short-lived terrestrial fossil plants were used5,7,8. To synthesize the timing of the LIS advancing to its outermost margins, we selected AMS 14C dates from fossil plants, caddisfly nests, the outer 2 tree rings and tree-ring α-cellulose of a tree stump, from proglacial lakes derived from recently evaluated 893 14C dates7; all these proxies are associated with the timing of the development of southernmost Bloomington/Shelby Moraines at 39.4 °N; 88.7 °W. Here, we report that the tree stump was positioned vertically in the proglacial lake, truncated by the advancing Bloomington/Shelby ice margin and preserved well in sediment. These accurate AMS 14C dates reveal that the LIS reached its maximum size during the late HS2 intra-phase at ~24.5−23.5 ka; this is an important finding for global paleoclimate studies (Figs. 1b, 2, Supplementary Fig. 1; Supplementary Table 1)5.

Fig. 2: Probability analysis of AMS14C dates and implication for LGM climate phase changes.
figure 2

a Cumulative probability curve of 19 AMS14C dates with the 2σ values derived from plant leaves in ice-walled lakes on Marseilles and Valparaiso Moraines5,7,13,14 shows that the LIS advanced during the early LGM and early HS1. The absence of ice-walled lacustrine sediments without 14C dates indicates ice-margin decay during the late LGM and GI2. During the HS2a, the LIS expanded to its maximum size, known as the local LGM. The GI2 ice retreat represents deglaciation by default because LIS has never reached the same size since that time. b The probabilities of 16 AMS and radiometric 14C dates obtained from fossil plants, wood fragments, and soil organics with 2σ values associated with meltwater deposits in the Illinois and Mississippi River valleys (Fig. 1a) show that the LGMa-phased meltwater flood occurred from 19.4 to 18.3 ka (red dashed line). c LIS thickness in meters from the ICE-6G reconstruction15 at 23 ka and (d) at 16 ka. The black box indicates the mask used in our ice volume changes estimation in 23-16-ka slices along the southern LIS margins, excluding any ice volume changes in other regions. e Soil brGDGT temperature reconstruction indicates a cold LGMb, warm LGMa, and less-cold HS1b phase at the Keller Farm loess-paleosol exposure, in which 2 incipient paleosols in the LGMb unit agree with 2 minor ice retreats during the LGMb (Fig. 1b). f Ice margin positions in the Great Lakes and Northeast regions from 21.2 to 18.8 ka3,16. An LIS decay of approximately 900,000 km2 (orange). Meltwater in Glacial Lake (L) Milwaukee drains into the Illinois River via the Chicago outlet, from the Lake Wisconsin into the Fox River, and from the Lakes Saginaw and Leverett partially into the Milwaukee and Kankakee Rivers to the Gulf of Mexico. Meltwater in Lakes Saginaw and Leverett partially drained to the North Atlantic Ocean via the Hudson River.

Full size image

To synthesize the timing of the LIS advancing to its 2nd outermost margin (which it did several times), we selected AMS 14C dates derived from terrestrial fossil plants growing in ice-walled lakes on glacial moraines and/or on the last glacial tillplain from our database7. The results reveal that the Eureka-Arlington-Minonk-Marseilles morainal systems advanced to their 2nd-outermost margins (several times) at approximately 40.1 °N; 88.2 °W during the early LGM from 22.5 to 20.5 ka. It is noticeable that these moraines display 3 northeast‒southwest and half circles ~30 km apart, importantly suggesting that during the LIS advances within the cold early LGM phase, at least 2 short lukewarm intervals occurred. Although our AMS 14C dates cannot significantly differentiate the timing of the 3 ice margins, the large spatial patterns should reflect climate changes at a higher resolution, e.g., on the intra D-O/Heinrich scales (Figs. 1b, 2, Supplementary Fig. 1; Supplementary Table 1)5.

To synthesize the timing of the LIS advancing to its 3rd-outermost margins, we also select AMS 14C dates derived from terrestrial fossil plants growing in ice-walled lakes on the Tinley-Deerfield-Valparaiso moraines; this proxy reveals that the Valparaiso ice margin advanced to its 3rd-outermost margins at 41.5 °N; 87.1 °W during the early HS1 phase from ~17.8−16.8 ka (Figs. 1b, 2, Supplementary Fig. 1; Supplementary Table 1)5.

However, it is challenging to provide an accurate and reliable 14C chronology to constrain ice margin retreat-induced meltwater events due to the lack of datable materials. Here, we provide a new approach to solve this historic problem. When the southern margins of the LIS advanced and overrode the Illinois River valley during the late HS2 and early LGM phases, two till stratigraphic units formed across the valley, namely one from HS2 and the other from early LGM. However, while both till units are found in the east and west valley walls, they are absent from the river valley itself9 (Fig. 1b). It thus becomes important observation to determine the timing of the LIS decay and associated meltwater release. We suggest that the meltwater floods that eroded the 2 till units in the Illinois River valley must have occurred after the 2nd till unit was deposited, which was known as Eureka-Arlington-Minonk-Marseille morainal system in the early LGM-phase9. Because the early HS1-phased Tinley-Deerfield-Valparaiso ice margins readvanced to south of the Glacial Lake Chicago and shut down the Chicago Outlet, stopping meltwater from being released into the Illinois River, we predict that meltwater-induced outwash deposition during LIS decay should have ended by the onset of the HS1 phase (Fig. 1b). Thus, the chronology of LIS decay-induced meltwater floods should be constrained in the late LGM warming phase from 20.5-17.8 ka. The concept of this warm late-LGM phase was also corroborated by a formation of a famous paleosol stratigraphy known as “Jules Soil” in the region, which cannot be missed during the Quaternary mapping project in the Midwestern U.S.A10. The “Jules Soil” is dated at ~16 14C ka, equivalent to 19 cal ka BP10,11,12. Although the early-stage 14C dating uncertainty is large, the timing of the “Jules Soil” formation falls clearly in the 20.5-17.8-ka window.

However, it is notoriously difficult to constrain any type of large-scale flood deposit due to the lack of datable materials. It is not surprising that the timing of the late LGM warming-induced meltwater floods has been debated for nearly a century. Ekblaw and Athy13 hypothesized that a large meltwater flood occurred near the southernmost margins of the LIS, known as the “Kankakee Torrent”, in the Kankakee River valley, a second tributary of the Illinois River (Fig. 1b). The meltwater outwash terraces can be found ~30 m below and ~20 m above the modern floodplain surface in the middle Illinois River valley and are traditionally attributed to the “Kankakee Torrent”14,15. The high outwash terrace can be traced to a massive outwash floodplain that ran down the Illinois River where the valley broadens from 5 km to > 30 km, implying the release of a substantial amount of meltwater (Fig. 1a). To date, the minimal age of the “Kankakee Torrent” has been determined to be 18.93 ± 0.11 ka BP from 4 AMS 14C dates on terrestrial tundra plant leaves archived in a small lake basin in the Oswego channel that was formed by overflow scour through the Marseilles Moraine system (Fig. 1b)15,16. Here, we provide a probability density analysis of accurate AMS14C dates derived from terrestrial fossil plants in the Oswego Channel and Bryn Mawr lacustrine deposits with other representative and reasonable 14C dates (Fig. 1a; Supplementary Table 2) to synthesize that the meltwater event occurred specifically within the 19.4−18.3-ka-IntCal20 interval (Fig. 2b). Again, the AMS and conventional 14C dates clearly link the meltwater release timing to the late LGM warming phase.

Ice volume changes

The consequent question is the extent to which the ice volume changed when the LIS retreated during the late-LGM warming phase17,18. Did meltwater released from the decaying LIS change the AMOC? These important concerns must be reconciled before discussing intra D-O/Heinrich-scale climate changes over the NH. In the midwestern and northeastern sectors of the U.S.A., it is historically understood that the Lake Michigan, Green Bay, Huron, Erie, Ontario, and Miami-Scioto ice sheet lobes all decayed from ~20−18 ka19,20,21,22,23. These lobes decayed slightly faster than the James, Des Moines, and Chippewa ice sheet lobes in the Great Plains region in the western sector of the U.S.A.3,24−26. The 10Be exposure dates measured on the Green Bay, Chippewa, and New England lobes along with the 14C chronology on the James, Des Moines, Lake Michigan, and Miami-Scioto lobes have indicated that the entire LIS outer margins retreated substantially from 20 to 18 ka with a peak retreat rate at ~19 ka20,25. A large ice volume appears to have melted and drained into the ocean at this time. Here, we estimate the ice volume changes of the LIS and their equivalent contribution to global sea level rise by using a geophysically constrained ICE-6G reconstruction of the LIS, in which Global Positioning System measurements of vertical crust motion are applied to constrain the local ice cover thickness17. The results show that from 22.5 to 18.0 ka, the entire southern margins of the LIS decreased from 9.3 to 7.7 m sea-level equivalent density volume (Supplementary Table 3). The ICE-6G reconstruction estimates a maximum sea level contribution of ~1.6 m from the decaying LIS southern margins during the late LGM phase at ~20.5−18 ka (Supplementary Table 3; Figs. 2c-2d, 2f). Apparently, such a small contribution to global sea level rise would not shut down the AMOC but would instead sustain it27,28, which is consistent with the strong AMOC mode observed from ~20−18 ka29,30.

Paleotemperature construction

Reconstructing annual temperature variations near the southern LIS margins is important for discussing the climate transitions from the cold to warm LGM phases and to the less-cold early HS1 phase. Such reconstructions can help researchers assess whether North American air temperature changes were more sensitive to radiation forcing than the global average temperatures were during the broad sense of the LGM due to the presence of the world’s largest ice sheet31. Here, we analyze branched glycerol dialkyl glycerol tetraethers (brGDGTs) obtained from the loess-paleosols exposed at the classical Keller Farm site in the central USA. The Keller Farm loess deposits, located ~120 km from the southern margins of the LIS at 38°45′ N; 90°00′ W, show 13 intra D-O/Heinrich-scale climate units (Supplementary Fig. 1)5. Because changes in brGDGTs are functions of the growing-season soil temperature, which is influenced by the standing plant coverage32,33,34, the Keller Farm loess−paleosol brGDGTs should indicate changes in nonwinter soil temperatures. The results reveal that the cold early LGM phase from ~22-20 ka had an average growing-season soil temperature of 10.4 ± 0.7 °C with the lowest temperature observed at 9.5 °C, the warm late LGM phase from ~20-18 ka had an average growing-season temperature of 11.8 ± 1.1 °C with the highest temperature at 13 °C, and the cool early HS1 phase from ~17.8-16.0 ka had an average growing-season temperature of 11.4 ± 0.8 °C with the lowest temperature at 11 °C (Fig. 2e). In fact, the Keller Farm loess sequence clearly shows physical changes in climate stratigraphy from a coarse loess layer to paleosol horizon and to a less-weathered loess unit, thus independently indicating multiple instances of climate-mode shifts between warm/wet and cold/dry conditions near the ice margins during the study interval (Fig. 2e). Because the brGDGT-based soil temperatures vary 0−6 °C higher than air temperatures, depending on the different proportions of vegetation coverage34, the air temperatures near the southern LIS margins could vary in relatively smaller ranges due to the low vegetation cover under generally relatively cold/dry glacial conditions. The variance in growing-season soil temperatures with a maximum of 3.5 °C (Fig. 2e) in the early to late LGM at places ~120 km from the southern margins of the LIS (Fig. 1) potentially reflects a greater change in atmospheric temperatures than in soil temperatures. Importantly, as global mean surface temperatures varied by only ~0.5 °C during the LGM31, our reconstructed paleotemperature variations suggest that North America was more sensitive to radiative forcing than the rest of the world at this time.

Atmospheric response to LIS retreat

Although late LGM warming did not release sufficient meltwater to change the AMOC behavior, we argue that the retreat of the LIS changed the ice sheet morphology (Fig. 2f), causing it to shift westerly from zonal flow to meridional flow across the North Atlantic, transporting additional heat/moisture northward35. Here, we compile unpublished data from a fully coupled ocean-atmosphere climate model to examine how North American westerlies responded to a reduction up to 20% in the LIS under LGM climate boundaries35. To apply a LIS reduction to its maximum scenario was based on 2 observations. (1) the global ice sheets were relatively thin at ~21 ka6 (20.5 ka in this study), and (2) the outermost LIS margins at Bloomington, Illinois, at 39.4 °N retreated to Milwaukee, Wisconsin, at 43.4°N by at least 400 km. We estimate that 900,000 km2 ice margins melted away between the late HS2 and late LGM phases (Fig. 2f). The experiments (LIS_0.8 versus full LGM) indicate that North American westerlies shifted up to 65 °N meridionally, transporting heat from the subtropical west to subpolar northeast of the Atlantic Basin and modulating subpolar and subtropical gyre systems (Fig. 3). The results also reveal that anomalous southerly-westerly winds in response to a maximum reduction in LIS height were intensified more in winter (Fig. 3a) than in summer (Fig. 3b), indicating that rising boreal summer insolation also fueled winter atmospheric circulations, thereby decreasing the winter sea ice cover in the northeast North Atlantic. Anomalous southern westerly winds brought heat and moisture from the Gulf of Mexico northward to the southern LIS margins, and this heat and moisture provided the energy required to further decay the LIS margins (Fig. 3c).

Fig. 3: Wind anomaly fields at an 850-hPa geopotential height.
figure 3

It is simulated by 80% minus 100% of the LIS height under the full LGM boundaries in (a) January, (b) July and (c) annual mean. Vectors (m/s) indicate anomalous wind fields. Shading represents the zonal wind anomaly. The black bold contour is the ice sheet margin. It appears that zonal wind fields over North America and the North Atlantic (160-0 °W; 30-70 °N) show that westerlies (red) flowed meridionally across the North Atlantic, delivering heat from the west to the Northeast Atlantic; we propose that this mechanism melted the southern EIS margins and thinned the ice sheet. Meanwhile, anomalous southerly winds delivered heat from the Gulf of Mexico to the LIS, which we propose facilitated the melting of the southern margins of the LIS and the thinning of the ice sheet in general. The simulation identifies westerlies in the 55-70 °N and easterlies in the 40-55 °N (blue) anomaly fields, implying that the full-height LIS boundaries during the early LGM could also sustain two separate atmospheric-oceanic systems (see text).

Full size image

Discussion

In phase with the advance, retreat, and readvance of the LIS in the early and late LGM and early HS1 phases, respectively, the high-latitude boreal temperature variations, as indicated by the Greenland ice δ18O record, clearly indicate cold-warm-cool climate patterns on the intra D-O/Heinrich scales (Fig. 4a)36. Our soil temperature reconstruction in the loess-paleosol sequence near the LIS margins confirms these climate patterns (Fig. 4a). The European Ice Sheet (EIS) also advanced in the early stage and retreated as meltwater was released into the North Atlantic in the late LGM stage before readvancing in the early HS1 phase (Fig. 4a)37. The key difference is that the LIS advanced to much lower latitudes than the EIS, which could have played an important factor in global climate feedbacks. On the Asian continent, records, despite large chronology uncertainties, show that Tibetan Mountain glaciers advanced and retreated in the early and late LGM, respectively38,39,40; by the late LGM, 36,000 km2 pan-Tibetan Lakes had recharged to a high level (Fig. 4a)41,42,43,44. The Indian summer monsoon rainfall intensity in the Indian subcontinent45, south of the Tibetan Plateau, southwest China46,47, and the composite record of the East Asian summer monsoon intensity from the Hulu, Sanbo, and Dongge caves in southeast China48 all demonstrate a pervasive warm climate chronozone during the late LGM (Fig. 4a). Moreover, it has long been recognized that early humans in southern Europe occupied LGM niches that were as warm as they are at present, interpreted as LGM refugia49,50. High-resolution pollen records of Alborán Sea core MD95-2043 confirm that the landscapes south of the EIS shifted from semidesert shrublands to forests in response to the climate transition51. Circulations in the western Pacific52,53 and tropical Indian Ocean54 seem to have shifted from weak to strong and then weak again across the LGM toward the early HS1 phase. The Lake Qinghai record on the northeastern Tibetan Plateau55, stalagmites along the lower tributaries of the Yangtze River56, stalagmites and maximum shoreline elevations in New Mexico, southern U.S.A57,58, and stalagmites of Middle East Soreq cave59 records all show monsoon and westerly rainfall intensity shifts from weak to strong and then weak again throughout the LGM to the early HS1 (Fig. 5). We find compelling records that late-LGM warming widely prevailed over the NH, and many proxies indicate that this warming was sandwiched by early LGM and early HS1 cooling records.

Fig. 4: Northern Hemisphere (NH) climate records for LGM phase changes.
figure 4

a NH climate records for very-cold early LGM, fairly-warm late LGM, and a less-cold early HS1 phases: GISP2 δ18O for boreal temperatures36; paleosol GDGT temperatures near the southern LIS margins (this study); Channel River εNd index indicating European Ice Sheet retreat in LGMa (5-point moving average on younger age variations)37; Tibetan pan-lake level; Indian BT cave δ18O45; lake Tengchongqinghai (TCQH) leaf-wax δD record, showing the Indian summer monsoon rainfall intensity in SW China47; and composite of Hulu, Sanbo, and Dongge cave δ18O records indicating East Asian summer monsoon properties48. b Mechanism: global sea level4; global CO258; AMOC proxy of 231Pa/230Th ratio27,28; north‒south migration index of North American westerlies5; GDGT temperatures of loess-paleosol sequences near LIS margins; advance-retreat distance of the Lake Michigan lobe (LML) along the southern LIS margins; and summertime (red) and wintertime (blue) insolation at 40° N75. The large dots label low summer insolation in LGMb, increasing summer insolation in late LGM, and decreasing winter insolation in HS1b, along with the associated climate highlights.

Full size image
Fig. 5: World speleothem δ18O records.
figure 5

a Carbonate content of Lake Qinghai on the NE Tibetan Plateau, a proxy of the East Asian monsoon intensity over the Eastern Hemisphere; the findings suggest that more precipitation fell in the LGMa phase than in the LGMb and HS1b phases55; (b) Hulu cave δ18O records (PD + MSD), a proxy of the East Asian Summer Monsoon behavior over the lower Yangtze River valley in East Central China; the findings suggest intensified summer monsoon circulation during LGMa56; (c) fort Standon cave, New Mexico δ18O record, a proxy of the Westerly Jet Stream in North America over the Western Hemisphere57, indicating warm/wet signals during the LGMa phase, fully agreeing with the maximum Lake Estancia level historically recorded in New Mexico58; (d) Soreq cave, Jerusalem δ18O record, a proxy of Mediterranean rainfall behaviors in the Mideast, also revealing a relatively warm LGMa phase59. All records show e ither a tendency of substantial warming during the LGMa phase or a warm phase sandwiched by two cold phases.

Full size image

Greenhouse gas forcing played a limited role in late LGM warming because the atmospheric CO2 concentration increases after the late traditional LGM, while its acceleration seems counterintuitive to the early HS1 cooling (Fig. 4b)60. Second, the limited global sea-level rise4 with a strong AMOC throughout the LGM27,28 suggests that the AMOC played some role in late-LGM warming but not in early LGM cooling (Fig. 4b). Our estimation of the small contribution of the 1.6-m global sea-level rise resulting from LIS decay agrees consistently with these records. Last, the impact of boreal solar insolation is also a paradox: early traditional LGM cooling corresponded to low summer but high winter insolation, while late-LGM warming corresponded to increasing summer but decreasing winter insolation and early HS1 cooling corresponded to high summer but low winter insolation. Here, we focus on the evolution of the global westerly circulation coupled with LIS fluctuations to propose that in the early LGM phase, the equatorward shift of North American westerlies and the LIS advance amplified the low summer insolation driving NH cooling (Fig. 2f; 3; 4b). In the late LGM phase, the poleward shift of westerlies and the LIS retreat with a strong AMOC amplified the small increase in summer insolation, steering NH warming. In the early HS1 phase, the equatorward shift of westerlies and the readvance of the LIS with a weak AMOC, against increasing summer insolation and rising atmospheric CO2, collectively amplified the small decrease in winter insolation, thereby driving NH cooling (Fig. 4b). We predict that the very cold and fairly warm early and late LGM climate phases were characterized by frequently occurring cold and warm summer anomalies, respectively, while the less-cold HS1 phase was characterized by frequently occurring cold winter anomalies; this finding needs future investigation and testing.

The underlying mechanism for this hypothesis is that the advance/thickening and retreat/thinning of the LIS were sensitive to NH thermal forcing because of the position of its southernmost margins at the relatively low latitude of 39.4 °N. Specifically, the advancing LIS resulted in formidable obstacles for the meridional migration of subtropical jet streams5,31,61,62,63,64,65,66,67, thus enhancing the southward bulging of cold polar jet streams68, expanding the sea ice on the northern North Atlantic, shifting the Intertropical Convergence Zone equatorward, increasing the surface albedo, and reducing shortwave absorption31,69,70. Fast-moving and zonally flowing cold/dry westerlies facilitate the North American-Eurasian atmospheric teleconnection6 by reorganizing atmospheric pressure systems, AMOC modes, and Rossby waves5,31,68,69,71,72. The polar jet streams with potential vorticity over Eurasia also migrated, drifted, and bulged equatorward, causing cold air to sink along NW coastal Europe and move radially over the Tibetan Plateau. The EIS southern margins advanced to as low as 52 °N, although this latitude was much higher than that reached by the LIS37. Tibetan Glaciers advanced to lower elevations at the same time38,39,40, while the troughs of Asian jet streams likely invaded equatorward, thus hindering meridional heat and moisture transport around the Tibetan Plateau and stagnating the Indian and East Asian Summer Monsoon rainbelts over relatively low middle latitudes.

The retreat of the LIS reduced the presence of blocking barriers, allowing warm subtropical jet streams to shift poleward and bringing warm rainfall from the Gulf of Mexico and Pacific and Atlantic Oceans to the ice margins (Fig. 3)5,8,31,63,64,65,66,67,73. Fast-moving and meridionally flowing warm/wet westerlies reduced sea ice in the upper North Atlantic, shifted the Intertropical Convergence Zone southward, decreased albedo, and increased shortwave absorption5,31,69. The propagation of meridionally flowing subtropical jet streams across the North Atlantic corresponded to EIS retreat (Fig. 3)31,37. Because westerly poleward displacement created a “vacuum” spatiotemporally over relatively high middle latitudes, the warm/moist monsoon circulations from the Indian and West Pacific Oceans rushed in and produced persistent warm rainfall events74. Tibetan Mountain glaciers begin to decay, and the pan-Tibetan Lake level over 36,000 km2 rose to 60–100 m44. We argue that these intra D-O/Heinrich-scale climate oscillations represented the NH climate system in response to global atmospheric boundary75 and LIS morphology changes after the world’s largest ice sheet, the LIS, began to decay. Importantly, the LGM climate shifted from very cold to fairly warm phases over many NH regions, which should be highlighted in future studies.

Methods

AMS 14C dating of fossil plants, wood fragments, and alpha-cellulose

Plant fossils were picked from washed residue in cores and identified under a trinocular dissecting scope. Digital images were made of all the macrofossils chosen for AMS 14C dating (Supplementary Fig. 1). To prepare the samples for AMS 14C dating, a weak acid (0.2 M HCl) was used to pretreat the leaf, needle, and wood fragments. The same pretreatment was also applied to laboratory reference materials, including ISGS 14C-free wood, IAEA C5 (Two Creek Forest) wood, and FIRI-D (Fifth International Radiocarbon Inter-Comparison D) wood. Alpha-cellulose was extracted from whole wood using Soxhlet extraction methods. Each sample was first boiled in DI water for 6 hours to remove soluble sugars. Then, the sample was extracted to holocellulose, removing lignin and proteins using an acidification/bleaching process that utilized glacial acetic acid and sodium chlorite in a water bath at 70 °C under sonication. The acidification and bleaching process continued for 6 h with sequential acetic acid and sodium chlorite additions each hour. The sample was then thoroughly rinsed with DI water before the holocellulose was converted to alpha-cellulose by a reaction with 10% (w/v) NaOH in a water bath at 80 °C and sonication for 45 minutes followed by rinsing and a 17% (w/v) NaOH treatment in a water bath at room temperature for an additional 45 minutes. Each sample was thoroughly rinsed with DI water, 1% HCl solution, and finally more DI water until the effluent had a neutral pH. The sample was white in color at this time. Approximately 1–5 mg of dried sample, deadwood background and secondary working wood standards, and approximately 35 mg NIST Oxalic Acid II primary standard were placed into preheated quartz tubes for sealed combustion at 800 °C with a minimal amount of preheated CuO granules. Purified CO2 was collected cryogenically and submitted to the Keck Carbon Cycle AMS Laboratory at the University of California-Irvine for analysis using the hydrogen-iron reduction method on prepared graphite. All results were corrected for isotopic fractionation. The sample preparation backgrounds were subtracted based on the measurements of 14C-free wood blanks. The results of the IAEA C5 and FIRI D working standards fell within 1σ of the consensus values. This method has been consistently applied to all samples processed in the ISGS Radiocarbon Dating Laboratory, including those applied in previously published works and discussed herein.

Geophysically constrained ICE-6G reconstruction of the LIS ice volume

The ice thickness is multiplied by the grid cell area and then summed over the deglaciation region from 23 to 17 ka with 1.0- and 0.5-ka time slices. The total contribution to sea level rise was calculated as the LIS ice volume divided by the surface area of the global ocean (3.61 × 1014 m2) and corrected by the freshwater density (Supplementary Table 3). We then used a simple square mask to calculate the ice volume on the southern margin within the southern topographic sector of the ice sheet (Figs. 2c-2d, 2f), which allowed for the calculation of the LIS volume change at the southern margins throughout time while excluding ice volume changes occurring in other sectors of the LIS.

GDGT analysis

As described in a previous study34, pulverized loess and paleosol samples of approximately 15 g were treated with dichloromethane (DCM): methanol (9:1) using an accelerated solvent extractor (ASE 350) at 100 °C and 1500 psi in 3 cycles with 5-min heating followed by 5-min static extraction. The extracts were dried under an N2 environment and isolated with a silicone gel column with 25 ml of DCM:methanol (1:1, V/V). The dried extracts were then redissolved in hexane:isopropanol (99:1, V/V) and filtered over a 0.45-μm polytetrafluoroethylene filter. The brGDGTs were analyzed using a Shimadzu liquid chromatography triple quadruple mass spectrometer (LC‒MS 8030) with an autosampler and Labsolutions manager software. Detection was achieved in a chemical ionization chamber under atmospheric pressure with ion monitoring at m/z 1050, 1048, 1046, 1036, 1034, 1032, 1022, 1020, and 1018. The brGDGTs were quantified from integrated peak areas of the [M + H]+ ions and compared with the C46 internal standards. The source conditions were as follows: interface voltage 4500 V, interface temperature 350 °C, drying gas (N2) flow 5 L/min, Neb gas flow 2.5 L/min, and heat block temperature 250 °C. The injection volume was 50 μL. The separation of brGDGTs was achieved with two coupled Inertsil silica columns (250 mm × 4.6 mm × 3 μm). GDGTs were separated isocratically for 85 min with 95% n-hexane and 5% isopropanol at a flow rate of 0.6 mL/min. After each analysis, the column was cleaned by flushing using 10% n-hexane/90% isopropanol for 20 min. To derive land temperatures from brGDGT distributions in loess, the mean annual temperature (MATmr) was calculated over the past 24 ka using the following equations according to the most recent proxy calibration as follows:

$${\rm{MATmr}}=7.17+17.1\times \left[{\rm{Ia}}\right]+25.9\times \left[{\rm{Ib}}\right]+34.4\times \left[{\rm{Ic}}\right]-28.6\times \left[{\rm{IIa}}\right]\left({\rm{n}}=222,{{\rm{r}}}^{2}=0.68\right)$$

where the fractional abundances of Ia, Ib, Ic, and IIa are relative to the sum of all brGDGTs. Owing to the very low concentrations of brGDGTs in some samples, MAT was calculated using the 3 most abundant and omnipresent brGDGTs, Ia, IIa, and IIIa, in the form of the multiple linear regression simple index (MATmrs) over the past 24 ka as follows:

$${\rm{MATmrs}}=5.58+17.91\times [{\rm{Ia}}]{{\mbox{-}}}18.77\times [{\rm{IIa}}]({\rm{n}}=231,{{\rm{r}}}^{2}=0.62)$$

where the fractional abundances of Ia and IIa are relative to the sum of Ia, IIa, and IIIa. Over the past 24 ka, the two calibrations give very similar results. The average analytical reproducibility of the MATmrs based on duplicate injections of a selected set of loess–paleosol samples is ca. 0.5 °C. A laboratory internal standard (lacustrine sediment) was also injected after every 30 samples to check the repeatability of the sample test, and the average analytical reproducibility of the MATmrs was <0.3 °C.

COSMOS (ECHAM5-JSBACH-MPIOM) model

Here, we use a comprehensive fully coupled atmosphere–ocean general circulation model (AOGCM), COSMOS (ECHAM5-JSBACH-MPI-OM), to mimic climate change caused by the lowering of the LGM Laurentide ice sheet (LIS) by up to 20%. The atmospheric model ECHAM576, complemented by the land surface component JSBACH77, was used at the T31 resolution (∼3.75), with 19 vertical layers. The ocean model MPI-OM78, including sea-ice dynamics formulated using viscous-plastic rheology79, has a resolution of GR30 (3 ×1.8) in the horizontal direction, with 40 uneven vertical layers. This climate model has already been used to investigate a range of paleoclimate phenomena spanning glacial-interglacial climate backgrounds31,80,81,82,83.

Regarding the LGM control experiment, the external forcing and boundary conditions were adopted from the PMIP3 protocol for the LGM (available at http://pmip3.lsce.ipsl.fr/). The respective boundary conditions for the LGM comprise orbital forcing, greenhouse gas concentrations (CO2 = 185 ppm; N2O = 200 ppb; and CH4 = 350 ppb), ocean bathymetry, land surface topography, runoff routes according to PMIP3 ice sheet reconstructions and increased global salinity ( + 1 psu compared to the modern value) to account for a sea level drop of ~116 m. Regarding the sensitive run, we kept all the boundary conditions identical to the LGM conditions while reducing only the LIS by 20% (LIS_0.8) to quantify the impacts of the lowering of the southern LIS margins on the northern westerlies.

Seasonal insolation changes

To assess the average summer insolation, based on a previous study75, we calculated this term in the low range of 464 W/m2 during the early LGM, 466 W/m2 during the late LGM (+0.4%), and 476 W/m2 during the early HS1 (+2.2%) at 40 °N. The average winter insolation was calculated in the range of 185 W/m2 during the early LGM, down to 181 W/m2 during the late LGM (−2.2%) and to 176 W/m2 during the early HS1 (−2.8%) phases. During the late LGM, small-scale rising summer insolation was modeled to trigger deglacial warming30,84,85,86,87,88,89,90,91, but at the same time, winter insolation decreased by a relatively high magnitude. During the early LGM, the relatively low summer insolation and relatively high winter insolation corresponded to Northern Hemisphere cooling, while during the early HS1 phase, the further increasing summer insolation and decreased winter insolation with a similar magnitude corresponded to Northern Hemisphere cooling (Fig. 4b).

North American westerly axis proxy

Climate stratigraphy of loess − paleosol, dune − wetlands marl − lacustrine records and LIS margin advance − retreat records along the retreating trajectory of the ice sheet were selected to develop the NA westerly proxy5. This was because the poleward migration of subtropical westerlies brought warm rainfall from the Pacific Ocean, Atlantic Ocean, and Gulf of Mexico toward the ice sheet margins, facilitating soil and fresh watershed formation, while the equatorward migration of polar westerlies favored the formation of loess, dune sands, and salty watersheds. A 13 unit of intra D-O/Heinrich-scale climate stratigraphic framework was established using the composite records, and the longest loess-paleosol succession was measured on the L* value variations with 1σ uncertainty to indicate a shift in the westerly NA axis5.