Multidecadal behavior of the North Atlantic Oscillation during the last millennium

The North Atlantic Oscillation (NAO) is a major source of atmospheric variability in the Northern Hemisphere, affecting temperature, precipitation, and storm tracks across North America and Eurasia. Understanding NAO variability on multidecadal to centennial timescales requires paleo-reconstructions, but previously published reconstructions disagree on the magnitude of low-frequency NAO variability over the last millennium. Paleoclimate proxies for the oxygen and hydrogen isotope composition of meteoric waters have thus far been under-utilized in published NAO reconstructions. Here, we investigate multidecadal variability in the reconstructed NAO over the last millennium using 94 high-resolution NAO-sensitive records from the Iso2k database, a collection of globally distributed water isotope-based paleoclimate proxies. We find significant multidecadal to centennial scale variability, which we also highlight in other independent reconstructions of the NAO. Critically, however, the strength of the low-frequency signal has not been consistent throughout the last millennium. Isotope-enabled model simulations did not reproduce the low-frequency signal in the NAO reconstructions and thus it may be necessary to account for low-frequency variability when projecting the impacts of the NAO on temperature and precipitation under future climate scenarios.

1 Introduction

Weather and climate in Europe, Greenland, and eastern North America are heavily influenced by the North Atlantic Oscillation (NAO), which is associated with the dominant pattern of sea level pressure variability in the North Atlantic. The NAO is responsible for variability and extremes in both temperature and hydroclimate across the region (Casanueva et al., 2014; Hurrell, 1995; Hurrell and Deser, 2009; Trigo et al., 2002; Villarini et al., 2011; Zanardo et al., 2019), so better understanding the NAO and its connections to temperature and precipitation is critical for contextualizing past and future climate changes.

The NAO index is typically defined as the difference between the normalized winter sea-level pressures in Stykkishólmur, Iceland, and Gibraltar (Jones et al., 1997; Vinther et al., 2003a). The index is a proxy for coherent large-scale patterns of atmospheric circulation. For instance, during positive NAO conditions the North Atlantic Jet is deflected northwards, driving warm and wet winters over northern Europe and eastern North America contrasted by cool and dry winters in the Mediterranean. The opposite patterns are observed during negative NAO conditions (Hurrell, 1995; Hurrell and Deser, 2009; Woollings et al., 2010, 2018; Woollings and Blackburn, 2012).

Europe is home to some of the longest continuous instrumental records of sea level pressure, extending back through 1823 CE (Jones et al., 1997; Vinther et al., 2003a). However, even these long records provide only a limited perspective on the low frequency (multidecadal to centennial timescale) variability of the NAO, timescales that are critical for contextualizing and understanding its response to anthropogenic emissions (Cook et al., 2019; Woollings et al., 2015). Additionally, the instrumental records underlying indices for the NAO are focused on Europe and the eastern North Atlantic, and thus are more limited in capturing NAO impacts further up- (e.g., North America and Greenland) and down-stream (e.g., Eurasia) (e.g. Felis et al., 2000; Liu et al., 2018; Xiaoge et al., 2010). Instrumental data for the NAO outside of the North Atlantic region is especially limited in duration and coverage, so paleoclimate proxy-based reconstructions are a valuable tool for characterizing the NAO and its impacts on timescales that extend beyond the instrumental record, and outside the regions with long instrumental records of sea level pressure.

Substantial disagreements still exist between published reconstructions of the NAO, especially on the nature of multidecadal to centennial scale variability during the last millennium (Cook et al., 2019; Schmutz et al., 2000). Geologic archives such as tree-ring width chronologies, speleothems, glacier ice, and lake sediments have been employed to investigate paleo-signals of the NAO (Baldini et al., 2008; Brittingham et al., 2019; Cook et al., 2002; Hernández et al., 2020; Kozachek et al., 2017; Kuhnert et al., 2005; Lehner et al., 2012; Michel et al., 2020; Ortega et al., 2015; Sánchez-López et al., 2016; Sorrel et al., 2007; Trouet et al., 2009; Vinther et al., 2003b). However, these published reconstructions suggest a wide range of NAO variability across timescales from higher relative variance at low frequencies (e.g., Ortega et al., 2015) to nearly equal variance at all frequencies (e.g., Cook et al., 2019). Properly constraining low-frequency variability in the NAO is critical not only for future projections, but also for establishing the mechanisms and drivers of the NAO across timescales.

Understanding discrepancies in the published reconstructions is complicated by the different proxy types used as predictors, which may bring their own inherent spectral biases. Ice cores, tree rings, and speleothems, for instance, have each been shown to exhibit unique spectral biases (Dee et al., 2017; Franke et al., 2013). Drought-sensitive tree-rings, which have been used extensively in reconstructions of the NAO, can exhibit substantial autocorrelation (Alley, 1984; Ault et al., 2013). Only a small fraction of the proxy records included in published NAO reconstructions are based on the stable oxygen or hydrogen isotopic compositions of environmental waters like precipitation, seawater, lake water, or soil and groundwater (hereafter called water isotopes). Water isotopes provide information about the NAO on broader spatial and temporal scales than other proxies because they integrate information on basin-wide to hemispheric scales that is otherwise overprinted by regional or site-specific noise (Moerman et al., 2013). Water isotope-based proxies therefore offer a valuable new perspective on the NAO during the last millennium.

Water isotopes have long been used as integrative tracers of the modern water cycle because they are modified by fractionation processes in which the rare heavy isotopologues of water (e.g., ¹H${}_{\mathrm{2}}^{\mathrm{18}}$O, ¹H²H¹⁶O) fractionate from their lighter, more common counterpart (¹H${}_{\mathrm{2}}^{\mathrm{16}}$O) during evaporation, condensation, and other phase changes (Bowen et al., 2019; Dansgaard, 1964; Galewsky et al., 2016; Rozanski et al., 1993). The ratios of the heavy to light isotopes, relative to Standard Mean Ocean Water, are represented by δ¹⁸O for oxygen and δD for hydrogen. We focus on δ¹⁸O because the records in our analysis are predominantly based on oxygen isotope measurements. These water isotope ratios can be used to track phase changes as the water moves through and among oceans, the atmosphere, and land (e.g., Bowen et al., 2019; Galewsky et al., 2016; Rozanski et al., 1993). Water isotopes are more closely correlated with major modes of variability than precipitation amount, for example, because the majority of diurnal to interannual rainfall δ¹⁸O (δ¹⁸O_precip) variability originates from regional scale hydrological processes like convective intensity, moisture transport history, and moisture source changes (Moerman et al., 2013). Isotope values at a particular location are often linked to teleconnections far afield from the measurement site, highlighting the value of water isotope proxy records for evaluating regional to continental scale phenomena (Puntsag et al., 2016; Vachon et al., 2010; Vuille and Werner, 2005). Therefore, water isotope observations and proxy data may provide useful information about the NAO, even outside the immediate centers of action in the North Atlantic.

Within Europe, previous studies demonstrate that the NAO is responsible for changes in δ¹⁸O_precip by comparing the NAO index with instrumental δ¹⁸O_precip measurements (Baldini et al., 2008; Comas-Bru et al., 2016; Deininger et al., 2016; Langebroek et al., 2011). For example, the continental effect, or the “distance from the coast” effect, is the gradient in amount-weighted δ¹⁸O_precip across the continent, where δ¹⁸O_precip becomes more negative further inland (Dansgaard, 1964). The slope of this continental effect depends on the winter NAO as a function of atmospheric temperature and precipitation (Comas-Bru et al., 2016; Deininger et al., 2016; Hurrell, 1995). Specifically, steeper air temperature gradients across the continent, in combination with decreased atmospheric water content in the continental interior relative to the coasts, drive steeper observed δ¹⁸O west-east gradients across northern Europe during negative winter NAO phases (Deininger et al., 2016; Trigo et al., 2002).

Changes in the moisture source or transport pathway may also influence δ¹⁸O_precip by altering patterns of rainout. In Greenland, for example, past studies apply moisture source and transport models to diagnose the effects that changes in atmospheric circulation and moisture source locations have on δ¹⁸O_precip during different phases of the NAO (Sodemann et al., 2008). Simulated water isotope ratios of Greenland precipitation have higher (lower) δ¹⁸O values during a positive (negative) NAO, as moisture is primarily drawn from the seas north of Iceland (the southwestern North Atlantic). Ice core tops collected from the central region of the ice sheet support these results, although the moisture source and transport models suggest spatial heterogeneity that may not be captured by the relatively limited distribution of ice cores (Sodemann et al., 2008). These initial studies demonstrate mechanistic links between the NAO and water isotopes. Since then, ice-core water isotope records from Greenland have been incorporated into a range of NAO reconstructions and dynamical interpretations (e.g., Sjolte et al., 2018, 2020). However, important uncertainties remain regarding the spatial heterogeneity and stationarity of NAO–δ¹⁸O relationships. Further work is needed to better quantify and reconcile NAO–isotope relationships across sites and proxy archives.

Here we investigate multidecadal variability in the reconstructed NAO over the last millennium (1000–2000 CE) using δ¹⁸O and δD proxies from the PAGES Iso2k database, a global compilation of water isotope-based records spanning the Common Era (Konecky et al., 2020). The global coverage of this database allows us to identify water isotope relationships with the NAO in its centers of action (Greenland and western Europe) and further afield (e.g., the Himalayas and Tibetan Plateau). We evaluate the reconstruction with the instrumental NAO index (Vinther et al., 2003a), observed δ¹⁸O_precip in the Global Network of Isotopes in Precipitation database (GNIP; IAEA/WMO, 2023), and climate model output from the isotope-enabled Community Earth System Model Last Millennium Ensemble (iCESM iLME; Brady et al., 2019). In particular, the isotope-enabled simulations allow us to infer how NAO-induced changes in temperature, precipitation, and moisture transport pathways translate to δ¹⁸O_precip, and therefore how GNIP stations and water isotope proxies in the Iso2k database record the NAO.

2 Methods

2.1 Datasets

We used the following four datasets: the Vinther et al. (2003a) instrumental NAO index (NAO_Vinther), the PAGES Iso2k database (Konecky et al., 2020), the GNIP database containing modern isotope measurements in rainfall (IAEA/WMO, 2023), and the iCESM iLME (Brady et al., 2019).

The Iso2k database is a collection of previously published proxy δ¹⁸O and δD records including ice cores, speleothems, and wood cellulose, that span the Common Era (0–2000 CE). Although the database contains proxies for oceanic conditions, we used only terrestrial Northern Hemisphere records as the NAO is predominantly an atmospheric mode of variability. This approach allowed the most direct comparison with modeled atmospheric variables associated with the NAO. Furthermore, we included only records that had at least annual resolution and 33 % data coverage in the calibration interval between 1823 and 2000 CE. We used the 33 % data coverage threshold based on extensive experiments using different time intervals and data coverage thresholds in attempts to balance the inclusion of records in the analysis while minimizing the influence of missing data. Increasing the data coverage requirement above 33 % dramatically cut the number of included records and limited the spatial coverage of the analysis, while lower thresholds permitted inclusion of records with large sections of missing data. Adjusting the bounds of the time interval also changed the records included in the analysis. Specifically, few Iso2k records have data extending into years more recent than 2000 CE (Konecky et al., 2020), and the NAO_Vinther index is less constrained by observational data before 1823 CE.

The GNIP database contains measurements of both δ¹⁸O and δD in precipitation. We only used δ¹⁸O for our analyses as δD data from stations was less commonly available or discontinuous. Additionally, we only used sites with at least 10 winters of data (not necessarily consecutive), based on the general statistical principle that no fewer than 10 data points should be used to establish a linear relationship (Pearson, 1901). We converted any sub-monthly data to monthly means following Putman and Bowen (2019), and averaged the monthly δ¹⁸O boreal winter (December–January–February – DJF) seasonal means after amount-weighting using the precipitation amount data from each site.

We used three fully forced (i.e., with all transient external forcings) iCESM iLME simulations, each spanning 850–2005 CE, to evaluate the spatial relationships across timescales of δ¹⁸O_precip, temperature, and precipitation with the NAO. The iLME used a ∼2° atmosphere and land and ∼1° ocean and sea ice version of the respective model components. We analyzed monthly mean values for sea level pressure, surface air temperature, total precipitation, wind vectors, specific humidity, and δ¹⁸O_precip. We calculated vertically integrated water vapor transport (IVT) by weighting the monthly mean wind vectors by monthly specific humidity and summing over the full atmospheric column. We defined the winter (DJF) NAO index in the iLME as the difference between the averaged sea level pressure in grid cells containing Reykjavik, Iceland (65° N, 22° W) and Ponta Delgada, Spain (38° N, 25° W). We then regressed this index (NAO_iLME) against DJF sea level pressure, surface air temperature, total precipitation, and δ¹⁸O_precip in the iLME to produce spatial correlation maps of these variables with the NAO. We defined years with positive and negative NAO events as those with winter NAO_iLME index values greater and less than one standard deviation from the mean, respectively.

Finally, we used two published NAO reconstructions for comparison with the Iso2k reconstruction. The first was the calibration-constrained NAO reconstruction from Ortega et al. (2015), based on a network of annually resolved proxy records including ice cores, tree rings, and speleothems. The second was the NAO reconstruction ensemble median from Cook et al. (2019), derived from a tree-ring network. Both reconstructions spanned the last millennium and provided independent estimates of past NAO variability against which the Iso2k reconstruction was evaluated.

2.2 NAO signals in Iso2k

We first annualized the Iso2k records by averaging any sub-annual data into a standard January–December year. We standardized each record to anomaly units by subtracting the mean and dividing by the standard deviation of the full series. This allowed comparisons between oxygen and hydrogen isotope records, as well as between records collected from archive types with differing climate sensitivities.

Age uncertainty was an important consideration when integrating diverse annually resolved proxy records into a large-scale reconstruction. Even within layer-counted archives such as many Greenland ice cores, cumulative counting errors can accrue with depth, and known issues in the Greenland Ice Core Chronology 2005 (GICC05) timescale prior to ∼1200 CE introduce age deviations on the order of several years (Adolphi and Muscheler, 2016; Sinnl et al., 2023). Other ice core chronologies, such as those from Svalbard and some alpine sites, rely on ice flow modeling anchored by age markers (e.g. Isaksson et al., 2005), which increases the age model uncertainty envelope. Speleothem chronologies also include uncertainties associated with uranium-thorium dating precision and age model construction. Because these uncertainties differ in both magnitude and structure across proxy types, their potential to influence reconstructed variability arises not only through absolute age offsets but also through spectral biases.

We restricted the Iso2k dataset to annually resolved, well-dated records to avoid the larger uncertainties associated with radiocarbon-based chronologies. However, we also included a small number of speleothem and flow modeled ice core records, consistent with common practice in other multiproxy reconstructions (e.g., Mann, 2007; PAGES 2k 2013; PAGES 2k 2019; Falster et al., 2023; Table S1 in the Supplement).

We applied a composite-plus-scale (CPS) reconstruction approach to extract the signal of the NAO from the Iso2k records. We produced the CPS ensemble using methods mirroring those applied by previous studies with similar goals (Mann et al., 2008; Neukom et al., 2014). We treated all proxy records uniformly with respect to temporal resolution and seasonal sensitivity. Only three of the 94 Iso2k records possessed subannual resolution, so we used annually averaged values for all proxies to ensure consistent treatment across the network. This procedure followed established precedent in comparable reconstruction studies (e.g., Konecky et al., 2023; Ortega et al., 2015; PAGES 2k Consortium, 2019). Although we incorporated the proxies at an annual timescale, our approach accounted for the seasonal character of the NAO by correlating each proxy series with the instrumental DJF NAO index when weighting the proxy contributions during calibration. This process captured the extent to which each annual proxy series reflected winter NAO variability, resulting in a reconstruction that represented the aggregated isotopic response to the seasonal NAO signal. For each member of the 1000-member ensemble, we randomly removed 15 % of the Northern Hemisphere Iso2k records (Sect. 2.1). We then performed the CPS reconstruction by normalizing and scaling the remaining records according to their correlation with the NAO_Vinther index in the calibration interval from 1873–2000 CE and computing the mean of the records at each timestep. We then normalized mean and variance of the composite to match that of the NAO_Vinther index. The resulting reconstructions should not be interpreted strictly as that of the NAO atmospheric pressure index itself, but rather a reconstruction of the integrated, long-range, and persistent impacts of NAO variability on atmospheric circulation and regional climate as recorded by the isotopic proxy network. We validated this relationship by testing the correlation of the reconstruction ensemble members against the first 50 years of the NAO_Vinther index in the validation interval from 1823–1873 CE.

Note that the choice of calibration and validation windows inherently influenced the number of proxy records available for reconstruction because not all Iso2k records had complete coverage through the 20th century. Using alternative calibration/validation configurations shifted these windows and altered the reconstruction performance by changing proxy availability during the calibration interval. For example, calibrating on an earlier segment of the industrial period (e.g. 1823–1923 CE) and validating on the later period (1924–1999 CE) resulted in the exclusion of substantial portions of shorter records. Because proxy availability varied through time, these tests underscored that the calibration interval must balance record length, data quality, and coverage to ensure a stable and interpretable set of validation statistics. We estimated the impact of changing record coverage throughout the millennium by repeating the CPS validation method after restricting included records according to their temporal coverage. To accomplish this, we performed the reconstruction using only records with start years no later than 1500 CE, and again using only records starting no later than 1000 CE.

We developed a null hypothesis by constructing an ensemble of first-order autoregressive (AR(1)) noise records corresponding to the Iso2k records above and then composited and scaled using the same CPS methods. We then tested the correlation of the resulting composite against the NAO_Vinther index. We repeated this process one hundred times to produce a range of correlations that might be expected from one hundred “reconstructions” of the NAO using autocorrelated noise.

To test the relative influence of records outside of the NAO centers of action on our reconstruction, we performed the same CPS procedures on only records from Greenland and Europe and separately on records from the Indo-Asian monsoon region (Fig. 1, blacked dashed regions). We compared the results with the full reconstruction, and, to ensure the high skill of the full reconstruction was not simply the result of having more records, performed the CPS procedure using records from the regional subsets after adding autocorrelated noise records to bring the total number of records equal to the full reconstruction (“noise-padded”).

Figure 1 Map of Iso2k records included in the analysis. Shapes correspond to archive type and colors denote the Pearson's r-value of each record with NAO_Vinther (1823–2000 CE). Bold outlines signify significant correlation (p<0.1). Regional subsets used for validation are indicated by dashed boxes over Greenland and southeast Asia.

In addition to the spatial correlation analysis, we calculated a model-based “pseudo-reconstruction” of the NAO using the iLME δ¹⁸O_precip field. To mimic the proxy reconstruction framework, we extracted the modeled δ¹⁸O_precip time series from each grid cell corresponding to the terrestrial Iso2k proxy locations. We then weighted each modeled δ¹⁸O_precip series by its correlation with the modeled NAO SLP index over the last millennium (1000–2000 CE). We composited the weighted series to yield a single δ¹⁸O-based pseudo-reconstruction of the NAO. We repeated this process repeated for 1000 ensemble members, each time removing 15 % of the modeled δ¹⁸O_precip grid cells corresponding to Iso2k proxy locations to mimic the data sensitivity experiments in the Iso2k NAO reconstruction. This method provided a process-based test of how the spatial sampling and covariance structure of the proxy network can recover NAO variability in the simulations. It also provided a more direct comparison between the proxy-based reconstruction of the impacts of the NAO with the modeled impacts of the NAO on δ¹⁸O_precip.

We performed spectral decomposition on the NAO_Vinther index, Iso2k NAO reconstruction, NAO_iLME, as well as the NAO reconstructions from Cook et al. (2019) and Ortega et al. (2015) using the multi-taper method (MTM; Thomson, 1982) through the R package geoChronR, which in turn leveraged the astrochron R package (Mckay et al., 2020; Meyers, 2012). We chose to focus on the calibration-constrained reconstruction from Ortega et al. (2015) because the reconstruction methodology most closely matched the methods we used for the Iso2k reconstruction. We performed the spectral analysis of the Iso2k NAO reconstruction and NAO_iLME across the reconstruction ensemble to produce a range of possible MTM results. We compared the resulting decompositions against an AR(1) null hypothesis (e.g., Zhu et al., 2019). We performed wavelet analysis on the NAO reconstructions using the R package “WaveletComp” (Roesch and Schmidbaur, 2018). This method allowed for an investigation of changes in the frequency of signals through time. We also performed wavelet analysis on reconstructions that included only a single proxy type to investigate the impact that record availability and proxy type had on the spectral characteristics of the reconstructions.

We used cross-wavelet power and wavelet coherence to assess the consistency of time domain features between NAO reconstructions. Cross-wavelet power identified regions of shared variance in time-frequency space, while wavelet coherence quantified the strength of the spectral relationship between reconstructions. We evaluated statistical significance against a red-noise (AR(1)) background spectrum, estimating 95 % confidence levels following the implementation in the WaveletComp package. We also compared between reconstructions using an 11-year running mean to investigate decadal-scale variability. We assessed the statistical significance of correlations between the decadally smoothed NAO reconstructions using the Ebisuzaki (1997) phase-randomization approach, which accounted for serial autocorrelation by comparing the observed correlation against a distribution of time series with identical spectral properties (Ebisuzaki, 1997).

3 Results

3.1 NAO-correlated Iso2k records and CPS reconstructions

We extracted a total of 94 records from Iso2k from the Northern Hemisphere with at least annual resolution and 33 % coverage between 1823 and 2000 CE. This included 2 speleothems, 37 wood cellulose, and 55 glacier ice records, with all records based on δ¹⁸O except three glacier ice records of δD. Pearson correlation r-values of these records with the NAO_Vinther index ranged from −0.37 to 0.29. Thirty-seven of the original 94 proxy records were significantly correlated (p-value < 0.1) with NAO_Vinther in the calibration interval including 15 wood cellulose, and 22 glacier ice records. Records with significant correlations spanned the full spatial distribution of the original 94 records, including glacier ice records in Greenland, the Alps, and the Himalayas, and wood cellulose in Europe, northern India, and the Tibetan Plateau (Fig. 1, bold symbols). A total of 20 records spanned the full last millennium including 15 glacier ice, 3 wood cellulose, and 2 speleothem records (Fig. 2b), of which only 7 glacier ice records were significantly correlated with the NAO_Vinther index (Fig. S1a).

Figure 2 (a) NAO reconstruction ensemble members (grey shading) and ensemble median (black line). Horizontal dashed line shows zero. Vertical dashed lines indicate calibration and validation cutoff years at 1823 and 1873 CE. (b) Record coverage through time for the full Northern Hemisphere Iso2k dataset including glacier ice (blue), wood cellulose (green), and speleothem (red) records. Solid black line indicates the total number of available records in each year.

The median of the ensemble of reconstructions ranged from approximately −3.5 to 3.4, while the full ensemble spanned from −4.2 to 4.1 (Fig. 2a). The ensemble members were all significantly (p<0.001) positively correlated with the NAO_Vinther index during the validation interval from 1823–1873 CE, with correlation values ranging from 0.38 to 0.57 and a median of 0.48 (Fig. 3, red).

Figure 3 Box and whisker plots of validation (1823–1873 CE) metrics for the reconstructions using first-order autocorrelated noise (grey), the Northern Hemisphere Iso2k dataset seen in Fig. 1 (red, n=94), only records from Greenland and Europe (blue, n=47), only records from southeast Asia and the Tibetan Plateau (black, n=27), and the combination of Greenland, Europe, southeast Asia, and the Tibetan plateau (green, n=74). Regional subsets correspond with dashed boxes in Fig. 1. Yellow distributions show validation scores when restricting the reconstruction to only records that extend back to 1500 CE or earlier (n=34), and 1000 CE or earlier (n=20). Individual colored squares indicate validation scores for other published NAO reconstructions.

Restricting the spatial extent of records included in the reconstruction to only Greenland and Europe produced a substantial decrease in correlations with the NAO_Vinther index during the validation interval, with correlation values between 0.3 to 0.45, and a median of 0.39. (Fig. 3, blue). Using only records from northern India and the Tibetan plateau in the reconstructions produced a median validation correlation value of 0.44 (with values ranging from 0.27 to 0.51) that exceeded the Greenland and Europe-only reconstructions (Fig. 3, black). Using the two regions together produced the reconstructions with the highest median correlation at 0.5, with individual ensemble values ranging from 0.39 to 0.58 (Fig. 3, green), slightly higher than even reconstructions using all 94 records (Fig. 3, red). Padding the regional reconstructions with autocorrelated noise did not improve their validation correlations, suggesting the skill of the full Northern Hemisphere reconstruction was not simply a function of the large number of records (Fig. S2). The high skill of the two regions combined highlights the importance of including records outside the center of action of the NAO index when reconstructing the NAO.

Reconstruction skill was also dependent on temporal data availability as shown by reconstructions that were restricted to only records extending back to 1500 or 1000 CE (Fig. 3, yellow). Restricting to only records extending to 1500 CE produced validation correlations ranging from 0.06 to 0.47, with a median of 0.32, while records extending to 1000 CE produced reconstruction validations ranging from 0 to 0.4 with a median of 0.28. These validation scores highlight a limitation common to paleoclimate reconstructions, that validation skill during observational intervals with good proxy coverage may not be representative of reconstruction fidelity further back in time when data availability is reduced. Note, however, that because the number of records spanning the full last millennium was lower, the reconstruction validation scores were more sensitive to the random removal of some records. This was reflected in the larger range of validation scores in the distribution, resulting in some subsets of the available data producing validation scores nearing the range of the full reconstruction.

3.2 NAO fingerprints in GNIP isotope stations

Of the 1005 stations extracted from the GNIP database, 146 stations met conditions of at least 10 non-consecutive winter (DJF) seasons of data (Fig. 4 symbols) and of those, 50 had a significant (p<0.1) correlation with the NAO_Vinther index (Fig. 4, symbols with bold outlines). Of those 50 stations, 47 were positively correlated with the NAO, with Pearson's R values ranging from 0.29 to 0.86, while 3 were negatively correlated, ranging from −0.34 to −0.53. The positively correlated stations were concentrated in central Europe, although stations in Finland, the Azores, Barbados, and Greece were also positively correlated.

Figure 4 Winter NAO_iLME correlations (color contours) with iLME winter δ¹⁸O_precip from 1823–2000 CE. Stippling denotes p <  0.01. Symbols show GNIP stations with ≥10 winters colored by δ¹⁸O_precip correlations with NAO_Vinther, bold symbols show significant (p <  0.1) correlations.

The negatively correlated stations were in Iceland, southwest Russia, and Turkey. The correlation dipole between the Iceland station and the stations in Europe was broadly consistent with the sea level pressure (SLP) dipole characteristic of the NAO, in which the NAO was negatively correlated with SLP over Iceland and positively correlated over Europe (Fig. S3). This pattern was also consistent with past studies of δ¹⁸O correlations with the NAO (Deininger et al., 2016).

3.3 NAO fingerprints in the iLME

The winter (DJF) NAO_iLME index values ranged from −5.29 to 3.15 with a standard deviation of 1.79 through the period 1823–2000 CE (Fig. S4a). As with the NAO_Vinther index, the winter NAO_iLME was slightly positively skewed, but with more extreme values during negative years (Fig. S4b–c). The NAO_iLME captured the main spatial features of observed NAO impacts in the North Atlantic region, with similar correlations to SLP, temperature, and precipitation (Compo et al., 2011; Hernández et al., 2020). Correlations between NAO_iLME and winter precipitation amount were less spatially extensive than with temperature (Fig. 5a–b). However, the precipitation correlations with the NAO_iLME showed an east/west contrast over Greenland, with positive correlations in the east and negative correlations spanning western Greenland and the northeastern Canadian Arctic. This differed notably from 20th Century Reanalysis correlations which showed no significant NAO-precipitation correlations over Greenland (Hernández et al., 2020). Positive winter NAO_iLME years were characterized by greater IVT from the North Atlantic and stronger anticyclonic flow in the subtropical Atlantic, which increased moisture transport from the North Atlantic into Northern Europe and Scandinavia relative to negative winter NAO_iLME years, consistent with a stronger and more northern jet (Fig. 5c).

Figure 5 Winter NAO_iLME correlations (color contours) with iLME winter surface air temperature, precipitation amount, and composite IVT from 1823–2000 CE. Stippling denotes p <  0.01. (a) NAO_iLME correlation with winter surface air temperature. (b) NAO_iLME correlation with winter precipitation amount. (c) Difference between winter IVT during positive and negative winter NAO_iLME years, defined as NAO_iLME greater or less than one standard deviation from the mean, respectively, over the period 1823–2000 CE.

Winter NAO_iLME correlations with amount-weighted δ¹⁸O of winter precipitation (δ¹⁸O_precip) were negative over western Greenland and Iceland, and positive over Europe, Scandinavia, and the mid-latitude North Atlantic (Fig. 4, shading). The NAO_iLME index was also positively correlated with winter δ¹⁸O_precip over the tropical Pacific, west Africa, and the Tibetan plateau. Although GNIP δ¹⁸O_precip observations are limited, and nonexistent over the oceans, the correlations were consistent with the available GNIP δ¹⁸O_precip data (Fig. 4, symbols).

3.4 Spectral analysis of the NAO reconstruction and the NAO in the iLME

Wavelet decomposition of the NAO reconstruction found significant (p<0.05) power in four main “bands” throughout the last millennium, which here we refer to as interannual (1–10 years), decadal (10–30 years), multidecadal (30–100 years), and centennial (100+ years; Fig. 6a). The multi-taper method spectral decomposition of the full reconstruction over the last millennium showed higher power at longer periods (Fig. 6b). Peaks that rose above the 95 % confidence level occurred in the interannual band (2–6 years), as well as at decadal, multidecadal, and centennial scales (11, 23, 94, and 300+ years). The wavelet analysis of the NAO reconstruction showed significant power in the interannual and decadal bands from the early 1700s through present (Fig. 6a). Power in the interannual band was weaker prior to 1550 CE. Significant decadal and multidecadal power also persisted from ca. 1000–1300 CE, despite being mostly absent by the mid-millennium. Significant multidecadal to centennial power extended throughout the full reconstruction, although with lower power than the higher frequency bands.

Figure 6 (a) Wavelet decomposition of the NAO reconstruction from 1000–2000 CE. Color contours are scaled by power quantiles. Black contours bound statistically significant power (p<0.05) and white curves show ridges (power maxima). White envelope shows the region where padding may influence estimates. (b) The multi-taper method spectral decomposition for the NAO reconstruction. Red line shows 95 % confidence level based on an AR(1) null. Grey shading shows ensemble range. (c) and (d) are as in (a) and (b) but for the NAO pseudo-reconstruction from the iLME δ¹⁸O_precip field.

We repeated the wavelet and multi-taper analysis twice to investigate the impacts of proxy type on the spectral character of the reconstruction: First using a reconstruction based only on glacier ice records, and again with a reconstruction using only wood cellulose records (Fig. S5). Like the full reconstruction, both proxy-specific reconstructions exhibited significant interannual to decadal power, but glacier ice records had higher multidecadal to centennial power than wood cellulose records.

The reconstruction using glacier ice records had consistent significant power across scales from the late 1700's through present and from 1000 to about 1250 CE (Fig. S5a). Nevertheless, the mid-millennium had only intermittent significant power on interannual to decadal scales. By contrast, power on multidecadal and centennial scales was notably weaker in the later part of the millennium in the reconstructions using wood cellulose records. While these reconstructions had significant power in the interannual to decadal band from about 1600–2000 CE (Fig. S5b), limited availability of wood cellulose records prior to 1600 CE means the wavelet decomposition was based only on a small number of records and likely not reliable (Fig. 2b).

Ensemble-based MTM analyses showed that variability in record selection (15 % removed per realization) could cause some ensemble members to fall below the AR(1) 95 % confidence threshold in the multidecadal band, suggesting that the magnitude of multidecadal variability is partially sensitive to network size and composition (Fig. 6b). However, the wide diversity of the network reduced the influence of any single record, and the shared climate signal could still be robustly identified despite heterogeneity in age uncertainty (Ortega et al., 2015). Sensitivity analyses confirmed this, as removing the non-layer-counted records either individually or collectively produced only minor changes in the reconstruction, with the most notable being a slight reduction in the multidecadal spectral peak when all were excluded (Fig. S6b). The low frequency variability in the wavelet remained consistent with the full reconstruction (Fig. S6a). Given the small age uncertainty introduced by the updated Greenland chronology ($\sim \pm \mathrm{5}$ years between 1000–1200 CE), and the limited leverage of the relatively few non-layer-counted records, we conclude that age-model uncertainty is unlikely to be a dominant contributor to the reconstructed low-frequency variability. Nonetheless, we highlight it as an important dimension of proxy heterogeneity that influences data synthesis efforts.

The wavelet decomposition of iLME pseudo-reconstruction of the NAO exhibited characteristics that were largely distinct from any of the reconstructions. Specifically, interannual variability dominated the pseudo-reconstruction over the last millennium, with the only significant (and still relatively weak) decadal power early in the millennium (Fig. 6c). The multi-taper method decomposition of the pseudo-reconstruction confirmed the lack of multidecadal and centennial scale power (Fig. 6d) and the generally “white noise” character of pseudo-reconstruction with several significant inter-annual peaks (2–8 years), one decadal peak (∼11 years), but overall low power at lower frequencies.

3.5 Spectral comparisons with other NAO reconstructions

We assessed low-frequency agreement between the Iso2k reconstruction and the NAO_Vinther index through decadal smoothing, cross-wavelet power, and wavelet coherence on the portions of the reconstructions that covered the instrumental interval. The 11-year centered running means showed broadly coherent decadal variability between the Iso2k reconstruction and the NAO_Vinther index (r=0.59, Ebisuzaki p=0.004, Fig. 7a). In the NAO_Vinther index, low frequency variability was present throughout, including a prolonged positive phase during the early- to mid-twentieth century as well as a shift toward positive values in the final ∼2 decades. The Iso2k reconstruction generally captured these decadal variations. Decadally smoothed reconstructions from Ortega et al. (2015) and Cook et al (2019) also showed moderate positive correlations with NAO_Vinther (Fig. 7b, r=0.45 and Fig. 7c, r=0.46, respectively), although the Cook et al. (2019) correlation was significant (Ebisuzaki p=0.024), while the Ortega et al. (2015) correlation was not (Ebisuzaki p=0.11). Cross-wavelet power analysis between the Iso2k reconstruction and the NAO_Vinther index identified statistically significant shared variance at interannual to decadal periodicities (∼5–10 years), with some significance at longer periods (∼20 years) during the late nineteenth and early twentieth centuries (Fig. S7a). Wavelet coherence also showed high coherence at (∼10–20-year periodicities across much of the record (Fig. S7b).

Figure 7 Decadally smoothed NAO reconstructions during the observational period compared with the decadally smoothed NAO_Vinther index (black lines). Series are smoothed using an 11-year centered running mean. (a) Iso2k CPS reconstruction (this study; orange) vs. NAO_Vinther. (b) Ortega et al. (2015) calibration constrained median (teal) vs. NAO_Vinther. (c) Cook et al. (2019) NAO_med (blue) vs. NAO_Vinther.

The spectral characteristics of the Ortega et al. (2015) and Cook et al. (2019) reconstructions generally resembled those of the Iso2k NAO reconstruction (Fig. 8). The wavelet decomposition of the Ortega et al. (2015) reconstruction showed significant power in the multidecadal to centennial band (30–200 years), particularly during ca. 1100–1600 CE, with weaker interannual to decadal power (Fig. 8a). The MTM power rose above the 95% confidence level at multidecadal to centennial periods, as well as on interannual timescales (2–5 years; Fig. 8b). The Cook et al. (2019) reconstruction showed significant multidecadal variability from ca. 1000–1400 CE (Fig. 8c), which was consistent with the variability identified in the Ortega et al. (2015) reconstruction and the Iso2k reconstruction during this period. The Cook et al. (2019) MTM showed elevated power at low frequencies but did not rise above the 95 % confidence level (Fig. 8d), suggesting that low-frequency variability was present but not as extensive throughout the last millennium as in the Iso2k or Ortega et al. (2015) reconstructions.

Figure 8 As in Fig. 6 but for the NAO reconstruction ensemble means from Ortega et al. (2015) (a, b) and Cook et al. (2019) (c, d).

Over the last millennium the decadally smoothed Iso2k NAO reconstruction was positively correlated with decadally smoothed medians from the Ortega et al. (2015) reconstruction (r=0.33, Ebisuzaki p<0.001, Fig. S8a) and the Cook et al. (2019) reconstruction (r=0.29, Ebisuzaki p<0.01, Fig. S8d). Cross-wavelet comparisons showed non-stationary variability in the multidecadal to centennial band (Fig. S8b, e), and wavelet coherence showed enhanced variance at periods of ∼30–100 years, particularly during ca. 1200–1500 CE and post-1800 CE portions of the reconstructions (Fig. S8c, f). These intervals overlapped with periods of elevated multidecadal power in the single-series wavelet decompositions (Figs. 6, 8).

4 Discussion

4.1  Regionally variable drivers of NAO-δ¹⁸O_precip relationships in Europe

There is a clear imprint of the NAO on records in the Iso2k database. Opposing positive and negative correlations in northern Europe and southern Greenland, respectively, are consistent with the observed contrast in temperature and precipitation associated with the winter NAO (Fig. 1; Hernández et al., 2020). Correlations between the NAO and observed and modeled δ¹⁸O_precip suggest that the temperature effect is a major control on these relationships, especially at high latitudes in western Greenland, northern Europe, and Scandinavia (Figs. 4, 5a). Specifically, warm air temperatures over northern Europe during a positive NAO cause weaker temperature-dependent Rayleigh fractionation during condensation and decreased rainout along the westerly path of distillation, while cool air temperatures during negative NAO cause stronger Rayleigh fractionation and increased rainout (Dansgaard, 1964), all of which should produce a positive δ¹⁸O_precip/NAO relationship. However, temperature is clearly not the only mechanism by which the NAO influences δ¹⁸O_precip. Globally, in locations where a precipitation “amount effect” is observed (particularly, though not exclusively, at low latitudes), the amount of precipitation is inversely correlated with δ¹⁸O_precip (Nusbaumer et al., 2017). The NAO impacts precipitation amounts in both Europe and Greenland by redirecting storm tracks in the Atlantic (Hurrell, 1995; Hurrell and Deser, 2009; Woollings et al., 2010, 2018; Woollings and Blackburn, 2012). Unlike temperature, precipitation amount-driven impacts of the NAO would therefore yield a negative correlation with δ¹⁸O_precip, with more negative δ¹⁸O_precip in areas experiencing increased precipitation and vice versa; in this way, the precipitation amount effect opposes the temperature effect. In observations and in the iLME, negative NAO-precipitation amount correlations corresponded with positive NAO-δ¹⁸O_precip correlations over the Azores, the subtropical Atlantic, and the tropical Pacific, consistent with an amount effect (Figs. 4, 5b). Changes in local precipitation amount as well as upstream fractionation during anomalous transport (Fig. 5c) can therefore explain the far reach of the NAO in these lower-latitude regions.

Interestingly, the NAO-δ¹⁸O_precip correlation over Europe is uniformly positive, despite the canonical north-south contrast in the NAO's influence on temperature and precipitation (Fig. 5a–b). This likely reflects a latitudinal difference in the primary controls on δ¹⁸O_precip, and thus the NAO-δ¹⁸O_precip relationship, across the continent. A positive, direct relationship with temperature can explain NAO-driven δ¹⁸O_precip in northern Europe and Scandinavia: Positive NAO phases bring warm conditions to this region, which translates to higher δ¹⁸O_precip (Figs. 4, 5a). Conversely, the negative relationship with precipitation amount can explain NAO-driven δ¹⁸O_precip in southern Europe. There, positive NAO phases reduce precipitation, resulting in higher δ¹⁸O_precip values (Figs. 4, 5b). In both regions, therefore, positive (negative) NAO conditions produce higher (lower) δ¹⁸O_precip values, resulting in a positive correlation between the NAO and δ¹⁸O_precip across the continent. Previous studies have relied on north-south gradients in European precipitation to reconstruct the NAO (e.g., Trouet et al., 2009). While this is valid for precipitation-driven proxies, our results indicate that studies utilizing δ¹⁸O_precip-based proxies must be careful to not presume a north-south gradient in δ¹⁸O_precip, because the drivers of NAO-δ¹⁸O_precip relationships are spatially variable across Europe.

4.2 NAO teleconnections with South and Central Asia

The positive correlations between the NAO_Vinther index and wood cellulose records from the Himalayas and Tibetan Plateau are strong enough that our NAO reconstructions perform best when they are included. We suggest that the winter NAO-δ¹⁸O_wood relationship arises from the combined influence of temperature, precipitation, and changes in moisture source and transport as they are integrated across seasons in the wood cellulose. During winter, the westerlies intensify and flow south of the Tibetan Plateau, bringing moisture to the region from the eastern Mediterranean and western Central Asia (Liu et al., 2017; Schiemann et al., 2009). However, because the eastern Mediterranean and western Central Asia are less humid during positive NAO years, upstream winter IVT is weaker, bringing less winter moisture to the Himalayas and Tibetan Plateau from these distal, westerly moisture sources (Fig. S9). Compared with local sources, moisture imported from distal sources should have lower δ¹⁸O from rainout along the transport path. We speculate that the reduction in imported moisture during positive NAO years would result in higher δ¹⁸O of winter precipitation, contributing to positive correlations between the winter NAO and the δ¹⁸O in wood cellulose records over long timescales.

The Tibetan plateau generally lacks observational δ¹⁸O_precip data that is long and continuous enough to rigorously evaluate the regional NAO-δ¹⁸O_precip relationships, but one station from Lhasa with discontinuous coverage between 1994–2006 CE shows the NAO is significantly negatively correlated with winter temperature and significantly positively correlated with winter precipitation amount (Yao et al., 2013; Table S2). At Lhasa, colder winter temperatures are associated with more negative δ¹⁸O_precip, whereas increased winter precipitation is associated with more positive δ¹⁸O_precip (opposite to the traditional “amount effect”; Dansgaard, 1964) (Table S2). Thus, during the colder and wetter winters of positive NAO years, the competing influences on δ¹⁸O_precip result in a weakly positive, though not significant, correlation between the winter NAO and δ¹⁸O_precip. While additional, continuous station data would be required to better evaluate NAO-δ¹⁸O_precip relationships in the Tibetan Plateau, multiple studies from central Asia have observed NAO correlations with δ¹⁸O_wood (Liu et al., 2009, 2015; Xu et al., 2021, 2019), supporting the connection which we believe is best explained through moisture source and transport changes.

Like observations, the iLME also shows significant negative correlations between the NAO and precipitation to the northwest (upstream) of the Tibetan Plateau, suggesting drier upstream conditions during positive NAO years, as well as significant positive correlations with δ¹⁸O_precip (Figs. 4, 5c). However, the relationship between the winter NAO_iLME and surface air temperature over the Tibetan Plateau is opposite that of observations (Fig. 5a). The Himalayas and the Tibetan Plateau are notoriously challenging to resolve in lower-resolution climate models such as that used for the iLME, and contribute to large local biases in precipitation amount, seasonality, moisture source (Li et al., 2022; Lin et al., 2018), and winter δ¹⁸O_precip (Gao et al., 2011). We must therefore use caution in interpreting the iLME results local to the Tibetan Plateau and ultimately suggest that reduced influx of depleted moisture from distal sources leads to overall higher δ¹⁸O values across the Tibetan Plateau during positive NAO years, contributing to the relationship observed in the wood cellulose δ¹⁸O records. Longer, continuous station data (>10 years) paired with higher-resolution simulations and proxy system models are needed to further test these mechanisms.

We note that these mechanisms are distinct from previous interpretations that the NAO-δ¹⁸O_wood relationship arises from teleconnections between the summer NAO and the Indian Summer Monsoon (ISM) (Sano et al., 2013; Wernicke et al., 2017; Xu et al., 2018). The positive NAO-ISM teleconnection produces a negative NAO-δ¹⁸O_wood relationship largely via the amount effect (Bamzai and Shukla, 1999; Rajeevan, 2002; Raman and Maliekal, 1985; Robock et al., 2003), which cannot explain the positive correlation that we observe. The positive winter NAO_Vinther-δ¹⁸O_wood correlation that we observe, and the Lhasa station data suggest that central Asian δ¹⁸O_precip increases following positive winter NAO (Figs. 1, 4). We therefore conclude that while the NAO-ISM teleconnection is important for summer δ¹⁸O_precip, the winter NAO's influence on winter δ¹⁸O_precip strongly influences δ¹⁸O_wood in locations where snowmelt plays a substantial role during the growing season.

4.3 Low frequency variability in the NAO reconstructions

In our reconstruction, we see strong evidence for multidecadal variability of the NAO (Fig. 6a), which has previously been debated, but also that there is nonstationarity in its strength over the last millennium. Significant decadal and multidecadal variability persists during the Medieval Climate Anomaly (MCA) from ca. 1000–1250 CE and then weakens during the early Little Ice Age (LIA), especially between ca. 1300–1700 CE, before increasing again in the final 300 years of the millennium. We argue that water isotope-based proxy records used in our reconstruction integrate climate processes on broad spatiotemporal scales, and therefore are well-poised to capture decadal or longer scale signals that can be missed by proxy systems like tree ring width and density, which are more sensitive to local high frequency variability (Ault et al., 2013; Dee et al., 2017; Moerman et al., 2013).

The nonstationary multidecadal variability in our reconstruction is not simply an artifact of changing proxy types and locations through the last millennium. Changes in proxy types could in theory be particularly important after 1600 CE, when many wood cellulose records become available; wood cellulose records have been suggested as biased towards high-frequencies (Dee et al., 2017). However, multidecadal variability is also evident in the wavelet decompositions performed on proxy-specific NAO reconstructions (Fig. S5). Furthermore, the wood cellulose records that do have coverage between ca. 1150–1300 CE (Fig. S5b) also share the increased multidecadal variability reflected in the ice core-based (Fig. S5a) and full NAO reconstruction (Fig. 6a). Moreover, NAO variability is strong across scales during the MCA when changes in record availability are gradual, and, more generally, the timing of nonstationarity does not coincide with changes in record availability (Figs. 2b, 6a). Nevertheless, accumulating chronological uncertainty in layer-counted ice core records may contribute to the spectral characteristics of the reconstruction. As dating errors grow back in time, high-frequency variability may be slightly dampened and redistributed toward lower frequencies, potentially contributing to the reduced interannual power prior to ∼1700 CE (Fig. S5a).

Nonstationarity in the variability of the NAO across scales through the last millennium is also evident in other NAO reconstructions, but the robustness of low frequency variability has been debated. Landmark reconstructions including Ortega et al. (2015) and Trouet et al. (2009) showed extensive low frequency variability throughout the last millennium (e.g. Fig. 8a, b), while Cook et al. (2002, 2019) argued that NAO variability during the last millennium is more consistent with “white noise”. That said, closer inspection of the power spectrum presented in Cook et al. (2019) also shows similar nonstationarity to our NAO reconstruction, with high spectral power in multidecadal bands during the MCA (ca. 950–1300 CE), which then decreases during the mid-millennium and LIA (Fig. 8c).

Significant decadal variability is present in both the instrumental NAO_Vinther index and the NAO reconstructions, especially for Iso2k and Cook et al. (2019) (Fig. 7a, c). The weaker, non-significant correlation with Ortega et al. (2015) likely results from the shorter overlap and pronounced divergence from the instrumental record after ca. 1960 CE (Fig. 7b). Over the last millennium the decadally smoothed series and individual wavelet analyses show moderate correspondence at decadal to centennial timescales, particularly during intervals of enhanced variance such as the MCA, indicating that the Iso2k NAO as well as the Ortega et al. (2015) and Cook et al. (2019) reconstructions capture similar low-frequency behavior (Figs. 6, 8). Notably, the cross-wavelet power and wavelet coherence analyses share strong signals in the multidecadal to centennial bands during the early millennium and the MCA, highlighting the shared signal across reconstructions at that time (Fig. S8b–f). However, this relationship is nonstationary, with significant coherence occurring only intermittently in time and across limited frequency bands. This highlights variability in the phase and amplitude relationships between reconstructions and through time, which is unsurprising given the differences in proxy networks, reconstruction methods, and chronological uncertainties. In this context, the wavelet coherence results should be interpreted as evidence for a temporally evolving relationship between reconstructions.

All of the reconstructions included here agree that the low-frequency variability in the NAO may be inherently nonstationary. We hypothesize that it arises from the NAO's interactions with changing latitudinal gradients in Northern Hemisphere temperature, and specifically in the state of the Atlantic Ocean. Latitudinal gradients in sea surface temperature between the tropical and high latitude North Atlantic determine the strength and position of the Azores High and Icelandic Low. Warmer Northern Hemisphere temperatures during the MCA, for example, corresponded with a weakened latitudinal temperature gradient (Li et al., 2013), which may have enhanced decadal to multidecadal variability in the NAO (e.g., Wang et al., 2017). More broadly, it is difficult to reconcile climate shifts across Europe and North America between the MCA and LIA (Edwards et al., 2017; Goosse et al., 2012), and explanations often invoke minor changes in solar and volcanic forcing that were then amplified by internal variability (Goosse et al., 2005; Graham et al., 2011). Better understanding the role for these forcings in multidecadal variability of the NAO, and its nonstationarity from the MCA to LIA, should be a priority (Birkel et al., 2018; Fernandez et al., 2025; Otterå et al., 2010).

It should be noted that the mechanisms underlying climate variability during the majority of the last millennium are quite different from those responsible for modern anthropogenic climate change (Goosse et al., 2012), as the magnitude of anthropogenic greenhouse gas forcing far outweighs the magnitude of natural forcings (Kuzmina et al., 2005; Stephenson et al., 2006). The positive trend and increase in significant spectral power across interannual, decadal, multidecadal, and centennial timescales in the industrial era of our reconstruction are suggestive of amplification of variability in Northern Hemisphere atmospheric circulation as a response to anthropogenic climate change (Gillett et al., 2002, 2003; Moore et al., 2017).

The multidecadal variability in the proxy reconstruction is not reproduced in the pseudo-reconstruction of the NAO derived from the iCESM iLME simulations. This discrepancy aligns with previously documented mismatches between the spectral characteristics of model simulations and those of proxies (Ault et al., 2013; Dee et al., 2017; Franke et al., 2013). Recent work does find multidecadal variability in a spectral analysis of the North Atlantic region in the Last Millennium Ensemble, though high power at multidecadal timescales was limited to periods with major volcanic eruptions and absent from purely internal variability (Fernandez et al., 2025). The fundamental cause of this discrepancy between reconstructions and models – whether due to inherent biases in the models, limitations of the proxies themselves, or a combination of both – remains unresolved.

5 Conclusions

Water isotope-based proxies capture the broad range of spatial and temporal impacts of the NAO because they integrate temperature responses across Greenland and Northern Europe, local precipitation amount responses over the Iberian Peninsula and Southern Europe, and moisture source and transport responses in teleconnected regions like South Asia. These relationships explain the skill of our water isotope-based NAO reconstruction, especially when leveraging proxy data from teleconnected regions outside the North Atlantic such as Asia. Our NAO reconstruction reveals strong but nonstationary multidecadal variability across much of the last millennium, which is also present in other independent reconstructions of the NAO, and future work should evaluate the mechanisms by which external and internal mechanisms affect the timescales of variability in the NAO. Importantly, iCESM (like other models) underestimates the NAO variability on multidecadal timescales in the reconstruction, suggesting that multidecadal variability could also be underestimated in projections of future climate change. Accounting for low-frequency variability in the NAO as illustrated here will be critical for assessing the risk of NAO impacts to temperature and precipitation across the North Atlantic region in the future.