Research ArticleCLIMATOLOGY

Western U.S. lake expansions during Heinrich stadials linked to Pacific Hadley circulation

See allHide authors and affiliations

Science Advances  28 Nov 2018:
Vol. 4, no. 11, eaav0118
DOI: 10.1126/sciadv.aav0118


Lake and cave records show that winter precipitation in the southwestern United States increased substantially during millennial-scale periods of Northern Hemisphere winter cooling known as Heinrich stadials. However, previous work has not produced a clear picture of the atmospheric circulation changes driving these precipitation increases. Here, we combine data with model simulations to show that maximum winter precipitation anomalies were related to an intensified subtropical jet and a deepened, southeastward-shifted Aleutian Low, which together increased atmospheric river–like transport of subtropical moisture into the western United States. The jet and Aleutian Low changes are tied to the southward displacement of the intertropical convergence zone and the accompanying intensification of the Hadley circulation in the central Pacific. These results refine our understanding of atmospheric changes accompanying Heinrich stadials and highlight the need for accurate representations of tropical-extratropical teleconnections in simulations of past and future precipitation changes in the region.


The winter precipitation–dominated portion of the southwestern United States is a water-stressed region in which ongoing climate change is projected to reduce water availability (expressed as PE, precipitation minus evaporation) in the coming decades (1). Substantial uncertainties in these projections remain, however, because of uncertainty in the circulation responses to future climate change (1, 2). For over a century, ancient shoreline deposits in hydrologically closed basins within the U.S. Great Basin, spanning southern California through northern Utah, have been understood to indicate that the region experienced much wetter conditions in the past (3, 4). Water budget analyses have identified that these lake highstands required substantially increased precipitation (a factor of ~1.7 to 2.4 relative to the present) in addition to reduced evaporation (57). These marked past changes in precipitation provide an important opportunity for identifying the atmosphere’s response to past climate changes, testing whether climate models produce realistic changes in precipitation when run with paleoclimate boundary conditions (811), and better understanding the dynamical mechanisms involved in present and future precipitation changes in the region (12).

Because of the correspondence of lake highstands with the end of the last glacial period, early work to understand precipitation maxima in the southwestern United States focused on the impact of North American ice sheets on the jet stream (13, 14). More recent work has highlighted the importance of both circulation and temperature changes associated with ice sheets in increasing moisture availability in the Great Basin during glacial periods (11). Improved chronologies for both ice extent and lake highstands, however, have demonstrated that the wettest conditions in the Great Basin did not occur during the time of maximum North American ice extent [the Last Glacial Maximum (LGM), ~19 to 23 thousand years (ka)]. Instead, most lakes reached their maximum extents during Heinrich stadial 1 (HS1; ~18 to 14.7 ka) (15, 16), a period of intense winter cooling of the North Atlantic and a southward shift of the intertropical convergence zone (ITCZ) (17, 18) at the beginning of the last deglaciation. As reviewed below, evidence indicates that previous Heinrich stadials were also marked by increased precipitation in the Great Basin, suggesting a consistent pattern.

Previous studies have suggested a range of changes in atmospheric circulation that could have increased precipitation in the southwestern United States during Heinrich stadials, including a southward shift of the upper tropospheric North Pacific jet (19), increases in jet intensity (10, 20, 21), deepening of the Aleutian Low (22), and increased summer precipitation (23). These results paint a confusing picture, with each pointing to different aspects of the North Pacific atmospheric circulation. Moreover, these studies disagree as to whether the proposed circulation changes were driven by high-latitude cooling directly (19), changes in the Pacific ITCZ and Hadley circulation (21), or teleconnections to tropical Atlantic convection (24).

One factor in this disagreement is that simulating the correct dynamical response to the forcings during Heinrich stadials remains a challenge for coarse-resolution global general circulation models (GCMs). For example, as discussed below, the TraCE-21ka simulation of the last deglaciation (which attempts to reproduce the transient climate response to changes in greenhouse gases, ice sheets, and Earth’s orbit) does not show a substantial precipitation increase in the Great Basin during HS1 with respect to the LGM. To shed light on the dynamics driving precipitation increases in the southwestern United States during Heinrich stadials, we draw on an updated compilation of regional data, dynamics observed in modern interannual variability, and an ensemble of climate model simulations under a range of forcings to provide a physical mechanism that integrates these previous results and clarifies the drivers and patterns of precipitation and atmospheric circulation changes in the southwestern United States during Heinrich stadials.


Western U.S. precipitation changes during Heinrich stadials

We compiled Great Basin lake highstand ages (i.e., estimates of the times at which lakes reached their maximum extents) from the last deglaciation to provide constraints on the mechanisms of precipitation increase during HS1 (Fig. 1 and Supplementary Text). Highstands occurred synchronously along southwest-northeast trends and progressed through time from southeast to northwest. Basins in the southwest, center, and northeast of the Great Basin attained their highstands from 16.0 to 17.5 ka, with many clustering around 16.5 to 17 ka (15). To the northwest of this southwest-northeast band, Lake Russell and Lake Lahontan reached their highstands slightly later, at 15.6 to 16.0 ka (15). Moving farther to the northwest, the highstand in the Chewaucan Basin (southeastern Oregon) occurred after 14.6 ± 0.3 ka (25). For basins with well-documented LGM lake levels, these deglacial highstands represent 49 to 82% increases in surface area relative to the LGM (Supplementary Text).

Fig. 1 Extents and ages of lake highstands in the U.S. Great Basin during the last deglaciation.

(A) Lakes at their greatest extents of the last glacial cycle, with colored circles denoting timing of wettest conditions during the last deglaciation (Supplementary Text). Blue arrow shows inferred direction of anomalous moisture transport during HS1. Gray line is the outline of the Great Basin. Percentages indicate the magnitude of lake area increase from the LGM to HS1 (Supplementary Text). (B) Age estimates with 95% confidence intervals for wettest conditions during the deglaciation in each basin as a function of distance from a line connecting the Panamint and Bonneville basins [dashed line in (A)]. Highstand ages are progressively younger with greater distance northwest from this line. Paleolakes: B, Bonneville; Ch, Chewaucan; Cl, Clover; F, Franklin; J, Jakes; L, Lahontan; P, Panamint; R, Russell; S, Surprise; W, Waring. Lake highstand map was adapted from (16).

This pattern suggests that anomalous moisture supply was derived from the southwest and transported toward the northeast, consistent with a previous compilation of paleo-data spanning the deglaciation (23) while adding improved chronological control and newer records. The data suggest an amplification of southwesterly “atmospheric river” moisture transport identified in LGM climate model simulations (9), but highstand ages are oriented orthogonally to the northwest-southeast steering of LGM storms suggested by Oster et al. (8). Northwestward progression of lake highstands during and immediately after HS1 could reflect diminishing ice sheet topography as the deglaciation progressed (10, 20).

Data for other Heinrich stadials are consistent with the patterns identified in HS1 (fig. S1). Oxygen and uranium isotope data indicate increased precipitation in northern Utah’s Bonneville Basin during HS2 (26). Sediments in the Manix Basin of the Mojave Desert (southern California) show lake-level increases during HS3, HS4, and other stadials between 43 and 25 ka, with lower lake levels during interstadials (27). Speleothems in Arizona and New Mexico show decreases in δ18O values during stadials between HS2 and HS5, consistent with increases in winter precipitation and/or decreases in summer precipitation during these events (19, 22). Tracers of local infiltration (growth rate, trace element concentrations, 87Sr/86Sr ratios, and/or δ13C values) from stalagmites in eastern Nevada (28) and the central Sierra Nevada (29) indicate wet conditions during HS11 and HS6, respectively, followed by rapid drying at the end of each stadial. In contrast, HS2 to HS5 appear to be marked by drying in the Chewaucan Basin of southeastern Oregon (30). Although more work is certainly needed to document the region’s pre-LGM hydrological history, these findings suggest that stadials before HS1 were also marked by greater winter precipitation in the southwestern United States, with drying in the northwest Great Basin, consistent with the southwest-northeast orientation of anomalous moisture transport.

An exception to this consistency are findings from Pyramid Lake (western Nevada) and Owens Lake (south-central California), which have been interpreted as indicating dry conditions during stadials between 27 and 50 ka (31). If the chronologies for these records are accurate, then these findings would suggest that the boundary between wet and dry conditions during stadials was located farther south in stadials between 27 and 50 ka than during HS1 (27). As noted by the authors, however, chronological uncertainties in these records make the absolute phasing between lake-level oscillations and stadial-interstadial variability difficult to determine.

Simulated precipitation changes due to North Atlantic cooling

We investigate the relevant physical mechanisms using four freshwater hosing experiments performed with the fully coupled Earth system model CM2Mc (Materials and Methods) (32). Given that all models have shortcomings (particularly in their representations of atmospheric physics, orography, and Pacific mean state), we expect that the simulated response to hosing under a given set of boundary conditions will always be incorrect to some degree. An instructive example is the TraCE-21ka transient simulation of the deglaciation, which is the only publicly available simulation using relatively realistic orbital, greenhouse gas, and ice sheet boundary conditions for HS1 (33). While it has been used previously in investigations of deglacial moisture delivery to the southwestern United States (20), it shows negligible changes in regional precipitation from the LGM to HS1 (see below), suggesting that it does not capture important elements of connections between North Atlantic cooling and western U.S. precipitation. To address this hurdle, we present an ensemble of simulations for which the atmospheric CO2 concentration, ice sheet topography, and sea level are prescribed to approximate those of the LGM (32), under a range of four orbital boundary conditions. The four orbital configurations, with opposite phases of obliquity (22.5° and 24°) and precession (precession angles of 90° and 270°), produce varied background climate states before hosing, which, in turn, produce different responses to hosing. These orbital configurations bracket those of HS1, which was characterized by intermediate values of both precession and obliquity. As shown below, the different responses then provide insights into the dynamics that drive responses similar to those observed in paleoclimate data. By testing across a wider range of conditions, our multisimulation approach helps make clear the robust dynamical responses among the ensemble of simulations—even if the response may be biased under any particular set of boundary conditions.

We focus on DJF precipitation, as this is the dominant control on stream flow and lake level in the Great Basin (34, 35) and because speleothem (19, 22) and pollen (36) data are inconsistent with the suggestion that increased summer precipitation drove lake highstands (23). Annual mean evaporation changes in our experiments are negligible, consistent with expectations for LGM-to-HS1 evaporation changes (see below). All hosing experiments show DJF precipitation increases extending over the central North Pacific toward western North America, reaching maximum values at the coast (Fig. 2). The precipitation increase is more than threefold larger and extends farther inland in the two simulations with 270° precession angles [high Northern Hemisphere (NH) seasonality] than in those with 90° precession angles; similar relative changes are observed in annual mean precipitation. The effect of obliquity is small by comparison. We focus on the two simulations with the maximum and minimum southwestern U.S. precipitation responses, hereafter termed “Strong” (270° precession angle and 24° obliquity) and “Weak” (90° precession angle and 24° obliquity); results from all experiments are shown in the Supplementary Materials.

Fig. 2 North Pacific precipitation and atmospheric circulation anomalies in hosing experiments.

(A) Anomalies in the DJF precipitation (in mm/day; shading), near-surface (10 m) wind (in m/s; vectors), and sea-level pressure (in Pa; purple contours) for the Strong hosing experiment with respect to its control. The black box shows the area plotted in (C) and (D). (B) As in (A), but for the Weak simulation. (C and D) Insets showing anomalies over the eastern Pacific and western United States for the Strong (C) and Weak (D) simulations. Note that the shading color scale is changed between the top and bottom panels. Black box in (C) shows the area of the map in Fig. 1A.

The precipitation increase is paralleled by near-surface (10 m) westerly and southwesterly wind and water vapor transport anomalies in the eastern subtropical North Pacific that are greatest in the Strong experiment (Fig. 2; water vapor transport anomalies are not shown). These patterns suggest that increased southwestern U.S. precipitation in the Strong experiment is related to increased atmospheric moisture transport from subtropical latitudes, which may include enhanced moisture transports by atmospheric river events (9), similar to anomalies associated with high-precipitation winters today (37, 38). The southwesterly orientation of wind anomalies in the northeastern Pacific matches the spatiotemporal pattern identified in deglacial lake highstands in Fig. 1, suggesting that the simulations broadly represent the dynamics of regional precipitation changes during Heinrich stadials.

Evaluation of the magnitude of simulated versus observed precipitation anomalies

To compare precipitation anomalies in our Strong and Weak simulations with the magnitude of precipitation anomalies required to explain hydrologic changes during Heinrich stadials, we examine lake-level changes between the LGM and HS1. We emphasize that, because of the coarse model resolution, model biases (39), the idealized nature of the modeling experiments, and the lack of inclusion of feedbacks like lake-effect precipitation (40), we do not expect perfect agreement; instead, our aim is to test whether the observed precipitation anomalies in our Strong simulation represent a substantial fraction of the precipitation changes needed to explain lake highstands and thus whether our experiments may be capturing an important part of the dynamics of Heinrich stadial (HS)–related precipitation changes in the southwestern United States.

As documented in the Supplementary Text, there are four Great Basin lakes with well-resolved lake-level histories spanning the LGM and HS1: Franklin, Lahontan, Russell, and Surprise. These records show lake area increases from the LGM to their HS1 highstands ranging from 49 to 82%. We leave out Bonneville despite its well-documented lake-level record because it was no longer hydrologically closed at its highstand.

We use these lake area increases to estimate the required precipitation changes using a standard steady-state closed-basin lake water balance equation (6, 4144)Embedded Image(1)where P is the precipitation averaged over the basin, AB is the area of the basin, AL is the area of the lake, EL is the evaporation from the lake, and k is the fraction of precipitation that falls in the basin outside of the lake that enters the lake as river discharge or groundwater, also known as the “runoff ratio.” The equation relates water inputs to the lake (on the left) to water outputs as evaporation (on the right), assuming that groundwater losses and changes in groundwater storage are small. Following Budyko’s framework (45, 46), k is a function of P, potential evapotranspiration (PET), and an exponent ω that determines the nonlinearity of runoff increases with increasing moisture availability (Supplementary Text)Embedded Image(2)

Following (43), we rearrange Eq. 1 to solve for the ratio of the area of the lake to the area of the basinEmbedded Image(3)

We now consider steady-state areas for a lake in two different climates, AL1 and AL2. By dividing AL2 by AL1, AB cancels and we obtainEmbedded Image(4)

To simplify the analysis, we make two assumptions. First, we assume that changes in evaporation between the LGM and HS1 are small: EL1EL2. This assumption is warranted on the grounds that decreases in NH temperature between the LGM and HS1 occurred primarily in the winter, with only small temperature changes in summer, the primary evaporation season in the western United States (47). Consistent with this assumption, our simulations show negligible (+3% to −2%) changes in annual mean evaporation in the southwestern United States in response to hosing. We relax this assumption and explore its influence in the Supplementary Text.

Second, we assume that EL ≈ PET. PET is meant to express the evaporative losses from a vegetated surface with unlimited water availability, and some expressions for PET are equivalent to evaporation from a lake surface (48). Offsets between EL and PET may occur seasonally because of heat storage in the lake or water-advected heat (e.g., due to seasonal runoff inputs), but in the annual average, these can be neglected (49). EL and PET may also differ in that many ways of calculating PET incorporate the conductance of water through vegetation. We explore the influence of varying EL/PET in the Supplementary Text but, for now, assume that EL ≈ PET.

The benefit of these assumptions is that the ratio of the lake areas can be expressed as a function of only three variables: P2/P1, PET/P1, and ωEmbedded Image(5)whereEmbedded Image(6)andEmbedded Image(7)

Model simulations included in the PMIP3 (Paleoclimate Modeling Intercomparison Project, version 3) ensemble indicate that PET/P in the southwestern United States was between 1 and 2 at the LGM [see (43)]. As we will demonstrate, this term has very little influence on the change in lake area between LGM and HS1. In the canonical “Budyko curve,” ω = 2.6, but there is a wide variability across the United States (46). Values of ω calculated from LGM simulations for the southwestern United States [see (43)] and for our control experiments (Materials and Methods) center around 1.8. Here, we report results for ω varying between 1.5 and 2.6.

In Fig. 3, we plot AL2/AL1 as a function of P2/P1, with different lines drawn for varying PET/P1 and ω. Note that PET/P1 has a negligible influence on the relationship between precipitation changes and lake area changes, while the choice of ω has a substantial influence. For values of ω similar to those estimated for our simulations (1.8), the precipitation change seen in the Strong simulation produces lake area changes that are up to 80% of the change required to satisfy the range of estimated LGM-to-HS1 lake area changes. With higher values of ω (>2.1), the Strong simulation predicts precipitation changes that are sufficient to satisfy the lake area changes. Thus, the dynamics captured by these simulations are potentially of the appropriate magnitude to explain the observed changes, despite the low resolution of the model and consequent biases. For all values of ω, precipitation changes in the Weak simulation produce lake area changes that are much smaller than observed.

Fig. 3 Simulated lake area changes (AL2/AL1) in response to annual mean precipitation changes (P2/P1).

Precipitation and lake areas with “1” subscripts represent pre-hosing values, while those with “2” subscripts represent values after hosing. The solid lines correspond to ω values ranging from 1.5 to 2.6, where ω is a basin-specific factor expressing the relationship between precipitation and runoff changes. The value estimated from LGM model simulations (43) and our control experiments is 1.8. All solid lines use PET/P1 = 1.5. The two dashed black lines are calculated for ω = 1.8 and PET/P1 = 1 and 2; note that there is almost no effect of varying PET/P1. All lines are calculated assuming no change in PET or EL and assuming that PET = EL. Bars at the top indicate precipitation changes in our Strong and Weak hosing simulations and between the LGM and HS1 in the TraCE-21ka experiment. The blue region indicates reconstructed lake area changes between the LGM and HS1 (see the Supplementary Materials).

We also examine precipitation changes from the LGM to HS1 simulated in the TraCE-21ka experiment (33). We computed precipitation anomalies for 18 to 17 ka, 17 to 16 ka, and 16 to 15 ka relative to 22 to 21 ka in the same region as for our experiments. For all three periods, precipitation anomalies relative to 22 to 21 ka are 5% or lower. As for our Weak simulation, these precipitation changes are much smaller than those required to explain observed lake area changes from the LGM to HS1, suggesting that key dynamics for western U.S. precipitation are not adequately represented in the TraCE-21ka experiment.


Dynamical connections

The wind anomalies in the central and eastern North Pacific track along the southeastern margin of a large zone of negative sea-level pressure anomalies in the northeast Pacific, which represent the deepening and southeastward shift of the Aleutian Low in the hosing simulations compared to the controls (Fig. 2). These surface changes are accompanied by an intensification and slight southward shift of the jet stream over the North Pacific, as illustrated by the positive anomalies in the zonal wind at ~300 hPa (Fig. 4 and fig. S2). For both the Aleutian Low and the jet, the changes are greater in the Strong experiment than in the Weak experiment.

Fig. 4 Central Pacific atmospheric circulation anomalies in the hosing experiments.

Anomalies in the DJF divergent wind and isobaric vertical velocity (vectors; in m/s and Pa/s, respectively) and in DJF zonal winds (shading; in m/s) for the Strong (A) and Weak (B) hosing experiments with respect to their corresponding controls for the Pacific region between 150°E and 140°W. Isobaric vertical velocities are multiplied by 100 for a better view of the circulation vectors. Note the southward-shifted and larger Hadley circulation anomalies and greater jet anomalies in the Strong simulation.

Two mechanisms can potentially explain the wintertime jet intensification and Aleutian Low changes in the hosing experiments. On the one hand, anomalous NH high-latitude winter cooling due to shutdown of the Atlantic Meridional Overturning Circulation (AMOC) increases the equator-to-pole temperature gradient, enhancing the upper tropospheric jet via thermal wind adjustment. This thermal wind mechanism dominates in the North Atlantic basin, where colder temperature anomalies are accompanied by a stronger jet stream (figs. S2 and S3). In the North Pacific, by contrast, the strongest jet anomalies are found in the Strong experiment, although high-latitude cooling and changes in the meridional temperature gradient are smaller in this simulation than in the Weak experiment. Jet anomalies also do not correlate with the latitude or strength of the jet before hosing. A previous study showed that deepening of the Aleutian Low during Heinrich stadials is also unlikely to be directly linked to high-latitude cooling, which, in isolation, would raise surface pressures throughout the mid- and high latitudes (24).

Here, we propose that the mechanism for the observed jet and Aleutian Low changes instead involves the poleward transport of angular momentum by the Pacific Hadley circulation. In the hosing experiments, the AMOC shutdown leads to an interhemispheric heating contrast with anomalous winter cooling in the NH. To enhance the northward cross-equatorial atmospheric heat transport and thereby compensate the anomalous heating contrast, the ITCZ and the Hadley circulation shift southward (50). The ITCZ shift intensifies the NH winter Hadley circulation and associated trade winds (51), both globally and in the North Pacific (Fig. 4). This strengthening enhances the momentum convergence at the latitude of the descending limb of the North Pacific Hadley circulation, thereby accelerating the subtropical jet.

In tandem with jet intensification, an enhanced divergent component of the meridional wind, directed across strong angular momentum gradients in the subtropical central North Pacific, acts as an effective Rossby wave source (52). This source alters stationary wave patterns to the north and east, leading to a deepening and southeastward shift of the Aleutian Low.

The jet and Aleutian Low responses described here bear instructive similarities to the North Pacific winter response during El Niño events. In both our Strong simulation and El Niño events, the central Pacific NH winter Hadley circulation intensifies, leading to intensification of the subtropical jet (53, 54). Likewise, El Niño events are marked by deepening and southeastward shifts of the Aleutian Low due to wave propagation related to strengthening of the central Pacific NH winter Hadley circulation (55). While Heinrich stadials differ from El Niño events in having opposite atmospheric responses in the NH and Southern Hemisphere (SH) rather than hemispheric symmetry (51, 54), both can intensify the North Pacific winter Hadley circulation (21, 53, 55), the central element in the mechanism described here.

As summarized in Fig. 5A, the responses of the jet and Aleutian Low to a stronger Pacific winter Hadley circulation increase southwesterly moisture transport into the western United States: Increased westward momentum in the jet is transported downward to the surface by midlatitude eddies, accelerating the low-level westerlies in the North Pacific (Figs. 2 and 4), and the deepening and southeastward shift of the Aleutian Low drive southwesterly surface wind anomalies and increased water vapor transport into the southwestern United States, as seen in modern interannual variability (37).

Fig. 5 Dynamical links between central Pacific ITCZ shifts and increased southwestern U.S. precipitation during Heinrich stadials.

(A) Proposed surface (black) and upper-level (gray) Pacific atmospheric circulation changes during Heinrich stadials. The southward shift of the central Pacific ITCZ (1) is accompanied by intensification of the NH winter Hadley circulation (2; streamlines), which intensifies the subtropical jet (blue contours, blue arrow) and initiates a Rossby wave response (3; dashed arrows, gray contours showing pressure anomalies) that causes the Aleutian Low to deepen and shift southeastward (black dashed contours). Together, these changes increase southwesterly moisture transport into the southwestern (SW) United States (4; black arrows). (B and C) Plots of (B) Pacific jet anomalies and (C) northeastern (NE) Pacific sea-level pressure (SLP) anomalies versus the central Pacific ITCZ shift in each experiment, showing greater jet intensification and a deepening and eastward shift of the Aleutian Low in the experiments in which the ITCZ shifts farther south. (D) Scatterplot of southwestern U.S. precipitation anomalies versus northeastern Pacific SLP anomalies, showing greater precipitation anomalies in the experiments with greater SLP responses and thus in the experiments with greater ITCZ and jet changes. Teal/brown symbols indicate precession angles of 270°/90°, and circles/triangles indicate obliquity angles of 24°/22.5°. All anomalies are calculated for DJF.

In contrast to the thermal wind adjustment mechanism related to high-latitude temperatures, the Hadley circulation mechanism explains the differences in the anomalies in the North Pacific across the hosing experiments. In the Strong simulation, the central Pacific ITCZ shifts farther south, and the Pacific winter Hadley circulation intensifies more (Figs. 4 and 5), converging more momentum into the subtropical central North Pacific (fig. S4). The strengthened circulation transfers more momentum to the jet and elicits a stronger wave response in the northeastern Pacific, explaining why the larger ITCZ shift in the Strong experiment is matched by greater jet and pressure anomalies (Fig. 5, B and C). The jet strengthening and the changes in the position and intensity of the Aleutian Low, in turn, lead to greater near-surface westerly and southwesterly wind anomalies and, ultimately, to the largest precipitation anomalies over western North America (Fig. 5D).

The hypothesized connection between the central Pacific and western U.S. lake levels is supported by recent reconstructions indicating large (~5°) southward shifts of the central Pacific ITCZ during the Heinrich stadials of the last two deglaciations (HS1 and HS11) (56, 57), substantially larger than estimates of global mean ITCZ displacements (17). Our analysis suggests that these southward shifts would have driven jet and Aleutian Low responses that, superimposed on the circulation effects of large remnant North American ice sheets, could have provided the winter precipitation increases necessary to fill basins to their highstand levels. Our results offer strong support for the hypothesis that subtropical jet strengthening related to southward ITCZ shifts drove precipitation increases in the southwestern United States during Heinrich stadials (21), although our findings indicate that the Aleutian Low response also plays a key role in the resulting precipitation changes. Our study agrees with the finding of consistent deepening of the Aleutian Low in response to hosing (24), but we link Aleutian Low changes to convection anomalies in the tropical Pacific rather than in the western tropical Atlantic (Supplementary Text).

A remaining puzzle is why the central Pacific ITCZ shifts farther south in the Strong simulation than in the Weak simulation. One potential explanation is that the central Pacific winter ITCZ begins farther north before hosing in the Strong experiment (fig. S5), likely because of low DJF insolation. Because northward cross-equatorial heat transport intensifies as the ITCZ shifts farther south (58), a southward-shifted ITCZ in the Weak control simulation may migrate a shorter distance in response to perturbations such as hosing (fig. S5). We also note that the greatest precipitation anomalies are seen in the experiment in which orbital parameters favor the greatest NH seasonality (boreal summer perihelion and high obliquity), offering a connection to observations that Heinrich stadials are characterized by extreme winter cooling and seasonal amplitude in the NH (47).

An important implication of our results is that the central Pacific ITCZ position needs to be accurately simulated to capture the magnitude of western U.S. precipitation changes in response to climate perturbations. Climate models consistently show a southward bias in ITCZ position, particularly in the Pacific basin, under present-day conditions (59). The mechanism identified here suggests that, in isolation, this ITCZ bias should cause models to overestimate winter precipitation in the U.S. Great Basin and southwest before hosing, consistent with the positive precipitation bias identified in the region in simulations of the modern climate (1). Simulations of the NH cooling may underestimate the magnitude of North Pacific atmospheric changes and western U.S. winter precipitation anomalies if the ITCZ and Hadley circulation begin with a southward bias. This persistent bias is an important motivator of our modeling approach, as the different orbital configurations produce varied Pacific ITCZ positions before freshwater forcing, which, in turn, appear tightly linked to the magnitude of the response (Fig. 5 and fig. S5). Further work should test this finding using hosing experiments in additional models.

A limitation of the model results presented here is that they do not precisely match the precipitation response from the LGM to HS1. HS1 was characterized by intermediate orbital conditions; thus, a simulation with this model using HS1-specific orbital conditions might be expected to produce precipitation changes between the Strong and Weak results, smaller than the values suggested by observations. The precessional parameter in our Strong experiment is more closely matched to conditions during HS4, HS6, and the end of HS11, and the precessional parameter in our Weak experiment is more closely matched to conditions during HS2 and HS5. This underestimation of the HS1 response might be anticipated because of the tropical Pacific biases that are persistent in this and other models; this offset motivates our exploration of a wider range of climates (and particularly Pacific mean states) to uncover the dynamics involved in producing lake-level maxima in the southwestern United States. In addition, including feedbacks related to lake-effect precipitation or vegetation changes might also be crucial to matching observation during HS1.

A second limitation is that freshwater hosing is an inherently imprecise way of simulating Heinrich stadials, in that stadials appear to begin with sea ice advances rather than with freshwater inputs (60) and actual freshwater inputs during stadials are poorly constrained. We note that this model also shows spontaneous abrupt North Atlantic coolings and weakenings of the AMOC under some boundary conditions. Precipitation responses are similar in hosed versus unhosed abrupt North Atlantic coolings in this model, although the responses are larger in response to hosing (32).

Despite these limitations, the mechanism we present is consistent with the southwesterly moisture transport suggested by lake-level reconstructions, is supported by data suggesting large southward shifts of the central Pacific ITCZ during Heinrich stadials, draws on straightforward consequences of angular momentum transports in the Pacific Hadley circulation, and bears distinct similarities to the dynamics observed during wet years in the southwestern United States today. The tropical-extratropical linkages documented here highlight the fact that changes in the North Atlantic can have global effects transmitted through the tropical Pacific, and they contribute to a global picture of atmospheric reorganizations during Heinrich stadials (18). These linkages also have important implications for future climate changes, as changes to the interhemispheric distribution of atmospheric energy inputs in boreal winter (e.g., due to changes in AMOC, sea ice, or clouds) may produce shifts in the Pacific ITCZ that would have pronounced effects on southwestern U.S. precipitation (18). Last, the mechanism presented here provides an important template for future data collection in the western United States, as the development of new well-dated hydroclimate records can offer additional tests of the spatial and temporal patterns and magnitudes of change we identify in the presently available data.


Climate model experiments and data analysis

The four sets of simulations were performed with the full-complexity coupled ocean-atmosphere CM2Mc model at a resolution of 3°. The simulations used obliquities of 24° and 22.5° and precession angles of 270° and 90° between perihelion and autumnal equinox, with 270° corresponding to maximum NH seasonal intensity/minimum SH seasonal intensity and 90° corresponding to the opposite. The eccentricity was held constant at 0.03. The hosing was applied as a uniform 0.2-sverdrup freshwater input in the North Atlantic for 1000 years, evenly distributed over the region 40° to 60°N and 60° to 12°W, similar in magnitude to freshwater fluxes used in other hosing experiments (33, 61). The aim of this input is to shut down the AMOC, as commonly done to simulate Heinrich stadials (33, 61). For each experiment, we took the mean state of the past 100 years of hosing and compared it with 100 years of the corresponding control simulation that was run under the same boundary conditions without hosing. More details about these experiments can be found in (32, 39).

Variables in the scatterplots in Fig. 5 are defined as follows: The Pacific ITCZ position is the centroid of the precipitation rate between 20°S and 20°N, zonally averaged between 120°E and 140°W; the Pacific jet stream speed is the maximum of the zonal wind at 300 hPa averaged between 150°E and 140°W; northeastern Pacific sea-level pressure was averaged over 30° to 60°N, 120°W to 180°; and southwestern U.S. precipitation was averaged over 32° to 45°N, 110° to 130°W.

Statistical significance of anomalies in Figs. 2 and 4 was calculated via a two-tailed Student’s t test, in which effective degrees of freedom and serial autocorrelation were taken into account (62). Thus, after calculating DJF means from five-day averages over 100 years in both the controls and the hosings, we applied the test to detect whether an anomaly value lies outside the 5th or 95th percentiles of a distribution given by the controls’ climatology and hence is statistically significant. For the sake of clarity in the main figures, significances are shown in figs. S6 to S8 for each experiment and variable. There, stippling (precipitation) and gray arrows (winds) indicate anomalies that are not significant at the 5% level. For experiment 24°/270°, significances could not be calculated for the zonal wind and the divergent wind (which was computed from both the zonal and meridional winds) because only monthly and yearly climatological means were properly stored. In these cases, we considered that anomalies in this experiment are significant if they are so in all three other experiments, which is a conservative estimate given the fact that responses to hosing are largest in this experiment. Significances of sea-level pressure and vertical wind changes could not be computed either, as only the monthly and yearly climatological means were archived for these variables.

Values for ω were estimated for each control experiment following (63). Briefly, annual mean values of net surface radiation (incoming shortwave and longwave radiation minus outgoing shortwave and longwave radiation) were used to compute an energetic limit for evaporation at each grid cell in the southwestern United States (E0, which is effectively EL). These values were then compared to actual evaporation and precipitation for all grid cells in which precipitation was greater than evaporation (n = 9 for each experiment). The relationships between E/E0 and P/E0 were then compared to the estimates with varying values of ω. For each experiment, the best fit was obtained with ω between 1.7 and 1.8, similar to the values estimated from LGM climate model simulations analyzed by (43).


Supplementary material for this article is available at

Supplementary Text

Table S1. Lake highstand ages used in Fig. 1.

Fig. S1. Well-dated proxy evidence of hydrologic changes during Heinrich stadials before the LGM.

Fig. S2. Temperature and wind anomalies in each hosing experiment.

Fig. S3. Relationship between jet changes and meridional temperature gradient changes.

Fig. S4. Anomalies of upper-level divergence in response to hosing.

Fig. S5. Comparison of central Pacific DJF ITCZ shift in response to hosing to its position in the corresponding glacial control simulation.

Fig. S6. DJF precipitation and wind anomalies for all four hosing experiments.

Fig. S7. DJF central Pacific zonal wind anomalies for all four hosing experiments.

Fig. S8. DJF central Pacific circulation anomalies for all four hosing experiments.

Fig. S9. Ages of terminal lake highstands in the U.S. Great Basin during the last deglaciation, including ages excluded by quality control criteria.

Fig. S10. Simulated lake area changes (AL2/AL1) in response to precipitation changes (P2/P1), evaporation changes, and differences in PET/EL.

References (6492)

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. Adams for assistance in preparing Fig. 1A and conversations about the Great Basin lake histories. We also thank W. S. Broecker, D. E. Ibarra, and two anonymous reviewers for thoughtful comments that improved the manuscript. Funding: D.M. acknowledges support from NSF award EAR-1103379. J.M. and E.M.-C. acknowledge support from NOAA award NA16OAR4310177. E.D.G. acknowledges support from Compute Canada award ayu-503 and Canadian Foundation for Innovation award 25402 and the Spanish Ministry of Economy and Competitiveness, through the María de Maeztu Programme for Centres/Units of Excellence in R&D (MDM-2015-0552). Author contributions: D.M. conceived the project, compiled paleo-data, and conducted lake water balance modeling. E.M.-C. analyzed the model output. E.D.G. conducted the model simulations. D.M. and E.M.-C. wrote the manuscript with contributions from all authors. All authors contributed to interpreting the results. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All paleoclimate data used in this study are previously published in the cited manuscripts. All model output analyzed in this study can be downloaded at All code used for the analysis of model output is available upon request from the corresponding author.
View Abstract

Navigate This Article