Patterns and drivers of Holocene moisture variability in mid-latitude eastern North America

Domain: earth_science
Published: 2025-04-15
Editorial outcome: major_revision

Abstract
Proxy data for eastern North American hydroclimate indicate strong and persistent multi-millennial droughts during the Holocene, but climate model simulations often fail to reproduce the proxy-inferred droughts. Diagnosing the data–model mismatch can offer valuable insights about the drivers of hydrological variability and different regional sensitivities to hydroclimate forcing. Here we present a proxy–modeling synthesis for Holocene climates in the eastern North American mid-latitudes, including machine-learning-based water balance reconstructions and high-resolution climate simulations. These data-model results resolve prior-generation inconsistencies, show consistent spatiotemporal patterns of Holocene hydroclimate change, and enable assessment of the driving mechanisms. This agreement suggests that the secular summer insolation trend, combined with the Laurentide Ice Sheet deglaciation and its effect on atmospheric circulation, together explain the extent and duration of drier-than-present climates. In addition, our high-resolution proxy data and transient simulations reveal clear multi-centennial climate variability. In our simulations, temperature-driven increases in evapotranspiration exceed regional precipitation gains, drying much of the region during the mid Holocene. This suggests that the mid-Holocene multi-millennial drought was driven by similar processes compared to the drying trajectory projected for mid-latitude North America over this century, which is also primarily driven by warming.

Introduction
Changes in water availability have a major impact on both regional-scale ecosystem processes1,2 and societies3. In the North American mid-continent, climate model projections for the 21st century indicate a trend towards increased atmospheric and hydrological aridity with rising greenhouse gas emissions4,5. Regional water scarcity combined with diminishing groundwater reserves mean that both wildland and agrarian ecosystems are sensitive to changes in water availability6,7. Because the recent decades of anthropogenic warming have already pushed the global climate system beyond the limits of historical record5, paleoclimate archives can offer valuable insights. America is unusually wet relative to the Holocene, with most of the Holocene characterized by drier-than-present conditions for most of the United States (except southwest) and southern Canada9. Current climate models poorly simulate these Holocene moisture variations, e.g. showing greater-than-present precipitation in mid-latitude North America during the early-to-mid Holocene, a stark contrast to the peak aridity suggested by the proxy data9,10 and underpredicting the magnitude of mid-Holocene lake-level drawdown in the mid-continent11. Moreover, many geomorphic and ecological systems show high climatic sensitivity and non-linear threshold responses to the modest forcings over the Holocene, with wide variations among sites in the timing and pace of drying, suggesting the possibility of unpredictable and potentially abrupt local responses to regional aridification8.

Late-Quaternary hydroclimatic variations in eastern North America have been inferred from many sources8,9,12,13, including fossil pollen, eoliandeposits (dune mobilization, loess deposition), carbon isotopes (as a measure of water balance of C4 vs. C3 grasses6), physical sedimentology14,11,12,15, oxygen isotope and calcite-aragonite ratios6, and diatom indicators of salinity16. Among these proxies, fossil pollen records have been a backbone of Holocene paleohydrological reconstructions in central and eastern North America17–20, because of the widespread availability of lake archives, a demonstrated sensitivity of vegetation composition (and consequently, pollen assemblages) to variations in water availability and temperature, and orthogonal regional gradients in temperature and moisture. The pollen-based quantitative paleoclimatic reconstructions typically used methods such as the modern analog technique (MAT12,21), response surfaces (a smoothed form of MAT18), or other classical transfer functions such as weighted averaging-partial least squares13,22. Other studies have employed semi-quantitative paleoclimate measures such as z-scores of proxy variability9, or quantitative reconstructions of variables indirectly related to moisture balance, including salinity23 or changes in vegetation cover24.

Key features of prior paleoclimatic reconstructions include long-term increases in moisture availability that produced wetter conditions today across much of mid-latitude North America compared to the early- and mid-Holocene, as well as abrupt millennial and centennial departures from these trends. Regionally, while the mid-continent was likely wet during the early Holocene, the eastern coastal areas were substantially drier than today25. The pattern reversed around ca. 8 ka, when large portions of the mid-continent dried8,15,19,23,26, in contrast to the eastern areas that experienced ~20% increases in effective precipitation1,27,28. A leading hypothesis suggests that the hydroclimate reorganization reflects the rapidly diminishing influence of the Laurentide Ice Sheet towards the end of the early Holocene12,15. After around 8 ka, a long-term increase in effective moisture, amplified by a rapid moistening at around 5.5 ka, affected most areas, possibly in response to long-term insolation trends and superimposed millennial variability12,29. These millennial and multi-centennial hydroclimate influenced eolian activity in the Great Plains30,31, altered water levels in lakes and wetlands14,32, and transformed forest composition1.

In this study, we present an analysis and synthesis of hydroclimate trends in the eastern North America, drawing on recent advances in the available proxy data, quantitative climate reconstruction methodology, and paleoclimate simulations. First, machine-learning based approaches have been adopted in proxy-based paleoclimatic reconstructions. These include ensemble models of regression trees (e.g., random forests, boosted regression trees) and artificial neural networks33. Here we focus on regression tree ensembles, which were first applied in biogeography to contemporary ecosystems34,35, and are now increasingly applied to paleoclimatic reconstruction (reviewed in refs. 33,36,37). A key strength of regression-tree based models is their ability to detect comparatively weak climatic signals in datasets, such as those associated with secondary (or tertiary, etc.) environmental drivers34. Regression trees achieve this by first screening the training data for indicator taxa showing strong responses to each climatic variable, and then selectively focusing on these taxa in prediction36,37. This characteristic of regression trees is useful in moisture reconstruction because North American pollen datasets primarily respond to temperature (with summer temperature carrying more predictive power than winter temperature), but with a clear secondary signal of moisture, with certain pollen taxa responding strongly to moisture-related variables17,20,37,38. Second, a suite of high-resolution and well-dated pollen records from the last two decades (Supplementary Table 1) now permit assessment of centennial-scale climate variability at local to regional scales. Third, an increasing number of transient Earth system model simulations for the Holocene is now available39, with higher spatial model resolution and various forcings, including solar and volcanic forcing40 and dynamic vegetation changes41 in addition to changes in orbital and greenhouse gas forcing.

Here we apply the boosted regression tree (BRT)-based proxy–climate calibration models37 to prepare synthesis reconstructions of annual water balance (calculated as annual precipitation minus evapotranspiration42) and July mean temperature (Tjul). These reconstructions are then integrated with an ensemble of four climate model simulations for the Holocene. Our paleoclimate reconstructions are based on fossil pollen sequences from three regions spanning the eastern North American mid-latitudes: Midwest (MW), Great Lakes (GL) and Northeast (NE), which collectively extend from the prairie–forest ecotone to the eastern seaboard (Fig. 1). While the general timing of multi-millennial drought is well established by available multi-proxy evidence, we explore open questions that especially require quantitative moisture reconstructions. First, we assess the magnitude and spatiotemporal progress of Holocene droughts, facilitating a direct comparison with trends from the Holocene simulations. Second, using a high-resolution subset of the data, we examine the shorter, centennial-scale moisture variations, including their prevalence, geographic distribution, periodicities, and magnitudes. Finally, we test the hypothesis that the spatiotemporal patterns of Holocene aridification in eastern North America are linked to atmospheric circulation changes driven by the retreat of the Laurentide Ice Sheet (LIS).

Results
Reconstructed millennial trends in moisture and temperature
Water balance reconstructions using the BRT calibration models from all three regions (Fig. 2a) show multi-millennial, early-to-mid Holocene relative drought with virtually the entire error margins of the synthesis curves under the preindustrial value. There is a notable longitudinal gradient in the timing of peak drought, which commences first in the NE at interglacial onset, with a sharp Holocene water balance minimum of ca. −180 mm reached at around 11 ka. Conversely, in the MW, peak aridity is reached later in mid Holocene (ca. 7 ka), with a reduced amplitude of ca. −110 mm. In the GL, the timing of the the water balance minimum at 10–6 ka falls between those in NE and MW and has a similar amplitude with MW at −100 mm. The fossil pollen samples generally have good modern analogs, with the compositional distance (squared-chord distance) to best-matching modern pollen sample largely around 0.1–0.15, and only 288 fossil samples out of 5257 (5.5%) exceeding a threshold distance of 0.25 recommended for a good analog38 (Supplementary Figs. 1–10). In the MW and NE, the analog quality remains good through the Holocene suggesting a lack of biases in the reconstructions due to poor modern analogs. However, in the GL analog quality is reduced through 12–10 ka (Supplementary Fig. 5). In consequence, in the GL, the timing of the onset of negative water balance anomalies is left uncertain, as the wide error band allows for trajectories ranging from an NE-like early onset to one identical to MW.

The Tjul reconstructions (Fig. 2b) show a Holocene temperature maximum (HTM) with a broadly uniform timing in all regions, with the millennial-bandwidth smoothers indicating a plateau at ca. 7 ka (Fig. 2b, lower panels), however with the error bands in the NE allowing for the possibility of an earlier temperature maximum at 11–8 ka. The amplitude of peak warmth decreases towards the Atlantic (ca. 0.9–1.0 °C in the MW, 0.8–0.9 °C in the GL, and 0.6–0.7 °C in the NE).

Some aspects of the proxy reconstructions largely confirm prior findings about the broadscale millennial moisture and temperature trends. In the well-studied MW region, the temporal evolution of the multi-millennial Holocene drought matches that identified in previous studies based on shifts in pollen percentages of forbs and moisture-sensitive tree types, corresponding pollen-based precipitation reconstructions, and lake-level reconstructions, with drying commencing at ca. 10 ka and with peak drought conditions reached by ca. 7 ka15,17,19,20 (Fig. 2a). In contrast, NE shows an early-Holocene drought followed by an increasing water balance and multi-century variability since 11 ka (Fig. 2a), which aligns with recent lake-level and pollen data15,32 and derived quantitative moisture reconstructions from this region1,20,29. However, fewer prior data exist for the intervening GL sector, where our water balance reconstruction shows that the timing of the lowest water balance in the early Holocene occurs between the timing in MW and NE (Fig. 2a).

The pacing of reconstructed water balance vs. Tjul changes is roughly synchronous in MW, while in GL and NE, the earlier onset of drought results in the water balance minimum preceding the temperature maximum (Fig. 2a, b), supporting the notion that independent signals of Tjul and water balance can be extracted from eastern North American fossil pollen data20,37,38. Our Tjul reconstructions could be slightly biased by postglacial vertical land movements (forebulge collapse), with subsidence in the southern part of the GL region estimated at close to 100 m since 10 ka43. At present, the rate of subsidence is highest in the MW and GL, while for NE the rate is smaller or transitions to glacio-isostatic uplift in the northern part of this cluster44. The climatic impact of the vertical land movements could be regionally up to ~ 0.6 °C based on a lapse rate of 6.4 °C/km, which falls within the uncertainty intervals of our reconstructions. However, the expected effects of the subsidence are to negatively bias the HTM amplitudes, due to higher mid-Holocene site elevation in forebulge collapse areas, and to delay the reconstructed timing of the HTM because the negative temperature bias increases towards early Holocene and thus delays the arrival to the mid-Holocene temperature maximum. Thus the bias would be strongest in the mid-continent and so would tend to mute the reconstructed gradient in HTM amplitude, which in our reconstructions increases from the Atlantic towards the mid-continent. Hence, correcting for the glacio-isostatic effect would reinforce the observed pattern.

A principal components analysis (PCA) of the water balance reconstructions covering all three site clusters reveal two dominant modes of variability (Fig. 3a), with the first component (PC1) explaining 37.0% and the second component (PC2) 22.5% of the total variance, and the remaining components not exceeding 6.1%. As shown by the positive PC1 site loadings at nearly all sites (50 out of 53; Fig. 3b), PC1 captures the significant shared water balance trend of all sites with an early-to-mid Holocene drought transitioning into a late-Holocene wetting pattern. This trend is evident in PC1 sample scores, which show a decline until 8 ka, followed by a rising trend towards the present (Fig. 3c). Based on the PC2 site loadings (Fig. 3d), showing clusters of positive loadings in MW and negative loadings in NE, PC2 represents the temporal variation in longitudinal water balance contrast. The PC2 sample scores (Fig. 3e) reach a maximum in the Early Holocene (ca. > 8.5 ka), representing the period of maximum contrast in water balance anomalies due to peak multi-millennial drought in NE before the early wet conditions ended in MW (Fig. 2a). However, the northernmost sites in NE showing positive PC2 loadings (Fig. 3d) do not reflect this trend, instead covarying with the mid-continental sites.

Comparison with climate model simulations
Overall, the simulations with CESM1, MPI-ESM and EC-Earth agree well with the evolution of water balance for the past 12 ka as reconstructed with our BRT-based calibration models (Fig. 2a). This agreement signifies an improvement over earlier simulations by TraCE-21ka and gives confidence in both the reconstructions and the numerical models. The initial wet conditions during the deglaciation period (12 to 11 ka) in the MW and GL region, and the following drying trend towards 9 ka in reconstructions is well captured by the CESM1 snapshot simulations. Conversely, the transient simulation with TraCE-21ka fall largely outside the uncertainty band of reconstructions throughout the entire period, with almost no changes in water balance.

The long-term wetting trend over the past 8 ka is well captured by MPI-ESM and EC-Earth simulations, separating these models from the more inconsistent hydrological trends presented by TraCE-21ka for the Holocene. Despite these advancements, some data–model discrepancies persist, such as the underestimation of early-Holocene dry conditions in the NE by CESM1, which may result from a poor spatial coverage of coastal areas in the model. Meanwhile, MPI-ESM and EC-Earth simulate slightly drier conditions than reconstructed from 8–7 ka for the MW, resulting in a stronger wetting trend than seen in the reconstructions.

The CESM1 simulations also show better agreement with the proxy reconstructions for Tjul than the transient simulations from TraCE-21k, MPI-ESM, and EC-Earth (Fig. 2b). While the MPI-ESM and EC-Earth simulations overcome the previous noted inconsistency in water balances seen in TraCE-21k, their simulated Tjul trends for the last 8 ka are significantly higher than in the reconstructions with up to +4.5 °C in the MW for EC-Earth and +3.5 °C for MPI-ESM) (Fig. 2b). Such deviations underscore the challenge of accurately reproducing Tjul trends in the transient simulations.

Sensitivity of vegetation and pollen to changes in water balance grounds our reconstruction, but the simulations enable us to disentangle the changes in precipitation versus evapotranspiration, which the proxy record cannot separately resolve. Figure 4 shows the simulated anomalies in annual precipitation (ΔP) and potential evapotranspiration (ΔPET) relative to preindustrial levels, which together inform the observed changes in water balance changes (Δ[P–PET]), which appear broadly consistent with the reconstructions in Fig. 2. Based on these hydrological anomalies, the wet conditions at 12 ka and 11 ka resulted primarily from reduced PET due to cooler and shorter than preindustrial summers, which compensated for a strong reduction in precipitation during this period (Fig. 4a). As the early Holocene climate warmed, however, the persistence of low precipitation favored drought conditions from 11 ka onward. For the mid-Holocene ~8 to 5 ka (Fig. 4b), the most severe drought conditions were driven by high PET during the substantially warmer than preindustrial summers. Conversely, the onset of neoglacial cooling and increased precipitation around 4 ka alleviated the drought conditions and allowed water balance to increase by 2 ka. It is noteworthy that the shift in PET from negative values before 8 ka to positive, but declining, afterward predominantly shaped the evolution of the annual water balance curves in our simulations, as the transient dynamics in PET overprinted the roughly linear increase of annual precipitation from 12 to 4 ka. This pattern underscores the dominant influence of PET, driven by thermal responses to ice sheet configuration and orbital forcing, on the regional water balance.

To investigate the large Tjul difference in simulations for the mid-Holocene, we also examine the data–model agreement for growing degree days (GDD), a measure integrating the biologically available heat across the growing season. For GDD5 (5 °C temperature threshold), transient simulations show a good agreement with reconstructions for the past 8 ka in GL and NE, while MPI-ESM and EC-Earth are again too warm over the MW (Supplementary Fig. 11). In contrast to the good agreement with Tjul, CESM1 performs worse for GDD5, indicating low GDD5 values for 12 to 9 ka, albeit following a comparable trend with reconstructions. There is an apparent paradox in the models showing consistent water balance, despite indicating much higher Tjul compared to proxies which would suggest a significant disagreement in the PET portion of the water balance equation. By design, water balance includes compensating effects between PET (driven by warm season temperatures above freezing) and the amount of annual precipitation. While all components are known in climate models including the warm season length contributing to PET, no direct proxy evidence is available for PET and annual P, making the paradox difficult to fully resolve. However, for GDD5, which characterizes the overall warm season instead of Tjul only, the agreement between proxies and the MPI-ESM and EC-Earth simulations improves considerably (Supplementary Fig. 11), which increases our confidence in the realism of the PET and precipitation contributions to the water balance of the past 8 ka indicated by MPI-ESM and EC-Earth (Fig. 4b). We hence consider it more likely that the water balance in the proxies and the models are correct and rather doubt the accuracy of the Tjul inferences. This also applies to CESM1, with the difference that here Tjul appears more similar with proxies than the warm-season quantity GDD5.

Sub-millennial climate events and periodicities
Drawing on the subset of well-dated and high-resolution fossil pollen records available for the MW and NE regions, our reconstructions confirm the presence of submillennial climate events and periodicities in temperature and moisture levels superimposed on the multi-millennial climate trends. At centennial smoothing bandwidths, both water balance and Tjul in all regions show significant variations (lower panels in Fig. 2a, b) which are, however, generally not coherent among regions (Fig. 5a). Also, the amplitude of the centennial shifts (up to ca. 50 mm; Fig. 5b) is smaller than that of the multi-millennial aridity trend. The outstanding, spatially coherent centennial event occurs at ca. 5.5–5.0 ka, initiating with a significant centennial-scale rise in water balance in all regions at ca. 5.5 ka (Fig. 2a) and leading to a transient multi-century moisture maximum (Fig. 5a). While we cannot rule out that this is only a coincidental alignment of the continuous sub-millennial variations seen in all clusters, the 5.5–5.0 ka event matches a wet period at 5.4 to 4.8 ka previously reported in the multiproxy data from the varved sequence of Elk Lake, Minnesota19, which interrupts the mid-Holocene period of aridity. Shuman29 also identified a period of increased moisture at 5.6–4.5 ka spanning the eastern North American mid-latitudes as the outstanding centennial-millennial scale deviation from the Holocene multi-millennial trends in this region. The amplitude of centennial moisture changes is overall stable through the Holocene in all regions (Fig. 5b). The one exception is a sharp apparent reduction in water balance in MW over the past 200 years (Fig. 5a), which is likely a spurious effect driven by the concurrent and probably anthropogenic upswings in Artemisia and Amaranthaceae pollen percentages (Supplementary Figs. 12–13).

Numerous significant periodicities in water balance and Tjul were found in wavelet analyzes, which were run on both the transient simulations and a subset of proxy reconstructions based on high-resolution pollen sequences available in the MW and NE clusters. These results are summarized in Fig. 6 as kernel density estimates fitted to the detected significant periodicities in the reconstructions (Fig. 6a, b) and the simulations (Fig. 6c, d). The periodicities are broadly similar for water balance (dotted lines) and Tjul (dashed lines), suggesting the influence of summer temperature variations on moisture dynamics, even at sub-millennial timescales (compare with Fig. 4). A notable clustering of periodicities around 0.2 ka is observed for both reconstructions and simulations in both MW and NE. However, the NE reconstructions also show a larger cluster between 0.4 and 0.5 ka, a feature not captured in the simulations. On the other hand, the simulations reveal significant multidecadal periodicities, which cannot be evaluated in the reconstructions given the insufficient temporal sampling of the fossil pollen datasets.

The ~200-year periodicities are also found in earlier simulation experiments with and without solar forcing, and may either reflect the role of solar forcing on North American hydroclimates, since cycles around 210 years has been linked to Suess/de Vries cycles, or represent internal variability39,45,46. In control simulations with the same EC-Earth model as used here for the past 8 ka, internal feedbacks are enough to sustain multi-centennial variability of the Atlantic meridional overturning circulation47 without any changes in external forcing. Simulated and reconstructed low-frequency variability is hence likely the result of combined internally and externally forced variability. The diminishing series of peaks at 0.5–1.2 ka periodicities in the MW proxy data (Fig. 6a) may represent subharmonics of the 200-year period. In NE reconstructions, however, the significant periodicities cluster around 400 and 500 years, especially for Tjul (Fig. 6b). Similar periodicities were observed in the marine sediment cores retrieved from the North Atlantic and Arctic oceans48–50 and also in earlier analyzes of lake-level reconstructions as well as the pollen data in the NE region29. Because the lakes of the NE cluster are located closer to the coastline, these periodicities (400–500 years) might be the effect of variations in ocean circulation on top of any solar forcing signal51, a feature that is not captured by proxies in the more continental MW region. An alternative explanation for the comparatively weaker 200-year periodicity in NE is the lower average resolution in NE fossil pollen sequences (mean gap between samples of 73.7 years over 8–0 ka) compared to MW (63.2 years), which could limit the detection of the 200-year cycle and allow the dominance of the 400–500 year cycle caused by oceanic forcing and/or arising as a subharmonic of the 200-year cycle.

Causes of Holocene millennial drought in eastern North America: comparison with prior studies
In our water balance reconstructions, the patterns for PC1 and PC2 sample values across time (Fig. 3c, e) parallel those found earlier for a dataset of pollen and alkenone based summer temperature reconstructions, lake-sediment stable isotope records, and lake-level and dust-flux moisture proxy data, spanning the mid-latitude US from the Rocky Mountains to the northeast Atlantic margin12. These PC1 and PC2 patterns closely correspond to two key forcings affecting this region: summer insolation (PC1; Fig. 5c) and the LIS deglaciation, which is closely associated with an increasing east-west moisture contrast (PC2; Fig. 5d)12. However, PC1 deviates from summer insolation before 8 ka due to the presence of a sizable LIS (Fig. 5c), while the PC2 trend reverses after 5 ka (Fig. 5d) due to a greater moisture increase in MW compared to NE, as the MW recovers from the later drought maximum in the mid-continent (Fig. 5a), which now starts to decrease the longitudinal moisture gradient captured by PC2.

Our reconstructed drought pattern, with the mid-continent drying over 10–8 ka while the NE gets wetter, also aligns with earlier reconstructions based on lake-level reconstructions, drought-indicating depositional hiatuses, and pollen abundances of drought-resistant and mesic taxa15. This spatiotemporal pattern has been hypothesized to be explained by LIS deglaciation, and a following waning of the glacial anticyclone and an increased influence of the Bermuda subtropical high, leading to a rerouting of the northward moisture advection towards the NE region and away from the mid-continent1,8,15,52. However, earlier climate simulations (represented here by TraCE-21ka; Fig. 2a) have shown inconsistent patterns of hydroclimatic change compared to proxy data9,10, and notably, indicating greater than modern precipitation and moisture balance (ratio of actual vs. potential evapotranspiration) at 6 ka through the U.S. mid-latitudes53. This poor data-model agreement has made it difficult to definitively attribute reconstructed hydroclimate variations to specific atmospheric dynamics.

Now, our climate reconstructions for the eastern North-American midlatitudes show good agreement with the latest generation model simulations regarding the first-order patterns of drought initiation and progress over the early and middle Holocene. Thus, we are now better positioned to discuss the changes in water balance and their drivers. As the data–model convergence resolves the fundamental mismatch seen in earlier studies, we suggest that the major spatiotemporal moisture patterns are adequately explained by the summer insolation trend combined with the Laurentide Ice Sheet deglaciation and its effect on atmospheric circulation.

Role of the Laurentide Ice Sheet (12–9 ka)
Based on the CESM1 simulations for the early Holocene, the distinct wet conditions at 12 to 11 ka can be largely attributed to a substantial reduction in PET, by ~30 to 40% (~180 to 230 mm) at 12 ka relative to preindustrial levels. This reduction in PET exceeds the decline in P during the same period, estimated at ~15 to 25% (~120 to 180 mm) (Fig. 4a). The strongly negative ΔPET is consistent with dominating cold and dry winds from northwest to east due to anticyclonic blocking over the LIS at 12 ka and 11 ka (Fig. 7a, d), where strong summer cooling with low PET over compensates the significantly reduced precipitation, resulting in positive water balance anomalies (Fig. 7b, e). This is consistent with several degrees colder Tjul in both the reconstructions (Fig. 2b) and the CESM1 simulations (Fig. 7c, f).

The dipole of very wet conditions in MW vs. dry in NE in the reconstructed water balance around 12 and 11 ka BP (Fig. 5a) can be explained by dynamical changes in atmospheric circulation linked to anticyclonic blocking over the ice sheet (Fig. 7a, d). To further explore the linkages between the regional water balance and large-scale atmospheric circulation, we performed a canonical correlation analysis (CCA) between large-scale variations of summer sea-level pressure (SLP) over North America and regional variations in water balance (P−PET) seen in the CESM1 runs for 12, 11 and 9 ka (Fig. 8). In general, CCA can identify pairs of linear combinations (canonical variates) from the two sets of variables that are maximally correlated with each other (see “Methods” for details). Applied to water balance in our study region, the CCA identifies strongly correlated patterns between variations in large-scale SLP and regional P−PET, suggesting that shifts in summer SLP are integral to the observed changes in water balance. As suggested by a decomposition of hydrological components of the water balance in model simulations (Fig. 4), the annual water balance variations can be expected to be linked to variations in summer and hence summer SLP, which is consistent with a comparable CCA for Europe54.

In our CESM1 simulation, the reconstructed water balance dipole around 12 ka BP is consistent with the pattern identified by CCA1 (7.4% explained variance) with SLP. The contrast between the wet MW and the dry coastal NE is caused by a meridional SLP anomaly pattern identified by CCA1 (9.5%) with positive anomalies over the ice sheet representing the anticyclonic blocking (Fig. 8a). Even CCA2 (8.4%) for the water balance yields such a dipole pattern, which is now linked to a meridional pressure pattern of SLP in CCA2 (17.1%) (Fig. 8e) that effectively steers humid air into the MW while dry air deflects precipitation in the NE (Fig. 7b). The west–east water balance dipole weakens at 11 ka, as the blocking anticyclone retreats northward with the retreating ice sheet, allowing more humidity to reach the Northeast (Fig. 8b). Along with the anticyclonic pattern, the meridional SLP pattern also moves northward but keeps humidity low to the south and southwest of the ice sheet (Fig. 8f).

By 9 ka BP, the anticyclonic variability pattern is replaced by more zonal SLP patterns due to the almost vanished ice sheet (Fig. 8c). The persisting dry conditions in NE seen in reconstructions (Fig. 5a) appear to be caused by the northward-shifted meridional SLP variability (indicated by CCA2 of SLP in Fig. 8g) compared to the modern climate (Fig. 8h). Note that even at 9 ka the co-variability pattern between SLP and water balance is very different compared to preindustrial, where the water balance is dominated by the cyclonic activity over north-central Canada (Fig. 8d). With the last remaining parts of LIS retreating to the far north at 9 ka, the climate over northeastern North America rapidly switches from the blocked deglaciation pattern seen at 12 and 11 ka (Fig. 7a, d) towards a regime driven by summer orbital forcing. This period marked the onset of increasing aridity, as summer warming (Fig. 7i) leads to an increase in PET that now exceeds P (Fig. 7h), intensified by the shift from cool maritime winds to warm continental westerlies due to the LIS retreat that removes the blocking anticyclonic flow from the region (Fig. 7g). Overall, these CCA results based on CESM1 confirm the hypothesis1,8,15,52 that the dipole in water balance in the early Holocene is indeed driven by anticyclonic blocking over the ice sheet.

Insolation-driven drought (8 ka onward)
The disappearance of the LIS by 8–6 ka (Fig. 5d) marked an end of the influence of the ice sheet on the spatiotemporal hydroclimate patterns, and from the mid-Holocene drought and its easing towards the present were largely driven by long-term changes in orbital summer insolation (Fig. 5c). Warming and drying rapidly increased to maximum drought conditions at 8 ka BP where the strong warming caused higher PET than today, while annual precipitation, despite a gradual rising trend, remained almost as low as during the late deglaciation (Fig. 4). This implies that drought is driven by both the higher PET, due to much warmer summer temperatures in response to orbital forcing, and clearly lower annual P than under preindustrial conditions.

The (relative) drought conditions through much of the Holocene, as delineated by both our reconstructions and the MPI-ESM and EC-Earth simulations, suggest that the recent and modern climate is unusually wet while drier conditions seem to be the norm during most of the Holocene. This period of enhanced aridity was largely due to elevated summer insolation, temperature, and evaporative demand. Future simulations project that, as greenhouse gas concentrations rise and global temperatures increase, mid continental North America will experience an increase in precipitation yet a decrease in plant-available soil moisture, because of enhanced evaporative demand5. It is hence plausible that the reversal of the natural neoglacial cooling observed over recent millennia through anthropogenic global warming might cause a return of mid continental aridity in eastern North America, for which the early to mid Holocene serves as a reasonable natural analog.

Methods
Pollen data
For pollen–climate calibration data, we use an eastern North American dataset37 derived from the North American Modern Pollen Database55,56, with the addition of 165 modern samples from the Neotoma Paleoecology Database57 (https://neotomadb.org/) and originating from more recent work22,58,59, to improve the coverage in the western prairie. The dataset (Fig. 1; Supplementary Data) includes 2419 surface pollen samples with climate means extracted from CRU CL v. 2.0 climate grids for 1961–199060. Correlation between water balance and Tjul is low (r=–0.11; Supplementary Fig. 14) due to the near-orthogonal gradient directions (north-south for Tjul, east-west for water balance), facilitating the independent modeling of temperature and moisture signals in the modern pollen data. The pollen taxonomy follows the 64-type list of ref. 38 for eastern North America, except for Pinus (pine) which is combined in a single column because the fossil pollen data used do not consistently differentiate between diploxylon and haploxylon pine. The calibration dataset excludes the region dominated by the ecologically distinct southeast pine species18,20. Importantly, this taxonomy does not include Ambrosia (ragweed), which is an important Holocene prairie taxon and indicator of hydroclimate variability26. However, Ambrosia has a strong human impact on its present distribution17 potentially biasing pollen–climate calibration models or leading to challenges in interpretation of paleoclimate reconstructions22.

To study the spatio-temporal patterns in drought, while pooling sufficient data to establish robust paleoclimate signals, we assembled fossil pollen sequences for the three regions (Fig. 1) from Neotoma. The following criteria were used to filter for acceptable fossil pollen sequences: at least 30 pollen samples and five absolute datings, to enable the detection and alignment of sub-millennial climate signals, a modern water balance within the range of the calibration data (−348 to 1343 mm) by a margin of at least 150 mm (Supplementary Fig. 14), to ensure that past variation relative to modern conditions can be reconstructed, and bottom age of at least 10 ka. The number of sites meeting these criteria was 20 for MW, 18 for GL, and 28 for NE (Supplementary Table 1). The MW and NE clusters each contained six sequences with considerably higher than average number of pollen samples (minimum 125) and 14C dates (minimum nine, excluding dates rejected by the authors), and we prepared separate reconstructions from these sets of well-dated high-resolution sequences. In MW, all six high-resolution sites use AMS 14C dates from plant macrofossils or charcoal, while the NE high-resolution sites largely use AMS dates from bulk sediment.

Paleoclimate reconstruction
The BRT models for pollen–water balance, pollen–Tjul, and pollen–GDD5 calibration were created using the R61 library gbm62 with the following settings: maximum number of trees = 3000, learning rate = 0.025, tree complexity = 4, bagging fraction = 0.5. For further details on these models, see ref. 37.

To align reconstructions from different modern climatic settings, all reconstructions were expressed as deviations from the site-specific mean over the past 4 ka. To equalize the impact of fossil sites on the reconstructed moisture and temperature levels regardless of their sampling resolution, we then interpolated all site-specific reconstructions at 50-year timestep. Synthesis reconstructions were then calculated for each region as the mean of all reconstructions, repeated 1000 times using bootstrap samples of all fossil sequences to calculate the 95% errors (2.5th and 97.5th percentiles) of the mean. SiZer maps63 (implemented using the R library SiZer64) were calculated based on the synthesis curves to identify significant rising and falling trends in the reconstructions when smoothed at a range of decadal-to-millennial bandwidths. To further analyze spatiotemporal patterns in the water balance reconstructions, we performed a principal components analysis (PCA) on the reconstructions from all sites covering the 11–0.25 ka timespan (53 sites out of 66)12.

In our reconstructions for the MW and NE regions, we primarily use the reconstructions based on the high-resolution subsets of sites. While the reconstructions using all sites are broadly similar and would not change the key interpretations regarding the patterns and drivers of multi-millennial moisture change, the all-site reconstructions appear to truncate the Holocene range of water balance variation, while the high-resolution reconstructions also reveal deeper sub-millennial moisture and temperature variations between the sites (Supplementary Fig. 15). However, the all-site reconstructions for MW and NE may be preferable for future studies concerned with multi-millennial moisture or temperature levels, which are likely to be sufficiently captured by all sites which passed our initial filtering, and we have thus also included the full results using all sites in the Supplementary Data.

To explore the effect of calibration model selection, we also prepared the water balance and Tjul reconstructions using MAT21, the most commonly used approach in pollen-based climate reconstructions in North America12,19,38. The MAT was implemented using the R library rioja65 and a weighted mean of five best analogs. The timing of the early-mid-Holocene multi-millennial droughts and the temperature maximum is broadly consistent between reconstructions prepared with BRT and MAT (Supplementary Fig. 16), showing that the qualitative patterns underpinning our main conclusions are robust to calibration method selection. However, the BRT reconstructions show a stronger peak warming in all sectors, while the broad scale Holocene drought is deeper in MW and GL but weaker in NE compared to MAT-based reconstructions. We primarily employ the BRT-based paleoclimate reconstructions in this work, because in cross-validations using a robust h-block scheme37, the BRT model performs considerably better for water balance (coefficient of determination 0.65 [BRT] vs. 0.56 [MAT]) with a smaller advantage in Tjul (0.88 [BRT] vs. 0.86 [MAT]) (Supplementary Fig. 17). Further, earlier h-block cross-validation experiments using a wide range of h values show that the performance advantage of BRT over MAT widens with large h, suggesting more robust performance of BRT with fossil samples with poor modern analogs.37

The BRT models gain additional support from an analysis of the relative contributions of the most important predictor taxa (Supplementary Figs. 12, 13) in the pollen–climate calibration data (Fig. 1). This analysis confirms the ecological realism of the BRT-based calibration models, with the water balance model relying heavily on the prairie forbs Artemisia (sagebrush) and Amaranthaceae (amaranth family). The five most important water balance predictors also include moisture-sensitive trees, including the comparatively drought-resistant Pinus (pine) but also the moisture-demanding Betula (birch) and Abies (fir)18,20, which the BRT models employ in this order as indicators of progressively wetter conditions (Supplementary Fig. 12). By comparison, the Tjul model uses largely different predictors (Supplementary Fig. 13) with a predominant driver in eastern Quercus (oak), a well-understood summer temperature-sensitive tree in eastern North America17,18,20,56. The water balance and Tjul calibration models thus have distinct structures (i.e., rely on different paleobotanical signals) and are consistent with prior ecological knowledge on indicator species for summer temperature and moisture availability.

Paleoclimate model simulations
We compare our reconstructions for water balance, Tjul, and GDD5 with a set of climate model simulations. The model outputs were extracted using boxes coinciding with the pollen site clusters (MW: 42–51°N, 92–100°W; GL: 40–47°N, 85–90°W; NE: 40–47°N, 69–75°W, excluding ocean), shown in Fig. 1. As a reference simulation covering the whole period, we reconsider the (quasi-)transient, fully coupled atmosphere-ocean simulation TraCE-21k that covers the full period from 22 ka to the modern climate66,67. The simulation incorporates changes in orbital and greenhouse gas forcing, adjustments of continental ice sheets, and sea-level rise, as well as meltwater fluxes triggering abrupt climate shifts in agreement with geological records. The CCSM3 model used for TraCE-21k has a relatively coarse spatial resolution of 3.75°×3.75° (~400 km) and was found to produce an inconsistent evolution of North American water balance in a previous study for the Holocene9,10.

To evaluate whether a higher spatial resolution alters the representation of atmospheric flow and related climate in response to the LIS during the period 12 to 9 ka, we analyze additional 100-year snapshot simulations conducted with CESM1.0.5 with a ~4× higher horizontal resolution of 0.9°×1.25° (~100 km) for the periods 12 ka, 11 ka, and 9 ka. The model setup and boundary conditions are adjusted to these periods for changes in radiative forcing, ice-sheet configurations, and sea-level change (described in more detail for 13 and 12 ka in ref. 69). The ocean and sea-ice states for these periods are prescribed from TraCE-21k using a mean annual cycle climatology calculated as ± 50 years around the target period. The high spatial resolution leads to fundamental differences in response to the presence of ice sheets over Europe68,69 compared to its parent model CCSM3 despite having identical radiative forcing and ocean states. Simulations with CESM1 for the deglaciation show a higher sensitivity for changes in the hydrological cycle and tropical convection with more extreme changes in seasonality compared to CCSM370,71.

For the remaining transient evolution of the Holocene, we make use of simulations covering the period 8 ka BP to modern conducted with the MPI-ESM1.2 model40,41 and EC-Earth3-veg-LR72,73. The MPI-ESM (Earth System Model of the Max Planck Institute for Meteorology) is forced by variations in orbital forcing, greenhouse gases and solar irradiance40, with an updated volcanic forcing as well as land-use and dynamical vegetation changes41 with a horizontal resolution of 1.875°×1.875° (nominal resolution of 1.5° for ocean and sea-ice). EC-Earth is forced by orbital and greenhouse gas changes at a horizontal model resolution of 1.125°×1.125° (1° for ocean and sea-ice).

For comparison with our reconstructions, we use simulated annual water balance, Tjul, and GDD5. Annual water balance is calculated from simulations as P−PET [mm/a] where P = annual total precipitation [mm/a] and PET = potential evapotranspiration. To be consistent with the calculation used in the pollen–climate calibration data37, we calculate PET as 58.93 × MABT [mm/a] where MABT is the mean annual biological temperature derived as the annual average of monthly means >0 °C. Similarly, GDD5 was calculated from model outputs using the same formula as for the pollen–climate calibration data, with mean daily temperatures interpolated from a sine regression fitted to the monthly mean temperature values37. Anomalies of water balance, Tjul, and GDD5 are calculated relative to the preindustrial climate (1850 AD) for CESM1 and relative to 0.25–0.75 ka (as done for reconstructions) for transient simulations of CCSM3, MPI-ESM and EC-Earth. Additional model output from CESM1 is used for dynamical analysis for sea-level pressure [hPa] and near-surface wind at the terrain following lowest sigma-hybrid coordinate level 992 (~20–50 m above ground). The pressure data for past climate states were corrected for systematic offsets by subtracting the difference of the global mean sea-level pressure between climate states in the past and preindustrial to avoid spurious offsets caused by non-dynamical changes of lower sea-level and glacio-isostatic land movements and ice sheets.

To investigate whether the presence of the retreating LIS in the early Holocene is indeed a controlling factor for regional changes in water balance12,15, we perform a Canonical Correlation Analysis74 (CCA) of simulated water balance P−PET variations with the large-scale atmospheric circulation over North America for the periods 12, 11 and 9 ka BP (Fig. 8). As annual precipitation P changes roughly linearly across the 12,000 years and evapotranspiration PET is ~0 during months ≤0 °C and hence dominated by summer changes (Fig. 4), the CCA is performed on regional annual P−PET vs. large-scale summer SLP using monthly data from CESM1.

The CCA results in two patterns showing maximum correlation among two geophysical fields X and Y. In our example X relates to SLP and Y relates to P−PET. Both variables are derived from CESM1 simulations for different periods of time, 12, 11 and 9 ka BP. The hypothesis is that basic (tele-) connection and structure between the different periods changed because of (profound) changes in background surface boundary conditions regarding continental eastern Laurentide ice sheets. These changes should then also be reflected in the canonical patterns between SLP (XCCA) and P−PET (YCCA).

A technical pre-processing of the original fields X and Y relates to an empirical orthogonal function (EOF)-truncation prior to carrying out the CCA. This simplifies the mathematical structures. In our case, five EOFs representing ~60% of total variance in P−PET and ~72% in SLP are used as input to generate the CCA patterns. One of the features of the CCA patterns is that their time coefficients are orthogonal, i.e. they do not show any temporal correlation.

Concerning the interpretation of the CCA results, the leading CCA1 patterns share the highest correlation among the variables X and Y, whereas the higher indexed CCA patterns show lower correlations. Therefore, higher indexed CCA patterns should be interpreted with care, especially when their canonical correlation indicates low values. Here, we focus on the first two leading patterns.

Analysis of periodicities
To assess the presence of significant periodicities, wavelet analyzes (Morlet wavelet) were performed for both paleoclimate reconstructions and the transient model simulations using the R library WaveletComp75. The wavelet analysis was performed individually on each of the water balance and Tjul reconstructions prepared from the subsets of high-resolution sites available for the MW and NE site clusters (12 sites total; Supplementary Table 1) using the BRT calibration models. For comparison, we also performed the wavelet analysis on the simulated water balance and Tjul anomalies from each transient model simulation (N=3) for the spatial domains of the MW and NE site clusters. In all wavelet analyzes, we considered the 8–0 ka timespan covered by all reconstructions and transient simulations. Prior to the wavelet analysis, the transient model simulations data were averaged using a 50-year running mean to mimic the 50-year resolution of the paleoclimate reconstructions, using the R library zoo76. To remove oscillations on the millennial timescale, all data (paleoclimate reconstructions and the transient model simulations) were filtered using a high-pass red noise Butterworth filter using the R library signal77 (cutoff = 1000 years). Finally, the significant periodicities (p<0.05) relative to red noise observed in the wavelet analysis were summarized in a kernel density plot.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Dear Editor,\nI would like to first apologize to the authors for delays in my submission of the current review. I have read the manuscript\nfrom Salonen and al. with great interest and found it interesting and worthy of publication in your journal after revision. While\nI have reservations on some methodological aspects described below that have to be addressed, I do not think it would\nfundamentally alter the conclusions of the paper, which rely on several lines of evidence. I have thus a few general\ncomments on the content and structure of the manuscript, and a list of detailed comments with line numbers.\nStructure wise, I think that the manuscript could be generally condensed a little both in terms of text and figures. I also note\nthat it seems to me that new results are integrated in the discussion. I don’t think the classic results-discussion separation is\nnecessary and that the sections could be restructured with more descriptive titles instead.\nMy general comments have to do mainly with (i) the handling of the reconstructions, specifically relating to the stacking and\nresolution, and (ii) with the analysis of periodicities.\n(i) If I understand well the stacking methodology, running 250-year windows at 50-year intervals were used to average all the\navailable data points. This means that a high-resolution series with say 5 points in the window would contribute 5 times\nmore than that of a low-resolution series with only1 point in the window, right? Not to mention that some series might not\neven contribute to some windows if some of their points are more than 250-year apart. This can lead to apparent non-\nstationarity in the variability as more estimates are available in the more recent period and thus, on one hand more noise is\naveraged out, and on the other hand, more smoothing occurs due to time-uncertainty. In both cases, we thus expect less\nvariability (and not more) for stacks of higher resolution series, unless the signal-to-noise ratio is high and the time-\nuncertainty is low. I would thus suggest to at least perform a temporal interpolation before averaging the series as was\ndone for example by https://doi.org/10.1038/nature25464.\nIs there evidence that BRT improves the signal-to-noise ratio (SNR) within clusters compared to MAT? It seems clear that\nthe variability of the BRT-based stacks is higher than the MAT-based ones, but is this due to higher SNR with BRT, or\nsimply because the variance of BRT reconstructions are generally higher than that of MAT ones? In other words, is the\nvariance ratio from individual series in clusters to their stack the same for both methods, or is it really higher (lower) for MAT\n(BRT)?\n(ii) While the analysis of periodicities seems to be computationally correct, I have reserves on the methodological choices.\nThe methods used have however been used quite frequently in palæoclimate studies and therefore my criticism is not\nnecessarily specific to the author’s work. First, I have reserves on the usage of wavelets to obtain time-dependant and\ntimescale-dependent variability estimates from palæoclimate data, the issue being that we are extracting a lot of information\nfrom relatively short and noisy timeseries, and changes in periodicity through time is likely random rather than expressing a\nreal change in the state. For such Holocene timeseries, resolving about one order of magnitude in frequencies, more robust\nresults would be obtained by looking at the wavelet (or regular power) spectrum for the whole time period (e.g. Fig.1 in\nhttps://doi.org/10.1038/s41561-022-01056-4).\nWavelet spectrogram can still be qualitatively informative in this case to see how the variability is distributed across time and\nfrequency, but to infer quantitative information requires some additional considerations. The null hypothesis has to describe\nthe “background variability” correctly. While an AR1 red noise is commonly implemented and used in palæoclimate analysis,\nit almost never describe the background variability satisfactorily and leads to spurious estimation of peak significance (e.g.\nhttps://doi.org/10.1175/JCLI-D-22-0011.1). Therefore, most, if not all, of the periodicities are likely not significant (and with\np<0.05 we still expect 5% of false positives). I also note here that I wonder how the periodicities were identified from the\nspectrogram as they are not stationary time, was that from a mean wavelet spectrum?\nI also have reserves with respect to the frequencies that could be resolved given the low- and high-pass filters that were\napplied. The pollen data is computed with 250-year sliding windows, which corresponds to a convolution with a rectangular\nwindow, or with the multiplication of a sinc function in the Fourier domain. What this means really, is that the cutoff frequency\naffected by the sliding window, effectively a not so ideal lowpass filter, is greater than 1/(500 years). Then on top of that, a\nhigh-pass filter is used to kill the power below 1/(1000 years), and thus we have a very narrow window between 1/500 and\n1/1000 years where variability estimates are less affected by the smoothings (albeit the running mean filter has a sizeable\nimpact even on the (1/1000 years) frequency). What is then the sense of the AR1 null hypothesis when the shape of the\nspectrum does not represent the signal but the smoothing filters right? In the end it looks like the identified periodicities are\nquite randomly distributed over the range of frequencies and I’d suspect that any similarity, albeit I fail to see the described\nbroad agreement, would be more likely an interaction between the shape of the filters and that of the null hypothesis.\nI think that in order to evaluate whether the water balance and July Temperature are linked, I would suggest a cross-spectral\nanalysis. If common periodicities do exist, they should be clearly evidenced by such analysis as significant levels of\ncoherency should be attain at the relevant frequencies.\nBelow are detailed comments with line and figure references:\nLine 72-76: A bit of a weird sentence, you say MAT and responses surfaces were used before, and that not MAT and other\nclassical methods are used. I'm not sure what the logic is, except to say that response surfaces might not be used so much\nanymore? I just think the sentence could be simplified to say \"different methods such as X and Y have been developed and\nused.\nLine 78: Last two sentences could be combined and shortened, and a summary sentence could be added to the paragraph\nwith the key takeaway point.\nLine 82: Is it really abrupt if it occurs over a millennium? Abrupt kinda implies a departure from more \"normal\" variability to\nsomething more extreme, I'm not sure if that's really the case, at least not when we look at the average of several\nreconstructions (which might be overly smoothed) but on the other hand single reconstructions are likely containing lots of\nnon-climatic noise.\nLine 94: The time is always ripe for a new and better analysis. I think a \"Why\" could be included here too already, \"because\nof new methods and data\".\nLine 95-96: \"have arisen\" sounds like they just appeared on their own, they were developed (by people).\nLine 118: transient climate model simulations? I would maybe use the phasing: \"an ensemble of N Holocene transient\nsimulations across 4 climate models\"\nLine 125: No need to repeat there are four sets of model simulations here, just \"the Holocene simulations\"\nLine 128: Finally instead Fourth?\nLine 129: Maybe skip the \"concerning the spatiotemporal patterns.... America\" and just say directly \"we test the hypotheses\nthat the patterns are linked to...\"\nFigure 1: Just wondering what the reason is for the lack of records between the Great Lakes cluster and the Northeast\ncluster? Lack of high resolution records for these areas or because the signal is too mixed there and not useful for the\nanalysis? I would also suggest to move the map to the extended data figures.\nLine 150: Typo \"the MD\" maybe \"in the MD\"?\nLine 151: Typo \"at a 10-7 ka\"\nLine 159: Unclear what means that the progress of the drought is well tracked by proxy data. Do the authors simply mean\nthat negative water balance anomalies are reconstructed for all regions from 10ka onward?\nFigure 2: Wondering if there are ways to simply the figure and caption. Could we show the average of MPI-ESM and EC-\nEarth together and maybe only show MAT in the extended data fig? Please indicate the bandwith h with units and without\nthe logarithm on the axis (the spacing can remain logarithmic, but the labels would be clearer if given directly in ka). On the\nother hand, are those SiZer maps really necessary? What do they tell us except what we can already see on the plot?\nThe GDD5 is maybe not necessary to be shown? Was there an additional CESM1 pre-industrial snapshot that was considered\nto calculate the anomalies? Could the models be combined as on Fig 4, and maybe keep TraCE-21ka for the supplement?\nFigure 2: Presumably the temperature and water balance ratio are relatively well correlated in the model simulations.\nTherefore, if water balance agrees with the reconstructions but not July temperature, I wonder whether a model that would\nsimulate correctly the July temperature would then end up overestimating the water balance?\nLine 190: Can the details of the comparison between BRT and MAT be moved to the methods or supplement to improve the\nreadability flow?\nLine 196: Maybe specify that the PCA is done across all the clusters.\nLine 198: Maybe write \"with the remaining components not exceeding 6.1%.\" rather than using a symbol \"< 6.1%\".\nLine 202-207: Could the description of PC2 be better and more simply explained? I just wonder a bit, what is the longitudinal\nwater balance contrast, or the period of maximum contrasts due to peak multi-millennial drought. I think this could be\nsaid more simply. I don't know if the PC loadings can be said to be moving with longitudes, I see mostly two clusters with positive\nloadings in the interior and negative on the coast, and not so much transition in between.\nLine 209: Are the two loadings covarying with the midwest sites doing so significantly? Given that the original series might\nhave different resolutions (and I imagine the PCA is performed after interpolation), high amplitude loadings could still be\nrandom or dominated by noise.\nFigure 3: What does ”scores for each sample of water balance reconstructions” mean? Revise and clarify sentence.\nLine 272: The shift from negative to positive PET anomalies looks like model bias, it's doubtful that a continuous curve could\ngo through. Are there no simulations of CESM1 that overlap with the other runs to evaluate this? Again I wonder about how\nthe CESM1 anomalies were computer and whether this could be related.\nFigure 4: It's interesting to note that the CESM1 simulations show a wetter midwest while the transient ones show a wetter\nnortheast. I wonder what we learn from the difference regional curves here other than this reversal between the models\nwhich doesn't seem to be interpreted. If we cannot trust the regional differences from the simulations, I would show instead\none curve for the average across all three regions. Again, it would be interesting if an addition CESM1 simulation at 8ka and\nthus overlapping with the transient simulations was performed in order to see if the discrepancy is a model bias. I think the\nfigure could be improved in general better evidence the key takeaway.\nFigure 5: Revise figure. Either show the average spectra or the coherency to show common periodicities.\nLine 344: Yet to prove with a coherency analysis.\nLine 353: Missing closing bracket.\nLine 362: Again unsure about the periodicity argument.\nFigure 6: It seems that all the data shown here has been shown before on different panels, it would be nice if somehow we\ncould condense the storyline and avoid duplication. (a) Seems like the Great Lakes and Northeast might have some anti-\ncorrelation, might show up on coherency analysis. (b) Why the shading under the curve? What high-pass filtering was done?\nLooks like the cutoff frequency might be too high, I would favour only detrending the very slow trend (could be done by\nsubtracting a loess smoothing with a long span (on the order of a couple thousands of years). (c) Seems that the\ndeglaciation could explain the first part of both curves, and that the second part would be the insolation trend.\nLine 387: How were the climatic impact of vertical land movements estimated? Why 0.5 K? Unclear to me how vertical land\nmovements affect the climate.\nLine 388: Why are they expected to negatively bias the HTM amplitude and delay reconstructed HTM?\nLine 397: What is the contradiction?\nLine 399: So are you suggesting that the July temperature reconstructions are of poor quality? Then why use them at all?\nLine 404: What is meant by \"significant\" variations? Significant should be used in the context of testing.\nLine 405: Are they however coherent between regions? I imagine that there is a high level on noise on single\nreconstructions, but still wondering on the robustness of centennial scale variations.\nLine 412: Could the reason that the 5.5-5.0 ka wettening is visible across the three regions simply that it has higher\namplitude and thus creates a stronger signal above the noise in all three regions, whereas elsewhere the centennial\nvariability is masked by proxy noise?\nLine 424: Were the series first interpolated before the average? The higher number of samples in the recent part means that\nmore noise will be averaged out and therefore, the difference in variability might be an indication that more proxy noise\nremains in the earlier part. In this sense, the higher data resolution should average more noise and decrease variability.\nLine 429: 200-year is below the effective resolution of the series, which were produce with 250-year sliding windows. It does\nnot make sense to interpret that timescale.\nLine 430: Or maybe the significant periodicities around the 200-year frequencies correspond to a transition timescale of the\nAR1 null hypothesis where a mismatch is more likely to occur and lead to spuriously significant periodicities? That would be\nmy guess.\nLine 474: Period missing.\nLine 493: What is meant by comparing large scale SLP and regional water balance? Do I understand right that in one case\nthe CCA is performed over a region from longitude 170W to 40W and Latitude 20N to 75N and the other 100W to 40W and\n32N to 58N only over land? These seem to be very arbitrary choices to me. Of course this type of analysis can be sensitive\nto the data that is included and choosing different domains can change the dipole spatial patterns observed.\nFigure 8: Are we really learning something from all those figures? Could it be simplified? I am not so familiar with the details\nof CCA, but I wonder if the correlations are significant? I think the Methods on CCA could better describe the procedure.\nLine 539: I struggle with this statement given that the rise in PET from 9ka to 8ka in model simulations was inferred from a\nchange from CESM1 to transient simulations.\nLine 616: Using 250-year windows means that the centennial scale variations, which were analyzed and interpreted, will be\nvery smoothed. I think any variations below that should not be interpreted, at least in terms of their amplitude as they will be\ndamped by construction. I think it would be better to interpolate to 50-year resolution all the series and stack them. The\nseries could then be smoothed when plotting timeseries to evidence slower variability, but unsmoothed series should be\nused for wavelet/spectral analysis.\nMethods\nLine 572: comma before \"which\".\nLine 573: Split sentence before \"however\".\nLine 584: How were the cluster limits chosen? I wonder whether the gap between the GL and NE areas are a choice or a\nlack of records.\nLine 587: Well-dated is good, but what does that mean concretely in terms of time uncertainty? Generally, I’d still expect a\ntime uncertainty on the order of hundreds of years, meaning that stacking nearby sites would lead to a smoothing of the\nsignal due to ‘misalignment’. The effect can be seen on surrogate data created from model simulations with realistic time\nuncertainty and the resolution of pollen records on Ext. Data Fig 1 in https://doi.org/10.1038/s41561-022-01056-4 . Of course\nthe effect should be smaller in this case given the selection for high quality records, but I still suspect it to be important for the\ncentennial variablity. I also wonder here the impact of sediment mixing on the effective resolution of the high-resolution\nrecords\nLine 616: Were the reconstructions interpolated before taking the rolling mean, as was done for example by Marsicek et al\n2018, albeit here there would be no spatial gridding since we are only interested in stacking the clusters.\nLine 625: Well it wouldn't be suprising that the stack of all records and high-resolution records show similar features as they\nare both using the high-resolution data, and if there is not interpolation before stacking done, then the high-resolution one\nwill contribute more samples to the rolling mean and have a higher impact on the all record stack.\nLine 642: How long were there simulated snapshots?\nLine 657: \"gasses\"\nLine 675: So the adjustment is a constant offset over the entire field? Meaning that it doesn't affect the results in any way\nexcept shifting the colour scale on the maps a little?\nLine 686: I think it would be useful to add a sentence describing what the CCA is about, broadly.\nLine 696: Does that then mean that the explained variance on figure 8 are with respect to 60% and 72% instead of 100%?\nLine 703: Are the CCA2 really worth interpreting then? It wouldn't hurt to simplify the analysis and figure if possible.\nLine 709: My general take on wavelets is that I don't think the quality of the palæo data here is good enough (even though it\nis indeed high quality pollen data) to estimate time-dependant and timescale-dependant variability estimates at once.\nLine 718: Specify the cutoff frequency of the Butterworth filter. Unclear to me why high-pass filtering was necessary for such\nan analysis. On one hand, one can just ignore the spectral estimates for longer than millennial frequencies if one does not\nwish to study them, and on the other hand, aren't millennial scale variations better resolved by the reconstructions and more\ninteresting than the more uncertain faster variability? I would perform such analysis without the high-pass filtering and\ndetrending instead, which could be linear detrending or sinusoidal detrending for precession, or no detrending (e.g.\nhttps://doi.org/10.1038/s41561-022-01056-4).\nSupplementary Figures 9-16: Unclear to me why the absolute amplitude of the wavelet power levels between simulations\nand reconstructions are generally an order of magnitude apart from each other. I think it’s because wavelet coefficients are\nnot normalized by the frequency right, so the discrepancy would come from the different resolutions of the input timeseries?\nIn any case, being able to compare the amplitudes would be nice.", "reviewer_2": "Patterns and drivers of Holocene moisture variability in mid-latitude eastern North America\nSubmitted to Nature: Communications\nThis manuscript presents a thorough and creative evaluation of Holocene hydroclimate patterns represented by fossil pollen\nassemblages and the latest generation of model simulations mid-latitude eastern North America. The manuscript helps\naddress the longstanding disagreement between model simulations and paleoclimate proxy reconstructions over the course\nof the Holocene in North America. I expect that the findings will be of keen interest to the paleoclimate community. The\nauthors use of a variety of methods and types of data (boosted regression tree proxy calibrations, time slice model\nsimulations, transient model simulations) is a strength of the paper.\nThe paper would benefit from better organization, especially in the Discussion, a s well as a more detailed description of\ntheir rationale for carrying out certain analyses and a more thorough explanation of their methodology. I outline the main\nissues with the Discussion and Methods below. I also include some more minor questions/suggestions for the Figures\nand for specific lines.\nDiscussion\nI recommend adding sub-headers to the Discussion section to help organize the various pieces of analysis.\nA more detailed discussion of how the new BRT-based proxy reconstruction compares to that of the traditional MAT\napproach would be welcome. Does this study indicate that BRT should be widely adopted instead of MAT?\nThe authors note that the latest generation of transient climate models do not agree with the reconstructions when it comes\nto TJul but say that the data-model match is much-improved for the entire growing season heat sum (GDD5) (which is also\ntrue for TraCE-21ka). This appears to be truer for the Great Lakes and Northeast region than it is for the Midwest (Figure 2c),\nwhich should be acknowledged/explained. Additionally, it would be worthwhile for the authors to describe why it is that the\ntransient simulations produce better agreement for GDD5 than they do for TJul, and why this observation does not hold true\nfor the CESM1 time slice experiments. Finally, is there a reason that the GDD5 values are not included in Supplementary\nTable 2? (Lines 394 – 402)\nThe trend reversal in PC2 after 5ka really stands out (Figure 3, Figure 6). Can the authors clarify what drives PC2? Is it the\nextent of the LIS up until it is gone at ~5 ka at which point the main driver becomes the east-west moisture contrast? Or is it\nthe east-west moisture contrast across the entirety of the record and PC2 aligns so closely with LIS extent from 11 to 5 ka\nbecause the LIS extent controls the east-west moisture contrast? (Line 452 – 454)\nIt is somewhat surprising to the reader that the authors focus the second half of the discussion entirely on the time slice\nexperiments from CESM1. I understand that they focus on CESM1 because they are probing the dynamical relationship\nbetween hydroclimate and the LIS, but the authors could do a better job of outlining the rationale for this focus in Lines 456 –\n469. Additionally, it would be interesting to see the same dynamical treatment shown in Figures 7 and 8 given to time slices\nextracted from the transient simulations. These analyses carried out for the transient simulations (individually or as an\nensemble) for 8, 7, and/or 6 ka would be useful for further elucidating the influence of the LIS as it disappears. (Lines 490 –\n537).\nThe CCA analysis is confusing and would benefit from a more fleshed-out description. It would be helpful is the authors\nexplicitly stated what the percentages in the Figure 8 sub-headers mean in terms of the CCA.\nMethods\nWhat criteria do the authors use to subset the “high-resolution” fossil pollen records?\nThe authors should outline which model grid cells are used to calculate mean values for each of the three regions. They\ncould also mark the boundaries of these regions (as defined for the spatial extents in the model output) in Figure 1.\nThe authors should describe how they calculate the water balance synthesis curves in Figure 6.\nFigures\nFigure 4: I’d recommend showing the points for 12,11, and 9 ka snapshots from CESM1, and potentially doing away with the\nlines connecting these points, to underscore that these are not transient simulations.\nWhy combine MPI-ESM and EC-Earth for this figure while keeping them separate in Figure 2? Does this model ensemble\noffer something that the individual models do not? The specifics of how this ensemble is calculated should be described in\nthe text. For example, were the models interpolated to the same spatial resolution prior to averaging? Would it be worthwhile\nto also include TraCE-21ka in Figure 4?\nWhat explains the big difference in sign of the precipitation anomaly between CESM1 at 9 ka and the mean of MPI-ESM1\nand EC-Earth at 8 ka? Is this due to differences in the models (CESM1 vs the transient simulations) or interpreted as an\nactual feature of the climate system?\nCan the authors include some metric of uncertainty in this figure, especially for the average of MPI-ESM and EC-Earth.\nFigure 7: The authors should consider graphically representing areas where the difference between the time slice\nexperiment (12, 11, or 9 ka) and the PI are not significant. This can be done using a Student’s T Test and shown on the\nfigures by graying or whiting out the anomaly (usually around the transition from positive to negative).\nFigure 8: I am confused as to what exactly these anomalies are relative to. Can the authors describe how the data in this\nfigure are calculated? Are all these anomalies statistically significant / is there a way to test and communicate this? I\nunderstand that CCA1 doesn’t always have highest explained variance for P-PET or SLP individually, but is it strange that\nthere are times when the values for both P-PET and SLP are lower than those for CCA2?\nLines\nLine 238: What is meant by “a poor spatial coverage of coastal areas in the model.”?\nLine 353: Add closing parenthesis. Add comma after “In contrast”\nLine 432: Can the authors briefly elaborate on the “internal variability” that may be contributing to the ~200-year periodicities\nthat are observed in both reconstructions and simulations?\nLine 495: Refer to Figure 8 at end of sentence introducing CCA.\nLine 637: Reword the beginning of the sentence about CCSM3.", "reviewer_3": "This paper presents an impressive synthesis of pollen records and model analyses, leveraging new regression tree methods\nto attempt to distinguish moisture balance from temperature signals in the eastern US. The boosted regression tree\napproach appears to provide new insights not available via the common modern analog technique. The authors then\ncompare these results to transient climate model simulations to speculate about the mechanisms responsible for past\nchanges. Overall, I really like this paper, but I think it could benefit from tightening the story and honing the ‘big picture’\nmessage of the work. I think these could be tackled via major revisions. My overall feeling is that the main strength of this\npaper comes from the interesting new reconstruction. The model analyses are interesting, but also raise many questions and\nsome confusion in those analyses detracts from the overall strength of this work. I would not be afraid to restructure this\npaper as a ‘data driven’ paper.\nBelow are some overall points to consider as well as detailed comments on the text.\nFirst, I think the organization of the paper would benefit from a combined results and discussion section rather than having\nthese sections distinct. As it stands, there are topics that should traditionally be included in the discussion in the Results,\nand the organization gets a little bit messy. Especially for a short format paper, combining these would be valuable.\nSecond, the overall conclusion of the paper could be stronger. The abstract reads more like a literature review rather than\nthe abstract of a short format paper with a ‘punchy’ conclusion, and a similar comment could be applied to the paper’s\ndiscussion. My personal thought is it might be better to focus either on the impact of long-term forcing on the PCA/water\nbalance reconstruction, or on the question of periodicities. As it stands, I don’t know that the model analysis is ‘enough’ to\nsupport the conclusion that summer insolation and the Laurentide are driving the patterns you see. But, it is worth noting that\nTrace21ka has single forcing simulations (e.g. with ONLY orbital changes, or with ONLY ice volume). :\nhttps://www.earthsystemgrid.org/project/trace.html\nIt is also worth taking a step back and starting with evaluating the PI model simulations’ climatology. Which ones are too\nrainy over your regions of interest? Too hot? Are the winds too strong relative to observations? This might provide some\nadditional insights about why model performance varies. It would also be good to make a table as to what forcing datasets\nwere used in each set of transient runs - the change in different generations of ice sheet reconstructions (e.g. 5G, 6G, etc)\nwill affect simulations even of the same model.\nFinally, there are other ‘data’ products out there for deglacial climate variables (see Osman et al 2022 in Nature for the data-\nassimilation approach). My sense is that that product performs well for temperature but might not perform well for\nhydroclimate-related variables. It would be worth comparing your reconstructions (real world data! So more likely to\nrepresent the truth) to the DA product for P and P/ET during the deglaciation. It might help emphasize the innovative nature\nof the reconstructions that you’ve produced.\nIntroduction\nLine 126-131: Instead of just reviewing the methods and approaches you used in the paper, provide a signpost of the overall\nconclusion of the paper.\nResults:\nFigure 2: Minor point, but I find it hard in some panels to distinguish the lines for EC-Earth and Trace 21ka. Make sure you\ncheck colorblindness compatibility of the figure, and perhaps use line weights/dash patterns to increase understandability.\nLine 183: This applies to the paper overall, but it is worth defining what you mean by ‘drought’ - for paleo records like this,\nwe tend to define drought as just below average P or P-E, while the term can mean something very different in the modern\nliterature (e.g. below a certain percentile of precipitation).\nLine 196: I think the PCA presents one of the most interesting results in this paper, and it is worth really strongly pinning\ndown what the driver of this pattern is. I am not sure I understand why PC2 is Laurentide deglaciation if it goes ‘up’ again\nafter 5 ka?\nLine 225: I am not sure I agree that there is a good comparison with the evolution of water balance between the models and\nthe data. CESM1 seems to overlap the water balance reconstructions, but seems to be very off in terms of GDD so this to me\nsuggests that it may be getting the ‘right’ answer in terms of hydroclimate change for the wrong reasons. I am therefore not\nsure how much to trust the CCA. The logic behind relying on CESM1 for interpretation needs to be clearer. An\nalternate/complementary analysis might involve comparing Trace to MPI and EC Earth for an interval like 7 ka, when the\nformer is not dry enough and the latter two are too dry. Difference fields between these runs might reveal why the models\ndiffer in skill at simulating deglacial climate change. See comments above about some basic diagnostic analyses as well.\nFigure 6: In panel d, if this mode reflects Laurentide deglaciation, why does it go back ‘up’ slightly after 5 ka?\nLine 305: I am not sure if the section on periodicities is a good idea for this paper. It seems unrelated to the overall story of\nthe paper, and could be a separate, more robust analysis in and of itself.\nDiscussion:\nSuggest combing with results section, much of this discussion could be integrated through the prior section\nLine 394: See comments above about model-data comparison.\nConclusions:\nThe concluding paragraph of your paper focuses on the MPI, EC Earth, and new reconstructions. If this is the focus, it would\nbe important to then focus on these model simulations in Figures 7 and 8.\nI would also check about the length limits for this journal - it seems like there are a lot of figures, and the text is quite long for\na short format journal."}, "author_response": "PATTERNS AND DRIVERS OF HOLOCENE MOISTURE VARIABILITY IN MID-\nLATITUDE EASTERN NORTH AMERICA\nResponse to Referees\nWe quote the reviewers’ full comments and describe the changes done to address each point. We\nare grateful to all reviewers for their thoughtful and detailed comments and we believe that the\nmanuscript has been strengthened by addressing them.\nReviewer 1: Dear Editor,\nI would like to first apologize to the authors for delays in my submission of the current\nreview. I have read the manuscript from Salonen and al. with great interest and found it\ninteresting and worthy of publication in your journal after revision. While I have\nreservations on some methodological aspects described below that have to be addressed,\nI do not think it would fundamentally alter the conclusions of the paper, which rely on\nseveral lines of evidence. I have thus a few general comments on the content and structure\nof the manuscript, and a list of detailed comments with line numbers.\nStructure wise, I think that the manuscript could be generally condensed a little both in\nterms of text and figures. I also note that it seems to me that new results are integrated in\nthe discussion. I don’t think the classic results-discussion separation is necessary and that\nthe sections could be restructured with more descriptive titles instead.\nWe have changed the manuscript to use a combined “Results and Discussion”, with thematic\nsubheadings.\nThe switch to Results and Discussion eliminated about 300 words due to less repetition,\nand we have mainly used this space to better describe two important elements of the paper: (A)\nthe significance of the difference in data-model match for T and GDD5, and (B) the CESM1\njul\nequilibrium runs and their CCA analyses at the end of the Results and Discussion.\nReviewer 1: My general comments have to do mainly with (i) the handling of the\nreconstructions, specifically relating to the stacking and resolution, and (ii) with the\nanalysis of periodicities.\n(i) If I understand well the stacking methodology, running 250-year windows at 50-year\nintervals were used to average all the available data points. This means that a high-\nresolution series with say 5 points in the window would contribute 5 times more than that\nof a low-resolution series with only 1 point in the window, right? Not to mention that some\nseries might not even contribute to some windows if some of their points are more than\n250-year apart. This can lead to apparent non-stationarity in the variability as more\nestimates are available in the more recent period and thus, on one hand more noise is\naveraged out, and on the other hand, more smoothing occurs due to time-uncertainty. In\nboth cases, we thus expect less variability (and not more) for stacks of higher resolution\nseries, unless the signal-to-noise ratio is high and the time-uncertainty is low. I would thus\nsuggest to at least perform a temporal interpolation before averaging the series as was\ndone for example by https://doi.org/10.1038/nature25464.\nThank you for this and other comments; we appreciate the careful review. We have done the\ntemporal interpolation of the individual site-specific reconstructions as requested by Reviewer 1,\nusing a 50-year time step.\nImportantly, we note that in our original submission, we already used the data with\ninterpolated values at 50-year time steps for both the PCA and the wavelet analyses. Thus, the\nPCA and wavelet analyses are entirely unaffected by the new synthesis curve methodology.\nthe methodology of the paper.\nWith the new stacked reconstructions based on the interpolated site reconstructions, we\ncontinue to use bootstrapping to calculate the 95% errors for the ensemble mean, using the values\ninterpolated for the individual sequences for each time point. To summarize multi-centennial\ntrends, we also plot the running 5-point mean of the ensemble means.\nFrom a methodological angle, we agree with Reviewer 1 that the benefit of the\ninterpolation of the individual sequences is that it equalizes the impact of higher and lower\nresolution fossil sites on the synthesis curve. This can be desirable especially in the study of\nmillennial trends, because (A) even the lower-resolution sites should be able to resolve the\nmillennial trends and (B) because individual sites may be subject to persistent biases in pollen\ntaxon percentages due to taphonomic factors such as basin size. When considering centennial\nfeatures, one might counter that ultimately a higher-resolution sequence provides more\nobservations on fossil pollen variation and should be given proportionally more weight.\n(Anecdotally, in another ongoing study we have faced a situation where interpolating the data\nsequences destroys a rather beautiful record of the 8.2k event provided by the highest-resolution\nfossil sites.) However, because our largest focus here is on millennial spatiotemporal patterns,\nand because our Wavelet analyses used to study periodicities employ the individual sequences\nanyway (instead of stacked reconstructions), on balance we do see the benefit of the interpolation\nfor this study.\nThe figure below compares the new vs. old reconstruction synthesis curves:\nTo summarize, these are the major changes in the synthesis curves:\n● Overall, the synthesis curves do not change greatly, and there are no changes to key\ninterpretations of the paper.\n2\n● Some changes are seen in the millennial moisture levels, especially in MW and NE where\nthe use of the subsets of only 6 high-resolution pollen sequences makes the reconstructed\nmoisture level sensitive to cases where some comparatively high- resolution sequence\nfalls towards the edges of the ensemble. In particular, a deepening early-Holocene\nmidcontinent-coastal contrast is seen due to wetter conditions in MW and a deeper\ndrought in NE. This reinforces one of the major observations in the paper.\n● Over the past 8 ka, the fit between water balance reconstructions and simulations (shown\nin the revised Fig. 2) improves further for both GL and NE, with GL showing deeper\nmid-Holocene drought, and NE somewhat shallower drought, which agrees with the new-\ngeneration transient models (EC-Earth and MPI-ESM).\nWe have also recalculated the synthesis curves of the transient climate simulations shown in Fig.\n2, using 50-year binning of the annual values instead of the 250-year binning used in our original\nsubmission, to remain as consistent as possible with our new reconstruction synthesis curves.\nThis change mainly affects the width of the 95% error bands and is not important for the paper.\nReviewer 1: Is there evidence that BRT improves the signal-to-noise ratio (SNR) within\nclusters compared to MAT? It seems clear that the variability of the BRT-based stacks is\nhigher than the MAT-based ones, but is this due to higher SNR with BRT, or simply because\nthe variance of BRT reconstructions are generally higher than that of MAT ones? In other\nwords, is the variance ratio from individual series in clusters to their stack the same for\nboth methods, or is it really higher (lower) for MAT (BRT)?\nOur choice of BRT for this paper draws on the cross-validation (CV) experiment shown in\nSalonen et al. (2019) which compared eight calibration models (including BRT and MAT) which\nuses a robust h-block CV scheme, indicating BRT as the best-performing moodel. That paper\nfurther showed that the BRT performance advantage over MAT would further increase with very\nhigh h values, which mimics the situation faced in paleo-reconstructions with poor modern pollen\nanalogs.\nWe acknowledge that even this approach cannot fully represent the problem of\nreconstruction from fossil data, where we suspect the reliability of methods remains, to a\nsignificant extent, unknowable. Hence we decided to also include the full MAT results to show\nthat the main conclusions of the paper are not sensitive to calibration method choice.\nFrom a practitioner viewpoint, our preference for BRT is also affected by the\nexceptionally good tools BRT provides to verify that the paleo-botanical signals (i.e. variation of\nspecific pollen types) used to reconstruct a given climate variable is consistent with ecological\nknowledge. In this study, such analysis of the BRT models gives strong ecological validation for\nthe BRT calibration models (Supplementary Figs. 12 and 13).\nWe have edited the reconstruction section of Methods to more clearly state this\nreasoning.\nReviewer 1: (ii) While the analysis of periodicities seems to be computationally correct, I\nhave reserves on the methodological choices. The methods used have however been used\nquite frequently in palæoclimate studies and therefore my criticism is not necessarily\nspecific to the author’s work. First, I have reserves on the usage of wavelets to obtain time-\ndependant and timescale-dependent variability estimates from palæoclimate data, the\nissue being that we are extracting a lot of information from relatively short and noisy\ntimeseries, and changes in periodicity through time is likely random rather than expressing\na real change in the state. For such Holocene timeseries, resolving about one order of\nmagnitude in frequencies, more robust results would be obtained by looking at the wavelet\n(or regular power) spectrum for the whole time period (e.g. Fig.1 in\nhttps://doi.org/10.1038/s41561-022-01056-4).\nWe calculated the wavelet for the past 8 ka, the same as in Hébert et al. (2022). We just removed\nthe Early Holocene because of the influence of the LIS in the lakes.\nhow the variability is distributed across time and frequency, but to infer quantitative\ninformation requires some additional considerations. The null hypothesis has to describe\nthe “background variability” correctly. While an AR1 red noise is commonly implemented\nand used in palæoclimate analysis, it almost never describe the background variability\nsatisfactorily and leads to spurious estimation of peak significance (e.g.\nhttps://doi.org/10.1175/JCLI-D-22-0011.1). Therefore, most, if not all, of the periodicities\nare likely not significant (and with p<0.05 we still expect 5% of false positives). I also note\nhere that I wonder how the periodicities were identified from the spectrogram as they are\nnot stationary time, was that from a mean wavelet spectrum?\nThe periodicities are indeed from the mean wavelet spectrum.\nIt is generally difficult, and for specific cases almost impossible, to know the “true”\nbackground variability which makes AR1 the most plausible choice. A better grip on the color of\nproxy noise might help in future studies (e.g. McPartland et al. (preprint):\nhttps://cp.copernicus.org/preprints/cp-2024-73/). Using many records as done in this study has\nthe advantage that we can be more confident about the robustness of the results if local records\nconsistently yield comparable frequencies independent of the type of significance testing.\nOut of curiosity, we tested 44 different methods/parameter choices (described below) of\ndescribing the “background variability” and, therefore, to calculate the significance of the\nobserved periodicities. Most of the periodicities were robust enough that they appear in most, if\nnot all, the methods. The tested methods included:\n● White noise; N = 1\n● \"shuffle\" (shuffling the given time series); N = 1\n● “Fourier.rand\" (time series with a similar spectrum); N = 1\n● “AR” (AutoRegression), with “p” (the number of past observations required to define the\nauto regression) varying between 1 and 5; N = 5\n● “ARIMA” (AutoRegressive Integrated Moving Average), with “p” (the number of past\nobservations required to define the auto regression) varying between 0 and 5, and “q”\n(the number of past observations required in the moving average) also varying between\n0 and 5; N = 36\nHowever, red noise (AR1) is usually implemented in time series analysis because its power\nspectrum is weighted toward low frequencies. Therefore, for the manuscript, we decided to use\nonly the tests based on red-noise methods, i.e., AR-1, ARIMA(101), ARIMA (102), ARIMA\n(103), ARIMA (104) and ARIMA (105). For the final figure in the manuscript, we kept only the\nperiodicities that were significant in more than 66% of the tested methods (i.e., in at least 4 out\nof 6 AR-based methods).\nWe have revised the Supplementary Table 3 to show all periodicities observed originally\nusing AR-1 method, and in addition, how many times (in %) the peaks were significant\nconsidering the 6 AR-based tested methods. From the results shown in Fig. 6, we have then\nfiltered out the periodicities which were significant in three or fewer of the tested methods. This\nresults in relatively minor changes in the kernel density plots fitted to the periodicities, and\ndoesn’t impact our main conclusions about the periodicities.\nReviewer 1: I also have reserves with respect to the frequencies that could be resolved\ngiven the low- and high-pass filters that were applied. The pollen data is computed with\n250-year sliding windows, which corresponds to a convolution with a rectangular\nwindow, or with the multiplication of a sinc function in the Fourier domain. What this means really,\nis that the cutoff frequency affected by the sliding window, effectively a not so ideal lowpass\nfilter, is greater than 1/(500 years). Then on top of that, a high-pass filter is used to kill\nthe power below 1/(1000 years), and thus we have a very narrow window between 1/500\nand 1/1000 years where variability estimates are less affected by the smoothings (albeit\nthe running mean filter has a sizeable impact even on the (1/1000 years) frequency). What\nis then the sense of the AR1 null hypothesis when the shape of the spectrum does not\nrepresent the signal but the smoothing filters right? In the end it looks like the identified\nperiodicities are quite randomly distributed over the range of frequencies and I’d suspect\nthat any similarity, albeit I fail to see the described broad agreement, would be more likely\nan interaction between the shape of the filters and that of the null hypothesis.\nFirst, we note an apparent misunderstanding represented by the statement “the pollen data is\ncomputed with 250-year sliding window” and we apologize for lack of clarity in the original\nsubmission. We did indeed perform our multi-site synthesis reconstructions in our original\nsubmission with the 250-year rolling window – although this approach was changed for this\nrevision, as discussed above. However, the wavelet analysis was not performed on these synthesis\ncurves, but on the individual, underlying pollen-based reconstructions, and even then, specifically\nonly on the 12 highest-resolution reconstructions available for the Midwest and Northeast\nclusters. As we note in the periodicity subsection of Results and Discussion, these pollen sites\nhave an average time resolution of about 60–70 years. These reconstructions had been further\ninterpolated at 50-year time step prior to the wavelet analyses; as we note above, this was also\ndone for the data incorporated in the synthesis curves in this revision.\nWe have edited the Methods section about the wavelet analysis, to be clearer on the\nsubset of data used for the analysis.\nReviewer 1: I think that in order to evaluate whether the water balance and July\nTemperature are linked, I would suggest a cross-spectral analysis. If common periodicities\ndo exist, they should be clearly evidenced by such analysis as significant levels of\ncoherency should be attain at the relevant frequencies.\nOur main focus was to identify if the water balance and the July temperature presented significant\nperiodicities. Because our intent was not to test if they varied simultaneously, we prefer to keep\nthe wavelet analysis of the two parameters separate. Additionally, a cross-spectral analysis would\nadd even more discussion to an already very complex paper.\nReviewer 1: Below are detailed comments with line and figure references:\nLine 72-76: A bit of a weird sentence, you say MAT and responses surfaces were used\nbefore, and that not MAT and other classical methods are used. I'm not sure what the logic\nis, except to say that response surfaces might not be used so much anymore? I just think\nthe sentence could be simplified to say \"different methods such as X and Y have been\ndeveloped and used.\nLine 78: Last two sentences could be combined and shortened, and a summary sentence\ncould be added to the paragraph with the key takeaway point.\nWe have edited this section, with a simplified list of reconstruction approaches as suggested by\nthe reviewer, and in the following sentence clarified our intended meaning, i.e. contrasting the\nquantitative moisture reconstructions with indirect measures of moisture variability.\nReviewer 1: Line 82: Is it really abrupt if it occurs over a millennium? Abrupt kinda\nimplies a departure from more \"normal\" variability to something more extreme, I'm not\nsure if that's really the case, at least not when we look at the average of several\nreconstructions (which might be overly smoothed) but on the other hand single\nreconstructions are likely containing lots of non-climatic noise.\nWe agree the word “rapid” is not needed here; removed.\nReviewer 1: Line 94: The time is always ripe for a new and better analysis. I think a \"Why\"\ncould be included here too already, \"because of new methods and data\".\nsynthesis and analysis of the eastern North American hydroclimate (because of these new\nadvances which we then describe).\nReviewer 1: Line 95-96: \"have arisen\" sounds like they just appeared on their own, they\nwere developed (by people).\nRephrased → “have been developed”.\nReviewer 1: Line 118: transient climate model simulations? I would maybe use the\nphasing: \"an ensemble of N Holocene transient simulations across 4 climate models\"\nRephrased → “an ensemble of four climate model simulations for the Holocene”.\nReviewer 1: Line 125: No need to repeat there are four sets of model simulations here,\njust \"the Holocene simulations\"\nDone.\nReviewer 1: Line 128: Finally instead Fourth?\nDone.\nReviewer 1: Line 129: Maybe skip the \"concerning the spatiotemporal patterns....\nAmerica\" and just say directly \"we test the hypotheses that the patterns are linked to...\"\nDone.\nReviewer 1: Figure 1: Just wondering what the reason is for the lack of records between\nthe Great Lakes cluster and the Northeast cluster? Lack of high resolution records for\nthese areas or because the signal is too mixed there and not useful for the analysis? I would\nalso suggest to move the map to the extended data figures.\nThe cluster limits reflect the (rather long) development history of this paper. We were originally\nworking with the Midwest cluster only, representing the much-studied forest-prairie ecotone and\na real clustering of fossil sites. After being satisfied that the BRT water balance models are also\nuseful within the wetter forest ecosystems, we introduced the Northeast cluster which also had a\nreal concentration of published high-resolution pollen data.\nThe Great Lakes cluster was the last addition. Here, the gap eastward before reaching the\nNortheast cluster is indeed real – after filtering the fossil pollen sites available on Neotoma for\nour required pollen sample and dating counts, no data were available to expand the Great Lakes\ncluster further eastward.\nWe agree Fig. 1 could potentially be in the Supplement. But given the relaxed article\nformat of Nature Communications, within which we also comfortably fit, we prefer to keep the\nmap figure, as it nicely grounds the introductory paragraphs about North American environmental\nhistory and the data used in this paper. Note that we have now also added boxes in the map used\nfor calculating the data regions from climate models (see another reply below, p. 17).\nReviewer 1: Line 150: Typo \"the MD\" maybe \"in the MD\"?\nFixed.\nReviewer 1: Line 151: Typo \"at a 10-7 ka\"\nReviewer 1: Line 159: Unclear what means that the progress of the drought is well tracked\nby proxy data. Do the authors simply mean that negative water balance anomalies are\nreconstructed for all regions from 10ka onward?\nRephrased to say that we mean the lack of specific biases due to poor modern analogs.\nReviewer 1: Figure 2: Wondering if there are ways to simply the figure and caption.\nCould we show the average of MPI-ESM and EC-Earth together and maybe only show MAT\nin the extended data fig? Please indicate the bandwith h with units and without the\nlogarithm on the axis (the spacing can remain logarithmic, but the labels would be clearer\nif given directly in ka). On the other hand, are those SiZer maps really necessary? What\ndo they tell us except what we can already see on the plot? The GDD5 is maybe not\nnecessary to be shown? Was there an additional CESM1 pre-industrial snapshot that was\nconsidered to calculate the anomalies? Could the models be combined as on Fig 4, and\nmaybe keep TraCE-21ka for the supplement?\nWe broadly agree with the above, and have done significant simplifications to Fig. 2:\n● The MAT reconstructions have been moved to supplement (new Supplementary Figure\n16)\n● The GDD5 reconstructions have also been moved to supplement (new Supplementary\nFigure 11).\n● We have done the suggested change to the SiZer panels, showing the smoothing\nbandwidth directly in ka. We have, however, opted to keep the SiZer maps in the figure.\nMethodologically, we’re quite fond of the SiZer because it simultaneously provides an\nassessment of significant trends on different timescales, while doing smoothing at a full\nrange of bandwidths and without requiring a subjective decision on the bandwidth to use\n(as is the case with e.g. LOESS).\nWith the overcrowding of the figure now greatly eased, we have opted to keep the transient\nmodels separated.\nAbout CESM1, there was indeed a separate PI control run; the results from this run are\nshown in Supplementary Table 2. We have added a note about the CESM1 PI run in the Fig. 2\ncaption.\nReviewer 1: Figure 2: Presumably the temperature and water balance ratio are relatively\nwell correlated in the model simulations. Therefore, if water balance agrees with the\nreconstructions but not July temperature, I wonder whether a model that would\nsimulate correctly the July temperature would then end up overestimating the water balance?\nThe water balance integrates different aspects of being a function of annual precipitation and\npotential evapotranspiration (PET). PET is generally a function of temperatures above freezing\npoint but includes both, how warm it is overall and how long the warm season is. As P and PET\nmay have compensating effects or not, the correlation is not necessarily obvious or stationary.\nSee also in another reply further below (p. 22) and how P vs. PET differs over time in Fig. 4.\nReviewer 1: Line 190: Can the details of the comparison between BRT and MAT be moved\nto the methods or supplement to improve the readability flow?\nAgreed. We have moved the BRT-vs-MAT comparison away from the main text. There is a new\nSupplementary Figure 16 comparing the BRT and MAT reconstructions. We run through the\ndifferences to MAT in the Methods, after presenting the BRT reconstruction approach.\nReviewer 1: Line 196: Maybe specify that the PCA is done across all the clusters.\nReviewer 1: Line 198: Maybe write \"with the remaining components not exceeding 6.1%.\"\nrather than using a symbol \"< 6.1%\".\nDone.\nReviewer 1: Line 202-207: Could the description of PC2 be better and more simply\nexplained? I just wonder a bit, what is the longitudinal water balance contrast, or the\nperiod of maximum contrasts due to peak multi-millennial drought. I think this could be\nsaid more simply. I don't know if the PC loadings can be said to be moving with longitudes,\nI see mostly two clusters with positive loadings in the interior and negative on the coast,\nand not so much transition in between.\nLine 209: Are the two loadings covarying with the midwest sites doing so significantly?\nGiven that the original series might have different resolutions (and I imagine the PCA is\nperformed after interpolation), high amplitude loadings could still be random or\ndominated by noise.\nWe agree that the PC2 loadings are better described as positive interior vs. negative coastal\nclusters. We have significantly edited this paragraph, to simplify the statement about the spatial\npattern in loadings. We have also rearranged the wording to first present the site loadings and\nthus the interpretation of each PC, and only then describe the temporal changes in sample scores.\nThe panels in the associated Fig. 3 have been reordered accordingly.\nReviewer 1: Line 209: Are the two loadings covarying with the midwest sites doing so\nsignificantly? Given that the original series might have different resolutions (and I imagine\nthe PCA is performed after interpolation), high amplitude loadings could still be random\nor dominated by noise.\nWhether a higher order PCA mode is meaningful is based on the explained variance of that mode\nand whether its interpretation can be linked to a “real” phenomenon. Based on the eigenvalues\nrepresenting the explained variances in Figure 3a, this is clearly the case for the first two modes.\nApart from this, PCA loadings are typically not tested for statistical significance as other\nstatistical measures. The loadings only say how strong the represented variability is within a given\nmode and the sign is used to indicate co-variability in phase or anti-phase. The lower the specific\nloadings, the less important within a given mode. Overall, PC1 and PC2 explain a significant\namount of the overall variability and can be physically explained as in Fig. 5\nReviewer 1: Figure 3: What does ”scores for each sample of water balance\nreconstructions” mean? Revise and clarify sentence.\nWe have significantly clarified the Fig. 3 caption.\nReviewer 1: Line 272: The shift from negative to positive PET anomalies looks like model\nbias, it's doubtful that a continuous curve could go through. Are there no simulations of\nCESM1 that overlap with the other runs to evaluate this? Again I wonder about how\nthe CESM1 anomalies were computer and whether this could be related.\nThere are two main problems here:\n1. Unfortunately, there are no transient simulations with sufficiently high resolution\nconnecting the early Holocene with the rest of the Holocene. While TraCE simulations do\ncapture the whole period, the very coarse spatial resolution and potentially other reasons\n8\nlead to fundamentally inconsistent results of that model with other models and proxies. We\nnow highlight the use of different models by splitting up Fig. 4 for the change in models.\n2. The transient high-resolution simulations of MPI-ESM and EC Earth starting at 8 ka BP\nstart from initial conditions that do not include any memory in the ocean from coming out\nof the ice age. This might explain the big jump relative to CESM1, however, there might\nalso be other reasons that could only be resolved with more transient simulations for the\nwhole period. For the specific region, the switch from still having an ice sheet (CESM1)\nto its demise (other models) may further contribute to the jump. The main motivation to\nuse the set of models included in the paper is to highlight the role of ice sheets (CESM1)\nin the early Holocene and the transient evolution afterwards (other models).\nReviewer 1: Figure 4: It's interesting to note that the CESM1 simulations show a wetter\nmidwest while the transient ones show a wetter northeast. I wonder what we learn from\nthe difference regional curves here other than this reversal between the models which\ndoesn't seem to be interpreted. If we cannot trust the regional differences from the\nsimulations, I would show instead one curve for the average across all three regions.\nAgain, it would be interesting if an addition CESM1 simulation at 8ka and thus\noverlapping with the transient simulations was performed in order to see if the discrepancy\nis a model bias. I think the figure could be improved in general better evidence the key\ntakeaway.\nBased on the proxy reconstructions, the water balance should differ regionally in the early\nHolocene as simulated by CESM1 in response to ice sheets. This west-east gradient changes with\nthe disappearance of the ice sheet in the proxies consistent with the transient climate models\nstarting at 8 ka. So the difference between CESM1 and the other models are real changes\nconsistent with the proxy data and are not dominated by the switch of models. Please note that\nthe model data is presented as anomalies relative to the pre-industrial simulation of each model.\nHence, systematic differences between models are largely already removed by this approach.\nReviewer 1: Figure 5: Revise figure. Either show the average spectra or the coherency to\nshow common periodicities.\nThe coloured (dashed and dotted) vertical lines in Fig. 5 (Fig. 6 in our revision) indicate the mean\nwavelet spectrum of each dataset. Because of the amount of analysed data (2 regions with 6 lake\nrecords and 3 model simulations per region, and 2 climate variables from each dataset or model),\nwe decided to present all observed significant periodicities as the vertical lines and integrate them\nusing a Kernel density plot. Presenting each average spectra individually would represent adding\na panel made of 36 different plots, which would be even more difficult to identify similar\nsignificant periodicities.\nReviewer 1: Line 344: Yet to prove with a coherency analysis.\nSee other reply below (p. 10); we have softened the statement on the coherency of the centennial\nevents.\nReviewer 1: Line 353: Missing closing bracket.\nFixed.\nReviewer 1: Line 362: Again unsure about the periodicity argument.\nThis passage was removed in combining the Results and Discussion.\nReviewer 1: Figure 6: It seems that all the data shown here has been shown before on\ndifferent panels, it would be nice if somehow we could condense the storyline and avoid\n9\nduplication. (a) Seems like the Great Lakes and Northeast might have some anti-\ncorrelation, might show up on coherency analysis. (b) Why the shading under the curve?\nWhat high-pass filtering was done? Looks like the cutoff frequency might be too high, I\nwould favour only detrending the very slow trend (could be done by subtracting a loess\nsmoothing with a long span (on the order of a couple thousands of years). (c) Seems that\nthe deglaciation could explain the first part of both curves, and that the second part would\nbe the insolation trend.\nWe would favor keeping this figure. First, we think the summary of water balance patterns in all\nclusters in a single panel is helpful for the Discussion and to summarize the paper. Also the\ncentennial variability and the comparisons vs. forcings do add significant new results.\nWe have done the requested change to the method of calculating centennial variability,\nusing the residuals of a 2-ka LOESS smoother.\nNotably, with the changes both to the method used to analyse the centennial variability\nand to the underlying climate curves (due to the new methodology described above), the results\ndo change quite significantly:\n● The 8.2k event doesn’t stand out as clearly, and we have removed to discussion regarding\nthis event.\n● The 5.5k wet shift still stands out as a centennial feature shared by all clusters, and the\ndiscussion regarding this is retained. However we now note the coherency of the 5.5k\nevent between regions could well be coincidental.\n● There’s no longer an evident shift in the degree of centennial variability over the\nHolocene. We note this briefly in the revised manuscript.\nThe disappearance of known events like 8.2k due to using a different filtering or interpolation\nmethod highlights that these choices should depend on the research question.\nIn sum, we are significantly shortening the discussion on centennial variability. This also\ngives us space to expand the discussion of the drivers of millennial change in a later part of the\nResults and Discussion (see below).\nWe have opted to not do the suggested coherency analyses. With the overall direction\nbeing towards expanding the discussion of the drivers of millennial change (and better discussion\nof Figs. 7 and 8) – while staying within Nature Communication’s length limits – we favor not\nexpanding but rather slightly shortening the discussion on centennial variability.\nSee comment below (p. 16) about the interpretation of the PC1 and PC2.\nReviewer 1: Line 387: How were the climatic impact of vertical land movements\nestimated? Why 0.5 K? Unclear to me how vertical land movements affect the climate.\nLine 388: Why are they expected to negatively bias the HTM amplitude and delay\nreconstructed HTM?\nThe idea here is that because of the isostatic readjustment, our proxy sites don’t capture a\npaleoclimate record from a stable elevation. Effectively, the isostatic forces carry the paleo-\nthermometer from one elevation to another over the millennia, and hence, the reconstruction\ndoesn’t represent the climate change that would have been observed at a site of stable elevation.\nRelevant for this work, in a region of forebulge collapse, the sites would have been at\nincreasingly higher elevation towards the early Holocene, leading to an increasing negative bias\nin reconstructed temperature towards the early Holocene.\nHence, in a region of forebulge collapse, the site would have also been at a higher\nelevation during a mid-Holocene temperature maximum, making the reconstructed temperature\namplitude vs. preindustrial lower than would’ve been the case at constant elevation. Also,\nbecause the negative bias in temperature increases towards the early Holocene, this delays the\narrival to the mid-Holocene temperature peak.\nThis quick sketch perhaps helps illustrate the idea:\nforebulge collapse in the Great Lakes region, which can be visually estimated at slightly under\n100 m based on a figure in the Clark et al. (2012) paper we cite. We thus come to a rough estimate\nof 0.6 ºC for the maximum bias – we have changed this estimate from the rougher estimate “~\n0.5 ºC” given in our original submission, based on the cited lapse rate.\nWe have done multiple edits to this paragraph which hopefully clarify our reasoning.\nReviewer 1: Line 397: What is the contradiction?\nWe have significantly expanded this paragraph to explain the implications of the proxy-data\nmatch achieved for T and GDD5.\njul\nReviewer 1: Line 399: So are you suggesting that the July temperature reconstructions\nare of poor quality? Then why use them at all?\n“Poor quality” is perhaps a harsh way of putting it, but we do consider that for a vegetation-based\nproxy a variable like GDD5 encompassing more of the growing season might be more optimal.\nIn practice, however, we’re often stuck with T /MTWA to enable synthesis with other proxies\njul\nand climate modelling, and also to allow comparison with the plethora of earlier studies working\nwith T /MTWA. The problem of using peak warm or cold temperatures goes beyond the current\njul\nstudy and touches upon the general problem of potential seasonal biases, i.e,. if changes in\nseasonality and/or its effect is not known.\nReviewer 1: Line 404: What is meant by \"significant\" variations? Significant should be\nused in the context of testing.\nThis sentence was removed in combining the Results and Discussion; however this referred to\nthe significant deviations of the derivative from zero, as determined by the SiZer analysis.\nReviewer 1: Line 405: Are they however coherent between regions? I imagine that there\nis a high level on noise on single reconstructions, but still wondering on the robustness of\ncentennial scale variations.\nThis was a miswording and meant to say “coherent between regions”. This is corrected in the\nrevised section on sub-millennial variations.\nthree regions simply that it has higher amplitude and thus creates a stronger signal above\nthe noise in all three regions, whereas elsewhere the centennial variability is masked by\nproxy noise?\nWith the new synthesis curves, we have softened the statement on the 5.5k event (see reply above,\np. 10).\nReviewer 1: Line 424: Were the series first interpolated before the average? The higher\nnumber of samples in the recent part means that more noise will be averaged out and\ntherefore, the difference in variability might be an indication that more proxy noise\nremains in the earlier part. In this sense, the higher data resolution should average more\nnoise and decrease variability.\nAs noted, we have revised the synthesis curves to use interpolated data.\nReviewer 1: Line 429: 200-year is below the effective resolution of the series, which were\nproduce with 250-year sliding windows. It does not make sense to interpret that timescale.\nLine 430: Or maybe the significant periodicities around the 200-year frequencies\ncorrespond to a transition timescale of the AR1 null hypothesis where a mismatch is more\nlikely to occur and lead to spuriously significant periodicities? That would be my guess.\nSee reply above (p. 5) about the nature of the data used for the wavelet analysis.\nReviewer 1: Line 474: Period missing.\nFixed.\nReviewer 1: Line 493: What is meant by comparing large scale SLP and regional water\nbalance? Do I understand right that in one case the CCA is performed over a region from\nlongitude 170W to 40W and Latitude 20N to 75N and the other 100W to 40W and 32N to\n58N only over land? These seem to be very arbitrary choices to me. Of course this type of\nanalysis can be sensitive to the data that is included and choosing different domains can\nchange the dipole spatial patterns observed.\nWe described the CCA in more detail in the methods sections which, however, comes only at the\nend of the manuscript. We now add two sentences to better describe the rationale of using the\nCCA and the choice of the regions in the revision version already briefly in the Results and\nDiscussion and refer for details to the Methods. As written in the latter, in general, CCA finds\npairs of linear combinations (canonical variates) from the two sets of variables that are maximally\ncorrelated with each other.\nIn our case, we use CCA to find the modes of maximum correlation between a local\nvariable – here the small domain of the water balance where the proxy sites are located – with\nthe large-scale atmospheric driver here represented by continental-scale SLP over whole North\nAmerica. So the domain of the water balance has to be local to be comparable with the proxies.\nOtherwise, the CCA would no longer explain the local variability of the target region of the small\ndomain. In general, the selection of the domains for the CCA is always somewhat arbitrary, but\nwe tried to account for the large-scale character of SLP changes using a larger domain and the\nlocal-scale character of the water balance using a smaller domain. The smaller domain was\nspecifically selected to appreciate the concrete location of the proxy sites used in the study.\nReviewer 1: Figure 8: Are we really learning something from all those figures? Could it\nbe simplified? I am not so familiar with the details of CCA, but I wonder if the correlations\nare significant? I think the Methods on CCA could better describe the procedure.\nabove. Despite the noisy nature of water balance, the correlations with the large-scale SLP\nvariability are generally quite high as shown by the numbers in Fig. 8. Interestingly, the\ncorrelations are systematically higher in the early Holocene (r>0.4) than during the PI, possibly\ndue to the role of ice sheets in the past which is absent today. We will add to the figure caption\nthat all correlations are statistically significant at the 99% level (e.g. the smallest r=0.34 in CCA2\nfor PI with p<<0.01 (p=0.00053825)).\nReviewer 1: Line 539: I struggle with this statement given that the rise in PET from 9ka\nto 8ka in model simulations was inferred from a change from CESM1 to transient\nsimulations.\nSee another reply above (p. 8). The general switch from low PET in the early Holocene results\nfrom still having an ice sheet and hence cooler temperatures nearby. This effect is gone when the\ntransient models start at 8 ka BP.\nReviewer 1: Line 616: Using 250-year windows means that the centennial scale variations,\nwhich were analyzed and interpreted, will be very smoothed. I think any variations below\nthat should not be interpreted, at least in terms of their amplitude as they will be\ndamped by construction. I think it would be better to interpolate to 50-year resolution all\nthe series and stack them. The series could then be smoothed when plotting timeseries to\nevidence slower variability, but unsmoothed series should be used for wavelet/spectral\nanalysis.\nAs noted above (p. 2), in the revision we have done the interpolation of the individual sequences\nto 50-year resolution before stacking the reconstructions.\nReviewer 1: Methods\nLine 572: comma before \"which\".\nDone.\nReviewer 1: Line 573: Split sentence before \"however\".\nDone.\nReviewer 1: Line 584: How were the cluster limits chosen? I wonder whether the gap\nbetween the GL and NE areas are a choice or a lack of records.\nAddressed above (p. 6).\nReviewer 1: Line 587: Well-dated is good, but what does that mean concretely in terms of\ntime uncertainty? Generally, I’d still expect a time uncertainty on the order of hundreds of\nyears, meaning that stacking nearby sites would lead to a smoothing of the signal due to\n‘misalignment’. The effect can be seen on surrogate data created from model simulations\nwith realistic time uncertainty and the resolution of pollen records on Ext. Data Fig 1 in\nhttps://doi.org/10.1038/s41561-022-01056-4 . Of course the effect should be smaller in\nthis case given the selection for high quality records, but I still suspect it to be important\nfor the centennial variablity. I also wonder here the impact of sediment mixing on the\neffective resolution of the high-resolution records\nWe agree on all counts. Overall, our revised manuscript deals somewhat less with the sub-\nmillennial features seen in the stacked synthesis reconstructions. Meanwhile, as noted above (p.\n5), the analysis of periodicities doesn’t use the stacked reconstructions at all, but the individual\nreconstructions from the subset of high-resolution sites.\nmean, as was done for example by Marsicek et al 2018, albeit here there would be no\nspatial gridding since we are only interested in stacking the clusters.\nLine 625: Well it wouldn't be suprising that the stack of all records and high-resolution\nrecords show similar features as they are both using the high-resolution data, and if there\nis not interpolation before stacking done, then the high-resolution one will contribute more\nsamples to the rolling mean and have a higher impact on the all record stack.\nAddressed above (p. 2); now doing the interpolation of reconstructions before stacking.\nReviewer 1: Line 642: How long were there simulated snapshots?\nAs described in Schenk et al. (2018) and Kuang et al. (2012), the simulations were run for 150\nyears but only the last 100 years are used for analysis. We have added this to the Methods text.\nReviewer 1: Line 657: \"gasses\"\nFixed → “gases”.\nReviewer 1: Line 675: So the adjustment is a constant offset over the entire field? Meaning\nthat it doesn't affect the results in any way except shifting the colour scale on the maps a\nlittle?\nYes, subtracting the systematic model differences (offset) in global mean sea-level pressure\nbetween past and pre-industrial climate state experiments helps avoid spurious results. For\naccurate dynamical interpretation and comparison of surface or sea-level pressure, it is essential\nthat the surface or sea-level heights are consistent. Given the significant changes in sea levels and\nice sheets in the past which are incorporated in CESM1, removing the systematic offset is\nnecessary to achieve this consistency.\nReviewer 1: Line 686: I think it would be useful to add a sentence describing what the\nCCA is about, broadly.\nWe have added this in the revised version.\nReviewer 1: Line 696: Does that then mean that the explained variance on figure 8 are\nwith respect to 60% and 72% instead of 100%?\nThe explained variances in Fig. 8 are calculated relative to 100% of the input data used in the\nCanonical Correlation Analysis (CCA). While the input data for the CCA represents a subset of\nthe total simulation data (due to the removal of \"noise\" by the Empirical Orthogonal Functions\n(EOFs)), the variances shown in Fig. 8 are based on the proportion of variance explained within\nthis subset. Therefore, the percentages (60% and 72%) reflect the explained variance relative to\nthe input data used in the analysis, not the total amount of variance represented by the original\nsimulation data.\nReviewer 1: Line 703: Are the CCA2 really worth interpreting then? It wouldn't hurt to\nsimplify the analysis and figure if possible.\nThe high canonical correlations and large explained variances of CCA2 clearly require it to be\nincluded. As written in the text, it is also physically-dynamically required because CCA1 reflects\nthe dominating large-scale dominance of pressure variations over the ice sheet while CCA2\nreflects the dependence on the strength and position of zonal pressure gradients.\npalæo data here is good enough (even though it is indeed high quality pollen data) to\nestimate time-dependant and timescale-dependant variability estimates at once.\nSee reply above (p. 5) about the nature of the data used for the wavelet analysis.\nReviewer 1: Line 718: Specify the cutoff frequency of the Butterworth filter. Unclear to\nme why high-pass filtering was necessary for such an analysis. On one hand, one can just\nignore the spectral estimates for longer than millennial frequencies if one does not wish to\nstudy them, and on the other hand, aren't millennial scale variations better resolved by the\nreconstructions and more interesting than the more uncertain faster variability? I would\nperform such analysis without the high-pass filtering and detrending instead, which could\nbe linear detrending or sinusoidal detrending for precession, or no detrending (e.g.\nhttps://doi.org/10.1038/s41561-022-01056-4).\nThe cutoff frequency was 1/1000 years, to remove the multi-millennial frequency. We have added\nthis information to the manuscript Methods.\nThe high-pass filter was necessary because the magnitude of the millennial time scale is\ndisproportionately larger than the sub-millennial time scale. We tested running the wavelet\nanalysis without removing the multi-millennial frequency, but all the sub-millennial oscillations\nwere masked. Therefore, to be able to resolve the not so well studied sub-millennial oscillations,\nwe decided to apply the high-pass Butterworth filter.\nReviewer 1: Supplementary Figures 9-16: Unclear to me why the absolute amplitude of\nthe wavelet power levels between simulations and reconstructions are generally an order\nof magnitude apart from each other. I think it’s because wavelet coefficients are not\nnormalized by the frequency right, so the discrepancy would come from the different\nresolutions of the input timeseries? In any case, being able to compare the amplitudes\nwould be nice.\nIndeed, the wavelet coefficients were not normalized by the data frequency. The pollen data has\na 50-year resolution, while the model simulations’ data has a 1-year resolution. This difference\nmagnifies the amplitude of the wavelet power levels from the model data when compared to the\npollen data.\nReviewer 2: This manuscript presents a thorough and creative evaluation of Holocene\nhydroclimate patterns represented by fossil pollen assemblages and the latest generation\nof model simulations mid-latitude eastern North America. The manuscript helps\naddress the longstanding disagreement between model simulations and paleoclimate proxy\nreconstructions over the course of the Holocene in North America. I expect that the\nfindings will be of keen interest to the paleoclimate community. The authors use of a variety\nof methods and types of data (boosted regression tree proxy calibrations, time slice model\nsimulations, transient model simulations) is a strength of the paper.\nThe paper would benefit from better organization, especially in the Discussion, as well as\na more detailed description of their rationale for carrying out certain analyses and a more\nthorough explanation of their methodology. I outline the main issues with the Discussion\nand Methods below. I also include some more minor questions/suggestions for the Figures\nand for specific lines.\nDiscussion\nI recommend adding sub-headers to the Discussion section to help organize the various\npieces of analysis.\nDiscussion in this revision. We have also added sub-headers to the Results and Discussion.\nReviewer 2: A more detailed discussion of how the new BRT-based proxy reconstruction\ncompares to that of the traditional MAT approach would be welcome. Does this study\nindicate that BRT should be widely adopted instead of MAT?\nWe have clarified our rationale for preferring BRT over MAT – i.e., mainly cross-validation\nperformance, though acknowledging its limitations when assessing paleo-reconstructions – in an\nearlier reply (p. 3). A more thorough assessment of reconstruction methods we consider beyond\nthe scope of this manuscript, rather on this we refer to our earlier experiments published in\nSalonen et al. (2019). Instead we have, following another comment (p. 7), moved the assessment\nof the MAT reconstructions entirely to the Methods and the supplement.\nReviewer 2: The authors note that the latest generation of transient climate models do not\nagree with the reconstructions when it comes to TJul but say that the data-model match is\nmuch-improved for the entire growing season heat sum (GDD5) (which is also true for\nTraCE-21ka). This appears to be truer for the Great Lakes and Northeast region than it is\nfor the Midwest (Figure 2c), which should be acknowledged/explained. Additionally, it\nwould be worthwhile for the authors to describe why it is that the transient simulations\nproduce better agreement for GDD5 than they do for TJul, and why this observation does\nnot hold true for the CESM1 time slice experiments. Finally, is there a reason that the\nGDD5 values are not included in Supplementary Table 2? (Lines 394 – 402)\nWe acknowledge the worse data-model fit for Midwest in the GDD5 paragraph of our revised\nResults and Discussion.\nIt is somewhat ironic that our aim was to resolve the model-proxy discrepancy of\nHolocene water balance (which we do) but now end up with a disagreement in the magnitude of\nwarmer July temperatures. The disagreement of July temperatures between the transient models\nand our proxy evidence would require more work including in other regions of North America or\nEurasia and using other proxies. This is beyond the scope of this study which has a focus on water\nbalance and we prefer to avoid speculation at this point. GDD5 was only a first test to see if the\nproblems are overall to warm summers or only too warm peak summer temperatures in July.\nGDD5 gives a first indication that peak temperatures are the problem rather than the overall\nsummer warming.\nAs noted in a reply to Reviewer 1 (p. 11), we have significantly expanded the discussion\nof the proxy-data match achieved for T vs. for GDD5. We have also added the GDD5 values to\njul\nSupplementary Table 2.\nReviewer 2: The trend reversal in PC2 after 5ka really stands out (Figure 3, Figure 6).\nCan the authors clarify what drives PC2? Is it the extent of the LIS up until it is gone at ~5\nka at which point the main driver becomes the east-west moisture contrast? Or is it the\neast-west moisture contrast across the entirety of the record and PC2 aligns so closely\nwith LIS extent from 11 to 5 ka because the LIS extent controls the east-west moisture\ncontrast? (Line 452 – 454)\nAbout the PC2 interpretation, based on the spatial pattern in the loadings (Fig. 3c) the early-mid-\nHolocene fall in PC2 represents the increase in moisture contrast as the MW, with a lower\nbaseline moisture level, further dries while the NE starts to get wetter following the early-\nHolocene drought there. (This the CESM1 simulations can link to the impact of the LIS\nanticyclone, as we discuss in the manuscript.) The reversal in PC2 trend after 5 ka follows from\nthe stronger increase in water balance in MW compared to NE (seen in Fig. 6a), which now starts\nto decrease the midcontinent-Atlantic moisture contrast.\nin eastern North America” section of the Results and Discussion. We have clarified the wording\nin that paragraph.\nReviewer 2: It is somewhat surprising to the reader that the authors focus the second half\nof the discussion entirely on the time slice experiments from CESM1. I understand that\nthey focus on CESM1 because they are probing the dynamical relationship between\nhydroclimate and the LIS, but the authors could do a better job of outlining the rationale\nfor this focus in Lines 456 – 469. Additionally, it would be interesting to see the same\ndynamical treatment shown in Figures 7 and 8 given to time slices extracted from the\ntransient simulations. These analyses carried out for the transient simulations\n(individually or as an ensemble) for 8, 7, and/or 6 ka would be useful for further\nelucidating the influence of the LIS as it disappears. (Lines 490 – 537).\nThe CCA analysis is confusing and would benefit from a more fleshed-out description. It\nwould be helpful is the authors explicitly stated what the percentages in the Figure 8 sub-\nheaders mean in terms of the CCA.\nThe climate models are used to help with the proxy interpretation where two time periods stand\nout:\n1. The early Holocene shows a west-east dipole which has been hypothesized to be caused\nby the ice sheet. CESM1 is used to evaluate this and confirms the interpretation (wind\npattern in Fig. 7, CCA in Fig. 8).\n2. The remainder of the Holocene shows long-term transient changes which are now\nconsistent with newer climate model simulations. By merging results and discussions,\nwe can now make these points clearer by using subheadings.\nGiven the length of the paper and that the remainder of the Holocene is only characterized\nby consistent long-term trends between models and proxies, we do not see any added value for\nCCA in this context. The CCA was specifically used to test the role of ice sheets which are absent\nafter ~8 ka BP. The CCA without the role of ice sheets is included in Fig. 8 representing PI.\nWe have improved the description of the CCA in the revised “Results and Discussion”\nand in the “Methods”. The explained variances in Fig. 8 indicate how much of the overall variance\nis explained by the given CCA1 or CCA2 pattern (similar to PCA where the eigenvalue can be\nused to calculate the explained variances).\nReviewer 2. Methods\nWhat criteria do the authors use to subset the “high-resolution” fossil pollen records?\nThe criteria were based on the fact that in the Midwest and Northeast, there were subsets of sites\nwhich clearly stood out in terms of resolution of pollen data and number of dating (see\nSupplementary Table 1). The criteria (minimum of 125 pollen samples and 9 dating) are\nmentioned in Methods. The authors also had good understanding of these sites from prior work,\nand confidence that these sites have no outstanding issues related to e.g. taphonomy,\nsedimentation or chronology that would not be obvious based on a simple screening on Neotoma,\nwhich increased our confidence in using these thresholds.\nReviewer 2. The authors should outline which model grid cells are used to calculate mean\nvalues for each of the three regions. They could also mark the boundaries of these regions\n(as defined for the spatial extents in the model output) in Figure 1.\nWe have added the coordinates of the boxes in Methods. We now also show these boxes in Fig.\n1.\ncurves in Figure 6.\nAdded some detail in the caption; in general, the methodology is exactly identical with the main\nresults shown in Fig. 2.\nReviewer 2. Figures\nFigure 4: I’d recommend showing the points for 12,11, and 9 ka snapshots from CESM1,\nand potentially doing away with the lines connecting these points, to underscore that these\nare not transient simulations.\nWe agree and have split Fig. 4 in two panels (a = CESM1 equilibrium runs for 12, 11 and 9 ka;\nb = MPI-ESM and EC-Earth transient runs for 8–0 ka) to underline the fundamental difference\nin the modeling for the two time windows.\nReviewer 2. Why combine MPI-ESM and EC-Earth for this figure while keeping them\nseparate in Figure 2? Does this model ensemble offer something that the individual models\ndo not? The specifics of how this ensemble is calculated should be described in the text.\nFor example, were the models interpolated to the same spatial resolution prior to\naveraging? Would it be worthwhile to also include TraCE-21ka in Figure 4?\nIn Fig. 2, we consider it necessary to show the raw outputs of all the simulations we have used.\nHowever, given the results seen in Fig. 2, namely that the new-generation MPI-ESM and EC-\nEarth give very similar results, while showing much improved fit with the reconstructions\ncompared to TraCE-21ka, we consider it sensible to only show in Fig. 4 two models (MPI-ESM\nand EC-Earth) which are consistent with the reconstructions. Further, because the MPI-ESM and\nEC-Earth results are so similar, we prefer to simplify Fig. 4 by showing their mean instead of two\nnear-overlapping sets of curves, because we consider Fig. 4 somewhat challenging to\ncommunicate as-is.\nReviewer 2. What explains the big difference in sign of the precipitation anomaly between\nCESM1 at 9 ka and the mean of MPI-ESM1 and EC-Earth at 8 ka? Is this due to differences\nin the models (CESM1 vs the transient simulations) or interpreted as an actual feature of\nthe climate system?\nAs explained above in response to Reviewer 1 (p. 8), the jump is likely a combination of a real\nchange as indicated by the agreement of simulated and reconstructed water balances across\nmodels and the switch to different model experiments. As explained further up, the transient\nmodels start from equilibrium initial conditions at 8 ka BP and miss out on ocean memory of the\nsystem coming out of an ice age.\nReviewer 2. Can the authors include some metric of uncertainty in this figure, especially\nfor the average of MPI-ESM and EC-Earth.\nWe would prefer to leave Fig. 4 unchanged in this respect, as it might look odd to add a standard\ndeviation for the two models in Fig. 4b and have none for CESM1 in Fig. 4a (where there is only\none simulation). As the individual model time series are shown in Fig. 2 and have very similar\ntrends, the errors would add little value. It could also make Fig. 4b quite difficult to read as there\nare already 9 time series shown, and the idea being communicated is already quite complex.\nReviewer 2. Figure 7: The authors should consider graphically representing areas where\nthe difference between the time slice experiment (12, 11, or 9 ka) and the PI are not\nsignificant. This can be done using a Student’s T Test and shown on the figures by graying\nor whiting out the anomaly (usually around the transition from positive to negative).\nSST and hence surface temperature over the ocean is prescribed as climatology and hence has\nzero variance in time precluding the calculation of significance for temperature. As we do not\nmake any formal comparisons of models or proxies where significance testing would be relevant,\nwe would prefer to use simple anomaly maps to highlight the large-scale changes independently\nof local grid point significance over land.\nReviewer 2. Figure 8: I am confused as to what exactly these anomalies are relative to.\nCan the authors describe how the data in this figure are calculated? Are all these\nanomalies statistically significant / is there a way to test and communicate this? I\nunderstand that CCA1 doesn’t always have highest explained variance for P-PET or SLP\nindividually, but is it strange that there are times when the values for both P-PET and SLP\nare lower than those for CCA2?\nWe have now provided a more detailed explanation of CCA in the text. CCA anomalies help\nidentify specific deviations from expected patterns based on the relationships between variables,\nhere local water balance and large-scale SLP. For example, in CCA1, a water balance anomaly\nof +50 mm indicates that the large-scale SLP linked to the water balance causes the water balance\nto be 50 mm higher than it would otherwise be. The CCA patterns represent anomaly patterns\nfrom the basis used to calculate the EOFs. The interpretation of the according anomalies within\nthe patterns needs to take into account the respective values of the Canonical time coefficients.\nFor instance, a value of 1 in the time coefficients indicates a typical anomaly SLP-local water\nbalance pattern represented in the individual canonical CCA patterns. Higher/lower values in\nCCA time coefficients scale accordingly and can even lead to inverse anomaly patterns in the\npresence of negative CCA time coefficients\nReviewer 2: Lines\nLine 238: What is meant by “a poor spatial coverage of coastal areas in the model.”?\nThe model resolution even at ~1° is still too coarse and does not capture the land-sea distribution\nof coastal areas well. While proxies from lakes near the coast would still represent conditions\nover land, the model would average land temperatures with SST grid points making a model vs.\nproxy comparison inconsistent for any grid point that is not 100% land fraction. To avoid this\nproblem, we set all grids to missing values that include fractional values from SST along the\ncoastline.\nReviewer 2: Line 353: Add closing parenthesis. Add comma after “In contrast”\nDone.\nReviewer 2: Line 432: Can the authors briefly elaborate on the “internal variability” that\nmay be contributing to the ~200-year periodicities that are observed in both\nreconstructions and simulations?\nThere is still an ongoing debate about the nature and origin of low-frequency variability. Some\ncontrol simulations with static forcing of some models do show such low-frequency variability,\nhence purely internal variability, while others don’t. As shown for the EC-Earth model version\nalso used in our study here, internal feedbacks can sustain multi-centennial variability in AMOC\n(Cao et al. 2023). Externally forced transient simulations including the ones used here typically\nshow variability around 200-year periodicities, consistent with proxies (e.g. Askjær et al. 2022).\nIn that case, multi-centennial variability can be a mixture of forced (e.g. episodes with volcanic\neruptions, solar forcing) and unforced internal variations. It is noteworthy that EC-Earth, with\nonly transient forcing changes of radiative forcing (GHG, orbital), shows comparable ~200 year\nvariability as MPI-ESM that in addition includes variations in volcanic and solar forcing. We\nhave added the Cao et al. (2023) reference to the manuscript.\nDone.\nReviewer 2: Line 637: Reword the beginning of the sentence about CCSM3.\nDone.\nReviewer 3: This paper presents an impressive synthesis of pollen records and model\nanalyses, leveraging new regression tree methods to attempt to distinguish moisture\nbalance from temperature signals in the eastern US. The boosted regression tree approach\nappears to provide new insights not available via the common modern analog technique.\nThe authors then compare these results to transient climate model simulations to speculate\nabout the mechanisms responsible for past changes. Overall, I really like this paper, but I\nthink it could benefit from tightening the story and honing the ‘big picture’ message of the\nwork. I think these could be tackled via major revisions. My overall feeling is that the main\nstrength of this paper comes from the interesting new reconstruction. The model analyses\nare interesting, but also raise many questions and some confusion in those analyses\ndetracts from the overall strength of this work. I would not be afraid to restructure this\npaper as a ‘data driven’ paper.\nBelow are some overall points to consider as well as detailed comments on the text.\nFirst, I think the organization of the paper would benefit from a combined results and\ndiscussion section rather than having these sections distinct. As it stands, there are topics\nthat should traditionally be included in the discussion in the Results, and the organization\ngets a little bit messy. Especially for a short format paper, combining these would be\nvaluable.\nThank you for these helpful comments and for your positive feedback. As noted above, we’ve\nshifted to a combined Results and Discussion (with subheadings).\nWe note that the main motivation to conduct this study was twofold: (1)\nopportunistically, the availability of new proxy methods and new model simulations and (2) the\npending mismatch between climate simulations and proxy evidence regarding water balance\ntrends across the Holocene. We hence re-investigate the issue by employing a range of classical\nand new proxy methods which are compared with new climate model simulations to resolve the\nissue. As shown in Fig. 2a, the mismatch is now mainly attributable to the older model\nsimulations with TraCE-21k, which is resolved by the newer model simulations rather than the\nproxy reconstruction. We hope this is clearer in our revised version.\nReviewer 3: Second, the overall conclusion of the paper could be stronger. The abstract\nreads more like a literature review rather than the abstract of a short format paper with a\n‘punchy’ conclusion, and a similar comment could be applied to the paper’s discussion.\nMy personal thought is it might be better to focus either on the impact of long-term forcing\non the PCA/water balance reconstruction, or on the question of periodicities. As it stands,\nI don’t know that the model analysis is ‘enough’ to support the conclusion that summer\ninsolation and the Laurentide are driving the patterns you see. But, it is worth noting that\nTrace21ka has single forcing simulations (e.g. with ONLY orbital changes, or with ONLY\nice volume). : https://www.earthsystemgrid.org/project/trace.html\nIt is also worth taking a step back and starting with evaluating the PI model simulations’\nclimatology. Which ones are too rainy over your regions of interest? Too hot? Are the\nwinds too strong relative to observations? This might provide some additional insights\nabout why model performance varies. It would also be good to make a table as to what\n20\nforcing datasets were used in each set of transient runs - the change in different\ngenerations of ice sheet reconstructions (e.g. 5G, 6G, etc) will affect simulations even of\nthe same model.\nFinally, there are other ‘data’ products out there for deglacial climate variables (see\nOsman et al 2022 in Nature for the data-assimilation approach). My sense is that that\nproduct performs well for temperature but might not perform well for\nhydroclimate-related variables. It would be worth comparing your reconstructions (real world data! So more\nlikely to represent the truth) to the DA product for P and P/ET during the deglaciation.\nIt might help emphasize the innovative nature of the reconstructions that you’ve produced.\nAs mentioned above about the scope of the paper, our motivation for the paper was both the\navailability of a new robust pollen–water balance calibration as well as the increasing body of\nwell-dated, high-resolution fossil pollen sequences. This paper was then deliberately written to\napproach the datasets from the various angles allowed by the data properties. We acknowledge\nthat the resulting paper is somewhat complex and our decision to submit to Nature\nCommunications, with its fairly relaxed manuscript format, was deliberate.\nWe acknowledge the problem with our initial submission’s long and convoluted\nDiscussion spanning many topics. Our revised manuscript’s Results and Discussion is now\ndivided in thematic subsections, and also includes some new introductory and concluding\nsentences for the thematic sections, which we hope make our main conclusions stand out better.\nUnless there is a specific reason, it is generally not common to re-evaluate the\nperformance of already tested state-of-the-art climate models regarding their systematic biases.\nInstead, these systematic model biases are usually removed by, e.g., subtracting the PI\nsimulations from past (or future) simulations. This standard procedure is applied here so that the\nfocus is only on relative changes over time. It should also be noted that biases in the PI do not\nnecessarily say much about biases in, e.g. CESM1, for periods with fundamentally different ice\nsheets and land-sea distributions, as is the case for the early Holocene. We hence do not see merit\nin speculating about how present-day biases – which we removed – might or might not influence\ntransient changes in the past. The model evaluation of these changes is here based on the\nagreement with proxy data.\nWe agree with the reviewer that single forcing runs like available for TraCE-21k would\nbe an interesting addition to explore temperature and water balance changes further. This is\nhowever beyond the scope of this already long study. In addition, the fundamental mismatch of\nthe full-forcing TraCE-1k result with proxies and the new transient simulations in this study (Fig.\n2) would, a priori, raise questions about its validity.\nReviewer 3: Introduction\nLine 126-131: Instead of just reviewing the methods and approaches you used in the paper,\nprovide a signpost of the overall conclusion of the paper.\nWith the revised structure of the Results and Discussion, with six thematically titled subsections,\ncombined with the Abstract, we hope our main conclusions are now easy to identify in the paper.\nWe’d prefer to not reiterate the findings and the end of Introduction because they are, due to the\nnature of the paper, somewhat wordy to express.\nWe have, instead, edited the first paragraph of Introduction for brevity, including\nremoving mention of the MAT test which is now out of the main manuscript text. The paragraph\nshould now read more as a concise list of the manuscript’s goals, and lead logically to the\nimmediately following Results and Discussion subsections.\nReviewer 3: Results:\nFigure 2: Minor point, but I find it hard in some panels to distinguish the lines for EC-\nEarth and Trace 21ka. Make sure you check colorblindness compatibility of the figure,\nand perhaps use line weights/dash patterns to increase understandability.\nthe Bang Wong color blind friendly palette (https://davidmathlogic.com/colorblind). The palettes\nin the equilibrium model simulated fields in Figs. 7 and 8 are unchanged, as they already avoid\nusing red and use white to separate negative and positive values.\nReviewer 3: Line 183: This applies to the paper overall, but it is worth defining what you\nmean by ‘drought’ - for paleo records like this, we tend to define drought as just below\naverage P or P-E, while the term can mean something very different in the modern\nliterature (e.g. below a certain percentile of precipitation).\nWe define “drought” here as annual water balance being under the preindustrial value. This is\narticulated in the first sentence of our revised Results and Discussion.\nReviewer 3: Line 196: I think the PCA presents one of the most interesting results in this\npaper, and it is worth really strongly pinning down what the driver of this pattern is. I am\nnot sure I understand why PC2 is Laurentide deglaciation if it goes ‘up’ again after 5 ka?\nSee earlier reply about PC2 interpretation (p. 16).\nReviewer 3: Line 225: I am not sure I agree that there is a good comparison with the\nevolution of water balance between the models and the data. CESM1 seems to overlap the\nwater balance reconstructions, but seems to be very off in terms of GDD so this to me\nsuggests that it may be getting the ‘right’ answer in terms of hydroclimate change for the\nwrong reasons. I am therefore not sure how much to trust the CCA. The logic behind\nrelying on CESM1 for interpretation needs to be clearer. An alternate/complementary\nanalysis might involve comparing Trace to MPI and EC Earth for an interval like 7 ka,\nwhen the former is not dry enough and the latter two are too dry. Difference fields between\nthese runs might reveal why the models differ in skill at simulating deglacial climate\nchange. See comments above about some basic diagnostic analyses as well.\nBy design, the water balance may include compensating effects between PET (which is driven\nby summer temperatures and incorporating the season length above freezing) and annual\nprecipitation. The relatively good agreement of simulated vs. reconstructions, i.e. regarding the\nconsistency of trends, is a clear improvement to previous studies. The additional comparison of\nreconstructed vs. simulated GDD5, which characterizes the overall warm season instead of peak\nJuly temperatures (now Supplementary Fig. 1) indicates a better agreement of GDD5 between\nproxies and MPI-ESM and EC-Earth, suggesting that the problem is rather in temperature and\nhence PET than the water balance itself. For CESM1, the problem is more the seasonal GDD5\nrather than July, when comparing to proxies. So in both cases, we would rather trust the water\nbalance than the details of peak July temperature. Hence, we’re confident that the CCA for water\nbalance is reasonable given the proxy-model agreement for water balance. The problem is rather\nin temperature (seasonal vs. peak July disagreement).\nAs noted in a reply to Reviewer 1 (p. 11), we have significantly expanded the paragraph\ndiscussing the significance of the GDD5 comparison between proxies and models.\nReviewer 3: Figure 6: In panel d, if this mode reflects Laurentide deglaciation, why does\nit go back ‘up’ slightly after 5 ka?\nSee earlier reply about PC2 interpretation (p. 16).\nReviewer 3: Line 305: I am not sure if the section on periodicities is a good idea for this\npaper. It seems unrelated to the overall story of the paper, and could be a separate, more\nrobust analysis in and of itself.\nSee earlier reply about our rationale for the scope of the paper (p. 21).\nSuggest combing with results section, much of this discussion could be integrated through\nthe prior section\nAs noted above, we have done this.\nReviewer 3: Line 394: See comments above about model-data comparison.\nWe refer back to that reply (p. 22).\nReviewer 3: Conclusions:\nThe concluding paragraph of your paper focuses on the MPI, EC Earth, and new\nreconstructions. If this is the focus, it would be important to then focus on these model\nsimulations in Figures 7 and 8.\nWe have added subheadings in the Results and Discussion which should clarify the reason for\nthe shift from CESM1 to MPI-ESM and EC-Earth, as we no longer deal with the period affected\nby the ice sheet (which we study using the CESM1 runs).\nReviewer 3: I would also check about the length limits for this journal - it seems like there\nare a lot of figures, and the text is quite long for a short format journal.\nWe’re currently within the length limits, which for Nature Communications are in fact\nsomewhat relaxed (5000 words excluding Methods, and 10 display items).\nIn addition to the changes described above, we have made minor edits for grammar and clarity.\nAlso, we note that we plan to improve the quality of the Supplementary Figures 1–10 and 12–\n13 for the final revisions before publication. This will include a better resolution, and a taller\nvertical scaling of the reconstruction panels in Supplementary Figures 1–10, making them\neasier to read. We’ve had technical issues exporting from R the dotted and dash-dot lines which\nare used in these figures for the individual fossil sites. We’ll sort this out for the final revisions.\nFinally, we wish to again thank the reviewers once more for their thorough and\nthroughtful comments.\nRespectfully,\nHelsinki, 6 December 2024\nDr. J. Sakari Salonen, on behalf of the co-authors\nE-mail: sakari.salonen@helsinki.fi\nReferences\nAskjær TG, Zhang Q, Schenk F, Ljungqvist FC, Lu Z, Brierley CM, Hopcroft PO, Jungclaus J,\nShi X, Lohmann G, Sun W, Liu J, Braconnot P, Otto-Bliesner BL, Wu Z, Yin Q, Kang Y,\nYang H (2022) Multi-Centennial Holocene Climate Variability in Proxy Records and\nTransient Model Simulations. Quaternary Science Reviews 296:107801.\nhttps://doi.org/10.1016/j.quascirev.2022.107801\nCao N, Zhang Q, Power KE, Schenk F, Wyser K, Yang H (2023) The Role of Internal Feedbacks\nin Sustaining Multi-Centennial Variability of the Atlantic Meridional Overturning\nCirculation Revealed by Ec-Earth3-Lr Simulations. Earth and Planetary Science Letters\n621:118372. https://doi.org/10.1016/j.epsl.2023.118372\nOverprinted by Ocean Temperature Fluctuations. Nature Geoscience 15:899-905.\nhttps://doi.org/10.1038/s41561-022-01056-4\nKuang X, Schenk F, Smittenberg R, Hällberg P, Zhang Q (2021) Seasonal Evolution Differences\nof East Asian Summer Monsoon Precipitation between Bølling-Allerød and Younger Dryas\nPeriods. Climatic Change 165:19. https://doi.org/10.1007/s10584-021-03025-z\nSalonen JS, Korpela M, Williams JW, Luoto M (2019) Machine-Learning Based Reconstructions\nof Primary and Secondary Climate Variables from North American and European Fossil\nPollen Data. Scientific Reports 9:15805. https://doi.org/10.1038/s41598-019-52293-4\nSchenk F, Väliranta M, Muschitiello F, Tarasov L, Heikkilä M, Björck S, Brandefelt J, Johansson\nAV, Näslund J-O, Wohlfarth B (2018) Warm Summers During the Younger Dryas Cold\nReversal. Nature Communications 9:1634. https://doi.org/10.1038/s41467-018-04071-5\nLATITUDE EASTERN NORTH AMERICA\nResponse to Referees\nWe quote below the reviewers’ full comments and describe the changes done to address each\npoint. We thank all the reviewers for their helpful comments.\nReviewer #1: Dear Editor,\nI am generally satisfied with the authors’ responses and have no major reservations\nregarding the manuscript. They have addressed my comments in many cases and have\ncarefully justified their choices in others. I appreciate their revisions, particularly the\nimproved methodology for stacking and the clarifications of other methodological aspects.\nOne aspect that still leaves me somewhat uncertain is the analysis of periodicities. I believe\nthe current framing places too much emphasis on the least informative aspect—wiggles in\nthe wavelet spectra that, after appropriate bandpass filtering, are interpreted as significant\nperiodicities. That said, I appreciate that the authors have tested the sensitivity of their\nperiodicity estimates by varying the null hypothesis, and the analysis is still informative of\nthe distribution of variabiltiy across frequencies.\nIf there are no restrictions on adding new supplementary figures, I would still recommend\nincluding the mean wavelet plots used for the periodicity estimates, as these would provide\nvaluable insights. Specifically, plotting the individual wavelet spectra for all sites within a\ncluster together would allow the reader to see the \"periodicities\" clustering and,\nimportantly, the overall shape of the wavelet spectra. Alternatively, the mean wavelet\nspectra could be shown for each region (i.e. a stack of the individual spectra, averaging in\nthe spectral rather than the time domain) with shading to indicate the spread among\nindividual records.\nOverall, I see no reason to further delay the publication of the manuscript.\nBest regards\nWe have added the mean power plots of all wavelet spectra (shown in Supplementary Figures\n18–25) in a new Supplementary Figure 26. This figure is cited, like Supplementary Figures 18–\n25, in the caption to Figure 6 which presents the wavelet results in the main manuscript.\nReviewer #2: The authors have submitted a response letter and updated version of the\nmanuscript that effectively addresses the comments raised by each of the three reviewers.\nThe reorganization of the Results/Discussion section in particular has strengthened the\nmanuscript. I recommend that the updated version of the manuscript be published and\nthank the authors for their detailed work to address reviewer comments.\nWe thank the reviewer once more for their helpful comments on this manuscript.\nReviewer #3: The authors have largely addressed my comments and those of the other\nreviewers. I could use some extra clarification on this text near line 428:\n\"This spatiotemporal pattern has been hypothesized to be explained by LIS deglaciation,\nand a following waning of the glacial anticyclone and an increased influence of the\nBermuda subtropical high, leading to a rerouting of the northward moisture advection\ntowards the NE region and away from the mid-continent\"\nI'd like a bit more detail - this sentence seems to suggest that after the Laurentide is gone,\n1\nthe upwards trend of PC2 post 5ka could reflect shifts in the strength of the North Atlantic\nSubtropical High? it would be nice to propose that explicitly in the text.\nThe rotational component of the western part of the Bermuda high generally helps to advect air\nand moisture from the (sub-)tropics northward over the US. This northward advection would be\nreduced while the LIS reaches far south. The role of the Bermuda subtropical high was mentioned\nas part of the original hypothesis in the literature but doesn’t imply a change in the high itself. It\nis rather the retreat of the ice sheets which allow the flow and moisture transport to move further\nnorth and by doing so increases the influence of the Bermuda high.\nWith the results presented in this paper, we have no means to infer potential changes in the\nBermuda high itself in terms of proxy data, and so we find it difficult to address this issue in the\ntext. The CESM1 simulations do not show any change in the wind or SLP pattern of the Bermuda\nhigh itself for the period 12, 11 and 9 ka (Fig. 7a,d,g).\nReviewer #3: Given the revised structure of the paper, it would be nice to have a conclusion\nsection that summarizes the overall argument of the paper and a few thoughts on future\ndirections. The last paragraph could be slightly expanded and serve as your Conclusion\nsection.\nWe have edited the last subsection of Discussion to start by recapping the fact that the early\nHolocene hydroclimate features showed a significant impact of the Laurentide Ice Sheet – and\nthen transition to discussing the temperature and evapotranspiration dominated regime of the mid-\nlate Holocene. With this, we feel the last subsection of Discussion adequately summarizes the\nmain features of the entire Holocene, while concluding with a forward-looking perspective.\nWe’re open to adding a separate Conclusions section, perhaps as bullet points, if the editor\nfeels this would be helpful. But on balance we’d favour going without, as we’re already at the\nword limit.\nReviewer #3: For the abstract, while it is true that the midHolocene features evaporatively-\ndriven drying similar to the future, it does not necessarily follow that they are analogs,\nwhich to me implies a symmetry of underlying mechanisms (and the radiative forcing of\nCO2 is different in character than seasonal changes in insoloation). Rather, it seems more\nfair to say that similar processes are at play. I would revise the abstract accordingly.\nAgreed; we have edited the Abstract to state that “...the mid-Holocene multi-millennial drought\nwas driven by similar processes compared to the drying trajectory projected…”\nIn addition to the changes described above, we have done the following changes\n● Fixed two cases of “analogue” as “analog”.\n● We have replotted Supplementary Figures 1-10 (paleoclimate reconstructions including\nresults for individual fossil sites) and Supplementary Figures 12 and 13 (diagnostic plots\nof most important predictor taxa in the calibration models) to make them more readable.\nThis does not change the data shown in the figures. The vertical scaling of the\nreconstruction panels in Supplementary Figures 1–10 has been increased, and\nSupplementary Figures 12 and 13 now use correct line symbology for the fossil sites.\n● We have added a Data Availability statement in the manuscript. We have uploaded the\nSupplementary Data file on Figshare (with a reserved DOI but not yet published).\nFinally, we wish to once more warmly thank the reviewers for their thoughtful comments during\nthis review process, which have greatly helped us improve the manuscript.\nRespectfully,\nHelsinki, 25 February 2025\nDr. J. Sakari Salonen, on behalf of the co-authors\nE-mail: sakari.salonen@helsinki.fi\n2"}]

Decision letter