Barents-Kara sea ice and European winters in EC-Earth

The potential link between decreasing Barents-Kara sea ice and cold winters in Europe is investigated using the enhanced resolution (horizontal atmospheric resolution of ∼80km\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 80\,\hbox {km}$$\end{document}) global, coupled climate model EC-Earth. Nudging sea ice only in the Barents-Kara Seas, five configurations of sea ice covers are used to assess the importance of the amount of sea ice in this region. Nudging in the coupled model is achieved by modifying the non-solar surface heat flux into the ice/ocean interface. The mean winter temperature response suggests a weak but statistically significant non-linear response with cooling over eastern Europe for moderate sea ice reductions in the Barents-Kara Seas, a weaker but still cold anomaly for minor reductions and warming for major reductions. However, this non-linear response is not reflected in the circulation. Instead, a negative mean sea level pressure anomaly over Barents-Kara Seas intensifies with sea ice reduction. In contrast to this, is the response in the coldest winters over central Europe: the larger the sea ice reduction, the stronger the Scandinavian pattern and the associated easterlies need to be to obtain cold winters over central Europe. The use of a coupled climate model is a potential explanation for the link between the intensified Scandinavian pattern and the cooling over Europe seen in this study, that is not observed in some atmosphere-only model studies.


Introduction
Arctic sea ice has declined over the last decades (Stroeve and Notz 2018). This decline has been observed over the whole Arctic, however, it is fastest in the Barents-Kara Seas (BAKA) (Onarheim et al. 2015(Onarheim et al. , 2018. At the same time a number of unusual cold winter events have been observed over Eurasia (Cohen et al. 2012;Liu et al. 2012). These apparent counter-intuitive developments have induced a number of hypotheses related to a possible causality between declining Arctic sea ice and cold extreme winters over Eurasia. In particular, how a rapid decline of BAKA sea ice could have an impact on the latter. Several studies based on observations and model results have suggested that there is a link between low sea ice in BAKA and cold winter conditions over Eurasia. According to Mori et al. (2014), the probability of severe winters in central Eurasia has more than doubled due to sea ice reductions in BAKA. Yang and Christensen (2012) analyzed the CMIP5 multimodel ensemble and found that despite an overall warming of Europe in the near future, cold winters (compared to the period  are still likely to occur. Some observational studies find a link between BAKA sea ice loss and cold Eurasian winters. Some of these even suggest Arctic sea ice loss as a driver of cold Eurasian winter. Koenigk et al. (2016) used reanalysis (ERA-Interim) and satellite sea ice data (OSI-SAF) to pinpoint autumn Barents Sea sea ice variations as the best predictor for the following winters NAO (low sea ice linked to negative NAO). The observed warming over BAKA has also been linked to cooling over East Asia (Kug et al. 2015). Outten and Esau (2012) showed a co-variability between the observed negative temperature trends over Eurasia and the Arctic sea ice loss. However, co-variability does not show causality.

Electronic supplementary material
The online version of this article (https ://doi.org/10.1007/s0038 2-020-05174 -w) contains supplementary material, which is available to authorized users.
Determining a potential driver of the Eurasian wintertime cooling was undertaken by Kretschmer et al. (2016). By using causal effects networks to analyze observational data, they showed that the BAKA sea ice reduction induces a weaker winter AO leading to wintertime cooling, thus being an important driver of the mid-latitude winter circulation. Further, Kretschmer et al. (2018) showed that cold winters in Eurasia and in North America are governed by different stratospheric patterns. Using cluster analysis, they showed that cold winters in northern Eurasia are linked to a weak polar vortex, upward propagating waves into the stratosphere leading to a negative NAO and resulting in cold spells over Eurasia.
Not only observational studies show this potential link and pathway. Using climate models, several model studies have investigated the effect of sea ice reduction on Eurasian winter temperatures. Petoukhov and Semenov (2010) used an atmosphere-only model (AGCM) at a coarse resolution ( 2.8 • × 2.8 • ) to investigate the effect of stepwise sea ice reduction in BAKA. They found a non-linear response to sea ice reduction, showing that a moderate sea ice reduction (40-80%) induced an enhanced anticyclonic anomaly over the Arctic Ocean and thus a strong easterly flow leading to cooling over Europe. In an AGCM study by Nakamura et al. (2015), sea ice reduction in BAKA induced a stationary Rossby wave leading to a negative AO/NAO like pattern, inducing prolonged Arctic cold air advection toward continental mid-latitudes.
Several mechanisms linking Arctic sea ice reduction, and in particular BAKA sea ice reduction, to cold Eurasian winter anomalies have been proposed. A typical proposed chain of events is: negative sea ice anomalies over BAKA in late autumn cause a larger surface heat flux from the ocean to the atmosphere, hence warming over BAKA; this induces upward propagating waves, leading to a weaker polar vortex, as well as weaker westerlies; downward propagation of the signal from the stratosphere back to the troposphere; the weaker polar vortex induces a negative AO/NAO like pattern linked to cold anomalies over Europe (Ruggieri et al. 2016;Garcia-Serrano et al. 2015;Kretschmer et al. 2018).
This potential link is still being widely discussed, as other studies have not been able to identify such a link. McCusker et al. (2016) compared the trend in BAKA sea ice decline to the trend of Eurasian winter temperatures for two time periods: 1979-1989 and 2002-2012. They found that although the Eurasian cooling from the first to the latter period was exceptional, the BAKA sea ice loss was not, hence suggesting that BAKA sea ice is not the main driver of Eurasian cooling. Instead, the observed cooling over Eurasia was most likely driven by internal atmospheric variability. The same conclusion was reached by Sorokina et al. (2016), based on ERA-Interim with no apparent correlation between Barents sea ice, turbulent heat fluxes over Barents Sea and the 'warm Arctic-cold Continents' pattern. Using an AGCM, Screen (2017a) reached a similar conclusion analyzing the link between NAO and Arctic sea ice loss and found that even though Arctic sea ice loss does intensify negative NAO events, this does not have a strong effect on the winter temperatures over Europe: the winter temperatures remain constant or even warm slightly. This 'missing' cooling is the result of the dynamical cooling due to the negative NAO being offset or exceeded by the thermodynamic warming due to the sea ice loss itself. Ogawa et al. (2018) also used AGCMs to investigate the impact of recent Arctic sea ice decline on winter climate, reaching the same conclusion as Screen (2017a) that the observed cooling over Siberia in winter is caused by internal atmospheric dynamics.
The discussion above clearly shows that there is no consensus among neither observational nor climate model studies as to whether there is a link between BAKA sea ice reductions and Eurasian wintertime cooling. The discrepancy between observations and model studies might be caused by an underestimation of the cooling in climate models (Mori et al. 2019). The discrepancy among modelling studies might be related to the use of AGCMs and coupled GCMs. As sea ice is often prescribed to AGCMs as new boundary conditions, atmosphere-ocean feedbacks are not included in the model. Further, the model sensitivity to sea ice perturbations may be offset due to circulation biases in the basic state. Deser et al. (2015) showed that including the oceanic feedbacks by coupling the atmosphere to the ocean results in global responses, as opposed to locally, poleward confined responses in AGCMs. Additionally, Screen (2018) states that coupled atmosphere-ocean models seem to agree more across models than atmosphere-only models. This suggests that coupled climate models are needed to obtain a more complete picture. As opposed to many of the studies presented above using AGCMs, we aim at using a coupled climate model, as this will capture the atmosphere-ocean feedbacks, thus capturing remote responses as well.
The aim of this study is to answer the following questions: 1. Using a coupled climate model and only modifying sea ice in the Barents-Kara Seas, are we able to find resulting colder winters in Europe? 2. Is the response, if identified, linearly related to the extent of the sea ice cover in Barents-Kara Seas? Or, is it nonlinear as Petoukhov and Semenov (2010) suggest? Screen (2018) discussed some of the methods used to modify sea ice in a coupled climate model. Sea ice can be modified by changing the sea ice concentration (SIC) in the model, overwriting the models own SIC using brute force. However, this is neither energy-nor (fresh-)water conserving. This can be improved by changing SIC over time, using assimilation and Newtonian relaxation (e.g. Tietsche et al. 2013). A different approach of energy and water conservation was used by Petrie et al. (2015) and Semmler et al. (2016). Both studies reduced sea ice thickness just before, or early in, the melt season, which induced a smaller summer sea ice extent. However, this approach only works for short time periods, as the system will quickly return to normal again. To avoid these constraints, some studies have used a 'ghost' flux to perturb sea ice in coupled climate models. Deser et al. (2015) adds to the ice module a seasonally, but not spatially, varying long-wave (LW) radiative flux (LRF), which is computed according to a target sea ice cover. They simulated scenarios with different LRF and from the resulting sea ice covers, estimated which magnitude of LRF results in the target sea ice cover. Oudar et al. (2017) compute the 'ghost' flux as the difference in non-solar (LW+turbulent) heat flux between a control simulation and target simulation. That is, they use the non-solar flux generated due to a specific sea ice cover. A slightly different approach is taken by McCusker et al. (2017). Their 'ghost' flux is computed as the amount of heat needed to melt (or create) a certain mass of sea ice. This heat is then added to/removed from the ice model. To form sea ice, heat is also removed from the ocean, in order to cool it down.

Methods for constraining sea ice
Our approach is similar to McCusker et al. (2017): the difference in sea ice concentration in a specific grid cell at each time step is used to compute a 'ghost' flux. However, our method differs in two ways: (1) instead of computing the heat needed to melt/form a certain mass of sea ice, we use the difference between the instantaneous SIC and a target SIC ( SIC = SIC model − SIC target ) as a weight on the magnitude of the instantaneous non-solar heat flux received by the ocean (HF): where A is an empirical constant. We have used a value of 5, as the resulting HF were too small to make a difference when using a smaller value. Using 5 ensures a large enough adjustment every time step, without causing the model to become unstable. If SIC exceeds a threshold value (here set to 5% ), the correction term is added to the instantaneous non-solar flux, HF.
(2) The heat added to melt sea ice is superimposed at the atmosphere/ocean boundary.

Experimental setup
In this study we use the fully coupled climate model EC-Earth V3.2. The atmosphere component of EC-Earth used here is the ECMWF's Integrated Forecast System (IFS) (1) HF = HF ⋅ SIC ⋅ A cycle 36r4 with a spectral resolution of T255 (equivalent to 80 km) and 91 vertical levels. The ocean component is the Nucleus for European Modelling of the Ocean (NEMO) version 3.6 on the ORCA1 tripolar grid ( 1 • × 1 • ) and with 75 vertical levels. NEMO is embedded with the sea ice component the Louvain-la-Neuve sea Ice Model (LIM3). IFS and NEMO are coupled through the OASIS-3 coupler. The older version EC-Earth V2.2 is presented and validated in Hazeleger et al. (2010) and. V2.2 differs from V3.2 by using IFS cycle 31r1 with horizontal resolution T159, NEMO version 2 and LIM2. These have all been updated since V2.2. One of the main differences is in the sea ice component LIM, as LIM3 has five sea ice thickness categories whereas LIM2 only has one. Validating EC-Earth V3.2 against ERA-Interim, similar to Hazeleger et al. (2012), show a general warm shift and reduction of the biases in most regions in V3.2, except over Sahara where the cold bias in V2.2 seems intensified in winter (DJF) (see the Supplementary figure S1). The systematic errors in the Arctic has shifted from a large cold bias in V2.2 to a slightly warm bias. The bias in summer (JJA) is also smaller in general than in V2.2, with a shift toward a warmer bias. Bias in the tropics have been reduced compared to V2.2. The cold bias over Greenland has decreased, as has the bias over Europe. A band across central Eurasia has a warmer bias in V3.2. The mean sea level pressure bias in V3.2 show similar patterns to that in V2.2 (Supplementary figure S2). In general, the bias shifted towards a more positive bias, i.e. higher pressure. To focus on a few areas, we note the bias is smaller over the North Atlantic in both summer and winter relative to V2.2. The high pressure bias has increased over Sahara. The Aleutian Low bias has shifted from being positive to negative. Hence, V3.2 have similar bias patterns to those of V2.2, with smaller, but slightly more positive values.
Six simulations were carried out, each 100 years long: one freely running model and five experiments, where sea ice concentrations were nudged. The freely running, unforced model simulation under present day conditions (GHG, aerosols etc) was run after an initial spin-up period of 350 years, hereafter referred to as CFree. This model simulation is used to generate the model sea ice concentration climatology ( CFree Clim SIC ) and the sea ice area climatology ( CFree Clim SIA ). There is still a drift (global mean surface (2m) temperature trend of ∼ 0.2 • C over 100 years) in CFree even after the spin-up. All model results are therefore linearly detrended before being analyzed. The five perturbed experiments are equal to CFree except the sea ice in BAKA (30 • -80 • E , 65 • -80 • N , see Fig. 1), where the SIC are nudged. Outside this region there is no nudging and the sea ice can evolve freely. CFree as well as the five perturbed simulations all show large centennial variability in the Atlantic Meridional Overturning Circulation (AMOC). It is assumed that the centennial variability of the perturbed experiments are similar to that of CFree, as they start from the same initial conditions and it is assumed that the nudged sea ice changes does not affect this. This is a potential limitation of the study, as only one member for each perturbation is performed and starting from a different initial state may result in different 100-year simulations. When nudging is applied, the SIC in BAKA is only reduced in the winter half year from November through April, as November sea ice have been shown to be correlated to the sign of the NAO (Koenigk et al. 2016). The rest of the year, SIC is nudged toward CFree Clim SIC . Following the setup from Petoukhov and Semenov (2010), SIC is stepwise reduced from 100% to 20% in steps of 20% in each simulation. This creates five perturbed sea ice covers for BAKA, all originating from the same initial conditions from CFree: C100 nudged toward 100% of CFree Clim SIC , C80 toward 80% (i.e. C80 = 0.8 ⋅ CFree Clim SIC ) and so forth. Comparing the long-term mean surface temperature for C100 and CFree, it is clear that only the local, perturbed region is affected by removing the internal SIC variability (not shown). Inside BAKA, C100 is up to 2 • C warmer than CFree for the winter (DJFM) mean, but outside this region the differences are ∼ ±0.5 • C . Comparing C100 SIC to the observed March SIC for the period 1979-2009 (OSISAF 2017) it can be seen that C100 replicates the observed sea ice cover quite well, although with a smaller sea ice extent in the Barents Sea, the Bering Sea and the Sea of Okhotsk, see Fig. 1. The monthly mean achieved sea ice area (SIA) for BAKA for each of the five perturbation experiments relative to the target SIA for C100 ( CFree Clim SIA ) are presented in Fig. 2. Generally, the achieved winter SIA is larger than the target, and the summer SIA is lower, resulting in an annual mean close to the target SIA. The large difference between target and achieved SIA in summer is related to the fact that BAKA is almost ice-free in summer. Therefore, differences between the two concentrations will be relatively larger in summer. The achieved SIA is not equal to the target, as the nudging method compares every grid point to the target and is only activated if the difference is larger than ±5% . Although the target SIA is not reached exactly, Fig. 2 shows a clear reduction going from C100 to C20. Reducing sea ice in BAKA warms the local ocean underneath, and increases the surface salinity in summer and autumn for all four perturbed runs relative to C100. In winter and spring, the BAKA ocean surface layers are fresher than C100 (not shown). This is related to the reduced melting/ forming of sea ice, leading to less fresh melt water release in summer and autumn, and less brine rejection in winter and spring, i.e. less salt being pushed out of the ice and into the ocean.
Changing water mass properties in the Nordic Seas, or changing the atmospheric circulation are important for the AMOC. Changing sea ice like we have done could possibly do the same, as suggested by several studies (Otterå et al. 2004;Zhang et al. 2011;Sévellec et al. 2017). The latter study found that a decline in Arctic sea ice was followed by a weakening of the AMOC. Figure 3 depicts the long time averaged annual mean AMOC and its ± standard deviation for the five experiments. There is large interannual to decadal variability of the same order of magnitude in the AMOC in all five experiments, while the differences in the long time average are small (maximum mean difference was 0.51 Sv between C100 and C60). However, statistical tests (two-sided Kolmogorov-Smirnov test) on both the annual mean time series and its decorrelated time series reveal that the AMOC in C60 and C20 are significantly weaker than in C100 at the 95% confidence level. As the AMOC time series . Colors indicate the five simulations: black (C100), dark blue (C80), light blue (C60), orange (C40) and red (C20). 100% is equal to the climatological BAKA SIA in CFree ( CFree Clim SIA ), i.e. the target for C100. Dashed lines indicate the target percentage were autocorrelated (AR(1) processes) they were decorrelated using the decorrelation time D (von Storch and Zwiers 2003, chapter 17). The fact that the AMOC in the other simulations do not show significant difference from that in C100 suggests the linkage between the slowdown of AMOC and the reduction of BAKA sea ice are not linear.
The atmospheric response of the reduced sea ice cover is assessed as the difference between the four experiments with reduced sea ice (C80, C60, C40, C20) and the control run C100. That is, the climatology of C100 have been subtracted from all experiments, so as to remove the seasonal cycle. Hence, results for C100 are presented as anomalies wrt the 100-year mean and C20 to C80 as changes of anomalies. As focus is on winter responses, we analyze the three bimonthly periods: early winter (December and January, DJ), mid winter (January and February (JF) and late winter (February and March, FM). By using bimonthly periods instead of a mean over all winter months (DJFM), a potential cold signal of shorter duration will not that easily be averaged out by a compensating warmer period. This is essential, since very cold events are not necessarily cold throughout the entire winter. The coldest winters are represented by the 5th percentiles temperatures ( T p5 ), which is computed for each grid point as well as for the area mean of Central Europe (CEU, They are based on the monthly values of the detrended temperature anomalies for each subseason, i.e. the subseason JF consists of 200 months (2 months for each of the 100 years).
Statistical significance of the mean and percentiles are assessed by bootstrapping using Monte Carlo simulations with resampling (von Storch and Zwiers 2003, chapter 5.5). Statistical significance is given at the 95% confidence level. The significance of the percentiles are computed by bootstrapping the percentiles of the C100 experiment, then testing the percentiles of the sea ice reduced experiments toward this.

Results
Results are discussed below for all three bimonthly periods (early, mid and late winter), but only results from mid winter are shown. In general, anomalies relative to C100 are shown (Fig. 4).

Mean atmospheric responses
We first focus on examining the change of the mean states as sea ice is reduced. Reducing sea ice in a certain area allows for more heat to escape from the ocean into the atmosphere. Therefore, the surface air temperature above and in the surrounding area is expected to increase. This is confirmed by the surface temperature distribution (Fig. 4). Comparing to the temperature mean of C100, it is seen that over the BAKA nudging region there is an increase in temperature of 3 • C or more. The more the sea ice is reduced, i.e. going from C80 to C20 , the larger this area becomes. For all four reductions and for each bimonthly period, the general picture is of a warming Northern Hemisphere (NH), largest increase in temperature over BAKA and adjacent areas of Eurasia and North America. Standing out is the simulation C60 with 60% of climatological sea ice, which shows significantly lower temperatures over an area of eastern Europe between Scandinavia and the Black Sea of about 1 • C as well as a larger area over the North Atlantic (around < 0.5 • C ), as seen in Fig. 4c. This implies a wider cold response. A smaller, but still significantly cold area over eastern Europe is present in late winter as well (not shown). For smaller sea ice reductions than this (C80), there is a lack of warming, but not statistically significantly different from the control run anywhere over Europe (Fig. 4b). Larger sea ice reductions (C40 and C20), does not exhibit colder areas over Europe, but rather statistically significant warming (Fig. 4d, e). It is interesting to note that the warmest winters appear to occur in C40, not in C20. The non-linearity of the surface temperature with BAKA sea ice reductions for all three bimonthly winter periods is visualized in Fig. 5, which shows the area mean temperature anomalies (i.e. Fig. 4) for two areas over Europe (areas shown in Fig. 1): CEU (black circles) and EU (green squares) which is the same area used in Petoukhov and Semenov (2010). The colder temperatures in C60 relative to C80, C40 and C20 are present in all three bimonthly periods and for both areas over Europe. For these area averaged values, only C20 and C40 are significantly different (warmer) from C100. However, the difference between the minimum and the maximum temperature (i.e. C60 and C40) is statistically significant for both EU and CEU in all three subseasons. Averaging over a larger area suppress the Fig. 3 The 100-year mean of annual AMOC (in Sv) for the latitudinal band between 45 • N and 48 • N (filled circles) for each of the five simulations: red (C20), orange (C40), light blue (C60), dark blue (C80) and black (C100). Errorbars indicate +∕− one standard deviation. C60 and C20 are significantly different from C100 on the 5% level using the Kolmogorov-Smirnov test Fig. 4 a Mid winter (JF) mean absolute surface (2m) temperature (in • C ) for C100. The difference in mean JF surface temperature compared to C100 for b C80 , c C60 , d C40 and e C20 . Black box indicates BAKA area. Dots indicate significance at 5% level, black dashed box indicate CEU region and green box indicate EU region Fig. 5 The area mean of the 100-year mean surface air temperature anomaly (in • C ) for two areas over Europe for each of the five perturbed simulations. The two areas CEU (black circles) and EU (green squares) are shown in Fig. 1. Filled markers indicate temperatures significantly different from C100 at a 5% level, i.e. only C20 and C40 are significantly different from C100. Significance was tested using bootstrapping with resampling, see Sect. 2.2. The temperature difference between C60 and C40 (i.e. minimum and maximum) is significantly different at a 5% level for both EU and CEU for all three subseasons relatively small statistically significant colder area over eastern Europe seen in Fig. 4 for C60. The mean absolute mean sea level pressure (MSLP) for the control simulation (C100) is shown in Fig. 6a. Similar to the observed winter MSLP pattern, two low pressure systems are identified over Iceland and the Bering Sea, while a high pressure system is located over Siberia and central Asia. The mean response to sea ice reductions for all four perturbed experiments are quite similar: a low pressure anomaly appear over the BAKA region and downstream thereof (Fig. 6b-e), with increase in strength as sea ice reduces. This tendency seems to saturate around a 40% sea ice reduction (C40), consistent with the maximum warming response in C40 (Fig. 4d). The low pressure anomaly is consistent with removing sea ice, which warms the surface and the air above, hence leading to expansion and rise of the air (a 'heat low'). It results in an overall southward displacement of the high pressure system over Siberia in conjunction with the development of the local pressure anomaly over and downstream of the forcing region. The linearity of the MSLP anomaly across reduction scenarios can be seen more clearly in Fig. 7. It shows the area mean root-mean-square-anomaly Fig. 6 a Mid winter (JF) mean absolute mean sea level pressure (in hPa) for C100. The difference compared to C100 for b C80 , c C60 , d C40 and e C20 . Dots indicate significance at 5% level (RMSA) of the MSLP for each perturbed experiment relative to C100 for the area north of 50 • N , encompassing the MSLP anomaly seen in Fig. 6b-e. In general, we see an increasing RMSA following the reduction in sea ice, going from 0.7 for C80 to 1.5 for C40 and C20 in JF. Hence, the MSLP anomaly in Fig. 6 becomes increasingly less like the C100 pattern as BAKA sea ice is reduced. Similar to the mean temperature response, the MSLP seems to saturate around C40.
The changes in the pressure systems are reflected in the wind patterns (not shown): the control run C100 shows the dominating westerlies over Europe, and all four perturbed experiments show a cyclonic wind anomaly over BAKA, associated to the low pressure anomaly discussed above. Also, as the low pressure anomaly increases (from C80 through C20 ), the wind anomaly over Europe increases. The results above indicate a potential link between cold winters and moderate sea ice reductions in BAKA.

Abnormally cold winters
After examining the mean atmospheric responses through all 100 years we will take a closer look at the coldest winters. There are several ways to investigate the response in the occurrence of the coldest winters. Here we first analyze the temperature distribution in each grid point, to identify geographical regions of interest. Afterwards, we focus at larger regions covering Europe.

Grid point wise percentiles
The coldest winters are shown as the 5th percentile temperature ( T p5 ) in each grid point in Fig. 8. Figure 8a shows T p5 for C100 ( T C100 p5 ). The difference between the control simulation and the four perturbed experiments are shown in Fig. 8b-e for mid winter. Clearly, large parts of the Northern Hemisphere experience warmer cold winters (orange colors), especially over the perturbed region (BAKA). The magnitude and spatial extent of these warmed areas increase as more sea ice is removed in BAKA (i.e. going from C80 to C20 ). This is similar to the 100-year mean atmospheric response seen in Fig. 4. However, there are also areas where T Cx p5 is colder than T C100 p5 . In C60 (Fig. 8c), Europe experi-ences a lack of warming, with small areas of statistically significant cooling over Eurasia, including parts of the area over eastern Europe that was significantly cold in the mean response (Fig. 4c). The temperature anomalies over eastern Europe reach 1-2 • C lower than that of the climatological run (C100). The 5% warmest winters demonstrate a similar pattern (not shown), with a lack of warming over Europe in C60. Combining this with the mean temperature response shown in Sect. 3.1 indicates that the temperature distribution for C60 is shifted toward lower temperatures compared to C100. This shows that C60 is statistically significant colder than C100 in some small areas and colder, but not statistically significant colder, in larger areas. This suggests that a moderate reduction of sea ice may link to statistically significant lower temperatures, but only for some small areas of Europe. For C80, there is not significant change over Europe, whereas C40 and C20 show a statistically significant warming over Europe. Quantifying the area averaged coldest temperatures over CEU and EU, the same non-linearity as for the mean response is implied in Fig. 9, but neither of the four reduction experiments are statistically significantly different from C100. However, the temperature difference between the minimum and the maximum (i.e. C80 and C40) is statistically significant for EU, but not for CEU, for JF and FM, indicating the non-linear behavior. The analysis in Fig. 8 is performed at each grid cell separately, meaning that the coldest temperature in neighboring grid points are not necessarily from the same winter months, thus not necessarily caused by the same weather systems. Lining up all the winter months chronologically and extracting the coldest month in each grid point reveals larger coherent areas, as depicted in Fig. 10. Further, looking at the second, third, fourth (and so on) coldest winter months results in a much messier picture with no large coherent areas showing up (not shown). Therefore, the above figures suggest that there might be a link between reduced BAKA sea ice and Eurasian winter temperatures, but it does not show why they occur.

Cold winter months composite of larger regions
Focusing on larger regions, the area mean winter temperature for CEU is used to determine the 5% coldest winter months, in each bimonthly period. These months are then Fig. 7 The root-mean-squareanomaly (RMSA) of the mean sea level pressure (in hPa) for each perturbed experiment relative to C100 over the area north of 50 • N used to construct composites of surface temperature and MSLP. The composite surface temperatures for cold CEU are shown in Fig. 11. The resulting composites of surface temperature anomalies for C100 demonstrates a cold anomaly of more than −5 • C in a large area of Europe and Western Asia, accompanied by a warm anomaly of above 2 • C over Greenland (Fig. 11a). Reducing sea ice in BAKA results in a general warming over most of the Arctic for the cold winter months for all three subseasons (Fig. 11b-e). However, it is worth noting that all four perturbed experiments have a similar, cold response in early and mid winter over the northern part of CEU and Scandinavia, as well as parts of Asia and North America, with respect to C100 (early winter not shown). The southern part of CEU warms relative to C100 for C40 and C20. The area mean temperature response over CEU show no statistically significant change as the BAKA sea ice is reduced (blue circles in Fig. 9).
The corresponding circulation responses are shown in Fig. 12. For the 5% coldest winter months in CEU in Fig. 12a, the overall MSLP is higher than the time average over the Nordic Seas and Scandinavia in C100, resembling the Scandinavian pattern (Barnston and Livezey 1987). Reducing the Fig. 8 a Mid winter (JF) composite of the surface temperature (in • C ) for the 5% coldest winters at each grid cell for C100. The difference compared to C100 for b C80 , c C60 , d C40 and e C20 . Orange/ brown colors indicate warmer 5th percentiles (warmer coldest winter months) and blue colors indicate colder 5th percentiles (colder cold winter months) in the perturbed experiments wrt to C100 sea ice in BAKA expands this high pressure anomaly northward while strengthening the anticyclonic anomalies, saturating around C40. This implies that as sea ice is reduced, cold winters in CEU are accompanied by an increasingly stronger Scandinavian pattern and hence stronger easterlies bringing cold air from Siberia into Central Europe. The strengthening of the Scandinavian pattern can be visualized by looking at the projection of the reduced experiments (C20-C80) MSLP anomalies onto the control simulation (C100) MSLP anomalies. Equation (2) shows the case for C80 where and are the longitudinal and latitudinal coordinates, MSLP C80 is the temporal mean of the detrended MSLP anomaly (i.e. Fig. 12b) and MSLP C100 is the same but for C100. The sum is taken over the region ( 40 • W -60 • E , 45 • N-90 • N ), the green box in Fig. 12a. The result for each experiment and for each of the 3 bimonthly winter periods can be seen in Fig. 13. For DJ and JF, P are generally above 1, indicating that the Scandinavian pattern intensifies when sea ice is removed in BAKA. For C40, the high pressure anomaly over Scandinavia has shifted northward and expanded (Fig. 12d), leading to a P value closer to 1 compared to C20, C60 and C80. This high pressure anomaly over Scandinavia is in contrast to the low pressure anomaly over BAKA observed for the mean response. This imply that even though the mean response is a low pressure anomaly, the 5% coldest winters (10 months out of the 200 months for each bimonthly period) occur when cold air is being advected from the North, associated with the Scandinavian pattern.

Discussion
The study presented here investigate the responses to the sea ice reduction by performing similar experiments as Petoukhov and Semenov (2010), but in an atmosphere-ocean coupled model system. Our mean non-linear temperature response was similar to that of Petoukhov and Semenov (2010). However, this non-linearity was not found when focusing on the 5% coldest winters. Neither was it found in the circulation responses. The difference in the mean responses between their results and ours could be model related: firstly, they used an AGCM with fixed SST patterns, whereas we used a coupled model. Although AGCMs have smaller biases relative to observations, due to AGCMs being forced by observed SST's, coupled models capture the atmosphere/ocean feedbacks. These feedbacks have been proven to be important for Fig. 9 The area and temporal mean surface temperature (in • C ) of the 5% coldest winter months in each grid cell [black circles (CEU) and green squares (EU)] and the mean of the 5% coldest winter months in CEU (blue circles). Filled circles/squares indicate significance at 5% level (see Sect. 2.2), i.e. none of the four perturbed experiments are significantly different from C100. The temperature difference between C80 and C40 (i.e. minimum and maximum) is significantly different at a 5% level in JF and FM for EU, but not for CEU Fig. 10 The month with the coldest surface air temperature occurring during 100 years at each grid cell for C60. The numbering is constructed by taking the four winter months (DJFM) and lining them up chronologically. Hence, month #1 is January of the first year, #2 is February, #3 is March and #4 is December etc capturing global wide responses to local changes in the climate system (Deser et al. 2015). Part of this may be related to a stationary Rossby wave driven feedback loop involving tropical SST's in the Pacific (not shown). An example of this feedback system is the non-linear response with a widespread cold signal over Europe and the North Atlantic seen in C60 and partly C20 for the air temperature. This non-linear response is manifested in the sea surface temperature (SST) as well, with significantly colder areas in the North Atlantic for C60 and a small, but significantly, colder area in C20 (Fig. 14). The SST changes are consistent with the AMOC changes seen in Sect. 2.2 and in Fig. 3, due to the linkage between the AMOC and North Atlantic SST (Schmith et al. 2014). At this point we cannot say that the AMOC changes are solely caused by the sea ice reduction. However, the AMOC in all five experiments has a similar level of inter-annual to decadal variability. This suggests that the sea ice changes, and not the inter-annual to decadal variability, is the main driver of the AMOC changes. The mean surface air temperature changes over Europe observed in C60 thus are consistent with the ocean response.
Secondly, our horizontal and vertical model resolution is higher (a horizontal resolution of 0.7 • and 91 vertical levels Fig. 11 Mid winter (JF) composite of the surface temperature (in • C) for the 5% coldest winters in CEU (black box). a Temperature difference in C100 ( T C100 − T C100,clim ) and the response as changes of anomalies in b C80 , c C60 , d C40 and e C20 . Orange colors indicate warmer than C100 and blue colors colder Fig. 12 Mid winter (JF) composite of mean sea level pressure (in hPa) difference wrt the climatology of C100 for the 5% coldest winters in CEU (black box) for a C100, b C80, c C60, d C40 and e C20. Green box indicate the region used to compute P in Fig. 13 Fig. 12a. P is unitless. The method for computing P is shown in Eq. (2). P = 1 indicate that the Scandinavian pattern in C20 to C80 is equal to the one in C100. P > 1 indicate a stronger Scandinavian pattern relative to C100, and P < 1 indicate a weaker Scandinavian pattern in T255 compared to 2.8 • in T42 and 19 vertical levels). In our setup, both aspects of the model configuration are consistent with a more adequate representation of physical processes.
Another potential reason for the discrepancies is the experimental protocol. The target sea ice reductions were not reached exactly, as would have been the case in an AGCM. Taking this into account, we still interpret our results to offer a more realistic sensitivity of real-world response to BAKA sea ice modifications.

Thermodynamic warming vs. dynamical cooling
Our results contribute to the ongoing discussion on whether there is a link between Arctic sea ice loss and cold Eurasian winters. Although we only found a weak but significant link for moderate BAKA sea ice loss and mean winter temperatures, this differs from most recent studies which does not find a link (i.e. Blackport et al. 2019), illustrating the discrepancies among AOGCMs. Among AGCMs the discrepancies are larger (see Sect. 1). This difference between AOGCMs and AGCMs might be caused by the balance between direct thermodynamic warming and indirect dynamical cooling (Screen et al. 2015b): melting sea ice leads to direct thermodynamic warming above the area but it might also induce a negative AO/NAO, leading to dynamical cooling. This battle have been used to explain the 'missing' cooling response to a negative NAO observed in several studies (Screen 2017a;McKenna et al. 2018), indicating that the thermodynamic warming overcomes the dynamical cooling. Both these studies used AGCMs. This potentially impacts the result, as the balance between thermodynamic warming and dynamical cooling might depend on the choice of experiment (observational, AGCM or AOGCM study) and to the location and magnitude of the sea ice loss, explaining the discrepancies between our results and other AOGCMs studies.

Observational studies, AGCMs or AOGCMs
The discrepancies between observational studies, AGCMs and AOGCMs might be related to the fact that they are observational studies, AGCM studies and AOGCM studies. Mori et al. (2019) showed that AGCMs systematically underestimated the Eurasian cooling due to a weaker 'warm Arctic-cold Eurasia' (WACE) pattern, relative to observations. They attributed the weaker WACE pattern to a too weak response to BAKA sea ice loss. Corrected for this bias by enhancing the WACE pattern, they conclude that ∼ 44% of the central Eurasian cooling observed from 1995 to 2014 was caused by BAKA sea ice loss. This might explain the 'missing' cooling over Eurasia found in some model studies (McCusker et al. 2016).
The discrepancies between atmosphere-only and coupled atmosphere-ocean model studies might be explained by the balance between thermodynamic warming and dynamical cooling as well. Smith et al. (2017) found that the response to sea ice reduction in AGCMs was dominantly a thermodynamic response to local warming, whereas the dominant response in coupled models was the dynamical response related to the AO/NAO. They suggest that the dynamic response does occur in AGCMs, but somewhat delayed compared to AOGCMs. Screen (2017a) and McKenna et al. (2018) did not observe the cooling over Europe, explaining this with the thermodynamic warming winning over the dynamical cooling. Following the above, using AOGCMs should result in a stronger dynamical cooling, leading to a cold signal as we found. However, most recent studies do not find this response. This suggests that this balance not only depends on the choice of model.

Location and magnitude of sea ice changes
Several studies have shown that the atmospheric response to sea ice loss depends greatly on the location of the sea ice loss (Koenigk et al. 2016;Pedersen et al. 2016). In particular, Screen (2017b) showed that some Arctic regions induce large-scale dynamical changes while other regions only induce a thermodynamic response from the local sea ice loss. One such region was the BAKA where sea ice loss induced a negative NAO in winter. Additionally, it was shown that the sum of the responses from each region is not the same as pan-Arctic sea ice loss response: the former leads to cooling over the mid-and high latitudes, while the latter leads to warming. They propose, the regional sea ice loss works on shorter time-scales and induces dynamical cooling, whereas the pan-Arctic sea ice loss is indicative of the longer time-scale sea ice loss resulting in thermodynamic warming winning over the dynamical cooling. As the BAKA is the region currently experiencing the largest sea ice loss, it goes in the category of short time-scales (Onarheim et al. 2018;Stroeve and Notz 2018). Hence, dominated by dynamical cooling.
The importance of the short time-scale regional sea ice loss (leading to dynamical cooling) will diminish as the longer time-scale pan-Arctic sea ice loss (leading to thermodynamic warming) takes over. Hence, the response to the sea ice loss observed until now is a transient response and not the same response as will be seen for Arctic sea ice loss in the future, i.e. it is a non-linear response, as suggested by several studies including Screen (2017b) and Peings and Magnusdottir (2014). In our study, the region of the sea ice reductions might be too small for the thermodynamic warming to compensate for the dynamical cooling, hence showing a weak link for moderate sea ice reductions.

Limitations of the experimental setup
Sea ice was only reduced in the BAKA, not across the Arctic. Although pan-Artic sea ice loss would be more realistic, the balance between warming and cooling might not be (see Sect. 4.1.2). Thus, by reducing sea ice across the Arctic, the warming might overwhelm the already too weak cooling. To avoid this, we only reduced sea ice in BAKA, as this region has the largest influence on European winter weather. Our choice not to simulate the pan-Arctic sea ice loss is most likely the reason we do not observe the intensification of the Aleutian Low observed for the six coupled models in Screen (2018). Due to this limitation, our results can 'only' be used to imply a link between reduced BAKA sea ice and a temporal mean Eurasian cooling.
The method used for nudging sea ice is another potential limitation on our results. Sea ice is reduced in each grid cell to a certain percentage of the climatology. Another possible method for changing the sea ice cover in a certain region is to change the overall volume of sea ice in each grid cell. This would allow for sea ice to cover a grid cell by 100% which is more physical than the entire region being covered by up to, and not exceeding, for example 60% sea ice. However, this would require removing sea ice in certain regions and not others, i.e. making assumptions about the location of the sea ice reduction. Thus, reducing the sea ice concentration in each grid cell is less realistic, but avoid further assumptions. Such additional sensitivity experiments are however well beyond the scope of the present work.

Conclusions
The coupled atmosphere-ocean model EC-Earth was used to investigate the potential link between BAKA sea ice loss and cold European winters. This was achieved by reducing sea ice in BAKA, keeping it everywhere else. Our 100-year mean results differed from those focusing on the 5% coldest winters over Europe. For the mean winter temperatures we found a weak but statistically significant non-linear response with cooling over eastern Europe for moderate sea ice reductions in the BAKA region, a weaker but still cold anomaly for minor reductions and warming for major reductions. This non-linear response was not identified in the circulation changes. Instead, a negative mean sea level pressure anomaly over BAKA intensified with sea ice reduction. In contrast to this was an intensification of the Scandinavian pattern and associated easterlies with increasing sea ice loss. This was needed to obtain cold winters over central Europe.
Our results do not clearly support those studies finding a link. Nor does it clearly support those not finding a link, as we see a weak but significant link, but only for the mean temperature. Whether or not studies find a link between BAKA sea ice reductions and cold European winters might depend on the battle between direct thermodynamic warming from heat release over formerly ice-covered areas and the indirect dynamical cooling from inducing a negative NAO. The choice of model (AGCM or AOGCM) or observational study affects this, as models in general and AGCMs in particular underestimates the dynamical cooling. This is counteracted by the impact from the location and magnitude of the sea ice loss, as regional sea ice loss favors dynamical cooling and pan-Arctic sea ice loss favors thermodynamic warming. Hence, finding a potential link between sea ice loss and cold European winters is highly sensitive to the choice of model and area of sea ice loss.