Simulations of Eurasian Winter Temperature Trends in Coupled and Uncoupled CFSv2

Over the last 30 years, global sea surface temperatures (SSTs) have substantially increased (Hartmann et al., 2013) and Arctic sea ice has declined (Serreze et al., 2007; Comiso et al., 2008). Decreases in Arctic sea-ice cover further increase near-surface temperatures in the Arctic due to the lower albedo of the open ocean compared to sea ice, allowing greater absorption of solar radiation. As a result, temperatures have warmed more in the Arctic than in the lower latitudes ——a phenomenon known as Arctic amplification (Kumar et al., 2010; Serreze and Barry, 2011). Numerous studies have attempted to link Arctic amplification with changes in the observed midlatitude synoptic variability, such as the frequency of extreme events or winter temperature decreases over Eurasia and North America, which have been extensively reviewed in (Cohen et al., 2014), (Gao et al., 2015), and (Overland et al., 2015). Because of the large uncertainty in the literature, questions remain regarding whether the recent observed changes in extratropical and high-latitude variability are a result of changes to the external forcing (Arctic sea-ice decline, low-frequency changes in SSTs, changes in greenhouse gas concentrations, volcanic forcing, and solar variations), or are a manifestation of natural variability (defined as changes in various aspects of the extratropical and high-latitude circulation on different timescales, which are unrelated to changes in external forcing). These natural circulation changes can also have a notable impact on Arctic sea-ice loss (Ding et al., 2017). However, for the purposes of this study, the long-term changes in SSTs and sea ice are represented as external forcing.

The hypothesized pathway for Arctic amplification to impact midlatitude climate, as suggested by (Francis and Vavrus, 2012), is that the increased Arctic warming creates a reduced latitudinal temperature gradient between the equator and North Pole, reducing the zonal upper-level winds through the hypsometric relationship. This in turn allows for a slower progression of upper-level atmospheric waves and a wavier jet stream (Francis and Vavrus, 2015), leading to more persistent weather patterns, thus increasing the probability for the occurrence of extreme events. These impacts on upper-level patterns are argued to have a variety of atmospheric consequences, both local and non-local. Some examples include a decrease in temperature variance in the mid and high latitudes (Screen, 2014), reduction in the strength of cold-air outbreaks in the midlatitudes (Ayarzagüena and Screen, 2016), and increased winter Eurasian snowfall due to increased moisture availability from melted sea ice (Liu et al., 2012).

Other studies have illustrated that Arctic sea ice can force a negative phase of the Arctic Oscillation (AO), which causes cooler temperatures over the midlatitudes (Nakamura et al., 2015). However, (Mori et al., 2014) showed that the cooling over Eurasia is independent of the AO and more related to autumn sea-ice loss in the Barents and Kara seas. (Petoukhov and Semenov, 2010) also proposed a pathway linking sea-ice loss in the Barents and Kara seas to cooling over Eurasia, but showed the relationship is highly nonlinear in terms of the temperature response to different magnitudes of sea-ice loss, while stressing the importance of the possible role of atmospheric teleconnection patterns. (Sato et al., 2014) deduced that a northward shift in the Gulf Stream, which results from a weaker Atlantic Meridional Overturning Circulation, reduces the sea ice in the Barents Sea, which then yields altered planetary wave patterns that promote Eurasian cooling. Several studies (Honda et al., 2009; Inoue et al., 2012; Zhang et al., 2012; Tang et al., 2013) have claimed that sea-ice decline ultimately leads to a stronger Siberian high, which enhances wintertime cooling over Eurasia. Other studies have explored the issue of conditional dependence (Overland, 2016), suggesting that for Arctic amplification to impact midlatitude weather, there must be some preexisting condition.

Linkages between Arctic warming and midlatitude cooling have also been contested. For example, (Barnes, 2013) and (Screen and Simmonds, 2013) highlighted that calculated differences in atmospheric wave patterns are sensitive to the parameters used to define them, and that under different definitions the linkages would not be as conclusive. (Screen et al., 2014) and (Perlwitz et al., 2015) demonstrated that, although sea ice plays a large role in the increase of near-surface Arctic temperatures, the signal is lost at higher altitudes, which would not support any feedbacks to the midlatitudes. In addition, (Barnes et al., 2014) argued that there has been no increase in the frequency of blocking patterns in the last three decades, and a link to Arctic sea-ice loss is nonexistent. (Screen et al., 2015) stated that the probability of colder winters over the midlatitudes will decrease due to sea-ice loss, while (McCusker et al., 2016) and (Sun et al., 2016) concluded that the cooler Eurasian winters are most likely due to atmospheric internal variability.

The present study attempts to determine the ability of an uncoupled (atmosphere-only) model and a coupled (ocean-atmosphere) model in producing the cooling over Eurasia, to provide further insight into the subject. The main questions that we address are: (1) Do uncoupled runs forced with either SST or sea-ice boundary conditions from 2005-14 simulate the observed December-February atmospheric changes relative to simulations forced with 1981-90 boundary conditions? (2) Are changes prevalent in the modeled variability of wintertime temperatures over Eurasia, and are there systematic atmospheric patterns associated with warm and cold winters seen within the different model configurations for the external forcing? (3) Can initialized coupled-model runs simulate the observed characteristics and, if so, is there a dependence on the forecast lead time? Answers to these questions will provide insights into the assessment of whether the observed changes are due to sea-ice changes, SST changes, or more related to atmospheric internal variability.

The Climate Forecast System, version 2 (CFSv2), is a fully coupled atmosphere-ocean-ice-land model developed at the National Centers for Environmental Prediction (Saha et al., 2014). For this study, the atmospheric component of CFSv2 is used for uncoupled Atmospheric Model Intercomparison Project (AMIP) simulations, starting at 0000 UTC 1 January 1985, with seasonal sea-ice and SST boundary conditions repeating annually, which are applied globally. Because only the atmospheric component of the model is used, the prescribed SSTs and sea ice will not be changed through direct model integration. The atmospheric component of CFSv2 is the Global Forecast System Model (Moorthi et al., 2001), which uses a T126 horizontal grid (∼100 km grid spacing) and finite differencing in the vertical grid with 64 sigma-pressure hybrid layers. Combinations of boundary conditions are used based on sea-ice concentration (SIC) and SST monthly mean 10-year average values for 1981-90 and 2005-14 from the Merged Hadley-NOAA/OI dataset (Hurrell et al., 2008).

The seasonal differences (2005-14 minus 1981-90) of SST and SIC from the Merged Hadley-NOAA/OI data are shown in Figs. 1 and 2, respectively. SSTs show a mostly uniform warming throughout the Northern Hemisphere, except for the west coast of the United States. Mean seasonal differences between the two periods for all points north of 30^°N are 0.29 K for March-April-May (MAM; Fig. 1a), 0.58 K for June-July-August (JJA; Fig. 1b), 0.56 K for September-October-November (SON; Fig. 1c), and 0.34 K for December-January-February (DJF; Fig. 1d). For SIC (Fig. 2) the changes are most pronounced in JJA (Fig. 2b) and SON (Fig. 2c), with large decreases seen over the Chukchi Sea extending along the Siberian coast and the Arctic Ocean. Modest decreases are still seen in MAM (Fig. 2a) and DJF (Fig. 2d) as well.

Figure 1.  SST change (K; top 10 m of ocean) between 1981-90 and 2005-14 from the Merged Hadley-NOAA/OI product: (a) MAM; (b) JJA; (c) SON; (d) DJF.

The average conditions of SST and sea ice for 1981-90 (2005-14) are denoted as SST1 and ICE1 (SST2 and ICE2), respectively. Configurations used for the runs are SST1ICE1, SST2ICE1, SST1ICE2, and SST2ICE2. SST1ICE1 utilizes both SST and SIC monthly means from 1981-90. SST2ICE1 uses SST from 2005-14 and SIC from 1981-90, with the SST changes only having an impact outside the sea-ice domain. Because the model is uncoupled, there is no interplay between the SSTs and sea ice; atmospheric impacts are controlled by fluxes from either the sea ice or ocean surface. Therefore, altering the SSTs will not have an impact on surface-to-atmosphere fluxes in ice-covered regions. SST1ICE2 uses SST from 1981-90 and SIC from 2005-14. In places where sea ice is removed, SST remains at the freezing point of sea water, as in the study by (Liu et al., 2012), with the advantage that the response to sea-ice loss is fully isolated. Other studies have perturbed SSTs to include the SST response that would occur due to sea-ice loss (i.e., Screen, 2014), but either method is considered valid, with each having advantages. Finally, SST2ICE2 incorporates SIC and SST from 2005-14, representing the combined impact of sea-ice loss and SST increase. In SST2ICE2, the SST changes will have an impact on a larger domain than in SST2ICE1, given there is less sea-ice coverage. For this reason, the results are not directly additive (adding SST1ICE2 to SST2ICE1 will not equal SST2ICE2).

Figure 2.  As in Fig. 1, but for SIC (%).

For each of the configurations, the uncoupled model is integrated for 101 years from the atmospheric initial state at 0000 UTC 1 January 1985 from the Climate Forecast System Reanalysis (CFSR; Saha et al., 2010). The first 11 months are excluded in the analysis. Paired differences of these runs are used to analyze the influence of changing SST, SIC, and both parameters. For example, differences between SST2ICE1 and SST1ICE1 are considered to be to the result of SST changes between the two 10-year periods, as that is the only parameter changed between the two simulations; while the differences between SST2ICE2 and SST1ICE1 are the result of impacts from both SST and SIC changes.

Starting with the first December of the 101-year runs, all DJF periods are averaged, creating a total of 100 seasonal means. DJF is chosen as the season of interest due to the large area of observed cooling seen over the Asian continent (Figs. 3a and d), and the goal is to see to what extent this can be duplicated in model simulations. For each configuration, the 100 DJF means are further divided into 10 sets of 10 individual DJF seasons and averaged together for comparisons with the 10-year reanalysis periods (1981-90 and 2005-14). All results from the simulations are based on differences with respect to the SST1ICE1 simulation. A dataset of differences is created using all possible combinations from the 10 sets of each configuration and the 10 sets of SST1ICE1, which equate to 100 possible combinations. The ensemble mean differences between two AMIP simulations highlight the responses from SST and/or SIC forcing, while the differences among individual runs within each simulation represent the contribution from atmospheric internal variability.

Another way to quantify the contribution of atmospheric internal variability is based on the analysis initialized predictions using hindcasts from the coupled CFSv2 model. If the observed anomalies have a contribution from atmospheric internal variability, then the capability of initialized predictions to retain these anomalies would diminish with lead time, with the impact of the changed boundary conditions becoming more apparent at longer leads. Therefore, using this approach, we assess the influence of internal variability in a coupled forecast setting, and on the time scale this influence is retained.

To provide insights into the impacts of initial conditions, the coupled hindcasts with CFSv2 are examined over the two periods that closely match those used in deriving the initial SST and SIC conditions in the AMIP simulations, with the caveat that the CFSv2 hindcasts are not available before 1982. Therefore, the periods compared are DJF 1982-90 and 2005-13. For these runs, the atmospheric model is coupled to the Geophysical Fluid Dynamics Laboratory (GFDL) Modular Ocean Model, version 4 (MOM4) (Griffies et al., 2003), and the GFDL sea-ice simulator. MOM4 contains 40 vertical layers with a resolution of 10 m at the surface that decreases to approximately 500 m at the bottom layer. The sea-ice model contains two equal layers of sea ice and one layer of snow, along with five sea-ice thickness categories (Saha et al., 2010). The model is integrated four times a day (0000 UTC, 0600 UTC, 1200 UTC, and 1800 UTC), every five days.

Mean DJF 2-m temperatures and 200-hPa heights are reconstructed from the hindcast simulations by averaging hindcasts with equal lead days and the four initial times from those days. For example, for the shortest lead time, which is approximately 5 days, the DJF mean for December, January and February would be substantially influenced by the initial conditions on 27 November, 27 December, and 26 January, respectively, and would also average the four initial times on those days for a total of 36 forecasts for each lead within each 9-year period. Leads out to 120 days are analyzed to determine if the observed temperatures and heights can be simulated in the model, and if there is a lead-time dependency. This approach was documented in (Chen et al., 2013). A plus/ minus five-day lead smoothing is also applied; so, as an example, a lead of 10 days would also incorporate the data for 5- and 15-day leads. This is done to reduce the large variations seen among consecutive lead times and produce more coherent results.

Figure 3 shows the mean DJF 2-m temperature and 200-hPa height differences (2005-14 minus 1981-90) for ERA-Interim (Dee et al., 2011) (Fig. 3a) and CFSR (Fig. 3b). The green outlined region in Fig. 3, and subsequent maps, signify the domain of interest for this study and encompass the region of the greatest cooling in ERA-Interim (40^°-60^°N, 30^°-150^°E). Geographically, this includes much of Eastern Europe and Central Asia. Both ERA-Interim and CFSR show pronounced cooling over this domain, not seen in other seasons (not shown). Both reanalyses show substantial warming in the Arctic.

The Eurasian cooling in the reanalyses is not simulated when comparing the means of the AMIP simulations, which all show either warming or no change (Figs. 4a-c) in 2-m DJF temperatures over Eurasia. The area-weighted mean 2-m DJF temperature change between the two periods over this domain in ERA-Interim is -0.27 K (Fig. 3a). For the uncoupled AMIP simulations (Figs. 4a-c), the mean 2-m temperature change is 0.72 K, 0.54 K and 0.32 K for SST2ICE2, SST2ICE1 and SST1ICE2, respectively —— all relative to SST1ICE1. In general, significant warming is confined to the Arctic for SST1ICE2, while the warming is more uniform throughout the Northern Hemisphere in SST2ICE1. The SST2ICE2 difference shows both the strong Arctic warming, and weaker but uniform warming at lower latitudes. Changes in 200-hPa heights are prevalent in ERA-Interim (Fig. 3a) and CFSR (Fig. 3b) and, albeit to a lesser degree, are still present in SST2ICE2 (Fig. 4a) and SST2ICE1 (Fig. 4b), particularly over the Arctic and extending south to Alaska. However, changes in the 200-hPa heights are more sporadic in the SST1ICE2 (Fig. 4c) runs.

Figure 3.  DJF 2-m temperature change (K; shaded) and 200-hPa geopotential height change (dm; contoured at intervals of 2 dm) between 1981-90 and 2005-14 from (a) ERA-Interim and (b) CFSR. The region outlined in green denotes the domain covering Eurasia used in this study.

Figures 4d and e show the zonal-mean 2-m temperature and 200-hPa height differences, respectively. For the reanalysis and the SST2ICE2 and SST2ICE1 simulations, the largest increases in 2-m temperature are in the Arctic, with smaller and more equal changes at lower latitudes. For 200-hPa heights, ERA-Interim shows the largest increase over the Arctic (consistent with Fig. 2a), with SST2ICE2 the next highest, followed by SST2ICE1, and finally SST1ICE2. Changes are also seen over lower latitudes for SST2ICE2 and SST2ICE1, but not for SST1ICE2.

Figure 4.  Mean DJF 2-m temperature (K; shaded) and 200-hPa geopotential height (dm; contoured at intervals of 2 dm) differences determined using all combinations of differences with respect to SST1ICE1 and using 10-year groupings of (a) SST2ICE2, (b) SST2ICE1, (c) and SST1ICE2. All plotted temperature differences are significant at the 95% confidence level. The region outlined in green denotes the Eurasian domain of interest. Zonal-mean 2-m temperature and 200-hPa geopotential height differences based on (a-c) (and Fig. 3a for ERA-Interim) are shown in panels (d, e), respectively.

Figure 5.  Mean DJF zonally averaged temperature (K; shaded) and geopotential height (dm; contours) differences: (a) ERA-Interim (2005-14 minus 1981-90); (b) SST2ICE2 minus SST1ICE1; (c) SST2ICE1 minus SST1ICE1; (d) SST1ICE2 minus SST1ICE1.

Figure 5 illustrates the zonal-mean vertical profiles of DJF temperatures and geopotential heights for ERA-Interim (Fig. 5a) and each of the AMIP configurations (Figs. 5b-d). A comparison of the two periods in ERA-Interim reveals warming throughout the Arctic troposphere, and this is also shown to some degree —— albeit with a weaker magnitude —— in the simulations that impose the SST forcing from 2005-14. Because the warming is vertically more extensive in ERA-Interim than in the AMIP runs, the 200-hPa geopotential height fields are under-simulated in AMIP relative to ERA-Interim, as seen in Fig. 4. The SST1ICE2 simulation only shows warming mainly in the lower Arctic troposphere, which is consistent with previous results (e.g., Kumar et al., 2010; Perlwitz et al., 2015). For Figs. 4 and 5, the results are not exactly additive, meaning that adding the impact of just the SST change (SST2ICE1) to the impact of just the sea-ice change (SST1ICE2) does not equal the result of the simulation where both SSTs and sea ice are changed (SST2ICE2). However, the combined result resembles the SST2ICE2 result (not shown), and small differences are anticipated due to the interplay between SSTs and sea ice, as well as because of issues related to sampling. However, the differences are generally smaller than those seen in the individual runs, and therefore nonlinearity is not expected to significantly impact the results.

To quantify the role of atmospheric internal variability, distributions of differences in the mean 2-m temperature relative to SST1ICE1 are shown in Fig. 6. For Fig. 6a, which takes the differences of SST1ICE1 with itself, a normal distribution is found, as expected, with a mean of 0 K. The standard deviation is 0.46 K. For SST2ICE2, SST2ICE1 and SST1ICE2, the standard deviation is 0.50 K, 0.49 K and 0.50 K, respectively, indicating that the amount of variability is similar in all configurations. The means shift towards warmer temperatures for all configurations (Figs. 6b-d), in accordance with the results shown in Fig. 2, but cooling is still seen in a small percentage of runs. The percentage of negative temperature differences, representing cooling, over Eurasia relative to SST1ICE1, is 45%, 11%, 17% and 24% for SST1ICE1, SST2ICE2, SST2ICE1 and SST1ICE2, respectively.

Figures 7 and 8 consider 2-m temperature (Fig. 7) and 200-hPa height (Fig. 8) patterns associated with the 10 warmest and coolest differences out of the 100 total differences over the domain of interest, with the goal to determine local and remote changes associated with warm and cold winters over Asia and, further, to assess if there is consistency among the different configurations.

For 2-m temperatures, the warmest Eurasian winters in the uncoupled simulations all show increased temperatures throughout the Eurasian domain (Figs. 7a-d). Arctic temperatures are increased in the configurations where sea ice from 2005-14 is used (Figs. 7b and d). The coldest DJF periods all show lower 2-m temperatures in the Eurasian domain (Figs. 7e-h), and have a slightly warmer Arctic than seen in the warm Eurasia differences for each configuration. However, significant warming is only noted near Greenland for all configurations when the warm differences are subtracted from the cold differences (Figs. 7i-l).

Analysis of the 200-hPa height patterns associated with the warmest and coldest DJF periods plotted in Fig. 7 reveals that the coldest DJF periods are associated with increased heights in the Arctic (Figs. 8e-h); specifically, significant increases are seen along the northern coast of Siberia for SST1ICE1 (Fig. 8i) and SST1ICE2 (Fig. 8l), and more expansively throughout the Arctic Ocean when SSTs from 2005-14 are used (SST2ICE2 and SST2ICE1, plotted in Figs. 8j and k, respectively). Due to the presence of these synoptic features in all configurations along with similar variability, it is plausible that the observed cooling over Eurasia is more related to internal atmospheric variability, rather than a change in SST or SIC forcing.

Figure 6.  Distribution of the area-weighted mean temperature change (K) over the Eurasian domain using the 100 combinations of differences relative to the control AMIP simulation: (a) SST1ICE1; (b) SSTICE2; (c) SST2ICE1; (d) SST1ICE2. The shorter vertical bars denote the mean and the horizontal black bars span 1 standard deviation from the mean.

Figure 7.  (a-d) DJF 2-m temperature differences (K) with respect to the control run using the 10 warmest DJF periods over the Eurasian domain based on the distribution in Fig. 5: (a) SST1ICE1; (b) SST2ICE2; (c) SST2ICE1; (d) SST1ICE2. (e-h) Mean of the 10 coldest DJF periods over the Eurasian domain: (e) SST1ICE1; (f) SST2ICE2; (g) SST2ICE1; (h) SST1ICE2. (i-1) Differences at the 95% confidence level between the two panels directly above (coldest minus warmest) for each configuration: (i) SST1ICE1; (j) SST2ICE2; (k) SST2ICE1; (l) SST1ICE2.

Figure 8.  As in Fig. 7, but for 200-hPa geopotential height differences (dm).

To further assess the contribution of atmospheric internal variability and coupled ocean-atmosphere interactions, the initialized predictions based on the coupled model are analyzed next. The background warming of SSTs and decline of sea ice between the two periods is present, with internal year-to-year variability represented within the forecast initial conditions. The coupled predictions show that the Arctic warming seen in ERA-Interim and CFSR is still replicated, but not the Eurasian cooling. Over Eurasia, 2-m temperature changes are small, but are generally positive at short leads and increase at longer leads. Figure 9 illustrates the DJF 2-m temperature difference based on reconstructed seasonal means at varying lead times. For all leads plotted, significant warming is seen at higher latitudes and, although Eurasian cooling is not seen, significant warming is more sporadic. For a 10-day lead, the area of significant warming is limited across the Eurasian domain (Figs. 9a and h), but becomes more prevalent at longer leads.

Figure 9g shows the evolution of the mean 2-m temperature difference over Eurasia as a function of lead time. The trend, based on least-squares fitting, is small (0.04 K per 30 days oflead), with the mean change for a 10-30-day lead being 0.88 K, and increasing to 1.04 K for a 90-110-day lead. However, the fraction of significant warming within the Eurasian domain increases with lead —— specifically, at a rate of 11.1% per 30 days of lead. The mean percentage of the domain that undergoes significant warming for a 10-day lead is 43%, but this increases to 94% at a 110-day lead. There is also a slight decline in spatial correlation with ERA-Interim differences, but consistent Arctic warming allows spatial correlations to remain modest, with a trend of -0.02 per 30 days of lead (Fig. 9i). Similar results are found when the CFSR data are used instead of ERA-Interim (not shown).

The 200-hPa geopotential height differences are more lead-time-dependent than the 2-m temperature differences. The greatest magnitude of significant Arctic height increases is seen at short leads (out to 30 days; Figs. 10a and b), which is similar to the differences between the two periods seen in the ERA-Interim and CFSR reanalysis products. For longer leads, the height increases become smaller (Figs. 10c-f). There is no discernible pattern noted across Eurasia. For the first 10-30-day lead period, the mean height increase over the Arctic (60^°N and polewards, represented by the circular domain in Figs. 10a-f), is 4.47 dm, but only 2.68 dm for the 90-110-day lead period, indicating a sharp downward trend of -0.6 dm per 30 days, which is significant at the 95% confidence level (Fig. 10g). Spatial correlation trends follow a similar pattern, with higher correlations within the first 10-30-day lead (0.49) to no relationship for the 90-110-day lead period (-0.09), resulting in a linear trend in the spatial correlation of -0.19 per 30 days. The downward trend and decreasing spatial correlation with lead time indicates a weakening of the impact of atmospheric initial conditions, implying that cooling, even if initialized, cannot be sustained. As with 2-m temperature, the results are similar for CFSR (not shown).

Figure 9.  DJF reconstructed mean 2-m temperature differences (K) using the CFSv2 coupled hindcasts (2005-13 minus 1982-90) at lead times of (a) 10 days, (b) 30 days, (c) 50 days, (d) 70 days, (e) 90 days, and (f)110 days. All differences in (a-f) are significant at the 95% confidence level and are smoothed with the last and next lead. (g) Area-weighted mean DJF 2-m temperature difference over the Eurasian domain as a function of lag. (h) Percentage of points within the Eurasian domain with a significant increase in temperature at the 95% confidence level. (i) Spatial correlation with ERA-Interim DJF 2-m temperature difference over the same period using all points polewards of 30^°N.

Figure 10.  DJF reconstructed mean 200-hPa geopotential height differences (dm) using the CFSv2 coupled hindcasts (2005-13 minus 1982-90): (a) 10-day lag; 30-day lag; (c) 50-day lag; (d) 70-day lag; (e) 90-day lag; (f) 110-day lag. All differences in (a-f) are significant at the 95% confidence level and are smoothed with the last and next lead. (g) Area-weighted mean DJF 200-hPa geopotential height difference over the Arctic (points 60^°N and polewards) as a function of lag. (h) Spatial correlation with ERA-Interim DJF 200-hPa geopotential height differences over the same period using all points polewards of 30^°N.

Results from the uncoupled and coupled simulations show that the Eurasian cooling in the reanalysis products is not well simulated as a response to changes in SST and sea-ice conditions. When sea-ice forcing is changed from 1981-90 mean conditions to the 2005-14 mean, the significant 2-m temperature changes are confined to the Arctic, with minimal geopotential height changes at 200 hPa and insignificant non-local changes —— a similar result to that of (Perlwitz et al., 2015). These changes in the Arctic do not propagate into the midlatitudes as suggested by other studies (i.e., Liu et al., 2012). For the SST forcing change, there is a significant increase in 2-m temperature throughout much of the Northern Hemisphere, along with a more robust change in the 200-hPa height field than when only the sea-ice forcing is changed. These results are in line with (Screen et al., 2012), showing that local sea-ice loss is the dominant cause of surface Arctic amplification but contributes much less to geopotential height increases at 200 hPa, which are driven by remote SST changes. Based on the results of Figs. 7 and 8, warm and cold Eurasian winters arise due to roughly the same large-scale atmospheric circulation conditions in all model configurations, implying internal variability as the dominant cause rather than changes in ocean surface forcing.

The hypothesized reasons for the inadequacy of models to simulate the Eurasian cooling fall into two main categories: (i) model results suggest that the cooling is a product of internal variability (McCusker et al., 2016) and a mean warming trend is expected to re-emerge (Sillmann et al., 2014; Sun et al., 2016), or (ii) models suffer from biases and deficiencies (Furtado et al., 2015; Handorf et al., 2015) that need to be addressed; for example, those related to stratospheric coupling, which will allow for better representation of surface patterns. Although the conclusion of (Mori et al., 2014) supported a relationship between sea-ice loss over the Barents and Kara seas with Eurasian temperatures, simulations using Coupled Model Intercomparison Project Phase 5 (CMIP5) models have revealed that the frequency of cool winters decreases despite continuing sea-ice decline. This was also supported by (Sun et al., 2016), who suggested that a warm Arctic and cold continents is only a temporary, transient phenomenon, and will become more unlikely under a warming climate.

In our simulations, each configuration of the model can produce cold winters on an individual sample basis, which are shown to be linked to higher Arctic heights resembling the AO pattern, as also shown in (Deser et al., 2004) and (Alexander et al., 2004). It is worth noting that the DJF AO averages more negatively over the 2005-14 period (-0.25) than over the 1981-90 period (0.00), with a greater frequency of stronger negative AO winters (e.g., 2009) in the later period but more strongly positive AO winters (1988 being the largest) in the earlier period. Because the variability of the Eurasian DJF temperatures in each model configuration is nearly the same, it cannot be concluded that the cooling is a response to SST or sea ice. Even in a warming climate, harsh winters will still occur-just less frequently, as highlighted by (Wallace et al., 2014). While the sea-ice forcing from 2005-14 does produce a larger probability of cooler simulations over Eurasia than the changed SST forcing, it is still warmer than the control simulation.

At short lead times, the coupled simulations produce warming and height increases in the Arctic and weaker to sometimes little warming at lower latitudes. At longer leads there is a gradual decrease in the ability of the coupled model simulations to match the reanalysis, especially in terms of 200-hPa heights, which show little correlation beyond a 30-day lead. Meanwhile, 2-m temperatures maintain some resemblance to the reanalysis products throughout the entire period, due to the stable, warmer Arctic and cooler midlatitudes pattern. Interannual variability patterns, such as the AO, Pacific-North America pattern, and the North Atlantic Oscillation, have shown a certain predictive ability at short lead times (Johansson, 2007; Riddle et al., 2013). As noted by (Cohen et al., 2012), due to a model's inability to predict the AO at longer leads, the lack of skill in predicting the 200-hPa height patterns beyond 30 days suggests that the changes at shorter leads are more related to the negative phase of the AO or the internal variability manifested within the initial conditions, rather than sea-ice or SST forcing, or air-sea interactions that exist in the coupled simulations.

The results of this study are valid on a seasonal scale only and do not take into account the frequency or magnitude of cold events within the DJF periods, which may or may not be influenced by the change in SST or sea-ice forcing. Another factor is that this study uses climatological-mean sea ice, as opposed to a single extremely low sea-ice year (e.g., 2012), which could have exerted more robust impacts on the atmosphere. However, because not every year will equate to 2012 levels, the climatological forcing used is more representative of our current climate. With the results of the uncoupled simulations showing that a subset of cooler Eurasian winters is simulated across all model configurations despite a mean warming, and the sharp drop-off in skill in predicting 200-hPa height changes in the coupled simulations after 30 days, the final message is that the observed changes in Eurasian winter 2-m temperatures may be more related to internal variability rather than any specific forcing, but with the caveats mentioned above taken into consideration.