Changes in Atmospheric Dynamics Over Dansgaard‐Oeschger Climate Oscillations Around 40 ka and Their Impact on Europe

Dansgaard‐Oeschger (D‐O) climate variability during the last glaciation was first evidenced in ice cores and marine sediments, and is also recorded in various terrestrial paleoclimate archives in Europe. The relative synchronicity across Greenland, the North Atlantic and Europe implies a tight and fast coupling between those regions, most probably effectuated by an atmospheric transmission mechanism. In this study, we investigated the atmospheric changes during Greenland interstadial (GI) and stadial (GS) phases based on regional climate model simulations using two specific periods, GI‐10 and GS‐9 both around 40 ka, as boundary conditions. Our simulations accurately capture the changes in temperature and precipitation as reconstructed by the available proxy data. Moreover, the simulations depict an intensified and southward shifted eddy‐driven jet during the stadial period. Ultimately, this affects the near‐surface circulation toward more southwesterly and cyclonic flow in western Europe during the stadial period, explaining much of the seasonal climate variability recorded by the proxy data, including oxygen isotopes, at the considered proxy sites.

Climate model studies of glacial climate states in Europe often focus on the last glacial maximum (LGM) (Kageyama et al., 2018(Kageyama et al., , 2021)).For the LGM, the boundary conditions are well constrained and many model studies have been conducted to disentangle the effects of the individual boundary conditions (e.g., orbital parameters, greenhouse gas concentrations, ice sheet extent) on atmospheric dynamics and thereby also on the glacial climate and weather patterns (e.g., Hofer et al., 2012;Ludwig et al., 2017;Shi et al., 2023).While the orbital configuration and the low greenhouse gas concentrations during the LGM result mainly in colder temperatures, the presence of large ice sheet.alters the location of the eddy-driven mid-latitude jet over the North Atlantic (Cao et al., 2019;Hofer et al., 2012;Pausata et al., 2011).Under interglacial conditions, three preferred modes of the jet position exist, while glacial conditions suppress this trimodality and favor a southward shift of a generally stronger and more zonally oriented jet instead (Löfverström et al., 2014;Merz et al., 2015;Woollings et al., 2010).
When addressing stadial and interstadial phases of MIS 3, climate model studies are less abundant.Many of them focus on the conditions of the North Atlantic ocean and variations in ocean circulation, given the importance of the AMOC in driving the transitions from stadials to interstadials (e.g., Pedro et al., 2022;Van Meerbeeck et al., 2011).Because D-O events occur relatively synchronously across the North Atlantic, Greenland and Europe, also an atmospheric transmission mechanism is hypothesized (Czymzik et al., 2020;Luetscher et al., 2015;Újvári et al., 2017, 2021).This role of atmospheric regime shifts was recently proposed for the explanation of abrupt climate changes occurred during the last deglaciation using transient simulations performed by a simplified Earth system model (Planet Simulator version 16) (Andres & Tarasov, 2019).It has been demonstrated recently that unforced, fully coupled global climate models are able to reproduce D-O-like climate oscillations (Armstrong et al., 2023;Fohlmeister et al., 2023;Kleppin et al., 2015;Klockmann et al., 2020;Peltier & Vettoretti, 2014), in response to gradual changes in the height of the Northern Hemisphere ice sheets (Zhang et al., 2014), CO 2 forcing (Vettoretti et al., 2022;Zhang et al., 2017) or salinity fluctuations (Armstrong et al., 2023;Peltier & Vettoretti, 2014) and explained by the internal feedbacks of the atmosphere-sea ice-ocean system (Li & Born, 2019).These models mainly investigate the underlying mechanisms of D-O variability, but are less concerned with the aforementioned reorganizations in atmospheric circulation and transmission mechanisms that may be responsible for the changes in the continental climate of Europe during stadial/interstadial periods of the last glaciation.
In the present paper, our main focus is on atmospheric changes over Greenland interstadial and stadial phases occurred in the North Atlantic and Europe for the period of 41.5 to 38.5 ka based on regional climate model simulations using GI-10 and GS-9/HE4 boundary conditions.We address the following three research questions: 1. How well do the paleoclimate simulations of GI-10 and GS-9 match the available proxy data? 2. What are the main differences between the interstadial and stadial large-scale atmospheric dynamics and regional circulations?3. What are the implications of these differences (if any) for proxy data interpretation?
For this purpose, we present a quantitative model-proxy data comparison in Section 3.1 and the simulated climate in Section 3.2.In Section 3.3, the simulations are analyzed by focusing on atmospheric circulation changes between the interstadial and stadial periods.The results are then discussed in the light of the available proxy data and concluded in Section 4.

Proxy Data
Stadials and interstadials of the Dansgaard-Oeschger cycles of MIS 3 are documented in various terrestrial paleoclimate archives including speleothems, lake sediments and loess deposits across Europe and western Asia (Voelker, 2002).However, the temporal coverage between 41.5 and 38.5 ka is limited and proxy data for temperature or precipitation is scarce, hampering quantitative paleodata-model comparisons for the GI-10 and GS-9 periods within our model domain.Speleothems provide the most precise chronologies, but their stable isotope signals (δ 13 C and δ 18 O) can only be qualitatively compared with simulated vegetation or temperature changes (Baker et al., 2011;Lechleitner et al., 2018).Of these, two prominent speleothem proxy records were chosen: (a) the Villars cave from France (Genty et al., 2003(Genty et al., , 2005)), and (b) the Sofular cave record from Turkey (Fleitmann et al., 2009).Both records show a rather muted δ 18 O signal and thus the variations in the δ 13 C records were emphasized and interpreted in relation to D-O variability in the original publications.These records were therefore used to compare the trends in simulated paleoclimate variations for the GI-10 and GS-9 periods with the regional vegetation/temperature changes detected in these speleothems.
Pollen grains, insects and especially chironomids preserved in lake sediments allow for the quantitative reconstruction of seasonal temperatures (Brooks & Birks, 2001;Heiri et al., 2011).However, only few such records Oxygen isotope and Ca 2+ data, the latter representing atmospheric dust variations, are from Seierstad et al. (2014).The paleotemperature reconstruction is from Kindler et al. (2014) and based on δ 15 N data combined with a firn densification and heat diffusion model.Greenland Stadial (GS) and Greenland Interstadial (GI) phases are indicated following Rasmussen et al. (2014).HE4 stands for Heinrich Event 4 (onset and termination after Guillevic et al. (2014)).All proxies are shown on the GICC05 timescale.
Loess deposits are the third type of archives allowing reconstructions of seasonal temperatures and mean annual precipitation on land for stadial/interstadial timescales using the stable carbon and oxygen isotope compositions of earthworm calcite granules (Prud'homme et al., 2022), and specific species of land snails (Radaković et al., 2023;Sümegi & Krolopp, 2002).Only one relatively well-dated loess record could be involved in our study, which provides quantitative paleoclimate information on GS-9: the Schwalbenberg loess sequence in Germany (Prud'homme et al., 2022).
The above-mentioned palaeoclimate archives are found between longitudes 0-32°E and latitudes 40-51.5°N in Europe and western Asia (see Figure 2), between 90 and 1,208 m above sea level, but predominantly below 500 m (Data Set S1).For the paleodata-model data comparison, we used the simulated data of the grid point in which the given proxy site is located.

Climate Model Simulations
Regional climate simulations are performed by nesting the Weather Research and Forecasting (WRF) model (Skamarock et al., 2021) into global runs of COSMOS (Klein et al., 2023;Zhang et al., 2013).These global runs have been chosen as the available data could be used as input for downscaling to regional scales.Several studies have demonstrated that increasing the model resolution through downscaling results in a more accurate representation of physical processes and orography (e.g., Giorgi, 2019;Ludwig et al., 2017Ludwig et al., , 2019)).Consequently, this approach provides added value for the representation of regional climate (Feser et al., 2011;Jacob et al., 2014) and has the potential to provide improved climate data for the comparison with proxy data.
In this study, the horizontal resolution across the model domain (including North America, the North Atlantic Ocean and Europe, see Figure 2) is refined to 1° (with 41 vertical layers) in both latitude and longitude compared to the coarser COSMOS data.WRF simulations for PI, GI and GS conditions are conducted for a period of 35 years each, where the first 5 years are reserved for model spin-up and the remaining 30 years are used for analysis.The WRF model uses 6-hourly boundary data from COSMOS simulations for PI to be compared with the current climate, as well as GI and GS boundary data for glacial conditions.The time step of the integration of the WRF model is set to 360 s and the output frequency is set to 6 hr.The global COSMOS model runs at T31 (triangular truncation at wave number 31) spatial resolution for the atmosphere (∼3.75° horizontal resolution) with 19 vertical layers and is coupled to an ocean model, including a thermodynamic sea ice model, employed on a bipolar curvilinear model grid with a formal horizontal resolution of 3.0° × 1.8° and 40 vertical layers.
Glacial boundary conditions in COSMOS are set to those at 40 ka, as described in Zhang et al. (2021).Specifically, trace gas concentrations are set to 193 ppm for CO 2 , 400 ppb for CH 4 and 240 ppb for NO 2 , while orbital parameters are set to 0.013146 for eccentricity, 23.6126 for obliquity and 358.167 for the longitude of perihelion with respect to the vernal equinox (Laskar et al., 2004).To realistically mimic GS conditions, the thermohaline circulation in the North Atlantic is weakened by freshwater release of 0.15 Sv to the Ruddiman Belt (40-55°N; 45-20°W), resulting in attenuated SSTs over the North Atlantic (Knorr et al., 2021).
Several modifications have successfully been implemented in the WRF model (Version 4.3), enabling research centered around glacial boundary conditions (Ludwig et al., 2017(Ludwig et al., , 2018;;Pinto & Ludwig, 2020;Schaffernicht et al., 2020;Újvári et al., 2022).In this study, the orbital configuration and greenhouse gas concentrations in the WRF model are the same as those used by COSMOS.The GI and GS simulations depict varying land ice coverage, ice sheet heights, and coastlines (see Figure 2) obtained from reconstructions by the ICE-6G-C (Peltier et al., 2015) available from the PMIP4 database (PMIP4, 2019).Based on various reconstructions, the sea level during the considered period (GI-11 to GS-9) varies between −55 m and −60 m (relative to today's sea level) for interstadial and −75 m to −90 m for stadial conditions (Clark et al., 2009;Lambeck et al., 2002;Yokoyama et al., 2007).
As there are no reconstructions available for periods of our interest, but to account for these variations, we utilized the ICE-6G-C reconstructions of ice sheet cover, topography, and coastlines at 17 ka for the regional GS simulation and the ICE-6G-C reconstructions at 14 ka for the regional GI simulation as best analogs, respectively.In the absence of comprehensive gridded data sets, we utilized the potential vegetation cover as outlined in (Ramankutty & Foley, 2010) for all regional simulations.

Evaluation Criteria of Model Simulation Data
First, an evaluation of the model data is compiled, comparing the PI simulation with ERA5 reanalysis data (Hersbach et al. (2020); see Text S1 in Supporting Information S1) and the paleoclimate simulations with the available proxy data.For the proxy-model comparison, we used the simulated data from the grid point where the corresponding proxy site is located and compare the respective long-term means with the standard deviation obtained from the 30 simulated years.
An overview of the simulated climate under PI, GI, and GS conditions is given with long-term annual, winter, and summer means of 2 m air temperature and daily precipitation distributions.The atmospheric dynamics and their changes between the stadial and interstadial phases are assessed via an in-depth analysis of the jet position and storm tracks.For the jet, the long-term annual and seasonal means of the wind speed in 200 hPa are taken into account, as well as a vertical cross-section of the zonal wind, zonally averaged over the North Atlantic region (60°W-0°W).The latitudinal jet position is then objectively identified following the method after Woollings et al. (2010).We interpolate daily data of the zonal wind in 850 hPa on a regular 1° × 1° grid and take a zonal average over the same North Atlantic region as for the vertical cross section (60°W-0°W).We apply a 10-day low-pass filter to remove synoptic systems such as individual storms crossing the North Atlantic.To finally determine the jet position, we use these data to locate the maximum wind speed between latitudes of 20°N and 75°N.
Storm track activity is a measure of synoptic activity and is identified using an Eulerian perspective by applying a 2.5-6 days band pass filter to 6-hourly mean sea level pressure data as described in Hoskins and Valdes (1990), Ulbrich et al. (2008) and Ludwig et al. (2016).
The regional atmospheric flow is analyzed based on the circulation weather type (CWT) approach after Jones et al. (1993).It has been applied for glacial conditions in for example, Ludwig et al. (2016) and Schaffernicht et al. (2020).With this method, the flow of each day is categorized into 10 CWTs, including eight directional weather types (north, northeast, east, southeast, south, southwest, west, northwest) and two rotational weather types (cyclonic and anticyclonic).For the classification, daily geopotential height in 850 hPa is interpolated on a regular 2.5° × 2.5° grid.A total of 16 grid points around a central point are considered to characterize the atmospheric flow at the central point.We classified the circulation for two locations in close proximity to the available proxy sites (see Figure 2): 50°N, 7.5°E near Schwalbenberg and 45°N, 0°E near Villars Cave.In our study, we combined all easterly circulation types (northeast, east, southeast) into a single category due to their infrequent occurrence.

Quantitative Model-Proxy Data Comparison
In this section, a quantitative comparison between proxy data and model data is presented.A direct comparison is performed in case the proxy reconstructions reveal absolute temperatures or precipitation values for either the GI-10 or the GS-9 phase.We consider a total of 17 data points for the different types of temperatures including mean annual air temperature, January temperature, July and summer temperature, and summer soil temperature.Ten values are related to the GS-9 phase, whilst seven values are associated with the GI-10 period.For precipitation, only two data points could be identified for each period.
The comparison of temperatures is shown in Figure 3.The magnitude of the different temperature values and their relative changes between the stadial and interstadial phases are well represented in our model.The Pearson correlation coefficient for all data points is p = 0.91, with a significance of 8.55.When considering only the GI or GS phase, the coefficient is p = 0.92 but with a slightly reduced significance (5.37 for GI, 6.88 for GS).The root mean square error (RMSE) amounts to 4.4°C (both periods), 3.2°C (GI) and 5.1°C (GS).These are reasonable values even for simulations at higher resolutions (e.g., Kageyama et al., 2021;Varga & Breuer, 2020).The model overestimates the coldest winter temperatures and underestimates the warmest summer temperatures.Nevertheless, soil temperatures are still within the estimated range of uncertainty of the proxy data.This leads to a lower seasonality in the simulations.
Comparison of the simulated precipitation with proxy data remains challenging due to the limited data availability.Applying robust statistics is not possible due to the scarcity of data points.The simulations are biased toward excessively wet conditions with an RMSE of 695 mm/year.This roughly corresponds with the amount that the proxy reconstructions reveal and is a known issue for glacial climate model simulations (Ludwig et al., 2017).Nevertheless, the transition from the warmer and wetter interstadial to the colder and drier stadial is captured, which is why the Pearson correlation coefficient (based on the four data points) is still good, with p = 0.95 and a significance of 4.22.

Simulated Present-Day and Glacial Climates in Europe
The simulated climate in Europe under PI conditions show largest seasonality in northeastern Europe, featuring cold winter and moderate summer temperatures (cf. Figure 4).In terms of precipitation amounts (cf. Figure 5) the highest levels are simulated over the North Atlantic and on the northeastern Iberian Peninsula during winter.Summer precipitation occurs mostly over land with higher values over the Alps, northern Great Britain and Scandinavia.
Under both stadial and interstadial conditions, the climate is much colder (by up to 8°C over land) compared to the pre-industrial (PI) era.The largest differences occur over the ice sheet itself, especially during winter.Despite the presence of ice sheets under glacial conditions, summer temperatures remain closer to PI values.Proxy-based temperature values are plotted against model temperature values of the associated grid point, both including the respective error bars: that denote one standard deviation derived from the original publication for proxy data or from model data, respectively.The shape of the data points denote the proxy type: circles represent data from lake sediments, squares represent data from peat bogs, and rotated squares represent data from loess.Colors denote temperature types: blue stands for January temperatures, orange for mean annual temperatures, different shades of red denote summer temperatures, July temperatures, and summer soil temperatures.Symbols with black outlines refer to the GI-10 period, symbols without outline refer to the GS-9 period.The black diagonal line denotes perfect model-proxy agreement.
Consistent with proxy data, the climate during stadial conditions is generally colder than during interstadial conditions.These differences are particularly large over the northern North Atlantic, the ice sheet and at the ice sheet edges primarily due to variances in ice sheet extent.In contrast, summer temperatures in western Europe remain similar in both phases.
The simulated precipitation amount is lower during glacial periods, particularly over land.Over the North Atlantic, the total precipitation is locally enhanced in both phases, but especially south of Iceland under interstadial conditions in winter.Europe is drier in summer, aside from the Black Sea, where rainfall is highest during stadial conditions.As the region south of Iceland receives particularly high amounts of rainfall during the interstadial in winter, the differences between the two glacial climates are most pronounced in this region.Both in winter and annually, the precipitation anomalies are divided in a north-south direction.While northern Europe and the northern North Atlantic are considerably drier, the southern North Atlantic and the Mediterranean region are wetter.In summer, precipitation over land is overall reduced.

Atmospheric Dynamics Under Present-Day and Glacial Conditions
To assess the large-scale atmospheric dynamics, we consider the jet position and storm tracks.The Circulation Weather Types for two central points near the available proxy sites in Europe (namely Schwalbenberg and Villars Cave) then bridge the gap to the regional circulation changes in Europe.
In winter, the wind speed at 200 hPa is zonally oriented under PI conditions, causing the winter jet to approach Europe at a latitude of the Iberian Peninsula around 45°N (see Figure 6).During summer, the jet weakens and a northward tilt is developed, resulting in a wider cross-section of the zonally averaged zonal wind field (see  d-f), GS conditions (g-i), and 2m air temperature differences between the GS and GI simulation (j-l).Thick black lines denote the respective coastlines, pink lines mark the ice sheet extent (GS in the difference fields).
Figure 7).This pattern shifts during glacial conditions.While the general behavior of a more tilted summer jet remains, a wave-like pattern develops.This pattern is even more pronounced during the colder stadial phase, leading to stronger wind speeds over the northern North Atlantic and weaker winds along the southern European coastlines.In winter, the jet is stronger and tilted northward under glacial conditions, with higher maximum wind speeds in the stadial phase.The change pattern between stadial and interstadial is reversed as the jet shifts southwards (see Figure 7k), resulting in stronger winds in southern Europe and weaker winds south of Iceland.Again, the differences between the two glacial periods are largest in winter.
Histograms depicting the occurrence frequency of daily jet latitude demonstrate the variability of the jet position over time as well as its changes under different climate conditions (see Figure 8).During all periods and seasons, there is a clear primary jet location around 45°N-most pronounced by high occurrence frequencies during the interstadial phase.An additional maximum is detected around 55°N for both glacial jets.This maximum is more pronounced during the interstadial in winter and in annual mean, whilst the southerly positions of the jet occur more often during the stadial phase.The more distinct wave-like pattern in the summer jet during the stadial phase results in more southern and more northern jet positions compared to the interstadial phase, while the peak located at 45°N is reduced and shifted southward.
Storms in the North Atlantic typically follow the northern edge of the jet stream and the respective sea ice boundaries (cf. Figure 9).Winter storm tracks are more intense and penetrate further into the continent, although the most intense parts remain over the ocean under PI conditions.Under glacial conditions, storms are deflected by the Fennoscandian Ice Sheet and a second maximum evolves downstream, resulting in higher storm track intensities over land in eastern Europe and western Asia.These second maxima are also visible in summer and annual mean, both with much lower intensities.d-f), GS conditions (g-i), and total precipitation differences between the GS and GI simulation (j-l).Thick black lines denote the respective coastlines, pink lines mark the ice sheet extent (GS in the difference fields).
10.1029/2023JD040247 9 of 16 In annual mean and in winter, the GS storm track is extended southward compared the GI phase.Regions where sea ice extends further south are also characterized by stronger storm tracks.This is particularly the case south of Iceland.However, storm track maxima are less intense during the stadial phase.
The prevailing circulation types and their changes are derived with a CWT analysis.Figure 10 illustrates the occurrence frequency in annual, winter and summer means during the different periods, alongside changes from the interstadial to the stadial phase.The westerly circulation type dominates winter flow, while south-and northwesterly flows also occur frequently, although less often than purely westerly circulation.The atmospheric circulation during PI period indicates more cyclonic and southwesterly flows compared to the glacial periods.The circulation close to the North Atlantic coast (Villars Cave) is more influenced by westerly wind conditions than Schwalbenberg.Here, the cyclonic and northwesterly circulations occur more often.In both locations, the anticyclonic and northwesterly flows are reduced during GS, while the cyclonic CWT and easterly (only Schwalbenberg) are enhanced.
In summer, the daily flow types exhibit a broader distribution, with a greater impact from anticyclonic circulation (Villars Cave) and more northerly and northwesterly flows (Schwalbenberg).Under stadial conditions, the westerly flow is suppressed, whereas the anticyclonic flow is intensified in both locations.This is in line with the wave-like jet structure observed for summer during both glacial phases, but more distinct during the stadial phase, when ice sheets were higher and extended further south.
The annual distribution of CWTs encompasses both winter and summer seasons, with prevailing westerly circulation types, and for Villars Cave, a frequent occurrence of anticyclonic circulation.

Discussion and Conclusions
The relative synchronicity of D-O events across the North Atlantic, Greenland and Europe implies a tightly coupled climate between these regions during the last glaciation.Obviously, the fast paleoenvironmental response Figure 6.Distribution of wind speed in 200 hPa in annual mean (left column), winter mean (middle column), and summer mean (right column) as simulated under PI conditions (a-c), GI conditions (d-f), GS conditions (g-i), and wind speed differences between the GS and GI simulation (j-l).Thick black lines denote the respective coastlines, pink lines mark the ice sheet extent (GS in the difference fields).
detected in these terrestrial archives suggests the existence of an atmospheric transmission mechanism of D-O variability.The most plausible hypothesis in this respect is that the stadial/interstadial phases are linked to a southward/northward shift of the eddy-driven polar jet stream (Czymzik et al., 2020;Luetscher et al., 2015;Újvári et al., 2017, 2021).Our timeslice experiments aimed to investigate this supposed connection by modeling atmospheric dynamics over Europe and the North Atlantic during Dansgaard-Oeschger oscillations around 40 ka.We conduct regional climate simulations with the WRF model forced by global data from the COSMOS model and simulate the PI as well as GI-9 and GS-10/HE4 periods for 35 years.
First, we evaluated our model performance using proxy data from mainland Europe.We found a very good agreement between the modeled temperatures and the proxy temperature data (Figure 3).The correlation is almost ideal (the Pearson correlation coefficient exceeds p ≃ 0.9) and the temperature estimation from the simulations is nearly unbiased with respect to proxy data.Remaining differences can mostly be attributed to orography differences.Very cold winter temperatures are often not well captured in model simulations (e.g., van Ulden & van Oldenborgh, 2006).The slightly underestimated summer temperatures resulted in smaller seasonal variations in simulations compared to those reconstructed from proxy data.Simulated precipitation shows larger differences from proxy records with values twice as large as reconstructed.Corresponding discrepancies have also been identified in other modeling studies (Ludwig et al., 2017;van Ulden & van Oldenborgh, 2006).Besides this large bias, the difference between stadial and interstadial precipitation amounts is well captured.We can therefore rely on the simulated variations between the two phases.
In general our simulations confirm that the European mainland climate was colder and drier during both stadial and interstadial periods with a more intense jet stream in comparison to PI (e.g., Merz et al., 2015;Sanchez-Goni & Harrison, 2010;Voelker, 2002).Similarly, but to a lesser extent, this holds also for the variations between the stadial and interstadial simulations.Here, the simulated GS climate is colder throughout the year, particularly in Figure 8. Kernel estimation of the PDFs of daily jet stream latitude (averaged over 60°W-0°W) for the whole year (a), for winter (b) and summer(c) as simulated under PI conditions (black lines), GI conditions (orange lines), and under GS conditions (blue lines), as well as GS minus GI differences (dashed gray lines).
Figure 9. Distribution of storm tracks in annual mean (left column), winter mean (middle column), and summer mean (right column) as simulated under PI conditions (a-c), GI conditions (d-f), GS conditions (g-i), and storm track differences between the GS and GI simulation (j-l).Thick black lines denote the respective coastlines, pink lines mark the ice sheet extent (GS in the difference fields), and cyan lines mark the 10% sea ice extent (GI dashed and GS solid in the difference fields).
winter.While temperature differences during summer are mostly negligible over (south-)western Europe, they reach −1.0 to −4.0°C in eastern Europe, British Isles and along the edge of the Fennoscandian Ice Sheet (see Figure 4).In line with the Clausius-Clapeyron relationship between temperature and atmospheric water holding capacity, precipitation decreases during the colder stadial phase, particularly in regions with distinct temperature differences.However, precipitation differences show strong spatial differences.The region south of Iceland shows the largest differences toward a considerably drier stadial winter in addition to a north-south aligned pattern with drier conditions north and wetter conditions south (Figure 5).This pattern is more pronounced in the winter season, thus being the dominant season in terms of temperature and precipitation changes for the simulated stadial and interstadial periods.
Next, we investigated the large-scale atmospheric dynamics and regional circulation in order to explain the observed climatic differences.Under stadial conditions, the eddy-driven jet stream intensified and its mean position is broadened and shifted southward during winter with a tilt toward the northeast (Figure 6).When considering the temporal variability of the jet position, it is evident that both stadial and interstadial periods show a clear preferred position around 44°N.However, during GS-9 the southern (40°N) positions are more frequent, while northern positions (50-55°N) are less prevalent (see Figure 8), leading to the general broadening of the jet stream distributions.Overall, the simulated storm tracks under glacial conditions are more intense and penetrate further into the continent in the winter season (Figure 9).During glacial periods, a second (annual and seasonal) maximum exists downstream of the Fennoscandian Ice Sheet, resulting in higher storm track activity over land in eastern Europe and western Asia.However, there are only minor differences in the annual mean and winter storm track differences between the GS and GI periods, which mainly occur south of Iceland and over the western Mediterranean.
In accordance with the shift and tilt of the jet stream during the stadial period, the near surface circulation exhibits fewer annual westerly and northwesterly flow types at Villars Cave and the Schwalbenberg loess site (see Figure 10).During winter, less northwesterly/anticyclonic and more southwesterly/cyclonic flow types are found at both proxy locations during the stadial period.This is in agreement with the colder/drier conditions indicated by higher δ 13 C/δ 18 O values in the Villars Cave calcite (Genty et al., 2003).These heavier δ 18 O compositions, primarily explained by temperature variations, can at least partly be accounted for by different, more southerly moisture sources as indicated by our CWT analysis.The stadial layers of the Schwalbenberg loess record show more negative δ 18 O of earthworm biospheroids, which is a summer season proxy (Prud'homme et al., 2022).The more negative values may dominantly be related to lower temperatures, but a moisture source effect also cannot be excluded considering more frequent appearance of northern flow types over the summer season at this site, suggesting a North Atlantic moisture source with potentially more negative δ 18 O values.Overall, our regional climate model simulations demonstrate significant atmospheric regime shifts over the last glaciation in comparison to the PI period.Additionally, the simulations show less pronounced but still substantial shifts in the position of the eddy-driven jet and associated storm tracks between the stadials and interstadials.We can confirm a more frequent southerly jet position during GS-9/HE4 compared to GI-10.Moreover, the CWT analysis indicate substantial regional circulation changes, with a prevalence of southwesterly and cyclonic CWTs at the locations of the Villars Cave and Schwalbenberg sites during stadial conditions.Further model simulations with other global and regional climate models are needed to verify the validity of the modeled jet characteristics over stadials/interstadials as described in this study.

Data Availability Statement
Proxy reconstructions for temperature and precipitation considered in this study are summarized in the Data Set S1 and are available from the corresponding literature.We used Speleothem data of the Villars Cave in France (Genty et al., 2003(Genty et al., , 2005) ) and of the Sofular cave in Turkey (Fleitmann et al., 2009), data from lake sediments, such as the Grand Pile peat bog insect record (Ponel, 1995), the chironomid record of Unterangerberg (Ilyashuk et al., 2022), and the pollen records of Les Echets and Grand Pile (Guiot et al., 1989), Lac du Bouchet (Thouveny et al., 1994), as well as Lago Grande di Monticchio (Allen et al., 1999).Further, the Schwalbenberg loess sequence in Germany (Prud'homme et al., 2022) is taken into account.
Regional climate model simulations were conducted with the WRF model.The model source code is freely available and can be downloaded from its webpage (https://www2.mmm.ucar.edu/wrf/users/download/get_source.htmlSkamarock et al., 2021).The simulations will be archived at DKRZ (German Climate Computing Centre) and are available from the corresponding author upon request.The ERA5 reanalysis data set (Hersbach et al., 2020) was used for validation of the PI simulation.Graduate School for Climate and Environment (GRACE) at the Karlsruhe Institute of Technology (KIT) for funding a research stay in Budapest and Sopron, Hungary.GU is funded by the NRDIO (OTKA ANN-142613 project) and PL is supported by the Helmholtz "Changing Earth" program.KHS and PL thank the German Climate Computing Centre (DKRZ, Hamburg) for providing computing resources within the DKRZ project 965 "Regional Palaeoclimate Modelling and Palaeoenvironmental reconstructions."The constructive comments of two anonymous reviewers are gratefully acknowledged.Open Access funding enabled and organized by Projekt DEAL.

Figure 1 .
Figure 1.Proxy records of δ 18 O ice , Ca 2+ and temperatures of the NGRIP ice core in central Greenland over the period of 43-38 ka, showing Dansgaard-Oeschger climate variability.Oxygen isotope and Ca 2+ data, the latter representing atmospheric dust variations, are from Seierstad et al. (2014).The paleotemperature reconstruction is from Kindler et al. (2014) and based on δ 15 N data combined with a firn densification and heat diffusion model.Greenland Stadial (GS) and Greenland Interstadial (GI) phases are indicated followingRasmussen et al. (2014).HE4 stands for Heinrich Event 4 (onset and termination afterGuillevic et al. (2014)).All proxies are shown on the GICC05 timescale.

Figure 2 .
Figure 2. Simulated domain and orography for the PI (a), GI (b), and GS (c) simulations.Zoom into Europe with PI orography, grid points show resolution of the simulations (d).Black crosses mark proxy sites used for a quantitative comparison, gray triangles with black edges mark speleothem proxy sites used for a qualitative comparison, and white circled crosses denote center points of the circulation weather type calculations.

Figure 3 .
Figure 3.Comparison of modeled versus proxy-based temperature data.Proxy-based temperature values are plotted against model temperature values of the associated grid point, both including the respective error bars: that denote one standard deviation derived from the original publication for proxy data or from model data, respectively.The shape of the data points denote the proxy type: circles represent data from lake sediments, squares represent data from peat bogs, and rotated squares represent data from loess.Colors denote temperature types: blue stands for January temperatures, orange for mean annual temperatures, different shades of red denote summer temperatures, July temperatures, and summer soil temperatures.Symbols with black outlines refer to the GI-10 period, symbols without outline refer to the GS-9 period.The black diagonal line denotes perfect model-proxy agreement.

Figure 4 .
Figure 4. Distribution of 2 m air temperature in annual mean (left column), winter mean (middle column), and summer mean (right column) as simulated under PI conditions (a-c), GI conditions (d-f), GS conditions (g-i), and 2m air temperature differences between the GS and GI simulation (j-l).Thick black lines denote the respective coastlines, pink lines mark the ice sheet extent (GS in the difference fields).

Figure 5 .
Figure5.Distribution of total precipitation in annual mean (left column), winter mean (middle column), and summer mean (right column) as simulated under PI conditions (a-c), GI conditions (d-f), GS conditions (g-i), and total precipitation differences between the GS and GI simulation (j-l).Thick black lines denote the respective coastlines, pink lines mark the ice sheet extent (GS in the difference fields).

Figure 7 .
Figure 7. Latitudinal distribution of zonally averaged zonal wind over height in annual mean (left column), winter mean (middle column), and summer mean (right column) as simulated under PI conditions (a-c), GI conditions (d-f), GS conditions (g-i).Wind differences between the GS and GI simulation (j-l; colored), contour lines show the absolute wind field under GS conditions.

Figure 10 .
Figure 10.Circulation Weather Types for Schwalbenberg (SN) and Villars Cave (VC) in annual mean (a), in winter (b) and summer (c), and the respective differences between CWTs during the stadial and during the interstadial phase (d-f).