North-south dipole in winter hydroclimate in the western United States during the last deglaciation

During the termination of the last glacial period the western U.S. experienced exceptionally wet conditions, driven by changes in location and strength of the mid-latitude winter storm track. The distribution of modern winter precipitation is frequently characterized by a north-south wet/dry dipole pattern, controlled by interaction of the storm track with ocean-atmosphere conditions over the Pacific and Atlantic Oceans. Here we show that a dipole pattern of similar geographic extent persisted and switched sign during millennial-scale abrupt climate changes of the last deglaciation, based on a new lake level reconstruction for pluvial Lake Chewaucan (northwestern U.S.), and a compilation of regional paleoclimate records. This suggests the dipole pattern is robust, and one mode may be favored for centuries, thereby creating persistent contrasting wet/dry conditions across the western U.S. The TraCE-21k climate model simulation shows an equatorward enhancement of winter storm track activity in the northeastern Pacific, favoring wet conditions in southwestern U.S. during the second half of Heinrich Stadial 1 (16.1–14.6 ka) and consistent with paleoclimate evidence. During the Bølling/Allerød (14.6–12.8 ka), the northeastern Pacific storm track contracted poleward, consistent with wetter conditions concentrated poleward toward the northwest U.S.

cool season (Oct-Mar) precipitation (P) and Jun-Nov Southern Oscillation Index (SOI); data from ref. 34 . Compilation of hydrologically sensitive paleoclimate records covering the Deglacial interval are plotted by climate indicator type and relative wetness during the HS1b or B/A periods (see legend). Full details and references for each numbered record are shown in the Methods and Table S1. Representative records lettered (d-h) are plotted in Fig. 2: (d) -Lake Chewaucan (this study, red box shows the map area of (b)), (e) -Lake Surprise 44 , (f) -Lake Franklin 52 , (g) -Lake Lahontan 53 , (h) -Lake Elsinore 4 . (b) Chewaucan drainage system. Drainage basin boundaries shown in black, approximate modern lake extents in light blue, B/A high shoreline extent shown in dark blue, LGM extent shown in green, early Holocene extent shown in yellow. 14 C sample locations are plotted by sample type. Named areas where samples were collected are labeled (e.g. the Abert Rim). Overflow point (1338 m asl) connecting the Summer Lake sub-basin and Abert Lake subbasin shown by yellow star. Location of shell 14 C age sample from ref. 73 that was omitted from our record is shown by red star, clearly well outside the high shoreline elevation (see Supplemental Information for details). Base hillshade was created using ArcGIS 10.4 using the 10 m-resolution digital elevation model from the USGS National Elevation Dataset 90 (available from https://nationalmap.gov). snowfall 23 . Since most paleoclimate archives considered here record changes in annual-average effective moisture, we focus on modern winter climate conditions to help interpret past hydroclimate changes.
Winter hydroclimate in the western U.S. is strongly controlled by the mid-latitude westerly storm track that transports moisture from the Pacific Ocean via extratropical cyclones. Along the western margin of North America, precipitation associated with atmospheric rivers-narrow corridors of strong horizontal moisture transport located in the warm sector of extratropical cyclones 24 -contribute up to 50% of total precipitation during the cool season 25 . The central and eastern Great Basin also receives ~30% of annual precipitation from extratropical cyclones that are not associated with atmospheric river conditions during winter, as well as cutoff and closed low pressure systems during spring 26 (March-June).
At interannual to decadal timescales winter hydroclimate often displays a north-south dipole in precipitation and runoff, wherein wet (dry) winters in the southwestern U.S. (SW, hereafter, corresponding to the SW dipole region in Fig. 1a) are mirrored by dry (wet) conditions in the northwestern U.S. (NW, corresponding to the NW dipole region) in a given year (Fig. 1a) 27,28 . This pattern is not always expressed, with some years showing anomalously wet or dry conditions across the whole region, but it is the leading mode of variability in regional winter precipitation at interannual and decadal timescales 22 . Extratropical cyclones with and without associated atmospheric river conditions contribute significantly to this dipole pattern, which characterizes the distribution of both total precipitation 22 and extreme events 29 . The dipole is partly driven by Pacific and Atlantic ocean-atmosphere conditions linked to lower frequency modes of climate variability such as the El Niño Southern Oscillation (ENSO) 22,30 , the Pacific Decadal Oscillation (PDO) 31,32 and the Atlantic Multidecadal Oscillation (AMO) 33,34 . These oscillations produce global rearrangements of sea surface temperatures, atmospheric pressure patterns, and tropical convection, which perturb the northern Pacific jet stream and storm track and ultimately influence the distribution and magnitude of western U.S. cool season precipitation. For example, during the warm (cool) phase of El Niño, as well as warm (cold) PDO and cold (warm) AMO conditions, moisture transport into the western U.S. is enhanced along an equatorward (poleward) track into the SW (NW) [32][33][34] . The modulation of the dipole by the PDO, which rather than a single dynamical mode is the collective oceanic response to coupled ocean-atmosphere interactions, appears to be related to changes in oceanic baroclinicity that shift storm tracks poleward (equatorward) during the negative (positive) phase 31 and constructively enhances the ENSO teleconnections 34 .
These lower frequency teleconnection patterns interact with Arctic sea ice conditions and higher frequency atmospheric oscillations to further influence precipitation distributions by altering the location of landfalling storms and triggering additional changes in convection and atmospheric Rossby waves 35 . Such high frequency oscillations include the Pacific North America Pattern 36 , the Arctic Oscillation 37 , the Madden-Julien Oscillation 36 , and internal atmospheric variability 38 . The transition zone between the modern SW and NW dipole regions displays appreciable interannual variability, but despite the multiple factors influencing precipitation distributions, this area has remained relatively stationary at ~40°N across the western U.S. for the past 500 years 39 . the Lake Chewaucan shoreline record. Pluvial Lake Chewaucan developed in a closed lake basin in the desert region of the NW, at the northwestern edge of the Great Basin (Fig. 1a). Modern annual precipitation in the basin peaks during the winter months, varying with elevation from 240 to 900 mm/yr 40 . Today, the closed drainage system terminates into two separate sub-basins, Summer Lake and Abert Lake, that were integrated at high lake level (>1338 m asl; Fig. 1b) to form Lake Chewaucan. Our lake level reconstruction for this system combines new shell and tufa 14 C ages with those from previous work [41][42][43][44][45] to date shoreline and nearshore lake deposits. Dated samples closely constrain lake surface elevation at the time of deposition (see Methods and Supplemental Information). The record spans intermittently from 21.3 ka to present. We interpret millennial-scale intervals with no ages to represent periods of lake regression below the elevation of our samples and/or subaerial exposure and erosion of deposits of that age (Fig. 2).
Three ages indicate Abert Lake level was moderately high ~21.0 ka, the lake receded from ~20.5-18.5 ka, then rose to nearly 1300 m asl from 18.4-17.8 ka. Two ages from the Summer Lake basin indicate a higher lake level of ~1324 m asl was achieved in this basin during the same ~21.0-18.0 ka interval. During the interval 17.8-14.6 ka, no ages were obtained on shoreline materials from either Summer or Abert Lake, suggesting a lake regression in both basins. This was followed by a rapid transgression to the high shoreline elevation (1350 m asl) at ~14.5 ka, as recorded by six overlapping tufa ages in both basins. Tufas of 14.5-13.4 ka age vary in elevation, indicating lake level may have fluctuated by ~30 m (Fig. S2), or that tufas formed at a range of elevations in the deep lake. However, in both sub-basins the deposits indicate lake level was mostly higher than the spilling threshold (1338 m asl), creating an integrated lake (Figs 1b and 2). All shell ages also fall within this interval. Shells collected near the highest shoreline elevations (1345-1350 m asl) agree well with the high shoreline tufa ages (~14.5 ka). All others have a narrow range of ages (~13.7 ka), consistent with the interpretation that the deposits represent a shell bed covering the lake bottom over a range of elevations, recording deep lake conditions in the Abert Lake sub-basin (see Methods and Supplemental Information).
Two ages on tufa, and one on a shell, indicate both lakes retreated below the spilling threshold by 12.9 ka, after which the lack of ages from 12.9-11.5 ka suggests the lakes receded below our sampled elevations. A lake transgression spanning 11.5-9.5 ka, reaching 1314 m asl, is recorded by four tufa ages on an extensive boulder shoreline exposed on the east side of the Abert sub-basin. No ages were obtained between 9.5 ka and 3.9 ka, indicating that Abert Lake receded to low lake level during the middle Holocene. In the Summer Lake sub-basin, the 7.6 ka Mazama tephra is found interbedded with dune deposits derived from deflation of the playa during the mid-Holocene 46 , in agreement with low Abert lake level. After this dry period, tufa ages from the playa along the Abert Rim indicate Abert lake level was similar to that of historical variability (<1300 m asl) 47 intermittently from 3.9 ka to present (Fig. 2).

Dipole-like hydroclimate change during the HS1-B/A transition.
A striking aspect of our shoreline reconstruction is that the timing of deep lake conditions for Lake Chewaucan are out of phase with most closed lake systems in the western U.S. (Figs 1a and 3). No 14 C ages fall within the HS1 interval ( Fig. 3a) 48 , including the HS1b period when many other Great Basin closed lakes were deepest ( Fig. 3e-g), and both wetlands 5 and lake sediment (Fig. 3h) 4,49 records indicate exceptionally wet conditions in the SW region. Instead, the deepest lake conditions at Lake Chewaucan developed abruptly at the beginning of the B/A and persisted throughout the interval (Fig. 3d). This pattern is more similar to records found in the NW, where a shift from dry climate to wetter conditions occurs in many pollen-based vegetation records north of Chewaucan at ~14.5 ka 50 ( Fig. 1a; Table S1). These climate conditions are also expressed in the diatom stratigraphy of Klamath Lake, just ~100 km to the west 51 . Overall, the records in our compilation (Fig. 1a) define a map pattern of drier (wetter) conditions during HS1b (B/A) in the NW, in contrast to wetter HS1b (drier B/A) conditions in most records south of 40°N. This pattern is similar in spatial extent to the winter precipitation dipole that characterizes interannual climate variability in the modern instrumental record (Fig. 1a) 22,34 . While the overall dipole pattern during this HS1-B/A shift is well supported by paleoclimate evidence, the transition zone exhibits a minor shift in latitude compared to the modern pattern. The deepest conditions observed for pluvial Lake Surprise coincide with HS1b, resembling the lakes in the SW region more than nearby Chewaucan (Fig. 3e) 44 . At 41.5°N, Lake Surprise, along with Lake Chewaucan, is located in the modern transition zone of the dipole pattern (Fig. 1a). Unlike Lake Chewaucan to the north (Fig. 3d, this study) and lakes to its south ( Fig. 3f, g) 52,53 , Lake Surprise remained >100 m deeper than the modern level throughout HS1, the B/A and YD, consistent with a transitional position. This suggests the present transition zone was shifted to the south of Chewaucan, but remained over Lake Surprise during the Deglacial. Despite small shifts in transition latitude, paleoclimate evidence strongly supports a persistent dipolar precipitation pattern similar in spatial extent to that expressed in modern interannual precipitation variability 34 .
In addition to the HS1-B/A interval, a dipole pattern of similar extent is visible in compilations of LGM paleoclimate records 8,9,54 (21 ka), indicating winter storms are displaced to similar northerly or southerly mean tracks under full glacial, Deglacial millennial-scale (this study), and modern conditions 39 . There is also some evidence a dipole-like pattern appears coincident with earlier stadial/interstadial climate shifts. Previous records have identified abrupt, millennial-scale moisture changes associated with Dansgaard-Oeschger (D-O) cycles during MIS 2-3, with moist (dry) conditions inferred in the NW during interstadials (stadials) 51,55 , balanced by moist conditions during D-O and Heinrich stadials recorded by wetlands and lake sediments in the SW 5,56 . However, other records have noted wet D/O interstadial conditions as far south as 34 °N closer to the west coast 57 , suggesting this regional pattern may differ from a strict north-south dipole during prior millennial-scale changes.
Ice sheet and AMoC forcing of Deglacial storm track position. It has been widely suggested that the westerly storm track shifted equatorward in mean winter position and/or strengthened, bringing more frequent or wetter winter storms to the western U.S. during the LGM and Deglacial periods 8,9,[58][59][60] . However, the relative importance of the forcing mechanism(s) causing these storm track shifts are less well constrained. Early hypotheses 61 and climate model simulations 58 focused on the atmospheric effect of the Laurentide Ice Sheet (LIS). Recent model/record comparisons have confirmed that persistent high atmospheric pressure over the LIS likely affected the storm track, favoring more southerly incidence of storms into the western U.S., and creating a dry-NW/ wet-SW dipole pattern relative to modern climate during peak LGM conditions (21 ka) 8,9 . Likewise, rapid melting and partial collapse of the LIS and Cordilleran Ice Sheet during the B/A (Fig. 3b) reduced this ice sheet perturbation, creating a poleward shift in the storm track 16,18 and moisture delivery 18 , contributing to the wet-NW/ Abert Shells (Licciardi, 2001) Abert Shells (Friedel, 2001) Abert tufas (this study) Summer tufas (Friedel, 2001) Summer tufas (this study) 10000 Abert tufa (Gehr, 1980 Figure 2. Lake Chewaucan lake level hydrograph, 0-25,000 cal yr BP. Abert Lake and Summer Lake playa elevations and threshold elevation connecting the Abert and Summer Lake sub-basins shown by dashed lines. Abert Lake (blue) and Summer Lake (gold) lake levels constrained by 14   www.nature.com/scientificreports www.nature.com/scientificreports/ dry-SW mode recorded by Lake Chewaucan and other records after the HS1-B/A transition (Figs 1a and 3). The LIS also likely continued to perturb atmospheric circulation into the early Holocene (>8 ka), maintaining wetter than modern conditions in the western U.S. despite peak Holocene summer insolation 62 . Ice sheet topography clearly was an important driver of storm track changes across the Deglacial interval. However, ice sheet changes do not explain the correspondence between the dipole-like shift in moisture in the western U.S. and North Atlantic-driven millennial scale climate changes (Fig. 3), which are known to have caused similar meridional shifts in circulation patterns in both tropical 63,64 and mid-latitude regions across the globe 7 .
Previous work 7 has suggested that variation in the interhemispheric meridional temperature gradient (centered at the thermal equator, or region of maximum mean average temperature), driven largely by Northern Hemisphere warming and Arctic sea ice change, may also influence the westerly storm track, and may explain the correspondence between millennial-scale dipole switches and global circulation shifts. The wettest conditions in the SW region records largely coincided with HS1b, displaying a dry-NW/wet-SW dipole pattern (Figs 1a and 3). During this time, massive freshwater discharge into the North Atlantic culminated in a major reduction of the Atlantic Meridional Overturning Circulation (AMOC; Fig. 3c) 65 . This led to millennial-scale cooling of the Arctic (Fig. 3a) 48 and southward shifts in precipitation-bringing circulation patterns in both hemispheres 7,63,64 . The shift to a wet-NW/dry-SW pattern at ~14.5 ka coincides with abrupt strengthening of AMOC, NH warming, and a northward shift in global precipitation patterns, consistent with a northward shift in the thermal equator 7 . This close correspondence suggests precipitation changes in the western U.S. form part of a larger global pattern in circulation changes that extend beyond ice sheet influence on the Pacific storm track.
To examine the global circulation dynamics that may have contributed to the observed HS1-B/A dipole shift, we compare storm track activity and zonally averaged vertical motions calculated from two 50-year time slices of daily output from the TraCE-21k transient climate model simulation (see Methods) 19 for AMOC shutdown conditions during HS1 (17.00 ka) and the resumption during the B/A (14.36 ka). Importantly, the B/A timeslice used here (14.36 ka) simulates conditions before the major ice sheet saddle collapse 18 , prescribed at 13.87 ka in the simulation 19 , to investigate the effect of AMOC-driven change on storminess. It has been noted that TraCE-21k underestimates the increase in precipitation required to create the deep lake conditions of the Deglacial 11 . However, it does reproduce the general spatial pattern of moisture change through the Deglacial interval seen in paleoclimate records 16,18 . The coarse resolution of global atmospheric models has been shown to introduce biases into important finer scale processes such as atmospheric rivers 66 . This renders the 3.25° horizontal resolution TraCE-21k insufficient to robustly assess moisture delivery in the present study. Nonetheless, as the only publicly available model with daily temporal resolution for the HS1 and B/A time slices, TraCE-21k still offers the ability to explore planetary-scale dynamics that accompanied abrupt climate change during the deglaciation.
Comparing global zonal-average atmospheric vertical motion (a measure of the location of major meridional circulation cells) in TraCE-21k shows a southward shift during HS1 relative to the B/A (Fig. 4a,c,e). This is consistent with a southward shift in the winter thermal equator under stadial conditions and a northward shift during the B/A interstadial 7 . Winter season (December-February) storm track activity (a measure of surface atmospheric pressure variance associated with the passage of extratropical cyclones) between HS1 and the B/A in TraCE-21k does not show a wholesale equatorward shift during HS1 (Fig. 4b,d,f). Instead, storm track activity is intensified and expanded equatorward during HS1 relative to the B/A across the north Pacific with the largest changes centered about the Aleutian Low region (Fig. 4f). This intensification and southward expansion of surface pressure variance in the Pacific basin is similar to the results of recent hosing experiments comparing LGM and HS1 conditions in the region 11 . In their simulation, AMOC shutdown conditions were characterized by a southeastward shifted and deepened Aleutian Low, and an intensified and equatorward shifted Pacific jet stream. These changes drove enhanced southwesterly wind anomalies and moisture transport from subtropical latitudes into the SW region. Moisture transport occurs primarily along the equatorward side of storm track activity maxima 9,67 , so the southward expansion and enhancement in storminess during HS1 (relative to the B/A time interval) appears consistent with an increased frequency of moisture-rich atmospheric river-type storms impacting the SW region 11 . Since the SW is significantly drier overall, the precipitation contribution of these intense events makes up a much greater proportion of total seasonal precipitation 25 and drives the interannual variability in this region 68 . While simulated changes in northern Pacific storm track activity do not demonstrate a uniform equatorward shift, our results suggest an increased frequency of intense storms in the SW region as indicated by the southward-expanded storm track. This is broadly consistent with the southward shift in the intertropical convergence zone identified in hosing experiments 11,60 and is captured by TraCE-21k (Fig. 4a,c,e). The increase in southwesterly moisture transport favored by the deepened Aleutian Low 11 (Fig. 4b) and stronger subtropical jet stream 59 during HS1 may also have facilitated more frequent inland penetration of atmospheric rivers into the interior SW region through the Mojave Desert 69 .
Relative to modern conditions, enhanced storm track activity during the B/A interval is also shifted equatorward over the North Pacific and southeastward over western North America consistent with a southward shift in zonal-average vertical motions during this time (Fig. 4g,h). This indicates that although the B/A average storm track was weaker and contracted poleward relative to HS1, it was still strengthened and further equatorward than that of the Holocene interglacial conditions. This likely reflects the effect of the large, but melting (Fig. 3b) 70 North American ice sheets on the storm track 18 , which favored more frequent and/or intense storms in the NW and promoted the exceptionally wet conditions (relative to modern) recorded by the highstand of Lake Chewaucan during the B/A.

Conclusions
The record of lake level fluctuations for Lake Chewaucan, located in the northern Great Basin, indicates peak wet conditions began abruptly at the start of the B/A period (14.6-12.8 ka) and persisted throughout the interval. This follows the HS1b period (16.1-14.6 ka) when other major closed basin lake systems indicate peak www.nature.com/scientificreports www.nature.com/scientificreports/ wetness, and occurs while lakes to the south dried. Based on a compilation of paleoclimate records for the region, this forms a pattern of north-south wet/dry conditions of similar spatial extent to the dipole observed in modern interannual winter precipitation. For the western U.S., the modern regional dipole pattern in winter precipitation is most clearly related to ENSO variability in the tropical Pacific 22 while being modulated by other modes of low and high frequency coupled atmosphere-ocean variability. However, our new lake level reconstruction for Lake Chewaucan, in conjunction with other paleoclimate records, clearly demonstrates that on millennial timescales, variations in global circulation associated with AMOC strength can favor the persistence of one average dipole state for hundreds to thousands of years. During the last deglaciation, these paleoclimate records clearly record exceptionally wet conditions in the western U.S., coupled with abrupt shifts in the regions  dashed in (a,c)). Statistically significant differences between the gridpoints at the alpha = 0.01 level are noted by stippling. Difference in vertical motions (g) and storm track activity (h) between the B/A from TraCE-21k output and the modern reanalysis observations 89 are displayed as (g,h) above for zonal averaged vertical motion and storm track activity, respectively. www.nature.com/scientificreports www.nature.com/scientificreports/ of wet and dry conditions corresponding with North Atlantic-forced millennial-scale climate changes. Global climate model simulations of abrupt climate events during the Deglacial are also consistent in suggesting a combination of atmospheric forcing by large North American ice sheets and changes in zonal-mean meridional circulations that altered North Pacific storm track and moisture delivery and created the hydroclimatic dipole signature. We recommend that continuing research on the millennial-scale Deglacial dipole focus on the development of lake level (or other proxy) records along the transition region (40 °N) as well as model studies performed at horizontal and vertical resolutions sufficient to capture key processes driving the magnitude and distribution of precipitation.

Methods
Field and laboratory sampling. Sample localities and stratigraphic sections in the study area were sampled during 2013-2015. Tufas were sampled in situ from bedrock outcrops or "in place" shoreline clasts using a hammer and chisel (Fig. S1). Mollusk shells were sampled in situ from natural exposures or hand-dug trenches (Fig. S3). The stratigraphy of the sediments was then described in detail (Fig. S4). Each sample location was recorded using a high-precision differential GPS unit. In the laboratory, we sampled tufas and mollusk shells following established methods 45 , targeting the densest tufa material and abrading the outer surface of samples that were exposed to subaerial weathering in order to minimize secondary contamination since tufa deposition. Single intact shells were used for 14 C dating in order to avoid combining material of different ages. 14 C dating. Each tufa or shell sample was reacted with 100% H 3 PO 4 and vacuum-extracted and purified 69 .
AMS and δ 13 C measurements were performed by the Arizona AMS Laboratory. Raw 14 C ages are reported with 1σ uncertainty and were calibrated using the OxCal 4.2 software 71 with the IntCal13 calibration curve 72 . All ages aggregated from previous publications [41][42][43][44][45] have been recalibrated using IntCal13 for consistency. All calibrated 14 C ages are reported in calendar years (cal yr BP) or calendar kiloyears before present (ka) as the median age of the calendar range plus or minus 2σ uncertainty.
Constructing the Chewaucan lake level record. The lake level record is based on 73 14 C ages of shorezone tufa samples collected on bedrock outcrops and boulder paleo-shoreline deposits, and shells of mollusks (Vorticifex effusa, Pyrgylopsis sp, and Valvata humeralis) that lived in shallow water habitats. We have omitted one previously reported shell age 73 , due to lack of information on the sample material and inconsistency between its reported location and elevation (see Supplemental Information).
For the lake level reconstruction, the depositional context (below water) of all sampled materials indicates sample elevation is a minimum estimate for lake surface elevation. Tufa forms at multiple paleo-shoreline levels on shoreline facies gravels in both sub-basins in areas with modern or paleo-spring discharge along faults. Their appearance on coarse shoreline gravels on high paleo-shorelines and on the modern shoreline of Lake Abert indicates they were largely deposited at only a few meters depth. V. effusa are freshwater pulmonate mollusks that feed on algae on rocky substrates like the Lake Chewaucan shorelines 74 , indicating they formed shells in shallow water. The other taxa live on aquatic vegetation, also in the near-shore environment 74 . All fossil shells were collected from the Abert Lake sub-basin, and of these, the majority were deposited in a 5-10 cm thick coquina in near-shore sandy deposits lower in elevation and distal to the shoreline gravels (Fig. S3). The coquina deposit is a conspicuous marker bed within the Abert Lake basin that comprises shells that were reworked from their natural habitat on the shorelines. Shells derived from this bed have also been reworked into gravels of lower shorelines, which are deposited unconformably on the coquina (Fig. S4, section CHL13-19). In both circumstances, significant reworking of shells is indicated. Based on these depositional characteristics, we regard the shorezone tufas to represent the best estimate of past lake level, while the shells provide minimum lake level estimates. Due to semi-continuous sedimentation into the basins since Deglacial time, all lake depths measured versus the modern playas are also regarded as minima. Normal faulting within the basin cuts shorelines and lake deposits of late Pleistocene age, creating scarps of 4-5 m height 75 , suggesting some variation in sample elevation may be imparted by tectonic uplift or drop since deposition, which may explain some variation in sample elevation within deep lake phases.
The reservoir effect within the lakes has been assessed by dating modern shorezone tufas and lake DIC, yielding bomb carbon-equilibrated ages (Table S2). In addition, previous work 70 indicates the carbon isotope composition of tufa is invariant regardless of past lake volume and consistent with inorganic precipitation from lake water DIC in equilibrium with pre-industrial atmospheric CO 2 . This finding indicates the lake epilimnion was well mixed with the atmosphere at all times in the past so that tufas and shells of shallow-water mollusk taxa should yield accurate 14 C ages. Shells are composed of 100% aragonite (based on XRD analyses of select samples) 70 indicating post-depositional weathering has not affected them. Ages from tufa and shell are in good agreement during the deepest lake interval, providing additional confidence in the age dataset. Based on this evidence, we applied no reservoir correction.
Compilation of paleohydrologic records. The compilation of 23 total paleohydrologic records (Table S1; Fig. 1a) were included for comparison to Lake Chewaucan and the TraCE-21k analysis. These include pollen-based vegetation records 1,50,76,77 , shoreline-based lake level records 44,52,53,[78][79][80][81][82][83][84][85][86] , lacustrine fauna (diatoms, ostracoda) moisture records 49,51,82 , other lake sediment-based hydrologic proxies 4,87 , and desert wetlands discharge records 5,15 . This compilation was chosen based on several criteria; the first being that they were sensitive primarily to effective moisture amount. Records that are less sensitive to effective moisture amount, such as glacier records (which respond strongly to temperature), and speleothem stable isotope records (which are also sensitive to temperature, moisture source, and seasonality) were omitted. This is with the exception of McLean's www.nature.com/scientificreports www.nature.com/scientificreports/ and Moaning Cave from the central Sierra Nevada, where trace element data were coupled with oxygen isotope data to provide a more robust hydrological signal 12 .
Records were also chosen based on their age model uncertainty, and the timespan covered during the Deglacial interval. Since the focus of our analysis was on the HS1b and B/A periods, continuous records (i.e. core reconstructions) needed to cover the entirety of both intervals. Continuous records also needed sufficiently low age model uncertainty (generally ~400 calendar years, Table S1) to clearly differentiate millennial-scale changes. For example, many pollen-based vegetation records produced before electronic datasets were widely available (so only the publication graphic of the record is available) and before routine calibration of 14 C age timescales, were omitted due to inability to translate those records objectively to a timescale comparable to more recent records. As the offset between calendar and 14 C years varies considerably between >1000 and >2000 years during the Deglacial interval 72 , these records have too much uncertainty to reliably define millennial scale changes. Age model uncertainty was estimated as the average 2 sigma uncertainty for all ages ( 14 C, U-Th series, or curve-matching) falling within the B/A and HS1 intervals (18.0-12.8 ka in total). Definition of relative moisture conditions during each millennial scale interval was based on the interpretations of proxy evidence by the original authors.
Inherently discontinuous records (i.e. shoreline records, wetlands sediment records) were chosen based on the occurrence of evidence of moist conditions at some part of the Deglacial interval, and required at least three ages within the interval of interest in order to be included (Table S1). Age model uncertainty was estimated identically to that for continuous records. Relative wetness during each millennial-scale interval with ages was defined by the interpretations of the original authors. For example, the lake level reconstruction for Jakes Lake only contains evidence of significant deep lake conditions during the HS1 interval, with the authors interpreting lack of lake deposits of other ages to indicate regression of the lake to lower levels.
TraCE-21k analysis methods. We used daily output from the 3.25° horizontal resolution and 16 isobaric level TraCE-21k global climate model simulation 19,20 to calculate storm track activity and vertical motions. Fifty-year time slices were used at 17.00 ka for HS1 and 14.36 ka for the B/A. Storm track activity (Fig. 4b,d,f) was calculated as the 24-hour, first-differenced sea level pressure variance and averaged for meteorological winter (December to February) 67,88 . Model-derived vertical motions at daily time steps during the 50-year time periods were used to estimate the mean winter position of the zonal-averaged meridional circulation and serve as a proxy for the position of the intertropical convergence zone (Fig. 4a,c,e). Statistically significant differences between the two time periods was estimated by requiring each grid point to satisfy both a two-sided Student's t-test for differences in mean values and a Kolmogorov-Smirnoff test for differences in distributions. An alpha level of 0.01 was used for both calculations. The modern (1981-2010) storm track activity and circulation, displayed for reference (Fig. 4b,d), was calculated using 24-hour, first-differenced sea level pressure variance from the 2.5° horizontal resolution and 26 isobaric level NCEP-NCAR reanalysis 89 . Simulated vertical motions were also obtained from the NCEP-NCAR reanalysis. NCEP-NCAR output was re-gridded to the resolution of TraCE-21k for calculation of differences between TraCE-21k vertical motions (Fig. 4g) and storm track activity (Fig. 4h) during the B/A and the modern (1981-2010) period.