The North American hydrologic cycle through the last deglaciation

While the climate evolution of North America during the last deglaciation has received considerable attention, few detailed model–data comparisons of the deglacial hydroclimate have been conducted at the continental scale. Here we use a transient climate simulation of the last deglaciation and a synthesis of hydroclimate proxies from across the continent to broadly assess the moisture budget and its evolution, including the primary components and mechanisms for changes, and evaluate the level of agreement between model and data in order to inform areas of major discrepancy for attention in future work. At the broadest scale, the simulation indicates that North America was wetter than modern through much of the last deglaciation, with the exception of the Pacific Northwest. This was principally due to increased moisture convergence by the mean flow over the west of the continent, and by transient eddies over the east. We find that the simulation demonstrates considerable skill in comparison to proxy records in the coastal Southwest and Pacific Northwest, with limited success in other regions. The highest disagreement occurs in the North-central region, where the model and data disagree on the sign of changes. In addition, we find that the areas of agreement coincide with regions whose hydroclimate is dominated by winter precipitation. In addition, regions of disagreement are those where summer precipitation is predominant, and where the model suggests dynamical changes are important. These results illustrate the impact of the ice sheet on the atmospheric dynamics, and therefore the importance of its accurate reconstruction, as Email address: juan.lora@yale.edu (Juan M. Lora) Preprint submitted to Quaternary Science Reviews October 1, 2019 well as that of model resolution. Our analysis and synthesis provide context for future transient climate simulations, and suggest numerous areas of priority for future research.


Introduction
The transition between the peak of the last glacial state and the beginning of the Holocene-the last deglaciation, approximately between 21,000 and 10,000 years ago (21-10 ka)-was a period of large and punctuated climatic changes (Denton et al., 2010;Shakun and Carlson, 2010;Clark et al., 2012). During the deglaciation, greenhouse gas concentrations (Marcott et al., 2014) and boreal summer insolation increased (Laskar et al., 2004(Laskar et al., , 2011, while large continental ice sheets retreated, leading to ∼100 m of sea level rise (Lambeck et al., 2014). Concurrent temperature  and atmospheric circulation (Toggweiler, 2009;Denton et al., 2010;McGee et al., 2014;Chiang et al., 2014) changes led to substantial modifications of the hydroclimate of various regions (Renssen et al., 2018), diagnosable through proxy records and climate model simulations. Thus, the last deglaciation offers a window into the evolution of regional hydroclimate during the most recent episode of important global warming.
The climate evolution of North America during the last deglaciation has been the focus of a substantial volume of research, as hydroclimate reconstructions from proxy evidence are ample and the proximity of the Laurentide Ice Sheet makes the region's changes stark and compelling. This is especially true in the western part of the continent, where dramatic deglacial water balance changes are recorded by multiple archives (e.g., Russell, 1885;Gilbert, 1890;Antevs, 1948Antevs, , 1954COHMAP Members, 1988;Bartlein et al., 2011;Oster et al., 2015a;Scheff et al., 2017;Ibarra et al., 2018), particularly in regions where future projections are uncertain since they fall in the transition zone between moistening extratropics and drying subtropics (Maloney et al., 2014;Neelin et al., 2013). Yet, few detailed model-data comparisons of the deglacial hydroclimate evolution have been carried out at a continental scale, particularly with a view towards the moisture budget and the mechanisms of its change.
Detailed analyses of the modern moisture budget-wherein net precipitation (precipitation minus evaporation; P − E) can be understood in terms of atmospheric moisture convergence by various components (e.g., Seager and Henderson, 2013)-have been performed over North America, based both on reanalysis products and simulations of modern and future climates Seager et al., 2014). These analyses indicate that, with global warming, changes of moisture convergence by the mean atmospheric flow will in particular impact net precipitation over the West, while changes in moisture convergence by transient eddies will contribute strongly over the Great Plains and Northeast. Furthermore, the changes in the mean flow convergence arise, in the context of anthropogenic warming, both as a result of the warming itself but also as a consequence of changes to the circulation.
Similarly, analyses of the moisture budget at the Last Glacial Maximum (LGM; 21 ka) show that large differences in P − E from the modern over North America are dominated by differences in the circulation, particularly over the West (Lora, 2018;Morrill et al., 2018). This coincides with the important impact of the Laurentide Ice Sheet on the circulation over the eastern North Pacific, North America, and the western North Atlantic, observed in numerous climate simulations (Manabe and Broccoli, 1985;COHMAP Members, 1988; Bartlein et al., 1998;Kim et al., 2008;Pausata et al., 2011;Ullman et al., 2014;Merz et al., 2015;Löfverström et al., 2016;Wang et al., 2018). In particular, the ice sheets forced a stationary wave that shifted the midlatitude jet and storm track south over the eastern North Pacific (Cook and Held, 1988;Bartlein et al., 1998;Kageyama and Valdes, 2000;Laîné et al., 2009;Lora et al., 2017), while intensifying and zonalizing the North Atlantic jet (Li and Battisti, 2008;Rivière et al., 2010;Merz et al., 2015;Löfverström et al., 2016). These changes strongly affected wintertime precipitation, especially over the West where moisture delivery by atmospheric rivers likely intensified into California . The impact of glacial conditions on summertime hydroclimate is less apparent, though simulations suggest enhanced P − E over the Great Plains and eastward from a strong low-level jet and increased storminess along the ice sheet margin (Bromwich et al., 2005;Lora, 2018). Simultaneously, the North American monsoon appears to have weakened because of northwesterly ventilation, a change driven primarily by ice sheet albedo (Bhattacharya et al., 2017).
Only a few studies have examined the North American hydroclimate during the deglaciation (e.g., Bartlein et al., 1998;Lora et al., 2016), in no small part due to the dearth of multi-millennial transient climate simulations. One such simulation, the Transient Climate Evolution of the past 21 kyr (TraCE21k; Liu et al., 2009;He, 2011), has been used to examine the impact of the receding ice sheets on the circulation and precipitation over western North America (Wong et al., 2016;Lora et al., 2016) and on the disposition of the North Atlantic jet (Löfverström and Lora, 2017). TraCE21k was run with the National Center for Atmospheric Research Community Climate System Model version 3 (NCAR CCSM3), with an atmospheric resolution roughly equivalent to 3.75 • in the horizontal and with 26 vertical levels. The simulation was initialized from an LGM state and forced with transient forcings including orbitally driven insolation changes, atmospheric greenhouse gases, prescribed meltwater fluxes, and prescribed land ice from the ICE-5G reconstruction (Peltier, 2004) updated at irregular intervals.
In TraCE21k, the receding ice sheets lead to a gradual zonalization of the North Pacific jet (Wong et al., 2016), punctuated by an abrupt shift northward around 13.8 ka due to a dramatic loss of ice, and accompanied by intensification of the circulation during meltwater forcing events (Wong et al., 2016;Lora et al., 2016). (The abrupt ice-sheet change at 13.8 ka is associated with Bølling-Allerød warming and Meltwater Pulse 1A, but imposed late in TraCE21k relative to these events, so the timing of deglacial episodes around 14 ka in the simulation should be interpreted with care (Lora et al., 2016).) The intensification and concomitant southeastward shift of the circulation, associated with a deepening of the Aleutian Low, leads to increased southwesterly moisture delivery to the coastal Southwest (and into the continental interior) during Heinrich events (Lora et al., 2016;McGee et al., 2018), similar to during the LGM .
The hydroclimate changes recorded by these archives allow for direct comparison to model output, particularly when quantitative hydroclimate variables can be reconstructed via proxy system modeling or transfer function approaches (Huang et al., 2002;Pribyl and Shuman, 2014;Maher et al., 2014;Voelker et al., 2015;Harbert and Nixon, 2018;Dee et al., 2018). While some compilations at the regional to continental scale have focused on the LGM (Thompson et al., 1999;Bartlein et al., 1998;Prentice et al., 2000;Bartlein et al., 2011;Oster et al., 2015a;Lora et al., 2017;Scheff et al., 2017;Sánchez Goñi et al., 2017;Oster and Ibarra, in press), few compilations have focused on compiling datasets at a continental scale for the deglaciation (Oster and Kelley, 2016;Shuman and Serravezza, 2017;Harbert and Nixon, 2018). Assembly of diverse hydroclimate datasets from the last deglaciation is important for 1) assessing the impact of climate change on the continent over the deglaciation (Lora et al., 2016), 2) validating future modeling work as new transient simulations become available (Ivanovic et al., 2016), and 3) providing context for the ecosystem changes (such as megafaunal extinction, C-4 grassland expansion and fire regime changes) that occurred during this climate transition (e.g., Gill et al., 2009;Bartlein et al., 2011;Briles et al., 2012;Cotton et al., 2016;Barnosky et al., 2016), as well as for the definitive presence of humans in the Americas during the deglaciation around 14 to 16 ka (Erlandson et al., 2007;Waters et al., 2011;Jenkins et al., 2012).
In this work, after contextualizing the hydroclimate, seasonality, and their deglacial evolution (Section 2), we analyze the North American moisture budget over the deglaciation as simulated by TraCE21k (Section 3), and then we provide quantitative comparisons to many of the available high-resolution hydroclimate records available from North America (Section 4). While not exhaustive, this work provides the first such synthesis of hydroclimate modeldata agreement at this scale.

Modern Climate
The seasonal character of the modern climate of North America varies dramatically by region, which we use as a starting point for our analysis in this paper. Figure 1 shows the modern mean annual temperature (M AT ) for North America, excluding the majority of Canada-which was glaciated for much of the last deglaciation-and Alaska. The most obvious features in modern mean annual temperature are the overall latitudinal gradient and the influence of topography on the western half of the continent; southern regions see mean temperatures that reach around 20 • C, while northern and higher altitude regions have average temperatures under 10 • C, and in some locations 0 • C.
These general temperature differences are accompanied by strong differences in precipitation, both in terms of overall yearly quantities as well as in their seasonal distribution. Figure 2 shows the seasonality of precipitation, which displays strong zonal gradients across the continent. Western North America, particularly near the coast, is dominated by wintertime precipitation; in contrast, the central part of the continent, on the eastern side of the North American Cordillera and the western central plains, is dominated by summertime precipitation. In other regions, precipitation falls more evenly across seasons, with a slight preference for the warm season in the Northeast relative to the Southeast. These seasonality differences reflect differing dominant mechanisms of moisture delivery. For example, the western regions are strongly impacted by wintertime extratropical systems and related landfalling atmospheric rivers, while the central southern regions sustain the North American monsoon and the Great Plains low-level jet from the Gulf 6 Figure 1: Mean annual temperature over North America. Modern average near-surface temperatures over North America from observational data (Willmott and Matsuura, 2001), showing the overall latitudinal temperature gradient. Side panels show the evolution of seasonal average surface temperature in each region through the deglaciation, as simulated by TraCE21k. Winter (DJF), spring (MAM), summer (JJA), and fall (SON) temperatures are shown respectively by blue, green, red, and orange curves. Black curves show only long-term trends, using a Gaussian low-pass filter (σ = 50 years). Note the different temperature ranges for each panel.
of Mexico. Therefore, these regions must be treated separately when considering their paleoclimatic evolution, because proxies may record different hydroclimate signals under different conditions, and models robustly simulate some (seasonal) mechanisms but not others.
Based on these regional contrasts, we define six overall regions for our analysis and subsequent proxy groupings that share climatic characteristics today: The Pacific Northwest, which is generally cool and wet, and receives mainly wintertime precipitation; the near-coastal Southwest, which is warmer and drier, and receives wintertime precipitation, especially closer to the coast; the central North, which is cold (with large seasonal temperature swings) and relatively dry, and receives primarily summertime precipitation; the central South (often referred to as the desert Southwest, though here we exclude CA/NV based on the seasonality), which is warm and dry, and receives mostly summertime precipitation; the Northeast, which is cool and wet, and sees large seasonal temperature contrasts but less seasonal variation of precipitation; and the Southeast, which is generally warm and wet, and receives ample precipitation through most of the year. Of course, variations exist on even more local scales due, for example, to topographic complexities, but for our purposes these six regions define useful divisions with starkly contrasting climates. Figure 1 displays the regionally-averaged seasonal temperature evolution, from 20 to 10 ka, as simulated by TraCE21k for the six regions of interest. In all cases, there is a long-term warming trend, which reflects the global increase in both temperature and atmospheric CO 2 , between the LGM at the beginning of the simulation and the start of the Holocene at the end, with significant fluctuations around 14 ka associated with that simulation's forced Bølling-Allerød (Liu et al., 2009). Other periods of the deglaciation have more subtle and more regional effects; for instance, coastal western and southern winter temperatures decrease during Heinrich Stadial I (HS1; ∼18-15 ka) despite slight increases in atmospheric CO 2 .

Temperatures through the Deglaciation
Other changes are also observable in these temperature trends, including regionally-varying changes in seasonality. In the three northern regions, seasonal temperature contrasts remain roughly constant, while in the coastal Southwest and, to a lesser degree the South, the difference between summer and winter temperatures increases through the deglaciation. The Southeast sees the opposite trend, with a slight contraction of this seasonal contrast. 8 Figure 2: Seasonality of precipitation over North America. Modern precipitation seasonality over North America from observational data (GPCC Schneider et al. (2011) and PRISM PRISM Climate Group (2010)), showing distinct regions dominated by wintertime precipitation (blue), summertime precipitation (red), or something in between. Side panels show the evolution of seasonal average precipitation in each region through the deglaciation, as simulated by TraCE21k. Winter (DJF), spring (MAM), summer (JJA), and fall (SON) precipitation are shown respectively by blue, green, red, and orange curves. Black curves show only long-term trends by using a Gaussian low-pass filter (σ = 50 years). Note the different precipitation ranges in each panel.
In addition, temperature changes from LGM to Holocene are much more extreme (a ∼20 • C range) on the eastern side of the continent, as a result of the stronger continental influence on the climate and the impact of the nearby receding ice sheet. In the Pacific Northwest, the damping influence of the Pacific mutes variations (to a ∼10 • C range) despite the Cordilleran Ice Sheet extending into the region during the LGM, while in both southwestern regions the largest seasonal LGM to Holocene warming is only about 5 • C.

M AP and E through the Deglaciation
All regions of North America undergo substantial shifts in their mean annual precipitation (M AP ) and evaporation rates through the deglaciation in the TraCE21k simulation, including in some cases dramatic changes in seasonal contrasts (Fig. 2). Large changes in the seasonality of evaporation, in particular, prompt some regions to experience drastic shifts in the seasonality and total annual magnitude of P − E (Fig. 3).
In the Pacific Northwest, fall and wintertime precipitation steadily increase after the LGM, from 18 ka to 10 ka, to average values above 4 mm d −1 . Springtime precipitation stays mostly constant at around 2.5 mm d −1 , and summertime precipitation dips slightly-from around half that value to below 1 mm d −1 -around HS1, and stays there for the rest of the deglaciation (Fig. 2). Seasonal evaporation rates are quite steady, with the exception of summertime, which increases around 14 ka and then again after 12 ka, bringing the summertime P − E to below 0 mm d −1 (Fig. 3a).
In contrast, the coastal Southwest sees a near halving of wintertime precipitation between 15 ka and 14 ka (see Lora et al., 2016), which occurs in two relatively abrupt stages straddling a slight temporary increase during the Bølling-Allerød (∼15-13 ka). Nevertheless, winter precipitation always accounts for the majority of total precipitation, as in the modern (Fig. 2). Spring and fall precipitation also decrease prior to 14 ka, though the latter first increases during HS1 relative to its LGM value, and summertime precipitation is essentially zero throughout the deglaciation. Evaporation rates for all seasons remain roughly constant in this region, such that P − E changes largely reflect the evolution of precipitation (Fig. 3b), which itself reflects changes in the atmospheric moisture transport over the North Pacific (Bartlein et al., 1998;Laîné et al., 2009;Wong et al., 2016;Lora et al., 2017). At the LGM and during HS1, P − E is positive, but around 14 ka spring and fall P − E both dip below zero, putting the annual average P − E close to zero thereafter. Over the North-central region, the annual total precipitation is dominated by summertime, as in the modern (Fig. 2). This summertime precipitation increases through the deglaciation, particularly after 12 ka. Precipitation in other seasons also increases slightly-with a step increase in wintertime totals around 14 ka-but accounts for much lower amounts than the summer. However, the seasonality of P −E experiences dramatic changes, which cause seasonal contrasts in P − E to vary substantially through the deglaciation (Fig. 3c). Summertime P − E drops from being highest at the LGM to lowest (averaging 0 mm d −1 ) just after 14 ka, and then recovers to nearglacial values at 10 ka; at the same time, wintertime P − E jumps around 14 ka to dominate the budget, and slightly declines afterward. These changes strongly reflect the presence of the ice sheet and its recession through the course of the simulation.
In southern North America, precipitation changes through the deglaciation are relatively small (and annual averages are roughly constant). Small LGM differences between winter, summer, and fall precipitation amounts disappear after 18 ka, while for springtime precipitation remains significantly lower than the rest of the year despite increasing slightly through the Younger Dryas (∼13-11.5 ka) (Fig. 2). These patterns are also reflected in the P − E of the region, which does not change considerably. There is a decreasing trend of fall P − E, which at the end of the deglaciation is roughly half of its glacial value; wintertime P − E dominates throughout, and the April-August averages are consistently negative (Fig. 3d). It should be noted that substantial fractions of this region's precipitation are delivered by the North American monsoon in late summer, but the low resolution of the TraCE21k simulation prevents this phenomenon from being properly simulated (Pascale et al., 2016). Therefore, for this region the summertime moisture budget (and by extension the total budget) presented here should be treated with caution.
In the Northeast, summertime precipitation is more than three times higher than wintertime at the LGM, and this difference is reduced through the deglaciation as the former decreases and the latter increases (Fig. 2). But similarly to the other ice-adjacent region to the west, these changes are accompanied by dramatic evaporation shifts that change the seasonality of P − E (Fig. 3e): The summer contribution shifts from most positive at the LGM to negative at the end of the deglaciation, while the winter contribution does not change much in magnitude but goes from lowest to highest (most positive). In annual averages, precipitation does not change substantially (it remains around 2 mm d −1 through the deglaciation), but P − E drops by roughly two thirds, to about 0.5 mm d −1 . The various stages of the deglaciation-HS1, the Bølling-Allerød, the Younger Dryasleave clear signatures in the precipitation and P − E trends in this region, perhaps more so than any other as a result of its proximity to the North Atlantic.
Finally, in the Southeast, summertime precipitation dominates, and, though it decreases between 18 ka and 14 ka, it largely recovers afterward (Fig. 2). September-May averages increase slightly during HS1, but then abruptly fall to below-glacial values around 14 ka and remain steady. The trends in P − E largely reflect these changes, except that the most positive values occur in winter due to reduced E, followed by fall and spring, and the summer values decrease from around 1 mm d −1 to less than zero around 15 ka, and remain approximately constant afterward (Fig. 3f). (Summertime evaporation increases from around 3 mm d −1 at the LGM to more than 4 mm d −1 by 10 ka.) In annual averages, precipitation increases and then decreases slightly, but P − E drops from more than 1 mm d −1 at the LGM to nearly 0 at 10 ka.

Moisture Budget at the LGM
Having established the overall character and evolution of the seasonality of precipitation and evaporation, in the following we turn to the net moisture budget of North America in the TraCE21k simulation, beginning with a brief evaluation of its LGM result. The LGM moisture budget of North America, including its cold-and warm-season components, has been previously investigated using the nine-model Paleoclimate Modelling Intercomparison Project, phase 3 (PMIP3) ensemble (Lora, 2018); those results form a baseline for comparison here.
The LGM distribution of annual-mean P − E is shown in Fig. 4a. As expected, P − E is weakly positive over the continent, while the surrounding ocean regions are neatly divided into subtropical areas of moisture divergence and extratropical areas (the storm track regions) of convergence. P − E is higher over the northwestern and eastern parts of the continent, and has its lowest values in the Southwest (negative net P − E over Baja California, for example, is an artifact of the model's low resolution). Figure 4b shows the differences between the average P − E at the LGM and the pre-industrial (PI; here defined using years 1750-1850 CE in TraCE21k).  Most importantly, the ice-free parts of the continent see increased P −E at the LGM, despite the globally lower temperatures, in agreement with the PMIP3 ensemble. This primarily results from increased cold-season precipitation in the west, and decreased warm-season evaporation in the east (Lora, 2018). Over the ice-sheet, LGM P − E is similar or lower than in the PI despite substantially decreased evaporation as a result of lower LGM precipitation totals over the region; the exceptions are the northern margin, which receives slightly more precipitation, and the southeastern margin, where lower evaporation exceeds the lower precipitation. The total atmospheric moisture convergence/divergence reflected in the distribution of P −E can be decomposed into contributions by the mean flow and transient eddies-that is, due to seasonally time-averaged flow (for example associated with semipermanent pressure systems) and to more rapidly changing disturbances or weather systems, like extratropical cyclones (see Appendix for details). Fig. 4c shows the mean flow component, which displays strong divergence over subtropical oceans, accompanied by convergence over the Pacific Northwest and the South and weak divergence along the icesheet margin and the Eastern seaboard. These patterns, and their LGM-PI differences (Fig. 4d), largely agree with the PMIP3 models. Notably, there is more moisture convergence over western North America at the LGM compared to PI in both PMIP3 models and TraCE21k (Lora, 2018;Morrill et al., 2018). In the western Gulf of Mexico, the high moisture convergence associated with northerly summertime flow ( Fig. 4c) is too intense and shows a southward bias compared to other models, though its relative LGM-PI difference is qualitatively similar to PMIP3. Figure 4e shows the mean moisture convergence by transient eddies at the LGM. Transient eddies generally lead to positive P − E over the continent, while divergence occurs over the oceans around 30 • N. As in the PMIP3 results, these patterns reflect decreases and increases in convergence over the continental west and east, respectively, relative to the PI (Fig. 4f). Thus, transient eddies tend to oppose the effects of the mean flow, both in terms of the annual means and their LGM-PI differences.

Evolution of the Moisture Budget
As can be intuited from the deglacial climate trends described in the previous section, various forcings were implemented in TraCE21k in order to produce a realistic transient simulation (Liu et al., 2009;He, 2011) responses (for example, P − E trends; see Figs. 3 and 5) are produced both as a direct result of nearby forcings-for example through the presence of ice sheets in certain regions-as well as through long-range effects mediated, for instance, by changes in the large-scale circulation. We next investigate the responses of the annual-mean moisture budget over North America through the deglaciation as a means of quantifying the simulation's hydroclimate evolution and identifying the dominant mechanisms of change. Figure 5 shows the trends in P − E for the six regions investigated here, and the component atmospheric moisture convergence trends that produce them. In the Pacific Northwest, the annual P −E is strongly positive throughout the time span considered, remaining roughly constant from the LGM through HS1 and into the Bølling-Allerød, then increasing to a maximum at the beginning of the Holocene. However, these small changes disguise significant changes in the balance of moisture convergence by atmospheric mean and transient components: During the LGM, convergence by transient eddies contributes more than three times more to P − E than convergence by the mean circulation, but by the end of the deglaciation, moisture convergence by transient eddies is slightly lower than that by the mean flow. Around 18 ka, the convergence by transient eddies and the mean slightly decrease and increase, respectively, maintaining a roughly constant total. Then, starting around 15 ka, abrupt and gradual changes through late HS1 and the Bølling-Allerød conspire to strengthen these trends, and after 14 ka to 10 ka the two components contribute roughly equally to the total P − E.
In the coastal Southwest, the annual-mean P − E is positive during the LGM and HS1-in agreement with copious proxy evidence (see, for example, Oster et al., 2015a;Lora et al., 2016)-and then zero after 14 ka, with a rapid transition between 15 and 14 ka (Fig. 5b). Throughout the LGM and deglaciation, transient eddies converge moisture while the mean flow diverges it on average, but the relative importance of the two components changes substantially. Between 22 and 15 ka, divergence by the mean flow is on average very small, while transient eddy convergence is around 1 mm d −1 . At the Bølling-Allerød transition, mean flow divergence and transient convergence both increase, but the former does so more strongly and thus shifts the total budget to lower P − E. Thus, water delivery to the region becomes more transient while drying by the mean flow increases. After 14 ka, these trends continue but much more gradually, and the balance between transient convergence and mean flow divergence is maintained. The conundrum of lower transient eddy activity during a time of higher baroclinicity in the LGM and early deglacial has been noted previously (Kageyama and Valdes, 2000;Li and Battisti, 2008), and an explanation (suggested for the North Atlantic but potentially also applicable to the Pacific) suggests that upperlevel storms seeding the strong baroclinic zones are weaker at the LGM than in the modern (Donohoe and Battisti, 2009).
Changes in the moisture budget along the west coast of North America are indeed related to the baroclinicity in the North Pacific, and along the western edge of the continent, which is relatively high between 21-15 ka ( Fig. 6a,b), coincident with the largest ice-sheet extents. The same glacial conditions also lead to relatively lower water vapor over the eastern North Pacific, which shows more meridional structure (Fig. 7a,b) than in the modern; in particular, early-deglacial annual mean column vapor is larger or equal at extratropical compared to subtropical latitudes along some meridians of the eastern North Pacific (∼130 • W), with a distinct band of high moisture extending southwest to northeast from the subtropics around the dateline to the coast of California. These patterns are also related to the fact that maximum temperature gradients are shifted southward relative to the modern in the central North Pacific. Thus, transient eddies converge moisture in the Pacific Northwest despite a southward-displaced storm track, while the mean flow dries the coastal southwest less than in the modern as a result of the climatological distribution of water vapor. In sum, these changes are in agreement with the interpretation that atmospheric rivers more efficiently channel water vapor-and therefore imprint on the mean convergence-between the subtropical ocean and southwestern North America during the peak glacial, despite lower overall transient activity, at the expense of moisture in the Pacific Northwest . At the Bølling-Allerød transition, North Pacific baroclinicity decreases relative to earlier stages, but the meridional temperature gradients over the southwest of the continent do not (Fig. 6c). At the same time, the annual-mean water vapor structure homogenizes zonally relative to the HS1 (Fig. 7c), indicating more variable moisture delivery to the southwestern part of the continent.
In the central North region, moisture budget trends are temporally smoother than those of other regions (Fig. 5c). The mean P − E is largely constant through the deglaciation, with a slight decrease during the Younger Dryas; however, the components of moisture convergence change dramatically. Convergence by transient eddies falls from a maximum of around 0.8 mm d −1 at the LGM to slightly negative by the Younger Dryas, decreasing at every stage of the deglaciation. Meanwhile, convergence by the mean flow increases from negative values at the LGM to a maximum around 0.5 mm d −1 around 13 ka, becoming positive around 15 ka. These changes reflect a decrease in transient activity in the region through the deglaciation resulting from the receding ice sheet and concomitant loss of baroclinicity in the region (Fig. 6).
The evolution of moisture convergence and its components is least substantial over the central South region (Fig. 5d). In this region of overall low average P − E, the mean flow diverges moisture while transient eddies converge it, and this occurs throughout 22-10 ka. There is a slight intensification of the divergence by the mean flow through this time-with the most rapid changes occurring around 18 and 12 ka-which results mainly from decreasing mean flow convergence to the north and east of the region. The net effect is to roughly halve the net P − E from close to 0.5 mm d −1 at the LGM to around 0.25 mm d −1 at 11 ka. Again, however, it is worth noting here that the TraCE21k simulation cannot resolve the North American monsoon, which is of primary importance in this region, and therefore these  results should be treated with caution. Previous work has suggested that monsoon rainfall was likely lower at the LGM as a result of northerly ventilation caused by the ice sheet (Bhattacharya et al., 2017), though decreases in evaporation may have also contributed to changes in P − E. The eastern part of the continent sees the most dramatic shifts in P − E through the deglaciation in TraCE21k, which are driven primarily by major shifts in evaporation caused by the presence of the ice sheet during the LGM (Lora, 2018), and its subsequent recession. During the LGM, transient eddies strongly converge moisture along the ice-sheet margin ( Fig. 5e; Lora, 2018), while there is relatively weaker divergence by the mean flow. At HS1, the mean flow divergence drops approximately to zero and the transient convergence decreases slightly, such that the total P − E, dominated by the latter, slightly increases. The simulation's forced Bølling-Allerød leads to large fluctuations between 15 and 13 ka, but a substantial decrease in convergence by transient eddies from the end of HS1 through 10 ka leads to late deglacial P − E totals that are less than one third of their early deglacial values.
Similar relative patterns produce comparable shifts in P −E in the Southeast (Fig. 5f). Average P − E drops from around 1.2 mm d −1 at the LGM to around 0.2 mm d −1 at 10 ka. In this region, both the mean flow and transient eddies produce moisture convergence from the LGM through the Bølling-Allerød, but convergence by transients decreases during HS1 (and continues decreasing afterward), while convergence by the mean flow roughly doubles during HS1 and then returns to LGM values around 14 ka. Large fluctuations also occur between 15 and 13 ka, but they are less significant than in the Northeast. And, unlike in the other five regions where variations about the long-term trends remain mostly steady, decadal variability in the components of moisture convergence increase considerably over the Southeast at the Younger Dryas.

Components of Mean Flow Moisture Convergence
Moisture convergence by the mean flow is important in all regions of North America throughout the deglaciation, in agreement with PMIP3 LGM results (Lora, 2018;Morrill et al., 2018). Therefore it is convenient to assess changes in its own components (mass divergent flow, advection across moisture gradients, and moisture transport along surface pressure gradients (Seager and Henderson, 2013;Seager et al., 2014;Lora, 2018;Morrill et al., 2018)), as well as investigate the mechanisms (thermodynamic versus dynamical) for changes (see Appendix for details).
In the modern, the mass divergent flow-upwelling or subsidence near the surface where humidity is high, associated with low or high pressure centersgenerally dries the atmosphere over the western part of the continent and moistens it in the region east of the Rockies, while producing relatively little convergence over the eastern regions (see Seager et al., 2014;Lora, 2018). We note that the TraCE21k simulation strongly overestimates the PI moisture convergence (moistening) due to this component in the central part of the continent relative to other models and reanalyses, likely as a result of the low resolution. Nevertheless, while there are slight changes in their magnitude, in all regions these overall patterns remain through the deglaciation (Fig. 8a shows these at the LGM). The biggest change is that drying by the mass divergence in the Pacific Northwest intensifies, particularly after 14 ka, and this is accompanied by increased moistening in the central North. In other words, these drying and moistening tendencies were lower in these two regions through much of the deglaciation relative to the modern climate.
Advection across moisture gradients, involving horizontal transport by winds, tends to moisten the Pacific Northwest, while drying the central regions, the coastal Southwest, and along the northeastern coast of the continent in the modern (Fig. 8b). As with the mass divergent flow, these patterns are not substantially altered through the deglaciation-at least at the resolution of the TraCE21k simulation-with the exception of intensified moistening in the Pacific Northwest that roughly compensates the increased drying by the mass divergence.
Lastly, the surface term (Fig. 8c), which describes convergence due to near-surface moisture transport along surface pressure gradients and is therefore important near topography, is generally much smaller than the other two components, and also remains close to constant (and of the same sign) in most regions through the deglaciation. The exception, again, is in the Pacific Northwest, where moistening by this component increases substantially after 14 ka, leading to increases in moisture convergence by the mean flow generally (see Fig. 5a and Lora et al., 2016), highlighting the particularly strong interaction between the mean flow and the retreating ice sheet in that region.

Mechanisms of Mean Flow Moisture Convergence Change
The evolution of the mechanisms of moisture convergence differences, relative to the PI, is more interesting and illuminating. As a starting point, Figure 8: Components of moisture convergence by the mean flow at the LGM. Maps of annual-average LGM moisture convergence due to (a) the mass divergent flow, (b) advection across moisture gradients, and (c) moisture transport along surface pressure gradients in the TraCE21k simulation. Thick black contours denote the LGM margins of the ice sheet. Fig. 9 shows the mechanisms of mean-flow moisture convergence differences between the LGM and the PI in TraCE21k.
The thermodynamic mechanism, involving changes in atmospheric humidity independent of circulation differences between the two climate states, leads to a weakening of the hydrologic cycle over the oceans at the LGM, with increased convergence over the subtropics and increased divergence in the extratropics (Boos, 2012;Lora, 2018). In addition, there is general LGM moistening across the unglaciated parts of the continent, except over Texas and northeastern Mexico; this latter divergence increase is substantially overestimated in TraCE21k relative to PMIP3 models (Lora, 2018). Changes in humidity also lead to dramatic moistening at the western ice-sheet margin, over the Pacific Northwest, and drying at the eastern margin, plus modest drying over the central part of the ice sheet.
The dynamical mechanism, involving changes in atmospheric circulation independent of humidity differences, also produces large LGM convergence differences, in agreement with the PMIP3 models (Lora, 2018). In particular, there is increased convergence over the eastern North Pacific and into the western part of the continent, with a band of higher divergence extending from the tropics into the west coast of Mexico. To the east, circulation changes lead to a large region of increased convergence over the coast west of the Gulf of Mexico, as well as north and northeast along the ice-sheet margin, surrounding an area of modestly increased divergence. A prominent feature of this mechanism is the substantial LGM drying over the bulk of the region over the ice sheets.
And, as in the PMIP3 simulations (Lora, 2018), the nonlinear term comprising changes in both humidity and circulation produces large convergence increases over the central region of the ice sheet and divergence increases along its southern margin, which extend south into the rest of the continent. In sum, the mechanisms lead to an LGM moisture budget featuring enhanced convergence by the mean flow across the west and south, and enhanced divergence in the east and over much of the ice sheet (Fig. 4d). Figure 10 shows the temporal change in the relative influence of thermodynamic, dynamical, and nonlinear mechanisms through the deglaciation, relative to the PI, in our six regions as simulated in TraCE21k. In the Pacific Northwest, humidity differences lead to increased moisture convergence from the LGM to about 14 ka, after which time these drop precipitously when the ice sheet recedes (Lora et al., 2016). Increased divergence due to circulation changes and the nonlinear term largely compensate, but the for- c Nonlinear term 2.0 f SE Figure 10: Trends of contributions to mean-flow moisture convergence changes. Trends in decadal-mean moisture convergence changes (relative to the PI) by the mean flow due to thermodynamic (red), dynamical (blue), and nonlinear contributions (gray) from the TraCE21k simulation. Panels correspond to regions as in Fig. 3. Black curves show only long-term trends by using a Gaussian low-pass filter (σ = 50 years). mer decreases substantially in magnitude at HS1. After 14 ka, these also collapse to near-zero.
In the coastal Southwest, all terms lead to increased moisture convergence, relative to the PI, during the LGM (Fig. 10b). At the start of HS1, the dynamical mechanism increases in importance while the thermodynamic and nonlinear mechanisms decrease. At the start of the Bølling-Allerød, the dynamical mechanism decreases abruptly, and does so again after 14 ka, while the thermodynamic mechanism decreases more monotonically through the deglaciation. Importantly, the timing of these changes indicates that both the hemispheric warming due to the Bølling-Allerød and the collapse of the nearby ice sheet contribute strongly to mean-flow moisture convergence changes.
In the central North region, both the thermodynamic and dynamical mechanisms contribute to relatively more mean-flow moisture convergence at the LGM (and, indeed, throughout the deglaciation), but the nonlinear term leads to dramatically more moisture divergence in the region from the

26
LGM to 14 ka compared to the PI (Fig. 10c). As a result, the whole region experiences more drying by the mean-flow, and lower P − E, until 14 ka, as a result of the proximity of the ice-sheet margin along which the nonlinear mechanism dominates. (It is important to note that the low resolution of TraCE21k may smear out this drying effect over a larger area than would occur at higher resolution; as a comparison, the nonlinear term effects in the PMIP3 LGM model ensemble are generally more localized around the icesheet margin, except south of the James/Des Moines Lobes (Lora, 2018).) During and after the Younger Dryas, both convergence changes by the dynamical mechanism and divergence changes by the nonlinear term, relative to the PI, decrease to near zero, leaving mean-flow convergence as the result of higher convergence by the thermodynamic mechanism.
In the central South region, the thermodynamic mechanism explains almost none of the differences between mean-flow moisture convergence during the LGM and the PI, while the dynamical mechanism accounts for relatively large increases in moisture convergence that are partially offset by increased convergence from the nonlinear term (Fig. 10d). Of course, these patterns should be interpreted with caution because this region straddles areas of both convergence and divergence increases by all three mechanisms (Fig. 9). Nevertheless, the influence of the dynamical and nonlinear mechanisms both decrease quickly into and during the HS1, such that mean-flow moisture convergence is largely unchanged in the region from 17 ka to modern, with some small excursions due to circulation changes at the Bølling-Allerød and the Younger Dryas.
Changes (relative to the PI) in moisture convergence by the mean flow through the deglaciation are comparatively small in the Northeast, where the evolving contribution by transient eddies more strongly alters P − E (Fig. 5e). Yet these changes result from large, but compensating, contributions from dynamical and nonlinear mechanisms (Fig. 10e). Unlike in other regions, the compensation is almost perfect, and indeed the influence of the thermodynamic mechanism is almost negligible until late in the deglaciation (Younger Dryas onward). Increased LGM moisture convergence from the dynamical mechanism and increased LGM moisture divergence due to the nonlinear term further intensify during HS1, then decrease, rebound, and further decrease through the Bølling-Allerød. It is therefore clear that, in this simulation at least, circulation changes are primary to the deglacial evolution of the moisture budget in the northeast of North America.
Finally, circulation changes are similarly key in the Southeast (Fig. 10f).
At the LGM, dynamical and thermodynamic mechanisms contribute roughly equally to increased convergence relative to the PI, and are partially compensated by the nonlinear term. But during HS1, the thermodynamic contribution shifts from contributing slightly more convergence to slightly less, while increased convergence by the dynamical mechanism triples in magnitude, before decreasing through the Bølling-Allerød in a pattern similar to that seen in the Northeast region. From 14 ka through the end of the deglaciation, the dynamical mechanism accounts for slightly higher convergence than the PI, while the thermodynamic and nonlinear terms are near zero. This analysis provides insights into the deglacial moisture budget of North America at the largest regional scales, with the caveats previously mentioned regarding the simulation's low resolution and other biases. First, the evolution of circulation changes is important to changing relative patterns of moisture convergence by the mean flow, particularly in the southern and eastern parts of the continent. In the northern and, in particular, western regions, the thermodynamic mechanism is also very important. Finally, contributions by the nonlinear term are entirely non-negligible, particularly for regions close to the ice-sheet margins.

Comparison to Hydrological Proxy Records
In order to evaluate the hydroclimate as simulated by TraCE21k, we next compare a compilation of high-resolution paleoclimate proxy records from across North America to results from the simulation. For this comparison, we interpolate simulated variables to the location of each record to maximize self-consistency, but acknowledge that the model's low resolution precludes capturing small-scale and local features that may have impacted the data, and thus some mismatch should be expected even if the large-scale hydroclimate were accurately simulated.
We use time series of relevant variables from the proxy records, deferring to the records' interpretations by their original authors. Many of the paleoclimate proxy records analyzed here are oxygen or hydrogen isotope records recorded in organic matter, lacustrine carbonate, groundwater deposits or speleothems. As records of continental precipitation, these datasets may be controlled by many variables (M AP , M AT , P − E, humidity, moisture source variability) and the relative contribution of these variables vary both spatially and temporally (Bowen and Revenaugh, 2003;Winnick et al., 2014;Jasechko et al., 2015;Bowen et al., 2019;Kukla et al., 2019;Oster et al., 2019) as the moisture budget and land surface of the continent evolves. Additionally, proxy specific processes related to the proxy type induce isotopic fractionation related to calcite precipitation (which is temperature dependent; Kim and O'Neil, 1997), lake evaporation (which is temperature and humidity dependent; Benson et al., 2011;Ibarra and Chamberlain, 2015), or other processes. Because of the equivocal nature of proxy records (isotope time series, grain size analyses, pollen abundance records), we compare to multiple hydroclimate variables from TraCE21k (M AP , P −E, M AT ). Note that because of the limited (monthly) temporal resolution of the TraCE21k output, we are not directly testing for moisture source variability (as inferred by a recent speleothem-based synthesis; Oster et al., 2019)). Nevertheless, we attempt a quantitative comparison when possible, using single and multiple linear regressions (least squares) and skill score calculations. In addition, we qualitatively compare the simulation to lake level changes compiled by Shuman and Serravezza (2017) : Fig. 11 shows the lake status (relative to late Holocene), as well as locations for other proxy sites in our analysis, for key periods (LGM, late Heinrich Stadial 1, Bølling-Allerød and Younger Dryas) along with the annual-mean P − E difference with respect to the PI.
Unfortunately the lack of significant spatial coverage of proxy data east of the North American Cordillera (e.g., Shuman and Serravezza, 2017;Oster et al., 2015a;Wong et al., 2016;Morrill et al., 2018;Oster and Ibarra, in press) limits the possible scope of proxy-model comparisons. Additionally, while many lake records exist and are included here (e.g., Benson et al., 2011;Kirby et al., 2013Kirby et al., , 2015Kirby et al., , 2018Heusser et al., 2015;Shuman and Serravezza, 2017), very few provide high resolution quantitative hydroclimate reconstructions. Nevertheless, we use the quantitative P − E reconstructions from Pribyl and Shuman (2014) in the North-central region and apply the transfer function(s) to M AT derived by Huang et al. (2002) for application to several hydrogen isotope records in the Northeast region.

Assessing Agreement Between TraCE21k and Proxy Records
To carry out our comparison we first manipulate the proxy-derived records and the interpolated TraCE21k time series in the following ways: 1) For all proxies we consider only data covering the deglaciation (20 to 10 ka); 2) at each record's age point we take a 200 year bin (±100 years) and average the variable of interest (M AP , P − E, and/or M AT ).
We rely on two quantitative statistical comparisons to determine agreement between the TraCE21k simulation and the proxy records: First, we 29 Figure 11: Compilation of proxy locations (colored circles), lake level status (colored squares), and sites used in a macrofossil reconstruction (smaller gray circles) compared to the TraCE21k P − E anomaly (relative to pre-industrial). Proxy data and lakes are coded by wetter (blue) and drier (orange) relative to pre-industrial, where lake data are from Shuman and Serravezza (2017) and the macrofossil reconstruction is from Harbert and Nixon (2018). Lake status data are categorized for the LGM (19 to 21 ka), Heinrich Stadial 1 (14.7 to 19 ka), Bølling-Allerød (12.9 to 14.7 ka), and Younger Dryas (11.7 to 12.9 ka) based on the Greenland ice core chronology (Rasmussen et al., 2006). Proxy data sources are in Table 1. compute both simple and multiple linear regressions; in the former case we report r to indicate when variables are correlated or anti-correlated (as is the case of many isotope-based proxy records). Second, we use a skill score (S) following Hargreaves et al. (2013): where m i , n i , o i and e i are the simulated results, reference results, the "observed" values (the proxy reconstruction), and the observation uncertainties, respectively. The simulated and observed values are normalized to 10 ka or the youngest age point in the proxy record. The reference state is an assumption of no change (n = 0) following Lora (2018). The skill score S allows us to compare the performance of the simulation in terms of both the sign of the trend and magnitude of change simulated, but only where direct comparisons are available. S values of 1 indicate perfect performance and S values of 0 indicate performance no better than the reference state of no change, while negative values indicate simulation inaccuracy worse than the reference state of no change.
In the following we discuss each region in sequence. First we discuss agreement between the paleoclimate records and the TraCE21k output in the winter precipitation-dominated (Fig. 2) Pacific Northwest and Southwest regions, followed by the South region. Then we compare the TraCE21k output to a regional compilation from Harbert and Nixon (2018) that covers the Southwest as well as portions of the Pacific Northwest and South. Finally, we discuss comparisons between Trace21k and the more sparsely available paleoclimate records from the North-central, Northeast, and Southeast regions.

Pacific Northwest
Pollen records in the Pacific Northwest and the TraCE21k simulation agree on the sign of hydroclimate changes (Fig. 12). Multiple regression between M AP and M AT with pollen time series from Carp Lake, WA (Whitlock and Bartlein, 1997) and Little Lake, OR (Worona and Whitlock, 1995) demonstrate that M AP is the strongest explanatory variable over the entire deglaciation. While both data sets are limited in their temporal resolution, they both show the most significant increases in hydroclimate change after 14 ka, with both records interpreted as reflecting increasing wet conditions through the deglaciation (Carp Lake and Little Lake are relatively close in latitude, at 45.9 • N and 44.2 • N, respectively). The Carp Lake record demonstrates an expansion of forest coincident with rising temperatures and precipitation. The record from Carp Lake has a strong correlation (Table 1) with precipitation (p = 0.09) between 14.5 and 10 ka coincident with a 0.5 mm day −1 increase in precipitation simulated by TraCE21k relative to pre-14.5 ka. In contrast, TraCE21k simulates reduced precipitation at Little Lake during the Bølling-Allerød and Younger Dryas relative to the LGM and HS1. These results potentially provide direct latitudinal constraints on the northward migration of the precipitation maximum over the deglaciation (Oster et al., 2015a;Lora et al., 2016Lora et al., , 2017Lora, 2018;Wong et al., 2016). The proxy data constrain the latitudinal position of maximum precipitation over western North America as inferred by Lora et al. (2016) who demonstrated that a northward shift (from approximately 39 • N to 46 • N) in maximum precipitation was coincident with the abrupt loss of a split jet (COHMAP Members, 1988;Oster et al., 2015a;Löfverström and Lora, 2017) over the continent. However, the lack of agreement between TraCE21k and the Little Lake record further south demonstrates a potential northward bias in the location of maximum precipitation simulated by TraCE21k over the west coast.

Coastal Southwest
In the coastal Southwest region there is substantial evidence for decreasing moisture availability associated with the deglaciation. Abundant proxy evidence from both coastal and inland sites has been previously compiled for the LGM Oster et al., 2015a;Lora et al., 2017;Kirby et al., 2018) demonstrating significantly wetter conditions relative to present. Furthermore, the bulk of PMIP3 models and TraCE21k simulate large precipitation anomalies for the LGM (Oster et al., 2015a;Lora et al., 2017;Lora, 2018;Morrill et al., 2018). Post-LGM lake highstands record the changing hydroclimate through the deglaciation, with asynchronous highstands occurring southeast to northwest (older to younger) over the Great Basin during mid to late HS1 Ibarra et al., 2014;Pribyl and Shuman, 2014;Ibarra et al., 2018;McGee et al., 2018;Morrill et al., 2018;Egger et al., 2018), and only Lake Chewaucan in the northwest corner of the Great Basin reaching its highstand during the Bølling-Allerød (Egger et al., 2018;Hudson et al., 2019). This is also evident in the lake level compilation shown in Fig. 11 (Pribyl and Shuman, 2014). ) and proxy data (grey) for the Pacific Northwest region. a) Percent forest pollen record from Carp Lake, WA (Whitlock and Bartlein, 1997) (see also Whitlock et al., 2000). b) Percent P inus record from Little Lake, OR (Worona and Whitlock, 1995) (see also Long et al., 2018). Left sub-panels are time series comparisons to hydroclimate variables, right sub-panels are cross-plots of proxy data and TraCE21k interpolated output (200 year bins) colored by age with simple linear regressions plotted for visual comparison. Multiple linear regression summary statistics are reported in Table 1.
Due to the lack of high resolution local moisture balance records in the northern portion of the Southwest, we focus here on high resolution records (Moseley et al., 2016;Kirby et al., 2013Kirby et al., , 2015Kirby et al., , 2018Heusser et al., 2015) from the southern Southwest. Fig. 13 shows good agreement between TraCE21k and the paleoclimate records from Devil's Hole, NV (Moseley et al., 2016), Silver Lake, CA , and Lake Elsinore, CA (Kirby et al., 2013(Kirby et al., , 2018Heusser et al., 2015). While these proxy records do not provide quantitative hydroclimate reconstructions, we calculate significant correlations between the TraCE21k hydroclimate variables and subaqueous calcite oxygen isotopes, lake sediment grain size, and the Juniperus pollen record. Devil's Hole, a regional groundwater record in subaqueous calcite thought to reflect regional water balance and wider temperature changes (see discussion and further references in Moseley et al., 2016;Winograd et al., 1992), is significantly anti-correlated with P − E and correlated with M AT from the deglaciation (multiple r 2 =0.93). The oxygen isotope data, while a relatively low-resolution time series, demonstrate a slow decrease in moisture availability and rising temperature, in agreement with the TraCE21k simulation. Note that if the isotopic composition of groundwater was constant, this positive relationship with temperature would be the opposite sign of the expected calcite-water temperature-dependent fractionation (Kim and O'Neil, 1997), suggesting instead that the isotopic composition of the groundwater changed over the deglaciation due to changes in regional P − E or the seasonality of precipitation (Moseley et al., 2016).
The Juniperus pollen record from Lake Elsinore (Heusser et al., 2015;Kirby et al., 2018) illustrates significant drying as a slow decrease in pollen between 16 and 13 ka, primarily occurring during the Bølling-Allerød. No distinct Younger Dryas wetting is observed in the pollen record. Both decreasing M AP and increasing M AT from the TraCE21k simulation strongly match (multiple r 2 =0.90) the Juniperus pollen record.
Finally, the grain size records from lake Elsinore (Kirby et al., 2013(Kirby et al., , 2018 and Silver Lake  both record decreased percent sand through the deglaciation, showing some structure in the later part and minimal change during the LGM and HS1. Notably, the Lake Elsinore record indicates the most significant drying during the Bølling-Allerød, whereas the Silver Lake record lags by ∼2 kyr. This could be due to the relatively coastal (Lake Elsinore) versus inland (Silver Lake) location of the two lakes. Grain size is thought to primarily be influenced by runoff (i.e., P − E or frequency and intensity of M AP ; Kirby et al., 2018), so we compare the individual variables separately since M AP and P − E are inherently correlated. The Lake Elsinore record shows significant agreement (r 2 values all >0.49 and p-values < 10 −10 ) over the deglaciation, and from 14.5 to 10 ka between both variables. In contrast, due to the delayed drying in the Silver Lake record, correlation to the TraCE21k P − E and M AP is weaker (r 2 values less than 0.4).

South
Relatively few high-resolution deglacial paleoclimate records exist for the South region, though pluvial lakes did persist in inward-draining basins during glacial periods (Allen and Anderson, 2000;Allen, 2005;Ibarra et al., 2018). Fig. 14 compares the TraCE21k simulation to two speleothem-based oxygen isotope time series from Cave of the Bells, AZ (Wagner et al., 2010) and Fort Stanton, NM (Asmerom et al., 2010), originally interpreted as indicative of changes in winter moisture delivery to the region. TraCE21k simulates slow warming over the deglaciation, coincident with decreasing P − E. Both speleothem time series record some structured change in hydroclimate, though the variability and temporal trends are not coherent between the two locations before ∼15 ka (Oster and Kelley, 2016;Oster et al., 2019). The Cave of the Bells record exhibits a large ∼2‰ increase in oxygen isotopes at the onset of the Bølling-Allerød followed by a slight decrease of ∼1‰ during the Younger Dryas; in contrast, the Fort Stanton time series records variable (−6 to −10‰) oxygen isotope values from the LGM through HS1, followed by increased values during the Bølling-Allerød and a slight decrease during the Younger Dryas. It is notable that the TraCE21ka simulation does not show substantial shifts through the deglaciation, apart from a distinct increase in temperature at the onset of the Bølling-Allerød).
Both records are anti-correlated with simulated P − E and positively correlated with M AT . If the isotopic composition of drip water was constant, this positive relationship with temperature would be the opposite sign of the expected calcite-water temperature-dependent fractionation (Kim and O'Neil, 1997), suggesting instead that the isotopic composition of the drip water recorded in these speleothems changed over the deglaciation due to changes in regional P − E, similarly to the Devil's Hole record. Alternatively, changes in the seasonality of precipitation could change the isotopic composition, as recently re-interpreted for Fort Stanton (Feng et al., 2014), though the original interpretation of this record was of wintertime precipitation (Asmerom et al., 2010). Multiple regression of M AT and P −E with the  (Moseley et al., 2016). Here we plot and compare to the updated DH2-D record (see Winograd et al., 1992;Moseley et al., 2016) for details. b) Grain size record from Silver Lake, CA . c) Grain size record from Lake Elsinore, CA (Kirby et al., 2013(Kirby et al., , 2018. d) Pollen record from Lake Elsinore, CA (Heusser et al., 2015;Kirby et al., 2018). Sub-panels as in Fig. 12. Multiple and single linear regression summary statistics are reported in Table 1.  (Wagner et al., 2010). b) Speleothem oxygen isotope record from Fort Stanton, NM (Asmerom et al., 2010). Sub-panels as in Fig. 12. Multiple linear regression summary statistics are reported in Table 1.
Fort Stanton and Cave of the Bells records show relatively low correlation (r 2 values less than 0.5). These results corroborate modeling work showing that simulated precipitation changes in this summer dominated (North American Monsoon) region are not particularly well resolved (further discussed in Section 5). As such, given these comparisons and the lack of coherent changes in the two records (Fig. 14), the interpretation of these records as being dominantly controlled by various hydroclimate variables (moisture source, temperature changes, regional or local P − E changes or precipitation amounts) remains an outstanding question.

Southwest Regional Compilation
A recent compilation and synthesis of macrofossil assemblage data (derived from packrat middens) from the western United States used a quanti-tative methodology (CRACLE: Climate Reconstruction Analysis using Coexistence Likelihood Estimation; Harbert and Nixon (2015)) to estimate regional changes in hydroclimate variables since 50 ka (Harbert and Nixon, 2018). Fig. 11 shows the distribution of the macrofossil localities included in the CRACLE reconstruction, which encompasses a large region of the Southwest, as well as parts of the South and Pacific Northwest regions, as defined in our previous analyses. Harbert and Nixon (2018) aggregated the geographically diverse data into anomaly-based time series of M AT , water balance (equivalent to P − E), and M AP (as well as other vegetation-sensitive variables that we do not consider). For our analysis, we averaged over a region defined by the entire suite of macrofossil sites (Fig. 11).
Despite the potential bias due to the inclusion of geographically diverse macrofossil localities, differences in the seasonality of precipitation, and large absolute differences in modern hydroclimate, there is good agreement between the regional change in M AT simulated by TraCE21k and the plant macrofossil reconstructions (Fig. 15). Temperature is correlated over the entire deglaciation (r = 0.62, p < 10 −10 ), though skill is poor (S = −0.13) due in part to a larger increase in regional temperature simulated by TraCE21k (∼5 K) compared to that recorded by the macrofossil reconstruction (∼2.5 K). In addition, the macrofossil reconstruction records a slight cooling over HS1, leading to further disagreement. Despite this, when considering just the post-Bølling-Allerød period (14.5 to 10 ka), both correlation (r = 0.92, p < 10 −11 ) and skill score (S = 0.78) improve significantly, suggesting that the degree of warming in response the model forcing is correctly captured over that interval.
Correlation and agreement in the trend and magnitude of change between the macrofossil reconstruction and TraCE21k for P −E and M AP are mixed. Over the deglaciation, the macrofossil reconstruction records a ∼200 mm yr −1 average decrease in water balance from the LGM through the Younger Dryas and earliest Holocene, whereas the TraCE21k simulation records a decrease of ∼125 mm yr −1 (equivalent to ∼0.35 mm day −1 in Fig. 15). This results in positive correlation (r = 0.66) and positive skill scores over the entire deglacial (S = 0.40) and pre-Bølling-Allerød (S = 0.53); however, the skill score for P −E decreases significantly (S = 0.18) post-Bølling-Allerød, due to drying in the macrofossil reconstruction not occurring until ∼13 ka. Finally, the greatest discrepancy occurs between the M AP simulated by TraCE21k and the macrofossil reconstruction. Most notable is an increase suggested by the reconstruction during the latest HS1 through the Bølling-Allerød, when the TraCE21k simulation suggests a significant decrease. This result could be due to the diverse geographic region being used to make the macrofossil reconstruction, and the lack of resolution in the simulation. This is also reflected in the relatively large error envelope in the reconstruction from 16 to 12 ka shown in Fig. 15. Furthermore, as demonstrated by the paleoclimate data from the Southwest (Moseley et al., 2016;Kirby et al., 2018) and the Pacific Northwest (Worona and Whitlock, 1995;Whitlock and Bartlein, 1997), as well as the TraCE21k output, the sign of the precipitation anomaly (relative to PI) varies strongly with latitude (Lora et al., 2016;Wong et al., 2016).

North-central
The North-central region sees the most significant disagreement in hydroclimate change between available hydroclimate reconstructions and TraCE21k, with complete disagreement on the trend and sign of change over the deglacial recorded by the paleoclimate records compared to the simulated hydroclimate. Two quantitative P − E reconstructions based on facies data from lake sediment records in Wyoming (Lake of the Woods and Little Windy Hill Pond) (Pribyl and Shuman, 2014) suggest increased P − E from 18 ka onward (Fig. 16). One record, Lake of the Woods, WY, records a slight decline at approximately 12 ka, possibly associated with the Younger Dryas. While lower resolution, these records both suggest over 200 mm yr −1 increase in P − E between HS1 and the latest Younger Dryas and earliest Holocene (Pribyl and Shuman, 2014); in contrast, the TraCE21k simulation declines by approximately 70 to 150 mm yr −1 over the deglaciation.
Similarly, a lacustrine carbonate oxygen isotope record from Blue Lake, UT, representing the period during desiccation of the main (northwest) subbasin of Lake Bonneville, records a slow increase in evaporatively-enriched isotope values. This includes the highest oxygen isotope values during possible warming associated with the Bølling-Allerød. Multiple regression indicates that M AT and P − E from the TraCE21ka simulation are not well correlated with this record (r 2 = 0.13). Additionally, numerous other lake level records from the North-central region show reduced levels through the deglaciation (relative to the late Holocene), contradicting the positive P − E anomalies simulated by TraCE21k (Fig. 11; Shuman and Serravezza, 2017).  Harbert and Nixon (2018), with error evelope shown by grey dashed lines. Sub-panels as in Fig. 12. Region where Macrofossil reconstructions are derived from are shown in Fig. 11, TraCE21k was averaged over the entire region which encompasses the Southwest as well as a portion of the South and Pacific Northwest, as defined here. Single linear regression summary statistics and skill scores are reported in Table 1.  Figure 16: Comparison of TraCE21k (P − E (green) and M AT (black)) and proxy data (grey) for the Northcentral region. a) Lake carbonate oxygen isotope record from Blue Lake Marsh, UT (Benson et al., 2011). This record represents the desiccation of the main (northwest) sub-basin of Lake Bonneville. b) P − E reconstruction from Lake of the Woods, WY (Pribyl and Shuman, 2014). c) P − E reconstruction from Little Windy Hill Pond, WY (Pribyl and Shuman, 2014). Sub-panels as in Fig. 12. Multiple and single linear regression summary statistics are reported in Table 1.

Northeast
We compare the TraCE21k output to several high-resolution isotope time series from lake sediments that record hydroclimate following the retreat of the Laurentide Ice Sheet from the region. Two of the three lake records, based on compound specific hydrogen isotope records (Hou et al., 2007;Shuman et al., 2001) provide direct comparison to temperature changes via transfer functions established by Huang et al. (2002). Additionally, given the likely influence of water balance (i.e., P − E changes) on lake systems, we compare these records through multiple regression to P − E. There is an abundance of lake level records available from the Northeast (Pribyl and Shuman, 2014) following the retreat of the ice sheet during HS1; however, while most suggest drier than modern conditions, several records in close proximity to the above time series and other lake records disagree. The Trace21k P − E anomalies show significantly wetter conditions, comparable to the LGM, during HS1, followed by a drier state during the Bølling-Allerød and Younger Dryas (Fig. 11).
Correlations with TraCE21k output are mixed for this region. Skill score calculations indicate favorable comparison with Crooked Pond, MA temperature changes (Shuman et al., 2001;Huang et al., 2002) (S = 0.25), whereas Blood Pond, MA (Hou et al., 2007) does not show skill (S = −0.65), though this is improved when only considering the post-14.5 ka period (S = 0.22). Similarly, multiple regression with the hydrogen isotope records and the carbonate oxygen isotope record from Fayetteville Green Lake, NY (Kirby et al., 2002) illustrate significant correlation with both P −E and M AT (r 2 = 0.60). In contrast, multiple regression against the hydrogen isotope record from Crooked Pond (r 2 = 0.60) suggests that M AT is the primary controlling variable, as expected given the skill score. Given that most records from this region are influenced by the presence and subsequent retreat of the ice sheet at 18 to 16 ka, future work to develop records from the entire deglaciation from the southern portion of this region will be important.

Southeast
There is limited high resolution paleoclimate data available for the deglaciation from the North American Southeast. We limit our analysis to two speleothem oxygen isotope records, from Cave Without a Name, TX (Feng et al., 2014), on the west edge of the region, and Abaco Cave, Bahamas (Arienzo et al., 2015(Arienzo et al., , 2017 (Hou et al., 2007). b) Palmitic acid hydrogen hydrogen isotope reconstruction from Crooked Pond, MA (Shuman et al., 2001;Huang et al., 2002). c) Lake carbonate oxygen isotope record from Fayetteville Green Lake, NY (Kirby et al., 2002). Sub-panels as in Fig. 12. Multiple linear regression summary statistics are reported in Table 1. level estimates shown in Fig. 11 (Pribyl and Shuman, 2014) indicate drier conditions relative to PI throughout the deglaciation.
The record from Cave Without a Name was interpreted as reflecting moisture source changes (resulting in changes in oxygen isotopes), rather than M AP , P − E, or M AT (Feng et al., 2014). In contrast,the records from Abaco Cave, Bahamas, were interpreted primarily as a temperature signal (Arienzo et al., 2015). We therefore compare both records with M AP and M AT , acknowledging the differing interpretations of the two records; indeed, correlations with TraCE21k are relatively weak. The oxygen isotope record from Cave Without a Name, which suggests a wetter Bølling-Allerød, is weakly positively correlated with M AP and negatively correlated with M AT . Multiple regression indicates that both variables are significant (r 2 = 0.52; Table 1). The Abaco Cave record is weakly negatively correlated with precipitation and weakly positively correlated with temperature-contradicting the original interpretation. Unfortunately, the lack of high resolution data from this region hinders further analysis. Thus, future work using lake sediments (for example, locations compiled in Pribyl and Shuman, 2014) will greatly improve the fidelity of similar comparisons.

Discussion
A clear outcome of the comparison between various proxy reconstructions of the deglacial hydroclimate over North America and its simulation in TraCE21k is that model-data agreement tends to be strongest in the western regions of the continent, and weakest toward the east and north (though a lack of quantitative reconstructions in the east may contribute to the apparent discrepancy). As has also been shown previously (Lora et al., 2016), the (low-resolution) simulation produces a climatic evolution through the last deglaciation that agrees with reconstructions of drying and wetting in the coastal Southwest and Pacific Northwest, respectively, and the modeled broad regional temperature trends also agree with the data. On the other hand, in the Northern regions near the ice sheet, the model appears to simulate climates that trend in the opposite direction to reconstructions through the deglaciation. Likewise, the simulation does not coincide well with records from New Mexico, Texas, or the Bahamas, though as discussed above this could be because oxygen isotope time series reflect factors other than those analyzed here (M AP , P − E and M AT ). ) and proxy data (grey) for the Southeast region. a) Speleothem oxygen isotope record from Cave Without a Name, TX (Feng et al., 2014). b) Speleothem oxygen isotope record from Abaco Cave, Bahamas (Arienzo et al., 2015(Arienzo et al., , 2017. Here the most continuous record over the deglacial derived from sample AB-DC-09 is plotted and analyzed. Sub-panels as in Fig. 12. Multiple linear regression summary statistics are reported in Table 1. Interestingly, these patterns of agreement are related to those of the seasonality of precipitation (Fig. 2) and, though less definitively, those of the mechanisms of mean moisture convergence changes. In other words, regions dominated by winter precipitation correspond to good model-data agreement, whereas regions dominated by summer precipitation correspond to areas where the model and data diverge. Recent climate models (CMIP5) better simulate wintertime moisture transport and precipitation over North America relative to observations (Sheffield et al., 2013), and our results therefore extend that finding to the deglacial climate simulated by TraCE21k. In addition, winter precipitation-dominated regions-the western portion of the continent-are also regions where the thermodynamic mechanism is especially important for changes in mean moisture convergence relative to the PI. In contrast, the eastern regions, where model-data agreement is weaker, are those where changes of the mean moisture convergence are more strongly affected by the dynamical mechanism. These patterns indicate the importance of accurately simulating the changing atmospheric circulation for understanding regional hydroclimate changes, as well as of differentiating between regionally-relevant hydroclimate processes and seasonality in interpreting climate reconstructions.
An additional inference from our analysis is the possible corroboration of substantially drier LGM conditions in the eastern part of North America indicated by pollen reconstructions , which are supported by drier conditions from other lake status reconstructions (Shuman and Serravezza, 2017) in the region through the deglaciation (Fig. 11), but in contrast to recent re-evaluation of pollen data suggesting very cold but mesic conditions (Izumi and Bartlein, 2016). Precipitation changes reconstructed from pollen lead to poor agreement with many PMIP3 models, which tend to simulate little or no precipitation declines at the LGM (Lora, 2018). These discrepancies have been suggested to result, at least in part, from inaccurate reconstructions that did not account for productivity decreases from lower atmospheric CO 2 concentrations (Cowling and Sykes, 1999;Prentice and Harrison, 2009;Prentice et al., 2011;Izumi and Bartlein, 2016;Scheff et al., 2017), but our results also seem to implicate inadequate simulation of changes in the atmospheric circulation. In TraCE21k in particular, but to a lesser extent likely also in the PMIP3 models, such deficiencies may result from uncertainties in the ice sheet reconstructions.
Many model-data inconsistencies are likely due to the size and retreat timing of the Laurentide and Cordilleran ice sheets included in the TraCE21k simulation, which were based on the ICE-5G reconstruction (Peltier, 2004). In comparison to the ice sheet elevations and extent to be used in PMIP4 (ICE-6G C and GLAC-1D; Tarasov et al., 2012;Peltier et al., 2015;Ivanovic et al., 2016), the TraCE21k ice sheet (shown as ice thickness and extent in Figs. 6 and 11) has several key differences. First, the timing of the large retreat implemented in TraCE21k at 13.8 ka likely occurs too late by several hundred years, and is therefore inaccurate and incongruent with the timing and implementation of meltwater pulses (Deschamps et al., 2012;Menounos et al., 2017;Ivanovic et al., 2018). Additionally, as shown in Fig. 11 there appears to be a discrepancy between the ice sheet extent and the presence of lake-based hydroclimate reconstructions (Shuman and Serravezza, 2017), possibly due to model resolution, during the deglaciation in the North-central and Northeast regions. Furthermore, the ICE-5G reconstruction used in TraCE21k is much thicker over central Canada and thinner over eastern Canada from the LGM through the Bølling-Allerød compared to the ICE-6G C reconstruction (see the comparison in Fig. S1 of Lora et al., 2016). Given that the height of the ice sheet crucially alters the atmospheric circulation, it is unsurprising that the North-central and Northeast regions along the margin of the ice sheet show a lack of model-data agreement. In addition, as noted earlier, the low resolution of the TraCE21k simulation results in the effects of the ice-sheet margin (for example, the non-linear term; see Section 3) extending over much larger areas than would otherwise be the case.
Separately, the largest change in the hydroclimate of western North America occurred during ice sheet retreat in the latest HS1, when northern Great Basin lake levels reached their maxima (Fig. 11;Oviatt et al., 1990;Adams et al., 2008;Munroe and Laabs, 2013;Reheis et al., 2014;Ibarra et al., 2014;Oviatt, 2015;Shuman and Serravezza, 2017;Egger et al., 2018;Ibarra et al., 2018;McGee et al., 2018;Hudson et al., 2019), and precipitation increased in the northwest, as recorded by the expansion of forest ecosystems during the Bølling-Allerød and through the Younger Dryas at Carp Lake, WA (Whitlock and Bartlein, 1997). In TraCE21k, the transition between increased and decreased LGM precipitation relative to modern occurs too far north (Oster et al., 2015a;Lora et al., 2017;Hudson et al., 2019). (Indeed, a similar problem occurs in the PMIP3 LGM simulation with CCSM4-a model of the same lineage as that used for TraCE21k, CCSM3-causing it to have the lowest skill relative to precipitation proxies in the west .) This northward bias of the transition zone is maintained as the deglacia-tion proceeds, as suggested by lake highstands at ∼16 ka at Lake Surprise Egger et al., 2018), 15.3 ka at Lake Lahontan (Adams et al., 2008;Adams and Rhodes, 2019), and at ∼14.5 ka at Lake Chewaucan (Egger et al., 2018;Hudson et al., 2019), along with the Pinus pollen record from Little Lake, OR (Worona and Whitlock, 1995). All of these indicate that wet conditions persisted through the Bølling-Allerød and into the Younger Dryas; in contrast, as we have shown, TraCE21k simulates earlier drying in comparison to the Little Lake, OR record (Fig. 12). Unfortunately, while several well-dated speleothem records straddle this transition zone (Oster and Kelley, 2016;Oster et al., 2019), at present the datasets are incomplete over the deglaciation (Vacco et al., 2005;Oster et al., 2009Oster et al., , 2015b and the oxygen isotope variations are difficult to interpret due to noise and minimal variations. Still, the carbon isotope and trace metal record from McLean's Cave in the Sierras (Oster et al., 2015b) and a composite isotope record from two inland Nevada sites (Lachniet et al., 2014) appear to correlate well with the lake records from the northern Great Basin and southern California (Oviatt et al., 1990;Benson et al., 1995;Adams et al., 2008;Reheis et al., 2014;Ibarra et al., 2014;Kirby et al., 2018;Feakins et al., 2019). At present, few proxy records provide quantitative constraints on the local moisture balance recorded on the continent. Further work developing proxy system modeling approaches, forward modeling or calibration of empirical transfer functions (e.g., Huang et al., 2002;Thompson et al., 2008;Bartlein et al., 2011;Maher et al., 2014;Pribyl and Shuman, 2014;Voelker et al., 2015;Ibarra et al., 2018;Dee et al., 2018;Harbert and Nixon, 2018), as well as the incorporation of improved boundary condition constraints-for instance, ice sheet thickness and timing, with uncertainties (Ivanovic et al., 2016), and inclusion of lakes boundary conditions (Pound et al., 2014)-will allow for more direct comparisons between geologic records and climate model simulations. For example, while the lake level status database compiled by Shuman and Serravezza (2017) is an excellent qualitative metric of hydroclimate change (catchment scale P − E changes), further work (e.g., Pribyl and Shuman, 2014) to quantify these changes robustly, across the entirely of the database, will allow for greater comparison to future climate model simulations. Additionally, recent work has highlighted the importance of revamping our interpretation of isotope records from leaf wax and speleothem archives (e.g., Bhattacharya et al., 2018;Oster et al., 2019;Feakins et al., 2019), as well as the development of complementary trace element and novel isotope system (such as calcium and uranium-series) datasets from speleothems and soil carbonates that relate to local moisture balance (e.g., Oster et al., 2012;Maher et al., 2014;Oster et al., 2015b;Owen et al., 2016;Steponaitis et al., 2015;Owen et al., 2016). Finally, greater spatial coverage and corroboration between proxy types in various regions is needed. This is particularly true in areas of sparse coverage (the Northeast, Southeast and North-central regions) and regions of model-data mismatch (Pacific Northwest and Southwest transition zone) important for further assessing model performance and skill (e.g., this study; Oster et al., 2015a;Lora et al., 2017;Lora, 2018;Oster and Ibarra, in press).

Conclusions
In this work we have analyzed the North American moisture budget over the last deglaciation using the TraCE21k simulation, and compared the results to a synthesis of hydroclimate records in six major regions of the continent defined based on modern seasonality. This model-data comparison yields the following insights: • At large scales, the LGM hydroclimate of North America simulated by TraCE21k agrees with that from PMIP3 models (Lora, 2018); relative to the PI, changes over ice-free regions consist of increased moisture convergence from the mean flow over the west and from transient eddies over the east, and increased moisture divergence of smaller magnitude from transient eddies over the west and from the mean flow over the east.
• Moisture convergence by transient eddies decreases through the deglaciation in most regions of North America, except in the coastal Southwest.
• Moisture convergence by the mean flow tends to increase through the deglaciation in northern regions and decrease in southern regions; in the south and eastern areas of the continent, these changes are associated with changes in the circulation, while in the west substantial thermodynamic changes are also important.
• The best agreement between TraCE21k and the hydroclimate records occurs in the western United States. While many of these records are not quantitative hydroclimate reconstructions, their qualitative spatiotemporal agreement with the TraCE21k simulation demonstrates the utility of such comparisons.
• These regions coincide with areas dominated by winter precipitation, and for which the thermodynamic mechanism is an important contributor to changes in the moisture convergence by the mean flow.
• By contrast, regions dominated by summer precipitation, and where the dynamical mechanism is strongest-in particular the North-central and eastern regions-yield the poorest model-data agreement.
• The magnitude of warming simulated by TraCE21k over the latter portion of the deglaciation is in agreement with the large regional southwest macrofossil reconstruction of Harbert and Nixon (2018), but inconsistent with the entire deglaciation possibly due to biases caused by spatial aggregation of macrofossil sites.
• We suggest the thickness and timing of the ice sheet retreat in TraCE21k strongly influences the lack of agreement with the hydroclimate records in North-central and Northeast regions.
In addition to providing the first detailed model-data analysis of the North American hydroclimate over the last deglaciation, our intention is that the synthesis presented here will provide context for upcoming higherresolution, transient simulations of the deglaciation through PMIP4 (Ivanovic et al., 2016). Importantly, future work in developing quantitative reconstructions and increasing the density of records from under-explored regions will be crucial for accessing the spatial and temporal evolution of hydroclimate over the deglaciation. Analysis of those results with respect to the continental moisture budget will enable a better understanding of the mechanisms and processes controlling hydroclimate changes during this last major episode of global warming.
In support of these efforts, attention to several developments in future simulations will be especially helpful: First and foremost, higher resolution modeling will enable more robust model-data comparisons that better capture the regional subtleties recorded by proxies; second, a quantification of the uncertainties, and associated dynamical effects, of ice sheet reconstructions will have a major impact on our understanding of regional changes particularly on the north and eastern parts of the continent; third, isotopeenabled modeling (especially in conjunction with proxy system forward modeling) will allow for more direct comparisons to high-resolution records, like speleothems, that have heretofore proven tricky to interpret quantitatively; and, lastly, the archiving and public release of high temporal-resolution (daily or better) simulation results will enable more detailed investigations of the evolution of and impacts of systems, like atmospheric rivers (Zhu and Newell, 1998;Newman et al., 2012), that are difficult to definitively categorize as part of the mean flow or transient systems. ∇ · ps 0 uqdp = ∇ · ps 0 u qdp − ∇ · ps 0 u q dp, where the first term on the right hand side represents convergence by the mean flow and the second term represents convergence by transient eddies.
The latter term can be computed as a residual, but therefore cannot be further decomposed. These moisture budget terms, and their evolution through the deglaciation in TraCE21k, are discussed in Section 3. In contrast to the transient eddy term, the mean flow convergence term (first term on the right hand side of Equation 3) can be broken into contributions that depend on mass divergence and horizontal advection, plus a surface term (e.g., Seager and Henderson, 2013): ∇ · ps 0 u qdp = ps 0 (q∇ · u + u · ∇q)dp + q s u s · ∇p s .
These terms are discussed in Section 3.3. Lastly, changes in the hydrologic cycle between two climate states can be examined in terms of changes to each of the components of the moisture budget discussed above. Again considering the convergence of moisture by the mean flow, ∆ ps 0 ∇ · (u q)dp = ps 0 ∇ · (u r ∆q)dp + ps 0 ∇ · (∆u q r )dp, where ∆ represent differences between the two climates, and the subscript r denotes the value for a reference climate. (Higher order terms in ∆ are ignored; Lora, 2018). Here, the first term on the right hand side represents the thermodynamic mechanism related to changes in the moisture field. The second term on the right hand side represents the dynamical mechanism related to changes in the atmospheric circulation. This decomposition allows for a diagnosis of the mechanisms of changes, which is presented for the TraCE21k simulation in Section 3.4, with the reference climate being the PI.