Observed winter Barents Kara Sea ice variations induce prominent sub-decadal variability and a multi-decadal trend in the Warm Arctic Cold Eurasia pattern

The observed winter Barents-Kara Sea (BKS) sea ice concentration (SIC) has shown a close association with the second empirical orthogonal function (EOF) mode of Eurasian winter surface air temperature (SAT) variability, known as Warm Arctic Cold Eurasia (WACE) pattern. However, the potential role of BKS SIC on this WACE pattern of variability and on its long-term trend remains elusive. Here, we show that from 1979 to 2022, the winter BKS SIC and WACE association is most prominent and statistically significant for the variability at the sub-decadal time scale for 5–6 years. We also show the critical role of the multi-decadal trend in the principal component of the WACE mode of variability for explaining the overall Eurasian winter temperature trend over the same period. Furthermore, a large multi-model ensemble of atmosphere-only experiments from 1979 to 2014, with and without the observed Arctic SIC forcing, suggests that the BKS SIC variations induce this observed sub-decadal variability and the multi-decadal trend in the WACE. Additionally, we analyse the model simulated first or the leading EOF mode of Eurasian winter SAT variability, which in observations, closely relates to the Arctic Oscillation (AO). We find a weaker association of this mode to AO and a statistically significant positive trend in our ensemble simulation, opposite to that found in observation. This contrasting nature reflects excessive hemispheric warming in the models, partly contributed by the modelled Arctic Sea ice loss.


Introduction
The second empirical orthogonal function (EOF) mode of Eurasian winter (December-February) surface air temperature (SAT) variability, known as Warm Arctic-Cold Eurasia (WACE) pattern, is closely related to Ural blocking or Siberian high (Mori et al 2014, 2019, Tyrlis et al 2019, Luo et al 2021).This WACE pattern of variability closely associates with the variations in Barents-Kara Sea (BKS) sea ice concentration (SIC) during the satellite era (Overland et al 2011, 2013, Cohen et al 2014, Horton et al 2015, Kug et al 2015).These associations have shown an upward trend after 2000 with a more persistent blocking regime (see figure 1(c)) with higher intensity (Li et al 2023).
Observational studies have suggested a possible role of the BKS SIC in the WACE pattern of variability and further in the cold Eurasian winter conditions (Petoukhov and Semenov 2010, Cohen et al 2013, Tang et al 2013, Mori et al 2014, Overland et al 2015, Kretschmer et al 2016, Kim and Son 2020, Xie et al 2020, Rudeva and Simmonds 2021, Wu et al 2022, Zhong and Wu 2022).However, the WACE pattern-related Ural blocking itself is shown to influence the BKS SIC variations on sub-seasonal to interannual time scale (Gong and Luo 2017, Blackport et al 2019, Peings 2019, Tyrlis et al 2019, Screen and Blackport 2019b, Tyrlis et al 2020, Blackport and Screen 2021, Komatsu et al 2022).Even a colder Eurasian temperature can bring Warm Arctic conditions by enhancing the Siberian high (Wu and Ding 2023).The frequent occurrence of the WACE pattern is also shown to be favoured by the warm phase of the Atlantic Multidecadal Variations (Luo et al 2017a, Jin et al 2020) and a negative phase of Pacific Decadal Oscillation (Luo et al 2022).Hence, it remains to be seen if any portion of the observed WACE mode of variations is essentially forced or modulated by the observed Arctic or BKS SIC variations.
Regarding the multi-decadal trend, the observational studies suggest that the BKS SIC loss-related warming in the Arctic reduces the tropospheric potential vorticity gradient over Eurasia, making the WACE related Ural blocking more persistent (Yao et al 2017, Luo et al 2017bLuo et al , 2018Luo et al , 2019b)).A dynamic and thermodynamic coupled view has also been proposed to explain the role of the warming Arctic in bringing cold Eurasian condition (Xie et al 2022).However, results from the previous climate model experiments have yet to establish a clear role of the recent Arctic sea ice loss on the observed Eurasian cooling.They are indicated to be mainly due to atmospheric internal variability (McCusker et al 2016, Sun et al 2016, Ogawa et al 2018, Wang and Chen 2022).Previous theory suggests that the nature of an internal mode of variability or circulation regime in the extratropics can be influenced by the external forcing (Palmer 1999).Hence, a comprehensive multimodel experimental causal analysis of the model's ability to precisely simulate the observed trend in the WACE mode of variability and the potential role of observed Arctic Sea ice loss behind such a trend is still to be explored and essential.Because through this trend in the observed WACE mode of variability, Arctic sea ice loss is suggested to relate to Eurasian cooling, apart from the contributions of the internal mode of variability (Mori et al 2014).
Apart from the WACE mode of variability, the first or the primary mode of observed Eurasian SAT variability is found to be closely related to Arctic Oscillation (AO), expressed as the first EOF mode of Northern Hemisphere (NH) (20 (Mori et al 2014).However, the simulated nature of this mode of variability and its association with AO remains unexplored, though important given its vital role in determining the long-term trend of observed Eurasian winter temperature (figure 2).Hence, to investigate the above-mentioned open questions, we analyse two sets of Atmospheric Model Intercomparison Project (AMIP)-type large ensemble of multi-model simulations with (ALL experiment) and without (SICclim experiment) observed Arctic SIC variations for the period 1979-2014 (see materials and methods) and compare with ERA5 reanalysis as a proxy for observations.

Data and methods
The observed winter season (December-January-February) monthly mean data for SAT, SLP and 500 hPa Geopotential height from 1980 (meaning December 1979 andJanuary-February 1980) to 2022 is taken from ERA5 reanalysis (Hersbach et al 2017(Hersbach et al , 2020)).SAT in ERA5 and the models denotes the air temperature at 2 m or reference height, depending on the available model output.The observed SIC data for the same period  is taken from the U.K. Met Office Hadley Centre SIC and SST version 1 (HadISST1) (Rayner et al 2003).Following Tyrlis et al 2020, the Ural Blocking Frequency (UBF) is calculated from the ERA5 daily data of potential temperature (θ) on 2PVU (potential vorticity unit, 1.0 × 10 −6 m 2 s −1 K kg −1 ) surface, which indicates the dynamical tropopause in the extra-tropics (Hoskins et al 1985).The blocking algorithm identifies wave-breaking areas, reversing the climatological meridional gradient of θ at the 2PVU surface (Tyrlis et al 2015(Tyrlis et al , 2019)).The UBF is defined as the number of Blocking episodes over 40   S1).For each of the eight models considered, various numbers of ensemble members ranging from 10 to 30 are produced, leading to a total of 145 members for ALL and SICclim.Following the previous research based on the same experiments (Liang et al 2020, 2021, Suo et al 2022), we treat the multi-model ensemble as a single-model large ensemble by giving the same weight to each member of each model.
For ERA5, the EOF1 and EOF2 of the SAT over Eurasia (20 • -90 • N, 0 • -180 • E) and SIC over the Northern Polar circle (60 • -90 • N, 0 • -360 • E) are calculated using the anomaly covariance matrix (North et al 1982).The EOF patterns remain consistent with rotation.Following (North et al 1982), a test of the sampling error in the eigen values is performed to confirm the absence of degeneracy between the Eurasian SAT EOF2 and EOF3 patterns.To represent the EOF patterns in the units of SAT (in Kelvin) or SIC (in %), the normalised EOF patterns are multiplied by the square root of their corresponding eigenvalues.The principal component (PC) time series are represented in normalised magnitude.In ALL and SICclim, each model's EOF1 and EOF2 SAT patterns are determined by performing the same analysis after concatenating the ensemble members for each model separately.The corresponding PC time series is then unpacked back to represent the PC time series for each ensemble member.Finally, the average appearance of EOF patterns in the models is shown by performing the multi-model means (MMM) of the EOF patterns.In ERA5-2014/2022, the portion of the entire data that is associated with a particular EOF pattern is constructed by multiplying the normalised EOF pattern with its PC time series.
To detect periodicity peaks, a power spectrum analysis is performed on the observed WACE mode, BKS SIC, SH and UBF Index time series.No smoothing is applied to the time series that are detrended and tapered in the time domain before calculating the spectrum.In addition, a theoretical Markov spectrum and its 95% (upper) spectrum are constructed using the lag-1 autocorrelation in the time series to identify the peaks in the spectrum that are statistically significant.The power spectrum is built for each ensemble member for the model experiments.Then an ensemble mean of those spectrums is constructed to show the persistent peaks in periodicities in the WACE mode time series.
The Theil-Sen nonparametric trend estimation method determines any linear trend of SAT fields in ERA5 (Sen 1968).The results are consistent and do not change using the least square method.The Mann-Kendall (M-K) non-parametric test (Mann 1945, Gilbert 1987) is performed on the Theil-Sen linear trend estimate to determine if the trend differs significantly from zero at the 5% level.For the linear trends in PCs, a value above 0.04 yr −1 turns significant at p <0.05 according to the M-K test.For the ensemble mean field, the statistical significance is determined using a one-sample, two-sided Student's t-test, comparing the ensemble mean value with zero at a 5% significance level.

The observed trend and sub-decadal variability of the BKS SIC and the WACE pattern
The first mode of observed Arctic SIC variations over the winter months (December to February) during the satellite period  shows a dominant variability over the BKS region (figures 1(a) and (c)).This region of the Arctic also shows one of the highest long-term trends in sea ice decline over this period (Simmonds andLi 2021, Thoman et al 2022).The observed WACE or the second EOF pattern of the Eurasian winter SAT variations also show prominent variations and a persistent upward trend, both closely associated with the BKS SIC variations (correlation = 0.82, p-value = 0.013) (figures 1(b) and (d)).This observed WACE mode of variability is also related to the observed SH index (correlation = 0.87, p-value = 0.002), confirming the close association of this mode with the key observed atmospheric circulation feature over Eurasia (figure 1(b)).
Here, we find that the observed BKS SIC exhibit prominent sub-decadal variations with a statistically significant spectral peak at around 5.5-6 years' time scale (figure 1(d)).A recent study reported a similar finding from an even longer record of the SIC data (Luo et al 2023).We also find the same statistically significant spectral peak in the observed WACE time series, SH index and UBF index (figure 1(d)).This affirms the WACE time series representing the most prominent variations in the Ural Blocking or SH index.Because, for all time series, the variations in this timescale hold most of the variance in the power spectrum analysis.Hence, now the question arises whether such sub-decadal variations in the WACE pattern could occur without the BKS SIC variations or if they essentially require feedback of the BKS SIC variations onto the WACE.

The role of WACE in shaping the observed Eurasian winter temperature trend
Apart from the prominent sub-decadal scale variability, the WACE time series has a statistically significant positive trend (figure 1(c)).The observed Eurasian cooling is suggested to be contributed by this trend in the WACE mode in association with the BKS SIC loss (Mori et al 2014).The Eurasian cooling trend has weakened in the recent year after 2014 (figure S2).This leads to a perception that the Arctic-Eurasia link during winter has waned in recent years (Blackport and Screen 2020).However, the trend in the WACE mode and BKS SIC association persists (figure 1(c)).Moreover, the portion of SAT trend contributed by the WACE-related SAT anomalies (for details, see data and methods) shows a stronger contribution to the central Eurasian temperature trend for 1980-2022 compared to 1980-2014 (figures 2(c) and (d)).This apparent paradox could be resolved if we examine the changes in PC1 of Eurasian SAT variability.
Our findings suggest that the recent reduction in the Eurasian cooling is a consequence of the change of phase in PC1 from negative to positive after 2014, by the change of phase in AO (compare black and blue curves in figure S1(c)).The PC1 trends in observations are not statistically significant either until 2014 Nevertheless, it should be noted that despite the reversal in the AO-related variability-induced trend from cooling to warming after 2014, the Eurasian winter temperature trend does not reverse (figures S2(a) and (b)).This shows a crucial role of the trend in the WACE mode of variability (figures 1(c) and 2(d)).Without this trend, we would not have been able to explain the observed central Eurasian temperature trend, especially in recent years under a positive phase of AO.If we remove the AO-related PC1 SAT trend from the entire SAT field, we find an enhanced residual Eurasian cooling trend in the SAT data over the last years (figures 2(e) and (f)).Moreover, the cooling in the residual data is statistically significant for a much larger area compared to the full data due to the removal of the AO-related SAT variations from the data (compare figures 2(e) and S2(a)).This residual Eurasian cooling can only be explained by the trend induced by the WACE (figures 2(c) and (d)).Hence, understanding the cause of the trend in the PC of the WACE mode of variability remains vital.

Simulated modes of winter Eurasian temperature variability
In both the ALL and the SICclim experiments, we calculate the EOF of the SAT over Eurasia (0 • -180 • E and 20 • N-90 • N) using all members concatenated separately for each model.The MMM EOF1 and EOF2 patterns, obtained by averaging the respective EOF patterns across the models, resemble those of ERA5 (pattern correlation ∼0.9 for each mode in MMM with a range of 0.8-0.9among the different models; figures 3(a), (e) and (b), (f) vs figures S3(a) and (b)).The associated atmospheric circulation, derived from the regression of the SLP fields on the respective normalised PCs, also resembles the observed structure.The monopolar central Eurasian (∼60 • N, 90 • E) warming pattern of EOF1 is associated with a low SLP centred around BKS (figure 3 The WACE pattern is primarily driven by the Ural blocking/Siberian high-related variability (figure 1(c)).It exists in SICclim (figures 3(f) and S3) but with a different characteristic than ALL.The extensive similarity in the WACE patterns for ALL and SICclim confirms that the WACE pattern is mainly driven internally to the atmosphere, ensuring previous studies (Mori et al 2014, Sorokina et al 2016, Mori et al 2019, Peings 2019).Nevertheless, the WACErelated positive SLP anomaly is systematically intensified in MMM and shifted northward towards the BKS in ALL compared to SICclim (figures 3(b) vs (f)).This enhancement of the WACE-related high-pressure anomaly under Arctic SIC loss in climate models is consistent with previous studies (Screen and Blackport 2019a).Moreover, the explained variances of the EOF2/WACE are slightly lesser in SICclim (range 12%-15%) than in ALL (range 12%-18%).
In MMM, the WACE SAT anomaly for ALL (figure 3(b)) and its positive centre over the BKS are similar to the observations (figure S3(b)), whereas the warm node is much weaker in SICclim (figure 3(f)).This suggests a link between the WACE and the BKS SIC variability-related temperature changes over the Arctic, which is missing in SICclim.The cold node of the WACE solely depends on the strength of the atmospheric circulation anomalies associated with WACE.For individual models, we see a consistent link between the strength of the simulated high-pressure anomaly and the magnitude of the central Eurasian cold anomaly in ALL and SICclim (figure S4).
In the case of the change in the circulation anomalies from ALL to SICclim associated with the changes in the EOF2 SAT anomalies, MMM depicts a colder northern Eurasia under the Arctic SIC forcing with an associated enhanced high-pressure anomaly (figure 4(a)).In the case of individual models, we continue to find consistent response patterns, with the models showing a stronger high-pressure anomaly bringing a larger cooling response over Eurasia (figures 4(b)-(i)).
The core of the cooling response differs among the models.Some models project the cooling in the north-western part of Eurasia, such as CMCC-CM2-HR4, EC-Earth3, and HadGEM3 (figures 4(a), (c) and (h)).Whereas there are also models with the core of the cooling response more towards the central to the eastern side of Eurasia, e.g.ECHAM6.3,CAM6-Nor, CESM2-WACCM6 (figures 4(d), (f), (g)).This difference is closely connected with the difference in the associated circulation structure, which confines more towards the west in the north-western Eurasian cooling-centric models.However, overall, a cooling associated with high pressure of varied strength can be seen in the model WACEs under Arctic SIC forcing.
Despite capturing the observed feature of the WACE pattern and the associated high-pressure enhancement, the variance explained by SAT EOF2 or WACE is systematically underestimated in the model experiments.The MMM of explained variances is 16% (with a range of 12%-18%) in ALL compared to 20% in ERA5.This is in line with the characteristics of the CMIP6 climate models, which show a systematic underestimation of the blockings over the Urals (Davini and D'Andrea 2020).

Low-frequency sub-decadal WACE variability under BKS SIC forcing
Apart from the WACE pattern, the multi-model ensemble mean PC2/WACE time series in ALL captures the response in the WACE time series forced by SST, SIC, and external forcing (as opposed to that driven internally by the atmosphere) and is highly similar to the sign reversed averaged observed BKS SIC time series with a correlation 0.93 (pvalue = 0.0008; red and blue curves, figure 3(d)).This relation remains important when using detrended time series (correlation 0.86, p-value < 0.0001).This close association of the BKS SIC variations and WACE mode is also found in observations (figure 1(c)).Further, in line with observations, the ensemble mean WACE time series in ALL has a prominent lowfrequency variation with a long-term positive trend.In contrast, such variability and trend are absent in SICclim (compare red curves, figures 3(d) vs (h)).This suggests a role of observed Arctic sea ice variations behind a prominent low-frequency variability and trend in the observed WACE time series.We perform power spectral analyses to affirm our inference regarding the low-frequency variations further.
The ensemble mean of power spectral density from all the WACE time series (see materials and methods) in ALL reveals a statistically significant spectral peak at around 5.5-7.5 years, consistent with the peak in observed BKS SIC variations (figures 5(a) and S5).Although, the corresponding spectral peak is missing in SICclim, which is not forced by the   observed Arctic/BKS SIC variations (figure 5(b)).We also find a statistically significant spectral peak at around four years in WACE from ALL and SICclim.This could result from other forcing factors in addition to Arctic SIC, e.g.external forcings or SSTs over the ENSO region or North Atlantic Ocean (Luo et al 2019a(Luo et al , 2023)).Overall, the main and only statistically significant difference in the power spectra between ALL and SICclim is the ∼6 year peak in ALL, suggesting that such periodicity in the WACE could not exist in the absence of the observed BKS SIC variations, which also shows a significant spectral peak at ∼6 year (figure 1(d)).
Apart from variability, the changes in the longterm trend of WACE and PC1 SAT under observed Arctic sea ice variations can be investigated by comparing their joint probability density functions (JPDFs) in ALL and SICclim with those in ERA5.The WACE or PC2 trends JPDF in ALL exhibits the ensemble mean of a statistically significant positive value of 0.040 yr −1 under observed SIC forcing, while the ensemble mean of the PC2 trends in SICclim is nearly zero (figure 6(a)), i.e. no statistically significant trend when driven by the climatological SIC and observed SST forcing.We find a statistically significant WACE trend of 0.046 yr −1 in observations, crucial in explaining the observed central Eurasian winter temperature trend (see figure 2 and its explanation in section 3.2).Such WACE trend could only be found in ALL and not in SICclim (figure 6(a)).Hence the experiments suggest that this important observed trend in WACE mode could only be found under observed Arctic/BKS SIC forcing.

Differences in observed and simulated first mode of Eurasian temperature variability
The modelled PC1 shows a striking difference from the ERA5 PC1 time series (figures 3(c) and (g)).The PC1 time series in ALL and, to a lesser extent, in SICclim have a consistently positive trend (grey lines), which we can also in their ensemble mean (red lines), while ERA5 shows no statistically significant PC1 trend (black line in figure 3(c)).The difference in the observed and simulated trends could be visualised in the JPDFs (figure 6(a)).The JPDF of the SICclim PC1 trends (blue contours, figure 6(a)) shows a small but statistically significant positive ensemble mean of ∼0.020 yr −1 (blue dot, figure 5(a)).The mean PC1 trend is enhanced in ALL under observed Arctic SIC forcing (shaded contours, figure 6(a)) with a statistically significant ensemble mean trend of 0.043 yr −1 (red dot, figure 6(a)).The observed PC1 trend from ERA5 is negative and not statistically significant at 5% (−0.015 yr −1 ).Not a single ensemble member displays such a negative PC1 trend.The observed PC1 trend being out of the spread of the model simulated PC1 trends suggest the differing nature of the trend in model PC1s from the observation.
The essence of this difference in the trend can be captured in the relationship between the PC1 and the AO.In ERA5, due to the close relation between PC1 and AO, the observed PC1 reflects the trend in the observed AO, which is not statistically significant (black star in figure 6(b)).Accordingly, in experiment ALL, the ensemble mean of the AO trend is close to zero due to an almost equal spread of the AO trends in both positive and negative sides arising mainly from AO-related internal variability.Further, the PC1 and AO trends also have a moderately close association in ALL (correlation = 0.5, p-value < 0.0001).However, unlike AO trends, all PC1s show positive trends and a statistically significant positive ensemble mean trend (red dot, figure 6(b)).The distribution of the PC1-AO correlations in ALL shows that most members show a relation weaker than that observed (figure S6).The distribution of the PC1-AO correlations shifts to a slightly higher value in SICclim compared to ALL, though the changes are not statistically significantly different.

Conclusions
We have shown that the Siberian high/Ural blocking related WACE pattern of variability exhibits a close association with BKS SIC variations in terms of the substantial ∼6 years variability and longterm trends both in the observation and the large ensemble atmospheric model simulations forced with the observed Arctic SIC variations.The statistically significant ∼6 years cycle and the long-term trends are no longer found when the models are forced with climatological mean Arctic SIC.These results demonstrate that although the WACE is an atmospheric internal mode of variability affecting BKS SIC variations, it has an influence or feedback from BKS SIC variations.We have also shown that the observed long-term trend in the WACE pattern continues to persist until winter 2022 and is essential in explaining the Eurasian SAT trend.
Our results mainly focus on the driving role of BKS SIC variations, as we cannot investigate the origin of the BKS SIC variability with our simulations.Previous works established that the atmosphere plays a major role in driving the observed BKS SIC variations in shorter time scales of weeks to months (Sorokina et al 2016, Peings 2019, Tyrlis et al 2020), together with oceanic processes on longer time scales (Jungclaus et al 2005, Sato et al 2014, Jung et al 2017, Yamagami et al 2022, Luo et al 2023).The observed BKS SIC variations imposed in ALL bear the signature of those influences.Hence the factors which influenced the BKS SIC variability and their association with Eurasian winter temperature through BKS SIC hold.Our results only show that BKS SIC variations have direct feedback on the WACE pattern, especially on sub-decadal timescales.And further, without BKS SIC variations and loss (in SICclim), the ∼6 years cycle and the long-term changes in the WACE pattern of variability are not found.Hence, the BKS SIC lowfrequency variations and trend are central in bringing the observed low-frequency association and the long-term trend to the WACE.Furthermore, we show a pivotal role of the multi-decadal trend in the WACE mode of variability for shaping the overall Eurasian winter SAT trend (see figure 2).Hence our results underscore a critical role of BKS SIC loss for shaping the observed Eurasian winter SAT trend.
In the simulations, we have shown that the root of the inconsistency between observations and the model for Eurasian winter is confined to the leading mode of Eurasian SAT variability (PC1), which is related to the AO.The PC1 in the model simulations show a statistically significant positive trend where the observed PC1 trend is found to be outside the spread of the model simulated trend.Additional research is necessary to understand better this disparity between the observed and model-simulated PC1 trend and the possible influence of Arctic sea ice loss on this divergence.Our results suggest that the observation and model differences in this aspect could be partly related to a too strong thermodynamics response reflecting in an overall higher land surface and hemispheric warming in the models and/or underestimation of the trend or the low-frequency component in AO-related dynamics in the models (Scaife andSmith 2018, Smith et al 2020).These inconsistencies could stem from the missing sea ice physics (Marcq andWeiss 2012, Lang et al 2017), missing Ocean-Sea ice-land-atmosphere coupled interaction (Smith et al 2017, Screen et al 2018), potential biases in sea ice-atmosphere interactions, atmospheric physics, Ocean-atmosphere teleconnections, two-way stratosphere-troposphere coupling.
Having the monopolar warm anomalies over almost the entire EOF1 domain, the modelled PC1 time series could capture the forced thermodynamic warming signals.A trend of the averaged NH SAT would represent such a warming level.We find that the modelled NH warming is much higher than the observed, which lies outside the spread of the modelsimulated trends in ALL (figure S7 The close relationship between the observed PC1 and AO explains why, during the first half of the 2010's, with the AO negative phase, the PC1 trend reinforced the WACE cold state over Eurasia.Thereafter (after 2014), the recently observed reduction in the Eurasian cooling trend is due to the shift to the AO positive phase, which has an opposing influence compared to the persistent WACE trend.
Finally, our finding of the link between the BKS SIC variations for the ∼6 years cycle in the observed WACE pattern of variability provides a potential source of sub-decadal predictability for Eurasian winter temperature.Further research should be directed to understand the possible drivers of such variations in the BKS SIC, the representation of those driving mechanisms in the coupled climate models and their changes under global warming.

Figure 1 .
Figure 1.For the period 1980-2022, (a) EOF1 of observed winter (December-February) Arctic SIC variability (60 • N-90 • N, 0 • -360 • E) in percentage from HadISST1, (b) from ERA5 reanalysis, the EOF2 of the observed winter Eurasian (20 • N-90 • N, 0-180 • E) SAT variability in Kelvin (in shading), known as Warm Arctic Cold Eurasia pattern.Black (grey) contours show the associated sea level pressure (500 hPa geopotential height) anomaly.The top right of the figures mentions the explained variance.(c) Associated time series of the PC1 of Arctic SIC variability in (a) (in light blue), the sign reversed area averaged winter Barents-Kara Sea (BKS) SIC anomaly (in blue, 65 • N-85 • N, 20 • E-90 • E, blue bold contour area in (a)), PC2 of the Eurasian SAT variability in (b) (in red), area-averaged winter Siberian high (SH) index (in black, 60 • N-85 • N, 10 • E-110 • E, black bold contour area in (b) and winter Ural Blocking Frequency (UBF) based on (Tyrlis et al 2020) averaged over 40 • E-100 • E at 61 • N.Each time series are normalised.(d) In the same colours, the power spectrum of the same five-time series is shown in (c) with the 95% (upper) confidence bounds of the associated Markov spectrum (in the same-coloured dashed lines).

Figure 2 .
Figure 2. From ERA5 reanalysis, the DJF SAT trend for the period (a) 1980-2014 and (b) 1980-2022 of the portions of the SAT data, which is associated with EOF1/PC1 of the Eurasian (0 • -90 • N, 0 • -180 • E) SAT.(c)-(d) as in (a)-(b) but for the portion of the SAT data associated with the EOF2/PC2 or WACE pattern of Eurasian SAT.The residual trend in the full SAT data over the box region is shown in figures (c) and (d) for the period (e) 1980-2014 and (f) 1980-2022 after the portion of the trend induced by PC1 SAT is removed.Stippled areas in (e) and (f) are showing regions where the residual SAT trend is statistically significant at 5% level.All units are in K/decade.
(a)).The variation of this SLP low is explained mainly by the Northern lobe of the AO(Mori et al 2014).The EOF2 or WACE pattern shows its warm centre over the BKS and cold centre over central-eastern Eurasia (40• -60 • N, figures 3(b), (f), 1(b) and S3(b)).The WACE pattern is associated with a high SLP centred around northern Eurasia/Siberia and related to Siberian high/Ural blocking(Mori et al 2014, Luo et al 2016, Gong and Luo 2017, Tyrlis et al 2019).

Figure 3 .
Figure 3.The multi-model mean (a) EOF1 and (b) EOF2 patterns in the winter (DJF) SAT (in K) over Eurasia (20 • N-90 • N, 0 • -180 • E) in ALL with daily year-to-year varying observed SST and SIC conditions for the period of 1980-2014.The black contours are the multi-model mean SLP (in hPa) fields associated with the EOFs, derived by regression of the SLP on the respective normalised PC time series for each model.(c) The normalised PC1 time series in grey for the 145 ensemble members of the ALL.The black line is the same but from ERA5 data for 1980-2014.The red line is the ensemble mean of the PC1 time series, and the brown lines are the 95% confidence intervals of the ensemble mean.(d) As in (c) but for the PC2, i.e.WACE mode.The other blue curve is the normalised sign reversed time series of the observed area averaged winter SIC anomalies in BKS from HadSST1 (65 • N-85 • N, 20 • E-90 • E, blue box in figure 1(a)).(e)-(h) are the same as (a)-(d) but for the SICclim experiment with daily climatological Arctic SIC forcing instead of year-to-year varying SIC.Top right corners of (a), (b), (e), (f) mention the multi-model mean of the percentages of the total variance explained by the respective EOFs.

Figure 4 .
Figure 4.The difference between ALL and SICclim for the EOF2 SAT pattern, i.e.WACE pattern (in shading), and its associated SLP pattern (in contours) for (a) multi-model mean (MMM) and in individual models, i.e.(b) CMCC-CM2-HR4 (c) EC-Earth3 (d) ECHAM6.3(e) LMDZOR6 (f) CAM6-Nor (g) CESM2-WACCM6 (h) HadGEM3 and (i) IAP4.Units are in Kelvin (K) and hPa.Stippled areas in (a) are the regions where the MMM differences are significantly different from zero at 95% confidence level relative to the spread of the difference in the 8 models.

Figure 5 .
Figure 5. Ensemble mean of the power spectrums from WACE/PC2 time series of 145 ensemble members in (a) ALL (Red) and in (b) SICclim (Blue), and the dashed lines in the same colours are the 95% (upper) confidence bound of their respective Markov spectrum.

Figure 6 .
Figure 6.(a) Scatter plot of the normalised PC1 and PC2/WACE trend of winter (DJF) Eurasian SAT for the 145 multi-model ensemble members (each marker represents the ensemble members of a particular model) in ALL (dots) together with the joint probability density function (JPDF) of the PC1 and PC2/WACE SAT trend in ALL (in colour shading) and in SICclim (in blue contour).The blue and red dots/lines are the respective ensemble means/ensemble standard deviations of the PC trends in SICclim and ALL.JPDF units are in percent.The grey dashed lines show the 95% confidence limit of the trend for ensemble members.(b) Scatter plot of the PC1 of winter (DJF) Eurasian SAT and AO trends in ALL for the 145 multi-model ensemble members with the same colour coding as in (a).A regression line (in black) indicates the relation between the trends in PC1 and AO.The black star in (a) and (b) is for the trends in the same for ERA5, and the respective black lines show the 95% confidence intervals.
(a)).Similar results are found in our central Eurasia region of interest (indicated by the black box in figures 2(c) and (d)), where the observed trend is outside the spread of the model simulated trends in both ALL and SICclim (figure S7(b)).This excessive warming reflects onto the PC1 positive trend in both experiments and is enhanced in ALL under observed Arctic SIC forcing (figure 6(a)).
• N during winter (Tyrlis et al 2020).The Siberian High (SH) index is the seasonal mean of area-averaged monthly 500 hPa geopotential height over 60 • N-85 • N and 30 • E-120 • E (black box region in figure 1(b)).The area represents the location where the 500 hPa geopotential height closely associates with the WACE pattern of variability (grey contours in figure 1(b)) and follows the region shown in a previous study as the footprint of SH at 500 hPa (Sun et al 2021).