# Millennial-timescale quantitative estimates of climate dynamics in central Europe from earthworm calcite granules in loess deposits

Charlotte Prud'homme, Peter Fischer, Olaf Jöris, Sergey Gromov, Mathias Vinnepand, Christine Hatté, Hubert Vonhof, Olivier Moine, Andreas Vött, Kathryn E Fitzsimmons

<figure><img src="image_1.png" /></figure>

## To cite this version:

Charlotte Prud'homme, Peter Fischer, Olaf Jöris, Sergey Gromov, Mathias Vinnepand, et al. Millennial-timescale quantitative estimates of climate dynamics in central Europe from earthworm calcite granules in loess deposits. Communications Earth & Environment, 2022, 3 (1), pp.267. 10.1038/s43247-022-00595-3. hal-03871063

HAL Id: hal-03871063

https://hal.science/hal-03871063

Submitted on 25 Nov 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Millennial-timescale quantitative estimates of climate dynamics in central Europe from earthworm calcite granules in loess deposits

Charlotte Prud'homme ID <fnref n="1" />,¹⁰✉, Peter Fischer ID ², Olaf Jöris ID ³, Sergey Gromov ID ⁴, Mathias Vinnepand², Christine Hatté ID ⁵,⁶, Hubert Vonhof ID ⁷, Olivier Moine⁸, Andreas Vött² & Kathryn E. Fitzsimmons<fnref n="1" />,⁹

Ice core and marine archives provide detailed quantitative records of last glacial climate changes, whereas comparable terrestrial records from the mid-latitudes remain scarce. Here we quantify warm season land-surface temperatures and precipitation over millennial timescales for central Europe for the period spanning 45,000-22,000 years before present that derive from two temporally overlapping loess-palaeosol-sequences, dated at high resolution by radiocarbon on earthworm calcite granules. Interstadial temperatures were 1-4 °C warmer than stadial climate, a temperature difference which is strongly attenuated compared to Greenland records. We show that climate in the Rhine Valley was significantly cooler during the warm season and overall drier with annual precipitation values reduced by up to 70% compared to the present day. We combine quantitative estimates with mesoscale wind and moisture transport modelling demonstrating that this region was dominated by westerlies and thereby inextricably linked to North Atlantic climate forcing, although ameliorated.

<footnote n="1">¹Research Group for Terrestrial Palaeoclimates, Max Planck Institute for Chemistry, Mainz, Germany. ²Institute of Geography, Johannes-Gutenberg-Universität Mainz, Mainz, Germany. ³MONREPOS Archaeological Research Centre and Museum for Human Behavioural Evolution, RGZM, Schloss Monrepos, Neuwied, Germany. ⁴Department of Atmospheric Chemistry, Max Planck Institute for Chemistry, Mainz, Germany. ⁵Laboratoire des Sciences du Climat et de l'Environnement, UMR 8212 CNRS CEA UVSQ, Université Paris-Saclay, Gif-sur-Yvette, France. ⁶Institute of Physics, Silesian University of Technology, Gliwice, Poland. ⁷Department of Climate Geochemistry, Max Planck Institute for Chemistry, Mainz, Germany. ⁸Laboratoire de Géographie Physique, Environnements Quaternaires et actuels, UMR 8591 CNRS Université Paris 1, UPEC Thiais, France. ⁹Department of Geosciences, University of Tübingen, Tübingen, Germany. ¹⁰Present address: Institute of Earth Surface Dynamics, Faculty of Geosciences and the Environment, Université de Lausanne, Lausanne, Switzerland. ✉email: charlotte.prudhomme@unil.ch</footnote>
Over the last glacial period, the climate of the Northern Hemisphere experienced numerous abrupt variations on millennial to centennial timescales. Our understanding of these changes is underpinned by the Greenland ice core¹ and North Atlantic marine² records; terrestrial responses and feedbacks are less well documented. Climatic variations include the Dansgaard-Oeschger cycles which comprise mild Greenland Interstadial events (GIs), defined by sudden decadal warming followed by gradual centennial cooling, bundled with Greenland stadial phases (GSs), characterized by abrupt cooling and more stable cold conditions¹,³,⁴, as well as short-lived, particularly cold Heinrich Events which are evidenced by coarse-grained layers of ice-rafted debris preserved within North Atlantic marine sediments⁵. Heinrich events are associated with a significant drop in sea-surface temperatures (SST)⁶, are coeval with the most prominent GSs, and typically occurred towards the end of sequences of GIs of decreasing amplitude (Bond-cycles)⁷. These short-term climate variations were likely triggered by a sudden slowdown of the Atlantic Meridional Overturning Circulation (AMOC) and are associated with increased North Atlantic sea ice cover⁸.

Within the terrestrial realm, these changes manifested in the landscape and ecology, yet quantitative reconstructions of these responses elude us still. On the European continent, evidence of millennial-timescale climate variations similar to the North Atlantic and Greenland records has been observed within loess-palaeosol sequences⁹⁻¹⁴, lake sediments¹⁵⁻¹⁸, continental margins¹⁹,²⁰, speleothems²¹⁻²³ and tufa²⁴, yet so far quantitative data on climate parameters such as temperature and rainfall are seldom available. These limitations constrain our understanding of the different influences of Atlantic air masses or the extent of ice sheets on continental Europe in response to climatic changes²⁵.

Here we provide quantitative estimates of temperature and precipitation for central Europe from terrestrial loess deposits at millennial-timescales comparable with the Greenland ice core and North Atlantic marine records. Compared with other terrestrial archives, loess deposits are widespread within the ca. 45–55° N latitudinal belt²⁶⁻²⁸ and characterized by high sedimentation rates within which climatic information is preserved. Loess accumulation rates of up to 1 mm yr⁻¹ ensure high resolution coverage of (sub-)millennial climatic cycles, expressed in the stratigraphy by alternating stadial primary loess and interstadial palaeosol units⁹,¹¹,¹²,²⁹,³⁰. The physical properties of loess sediments depend on the intensity and duration of interstadial and stadial conditions and sedimentation rates¹²,³¹; past climatic reconstruction is additionally facilitated by analysis of fossil mollusc assemblages³²⁻³⁸, land snail shell carbonate clumped isotope thermometry¹⁰, the isotopic signal of organic matter³⁹,⁴⁰, and lipid and plant wax biomarkers⁴¹⁻⁴³.

We make use of a singular material for quantitative, high resolution palaeoclimate reconstruction from loess sediments, which can be analyzed not only as palaeothermometers⁴⁴ and palaeo-rainfall gauges⁴⁵, but can also be directly dated using radiocarbon¹¹,¹²; fossil earthworm calcite granules (ECGs). ECGs are secreted daily at the soil surface⁴⁶ by various earthworm species (*Lumbricus terrestris* being the most productive⁴⁷), comprise rhombohedral low-Mg calcite crystals which are resistant to diagenetic alteration and experience limited vertical mixing within the loess sediment column⁴⁸,⁴⁹. The oxygen and carbon isotopic compositions of ECGs provide direct proxies for warm-season land-surface temperature (LSTws)⁴⁴ and precipitation⁴⁵.

Our quantitative dataset based on coupled stable isotope composition and radiocarbon dating of ECGs provides a record of high temporal resolution, spanning 45,000–22,000 years cal BP (late Marine Isotope Stage 3 and MIS 2), from two temporally

overlapping loess-palaeosol sequences. The Schwalbenberg¹²,⁵⁰ and Nussloch¹¹,²⁹ sequences are located c. 180 km from one another along a northwest-southeast transect in the non-glaciated zone in the Rhine River valley of western Germany and therefore reflect terrestrial responses to oscillating climates of the last glacial period for central Europe (Fig. 1). We embed our climate reconstructions within Bayesian age models for both sequences, facilitating precise correlation between both one another and with other archives from the North Atlantic region. We combine our results with continental-scale simulations of wind and precipitation trajectories for the last glacial maximum in order to more meaningfully understand potential teleconnections between continental Europe, the Fennoscandian ice sheet and the North Atlantic Ocean.

## Results

**Stratigraphy.** The Schwalbenberg profile is located on the left bank of the Rhine in its Middle Valley, at the confluence of the Rhine and Ahr rivers near the city of Remagen. The locality forms a gently south-eastward sloping plateau¹²,⁵⁰ (Fig. 1); sediment thickness reaches 30 m and spans the last Interglacial-Glacial Cycle¹². The aeolian sediments at Schwalbenberg are assumed to derive from adjacent river systems and the Rhenish Massif uplands. The 6 m-thick Schwalbenberg RP1 profile is exposed on the southeastern flank of the plateau and spans late MIS 3 and MIS 2, with a distinct unconformity separating the lower from the upper part of the section. Late MIS 3 is resolved in high resolution by multiple pedocomplexes (Calcaric Cambisols) and a tundra gley, overlain by a cyclic succession of MIS 2 loess and tundra gley horizons (Fig. 2)¹²,⁵⁰.

Nussloch is situated on the eastern margin of the Upper Rhine Valley (Fig. 1). The aeolian sediments of this locality are up to 17 m thick and are assumed to derive from local deflation of the Rhine alluvial plain to the west, accumulating at the junction between the Kraichgau and Odenwald uplands as NNW-SSE elongated ridges named “gredas”²⁹. A combination of topography and strong northwesterly wind regimes²⁹ facilitated high dust accumulation rates at Nussloch during MIS 2. The lower 5 m dates to late MIS 3 and comprises arctic and boreal brown soils alternating with primary loess and tundra gley horizons, overlain by 12 m of MIS 2 age alternating primary loess and tundra gley layers (Fig. 3)²⁹.

Both Schwalbenberg and Nussloch preserve evidence of periglacial and permafrost activity to varying degrees. They are situated between the Last Glacial Maximum (LGM) extents of the Alpine and Fennoscandian ice sheets. Tundra gley horizons are present at both sites and reflect short-term phases of formation and subsequent degradation of active permafrost layers²⁹. During GI 6, a Calcaric Cambisol developed at Schwalbenberg, coeval with tundra gleys at Nussloch, suggesting that permafrost features were significantly weaker at Schwalbenberg compared to Nussloch, possibly due to the leeward position and southerly aspect of the former.

**Age-depth models and temporal resolution.** Our chronologies for both Schwalbenberg (RP1, Fig. 2) and Nussloch (NUSS P8, Fig. 3) profiles were established using radiocarbon dating of ECGs (Supplementary Data 1). We undertook Bayesian age modelling using Bacon software⁵¹ (see “Methods”) (Fig. 4 and Supplementary Data 1).

Our chronology for Schwalbenberg RP1 is based on 22 dates and spans 39,963–38,162 to 22,150–21,426 cal BP (Fig. 2). We constrain the timespan of the unconformity between units 12 and 13 as a ~5 ka hiatus spanning 30,970–30,199 to 24,619–23,863 cal BP. Sedimentation rates range from ~0.5 to 0.6 mm yr⁻¹ for MIS
<figure><img src="image_4.png" /><figcaption>**Fig. 1** DEM-based map showing the Palaeogeography of western and central Europe during the Last Glacial Maximum (LGM ~21 ka). The distance between Schwalbenberg (50.56°N; 7.24°E) and Nussloch (49.18°N; 8.70°E) loess sites is ~180 km. Yellow patches represented areas of loess and loess derivates<sup>26,28</sup>. Beige-coloured areas represent the dry continental shelves<sup>28</sup>. The DEM is based on combined DEMs<sup>26</sup> of the GLOBE Project (1 km resolution DEM).</figcaption></figure>

3 and post 22 cal kBP, reaching a peak of 2.0 mm yr⁻¹ between 23 and 22 cal kBP.

We undertake Bayesian age-depth modelling and calculation of sedimentation rates based on published ECG data for the Nussloch P8 profile¹¹. Age-depth modelled calibrated dates range from 47,862–45,354 to 22,752–21,577 cal BP, showing that the Nussloch P8 profile predates Schwalbenberg RP1 by c. 7 ka (Fig. 3). Between 43 and 32 cal kBP, Nussloch sedimentation rates are comparable to those at Schwalbenberg RP1 (~0.5 to 0.6 mm yr⁻¹). Short-lived peaks (up to 1 mm yr⁻¹) in dust accumulation during MIS 2 were constrained to c. 30, 25.5, and 24 cal kBP.

**Stable isotope composition.** ECGs were collected from 5 cm thick vertical slices according to published methods, 23 samples from Schwalbenberg and 18 from Nussloch. From each sample we extracted 30 granules, measuring the carbon and oxygen isotope values on 690 individual ECGs from Schwalbenberg (Fig. 4, Supplementary Data 2, and Supplementary Fig. 1). The Nussloch record incorporates an already published stable isotope dataset (12 samples representing 360 analyses)⁴⁴,⁴⁵, to which an additional 180 analyses from 6 samples were undertaken (Fig. 4, Supplementary Data 2, and Supplementary Fig. 2).

At Schwalbenberg, ECGs extracted from the interstadial phases (9 samples, representing 270 analyses) yielded mean values of -12.5 ± 1.0 ‰ for δ¹³C and -4.5 ± 1.8‰ for δ¹⁸O. Results for samples from stadial phases (14 samples, representing 420 analyses) yielded mean δ¹³C and δ¹⁸O values of -12.2 ± 1.0‰

and -4.8 ± 1.8‰ respectively. ECG stable isotope results are statistically different for interstadial and stadial phases; T-test results yield *p*-values of 0.04 and 8.9 × 10⁻⁶ for δ¹³C and δ¹⁸O respectively.

The mean carbon isotopic composition values for ECGs from Nussloch are -12.4 ± 0.9 ‰ for interstadials (14 samples representing 420 analyses) and -11.7 ± 0.9 ‰ for stadial phases (4 samples representing 120 analyses). Mean oxygen values of ECGs during interstadials and stadials are -5.2 ± 1.9‰ and -6.3 ± 1.7‰ respectively. T-test values of 1.6 × 10⁻¹⁵ and 2.3 × 10⁻⁸ for our mean carbon and oxygen isotope ratio values also indicate significant statistical differences between interstadial and stadial phases.

**Quantitative estimates of terrestrial climate parameters.** Transfer functions were used to convert the ECG stable carbon and oxygen isotope compositions into estimates of annual precipitation⁴⁵ and warm-season land-surface temperature⁴⁴.

Since the main source of carbon in ECGs derives from a mixture of litter and soil organic matter⁵²,⁵³, the δ¹³C values of plants can be estimated by accounting for the fractionation factor (-11.7 ± 1.5‰) that results from the difference between the carbon ingested by the earthworms and that which is integrated into the ECGs⁵⁴. This fractionation factor includes the carbon isotope fractionation between amorphous calcium carbonate (ACC) and calcite produced by the earthworm (ε calcite-ACC = -1.2 ± 0.5‰)⁵⁵. The mean δ¹³C values of plants
<figure><img src="image_5.png" /><figcaption>Fig. 2 Age-depth model of the Schwalbenberg loess profile RP1. Age-depth model is based on Bayesian modelling using the software Bacon⁵¹ with 28 radiocarbon ages calibrated with IntCal20. The chronology of the profile is correlated with the δ¹⁸O values of NGRIP³. SU - sediment units. Tundra gley horizons are represented in grey and palaeosols (Cambisols) in brown.</figcaption></figure>

calculated from the carbon isotope composition of the ECG at Schwalbenberg (−24.0 ± 1.0‰) and at Nussloch (−24.0 ± 0.9‰) for the period from 45 to 20 cal kBP indicate that the main source of organic matter derived from C₃ plants⁴⁰. Following calculation of the δ¹³C values of plants, we used a published empirical equation assuming a linear relationship between C₃ plant leaf discrimination (Δ¹³Cplants) and mean annual precipitation (MAP) computed by ref. ⁵⁶ compiled from 480 modern sites⁵⁷ (Eq. 1).

$$ \text{Log(MAP} + 300(\text{mm.yr}^{-1})) = 0.092(\pm 0.004)\Delta^{13}\text{C}_{\text{plants}} + 1.148(\pm 0.074); $$

$$ (R^2 = 0.56; \sigma = 0.14) \text{ with } \Delta^{13}\text{C}_{\text{plants}} = \frac{\delta^{13}\text{C}_{\text{atm}} - \delta^{13}\text{C}_{\text{plants}}}{1 + \delta^{13}\text{C}_{\text{plants}}} \quad (1) $$

The Δ¹³Cplants is dependent on the δ¹³C of plants (‰ VPDB) and calculated from the δ¹³C of ECG and from δ¹³C values of atmospheric CO₂ (δ¹³Catm) which were estimated from air bubbles trapped within the Vostok ice core between 50 and 20 ka⁵⁸. We propagated the root mean squared error (σ = 0.14) through our calculations of MAP (Supplementary Data 2).

ECG oxygen isotope composition is related both to crystallization temperature and the isotopic composition of soil water⁵⁹. Based on an empirical equation using observations of modern earthworm response to temperature and taking into account the vital effect (Eqs. 2 and 3)⁵⁹, it then remains to identify the δ¹⁸O value of soil water in order to quantify the temperature.

$$ 1000\ln(\alpha_{\text{calcite-water}}) = [20.21 \pm 0.92](10^3 \cdot T^{-1}) - [38.58 \pm 3.18] \quad (2) $$

$$ \text{with } \alpha_{\text{calcite-water}} = \frac{[1 + \delta^{18}\text{O}_{\text{ECG}}]}{[1 + \delta^{18}\text{O}_{\text{mw}}]}; \quad (3) $$
T(K) and δ¹⁸O (‰V - SMOW)

Based on clumped isotope data, the δ¹⁸O values of soil water are not clearly biased toward seasonal δ¹⁸O values of meteoric water (δ¹⁸Omw) and are primarily controlled by the annual δ¹⁸Omw⁶⁰,⁶¹. In periglacial environment, δ¹⁸O values of soil water are likely dependent on freeze/thawing cycles. Therefore, as a first approximation, we assumed that the δ¹⁸O values of the granules
<figure><img src="image_6.png" /><figcaption>Fig. 3 Age-depth model of the Nussloch loess profile P8. Age-depth model is based on Bayesian modelling using the software Bacon⁵¹ with 29 radiocarbon ages¹¹ calibrated with IntCal20. The chronology of the profile is correlated with the δ¹⁸O values of NGRIP³. Tundra gley horizons are represented in grey and arctic/boreal brown soils in brown.</figcaption></figure>

are dependant of the annual δ¹⁸Omw values, which are themselves linearly dependent on mean annual air temperature at mid- to high latitudes⁶²,⁶³.

We established a linear equation (Eq. 4) relating mean annual weighted δ¹⁸Omw to mean annual air temperature (Tₐ) by compiling modern rainfall data from 80 mid-high latitude European stations from the Global Network of Isotopes in Precipitation (GNIP-IAEA) (Supplementary Data 3 and Supplementary Fig. 3):

$$ \delta^{18}O_{mw} = [0.68 \pm 0.08]Ta - [14.7 \pm 0.6]; R^2 = 0.80 \quad (4) $$

To better constrain our paleotemperature estimates, we used a linear equation relating daily ground temperature (Tₐ) to daily air temperature (Tₐ) under snow-free conditions at West Dock and Franklin Bluffs in Alaska, for the time period spanning from 1987 to 1992 as follows⁶⁴,⁶⁵:

$$ T_{soil} = 0.89 * T_{a} + 3.40 \quad (5) $$

This equation is in agreement with the assumption that air temperatures are ~3 °C colder than soil temperatures (Supplementary Data 2) as is typical for present-day periglacial environments⁶⁴,⁶⁶.

The time interval covered in this study (late MIS 3-2) oversaw substantial changes in hemispheric ice volume and, consequently in the oxygen isotope ratios of both sea and meteoric water. These needed to be accounted for in our calculations. During MIS 3 and MIS 2, a drop in sea level between ~60 and 120 m below the modern datum⁶⁷ contributed to a major reshaping of northwest European palaeogeography (Fig. 1)⁶⁸. The δ¹⁸O of seawater at the

MIS 2 lowest stand was 0.8 ± 0.1‰⁶⁹; we took MIS 3 δ¹⁸O values to be intermediate between MIS 2 and the Holocene (0 ‰; SMOW reference)⁷⁰,⁷¹ and assume a value of 0.5 ‰.

We calculated land-surface temperatures for the season of greatest ECG production by estimating the δ¹⁸O values of meteoric water (Eq. 4), taking into account the ice volume effect and by combining Eqs. (2), (3), (4) and (5). Since earthworms produced ECGs during the active layer thawing season⁴⁹ (May to September, the five warmest months of the year; data from Global Terrestrial Network for Permafrost http://gtnpdatabase.org), we propose that the oxygen isotope composition of ECGs records LSTws rather than mean annual land-surface temperature (LST). We reported the distribution of the 30 values for each sample of 5 cm thick at 1σ (Fig. 5).

At Schwalbenberg, LSTws range from 8.8 ± 4.7 to 12.9 ± 3.3 °C, and MAP values from 215 ± 172 to 518 ± 270 mm yr⁻¹ (Supplementary Data 2). Nussloch temperatures were cooler (5.0 ± 6.7 to 12.0 ± 5.5 °C) but recorded comparable annual precipitation (224 ± 175 to 490 ± 265 mm yr⁻¹) (Supplementary Data 2). The Rhine Valley is particularly cold with a mean value of 7.3 ± 4.7 °C during the GS 3 (26–25 ka), a stadial coeval with Heinrich Stadial 2. At Nussloch, the GS 5.1 is slightly cooler with a mean LSTws of 5.0 ± 6.7 °C. The timing of the warmest summer phases varies between the two sites; the warmest temperatures occurred at Schwalbenberg during GI 7 (~35 ka) and during GI 12 (~45.5 ka) at Nussloch, although it should be noted that GI 12 was not exposed at Schwalbenberg RP1. The timing of the driest conditions also varies between the sites: at Schwalbenberg these prevailed ~22 ka (during the LGM), and at Nussloch were coeval
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Stable isotope composition of ECGs from Schwalbenberg and Nussloch. **a** Histogram of oxygen and carbon compositions for both sequences; **b** Schwalbenberg (dark-coloured plain lines) and Nussloch (light-coloured dotted lines) δ¹³C (‰ VPDB) and δ¹⁸O (‰ VPDB) values plotted against the Bayesian age-depth model of the radiocarbon dated sequence using IntCal20 (black horizontal error bars represent 95% confidence levels for modelled ages). Error bars on stable isotope values represent the standard deviation (1σ) of 30 individual measurements per sample of 5 cm thickness.</figcaption></figure>

with the coldest phase (Henrich stadial 2). The wettest conditions occurred during the transition from GI 6 to GS 6 (~33 ka), and during GI 5.2 (32 ka), at Schwalbenberg and Nussloch, respectively. We note that apparent discrepancies in the timing of extreme values between sites reflects the lack of temporally overlapping data for the respective phases.

**Climate simulations: modelling of wind and precipitation regimes.** Our simulations of wind regimes and precipitation arriving at Nussloch and Schwalbenberg were undertaken for two time slices: the present day (1989–2018) and LGM (21 ka) (Fig. 6 and Supplementary Figs. 4, 5). Our results indicate that for both time periods, the two sites experience spatial differences in general circulation and the trajectories of air masses responsible for dust transport. At Nussloch in the present day, south-westerlies prevail throughout the year, with reduced wind intensities in the summer months and steady winds in the winters; during the

LGM, north-westerly and southerly winds dominated which were steady in summer. By contrast, Schwalbenberg presently experiences more unstable, stronger, predominantly south-westerly winds, with higher speeds in autumn and winter; by contrast, during the LGM only summer west/north-westerly winds were comparable with the present day in intensity and steadiness, with easterlies more pronounced in other seasons (Fig. 6 and Supplementary Figs. 4, 5). LGM precipitation derived from south-westerly air masses at Schwalbenberg, whereas at Nussloch, precipitation derived from westerlies.

## Discussion

Previous studies have demonstrated that the Schwalbenberg¹²,⁷²,⁷³ and Nussloch¹¹,¹⁴,²⁹,³⁴,⁴⁹,⁷⁴ sequences preserve stadial-interstadial climate variability at high temporal resolution. Our new age modelling and geochemical data, combined with wind regime modelling for the present day and LGM time slices, provide the first
<figure><img src="image_8.png" /></figure>

opportunity to meaningfully investigate terrestrial palaeoclimate
change in central Europe over millennial timescales. Our quanti-
tative climatic estimates are derived from in situ ECG response to
contemporary conditions at the land surface and are therefore
independent of sedimentation rates or direct wind transport vec-
tors. As such, they provide insights into terrestrial climate para-
meters at sufficiently high resolution for comparison with other
archives from Europe and the North Atlantic region (Fig. 5).

Our estimations for the period 45.5-22.5 cal kBP indicate
temperatures 3.3-11.8 °C cooler and significantly drier condi-
tions than experienced in the present day (Table 1). Modern
(1961-1990) warm season mean temperatures at Schwalben-
berg and Nussloch are 16.2 and 16.8 °C and modern rainfall is
654 and 784 mm yr⁻¹ respectively⁷⁵. During the last glacial
period, the northern site of Schwalbenberg yields higher LSTws
than Nussloch throughout, which may be attributed to local
**Fig. 5** Climate quantitative estimates and dust variation in terrestrial and glacial environments across the North Atlantic Region. **a** June insolation at 60°N<sup>108</sup>; **b** δ<sup>18</sup>O ratio and **c** Ca<sup>+</sup> concentration of NGRIP Greenland ice core<sup>3</sup>. Originally given in b2K ages, this latter chronological scale has been shifted by 50 years to allow for direct comparisons with the radiocarbon dating given in years BP (before 1950). **d** Iberian margin sea surface temperature (SST, MD01 2444)<sup>6</sup>, radiocarbon dating calibrated with IntCal13; **e** sedimentation rates from Nussloch (yellow) and Schwalbenberg (red) profiles, based on Bayesian age-depth modelling of the radiocarbon dated sequence using IntCal20; **f** mean annual precipitation estimates based on the δ<sup>13</sup>C values of ECGs for both loess profiles (Nussloch: light-coloured dotted lines and Schwalbenberg: dark-coloured plain lines; thick colour error bars represent SE), plotted against the Bayesian age-depth model of the radiocarbon dated sequence using IntCal20 (black horizontal error bars represent 95% confidence levels for modelled ages); **g** land-surface temperature of the warmest month estimated from the δ<sup>18</sup>O values of ECGs for both loess profiles (Nussloch: light-coloured dotted lines and Schwalbenberg: dark-coloured plain lines; thick colour error bars represent 1σ confidence levels), plotted against the Bayesian age-depth model of the radiocarbon dated sequence using IntCal20 (black horizontal error bars represent 95% confidence levels for modelled ages); **h** clumped temperature from land snail shells of *T. hispidus* (circle) and *S. oblonga* (square), error bars represent 68% confidence levels in the Dunaszekcső record (Hungary) referred to age-depth model of the radiocarbon dated sequence using IntCal13<sup>9,10</sup>; **i** bulk mass accumulation rates (BMARs) at the Dunaszekcső site<sup>9</sup>; **j** δ<sup>18</sup>O of meteoric water from the Alps based on the δ<sup>18</sup>O of the Sieben Hengste Cave (7H) speleothem<sup>21</sup>, U-Th dating in ka; **k** European ice sheet volume reconstruction (SIS)<sup>109</sup>. H2-H5 – Heinrich stadials (in blue). Palaeosols formed during long Greenland Interstadials (>1ka) are highlighted in reddish (GI 12, 11, 8) and tundra gleys developed during short Greenland Interstadials (<1ka) are in yellow (GI 10, 9, 7, 6, 5, 2, 1, 4, 3, 2, 2.1)<sup>31</sup>.

influences such as the south-facing aspect (i.e., insolation was more intensive at this locality). We observe comparable mean annual precipitation values and variability throughout the investigated time period at both Schwalbenberg and Nussloch, with Nussloch being even drier over the entire period covered (Fig. 5).

Overall, MIS 2 was c. 2.5 °C cooler than late MIS 3 at both sites and is interpreted to relate to the global cooling towards the LGM (Fig. 5). We are aware that uncertainties associated with our quantitative estimates of climate parameters are sizeable (Supplementary Data 2, “Methods”). Nevertheless, statistical tests (T-test) are showing that the mean values of LSTws between MIS 3 and MIS 2 are statistically different at both site (Supplementary Table 2). Our estimates are in agreement with other (lower resolution) quantitative estimates for past climate parameters in terrestrial central Europe based on palaeoecological transfer functions (malacofauna<sup>35,37,38</sup>, chironomids<sup>76</sup> and diatoms<sup>77</sup>, biomarker-based studies<sup>78,79</sup>), stable isotope composition of loess carbonates<sup>44,80</sup> and of large mammal remains<sup>81</sup> which provide mean temperatures for the warmest season ranging from 14 to 8 °C during MIS 3 and MIS 2. We suggest that the milder climatic fluctuations observed in the terrestrial mid-latitudes were driven at multiple scales by insolation and local topography. At local scales, direct insolation may have led to comparably high local summer temperatures, even during colder periods (MIS 3 stadials). Our mean annual precipitation estimates for the Rhine valley are consistent with the precipitation values modelled by the BIOME 4 inverse model simulation using δ<sup>13</sup>C of soil organic matter<sup>39,74</sup>, revealing driest (and coldest) values for the LGM. In comparison, the late MIS 3 was characterized by wetter conditions. Observations of palaeosol development indicate increased vegetation density at both sites at this time and agree with vegetation changes indicated by the Schwalbenberg and Nussloch mollusc fauna<sup>73,82</sup>, Nussloch biomarker-based studies<sup>43</sup> as well as palaeoecological reconstructions from the nearby Eifel maar lakes<sup>15,17</sup> and Bergsee<sup>18</sup>.

At millennial-timescale, our data show that the stadial phases in the Rhine Valley were on average 1–4 °C cooler during summer and significantly drier with less than half of the interstadial precipitation (Fig. 5, Supplementary Data 2). The lowest MAP estimates occurred during GS2.1 at Schwalbenberg and during GS3 at Nussloch with precipitation reduced by 67 and 71% respectively from present-day averages. The amplitude of ~10 °C difference between temperature from the coldest stadials (GS 3) and today in the Rhine Valley is in agreement with groundwater noble gas records which indicate a 7–10 °C terrestrial cooling at 40°N during the LGM<sup>83</sup>.

The interstadial between the Heinrich Stadials 2 and 3 (GI 3) recorded warmer temperatures (11.3 ± 5.0 °C) at Nussloch, a

phase which was characterized by abrupt warming in the Greenland region<sup>84</sup> (Fig. 5). We compare our estimates with growing season temperatures based on clumped isotope data from land snails in the Dunaszekcső loess profile of southern Hungary. The interstadials GI 5.1 and GI 3 were 7 °C warmer than at Nussloch<sup>10</sup>. A cooling between 4 and 7 °C was quantified for the stadial phases during 31–26.5 cal kBP<sup>10</sup>, suggesting a more continental climate in eastern Europe than central Europe during the last glacial period. This might be due to the more eastern and southern air mass influences in eastern Europe compared to central Europe which is mostly influenced by North Atlantic climate and the westerlies<sup>85</sup>. Continental response at the transition from stadial to interstadial condition in Europe was—within dating uncertainties – both rapid and synchronous, albeit not as pronounced as in the Greenland (+10 to 15 °C) and Mediterranean Sea regions (SST: +3 to 5 °C) (Fig. 5). During MIS 2 interstadials, high concentrations of mollusc shells combined with high juvenile mortality rates suggest increasing summer temperatures with abrupt temperature decreases during the transitions into the winter season reflecting contrast in seasonality during interstadials<sup>34,86</sup>. High contrast seasonality might also influence stadial phases in the Rhine Valley, as permafrost evidence in the stratigraphy<sup>12,29</sup> suggest particularly cold winter temperatures.

Increased precipitation in the Rhine Valley is recorded between 33 and 31 cal kBP (GI 6 and GI 5.2) and is coeval with the onset of the build-up of the British Island Ice Sheet (BIIS) and the covering of the Baltic Sea by the Fennoscandian ice sheet (FIS)<sup>87</sup>. Subsequent to the rainfall increase, stepwise drying is then recorded at Nussloch over the millennia from 31 to 23 cal kBP and correlates with both a weakening of the AMOC and synchronous and rapid growth of the FIS<sup>87–89</sup>, enhancing continental effective aridity on a global scale<sup>89</sup>. The isotope composition of organic carbon along with the litho- and pedostratigraphy of the Schwalbenberg REM3 key loess-palaeosol sequence<sup>12,72</sup> documents this period of overall increasing effective aridity in the formation of thick loess deposits reaching up to 8 m (superordinate sediment unit SSU E) within an interval of approximately 5,000 years duration. This unit (SSU E) is, however, strongly eroded and truncated in the Schwalbenberg RP1 section. During this time period, the Laurentide Ice Sheet reached its maximum extent (GS 3, from 27 to 23 ka) leading to significantly reduced precipitation in north-western Europe<sup>90</sup>, which is also reflected by the absolute minimum at Nussloch. This also resulted in a southward shift of the jet stream over southern Europe<sup>85,91,92</sup> immediately preceding the LGM, suggesting a weakening of the North Atlantic Oscillation (NAO)<sup>93</sup> associated with a strengthened Siberian High<sup>9</sup> and an increase of the sea ice cover<sup>94</sup> thereby leading to a reduction of the moisture amount over the North Atlantic with a weakening of the westerlies and a cooling of the air masses. These
<figure><img src="image_9.png" /><figcaption>**Fig. 6** Annual wind regime and precipitation statistics. The models are obtained for present day and LGM conditions at Nussloch (a, b) and Schwalbenberg (c, d). Shaded area presents the average precipitation flux at given wind direction/strength (percent of total). Sector line plots show wind speed (black), wind occurrence (red solid), precipitation rate (blue dashed) and total precipitation (red dashed, normalized) distribution about rhumbs. Thin solid/long-dashed line circles visualize the scalar/absolute winds and their directions, respectively; thin-dashed circle shows wind vector. SD. Note that wind direction is meteorological (from where wind blows). Further statistics tabled: Vs and VR — absolute and scalar wind speed [m/s], Rw and Rp — rhumbs of maximum wind and precipitation occurrence [percent], Q — wind steadiness (Vs / VR ratio), rxy — wind zonal and meridional components correlation, σr — wind vector SD [m/s], λ — wind absolute speed to SD ratio, Σ(P): total annual precipitation [mm].</figcaption></figure>

mechanisms have probably triggered pre-LGM aridification and cooling in the Rhine Valley.

The integration of sedimentation rate calculations (based on age-depth modelling) and wind regime reconstructions suggests that the two Rhenish sites provide contrasting information relating to atmospheric circulation. At Nussloch, dust accumulated almost continuously from 50 to 20 cal kBP and reach peaks during Heinrich stadials 3 and 2 from 30 to 23 cal kBP. The timing of Nussloch peak sedimentation rates was coeval with dust peaks in the Greenland ice cores and suggest periods of stronger winds and drier conditions at Nussloch, likely to indicate phases of significant transitions in Northern Hemisphere wind regimes (Fig. 5). This trend is also observed in eastern Europe from the Dunaszekcső loess-palaeosol profile⁹. During this time period, the sedimentation rates of both Nussloch and Dunaszekcső loess sequences correlated with the oxygen isotopic composition of the Sieben Hengste Cave (7H) used to reconstruct changes in moisture sources in the Alps²¹, implying that the increased dust deposition over western and eastern central Europe was caused by large-scale atmospheric re-organization⁹ (Fig. 5). By contrast, the

Schwalbenberg RP1 section located down slope at the rim of the plateau was prone to erosion at c. 25 cal kBP removing sediments that would cover the period 30–25 cal kBP (Fig. 5). However, it is important to note that this period is characterized by a significant increase of the sedimentation rates in the REM 3 core located in the interfluve position¹²,⁷².

Our combined geochemical and climate modelling results indicate that not only westerly but also southwesterly and northwesterly winds were responsible for LGM moisture transported into the Rhine Valley, and that the steadiness of the dominant winds varied significantly between the seasons (Supplementary Figs. 4, 5). Our wind and precipitation reconstructions for the present day and LGM time slices primarily reflect changes to the general circulation pattern: the hemispheric extent of the extratropical jet, and typical location and intensity of (anti)cyclonic systems. We assume that the general circulation pattern is important for identifying dust and atmospheric moisture transport pathways⁸⁵,⁹⁵ (Fig. 6 and Supplementary Figs. 3, 4). During the LGM, the main source of moisture to the Rhine Valley derived from the westerlies (SW
<table><caption>Table 1 Temperature and precipitation difference between modern and late MIS 3-MIS 2 at Schwalbenberg and Nussloch.</caption><thead><tr><th>Time period</th><th>LSTwS (mean °C)</th><th>Min/ max (°C)</th><th>Δ Mean T to 1961-1990 min- max (°C)</th><th>MAP (mm)</th><th>Min/max (mm)</th><th>Δ Mean MAP to 1961-1990 min- max (-%)</th></tr></thead><tbody><tr><td>Entire period covered</td><td>11.2</td><td>8.8/12.9</td><td>3.3-7.4</td><td>350</td><td>215/518</td><td>21-67</td></tr><tr><td>Overlap 40-31 cal kBP</td><td>11.8</td><td>10.4/12.9</td><td>3.3-5.8</td><td>369</td><td>252/518</td><td>21-61</td></tr><tr><td>Overlap 24-23 cal kBP</td><td>9.5</td><td>8.8/10.2</td><td>6.0-7.4</td><td>314</td><td>293/335</td><td>49-55</td></tr><tr><td>Interstadials</td><td>11.7</td><td>10.2/12.9</td><td>3.3-6.0</td><td>386</td><td>293/462</td><td>29-55</td></tr><tr><td>Stadials</td><td>10.8</td><td>8.8/12.4</td><td>3.8-7.4</td><td>326</td><td>215/518</td><td>21-67</td></tr><tr><td>Nussloch (1961-1990: 16.8 °C/MAP 784 mm)</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Entire period covered</td><td>9.2</td><td>5.0/12.0</td><td>4.8-11.8</td><td>342</td><td>224/490</td><td>38-71</td></tr><tr><td>Overlap 40-31 cal kBP</td><td>9.9</td><td>9.0/10.8</td><td>6.0-7.8</td><td>389</td><td>338/490</td><td>38-57</td></tr><tr><td>Overlap 24-23 cal kBP</td><td>9.3</td><td>7.3/11.3</td><td>5.5-9.5</td><td>260</td><td>224/297</td><td>62-71</td></tr><tr><td>Interstadials</td><td>10.3</td><td>7.7/12.0</td><td>4.8-9.1</td><td>360</td><td>297/490</td><td>38-62</td></tr><tr><td>Stadials</td><td>6.7</td><td>6.2/7.7</td><td>9.1-10.6</td><td>270</td><td>224/314</td><td>60-71</td></tr></tbody></table>

Modern climate parameters (land surface temperature from May to September, LSTwS and mean annual precipitation, MAP) are calculated from ref. 75 for the period 1961-1990 and past climate parameters are from Supplementary Data 2. The time-periods 40–31 and 24–23 cal kBP represent the intervals where the two sequences overlapped. Data between stadials and interstadials are reported in Supplementary Data 2.

and NW), which may have been weakened during phases of ice sheet volume expansion due to a reduction of evaporation over the North Atlantic Ocean and associated westward expansion of the high-pressure centres over the BIIS and FIS. Cooler and drier conditions during such episodes (GS) are also likely to have affected seasonality by shortening warm seasons (Figs. 5, 6 and Supplementary Figs. 4, 5).

To conclude, this study generates new millennial-timescale quantitative terrestrial climate estimates for the Rhine Valley, in central Europe. We see our approach as a first step towards setting up spatially widely distributed data on land-surface temperatures based on ECGs derived from loess-palaeosol sequences, which can be compared with other methods such as land snail shell carbonate clumped isotope thermometry. Our air circulation models for the present-day and LGM time slices help us contextualize the impact of extreme climate phases on terrestrial environments. Our dataset hints at significant correlation between changes in sea ice cover and extent of periglacial aeolian environments (such as the Rhine Valley deposits) over continental Europe during the second half of the last glacial period. Although the correlation between temperature and precipitation recorded within the Rhine Valley with higher latitude ice sheet dynamics suggests a passive, and more moderate response to hemispheric climate drivers, such environments provide an appropriate and quantitative illustration of changes occurring during abrupt climatic transitions. Overall, we observe comparable trends between loess records from western and eastern central Europe and isotopic proxies from the Greenland ice cores confirming that the Schwalbenberg and Nussloch and archives from central Europe were directly influenced by changes in the North Atlantic climate during the last glacial period.

## Methods

**Sampling strategy.** Continuous columns of 10 litres samples were taken throughout the 5.5 m-thick Schwalbenberg profile RP1 at Remagen and throughout the 17 m-thick profile P8 at Nussloch. Schwalbenberg RP1 is located immediately at the northern edge of a former archaeological excavation site¹²,⁹⁶. Nussloch P8 is situated about 50 m southeast of the reference profile P4²⁹. Before sampling, we refreshed and cleaned both sections by removing vertically about 50 cm of sediments. ECG were sampled at 5 cm resolution by collecting 10 l of sediments. After wet sieving at 0.5 mm, ECG with a diameter larger than 0.8 mm were extracted under a binocular microscope. For stable isotope analyses, we selected 23 samples throughout the profile RP1 at Schwalbenberg (Supplementary Fig. 1). For Nussloch profile P8, we analyzed 7 new samples (each 5 cm thick) that we added to the 11 previously published⁴⁴,⁴⁵ for this profile (Supplementary Fig. 2). ECGs were selected from tundra gley units and brown soil horizons and loess units.

**Radiocarbon dating.** Sample preparation and radiocarbon analysis protocol followed the procedure as described by ref. ¹¹. After cleaning the ECGs by ultrasonic bath, 100 µg (gas source) or 1000 µg (solid source) per dating sample were crushed in an agate mortar. All samples were then leached with 0.01 M HNO₃ at room temperature for 30 min and rinsed with Milli-Q water to remove superficial contamination and to oxidize any remaining organic matter. Samples are then introduced into the bottom of a two-fingers reactor and 1 cm³ of pure H₃PO₄ is added into the lateral reservoir. The reactor with the wet sample and H₃PO₄ is then connected to the semi-automated vacuum line⁹⁷. The sample is finally dried on the line, and the reactor is manually rotated to pour H₃PO₄ onto the samples. Subsequent steps (reduction for solid source and physical measurement) are either performed by the Laboratoire de Mesure du Carbone 14, Saclay⁹⁸ (LMC14, lab code SacA) on a traditional large AMS or by the ¹⁴C research group of the LSCE, Gif-sur-Yvette (lab codes GifA for chemistry and ECHO for physical measurement) on both solid and gas source of the ECHoMICADAS⁹⁹⁻¹⁰¹. To control the impact of chemical treatment and of the combustion, “blank” and international standards were run. Results are expressed in F¹⁴C as recommended by Reimer et al.¹⁰² and provided as ¹⁴C ages (BP) following Stuiver and Polach (1977)¹⁰³ convention (Supplementary Data 1). Replicates and triplicates from six samples of 5-cm thick in RP1 have been performed to check the reliability between the two methods, the gas source using 10 times less materials (Supplementary Data 1).
**Age-depth models and sedimentation rates.** The chronologies of Schwalbenberg and Nussloch profiles were determined using radiocarbon analyses on ECG (Figs. 2, 3 and Supplementary Data 1). We derived age-depth models using the software BACON 2.2⁵¹ based on Bayesian statistics. For Nussloch, we selected the radiocarbon ages (29 samples) of the profile P8 as age control points prior to running the Bayesian age model. For Schwalbenberg, the model is based on 23 age control points. On average, at Schwalbenberg, the resolution of one radiocarbon sample is about every 20 cm whereas at Nussloch the resolution is about every 60 cm. Probability distributions of calibrated radiocarbon ages were generated using IntCal20¹⁰⁴ in the Bacon code *cc*=1. In Bacon programme⁵¹, two prior distributions (gamma distribution for the accumulation rate and beta distribution for the accumulation memory) are constrained. We used the default values: *accumulation (acc).shape*=1.5 for the gamma distribution; *memory(mem).mean*=0.7 and *mem.strength*=4 for the beta distribution. The mean accumulation rate was changed to 50 yr cm⁻¹ according to information available from calibrated ¹⁴C dates and we run our sequences with 50 and 135 sections respectively for Schwalbenberg and for Nussloch. We choose a resolution of 5 cm for the age-modelling to obtain an individual age for each stable isotope sample. The age-modelling data are compiled in the Supplementary Data 1. The sedimentation rates for both profiles were calculated based on age modelling output at 5 cm resolution.

**Stable isotope analysis of earthworm calcite granules.** We followed the procedure for the extraction of the ECGs and the preparation protocol from refs. 44,45. After cleaning in an ultrasonic bath with deionized water (2 times 5 min), thirty ECGs were extracted from each layer for stable isotope analysis based on the absence of impurities and the well-preserved surface crystals. We performed 180 (30 times 6 layers of 5 cm thickness) and 690 (30 times 23 layers of 5 cm thickness) stable isotope analyses for Nussloch and Schwalbenberg, respectively. Individual granules were crushed, and an aliquot between 200 and 300 µm of carbonate powder was digested with phosphoric acid at 70 °C. The liberated CO₂ gas was analyzed using a Thermo Finnigan GasBench II preparation device interfaced with a Thermo Delta V advantage mass spectrometer at the Max Plank Institute of Chemistry in Mainz. All δ¹³C and δ¹⁸O values are reported on the Vienna-Peecee Belemnite (V-PDB) scale. The long-term standard deviation of a routinely analyzed in-house CaCO₃ standard was <0.1‰ (1σ) for both carbon and oxygen isotope ratios. The CaCO₃ standards are calibrated using NBS19, NBS18 and IAEA 603. The data are reported in Supplementary Data 2.

**Limits and uncertainties for the quantification of the climate parameters and temporal resolution of age models.** The uncertainties of our method are related to sampling, chronology, representativeness and the past precipitation and temperature estimations. The two profiles (Schwalbenberg RP1 and Nussloch P8) have been sampled at 5 cm resolution which represents between 50 and 150 years and the ECGs were extracted from 101 of sediments (~15 kg of sediment). By taking 15 kg of sediment, we reduce the uncertainties associated with the distribution of ECGs inside one sample of 5 cm thickness. Furthermore, ECG concentration is well correlated with the loess palaeosol pedostratigraphy in Northern France⁴⁸ and in the Upper Rhine Valley⁴⁹.

Age-depth modelling is based on 28 and 29 radiocarbon analyses for Schwalbenberg and Nussloch, respectively (Supplementary Data 1). The uncertainties for radiocarbon analysis are in the same order as the maximal counting errors in the NGRIP GICC05 chronology³.

Today earthworms daily produce granules when they are active. We assumed that during the last glacial period, earthworms were active during the warmer half of the year from May to September⁴⁹. The concentrations of the granules (per sediment volume) produced by *Lumbricus terrestris* during the last glacial period is low compared to present-day data⁴⁶. This might be due to either lower productivity of the earthworms because of climate conditions or lower earthworm populations or both. Experiments are still necessary to better constrain parameters that mostly influence the granule production.

We extracted and analyzed 30 individual granules per sample for the stable isotope analyses (one sample is thus made of 30 replicates) to be as representative as possible. To estimate the precipitation, we used a modern dataset between δ¹³C of C₃ plants and MAP so we assumed that this relation did not change during the last glacial period. We calculated and propagated the root squared error through our estimation of the MAP (Fig. 5 and Supplementary Data 2). The distribution of 30 replicates for each sample of 5 cm-thick are sizeable. However, the comparison between the mean values of MIS 3 and MIS 2 are significantly different at both sites (Supplementary Table 1).

To estimate the past temperature, we used a fractionation equation that took into account the vital effect of the earthworms⁵⁰. Based on clumped isotope data⁶⁰, we assumed that the δ¹⁸O value of soil water is not seasonally biased. Soil water evaporation is unfortunately not possible to quantify over past period, thus we did not take into account this parameter in our equation. Moreover, to estimate the δ¹⁸O value of the rain we used modern data and assumed that present-day δ¹⁸O-T gradients are also valid under glacial climate conditions (Supplementary Fig. 3 and Supplementary Data 3). We corrected the data for -0.5 and -0.8‰ according to the method developed by ref. ⁸¹. We reported the distribution of the 30 values for

each sample of 5 cm thick at 1σ. The uncertainties are sizeable (Fig. 5) however, the T-test indicated that the LSTws mean values are significantly different between Schwalbenberg and Nussloch and between MIS 3 and MIS 2 (Supplementary Table 1).

We should be aware that both methods give us a first approximation of climate parameter estimates.

**Wind dynamics model reconstructions.** Wind and precipitation statistics are derived from the instantaneous values of 10-metre winds and total precipitation sampled at 6-hourly (synoptic) intervals from continuous 30 years output of the respective atmospheric general circulation model products. For the present-day conditions, we use the output of the ERA INTERIM reanalysis¹⁰⁵ at 1° by 1° horizontal resolution covering the 1989–2018 period. Note that reanalysis precipitation data are based on forecasted fields. For the LGM conditions, we use the output of the MPI-ESM-P model (simulation r1i1p1-P, approx. 1.88° horizontal resolution) participated in the Climate Model Intercomparison Project, Phase 5 (CMIP5). MPI-ESM-P combines the coupled general circulation models for the atmosphere ECHAM6 and ocean MPIOM and the land surface model JSBACH with vegetation cover prescribed to pre-industrial reconstructions¹⁰⁶. LGM climate simulated in MPI-ESM-P compared well with to date available proxy data¹⁰⁷. Values are sampled in respective model grid cells without interpolation.

**Reporting summary.** Further information on research design is available in the Nature Research Reporting Summary linked to this article.

## Data availability

All data generated or analyzed during this study are included in this published article (and its Supplementary Information files). Supplementary Data files are available at https://doi.org/10.17632/2dsgypk45h1.1. We confirm that the geological materials samples were collected in a responsible manner and in accordance with relevant permits and local laws.

Received: 20 January 2022; Accepted: 19 October 2022; Published online: 21 November 2022

## References

1. Dansgaard, W. et al. Evidence for general instability of past climate from a 250-kyr ice-core record. *Nature* **364**, 218–220 (1993).

2. Bond, G. et al. Correlations between climate records from North Atlantic sediments and Greenland ice. *Nat. Lett.* **363**, 143–147 (1993).

3. Rasmussen, S. O. et al. A stratigraphic framework for abrupt climatic changes during the Last Glacial period based on three synchronized Greenland ice-core records: Refining and extending the INTIMATE event stratigraphy. *Quat. Sci. Rev.* **106**, 14–28 (2014).

4. NGRIP members. High-resolution record of Northern Hemisphere climate extending into the last interglacial period. *Nature* **431**, 147–151 (2004).

5. Heinrich, H. Origin and consequences of cyclic ice rafting in the Northeast Atlantic Ocean during the past 130,000 years. *Quat. Res.* **29**, 142–152 (1988).

6. Martrat, B. et al. Four climate cycles of recurring deep and surface water destabilizations on the Iberian margin. *Science* **317**, 502–507 (2007).

7. Broecker, W. S. Massive iceberg discharges as triggers for global climate change. *Nature* **372**, 421–424, (1994).

8. Menviel, L. C., Skinner, L. C., Tarasov, L. & Tzedakis, P. C. An ice-climate oscillatory framework for Dansgaard–Oeschger cycles. *Nat. Rev. Earth Environ.* **1**, 677–693 (2020).

9. Újvári, G. et al. Coupled European and Greenland last glacial dust activity driven by North Atlantic climate. *Proc. Natl Acad. Sci. USA* https://doi.org/10.1073/pnas.1712651114 (2017).

10. Újvári, G. et al. Stadia-interstadial temperature and aridity variations in east central europe preceding the last glacial maximum. *Paleoceanogr. Paleoclimatol.* **36**, e2020PA004170 (2021).

11. Moine, O. et al. The impact of Last Glacial climate variability in west-European loess revealed by radiocarbon dating of fossil earthworm granules. *Proc. Natl Acad. Sci. USA* **114**, 1–6 (2017).

12. Fischer, P. et al. Millennial-scale terrestrial ecosystem responses to Upper Pleistocene climatic changes: 4D-reconstruction of the Schwalbenberg Loess-Palaeosol-Sequence (Middle Rhine Valley, Germany). *Catena* **196**, 104913 (2021).
13. Zeeden, C. et al. Millennial scale climate oscillations recorded in the Lower Danube loess over the last glacial period. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **509**, 164–181 (2018).

14. Rousseau, D.-D. et al. Link between European and North Atlantic abrupt climate changes over the last glaciation. *Geophys. Res. Lett.* **34**, L22713 (2007).

15. Sirocko, F. et al. The ELSA-Vegetation-Stack: Reconstruction of Landscape Evolution Zones (LEZ) from laminated Eifel maar sediments of the last 60,000 years. *Glob. Planet Change* **142**, 108–135 (2015).

16. Bohncke, S. J. P., Bos, J. A. A., Engels, S., Heiri, O. & Kasse, C. Rapid climatic events as recorded in Middle Weichselian thermokarst lake sediments. *Quat. Sci. Rev.* **27**, 162–174 (2008).

17. Sirocko, F. et al. Muted multidecadal climate variability in central Europe during cold stadial periods. *Nat. Geosci.* **14**, 651–658 (2021).

18. Duprat-Oualid, F. et al. Vegetation response to abrupt climate changes in Western Europe from 45 to 14.7k cal a BP: The Bergsee lacustrine record (Black Forest, Germany). *J. Quat. Sci.* **32**, 1008–1021 (2017).

19. Sánchez-Goñi, M. F. et al. Contrasting impacts of Dansgaard-Oeschger events over a western European latitudinal transect modulated by orbital parameters. *Quat. Sci. Rev.* **27**, 1136–1151 (2008).

20. Fletcher, W. J. et al. Millennial-scale variability during the last glacial in vegetation records from Europe. *Quat. Sci. Rev.* **29**, 2839–2864 (2010).

21. Luetscher, M. et al. North Atlantic storm track changes during the Last Glacial Maximum recorded by Alpine speleothems. *Nat. Commun.* **6**, 6344 (2015).

22. Genty, D. et al. Precise dating of Dansgaard-Oeschger climate oscillations in western Europe from stalagmite data. *Nature* **421**, 833–837 (2003).

23. Fankhauser, A., McDermott, F. & Fleitmann, D. Episodic speleothem deposition tracks the terrestrial impact of millennial-scale last glacial climate variability in SW Ireland. *Quat. Sci. Rev.* **152**, 104–117 (2016).

24. Andrews, J. E. Palaeoclimatic records from stable isotopes in riverine tufas: Synthesis and review. *Earth Sci. Rev.* **75**, 85–104 (2006).

25. Osman, M. B. et al. Globally resolved surface temperatures since the Last Glacial Maximum. *Nature* **599**, 239–244 (2021).

26. Haase, D. et al. Loess in Europe-its spatial distribution based on a European Loess Map, scale 1:2,500,000. *Quat. Sci. Rev.* **26**, 1301–1312 (2007).

27. Fenn, K. & Prud'homme, C. *Treatise on Geomorphology* (eds Shroder, J. F., Lancaster, N., Sherman, D. J. & Baas, A. C. W.) (Academic Press, 2021).

28. Lehmkuhl, F. et al. Loess landscapes of Europe – Mapping, geomorphology, and zonal differentiation. *Earth Sci. Rev.* https://doi.org/10.1016/j.earscirev. 2020.103496 (2020).

29. Antoine, P. et al. Rapid and cyclic aeolian deposition during the Last Glacial in European loess: A high-resolution record from Nussloch, Germany. *Quat. Sci. Rev.* **28**, 2955–2973 (2009).

30. Rousseau, D.-D. et al. Dansgaard-Oeschger-like events of the penultimate climate cycle: The loess point of view. *Clim. Past* **16**, 713–727 (2020).

31. Rousseau, D. D. et al. (MIS3 & 2) millennial oscillations in Greenland dust and Eurasian aeolian records – a paleosol perspective. *Quat. Sci. Rev.* **139**, 99–113 (2017).

32. Ložek, V. Molluscs in loess, their paleoecological significance and role in geochronology - Principles and methods. *Quat. Int.* **7–8**, 71–79 (1990).

33. Rousseau, D.-D., Puisségur, J. J. & Lautridou, J. P. Biogeography of the Pleistocene pleniglacial malacofaunas in Europe. Stratigraphic and climatic implications. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **80**, 7–23 (1990).

34. Moine, O., Rousseau, D.-D. & Antoine, P. The impact of Dansgaard-Oeschger cycles on the loessic environment and malacofauna of Nussloch (Germany) during the Upper Weichselian. *Quat. Res.* **70**, 91–104 (2008).

35. Sümegi, P. et al. High-resolution proxy record of the environmental response to climatic variations during transition MIS3/MIS2 and MIS2 in Central Europe: The loess-paleosol sequence of Katymár brickyard (Hungary). *Quat. Int.* **504**, 40–55 (2019).

36. Sümegi, P. & Krolopp, E. Quatermalacological analyses for modeling of the Upper Weichselian palaeoenvironmental changes in the Carpathian Basin. *Quat. Int.* **91**, 53–63 (2002).

37. Moine, O., Rousseau, D.-D., Jolly, D. & Vianey-Liaud, M. Paleoclimatic reconstruction using mutual climatic range on terrestrial mollusks. *Quat. Res.* **57**, 162–172 (2002).

38. Rousseau, D.-D. Climatic transfer-function from quaternary mollusks in European Loess deposits. *Quat. Res.* **36**, 195–209 (1991).

39. Hatté, C. & Guiot, J. Palaeoprecipitation reconstruction by inverse modelling using the isotopic signal of loess organic matter: Application to the Nußloch loess sequence (Rhine Valley, Germany). *Clim Dyn* **25**, 315–327 (2005).

40. Hatté, C. et al. δ13C of Loess organic matter as a potential proxy for paleoprecipitation. *Quat. Res.* **55**, 33–38 (2001).

41. Zech, R., Gao, L., Tarozo, R. & Huang, Y. Branched glycerol dialkyl glycerol tetraethers in Pleistocene loess-paleosol sequences: Three case studies. *Org. Geochem.* **53**, 38–44 (2012).

42. Trigui, Y. et al. First calibration and application of leaf wax n-alkane biomarkers in loess-paleosol sequences and modern plants and soils in Armenia. *Geosciences (Basel)* **9**, 263 (2019).

43. Zech, M., Rass, S., Buggle, B., Löscher, M. & Zöller, L. Reconstruction of the late Quaternary paleoenvironments of the Nussloch loess paleosol sequence, Germany, using n-alkane biomarkers. *Quat. Res.* **78**, 226–235 (2012).

44. Prud'homme, C. et al. Palaeotemperature reconstruction during the Last Glacial from δ¹⁸O of earthworm calcite granules from Nussloch loess sequence, Germany. *Earth Planet Sci. Lett.* **442**, 13–20 (2016).

45. Prud'homme, C. et al. δ¹³C signal of earthworm calcite granules: A new proxy for palaeoprecipitation reconstructions during the Last Glacial in Western Europe. *Quat. Sci. Rev.* **179**, 158–166 (2018).

46. Canti, M. G. & Piearce, T. G. Morphology and dynamics of calcium carbonate granules produced by different earthworm species. *Pedobiologia (Jena)* **47**, 511–521 (2003).

47. Canti, M. G. Origin of calcium carbonate granules found in buried soils and Quaternary deposits. *Boreas* **27**, 275–288 (1998).

48. Prud'homme, C. et al. Earthworm calcite granules: a new tracker of millennial-timescale environmental changes in Last Glacial loess deposits. *J. Quat. Sci.* **30**, 529–536 (2015).

49. Prud'homme, C., Moine, O., Mathieu, J., Saulnier-Copard, S. & Antoine, P. High-resolution quantification of earthworm calcite granules from western European loess sequences reveals stadial-interstadial climatic variability during the Last Glacial. *Boreas* **48**, 257–268 (2019).

50. Schirmer, W. Rhine loess at Schwalbenberg II – MIS 4 and 3. *E&G Quat. Sci. J.* **61**, 32–47 (2012).

51. Blaauw, M. & Christeny, J. A. Flexible paleoclimate age-depth models using an autoregressive gamma process. *Bayesian Anal.* **6**, 457–474 (2011).

52. Briones, M. J. I., Ostle, N. J. & Piearce, T. G. Stable isotopes reveal that the calciferous gland of earthworms is a CO₂-fixing organ. *Soil Biol. Biochem.* **40**, 554–557 (2008).

53. Canti, M. G., Bronk-Ramsey, C., Hua, Q. & Marshall, P. Chronometry of pedogenic and stratigraphic events from calcite produced by earthworms. *Quat. Geochronol.* **28**, 96–102 (2015).

54. Canti, M. G. Experiments on the origin of ¹³C in the calcium carbonate granules produced by the earthworm *Lumbricus terrestris*. *Soil Biol. Biochem.* **41**, 2588–2592 (2009).

55. Versteegh, E. A. A., Black, S. & Hodson, M. E. Carbon isotope fractionation between amorphous calcium carbonate and calcite in earthworm-produced calcium carbonate. *Appl. Geochem.* **78**, 351–356 (2017).

56. Rey, K. et al. Late Miocene climatic and environmental variations in northern Greece inferred from stable isotope compositions (delta O-18, delta C-13) of equid teeth apatite. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **388**, 48–57 (2013).

57. Kohn, M. J. Carbon isotope compositions of terrestrial C3 plants as indicators of (paleo)ecology and (paleo)climate. *Proc. Natl Acad. Sci.* **107**, 19691–19695 (2010).

58. Hare, V. J., Loftus, E., Jeffrey, A. & Ramsey, C. B. Atmospheric CO₂ effect on stable carbon isotope composition of terrestrial fossil archives. *Nat. Commun.* **9**, 252 (2018).

59. Versteegh, E. A. A., Black, S., Canti, M. G. & Hodson, M. E. Earthworm-produced calcite granules: A new terrestrial palaeothermometer? *Geochim. Cosmochim. Acta* **123**, 351–357 (2013).

60. Kelson, J. R. et al. A proxy for all seasons? A synthesis of clumped isotope data from Holocene soil carbonates. *Quat. Sci. Rev.* **234**, 106259 (2020).

61. Tan, H., Liu, Z., Rao, W., Jin, B. & Zhang, Y. Understanding recharge in soil-groundwater systems in high loess hills on the Loess Plateau using isotopic data. *Catena* **156**, 18–29 (2017).

62. Dansgaard, W. Stable isotopes in precipitation. *Telltus* **16**, 436–468 (1964).

63. Rozanski, K. Deuterium and oxygen-18 in European groundwaters-links to atmospheric circulation in the past. *Chem. Geol.* **52**, 349–363 (1985).

64. Zhang, T., Osterkamp, T. E. & Stamnes, K. Effects of climate on the active layer and permafrost on the North slope of Alaska, U.S.A. *Permafr. Periglac. Process* **8**, 45–67 (1997).

65. Zhang, T. Influence of seasonal snow cover on the ground thermal regime: An overview. *Rev. Geophys.* **43**, RG4002 (2005).

66. Christiansen, H. H. et al. The thermal state of permafrost in the nordic area during the international polar year 2007–2009. *Permafr. Periglac. Process* **21**, 156–181 (2010).

67. Spratt, R. M. & Lisiecki, L. E. A Late Pleistocene sea level stack. *Clim. Past* **12**, 1079–1092 (2016).

68. Toucanne, S. et al. Timing of massive 'Fleuve Manche' discharges over the last 350 kyr: Insights into the European ice-sheet oscillations and the European drainage network from MIS 10 to 2. *Quat. Sci. Rev.* **28**, 1238–1256 (2009).

69. Schrag, D. P. et al. The oxygen isotopic composition of seawater during the Last Glacial Maximum. *Quat. Sci. Rev.* **21**, 331–342 (2002).

70. Lambeck, K. & Chappell, J. Sea level change through the last glacial cycle. *Science* **292**, 679–686 (2001).
71. Shackleton, N. J. The 100,000-year ice-age cycle identified and found to lag temperature, carbon dioxide, and orbital eccentricity. *Science* **289**, 1897–1902 (2000).

72. Vinnepand, M. et al. Combining inorganic and organic carbon stable isotope signatures in the Schwalbenberg Loess-Paleosol-sequence near Remagen (Middle Rhine Valley, Germany). *Front. Earth Sci.* https://doi.org/10.3389/feart.2020.00276 (2020).

73. Schiermeyer, J. Würmzeitliche Lößmollusken aus der Eifel (Universität Düsseldorf, 2000).

74. Hatté, C., Rousseau, D.-D. & Guiot, J. Climate reconstruction from pollen and δ¹³C using inverse vegetation modeling. Implication for past and future climates. *Clim. Past Discuss.* 5, 73–98 (2009).

75. Kaspar, F. et al. Monitoring of climate change in Germany – data, products and services of Germany’s National Climate Data Centre. *Adv. Sci. Res.* **10**, 99–106 (2013).

76. Engels, S. et al. Environmental inferences and chironomid-based temperature reconstructions from fragmentary records of the Weichselian Early Glacial and Pleniglacial periods in the Niederlausitz area (eastern Germany). *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **260**, 405–416 (2008).

77. Ampel, L. et al. Modest summer temperature variability during DO cycles in western Europe. *Quat. Sci. Rev.* **29**, 1322–1327 (2010).

78. Obreht, I. et al. The Late Pleistocene Beloturn section (southern Serbia) at the southern limit of the European loess belt: Environmental and climate reconstruction using grain size and stable C and N isotopes. *Quat. Int.* **334–335**, 10–19 (2014).

79. Schreuder, L. T., Beets, C. J., Prins, M. A., Hatté, C. & Peterse, F. Late Pleistocene climate evolution in Southeastern Europe recorded by soil bacterial membrane lipids in Serbian loess. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **449**, 141–148 (2016).

80. Zamanian, K., Lechler, A. R., Schauer, A. J., Kuzyakov, Y. & Huntington, K. W. The δ13C, δ18O and Δ 47 records in biogenic, pedogenic and geogenic carbonate types from paleosol-loess sequence and their paleoenvironmental meaning. *Quat. Res.* **101**, 256–272 (2021).

81. Lécuyer, C., Hillaire-Marcel, C., Burke, A., Julien, M. A. & Hélie, J. F. Temperature and precipitation regime in LGM human refugia of southwestern Europe inferred from δ13C and δ18O of large mammal remains. *Quat. Sci. Rev.* **255**, 106796 (2021).

82. Moine, O., Rousseau, D.-D. & Antoine, P. Terrestrial molluscan records of Weichselian Lower to Middle Pleniglacial climatic changes from the Nussloch loess series (Rhine Valley, Germany): The impact of local factors. *Boreas* **34**, 363–380 (2005).

83. Seltzer, A. M. et al. Widespread six degrees Celsius cooling on land during the Last Glacial Maximum. *Nature* **593**, 228–232 (2021).

84. Kindler, P. et al. Temperature reconstruction from 10 to 120 kyr b2k from the NGRIP ice core. *Clim. Past* **10**, 887–902 (2014).

85. Ludwig, P., Schaffernicht, E. J., Shao, Y. & Pinto, J. G. Regional atmospheric circulation over Europe during the last Glacial maximum and its links to precipitation. *J Geophys. Res.* **121**, 2130–2145 (2016).

86. Moine, O. Weichselian Upper Pleniglacial environmental variability in north-western Europe reconstructed from terrestrial mollusc faunas and its relationship with the presence/absence of human settlements. *Quat. Int.* **337**, 90–113 (2014).

87. Hughes, A. L. C., Gyllencreutz, R., Lohne, Ø. S., Mangerud, J. & Svendsen, J. I. The last Eurasian ice sheets - a chronological database and time-slice reconstruction, DATED-1. *Boreas* **45**, 1–45 (2016).

88. Toucanne, S. et al. The North Atlantic glacial eastern boundary current as a key driver for ice-sheet—AMOC interactions and climate instability. *Paleoceanogr. Paleoclimatol.* **36**, e2020PA004068 (2021).

89. Batchelor, C. L. et al. The configuration of Northern Hemisphere ice sheets through the Quaternary. *Nat. Commun.* **10**, 3713 (2019).

90. Liakka, J., Löfverström, M. & Colleoni, F. The impact of the North American glacial topography on the evolution of the Eurasian ice sheet over the last glacial cycle. *Clim. Past* **12**, 1225–1241 (2016).

91. Spötl, C., Koltai, G., Jarosch, A. H. & Cheng, H. Increased autumn and winter precipitation during the Last Glacial Maximum in the European Alps. *Nat. Commun.* **12**, 1839 (2021).

92. Stokes, C. R., Tarasov, L. & Dyke, A. S. Dynamics of the North American ice sheet complex during its inception and build-up to the last glacial maximum. *Quat. Sci. Rev.* **50**, 86–104 (2012).

93. Li, C. & Born, A. Coupled atmosphere-ice-ocean dynamics in Dansgaard-Oeschger events. *Quat. Sci. Rev.* **203**, 1–20 (2019).

94. Vettoretti, G. & Peltier, W. R. Last Glacial Maximum ice sheet impacts on North Atlantic climate variability: The importance of the sea ice lid. *Geophys. Res. Lett.* **40**, 6378–6383 (2013).

95. Schaffernicht, E. J., Ludwig, P. & Shao, Y. Linkage between dust cycle and loess of the Last Glacial Maximum in Europe. *Atmos. Chem. Phys.* **20**, 4969–4986 (2020).

96. App, V., Auffermann, B., Hahn, J., Pasda, C. & Stephan, E. Die altsteinzeitliche Fundstelle auf dem Schwalbenberg bei Remagen. *Berichte zur Archäologie an Mittelrhein und Mosel* **4**, 11–136 (1995).

97. Tisnérat-Laborde, N., Poupeau, J., Tannau, J. & Paterne, M. Development of a semi-automated system for routine preparation of carbonate samples. *Radiocarbon* **43**, 299–304 (2001).

98. Cottereau, E. et al. Artemis, the new 14C AMS at LMC14 in Saclay, France. *Radiocarbon* **49**, 291–299 (2007).

99. Ruff, M. et al. A gas ion source for radiocarbon measurements at 200kV. *Radiocarbon* **49**, 307–314 (2007).

100. Tisnérat-Laborde, N. et al. ECHoMICADAS: A new compact AMS system to measuring 14C for Environment, Climate and Human Sciences. In *International Radiocarbon Conference*, Dakar, Senegal 16–20 (2015).

101. Synal, H., Stocker, M. & Suter, M. MICADAS: A new compact radiocarbon AMS system. *Nucl. Instrum. Methods Phys. Res. B* **259**, 7–13 (2007).

102. Reimer, P. J., Brown, T. A. & Reimer, R. W. Discussion: Reporting and calibration of post-bomb 14C data. *Radiocarbon* **46**, 1299–1304 (2004).

103. Stuiver, M. & Henry, Polach. Discussion reporting of 14C Data. *Radiocarbon* **19**, 355–363 (1977).

104. Reimer, P. J. et al. The IntCal20 northern hemisphere radiocarbon age calibration curve (0–55cal kBP). *Radiocarbon* **62**, 725–757 (2020).

105. Dee, D. P. et al. The ERA-Interim reanalysis: Configuration and performance of the data assimilation system. *Quart. J. R. Meteorol. Soc.* **137**, 553–597 (2011).

106. Pongratz, J., Reick, C. H., Raddatz, T. & Claussen, M. Effects of anthropogenic land cover change on the carbon cycle of the last millennium. *Glob. Biogeochem. Cycles* **23**, GB4001 https://doi.org/10.1029/2009GB003488 (2009).

107. Hargreaves, J. C., Annan, J. D., Ohgaito, R., Paul, A. & Abe-Ouchi, A. Skill and reliability of climate model ensembles at the Last Glacial Maximum and mid-Holocene. *Clim. Past* **9**, 811–823 (2013).

108. Berger, A. & Loutre, M. F. Insolation values for the climate of the last 10 million years. *Quat. Sci. Rev.* **10**, 297–317 (1991).

109. Toucanne, S. et al. Millennial-scale fluctuations of the European Ice Sheet at the end of the last glacial, and their potential impact on global climate. *Quat. Sci. Rev.* **123**, 113–133 (2015).

## Acknowledgements

This work was supported by the German Science Foundation (Deutsche Forschungsgemeinschaft, DFG) in the framework of the TERRACIME-Project (FI 1941/5-1; FI 1918/4-1; VO 938/25-1) awarded to P.F., K.E.F., and A.V., and by an independent Max Planck Research Group awarded to K.E.F. S.G. acknowledges support through the project PalMod, funded by the German Federal Ministry of Education and Research (BMBF). We thank Sven Brömme for help in the stable isotope laboratory and Felix Hettwer for helping to extract the earthworm calcite granules. We thank everyone who assisted with fieldwork in Remagen, including Stefania Milano and Kristina Reetz. We are grateful to the three reviewers for their comments that helped us to improve the scientific content of this work.

## Author contributions

C.P., P.F., K.E.F., and A.V. designed the study and C.P., P.F., O.J., M.V., and O.M. conducted the fieldwork. C.P. and H.V. conducted the stable isotope analyses C.P. prepared the samples prior radiocarbon dating and C.H. performed the analyses S.G. performed the analysis of wind and precipitation regime reconstructions. All authors contributed to discussion, interpretation of the results and writing of the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-022-00595-3.

**Correspondence** and requests for materials should be addressed to Charlotte Prud'homme.

**Peer review information** *Communications Earth & Environment* thanks Thomas Stevens and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Ola Kwiecien and Clare Davis.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
<figure><img src="image_10.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022