Tracking westerly wind directions over Europe since the middle Holocene

The variability of the northern westerlies has been considered as one of the key elements for modern and past climate evolution. Their multiscale behavior and underlying control mechanisms, however, are incompletely understood, owing to the complex dynamics of Atlantic sea-level pressures. Here, we present a multi-annually resolved record of the westerly drift over the past 6,500 years from northern Italy. In combination with more than 20 other westerly-sensitive records, our results depict the non-stationary westerly-affected regions over mainland Europe on multi-decadal to multi-centennial time scales, showing that the direction of the westerlies has changed with respect to the migrations of the North Atlantic centers of action since the middle Holocene. Our findings suggest the crucial role of the migrations of the North Atlantic dipole in modulating the westerly-affected domain over Europe, possibly modulated by Atlantic Ocean variability.

The Mediterranean Basin has been a cradle of civilizations since the middle Holocene (8.2-4.2 thousand years before 1950 C.E., kyr BP). Today some 400 million people live in this region. As a "climate hot spot," the Mediterranean features large hydroclimate variability in response to global climate change 1 and has been experiencing exceptionally low rainfall over the past two decades 2 , resulting in significant impacts on ecosystems, human society, and the economy 3 . The drought has been attributed to an enhanced sea-level pressure (SLP) contrast between the Azores High and the Icelandic Low (the SLP dipole, hereafter; Supplementary Fig. 1), known as the positive phase of the North Atlantic Oscillation (NAO) 4 , which led to a northerly migration of the westerlies, transporting moisture away from the Mediterranean region ( Supplementary Fig. 2a). The location of the SLP dipole is not stable through time. Its spatial displacement determines the angle of the westerly tracks and the phase of the East Atlantic (EA) pattern [5][6][7][8][9][10][11] , defined by a pressure anomaly centered over the eastern North Atlantic (52°30'N, 27°30'W; Supplementary Fig. 2b) 12 . The variabilities of the SLP dipole and varying westerly tracks can therefore result in unstable NAO-correlated rainfall regimes (red and blue areas in Supplementary Movie 1) in Europe on decadal to multidecadal timescales, termed "non-stationary NAO" behavior 13 .
Climate proxy records are essential to develop an in-depth understanding of the westerlies on a range of timescales. For example, terrestrial and marine sediments in Iberia 14 , and Moroccan tree-ring data combined with Scottish speleothem records 15 suggest a prolonged positive NAO phase during the Medieval Climate Anomaly (1050-850 yr BP). However, Northern Hemisphere proxy assemblages and model simulations argued that proxy-based paleo-NAO reconstructions are sometimes contentious [16][17][18] , largely because of the complex interactions among Atlantic SLP. Only few proxy records reliably reflect Atlantic SLP patterns, rendering tracking the westerlies' position across mainland Europe challenging 19 . Since NAO-correlated rainfall regimes (Supplementary Movie 1) reflect the position of the westerlies [5][6][7][8][9][10][11] , past positions of the westerlies can be constrained by comparing a series of westerly-sensitive records in the North Atlantic region.
Here, we present a multi-annually resolved precipitation record from such a sensitive region (northern Italy). In conjunction with other regional westerly-sensitive archives, these data depict large-scale hydroclimate changes and the variability of the SLP dipole and the westerlies over the past 6500 years.  21 and cave records from central Italy ( Supplementary Fig. 9b, c) 22,23 are -0.42 (n = 31, p < 0.05) and 0.67 (n = 31, p < 0.05), respectively (Supplementary Table 2). The drought between 5 and 4 kyr BP was concurrent with dry periods recorded by speleothems from Turkey 27 , Lebanon 28 , Israel 29 , and Morocco 30,31 (Supplementary Fig. 10a-d). The following wet period during 4-3 kyr BP is in good agreement with records from circum-Mediterranean countries (Supplementary Figs. 9d-f, 10a, b, and 11c, d) [24][25][26][27][28]32,33 . These events likely had a profound impact on ancient Mediterranean human societies ( Supplementary Fig. 12). For example, the 5.2 kyr BP dry period has been suggested to have resulted in the demise of the Uruk period 29 , while the dry period at 4.8-4.1 kyr BP (4.2kyr event) possibly induced the collapse of the Old Kingdom 34 and the Akkadian Empire 35 .

Results and discussion
Hydroclimate changes in northern Italy are strongly affected by North Atlantic SLP variability ( Supplementary Fig. 2a, b) and westerly dynamics. Instrumental precipitation data from Genoa (Genoa PP) are positively/negatively correlated with the 850-mb zonal wind around Gibraltar/Iceland during the rainy season ( Supplementary Fig. 2c), suggesting that changes in the position of the westerlies affect Genoa PP. Genoa PP show a strong negative correlation with SLP variations centered over northwestern Europe ( Supplementary Fig. 2d), resembling the EA pattern 12 , implying that spatial changes of the SLP dipole dominate regional rainfall patterns. The SLP in northwestern Europe (averaged over [45-55°N, 20°W-0]) is also significantly anti-correlated with the 850-mb zonal wind around Toirano (averaged over [42-47°N, 5-15°E]), with a correlation coefficient of −0.6 ± 0.04 (n = 97; p < 0.01; 1920-2016 C.E.) during winter (December-February, DJF) based on NCAR/NCEP Reanalysis v3. This anticorrelation emphasizes that the westerlies dominate the hydroclimate in Toirano. The leading empirical orthogonal function (EOF) 1 of European-Atlantic (60°W-40°E, 20-80°N) precipitation represents an EA teleconnection ( Supplementary Fig. 13a) that is positively correlated with Genoa PP (r = 0.43 ± 0.04, n = 173, p < 0.01). Combined with the NAO pattern observed in EOF2 ( Supplementary Fig. 13b), our results show that the variability of the SLP dipole influences Toirano precipitation (Supplementary Fig. 2a, b). The precipitation patterns documented by the Toirano records hence reflect movements of the westerlies under the two climate modes, consistent with previous studies [36][37][38] .
On multi-decadal to multi-centennial time scales, Bàsura stalagmite-inferred paleo-precipitation is positively correlated (r = 0.69, n = 292, 99%, Methods) with a NAO reconstruction 16 of the past millennium (1049-1969 C.E.; Supplementary Fig. 14). In contrast to the negative correlation of modern Genoa PP with the NAO index 39 in winter (DJF; Supplementary Fig. 2a), the comparison suggests that Toirano was in a positive NAO-correlated region over the past millennium, where the westerlies were strong and precipitation was high during the positive NAO phase. When compared with lacustrine sediment records from western Greenland 40 ( Fig. 1d) that reflect NAOdriven local temperature changes and westerly dynamics, Bàsura Δ 18 O inferred precipitation patterns also show a positive correlation in the interval 2.2-1.2 kyr BP (r = 0.37, n = 64, 95%) and 0.8-0.3 kyr BP (r = 0.69, n = 41, 95%), but a negative correlation (r = -0.48, n = 84, 99%) between 4.2 and 2.2 kyr BP (Fig. 1d). In other words, the relationship of Toirano precipitation with the reconstructed NAO index 16 remained positive after 2.2 kyr BP, but was negative in the interval 4.2-2.2 kyr BP, revealing that the NAO-correlated regions over Europe have changed since the middle Holocene.

Non-stationary North Atlantic Oscillation
Over the past 150 years, correlations of Toirano precipitation with western Greenland temperature ( Fig. 1c; black line) and with the NAO index 39 ( Fig. 1c; blue line) have been unstable, reflecting the nonstationary NAO behavior (Supplementary Movie 1). Correlations between NAO index and Toirano precipitation were, for example, more significant for 1870-1900 and 1950-1980 C.E. than for 1905-1935 C.E. (Fig. 1c), suggesting that Toirano was located in negative NAOcorrelated regions in the late 19 th and late 20 th centuries ( Fig. 1a), but close to the boundary of (or even temporally in positive) NAOcorrelated regions in the early 20 th century (Fig. 1b). Such changes in NAO-correlated regions have been suggested to be modulated by the migrations of the SLP dipole (i.e., EA phases) [5][6][7][8][9][10][11] . Northward or counterclockwise rotated displacements of the SLP dipole could lead the westerly tracks towards a southwest-northeast tilt and shift the positive NAO-correlated regions towards higher latitudes 5,9 (e.g., Fig. 1a).
Variable correlations between Toirano precipitation and hydroclimate in the North Atlantic ( Fig. 1c; black and blue lines) match the trend of the SLP time series over southwestern England ( Fig. 1c; brown line) and the reconstructed EA index ( Fig. 1c; red line) 41 . This confirms that EA phase changes could affect the locations of the SLP dipole and hence the NAO-correlation regions over Europe, in accordance with previous studies 9 . For longer timescales, isotope-enabled general circulation models indicate that the multi-decadal correlation patterns of NAO-precipitation δ 18 O in Europe can be affected by the EA phases 42 . The Bàsura Δ 18 O also shows similarity with the 700-mb height pressure variation in the eastern Atlantic on centennial to millennial scales obtained from a transient simulation (TraCE-21ka) with all forcings using the Community Climate System Model version 3 (Fig. 2e). Accordingly, correlation changes between the NAO index and Bàsura Δ 18 O (Fig. 1d) over the past thousands of years can be attributed to the variability in the position of the SLP dipole, which modulates the NAOcorrelation regions.
The observation of the changing NAO-correlated regions demonstrates a centennial to millennial mixing effect of NAO and EA modes on European hydroclimate changes. The NAO mode in principle dominated the "dipole" precipitation pattern over mainland Europe (Fig. 1a, b) over the past 6500 years. The EA mode also played an important role in modulating the centennial to millennial positions of the NAO dipole. This phenomenon confirms the combined effect of NAO and EA on positions of the NAO dipole over the instrumental period 8,9 . For example, a concurrent positive NAO phase and negative EA phase during the period 5.4-3.5 kyr BP (Fig. 2g) forced the Azores High northeastward. Coinciding positive NAO and EA phases from 2.2-1.2 kyr BP, on the other hand, resulted in a southwestward shift of the Azores High.

Migration of the westerlies since the middle Holocene
SLP dipole movements could alter the route of the westerlies over the East Atlantic [5][6][7][8][9][10][11] . A northward displacement of the Azores High, such as prior to 2.2 kyr BP, favors a SW-NE direction of the westerlies over Europe. In this case, the pivot point of the SLP dipole (the NAO-zero correlation line) can also rotate into a SW-NE direction, restricting the positive NAO-correlated regions to northwestern Europe (Fig. 2g). After 2.2 kyr BP, the SLP dipole migrated southwards and the westerlies tilted less over Europe. The positive NAO-correlated regions shifted southward and northern Italy was thus positively correlated with the reconstructed NAO index (Fig. 2h).
A high degree of similarity between our Bàsura Δ 18 O record and a stalagmite δ 18 O series from Kaite cave, Spain (Supplementary Fig. 9d)   Their colors indicate wet/warm (red), cold/dry (blue), and conditions difficult to categorize (gray). Maps in g and h were generated using Ocean Data View. 1: Lake SS1220 40 Fig. 9d), therefore, suggest the co-occurrence of westward shifts of the Icelandic Low with eastward shifts of the Azores High. The induced counterclockwise rotation of the SLP dipole in turn can lead to a change in the NAO-correlated regions over Europe 8,9 . Combined with studies that suggest more frequent positive NAO phases from the middle to the late Holocene 40,56,57 , our results reveal the complex interaction between the NAO state and its dipole locations on centennial to millennial scales, with associated changes in the orientation of westerly routes across Europe.

Possible effect of latitudinal oceanic thermal gradient on the sea-level pressure dipole
The movements of the SLP dipole are likely modulated by the latitudinal gradient of sea-surface temperatures (SSTs) in the North Atlantic. On decadal to centennial scales, observational reanalysis and model simulations have shown that the multidecadal SLP response to Atlantic Multidecadal Variability (AMV) projects on migrations of the SLP dipole (i.e., changes in EA phases) 58,59 . Warm AMV phases that feature a warm North Atlantic Ocean correspond to a southward shift of the SLP dipole (positive EA). Over the past 2.2 kyr, SSTs have increased in northern Norway 60 (K23258, Fig. 3b) and eastern Greenland 61,62 (MD99-2269, Fig. 3d; MD99-2322, Fig. 3e), in line with elevated sea-surface salinity (SSS) levels around the Subpolar Gyre 63,64 ( Fig. 3a; RAPiD-12-1K, Fig. 3f; MD99-2227, Fig. 3g) and an increased proportion of Atlantic water-sensitive coccoliths in the Nordic Sea 65 (MD95-2011, Fig. 3c). These observations reveal that, over the past 2.2 kyr, more warm, salty, low-latitude Atlantic water has been delivered towards high latitudes compared to before, possibly reducing the latitudinal SST gradient and hence leading to a southward migration of the SLP dipole. The linkage of AMV and EA is also in accord with results of the preindustrial controlled Community Earth System Model 2 (CESM2) 66 (Methods), which shows a strong coherence between SLP anomalies over the East Atlantic and SST changes along the North Atlantic Current (Fig. 3a), suggesting that Atlantic Ocean dynamics play an important role in the migration of the SLP.
Our Δ 18 O series in combination with other westerly-sensitive records show multiple patterns of European westerly drift on decadal to millennial scales since the middle Holocene in response to the nonstationary behavior of the NAO that was possibly modulated by Atlantic sea-surface conditions. Our results underscore the impact of changing Atlantic Ocean circulations on the position of the westerlies over Europe that in the near future might be associated with variability of the meridional overturning circulation, increasing greenhouse gases, and/or varied aerosol concentrations 59 .

Methods
Stalagmite samples, U-Th dating, and age model Two stalagmites, 190-mm-long BA14-1 and 90 mm-long BA18-4, were collected in two cave chambers of narrow side passages of Bàsura cave, 500 and 800 m from the entrance (Supplementary Fig. 3a), in January 2014 and June 2018, respectively. The stalagmites were cut into halves and polished (Supplementary Fig. 4). High resolution scanning electron microscope (SEM) analysis shows that BA14-1 formed as needleshaped aragonite and BA18-4 as rhombic calcite crystals (Supplementary Fig. 4).
A total of 78 (BA14-1) and 18 (BA18-4) subsamples, 10-100 mg each, were drilled for U-Th chemistry 67 and dating 67,68 (Supplementary Data 1) at the High-Precision Mass Spectrometry and Environment Change Laboratory (HISPEC), Department of Geosciences, National Taiwan University. All U-Th isotopic measurements were conducted on a multi-collector inductively coupled plasma mass spectrometer, Thermo-Finnigan Neptune 68 . A gravimetrically calibrated 69 triplespike, 229 Th-233 U-236 U, and the isotope dilution method were employed to correct for mass bias and to determine U-Th isotopic compositions and contents. Half-lives of U-Th nuclides used for U-Th age calculation are given in ref. 69. Uncertainties in isotopic data and dates, relative to 1950 C.E., are given at the two-sigma (2σ) level or two standard deviations of the mean (2σ m ) unless otherwise noted. U-Th contents, isotopic compositions and ages are given in Supplementary Data 1. Age corrections for the initial 230 Th are 0-2 years, smaller than dating errors of ± 2-30 years. A Monte-Carlo-derived age-depth model was constructed using StalAge 70 techniques ( Supplementary Fig. 5).

Carbonate Sr/Ca
Powdered subsamples, 10-50 μg each, were micro-milled at 0.1-0.5 mm intervals along the growth axis for Sr/Ca determination on BA14-1 and BA18-4, with external matrix-matched standards for every 4-5 samples on an inductively coupled plasma sector-field mass spectrometer, Finnigan Element II, at the HISPEC, National Taiwan University 72 . The 2-sigma reproducibility is ± 0.5%.
BA14-1 Sr/Ca, with a mean of 1.2 (± 0.2, 1σ) mmol/mol, vary from 0.45 to 1.9 mmol/mol. BA18-4 Sr/Ca range from 0.035 to 0.25 mmol/ mol and the mean is 0.048 (± 0.017, 1σ) mmol/mol. The offset between the two Sr/Ca datasets is caused by different partition coefficients, 0.1-0.2 for calcite 73 and 0.8-2.0 for aragonite 74 . To combine the two Sr/Ca series, BA14-1 and BA18-4 Sr/Ca data were first converted to z-standard records using the mean and standard deviation of each dataset. The converted z-standard records (Supplementary Data 3) were zigzagged in the overlapping interval, 728-875 yr BP, based on the assumption that the precipitation changes influenced Sr/Ca of aragonite (BA14-1) and calcite (BA18-4) in the same direction (Supplementary Text 3).

Error propagation of correlation analysis
For correlations between data without age uncertainties, the upper/ lower limits for the correlation coefficient were added, derived from the probable errors 75 (P.E.) in correlations with P.E. = 0.6745 × (1−r 2 )/ N 0.5 , where r is the correlation coefficient and N is the number of observation pairs.
For the correlation of two time series with age uncertainties, the method of   76 was applied. This method embraces a Monte Carlo approach to simulate the best correlation between two series. The maximum correlation coefficient between two real series is presented in the text followed by the significance level. Specifically, more than two thousand artificial random time series that have the same characteristics as the original time series (e.g., variance, auto-correlation coefficients, data resolution and absolute ages with age uncertainties) were generated. By tuning these artificially generated time series, the distribution of the best correlations and the significant level can be obtained. For example, if 5% of tuned artificial series yield a correlation coefficient above 0.5, the best correlation of real series with a coefficient of > 0.5 is significant at the 95% significance level. The maximum correlation coefficient between two series is presented in the text followed by the significance level.
For Supplementary Table 2, the correlation analysis was conducted on centennial to multi-centennial scales. We first resampled each of the 20 proxy records, using 15-year time windows and averaged these values for consecutive 150-year periods, resulting in time series of X values (degrees of freedom = X-1), each representing a 150-year average. We then calculated correlation coefficients between each of these time series and the transformed Bàsura record.

Community Earth System Model Version 2
A Community Earth System Model Version 2 (CESM2) preindustrial simulation 66 was performed at approximately 2°× 2°resolution to represent natural climate variability. To verify that, the proposed processing holds on multidecadal to centennial time scales as in our records, a 5-yr running mean on this long preindustrial simulation was performed in a 300-year window. The sea-surface temperature data were correlated with the 700-mb height pressure variation in the eastern Atlantic (52˚30'N, 27˚30'W).