Physical and Unphysical Causes of Nonstationarity in the Relationship Between Barents‐Kara Sea Ice and the North Atlantic Oscillation

The role of internal variability in generating an apparent link between autumn Barents‐Kara sea (BKS) ice and the winter North Atlantic Oscillation (NAO) has been intensely debated. In particular, the robustness and causality of the link has been questioned by showing that BKS‐NAO correlations exhibit nonstationarity in both reanalysis and climate model simulations. We show that the lack of ice observations means nonstationarity cannot be confidently assessed using reanalysis prior to 1961. Model simulations are used to corroborate an argument that forced nonstationarity could result from ice edge changes due to global warming. Consequently, the observed change in BKS‐NAO correlations since 1960 might not be purely a result of internal variability and may also reflect that the ice edge has moved. The change could also reflect the availability of more accurate ice observations. We discuss potential implications for analysis based on coupled climate models, which exhibit large ice edge biases.


Introduction
Many studies have suggested that anomalous Barents-Kara sea ice (BKS) in autumn can trigger predictable shifts in the winter North Atlantic Oscillation (NAO), and hence midlatitude winter weather (Caian et al., 2018;Deser et al., 2007;Dunstone et al., 2016;García-Serrano et al., 2015;Kretschmer et al., 2016;Sun et al., 2015;Wang et al., 2017).This teleconnection manifests as a positive correlation between autumn BKS and the winter NAO, with a reduction in sea ice appearing to force a negative NAO.However, there remains considerable skepticism in the literature on the robustness and even causality of this teleconnection.
A key source of skepticism arises from modeling studies.Recent studies using large ensembles show that coupled climate models largely reproduce a positive BKS-NAO correlation over the satellite era (Blackport & Screen, 2021).However, the magnitude of the correlation is notably smaller in the models compared to estimates based on reanalysis (Blackport & Screen, 2021;Siew et al., 2021;Strommen et al., 2022), and there is considerable ensemble spread, with individual ensemble members simulating a wide range of positive and negative correlations (Blackport & Screen, 2021;Koenigk & Brodeau, 2017;Siew et al., 2021).Several studies have argued that this may suggest the stronger BKS-NAO link seen in reanalysis data is largely a result of internal variability (Koenigk & Brodeau, 2017;Warner et al., 2020), which we here understand to mean a combination of (a) random sampling variability and (b) unforced modes of decadal variability that can be triggered chaotically by noise; such internal variability is thus directly measured by ensemble spread.Idealized modeling studies also show weak signals with large internal variability (Smith et al., 2022;Ye et al., 2023Ye et al., , 2024)).It has additionally been argued that the weak links simulated by coupled models may mostly reflect atmospheric forcing on the ice (Blackport & Screen, 2021).However, it has also been argued that model biases lead to a systematic underestimation of Arctic-midlatitude teleconnections in model ensembles (Smith et al., 2022;Strommen et al., 2022), making it unclear if the large internal variability and weak signals in models are realistic (see also Mori, Kosaka, Watanabe, Nakamura, and Kimoto (2019); Screen and Blackport (2019); Zappa et al. (2021) and Mori, Kosaka, Watanabe, Taguchi, et al. (2019)).
In order to understand the relative roles of internal variability and model biases it is necessary to understand the extent of internal variability in the real world.Because the real world is just a single ensemble member, this cannot be done as straightforwardly as in modeling studies.However, Kolstad and Screen (2019) (hereafter KS19), argue that there is clear evidence of nonstationarity (i.e., variations in time) in the BKS-NAO link in reanalysis data spanning the twentieth century, with the recent period standing out as one of unusually high correlations.If this nonstationarity over the twentieth century is a result of internal variability, then KS19 and the aforementioned modeling studies could be jointly interpreted as suggesting that the BKS-NAO correlation seen in the satellite era does not actually reflect a large, robust, causal relationship.Indeed, KS19 conclude by cautioning against using BKS as a predictor of the NAO.
The purpose of this paper is to make two points concerning nonstationarity of BKS-NAO links, expanding on brief comments made in Strommen et al. (2022).Firstly, while it is well known that observations are sparser further back in time, the implications this has for how confidently one can assess nonstationarity in reanalysis do not appear to have been commented on.Secondly, KS19 suggested that one cause of the apparent nonstationarity could be a dependence of the BKS-NAO link on the mean state, citing decadal North Atlantic variability as a potential source of such mean-state dependence.However, they do not explicitly discuss the potential role of global warming, with their analysis only implicitly suggesting no forced trends.Here, we will argue.
1.That the lack of observations of autumn/winter sea ice means nonstationarity cannot be meaningfully assessed using reanalysis data extending further back than 1950, and is questionable even in the period 1950-1960.2. That changes to the ice edge over time, such as those due to global warming, could have induced some degree of physically forced nonstationarity since 1960.
In particular, while the real world BKS-NAO link necessarily exhibits some degree of internal variability, its magnitude cannot be easily measured, because forced nonstationarity is possible even in the absence of any unforced internal variability.In the Discussion we will also comment on the potential implications of point 2 above for analysis based on coupled models, known to exhibit considerable biases in their simulated ice edge.

Reanalysis and Observational Data
To demonstrate the effect the lack of sea ice observations has on ice-atmosphere links, we focus on ERA20C (Poli et al., 2016), which uses sea ice from HadISSTv2.1 (Titchner & Rayner, 2014).The sea ice in HadISSTv2.1 prior to 1972 comes from the Walsh and Chapman dataset (J.J. Walsh, 1978;Chapman & Walsh, 1991; J. E. Walsh & Chapman, 2001).The Walsh and Chapman dataset is the primary input for other reanalysis data covering the twentieth century, such as version 2 of NCEP's Twentieth Century Reanalysis (Compo et al., 2011) (which uses HadISSTv2.1),JRA-55 (Kobayashi et al., 2015) and ERA5 (Hersbach et al., 2020).Restricting to ERA20C is thus a representative example.More recently, Walsh and Chapman updated their original dataset (J.E. Walsh et al., 2017): this was used in KS19 and is also used in version three of the NCEP reanalysis (Slivinski et al., 2019).This update includes the use of "analog" methods to make a realistic looking guess for ice concentrations in regions/months where no observations are available.While this has the effect of making past sea ice variability look more realistic, the new variability is fundamentally unconstrained by observations and, as shown in KS19, does not produce notably different estimates of ice-atmosphere links than the original dataset.Further discussion on these points is included in Text S1 in Supporting Information S1.We also use ERA5 (Hersbach et al., 2020) for context and when assessing CMIP6 model biases.
To assess the number of available observations in the Barents-Kara region over time prior to the satellite era (approximately 1979 onwards) we consider two sources of observations.Firstly, we use a count of the number of HadISST ship observations of sea surface temperatures (SST).From this we computed the number of available observations in November anywhere within the Barents-Kara region (70-85°N, 30-90°E).This assumes that every ship visiting this region took a measurement of the ice, which is unlikely to be true.This is

Geophysical Research Letters
10.1029/2023GL107609 therefore best thought of as an upper bound on the true number.Secondly, we count the number of available ice edge charts from the Russian Arctic and Antarctic Research Institute (AARI) (A.R. Mahoney et al., 2008).We use the average number of charts available in the Barents-Kara region to measure chart availability.The only other source of Barents-Kara observations used by Walsh and Chapman were charts collected by the Danish Meteorological Institute and the Arctic Climate System Study (J.E. Walsh et al., 2017).However, both these sources only cover the summer months and so do not contribute to estimates of sea ice in October/November.Availability of ship observations and AARI charts therefore provide a reasonable picture of available sea ice observations.

Model Data
To assess how the ice-NAO link may depend on the ice edge mean state, we make use of an ensemble of coupled climate model simulations with stochastic ice and ocean parameterizations.This ensemble was introduced and studied in Strommen et al. (2022), and consists of six members spanning the period 1950-2015 using historical forcing data.The inclusion of stochastic parameterizations results in the model simulating consistently positive BKS-NAO correlations over the period 1980-2015 which are comparable in magnitude to that observed in reanalysis (Strommen et al., 2022).This close and consistent fidelity to observations is not observed in other model ensembles not making use of stochastic ice/ocean parameterizations (Blackport & Screen, 2021;Siew et al., 2021), making it a valuable resource for studying BKS-NAO links.Following Strommen et al. (2022), we refer to this ensemble as OCE.
Details about the model configuration can be found in Strommen et al. (2022).The model is based on a version of EC-Earth3 (Haarsma et al., 2020), itself based on the Integrated Forecast System (IFS) developed by the European Center for Medium Range Weather Forecasts (ECMWF).The ocean component is NEMO version 3.6 (Madec & the NEMO team, 2016) which includes the LIM3 sea ice model (Vancoppenolle et al., 2012).Three stochastic ocean schemes (Juricke et al., 2017(Juricke et al., , 2018) ) and one stochastic sea ice scheme (Juricke et al., 2013(Juricke et al., , 2014;;Juricke & Jung, 2014) are included.The IFS is run at a spectral resolution of T255 with 91 vertical layers.NEMO is run at a resolution of around 1°with 75 vertical layers.

Methods
We follow KS19 and define BKS as sea ice concentration averaged over the box 70-85°N, 30-90°E.We focus on November, rather than October as in KS19.This is motivated by the fact that the physical pathway from October sea ice to the NAO appears to be primarily via its influence on November sea ice (García-Serrano et al., 2015;King et al., 2016;King & García-Serrano, 2016;Santolaria-Otín et al., 2021), with viable atmospheric pathways from November sea ice being more widely documented and studied (García-Serrano et al., 2015;Sun et al., 2015).Furthermore, seasonal forecasts of the winter NAO, such as those issued by ECMWF or the UK Met Office, are initialized using November initial conditions, making November BKS more relevant for actual forecasts.Thus, unless stated otherwise, informal references to ice, sea ice or BKS always refer to November sea ice concentration.The comparison with October will be discussed: the choice of November versus October does not alter our main conclusions.
We define the NAO index in the OCE ensemble as the first principal component of 500 hPa geopotential height; a daily principal component timeseries has a seasonal cycle removed from it before DJF averages are taken.When correlating ice with the NAO in the OCE ensemble, we concatenate all six members back to back before computing the correlation.Both BKS and NAO indices show overall negative trends over 1950-2015 in OCE.
Since the focus here is on the interannual teleconnection, these trends are removed from both indices.See Text S2 in Supporting Information S1 for more information on the detrending procedure.
When determining statistical significance of correlations between ice and the NAO, our null-hypothesis models the DJF NAO as white noise and sea ice using an independent autoregressive model with a lag of 1 year.By fitting these models to the data and generating 1,000 random timeseries, we can estimate p-values for the null hypothesis.
All climate variables are regridded onto a regular 1°grid before analysis is carried out.The heatflux (=sensible + latent) sign convention is that "positive = upwards", that is, heat flowing from the surface to the atmosphere.

Unphysical Causes of Nonstationarity: Missing Data
We first examine the impact of data availability.Figure 1a shows the November BKS timeseries in blue.It is immediately apparent that there is a dramatic difference in variability before and after 1950, with essentially zero variability before 1950.Figure 1d shows the November sea ice variance in the modern period 1980-2010 at all gridpoints in the Arctic, while 1(e) shows the same over the period 1900-1949.It is clear that the variability has vanished almost everywhere, not just in the Barents-Kara region.That this has a big impact on assessments of iceatmosphere interactions can be seen already at the level of the local interaction between two-m temperature (T2M) and sea ice.The black line in Figure 1a shows correlations between BKS-averaged November T2M and November BKS for successive 30-year periods.A sharp discontinuity is apparent, with the correlations jumping from around 0.1 to 0.5 depending on whether the 30-year period includes years post-1950.The correlations drop again to ≈ 0.6 once the period includes years post-1979: note that this drop occurs whether or not trends in ice and T2M from 1990 onwards are removed, suggesting it is related to changes in variability.For completeness, Figure 1c shows 30-year BKS-NAO correlations in ERA20C (black line), reproducing the result from KS19; correlations using ERA5 (red;1950-2022) are shown in red for further context.Figure 1f shows correlations between the NAO and ice at every gridpoint.
To understand these differences, Figure 1b shows the availability of Barents-Kara observations over time.Prior to 1950 there are effectively zero observations in this region in November.Figure 1d suggests this lack of observations extends Arctic-wide and that the Walsh and Chapman dataset consequently use a climatological value for the sea ice.Indeed, the widespread use of climatological values in the Walsh and Chapman dataset pre-1953 is made explicit in the HadISSTv2.1 documentation (Titchner & Rayner, 2014), which contains a more extensive discussion on the limitations of ice data in the earlier twentieth century.In the period 1950-1970, ship observations start becoming available, but there are no AARI charts.A closer examination of the monthly BKS timeseries shows that the variability in the period 1950-1960 exhibits visibly unphysical periodicity, resulting in the November BKS being near-constant across multiple successive years (Figure S1 in Supporting Information S1).
We conclude that estimates of BKS-NAO correlations cannot be sensibly made prior to 1961 due to the collapse of variability owing to lacking observations, which leads to spurious unphysical effects in reanalysis already at the level of local ice-T2M links, particularly pre-1950.Because AARI observations are not widely available until 1970 onwards, we argue that caution should also be shown for the period 1961-1970.We note that 1(c) suggests non-stationarity is present also post-1960, with the correlations mostly increasing from roughly 0.1 to 0.4 (black line, ERA20C) and 0.2 to 0.4 (red line, ERA5).This is discussed further in upcoming sections.
The collapse of past sea ice variability is somewhat less dramatic in October, owing to better availability of AARI observations, but the difference before and after 1960 is still considerable and again leads to apparent nonstationarity in the ice-T2M link (see Figures S1b and S2 in Supporting Information S1).KS19 do state that changes in correlations may be partially influenced by the reduced variability in the early twentieth century: here we show that the extent and source of this reduction places serious limitations for how confidently nonstationarity can be assessed.Most fundamentally, while the NAO timeseries is well constrained in twentieth century reanalysis due to the large number of available observations over the North Atlantic, the BKS timeseries is very poorly (or not at all) constrained prior to 1961, and hence reanalysis/station-based timeseries cannot be expected to reproduce any real-world correlations that may exist between the two.

Physical Causes of Nonstationarity: A Changing Ice Edge
Next, we examine the role of a changing ice edge as a source of nonstationarity.Due to the limitations with reanalysis highlighted in the previous section, we do this using the OCE ensemble.Potential relevance to the real world is discussed in Section 5.
The position of the Arctic ice edge has changed over time, primarily due to global warming, which has led to a gradual retreat of the edge (Notz & Community, 2020;Stroeve & Notz, 2018).This is well reproduced by coupled climate models, including the OCE ensemble.Figure 2a shows the mean state of the OCE ensemble in the period 1950-1980, andFigure 2b the change between this period and the more recent period 1980-2015, demonstrating this ice edge retreat.Changes in the mean ice edge have immediate implications for changes to Arctic variability.This is because the interior of the Arctic is entirely frozen every November (≈100% sea ice concentration) and thus experiences zero interannual variability.Instead all interannual variability is concentrated at the ice edge.This is shown in Figures 2c and 2d, showing November variance in the period 1950-1980 and the difference in the recent period.As the ice edge changes, the regions experiencing considerable variability therefore also change.
Physical reasoning implies that changes to the ice edge variability will impact any teleconnections from the Arctic to the NAO that may exist, because such teleconnections are mediated by heatfluxes.A negative sea ice anomaly may result in comparatively warm Arctic waters being exposed to cold air aloft, and the resulting thermal contrast can trigger heatflux anomalies as high as 500 Wm 2 (Koenigk et al., 2009).These heatflux anomalies generate circulation anomalies that can propagate to the lower latitudes, via tropospheric (Deser et al., 2007) and/or stratospheric (Sun et al., 2015) pathways.Crucially, significant ice-induced heatflux anomalies can only occur at or near the ice edge, since (a) this is the only place where anomalous ice can expose or cover up the ocean, and (b) this is the only place where ice variability occurs at all.This is demonstrated using the OCE ensemble in Figure 2e, which shows the 1950-1980 interannual heatflux variability at gridpoints with a mean sea ice concentration of at least 5%.The heatflux variability is co-located with the 1950-1980 ice edge, and retreats in tandem with the edge under global warming (Figure 2f).
If Arctic sea ice really does exert a non-negligible forcing on the NAO, the above discussion implies that the exact region of the Arctic responsible changes over time.In fact, this is what seems to happen in the OCE ensemble.Figures 2g and 2h show correlations between the winter NAO and November sea ice at every gridpoint for the two time periods.In the earlier period 1950-1980, significant correlations are found in the Barents sea, Greenland sea, and the coast of Greenland more broadly.These correlations are co-located with the peak heatflux variability associated with the more extended ice edge of that period.No correlations are found in the Kara sea (70-85°N, 70-90°E), consistent with the fact that in OCE the Kara sea is almost permanently ice covered in November in the period 1950-1980, and thus experiences little/no ice or heatflux variability (Figures 2a and 2e).In the later period 1980-2015, correlations are still found in the Barents sea, but have largely vanished from around Greenland, consistent with the retreat of the ice and subsequent loss of heatflux variability there.On the other hand, the retreating ice means that the Kara sea has now become partially exposed, and thus begins to contribute non-negligible heatflux anomalies.The OCE ensemble now also shows significant correlations in this region.
The behavior in OCE is thus consistent with physically forced nonstationarity of Arctic-NAO links.What about the BKS-NAO link in particular?Quantitatively, the Kara-NAO correlation of the concatenated OCE ensemble (N = 210) increases from 0.06 (p ≈ 0.35) to 0.15 (p < 0.05).This is a significant change under our null hypothesis (p < 0.05), reflecting the shift from permanently ice-covered to partially exposed.The Barents-NAO correlations are approximately stationary (Figure 2).Thus the naive expectation is that the total forcing from the combined Barents-Kara sea should increase over the period 1950 to present.While this expectation is confounded by several issues which we discuss in Section 5 (including the relative magnitudes of atmospheric forcing from the two regions), it is nominally what seems to happen in the OCE ensemble: the BKS-NAO correlation of the concatenated members is 0.13 (p ≈ 0.05) in the earlier period and 0.24 (p ≪ 0.05) in the modern period.However, Table S1 in Supporting Information S1 shows that the correlation does not increase for all ensemble members individually, and there is considerable ensemble spread.Indeed the increase in the BKS-NAO correlation is not significant under our null hypothesis (p ≈ 0.25).Nevertheless, we note that (a) the BKS-NAO change is consistent with a statistically significant change in the Kara-NAO link having been smeared away by averaging with a stationary Barents-NAO link; (b) the individual BKS-NAO correlations are above and below zero in the early period but uniformly positive in the later period (Table S1 in Supporting Information S1); and (c) the simplicity of the physical mechanism is strong evidence that some change in the BKS-NAO link occurs (Shepherd, 2021).
We stress that our remarks on the relationship between the ice edge and heatfluxes are in no way novel, and many studies have emphasized that the Barents-Kara sea appears to be important by virtue of being where the maximum ice edge variability takes place (Deser et al., 2000;Koenigk et al., 2009;Vinje, 2001).However, the implications this has for nonstationarity of teleconnections do not appear to have been published before.We also note that our argument implies that changes to the ice edge would be expected to induce apparent nonstationarity in the BKS-NAO link even if the link is a result of a common atmospheric driver of both (Blackport & Screen, 2021;Peings, 2019;Siew et al., 2021).Finally, any BKS-NAO correlations must eventually decline and go to zero in a future where Barents-Kara becomes permanently ice free.The potential future weakening of Arctic-midlatitude links has been noted before (Kretschmer et al., 2020;Manzini et al., 2018;Peings & Magnusdottir, 2014).

Discussion
The analysis in Section 4 corroborates our proposed mechanism of forced nonstationarity in the OCE simulations, but we have made no attempt at determining if this mechanism is relevant in the real world, due to several confounding issues.Firstly, different regions of the Arctic may affect the NAO differently (Koenigk et al., 2016;Pedersen et al., 2016;Rinke et al., 2013;Sun et al., 2015), and such links can be sensitive to the background mean state (Deser et al., 2007;Smith et al., 2017;Strong & Magnusdottir, 2010); for example, differences in the climatological position of the jet may result in forcing from the same geographical region affecting the jet differently (Baker et al., 2017(Baker et al., , 2019)).This could lead to differences across models and between models and the real world, which may respond differently to changes in the ice edge.It also means that changes to BKS-NAO correlations may not be independent of other changes in the Arctic, such as those seen around Greenland (see Section 4).Secondly, unforced internal variability is clearly non-zero in the real world, but it is hard to reliably estimate the magnitude of this variability given the sample size limitations highlighted in Section 3. Thirdly, with reanalysis data there is an added complication that higher correlations in the satellite period compared to presatellite could be influenced by the availability of better observations.Figure 1c shows that in the more reliable period 1961-2022, correlations in reanalysis rise rapidly from ≈0.1 to ≈ 0.4; the ice edge variability changes over the same period (see Figure S6 in Supporting Information S1), much like seen in OCE (Figure 2).This is consistent with our proposed forced nonstationarity, but also with a hypothesis that more/better observations have resulted in the teleconnection being measured more accurately (see Text S1 in Supporting Information S1).It is plausible that all three factors (forced nonstationarity; internal variability; improved observations) play a role in modulating Arctic-NAO links over time in models/reanalysis, making it challenging to assess the role of any individual factor.We made no attempt at addressing this challenge here, including for the OCE ensemble.However, if the overall BKS-NAO link is weak, as several studies have suggested (Blackport & Screen, 2021; Geophysical Research Letters 10.1029/2023GL107609 Smith et al., 2022), then the forced nonstationarity would also be weak.It is also possible that the change in the real-world autumn ice edge has not yet been large enough to generate appreciable forced nonstationarity.
The fact that Arctic teleconnections may be linked to the location and variability of the ice edge is not just relevant for nonstationarity, but also has potential implications for model studies.It is well known that models exhibit considerable biases in both the mean and variability of Arctic sea ice (Gastineau et al., 2020;Khosravi et al., 2022;Koenigk et al., 2014;Notz & Community, 2020;Roach et al., 2018;Watts et al., 2021).It follows that the Arctic regions capable of forcing the NAO may vary from model to model.Most studies using models apply a predefined BKS region to both reanalysis and models alike (Blackport & Screen, 2021;Kolstad & Screen, 2019;Siew et al., 2021).While this avoids potential "cherrypicking", it also risks exaggerating the weakness of model signals.For example, given a model with a strong forcing from the Barents sea but no forcing from the Kara sea (e.g., due to the model simulating a permanently ice-covered Kara sea), a correlation based on the Barents-Kara sea could be misleading.Figure 3 demonstrates that CMIP6 models cannot be assumed to exhibit non-trivial sea ice variability in either the Barents or Kara seas.It seems of clear interest to assess how ice edge biases may affect model teleconnections, and furthermore to develop methods that allow for a more objective, physically motivated way to identify which Arctic regions may be forcing the NAO in a given model.
The behavior of the OCE model corroborated our physical reasoning for nonstationarity, but this might be a particular feature of OCE which other models do not replicate.While the OCE ensemble appears to be uniquely good at replicating the observed BKS-NAO correlations, the reasons for this are still poorly understood, and the ensemble does exhibit some biases, such as high sea ice concentration in the 50s (Figure S8 in Supporting Information S1). Figure S6 in Supporting Information S1 also shows that the spatial pattern of correlations is different between OCE and ERA5 in the pre-satellite period.The behavior of OCE may therefore not reflect the behavior of the real world.

Conclusions
We have shown that a lack of past observations from the Barents and Kara seas means that nonstationarity cannot be sensibly assessed using renanalysis prior to 1950, and great caution should be shown all the way up to at least 1961.We have also argued that simple physical reasoning implies forced nonstationarity is possible due to globalwarming induced changes in the ice edge, since the location and variability of the ice edge determines the location of the heatflux anomalies responsible for generating Arctic signals.In particular, there are three distinct "regimes" that an Arctic region can inhabit: (a) permanently ice-covered; (b) partially covered; and (c) permanently ice-free.
Teleconnections to the midlatitudes can only exist for a region in the second regime, and hence ice edge changes will generate nonstationarity by forcing different regions to transition between the three regimes.Our model ensemble does simulate changes to Arctic-NAO correlations which are consistent with such forced nonstationarity between 1950 and 2015, but the individual ensemble members were too variable to distinguish a statistically significant change in the BKS-NAO correlation over this period.We also emphasize that we have not attempted to estimate the magnitude of this effect in the real world, and it could be negligible relative to internal variability.
Importantly, while this does suggest that BKS-NAO correlations may have been lower in the past, it would be wrong to conclude from this that there is no robust and causal forcing on the NAO from BKS in recent decades.Rather, it may simply reflect that ice-edge changes mean some regions, like the Greenland sea, were more important in the past, and some were less important, like the Kara sea.These arguments were corroborated using climate model simulations which exhibit strong correlations from the Greenland-Barents region to the NAO in the period 1950-1980, which move along with the ice edge to the Barents-Kara region in the period 1980-2015.In other words, Arctic sea ice may have been causally forcing the NAO across the entire twentieth century, just not always from the same place.When restricted to the more reliable period 1961 onwards, ERA20C and ERA5 both exhibit steadily increasing BKS-NAO correlations, consistent with a forced response due to the increased exposure of the Kara sea and decreased Greenland sea ice.Improved availability of observations post-1979 may also contribute to these increased correlations.Thus changes in BKS-NAO links over time cannot be conclusively attributed to internal variability alone.
We conclude that the magnitude of internal variability of the BKS-NAO link in the real world remains unclear.
The possibility that models exaggerate the magnitude of unforced internal variability and underestimate the magnitude of forced signals from BKS cannot therefore be easily ruled out.Models are known to underestimate forced signals in the NAO on seasonal-to-decadal timescales (Scaife & Smith, 2018), lending further weight to this possibility.We also hypothesize that model biases in the ice edge may have obfuscated a clear understanding of Arctic-NAO teleconnections in multi-model studies which prescribe the Barents-Kara region upfront.This is a clear avenue of future research.

Figure 1 .
Figure 1.In (a): November BKS timeseries of ERA20C (blue) and successive 30-year correlations between November BKS and November T2M averaged over the same region (black).Each 30-year correlation is centered at the midpoint of the period (i.e., the point at 1965 corresponds to the period 1950-1980).In (b): total number of ship SST observations (red) and average number of AARI charts (black) over the Barents-Kara region.In (c): 30-year correlations between November BKS and the NAO for ERA20C (black) and ERA5 (red; 1950-2022).The vertical mauve line indicates the point at which the 30-year chunks contain no years pre-1961.In (d): November interannual sea ice variance of ERA20C.In (e), the same but over the period1900 -1949.In (f).In (f): correlations between the NAO and November sea ice at every gridpoint.The blue boxes highlight the Barents-Kara region.

Figure 2 .
Figure 2. In (a) the mean November sea ice across the OCE ensemble in the period 1950-1980, and (b) the difference in the ice mean between1980 -2015  and 1950 -1980.In (c) .In (c)  and (d): the same but for the variance rather than the mean.In (e) and (f) the same but for the heat flux variance.In (g): correlations between the DJF NAO and November sea ice concentration at gridpoints in the period 1950-1980 using the OCE ensemble.In (h): the same but over the period 1980-2015.In (g) and (h) all gridpoints outside the zero contour are significantly different from our null-hypothesis ( p < 0.05).The blue boxes highlight the Barents-Kara region.The heatflux sign convention is "positive = upwards".

Figure 3 .
Figure 3. Interannual November sea ice concentration variance over the period 1980-2015 for (a) the ERA5 reanalysis, (b)-(d) three different coupled CMIP6 models (historical forcing scenario).The models have been hand-selected for being illustrative.The Barents-Kara region is highlighted with blue boxes.