Extended North Atlantic Oscillation and Greenland Blocking Indices 1800–2020 from New Meteorological Reanalysis

: Based on newly-available meteorological reanalysis, we compile and present extended seasonal series of the North Atlantic Oscillation (NAO) and Greenland Blocking indices spanning 1800–2020, which we analyse for evidence of signiﬁcant trends. This represents a major backward extension of the previously available instrumental-/reanalysis-based Azores–Iceland and principal component-based NAO indices, and allows us to evaluate the potential effect of natural climate perturbations, especially the 1809 and 1815 major volcanic eruptions and ~1790s–1830 Dalton solar minimum, on North Atlantic atmospheric circulation. We ﬁnd that winters 1809/10 and 1816/17 mark positive NAO peaks, relative to several years before and afterwards, which is in accordance with the theory of volcanic forcing of climate. However, there is little evidence of a summer NAO volcanic signature. Overall, based on the signiﬁcantly longer new reanalysis time series, the new series presented here corroborate and extend our previous results of: (1) a signiﬁcantly more variable year-to-year NAO with a recent exceptional clustering of extreme events since 2000 for winter; (2) a signiﬁcant increasing trend in blocking over Greenland in summer. These trends have major repercussions for the probability of the occurrence of extreme weather events over northwest Europe and for the sensitivity and response of the Greenland Ice Sheet to global warming, especially if they continue as an integral part of anthropogenic climate change.


Introduction
The North Atlantic Oscillation (NAO) is a widely-used measure of North Atlantic polar jet stream variation (mainly the north-south position of the jet), which is an important control of climatic conditions and extreme weather fluctuations across the UK, Northwest Europe and eastern North America [1][2][3][4]. The NAO is variously measured using a south minus north mean sea-level pressure difference (typically Azores-Iceland but sometimes using Lisbon or Gibraltar as the southern station) or a principal component analysis of gridded pressure or geopotential height data over the North Atlantic [1, 2,[5][6][7]. Representing a complex, non-linear dynamical system, the NAO has a broad spectrum of temporal variations [8][9][10] and, apart from internal variability, is likely to be influenced by a number of jet stream drivers [11]. Potential drivers of the jet stream and NAO variability can be broadly grouped into several categories: cryospheric effects from variations in seaice extent and snow cover; oceanic effects from North Atlantic sea-surface temperatures, which may also, for example, influence the NAO teleconnection with Tibetan Plateau surface air temperatures [12]; tropical influences such as the El-Niño Southern Oscillation; stratospheric effects due to stratospheric circulation variability, solar variability, volcanic In previous work, we used contemporary meteorological reanalysis data to construct monthly and seasonal GBI time series spanning from 1851-2015, finding a significant increase in summer GBI since 2007 [30]. This Greenland Blocking increase was subsequently confirmed by other studies [21,22,32,34,40]. In a parallel analysis, we analysed trends and variability in NAO from 1899-2013 and reported a significant decrease in the summer NAO since the 1990s as well as a significant long-term increase in the year-to-year variability of early winter (December) NAO values, which culminated in a clustering of record high and low values during the last decade of 2004-2013 [7]. However, these analyses are limited by the available lengths of record, so it is important to test their results using newly-available longer and potentially more reliable climate datasets.
Recent efforts by the main climate modelling centres have produced new generations of meteorological reanalysis. These include products available back to 1806 (the extended Twentieth Century Reanalysis dataset, version 3 = 20CRv3 [42]), to 1600 (the EKF400v2 reanalysis [43]) and even for as long as the last 1000 years (The Last Millennium Reanalysis [44]). Here, we use one of these products, the EKF400v2 reanalyses, as a basis for constructing and analysing new long seasonal records of NAO and GBI covering the period 1800 to 2020. Extending our previously published NAO and GBI analysis back by 50-99 years keeps within the timeframe of at least some instrumental data being available (which is likely to make a reanalysis more reliable than using proxy data alone) and enables us to reassess the significance of North Atlantic atmospheric circulation changes in the last few decades, which are likely to have been strongly influenced by humandriven global warming, in a largely unperturbed (by anthropogenic greenhouse gases) early-mid-nineteenth century context. To bring our analysis up to date and use the best available composite product, we supplement EKF400v2 data with the European Centre for Medium-Range Weather Forecast's (ECMWF) recently available ERA5 reanalysis [45].
Extending the NAO/GBI records back to 1800 also allows us to assess the potential signature on North Atlantic air circulation patterns and variability of a couple of significant natural climate variations that occurred during the early nineteenth century, specifically, the Dalton solar minimum of 1790-1830 [46] and the April 1815 Mt. Tambora volcanic eruption, which was the most powerful such event in human history [47]. Previous modelling work suggests that, at a global or hemispheric level, volcanic effects from multiple eruptions during this period, rather than solar forcing, are thought to have mainly contributed to the global temperature dip during these decades [46]. Major tropical eruptions warm the lower stratosphere at low latitudes through the absorption of solar radiation, which enhances the stratospheric meridional temperature gradient. In turn, this strengthens the polar vortex and favours positive Arctic Oscillation (AO)/NAO conditions with milder, wetter winters in Northern Europe, contemporaneous with global mean surface cooling due to aerosol reflection, in the first couple of years following an eruption. On the other hand, tropospheric radiative cooling dominates in summer, with cooler conditions across much of Europe and a tendency towards a weaker Scandinavian high-pressure system [48][49][50], although the latter is not necessarily reflected through a negative NAO (e.g., [4]). Reference [51] suggests at first a negative and later a positive NAO arising from changes in the atmospheric dust veil distribution and meridional heating gradient following climatically-significant tropical eruptions. Meanwhile, periods of reduced solar output are thought to weaken the polar vortex and lead to weaker polar westerlies and a weaker North Atlantic polar jet stream in winter [52,53]. Despite earlier modelling studies, intriguing questions remain concerning the relative effects of the Dalton solar minimum and the Tambora volcanic eruption on the Euro-North Atlantic (and global) climate during the mid-late 1810s. Since previous purely instrumental-or reanalysis-based NAO and GBI values have not previously been available before 1821 and 1851, respectively (and 1899 for principal component-based NAO records), we also address this issue here.
Finally, it has recently been reported [54] that NAO-AO relations may break down under climate change due to surface temperature anomalies over Eurasia and the Pacific Atmosphere 2022, 13, 436 4 of 28 and changes in troposphere-stratosphere coupling, and our extended NAO record may be useful for addressing this and other questions in the future.

Reanalysis Datasets
We combine NAO and Greenland geopotential height time series from two recently available pioneering reanalyses: the EKF400v2 and ERA5. The EKF400v2 reanalysis uses the ECHAM5.4 global climate model with ensemble Kalman filtering to assimilate instrumental (pressure, temperature and precipitation) and proxy (tree ring and coral) records, to produce monthly means spanning 1601-2003 at 2 • latitude by 2 • longitude resolution with 32 vertical levels [43]. An ensemble of pre-computed transient model simulations are blended offline with the above-mentioned observational data. Full methodological details are provided in [43]. EKF400(v2) is the first paleo-reanalysis with multivariate monthly information on the state of the atmosphere for the past 400 years.  [43]. Here, we use EKF400v2 data for the period 1800-1949; 1800 is an appropriate early cut-off date, as indicated by significantly increased ensemble spread (divergence of individual model members) reflecting restrictions in input data before the late eighteenth century.
Our modern reanalysis of choice, ERA5, covering the period January 1950 to nearpresent, replaces the older ERA-Interim analysis which stopped production on 31 August 2019, and assimilates vast amounts of historical in situ and satellite data [45]. ERA5 is based on ECMWF's Integrated Forecasting System (IFS) Cy41r2 that was operational in 2016, uses an updated assimilation scheme, has much higher resolution in time and space, with hourly analysis fields on a 31-km grid (compared with 6-hourly and 79 km for the preceding ERA-I reanalysis), and has 137 vertical levels up to 80 km altitude. As with EKF400v2, ERA5 also ingests newly re-processed extra datasets. Another innovation relative to ERA-I is uncertainty estimates, although, because we are highly confident about reanalysed NAO and GBI values for the last few decades, we do not use these here.
For comparison with EKF400v2 results and benchmarking based on weather station data, we use 20CRv3, which is the fourth-generation (following v1, v2 and v2c) product of the Twentieth Century Reanalysis that covers the period 1806-2015, which represents a markedly longer period than 20CRv2c (1851-2012) [42]. The 20CRv3 reanalysis assimilates only surface pressure data into the coupled Global Forecast System land-ocean model with prescribed sea-surface temperatures and sea ice concentration, and comprises 3-hourly estimates of meteorological variations on a 75-km grid with 64 vertical levels. The previous version of this dataset (20CRv2c) has several significant issues, including a global sea level pressure bias in the mid-19th century and problems with ensemble-based estimates of confidence. 20CRv3 uses upgraded data assimilation methods and assimilates a greater number of pressure observations, as well as having a newer, higher-resolution forecast model that specifies dry air mass. These improvements have reduced the sea-level pressure bias and helped to constrain error estimates. However, our checking of International Surface Pressure Databank (ISPD, [55]) archives does not show the presence of Icelandic surface pressure data until November 1845, in both available and assimilated data ( Figure S1). This is somewhat surprising given the known presence of daily pressure observations from Southwest Iceland (Reykjavík and Stykkishólmur) back to 1820 [5,56]. Not using Icelandic barometric pressure data before 1845 could influence the representation of NAO and GBI in 20CRv3, by potentially suppressing variability in extreme seasons and increasing the ensemble spread due to poorer constraining of the North Atlantic atmospheric circulation patterns. Furthermore, Arctic sea ice switches from a repeat climatology (1806-1849) based  [58] dataset [42], which could also impact the representation of high-latitude North Atlantic climatology.
Uncertainty in these products is typically expressed by means of a model ensemble, of which there are, respectively, 30 and 80 ensemble members for EKF400v2 and 20CRv3. Here, we mainly use ensemble mean statistics but also use the EKF400v2 ensemble member series to help evaluate the significance of NAO and GBI trends (Section 3.1). We also use ensemble member data to help delimit the uncertainties of individual yearly GBI/NAO data for 1800-1830 (Section 3.2).

Station-Based Mean Sea-Level Pressure Series and Their Use in Evaluating Reanalysis Data
Very few mean sea-level pressure (MSLP) series are available for before the midnineteenth century. Here, we introduce backward extensions to the Southwest Iceland and Gibraltar MSLP series [5] and updated versions of the London and Paris MSLP series [59] to compare with/validate 20CRv3 and EKF400v2 reanalysis MSLP. These are four of the longest-running, nearly continuous MSLP station series, and between them cover the eastern North Atlantic/western Europe zone of interest.
The monthly version of the Southwest Iceland pressure series, which started in 1821, has been extended here to cover the period September 1807 to August 1814 (with a remaining gap from 1814-1821), based on observations made at Akureyri (65 • 40 N, 18 • 05 W; 5 m above sea level) by a coastal surveying team [60]. The original journals have not been found but they were all copied by a member of the team in the early 1830s. These copies (now to be found at the national library of Iceland) were copied again by the Danish Meteorological Institute in the late 19th century. These later copies were digitised by Dr. Trausti Jonssón (pers. comm. 2021) who applied necessary corrections and supplied accompanying metadata. The original French inches (27.07 mm) were converted to hPa. The barometer was located in an unheated room, and its temperature was not available but was estimated. The relationship between the ambient temperature at the barometer (ttt b ) and outside temperature (ttt o ) was based on comparative measurements at Nes near Reykjavík made in 1822-1829: It was assumed that in the case of very low outside temperatures (ttt o < −6.5 • C), the inside temperature did not fall below −3 • C. The Celsius freezing point was used as a reference for the correction which was 0.16917 hPa/ • C. Adjustments were made for height above sea level (+0.6 hPa) and gravity (+1.5 hPa). Based on high correlations (0.91 ≤ r ≤ 0.98) between monthly means of the pressure in Akureyri and Reykjavík for 1951-2020, we used mean monthly MSLP differences between Akureyri and Reykjavík MSLP data for this more modern period to correct the 1807-1814 Akureyri monthly MSLP values.
The Gibraltar MSLP monthly dataset from 1821-2020 ( [5], updated data) was extended back to 1799 using nearby Cadiz/San Fernando data: (1) from January 1799 to August 1801, 1804 to 1809 (with a few monthly gaps) and May 1812 to December 1813 digitised by Dr. Deborah Smith (pers. comm.), and (2) from September 1816 to December 1819 and October-December 1820 digitised by co-author MB. Similar to the old Icelandic pressure data described above, these newly unearthed San Fernando pressure data had corrections applied for altitude (~3 hPa) but, due to limitations in metadata, only part (1) was corrected for gravity (−0.88 hPa) and temperature (using the formula that fully corrected MSLP = partly corrected pressure−{1.65 + [temperature ( • C) − 10] × 0.16}).
We use monthly London and Paris MSLP updated from [59]. The London and Paris MSLP series used in this analysis are based on the daily series described in [61,62]. However, rather than being a combination of sub-daily observations and 24-h mean values, as was the case in the previous series, the data values in the new version are sub-daily throughout. This required the digitization of additional data to complete certain periods. In the case of London for the period 1843-1854, the 9 a.m. and 2 p.m. readings recorded in the private weather diary of John Henry Belville, taken at his home at various sites around Greenwich, were used; from 1854-1920 the manual barometer readings taken at Royal Observatory Greenwich at 9 a.m., 12 p.m. and 3 p.m. were used; from 1920-1949 (and to fill gaps in earlier periods) London observations from the UK Meteorological Office Daily Weather Report (DWR) were used [63]. Gaps in the Belville series were completed using observations kept by William Rogerson, as published in the Nautical Magazine. These new data supplement the existing series used in the earlier version of the London MSLP series. In the case of the Paris series, data from the Montsouris Observatory were used for the period 1872-present. To these data were added newly digitized data recorded at the Paris Observatory that had been extracted from the Journal de Physique publication prior to 1816 and from the Observatory's register for the period 1816-1871.
Adjustments were applied where necessary to the digitised barometric pressure series to reduce the readings to zero degrees Celsius, to standard gravity and to mean sea-level from their respective observation heights. The readings were also converted to the unit of hPa from the original units of observation. The data were quality-controlled by eliminating un-physical values, through the comparison of the series against each other as well as the contemporary sub-daily MSLP values recorded at De Bilt in the Netherlands, and against the 20CRv3 values after the method described by reference [64]. Large, rapid changes in daily values of MSLP were also manually examined and removed where necessary after reference [65]. The monthly data were then tested for temporal homogeneity through the comparison against the weighted average of the closest neighbouring stations using the Penalized Maximal t test [66] and without a reference series using the Penalized Maximal F test [67] from the RHtestsv4 software [68]. The reference series were the same as used in references [61,62]. However, given the improvements made through the digitization of the new data sources, the only adjustment deemed necessary was to account for drift in the London series prior to 1822 [61]. An adjustment of +0.5 hPa was applied to the Paris data before 1856 to correct for a step-change in the data, in addition to a correction of +3/4 line (2.3 hPa) before 21 May 1800. The sub-daily values were interpolated to a noon observation from the respective recording time, and monthly mean values were calculated from these noon values.
Comparative statistics and visuals of seasonal 20CRv3 and EKF400v2 reanalysis data against these four observed time series for 1806-1850 are presented in Table S1 and Figures S2-S5. They show a generally superior performance of EKF400v2 relative to 20CRv3 for this early-mid-nineteenth century period. The highest EKF400v2-weather station correlations (0.88 ≤ r ≤ 0.98) are found for London and Paris (Table S1), which is unsurprising, as these stations were used in the EKF400v2 reanalysis (Valler et al. 2021). 20CRv3-weather station correlations for these two sites are still high (0.67 ≤ r ≤ 0.89) for three seasons, although they drop off for summer (0.50 for London and 0.46 for Paris). The ISPD, version 4 [69], used as a basis for 20CRv3, contains time series for 1811-1880 for Paris but only 1815-1817 for London. Reanalysis-station correlations for SW Iceland are much higher for EKF400v2 (0.81 ≤ r ≤ 0.90) than 20CRv3 (0.23 ≤ r ≤ 0.55), which reflects the inclusion of early (1821-1844) SW Iceland MSLP data in GHCN-Monthly v2 [70], which is ingested into EKF400v2, and its lack of incorporation into 20CRv3 until 1845. Similarly, Gibraltar is not incorporated into 20CRv3 until 1852, compared with EKF400v2, which ingests Gibraltar from 1821 via the GHCN v2 (https://iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NCDC/.GHCN/.v2/, accessed on 2 March 2022). While EKF400v2 shows overall better performance than 20CRv3 when both are compared with observations for the Gibraltar grid-point, EKF400v2 performs poorly for summer correlations, while its mean errors are slightly greater in winter. However, overall, based on comparisons for these four locations and with generally much lower RMSE/MAE error values for EKF400v2 MSLP in the zone of interest compared with 20CRv3 (Table S1), and evidence of artificially suppressed year-to-year variability in 20CRv3 pre-1850 ( Figures S2-S5), we use EKF400v2 in the following analysis.

NAO and GBI Series
To maximise temporal consistency of the NAO and GBI time series over the whole time period, we spliced EKF400v2 and ERA5 mean sea-level pressure (NAO) and geopotential height (GBI) monthly datasets (i.e., adjusted the mean and standard deviation for the respective month so that they matched) against each other for the common EKF400v2/ERA5 overlap period of 1950-2003, before aggregating these monthly time series to seasonal before our analysis. Seasonal mean conditions and spatial anomalies of the patterns of 500 hPa geopotential height over the Greenland region and of North Atlantic mean sea-level pressure for different periods based on our spliced series are shown in Figures 1 and 2.
We use three versions of the NAO in this study to account for potential differences between methodologies and underlying data. The first NAO, hereafter NAO Fixed , represents a pseudo-station NAO using EKF400v2 and ERA5 sea level pressure from the reanalysis grid boxes encompassing Ponta Delgada, Azores and Reykjavik, Iceland. The Azores and Iceland time series are used to create the NAO for all 30 ensemble members of the spliced EKF400v2-ERA5 reanalysis. Both the Azores and Iceland MSLP seasonal time series are normalised using their respective means and standard deviations for 1951 to 2000 before the Iceland series is subtracted from the Azores. The Azores-Icelandic NAO is thought to better capture the annual migration of the NAO centres of action than using a continental southern station such as Lisbon or Gibraltar [5] (however, a British Isles-Greenland dipole is often identified as a more reflective Summer NAO spatial pattern [71]).
Atmosphere 2022, 13, x FOR PEER REVIEW 7 of 28 observations for the Gibraltar grid-point, EKF400v2 performs poorly for summer correlations, while its mean errors are slightly greater in winter. However, overall, based on comparisons for these four locations and with generally much lower RMSE/MAE error values for EKF400v2 MSLP in the zone of interest compared with 20CRv3 (Table S1), and evidence of artificially suppressed year-to-year variability in 20CRv3 pre-1850 (Figures S2-S5), we use EKF400v2 in the following analysis.

NAO and GBI Series
To maximise temporal consistency of the NAO and GBI time series over the whole time period, we spliced EKF400v2 and ERA5 mean sea-level pressure (NAO) and geopotential height (GBI) monthly datasets (i.e., adjusted the mean and standard deviation for the respective month so that they matched) against each other for the common EKF400v2/ERA5 overlap period of 1950-2003, before aggregating these monthly time series to seasonal before our analysis. Seasonal mean conditions and spatial anomalies of the patterns of 500 hPa geopotential height over the Greenland region and of North Atlantic mean sea-level pressure for different periods based on our spliced series are shown in Figures 1 and 2.  We use three versions of the NAO in this study to account for potential differences between methodologies and underlying data. The first NAO, hereafter NAOFixed, represents a pseudo-station NAO using EKF400v2 and ERA5 sea level pressure from the reanalysis grid boxes encompassing Ponta Delgada, Azores and Reykjavik, Iceland. The Azores and Iceland time series are used to create the NAO for all 30 ensemble members of the spliced EKF400v2-ERA5 reanalysis. Both the Azores and Iceland MSLP seasonal time series are normalised using their respective means and standard deviations for 1951 to 2000 before the Iceland series is subtracted from the Azores. The Azores-Icelandic NAO is thought to better capture the annual migration of the NAO centres of action than using a continental southern station such as Lisbon or Gibraltar [5] (however, a British Isles-Greenland dipole is often identified as a more reflective Summer NAO spatial pattern [71]).
To better capture the summer NAO pattern, the second NAO, hereafter NAOPC, uses the first principal component of detrended, monthly (aggregated to seasonal) SLP anomalies from the spliced EKF400v2-ERA5 dataset (1800-2020, anomalies using 1961-1990) covering −90° W-40° E, 20° N-80° N. We construct the NAOPC from individual ensemble members, then average the time series (inverting the time series if necessary so the Azores dipole, or British dipole during summer, is the positive one). The explained variance varies for each ensemble member, but the mean across all ensemble members (i.e., NAOPC) To better capture the summer NAO pattern, the second NAO, hereafter NAO PC , uses the first principal component of detrended, monthly (aggregated to seasonal) SLP anomalies from the spliced EKF400v2-ERA5 dataset (1800-2020, anomalies using 1961-1990) covering −90 • W-40 • E, 20 • N-80 • N. We construct the NAO PC from individual ensemble members, then average the time series (inverting the time series if necessary so the Azores dipole, or British dipole during summer, is the positive one). The explained variance varies for each ensemble member, but the mean across all ensemble members (i.e., NAO PC ) explains 34% (DJF), 32% (MAM), 30% (JJA) and 26% (SON) of the underlying SLP variability. NAO Fixed and NAO PC are generally well correlated but less strongly, and occasionally insignificantly, in summer compared with winter (Table S2), which reflects seasonal shifts in the NAO centres of action [71].
The final NAO, hereafter NAO Luter [72] is used as a benchmark for comparing our 1800-1850 NAO Fixed and NAO PC values. NAO Luter is an Azores-Iceland NAO series that was built based on proxy predictors from Eurasian instrumental (air pressure, precipitation and temperature) and documentary (cloud cover, snow and sea ice and phenological observations) sources. It is the only available (partly) instrumental-based Azores-Iceland record that extends before 1850, and provides a useful cross-reference in our evaluation of NAO seasonal values during 1800-1830 that may have been influenced by the 1815 Tambora volcanic eruption and/or the Dalton solar minimum.
For GBI, we use two measures, where GB 1 is simply the normalised area-weighted 500 hPa geopotential height over the region 60-80 • N, 20-80 • W [17,18] and GB 2 is the areaweighted 500 hPa geopotential height over this Greenland region minus that calculated over the entire hemispheric zonal band of 60-80 • N [34]. GB 2 is sometimes preferred to mitigate the possible effect of global warming on mean height levels, as discussed in Section 1. GB 1 and GB 2 are very highly correlated for all seasons and periods (Table S3). We normalise GB 1 and GB 2 values using the 1951-2000 period.

Statistical Analyses
We use a 21-point binomial filter to emphasize long-term decadal variability in GBI and NAO. Following Hanna et al. (2015), we also use an 11-year running standard deviation to examine potential long-term changes in GBI/NAO year-to-year variability. We use standard significance testing of linear trends, ascribing significance to results where p ≤ 0.05. Standard three-month meteorological seasons are used, where winter (DJF) corresponds to the year marked by the January. To evaluate statistical properties of our new NAO and GBI time series, we evaluate various periods throughout the series but with a focus on two recent periods 2000-2020 and 2007-2020, which mark the recent era of strong Arctic warming and, especially, the record run of low Arctic sea ice summers since 2007 (using two recent periods allows us to assess the time sensitivity of our results). We use the hypergeometric distribution to identify the significance of 21st century NAO and GBI record high/low events, as in [7] following [73]. The hypergeometric distribution is a discrete probability distribution describing the probability of k successes (random draws for which the object drawn has a specified feature) in n draws, without replacement (and is therefore different from the binomial distribution), from a finite population containing K objects with that feature, in which each draw is either a success or a failure. For example, in an N-year record for a particular season, we define 10 record years (extreme high and low GBI or NAO values), N-10 non-record years and 21 years (2000-2020) as sampling "draws". Due to testing multiple seasons and extreme value types (both high and low values), Bonferroni corrections (×8 scaling) are applied to reported p-values to mitigate the increased likelihood of rejecting the null hypothesis of no significant clustering of extreme GBI or NAO seasonal values in the last~20 years of record.

Key Features of 1800-2020 GBI/NAO Series
The winter seasonal time series of normalised GB 1 for the full period of records shows a record peak (2.94 standard deviations above the 1951-2000 winter mean) in 2010 ( Figure 3a). However, in the winter GB 2 series, 2010 is only the fifth highest value, while the three lowest values are in 1984, 2015 and 2018, in contrast to a more even distribution of extreme low values over time in the GB 1 series (Figure 3a). Despite these differences, there is no obvious overall trend in winter GB 1 or GB 2 values or variability over the whole period of records. However, the spatial geopotential height plots show moderate negative anomalies centred over southern and central Greenland, surrounded by moderate positive anomalies, for winter for the 1991-2020 period, with these mostly non-significant anomalies cancelling the wider GBI region (Figure 1). For spring, GB 1 shows the highest value in 2010, while GB 2 has its lowest value in 2011, but again there is no clear long-term trend or overall change in variability with time ( Figure 3b). Spatial anomalies for 1991-2020 for spring are significantly positive and greatest over northwest Greenland and towards the Canadian Arctic (Figure 1), which may reflect greater heating of the near-surface atmosphere from seasonally reduced sea ice cover. The summer GBI graph shows several highly unusual high values since 2007 and extreme year-to-year fluctuations between successive summer seasons, i.e., 2012 (high GBI)-2013 (low GBI) and 2018 (low GBI)-2019 (high GBI) (Figure 3c). These two recent summer pairs have by far the highest such interannual fluctuations in the entire period of the GB 1 record, and two of the three greatest year-to-year fluctuations in the GB 2 record. The overall increase in summer GBI in the last few decades is reflected by the 1991-2020 spatial anomaly map (Figure 1), which shows significantly high geopotential height anomalies over the whole of Greenland and its immediate surroundings but are greatest in west and northwest Greenland and in the Davis Strait. Autumn has the highest (lowest) GB 1 value in 1876 (1837) and the highest (lowest) GB 2 value in 1968 (1938) with no clear pattern in the temporal distribution of extreme values (Figure 3d). An underlying increase in autumn GBI from the late 1970s to the mid-2000s (also shown as significant positive anomalies from 1991-2020 in Figure 1), following a slight decrease since the midnineteenth century, generally reversed in the last 15 years of the records. The most frequent categories of GB 1 values (i.e., the numbers of years with mean GBI values for the respective period) shifted in summer from GB 1 < −0.5 during all the 50-year periods to 2000, to >0.5 in 2000-2020 and >0.5 and >1 (joint equal) in 2007-2020 (Figure 4). The only mean summer GB 1 values exceeding one standard deviation (with the year-to-year standard deviation for the respective period shown in brackets) were 1. 35 (1.28σ) for 2000-2020 and 1.74 (1.38σ) for 2007-2020 ( Figure 4). This pattern of a recent moderate increase in summer GBI is confirmed using GB 2 in Table S4, which suggests that the increase is not solely due to thermodynamics and there is a circulation component. There are no clear overall changes for the other seasons.
summer seasons, i.e., 2012 (high GBI)-2013 (low GBI) and 2018 (low GBI)-2019 (high GBI) (Figure 3c). These two recent summer pairs have by far the highest such interannual fluctuations in the entire period of the GB1 record, and two of the three greatest year-to-year fluctuations in the GB2 record. The overall increase in summer GBI in the last few decades is reflected by the 1991-2020 spatial anomaly map (Figure 1), which shows significantly high geopotential height anomalies over the whole of Greenland and its immediate surroundings but are greatest in west and northwest Greenland and in the Davis Strait. Autumn has the highest (lowest) GB1 value in 1876 (1837) and the highest (lowest) GB2 value in 1968GB2 value in (1938 with no clear pattern in the temporal distribution of extreme values ( Figure  3d). An underlying increase in autumn GBI from the late 1970s to the mid-2000s (also shown as significant positive anomalies from 1991-2020 in Figure 1), following a slight decrease since the mid-nineteenth century, generally reversed in the last 15 years of the records. The most frequent categories of GB1 values (i.e., the numbers of years with mean GBI values for the respective period) shifted in summer from GB1 < −0.5 during all the 50year periods to 2000, to >0.5 in 2000-2020 and >0.5 and >1 (joint equal) in 2007-2020 ( Figure  4). The only mean summer GB1 values exceeding one standard deviation (with the yearto-year standard deviation for the respective period shown in brackets) were 1. 35 (1.28σ) for 2000-2020 and 1.74 (1.38σ) for 2007-2020 ( Figure 4). This pattern of a recent moderate increase in summer GBI is confirmed using GB2 in Table S4, which suggests that the increase is not solely due to thermodynamics and there is a circulation component. There are no clear overall changes for the other seasons. Smoothed series of seasonal GBI ( Figure 5) show major peaks in the mid-late 1960s for winter and around 2010 for summer. The recent summer increase is statistically significant since 1951 in GB 1 but not GB 2 (Table 1). Peak spring GBI values occur around 1930, which is most pronounced for GB 2 . Autumn values are relatively high around the late 1870s and early 2000s (GB 1 only). Low points are seen around 1990 for winter and spring, the early 1920s and (GB 2 ) late 1930s for summer and 1832 for autumn. Peaks in interannual variability in GBI are shown for winter, and especially for spring and summer, in the last decade of the record ( Figure 6). However, linear trends (associated with these recent peaks and other fluctuations) of the overall GB 1 /GB 2 variability time series and selected sub-series are mostly not statistically significant, except for a significant increase in GB 2 variability for spring during 1951-2015 (Table 2). Eight (six) of the ten highest GB 1 (GB 2 ) summer values in the entire record have occurred since 2000 (Table 3). This is very highly significant (p ≤ 0.00045) according to the hypergeometric distribution, with a Bonferroni correction applied to the p-value to account for multiple sampling across four seasons and the two high/low extremes we sample for in Table 3. The top two GBI summers (combined GB 1 and GB 2 rankings) are 2012 and 2019, closely followed by 2011, 2008 and 2010. No similar (significant) recent clustering of extreme GBI years is seen for the other seasons. However, 2010 features in the ten highest GB 1 and GB 2 years for winter, spring (GB 1 only), summer and autumn, while 2018 is in the ten lowest winter and autumn GB 2 (but not GB 1 ) values, and 2013 has the third lowest summer GB 2 value.
Atmosphere 2022, 13, x FOR PEER REVIEW 11 of 28  Smoothed series of seasonal GBI ( Figure 5) show major peaks in the mid-late 1960s for winter and around 2010 for summer. The recent summer increase is statistically significant since 1951 in GB1 but not GB2 (Table 1). Peak spring GBI values occur around 1930, which is most pronounced for GB2. Autumn values are relatively high around the late 1870s and early 2000s (GB1 only). Low points are seen around 1990 for winter and spring, the early 1920s and (GB2) late 1930s for summer and 1832 for autumn. Peaks in interannual variability in GBI are shown for winter, and especially for spring and summer, in the last decade of the record ( Figure 6). However, linear trends (associated with these recent peaks and other fluctuations) of the overall GB1/GB2 variability time series and selected sub-series are mostly not statistically significant, except for a significant increase in GB2 variability for spring during 1951-2015 (Table 2). Eight (six) of the ten highest GB1 (GB2) summer values in the entire record have occurred since 2000 (Table 3). This is very highly significant (p ≤ 0.00045) according to the hypergeometric distribution, with a Bonferroni correction applied to the p-value to account for multiple sampling across four seasons and the two high/low extremes we sample for in Table 3. The top two GBI summers (combined GB1 and GB2 rankings) are 2012 and 2019, closely followed by 2011, 2008 and 2010. No similar (significant) recent clustering of extreme GBI years is seen for the other seasons. However, 2010 features in the ten highest GB1 and GB2 years for winter, spring (GB1 only), summer and autumn, while 2018 is in the ten lowest winter and autumn GB2 (but not GB1) values, and 2013 has the third lowest summer GB2 value.       Table 1 caption. In assessing statistical significance, degrees of freedom were scaled by the smoothing interval (i.e., number of years minus two were divided by eleven). Significant upward temporal overall trends, determined based on reference [74], with the corre-sponding number of significant ensemble trends, are emboldened in red.   NAO seasonal data (Figures 7 and 8) show no overall long-term trends in mean values but exhibit increased year-to-year variability in the last decade in winter (Figure 7a), but other seasons (Figure 7b-d) show the greatest variability at various points in the midlate nineteenth century and early-or mid-twentieth century. The North Atlantic MSLP winter anomaly plot for 1991-2020 shows significant positive anomalies centred around the Azores and significant negative anomalies close to the north of Iceland (Figure 2), which reflect relatively positive NAO values overall during the last few decades of the record (Figure 7a). There are signs of a summer NAO decrease since about 1990 (Figures 7c and 8), confirmed by the North Atlantic MSLP summer anomaly plot for 1991-2020 (Figure 2), which supplements the summer GB 1 increase and follows expectation based on the strong negative correlation between the NAO PC and GB 1 series (Table 4) Table S5.        Interannual variability of NAO seasonal time series ( Figure 10, Table 2) shows a significant increase in winter variability (p = 0.005) for NAOPC during 1901-2015, although, due to the high degree of smoothing, the NAOPC winter variability increase from 1951-  Interannual variability of NAO seasonal time series ( Figure 10, Table 2) shows a significant increase in winter variability (p = 0.005) for NAO PC during 1901-2015, although, due to the high degree of smoothing, the NAO PC winter variability increase from 1951-2015 is not quite significant (p = 0.07). Increases in NAO Fixed winter variability are not statistically significant, although NAO Fixed does become significantly more variable in spring during 1951-2015 (p = 0.021). Secular changes in the year-to-year variability of winter NAO include lower variability around 1850, 1910, 1950 and 2005) (Figure 10a), while spring variability peaks around 1920 and then decreases during the mid-twentieth century before once again increasing to reach a similar or slightly greater peak in about 2010 (Figure 10b). Summer variability peaks around 1820, 1900 and in the late 1950s (Figure 10c). The recent peak (~2010) is slightly higher than the earlier summer peaks for the NAO PC series and similar to or lower than the former peaks for the NAO Fixed record. Table 5 (Figure 10a), while spring variability peaks around 1920 and then decreases during the mid-twentieth century before once again increasing to reach a similar or slightly greater peak in about 2010 (Figure 10b). Summer variability peaks around 1820, 1900 and in the late 1950s ( Figure  10c). The recent peak (~2010) is slightly higher than the earlier summer peaks for the NAOPC series and similar to or lower than the former peaks for the NAOFixed record. Table  5     Detrended correlation coefficients between GB 1 and the two NAO indices are summarised in Table 4. These are nearly all significant and are strongest in winter and weaker in summer. Correlations are generally higher using NAO PC index, while GB 1 -NAO Fixed location correlations in summer are insignificant for 1971-2000. There are no clear long-term trends or variations in these correlations, except for weaker summer GBI and NAO Fixed values clustered in the mid-late twentieth century. Our combined results suggest this may reflect changing centres of action of the NAO during this period rather than a weaker influence of Greenland Blocking on the NAO.

Evaluation of NAO and GBI Values from 1800-1830
This period is of particular interest because it coincides with the Dalton solar minimum as well as the 1815 Tambora volcanic eruption. We therefore next examine summary statistics of NAO and GBI annual values and their variations. Figure 11 shows variations in NAO based on our NAO PC and NAO Fixed reconstructions and NAO Luter . Our NAO series derived from EKF400v2 are significantly correlated with NAO Luter for all seasons, with the strongest correlations for winter and when using NAO Fixed rather than NAO PC for the other three seasons (Table 6). EKF400v2 ensemble mean seasonal NAO mean values and their interannual variability for 1801-1830, together with yearly seasonal values for 1815-1817 (coincident with the immediate aftermath of Tambora) and their variability across ensemble members of EKF400v2, are shown in Table 7. Winters 1815 and 1816 have slightly negative or near-neutral NAO conditions; however, 1817 is the most positive value from the beginning (1800) of both EKF400v2-based NAO winter series until after 1850, with good agreement between the respective NAO indices (Figure 11a, Table 7). Box and whisker plots showing the spread of estimates between EKF400v2 ensemble members confirm a clear shift to a highly positive NAO in winter 1817 (Figure 11a). This is in agreement with the tendency for a major tropical volcanic eruption, such as Tambora, to favour a positive winter NAO [50], and this potential driver may have been augmented by solar activity being near the peak of the 11-year sunspot cycle in 1816 [75], although sunspot activity was then only midrange relative to a generally more active 11-year cycle outside the Dalton Minimum regime. A negative and then positive winter NAO in successive years following Tambora is in agreement with Reference [51]. Our results also show a similar positive NAO excursion in winter 1810, following another major volcanic eruption of unknown location in 1809, which was the third largest such event since 1500 with double the sulphur emission of the 1991 Pinatubo eruption [76] (Figures 11a and 12a).      Summer and autumn NAO values for 1815-1817 are not far from neutral, with the possible exceptions of summer 1817, which NAO Luter suggests was highly positive, and autumn 1816, which was strongly negative in NAO Fixed but not in the other two indices (Figure 11c). The NAO PC summer series indicates little apparent response to either the 1809 (unknown source) or Tambora 1815 volcanic eruptions. Available land-station pressure records and ships' logs indicate a weak Azores High and strong Icelandic Low regime in summer 1816 (the infamous "year without a summer" (e.g., [77])), with low pressure anomalies focused over the UK [50] and cyclonic north-westerly flow giving anomalously wet and cool conditions both there ( [78], pp. 339-340) and over central Europe; however, due to wider changes in the Euro-Atlantic pressure field, this circulation pattern might well correspond to a near-neutral rather than negative NAO, as indicated by [72]. Although fixed-location NAO indices can have problems capturing summer NAO changes (e.g., [2]), a near-neutral summer 1816 NAO value is suggested by all three NAO indices used in our study, including NAO PC , which is not affected by this potential issue. Moreover, summers 1821, 1822 and 1824 are noteworthy as they are highly negative in all three NAO series that we use (Figure 11c). On the basis of the prevalence of sulfate fallout in ice cores in Antarctica and Greenland [79], the Tambora dust veil seems unlikely to have persisted into 1821. However, reference ( [78], pp. 340-342) notes observations of notable coloured sunsets and twilights in London in summer 1821, attributing these to a volcanic dust veil effect, and confirms the impact of the highly negative summer NAO of 1821 on that year's "very cool summer" in the British Isles.
During 1800-1830 (coincident with most of the Dalton Minimum), winter mean NAO values from the three available indices are on average near-neutral (Table 7), and there are a few moderate but unexceptional negative winter NAO excursions during this period (Figures 7a and 11a). Although there is significant ensemble variability in EKF400v2 estimates of NAO for this period, clear differences between some consecutive years are still apparent (Figure 12), and the significant correlation with NAO Luter values ( Table 6) increases confidence in the EKF400v2 central estimates. Although there are no independent data that can help verify reanalysis estimates of Greenland Blocking for this period, our analysis suggests a number of negatively blocked (cyclonic Greenland) seasons during 1815-1817, i.e., winter 1817 GB 1 (GB 2 ) of −1.82 (−1.38), linked with the previously mentioned highly positive NAO, and a similar lack of blocking during the subsequent winter 1818 and in winter 1810; springs 1815 and 1817 were highly cyclonic, summer 1816 was moderately to strongly cyclonic, while autumns 1815 to 1818 were also cyclonic over Greenland.

Conclusions
We have presented newly developed Greenland Blocking and NAO indices based on recently released climate reanalyses, that exceed the span of previous similar records by 34% (221/165 years), including an insightful extension for much of the first few decades of the nineteenth century. Using the individual (n = 30) ensemble members to help assess trends has enabled us to overcome a major statistical limitation of using archived ensemble mean data. Our analysis confirms and extends (based on a longer record) previous results regarding a significant long-term increase in the year-to-year variability of the winter NAO PC , which was previously highlighted for 1899-2014 for the month of December by [7], and here culminates in a significant clustering of extreme low/high values during 2000-2020 in our NAO Fixed series. Whilst detailed physical interpretation is beyond the scope of this paper, we surmise that multiple factors acting on the winter polar jet to move it alternatively north and south appear to have been particularly active in recent decades: these factors are likely to include Arctic sea ice depletion and the recent deep solar minimum of~2010 related to the downturn in the~80-year solar cycle, as well as changes in North Atlantic sea-surface temperature patterns, tropical convection and ENSO [11,13,16].
Increased year-to-year variability of NAO in recent decades has important implications for seasonal-decadal weather prediction and severe weather planning, especially if these trends continue. The interannual variability increase in winter, with a simultaneously, positively-biased NAO in the future could both be sustained if meteorologically extreme events such as the 2010 and 2011 winter seasons interspace generally positive NAO winters, or will decrease in their absence. However, we also highlight that a positive trending NAO with climate warming is far from certain, and there was a predominately negative winter trend observed as recently as the 1991-2014 period [7].
Greenland Blocking, according to both of our extended GBI metrics used here, has had an unusual clustering of record high values in summer since 2007, accompanied by a significant summer increase in GB 1 from 1991-2020, corroborating previous results. The significant recent increase in GBI in summer is also a potentially powerful positive feedback on Greenland Ice Sheet melt that is currently mainly unaccounted for in global and regional climate models [29,34]. After two years of relatively low summertime mass loss (2017-2018), 2019 and then 2021 marked a return to exceptional melting [31,80,81]. Reference [82] suggested that the increased freshening of the North Atlantic that follows higher summer Greenland Ice Sheet mass loss events may lead to a positive NAO. Therefore, an increasingly positive (negative) summer GBI (NAO) could favour a positive winter NAO but further work is needed on a possible feedback. A weakly positive trending NAO is also an expected response to anthropogenic forcing in climate models [83], although the reliability of model projections for atmospheric circulation changes over the Greenland/North Atlantic region has been questioned [29,34].
Due to their magnitude, persistence and uniqueness, it seems unlikely that the significant trends in summer GBI and winter NAO variability and clustering of extreme summer GBI and winter NAO values from 2000-2020 that we report wholly reflect internal variability, and future work should therefore focus on process understanding and improving their representation in climate models. The latter are well attuned to capturing global climate averages but less so regional and temporal (e.g., year-to-year or decadal) variability.
The 1809 (unknown source) and 1815 Tambora volcanic eruptions are expected to have affected the NAO mainly in winter (a positive signal) but also potentially in summer (negative NAO), although here we find greater evidence of the projected winter signal. However, NAO values in summer are less reflective of the overall North Atlantic atmospheric circulation regime, e.g., [4]. We conclude that any Dalton Minimum solar signal, which would again be focused in winter in terms of its possible NAO response, appears to be limited in terms of Euro-Atlantic atmospheric circulation change, although future work should use analysis of NAO reconstructions at the monthly or even weekly timescale(s) to test this result. This aspect merits further attention based on the apparent influence of deep solar minima on Euro-North Atlantic surface air temperatures and circulation patterns since the late nineteenth century [52,84]. A clustering of significant (~−2σ) negative spring values in the Luterbacher et al. [72] NAO series from 1806-1816 inclusive, with positive NAO winters in 1810 and 1817, is reflected for the latter in our EKF400v2-based NAO series, and may well have been influenced by both solar and volcanic forcing during the depths of the solar minimum phase.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/atmos13030436/s1, Figure S1: ISPD observational coverage for selected years in 20CRv3; Figure S2. Seasonal time series of SW Iceland mean sea-level pressure (MSLP) from extended weather station record (grey) and 20CRv3 (orange) and EKF400v2 (blue) reanalyses: (a) winter; (b) spring; (c) summer; (d) autumn; Figure S3. Seasonal time series of Gibraltar mean sea-level pressure (MSLP) from extended weather station record (grey) and 20CRv3 (orange) and EKF400v2 (blue) reanalyses: (a) winter; (b) spring; (c) summer; (d) autumn. Table S1: title; Figure S4. Seasonal time series of London mean sea-level pressure (MSLP) from extended weather station record (grey) and 20CRv3 (orange) and EKF400v2 (blue) reanalyses: (a) winter; (b) spring; (c) summer; (d) autumn; Figure S5. Seasonal time series of Paris mean sea-level pressure (MSLP) from extended weather station record (grey) and 20CRv3 (orange) and EKF400v2 (blue) reanalyses: (a) winter; (b) spring; (c) summer; (d) autumn; Table S1. Correlation coefficients (a), Root Mean Square Error (RMSE) (b) and Mean Absolute Error (MAE) (c) between EKF400v2 (20CRv3) reanalysis and observations seasonal mean MSLP series for 1806-1850. Red highlighted are cases where 20CRv3 is superior; Table S2. De-trended correlation coefficients between seasonal NAOFixed and NAOPC indices for various periods. Both NAO indices are calculated from EKF400v2 and ERA5 reanalysis data spliced together, as discussed in the main text. The shaded value is *not* statistically significant at the p ≤ 0.05 level; Table S3. De-trended correlation coefficients between seasonal GB1 and GB2 indices for various periods; Table S4. GB2 means, standard deviations (stdev) and numbers of extreme values for selected periods. Means greater than 1σ and the highest frequencies for each period/season are emboldened; Table S5. NAOFixed means, standard deviations (stdev) and numbers of extreme values for selected periods. Means greater than 1σ and the highest frequencies for each period/season are emboldened.