Evaluating the decadal-to-centennial evolution of a new proxy-based NAO reconstruction during the Common Era Hernandez

3 The North Atlantic Oscillation (NAO) is the major atmospheric mode ruling European climate variability during winter and its significance is underpinned by the number of recent studies aimed at reconstructing past NAO variability across different time scales and temporal resolutions. We present a new 2000-year multi-annual, proxy-based, local NAO impact reconstruction, with associated uncertainties, employing a Bayesian approach. This new local NAO reconstruction is focused on the Iberian Peninsula, a geographical area not included in previous NAO reconstructions despite being a widely used region for instrumental-based NAO measurements. We also assess the main forcing drivers in order to unveil any discrepancies with previously published NAO reconstructions. Results show that, on a decadal scale, a low number of sunspots correlate to low NAO values. The comparison with other previously published NAO reconstructions allows us to observe the variability of solar influence on the NAO across a latitudinal gradient according to the position of the employed archives. Moreover, we highlight the potential role of other North Atlantic modes of variability (i.e., East Atlantic pattern) on the non-stationary behaviour of the NAO through the Common Era, likely via solar impact.

The NAO reconstructions available so far generally agree with each other as well as with instrumental NAO indices until as late as c. 1850 CE (Common Era). However, considerable discrepancies are evident further back in time (Pinto and Raible, 2012). Multiple factors have been put forward to account for these disagreements, ranging from the quality of instrumental data, chronological uncertainties, the use of diverse calibration periods, the rise of temperature linked to the anthropogenic climate change, differences in the sensitivity of the archives to climate variations and/or limitations on the stationary relationship assumed between the climate sensors and the NAO (Lehner et al., 2012;Michel et al., 2018;Schmutz et al., 2000;Sjolte et al., 2018;Timm et al., 2004;Wanner et al., 2001;Zorita and González-Rouco, 2002). The relatively large internal variability of the NAO has been commonly associated with its non-stationary behaviour (Lehner et al., 2012;Raible et al., 2006). More recently the amplitude of the internal variability has been ascribed to the influence of other North Atlantic modes of climate variability (i.e., the East Atlantic-EAand Scandinavian-SCA-patterns) that would modulate the variations in the strength and location of the NAO dipole from annual to multidecadal scales (Comas-Bru and McDermott, 2014;Mellado-Cano et al., 2019;Moore and Renfrew, 2012). Another important contributor to reduce the discrepancies amongst NAO reconstructions is to use more accurate definitions (i.e., local impacts vs regional reconstructions). For instance, a regional NAO assembling a large set of sub-decadal data from different archives and localities could provide much better constrains (Cook et al., 2019;Ortega et al., 2015;Sjolte et al., 2018) (Fig. 1a). However, each of those data sources may be recording the sign in a different way and this would result in a smoothed NAO signal. Many studies use just one or two records (Baker et al., 2015;Faust et al., 2016;Olsen et al., 2012;Trouet et al., 2009) (Fig. 1a), which provide a more accurate (avoiding over-smoothing) yet spatially limited reconstruction of the (local) NAO impact.  20CRv2c;Compo et al., 2011). Location of the proxy-based records (ice cores, lake sediments, speleothems, tree rings, and marine sediments) employed in this study, using symbols and colours to represent the different types of archives and the reconstructed climate variables. (b) Correlation distribution maps between the winter precipitation and temperature (wPre and wTmp) datasets and the NAO, for the boreal winters (December -February) between 1901 and 2016, calculated using the CRU-TS4.1 global climate dataset (Harris et al., 2014) and the NAO and EA indices from Comas-Bru and Hernández (2018). Positive Spearman rank correlations are shown in red and negative correlations are shown in blue. Location of record used for NAOIP is indicated.

6
Although the impact of external forcing mechanisms on the NAO is still a matter of debate (Swingedouw et al., 2017), it has been traditionally assigned to the volcanic activity (Wanner et al., 2015), as seen by a predominance of positive NAO phases during periods of increased volcanic activity Sjolte et al., 2018). The role of solar activity, however, is even more controversial, with contradicting conclusions from proxy-based reconstructions (Brugnara et al., 2013;Faust et al., 2016;Ortega et al., 2015;Sjolte et al., 2018). Modelling (Chiodo et al., 2019;Woollings et al., 2010) and observational (Chiodo et al., 2012;Ineson et al., 2011) studies also display contradictory conclusions between the 11-year solar cycle and the NAO relationship (Chiodo et al., 2019;Thiéblemont et al., 2015). The studies arguing for a solar impact on the NAO do it through a "top-down" mechanism related to the ultraviolet irradiance pattern ). An increase of UV radiation during periods of high solar activity warm up the middle atmosphere leading to an altered stratospheric circulation that propagates pole-and down-wards to affect tropospheric jet streams and thus atmospheric circulation (Ineson et al., 2011;Martin-Puertas et al., 2012). However, the response of the NAO to the solar cycle would not occur immediately but rather after a lag of c. 3 years. This is because the impact of solar heating accumulates for several years in the ocean, and because of a positive feedback between the ocean and atmosphere Scaife et al., 2013). Low sunspot activity results in a climate pattern very similar to the negative phase of the NAO (Lockwood, 2012) with longer lasting and more intense blocking episodes than during high solar blocking events (Barriopedro et al., 2008).
A recent proxy-based study (Faust et al., 2016) further supports a linkage between the Grand minima of solar activity and negative NAO phases that accompany cooling events (e.g. Little Ice Age-LIA-) at decadal-to-centennial timescales.
NAO variability only accounts for c. 40 % of the climate variance that is ultimately captured by palaeorecords (Pinto and Raible, 2012), and on this basis a perfect NAO reconstruction is far from achievable beyond the instrumental period. For this, the interpretation of available NAO reconstructions requires a careful re-examination. The sensitivity of different climate archives to the NAO may vary at different spatio-temporal scales. This highlights the need to choose localities where the NAO impact on local climate is likely to be (more) stationary and to select paleo-sensors-the seasonality of which coincides to the season where the NAO's impact on climate is larger (Lehner et al., 2012;Ortega et al., 2015).
Lakes are highly sensitive to climate through changes in air temperature, wind speed, precipitation, and solar irradiance (Margalef, 1983). These climatic controls on lake and catchment processes impinge on the water properties of the lake and result in its sediments recording climate variability over timescales from seasonal to decadal (Blenckner et al., 2007;Catalan et al., 2013;Pla-Rabes and Catalan, 2011;Williamson et al., 2009). Because the NAO controls the variability of a wide range of climate-related variables (e.g., temperature, precipitation, wind-speed, solar radiation), sedimentary records capturing lake dynamics may integrate NAO variability better than a single local weather record .
Most of these NAO reconstructions have been based on the use of some variant of regression models, often coupled with Principal Component Analysis (Michel et al., 2018). By contrast, richer models using Bayesian inference have been extensively used during the last decade for age-depth chronological building (Blaauw et al., 2007;Haslett and Parnell, 2008;Ramsey, 2008) as well as for climate and environmental reconstructions using biological proxies (Cahill et al., 2016;Haslett et al., 2006;Parnell et al., 2015;Tingley and Huybers, 2009). Nevertheless, neither Bayesian inference nor non-biological proxies have yet been used reconstruct the NAO. The Bayesian approach holds a major advantage over traditional methods, as it is conceptually simpler to build a complex model which quantifies the relationship between multiple proxy and climate variables simultaneously-rather than relying on individual coefficients to describe the relationship (Birks et al., 2010;Li et al., 2010). Furthermore, it is possible to model observations under all conditions (i.e., modern analogues) (Guiot and de Vernal, 2007). In contrast, the price to be paid for ''no modern analogue'' situations means considerably larger uncertainties, which can be however, accounted for in the resulting reconstructions.
Here we present a quantitative reconstruction of the NAO for the central Iberian Peninsula (IP) over the last two millennia, along with its uncertainties, applying a Bayesian approach.
We also assess the coherence between our new NAO reconstruction and previously published reconstructions from other locations, as well as potential external forcing mechanisms that would lead to disagreements between them as a result of the non-stationary spatial behaviour of the NAO.

2.1-Proxy data
We used the chemical composition of a lacustrine sediment core (CIM12-04A, 124.8 cm long) retrieved from an alpine lake located in the Central Iberian Range (Cimera Lake, 40°15'N -5°18'W, 2140 m a.s.l.) (Sánchez-López et al., 2016). The chemical composition of the sediments was obtained by continuous X-ray fluorescence (XRF) analysis using the XRF Avaatech core scanner located at the University of Barcelona (Spain). The XRF settings (working conditions) were: 2 mm of spatial resolution, 2 mA, 15 s count times and 10 kV for lighter elements, with 55 s and 30 kV for heavier elements. Thirty chemical elements were measured, but only ten light (Al, Si, K, Ca, Ti, V, Cr, Mn, Fe and Zn) and three heavy (Rb, Sr and Zr) elements had enough counts to be considered robust.
The chronology of the sediment deposition of the CIM12-04A core have been previously determined by Sánchez-López et al., (2016). It was derived using the activitydepth profile of 210Pb in the uppermost 9 cm of the core together with six AMS 14C dating.

2.2-Climate datasets
We used the NAO extended winter index (Jan-May) spanning the period 1824 -2012 CE was employed to produce the NAO influence reconstruction. The data were obtained from the Climatic Research Unit (CRU) at the University of East Anglia (UK) (https://crudata.uea.ac.uk/cru/data/pci.htm). This NAO index was defined by Jones et al., (1997) and modified by Vinther et al., (2003). It was defined as the difference between the normalized monthly sea level pressure anomalies recorded at Reykjavik (Iceland) and those observed at Iberia (Gibraltar/Cadiz). Precipitation and temperature datasets employed for figures were obtained from CRU-TS4.01 (Harris et al., 2014), whereas NAO and EA data for figures were acquired from Comas-Bru and Hernández (2018).

2.3-Bayesian Model
We follow the Bayesian modelling approach of Parnell et al., (2016Parnell et al., ( , 2015 to produce a reconstruction of NAO impact in the central IP. The relationship between proxy and climate is derived from a training data set for the instrumental/proxy calibration period and is expressed through a likelihood function. This function is combined with a prior probability density function containing parameter information to obtain a posterior probability distribution of the reconstructed NAO values using Bayes' theorem (Juggins and Birks, 2012). Whilst Parnell et al., (2015) based their framework on reconstructing multivariate temperature and moisture measurements from raw pollen data, the method is easily adaptable to other proxies and climate variables. Indeed, Cahill et al., (2016) used a similar approach to reconstruct sea level from foraminifera. In all cases the measures/counts of the proxy are required for a set of sediment layers (depths) in a core.
We summarize the mathematical details of the model in this section. Full technical details for the model fitting process are described in Parnell et al. (2015). We provide all the code used to create the reconstructions at www.github.com/andrewcparnell/NAO.
The notation we use is as follows: • NAO(t) is the North Atlantic Oscillation value at time t. The goal of our model is to estimate NAO(t), and its uncertainty, for a set of chosen times.
• We superscript both the quantities above with m and f so that XRF m refers to the modern XRF data set, with associated known NAO m , and XRF f refers to the fossil data set, for which we wish to estimate NAO f .
• ti refers to the age of the core at depth i. The ages in our core are all given in years BC/AD.
• θ refers to the set of parameters governing the relationship between the NAO and the XRF measurements, as well as the dynamics of how the NAO changes over time.
Our model proceeds by creating a Bayesian joint posterior distribution: The term on the left-hand side of the equation is the posterior distribution and represents the probability distribution of fossil NAO impacts given the observed data. The terms on the right-hand side represent respectively, the likelihood (the probability distribution of XRF f given NAO f ), the distribution of XRF m given NAO m , and the prior distribution of the parameters governing the relationship between XRF and NAO.
For the distribution of XRF given NAO, we standardise all the XRF values (by chemical element) and fit a multivariate normal polynomial regression model (MVN). This means, for the values k = 1,..., 13 chemical elements, we use: 1,...,μi,13] with μik = β0k + β1k NAO(ti) + β2k NAO(ti) 2 , and  is a covariance matrix which captures the extra dependence between elements not explained by differences in NAO.
We set the prior distribution on NAO f as a continuous-time random walk, which should reasonably match climate behaviour over the reconstructed time period (as in Haslett et al., (2006)). Other choices are available, such as long-memory or long-tailed stochastic processes Parnell et al., (2015) use a Normal inverse Gaussian process. Our prior distribution is: where σ 2 is a parameter representing the variance 297 of the NAO f increments for a unit of time.
Finally, we set uninformative prior distributions on the remaining parameters: where N and U represent Normal and Uniform distributions and I is the identity matrix.
The above model is computationally expensive to fit using the default tools for Bayesian model fitting due to the large number of parameters and the high data dimension. Instead, as stated above, we follow the approach of Parnell et al. (Parnell et al., 2015), which involves a computational approximation to fit the model in three steps. The first step involves fitting the model to the modern data only. The second step involves estimating NAO f for the fossil layers, and the third step involves constraining the estimated fossil NAO f values according to a random walk model.
In the first step of the model, the total overlapping period between modern data, XRF proxy and observed NAO index (i.e., NAO m and XRF m , respectively), extends from 1825 until 2012 AD. However, the XRF data from 1825 until 1930 AD has a lower resolution (i.e., decadal) than the 1930-2012 AD period, owing to the usual slight decrease in the age-depth model accuracy. Thus we restricted the overlapping fitting period employed in the analysis to 1930 -2012 AD. XRF m data were resampled with a yearly resolution for the overlapping period using the R function "approxTime" from the package "simecol" (Petzoldt and Rinke, 2007). We fit the model in R(R Core Team, 2016) and use the JAGS software (Plummer, 2003) (Just Another Gibbs Sampler).
The performance of the fitting algorithm can be determined by looking at the Brooks-Gelman Rubin (R̂) statistic (Brooks and Gelman, 1998;Gelman and Rubin, 1992).
We run the algorithm until all R̂ values are less than 1.05, which indicates satisfactory convergence of the algorithm to the posterior distribution.
We evaluate the performance of the model by testing the predictive performance of modern relationship between NAO m and XRF m (Step 1 as outlined above). We created predicted NAO values from the modern model using the modern NAO data (i.e, NAO index data). If the model is estimating NAO correctly there should be only 5% of observations outside the 95% interval, and 50% outside the 50% interval. Finally, the complete impact reconstruction is created using a 10-year time grid, and includes both 95% and 50% uncertainty intervals.

2.4-Statistical analyses
All NAO reconstructions have been converted to decadal time-scales to facilitate the comparison. For each reconstruction, all NAO values within the same decade have been averaged and we use that average value for its particular decade. The magnitude of the relationships between NAOs, sunspot numbers and volcanic eruptions were obtained according to Spearman's rank correlation coefficients (ρ) and associated p values. Unless otherwise stated, significance (p-value) is always considered at values of p < 0.01. The statistical treatment of the data was performed with R software (R Core Team, 2017).

3-A new local NAO reconstruction: Central Iberian Peninsula
The NAO has a great effect on winter climate in the Iberian Peninsula (IP) (Comas-Bru and Hernández, 2018;Hernández et al., 2015;Trigo et al., 2008). In particular, highmountain lakes from the IP are highly influenced by the NAO, because the cold and wet conditions during negative NAO phases control annual ice-cover dynamics (i.e., freezing and thawing) via the interaction between air temperature and precipitation (Sánchez-López et al., 2015). While the NAO is particularly relevant during the boreal winter, its impact on ecosystem and ice-cover dynamics is not restricted to this season (Gouveia et al., 2008); the NAO signal that is captured in lake records from this region spans from January to May (Sánchez-López et al., 2015).
A previous study (Sánchez-López et al., 2016) using geochemical and mineralogical data from the Cimera Lake record (40°15'N -5°18'W, 2,140 m a.s.l.) revealed that climatic conditions, highly influenced by the NAO, are preserved in the sediments via the frequency of intense runoff episodes caused by rain-on-snow events as well as by lake productivity, which is primarily governed by the ice-cover duration. Here we use this data to reconstruct the NAO impact for the central IP (NAOIP) for the last two millennia using a Bayesian modelling approach (see Methods).
The reconstructed local NAO impact ranges between -2 and 2 and represents the quantitative impact of NAO for the central IP (NAOIP; Fig. 2). Only 2.9% of the observations fall outside the 95% confidence interval (Fig. 3). These results indicate satisfactory performance of the model and validate the NAOIP reconstruction (see methods). The NAOIP has not been reconstructed for the period c. 1200 -1260 CE due to the lack of proxy data.
Though the model permits such interpolation, the uncertainties would be too wide to enable any reasonable interpretation.

4-Comparison with previous NAO proxy-based records
The comparison between different proxy-based NAO reconstructions published in the last two decades (Table 1) points out a number of periods with consistent signals as well as some with notable differences (Fig. 4). All reconstructions show a similar long-term evolution with high positive NAO values during the MCA and lower positive or negative NAO values during the LIA (Fig. 4). This coherence across regions suggests an hemispheric imprint of this climate mode at low frequencies compared to the local impacts that might be recorded by each reconstruction over shorter timescales.
To establish the extent to which the multi-decadal variability between positive and negative excursions in NAOIP during the last 2 ka compares to previously published reconstructions, we calculated Spearman's Rank Correlation Coefficients for decadal timescales (Table 2 Olsen et al., (2012;NAOOLS) or negative (NAOBAK and NAOFAU) -for those NAOs which cover these periods. These NAO differences might be attributed to the type, number, location and sensitivity of the archives used, the statistical approaches and external and/or internal forcings affecting each location in different ways.
Unlike gridded instrumental-based NAOs arising from datasets (e.g., atmospheric sealevel pressure) covering a vast geographical region (e.g., the North Atlantic region), the proxy-based NAO records are usually inferred using only a few climate proxies from restricted regions or locations (e.g., NAOTRO, NAOBAK, NAO by Sjolte et al., (2018;NAOSJO), or are based on a single paleorecord, as is the case of NAOOLS, NAOFAU and NAOIP. This restricted geographical constraint would imply that the reconstructed NAO signal is attributed to the local response of climate to this mode of variability rather than reflect a regional signal. Alternatively, the NAO by Cook et al., (2019;NAOCOOK) NAO by Luterbacher et al., (2001;NAOLUT) and NAOORT should be considered differently since they were reconstructed using a larger number of archives and have a wider spatial representation; i.e., a regional assemblage of local NAO impacts.
Nevertheless, we show here (Fig. 5) that the local impact of the NAOIP is capable of reflecting a much wider regional signal than what could be foreseen before. Thus, the reconstructed signal may be more accurate than in multi-archive reconstructions of regional NAOs. The reason for this is that the latter include different proxies from different archives that usually represent different seasons and climate variables, and, therefore, include a mixture of signals that will result in a smoother reconstruction.
Equally importantly, it is now clear that the NAO is far from unique in controlling the large-scale atmospheric variability in the North-Atlantic European region, acting simultaneously with other modes of variability that also play an important role in this region-namely the EA and SCA . In particular, Comas-Bru and McDermott (Comas-Bru and McDermott, 2014) previously showed that different NAO/EA and NAO/SCA combinations influence European winter temperatures and precipitation as a consequence of sea level pressure (SLP) dipole migrations in the North Atlantic. In particular, this is expressed by a more homogeneous spatial pattern in temperature and precipitation over the IP for periods with a predominance of the same sign for both the NAO and the EA (i.e., MCA and LIA) than for periods when these modes are of opposite sign (i.e., RP and EMA) (Sánchez-López et al., 2016). It is noticeable that the geographical displacement of the northern centre of action of the North Atlantic SLP dipole is relatively smaller than the concomitant migration of its southern pole during winters-where the NAO/EA have the same sign compared to years of opposite sign (Fig. 6). This could potentially explain the latitudinal variability observed for NAO reconstructions over the last millennia. Results also demonstrate that differences in precipitation between years of NAO and EA with the same and opposite signs (Fig. 7) are generally small for mid-latitudes (e.g., IP and Mediterranean basin) and larger for high-latitudes (e.g., Greenland, Ireland, and the UK). On the contrary, differences in temperature are larger for mid-latitudes and smaller for high-latitudes. These differences are almost inexistent in the western Atlantic sector (US and Canada) where the NAO impact is also significantly weaker. Hence, regional NAO reconstructions which use different archives (e.g., tree-rings, speleothems, ice-cores) recording different climate variables (e.g., precipitation, temperature) from different locations (high-vs mid-latitudes) could experience some kind of bias due to EA interactions.

5-Volcanic eruption impact
Despite some previous proxy-based NAO reconstructions show a robust positive NAO response in the 4 -5 years following the major eruptions of the last millennium Sjolte et al., 2018), a recent review of the impact of explosive volcanic eruptions on the main climate variability modes determined that no firm conclusions can be drawn at the moment (Swingedouw et al., 2017).
We compared NAOIP with volcanic eruptions (n=42) responsible for the largest accumulation of sulphate (> 20 kg·km -2 ) over Greenland according to Sigl et al., (2015) for the last 2 ka for decadal timescales (Table 3). We find positive reconstructed NAO values for c. 40% of the decades after these eruptions. This reinforces previous works Sjolte et al., 2018) that do not reach compelling conclusions on the relationship between NAO and volcanic activity. Nevertheless, if we relax the minimum emissions threshold and consider all Northern Hemispheric eruptions with an accumulation of sulphates larger than 3 kg·km -2 (n = 93), we find that decadal NAO values after these eruptions are positive in the 60% of the cases. Moreover, to explore the evolution of the volcanic impact on the NAOIP, and its relationship with other external forcings (i.e., solar activity), the studied interval (150 BCE /Present) was divided into two millennial subperiods (higher/lower variance of solar activity): S1 (150 BCE -1000 CE) and S2 (1000 -2012 CE). Decades with positive values after NH eruptions during S1 (n = 52) total 83%, whereas for S2 (n = 41) this decreases to only 34%. Finally, if we compare the NAOIP to decades with more than one eruption (n = 18), we find that 72% of them reconstruct positive NAO values. This suggests that the sum of the eruptions has more impact than a single large eruption. A further comparison including all the NAO reconstructions shows a wide range of percentages (37 -91%) of NAO+ decadal values agreeing with the largest NH volcanic eruptions (Table 3).
Hence, there is an apparent influence from the occurrence of NH volcanic eruptions and the preferred signal of the NAO pattern over decadal timescales.

6-Solar forcing modulation
We also compared the NAOIP with a decadal sunspot number reconstruction (Usoskin et al., 2016) through the CE. To avoid volcanic eruption interferences in the analysis, the decades corresponding to the fifteen largest volcanic eruptions globally were identified ( Fig.   2 and Table S1) and, considered to be outliers, were removed from the solar forcing analysis.
The linear correlation between the decadal NAOIP index and the sunspot number is nonsignificant (= 0.20;p > 0.1;DF = 196). However, the graphical representation of this relationship clearly displays an inflection point around the occurrence of 40 sunspots (Fig.   8). If we only consider sunspot numbers under 40, the linear correlation with NAOIP becomes significant (= 0.56; p < 0.01; DF = 103). Additionally, to explore the stationarity of the solar influence on the NAO evolution, the studied interval (150 BCE -Present) was divided into the same two millennial subperiods considered previously: S1 (150 BCE -1000 CE) and S2 (1000 -2012 CE). The sunspot reconstruction during S1 displays fewer oscillations than during S2, and only includes one Grand solar minimum (Fig. 2). NAOIP does not exhibit a significant correlation (= 0.29) for this period (S1). In contrast, during the S2, the sunspot exhibits a more oscillatory pattern with up to six reconstructed grand solar minima (Fig. 2).
For this period, the NAOIP vs Sunspot present a correlation value of = 0.58; p < 0.01; DF = 46 ( Fig. 2 and 8; Table 4). This correlation is evident in the last 600 years (Fig. 2), where the NAOIP shows a good agreement with the three latest Grand solar minima (i.e., Spörer, Maunder, and Dalton). To evaluate this relationship, we have reproduced the same analysis with all the NAO reconstructions (Table 4 and Fig. S1). Results are similar to those obtained for the NAOip exhibiting no significant correlations for the unfiltered data (i.e., using all sunspot numbers).
However, the sunspot numbers under 40 also present a tipping point (Fig. S1)  and NAO for mid-latitudes. These latitudinal differences of solar activity impact can be attributed, among other criteria previously exposed (e.g., type of proxy and methodology), to the fact that the solar variability signal is not uniformly distributed (Lean, 2017). Annual and decadal variations in solar activity have the largest impacts in the mid-latitudes (Haigh, 2011). Previous analyses of surface air temperatures (Lean and Rind, 2008) have already shown preferential warming in regions approximately in the range 30 -60° latitude for both hemispheres.

7-Discussion and conclusions
This study presents a new multi-annual proxy-based NAO reconstruction using a Bayesian To disentangle the reasons underpinning the observed disagreements between the NAO reconstructions throughout time, we analysed their regional distribution and the available records, as well as the distinct impact of potential external forcings (i.e., volcanic eruptions and solar activity).
Our results demonstrate a latitudinal solar activity impact on the NAO variability over Besides the influence of these two external forcing mechanisms we also assessed the role of an internal mechanism-namely the interaction of NAO with the second most important large-scale pattern of atmospheric circulation in the North-Atlantic European sector, i.e., the Eastern Atlantic (EA) pattern. It is well known that combinations of NAO and EA phases can change the geographical position of the NAO centres of action and affect the strength and latitudinal location of the dominant westerlies entering Europe from the Atlantic. As a consequence of this differential NAO impact on these climate variables according to the state of the EA, the chosen proxy and its location are crucial to more accurately reconstruct the NAO. Moreover, we are firmly convinced that this differential NAO impact could explain the differences between the available NAO reconstructions according to the geographical location of their proxies, the selected time intervals and the climate variables that they are sensitive to. Although a wide regional distribution of records could apparently yield better results, the antagonistic impact of the combined NAO and EA modes on some climate variables over the mid-and high-latitude records could be masking or, even, cancelling out the actual NAO pattern.
It appears obvious that further studies are required to better understand the NAO's behaviour and the disagreements between the continuously increasing number of available NAO reconstructions. Regional NAO reconstructions like the ones derived by integrating a grid of instrumental or proxy-based regional data can be considered more robust-and can aid the understanding of general climate dynamics-only if the records employed are sensitive to the same forcing and therefore can capture the same signal. In contrast, local NAO reconstructions would be more useful to determine NAO impacts of local meteorological variables, which in turn influence the local ecosystems (and natural archives), economy, and society. Therefore, local NAO reconstructions can help to develop better mitigation policies against problems derived from NAO climatic effects such as agricultural yield or water scarcity. We also argue that new applied statistical approaches (i.e., Bayesian) would help to improve the reliability of the results. Overall, we stand for a more adaptable concept of proxy-based NAO using appropriate distinctions, since it is impossible to understand the NAO as a single pattern. Rather, it should be regarded as complex system that is controlled by multiple factors, some of which are stochastic and are therefore difficult to constrain.