1 Introduction

Changes in surface air temperature (SAT) exert important influences on agriculture, ecosystem, and human lives (Beniston 2004; Feudale and Shukla 2011). Global coupled climate model simulations provide estimates of the climatic responses to anthropogenic and natural radiative forcings (IPCC 2013). However, there are large uncertainties in the simulated SAT changes on regional spatial scales over several decade and longer periods due to the superposition of the external forcing-induced change and internal climate variability (Dai and Bloecker 2019; Deser et al. 2012a, b; Thompson et al. 2015).

Internal climate variability stems from processes within the atmosphere, ocean, land, and cryosphere, and coupled interactions among them (Deser et al. 2014). In most terrestrial regions, a large proportion of internal variability in surface climate change arises from variations in atmospheric circulation (Deser et al. 2014, 2017a, b; Hu et al. 2019; Thompson et al. 2015). Deser et al. (2016) estimated that atmospheric circulation variability accounts for approximately one-third of the observed winter SAT trends during 1963–2012 over North America. Chen et al. (2019) showed that the projected spring SAT trends during 2006–2060 over the mid- and high-latitudes of Eurasia are subject to high uncertainty mainly arising from intrinsic atmospheric circulation variability.

One of the most robust ways to estimate the role of internal variability in simulated climate trends is to employ a large ensemble of simulations with a specific climate model in which each ensemble member is subject to the same external forcing, but is initiated at a slightly different atmospheric condition (Kay et al. 2015; Wettstein and Deser 2014). Since the forcing and model are identical in all ensemble members, the spread in climate trends among the ensemble members derives entirely from the internal variability in the model (Kay et al. 2015; Thompson et al. 2015).

Several studies have found that in addition to the internal atmospheric variability, thermodynamic processes may also contribute to uncertainties in the SAT trends (Chen et al. 2019; Deser et al. 2016; Ding and Wu 2020). Changes in the Arctic sea ice and Eurasian snow have been shown to exert important influences on Eurasian temperature variability (Chen et al. 2016; Chen and Wu 2018; Gao et al. 2014; Wu et al. 2014; Wu and Chen 2016; Ye et al. 2015). Ding and Wu (2020) suggested that, after the internal component of circulation-induced SAT trends has been removed, anomalies in sea ice and snow may contribute to the diversity of residual boreal spring SAT trends over the mid- and high-latitudes of Eurasia among the 30 ensemble members of the Community Earth System Model (CESM). However, Ding and Wu (2020) did not distinguish the cause and effect though there is a good correspondence between spring SAT trends and snow anomalies. While snow anomalies may induce SAT changes, changes in SAT may provide a favorable thermal condition for the formation of snow anomalies. In addition, the contribution of changes in sea ice and snow to SAT anomalies in other seasons has not been addressed. Given the fact that the sea ice and snow coverage displays large seasonal changes (Stroeve et al. 2012; Wu and Chen 2016), an issue worthy of investigation is whether the contribution of sea ice and snow to the internal component of thermodynamically induced SAT trends varies with the season.

This work compares the contributions of sea ice and snow to the uncertainty in multi-decadal trends of SAT during boreal autumn, winter and spring over the mid- and high-latitudes of Eurasia based on a large initial-condition ensemble simulation. Since the snow effect in boreal summer is negligible over most parts of Eurasia, the present analysis excludes summer. The present analysis is conducted for CESM simulations for comparison with Ding and Wu (2020) with an extension of the time period of analysis. The remainder of the paper is structured as follows. Section 2 describes the observational datasets, model simulations, and methods. In Sect. 3, the dynamical adjustment method is applied to remove the atmospheric circulation-induced internal variability from model simulations and observations and analyze the leading patterns of residual SAT trends in the three seasons. Section 4 investigates and compares the roles of sea ice and snow in the formation of residual SAT trends in the three seasons. Section 5 provides the conclusions and discussions.

2 Data and methods

2.1 Data

The observational monthly mean SAT data were derived from the Climatic Research Unit of the University of East Anglia (CRUTS), version 4.04 (Harris et al. 2020), available from 1901 to 2019, on a regular grid of 0.5° latitude by 0.5° longitude. Other SAT datasets, such as the NOAA Merged Land Ocean Global Surface Temperature Analysis (NOAAGlobalTemp), version 5 (Smith et al. 2008; Vose et al. 2012) and the Climatic Research Unit Temperature, version 4 (CRUTEM4; Osborn and Jones 2014), give similar results. For conciseness, only results obtained with CRUTS data are presented in this study. Monthly mean sea level pressure (SLP) data used in this study are derived from the Twentieth Century Reanalysis, version 3 (20CR; Compo et al. 2011) on a 1° latitude × 1° longitude grid from 1901 to 2015. The SLP data are extended to December 2019 by adding monthly SLP anomalies during 2016–2019 from ERA5 (Hersbach et al. 2020) of the European Centre for Medium-Range Weather Forecasts (ECMWF) to the 1901–2015 SLP climatology of 20CR (Lehner et al. 2017). Monthly mean sea ice concentration data are derived from the Hadley Centre Sea Ice and Sea Surface Temperature dataset (HadISST), which are available from 1870 to the present with a horizontal resolution of 1° × 1° (Rayner et al. 2003). The snow cover data are constructed from the Northern Hemisphere 25 km Equal-Area Scalable Earth Grid (EASE-Grid) 2.0 weekly Snow Cover and Sea Ice Extent, version 4 (Brodzik and Armstrong 2013). These EASE-Grid snow cover data are available from October 1966 to the present and are obtained from the National Snow and Ice Data Center (NSIDC). The raw weekly mean snow cover data have been converted to monthly mean on a regular 1° × 1° grid. In the conversion, the weekly snow cover value is first assigned to each day and then the monthly mean value is obtained as average of the daily snow cover values in the specific month. Monthly mean snow water equivalent data are obtained from the ERA5 at 0.25° horizontal resolution for the time period from 1979 to the present. In addition, a parallel analysis is performed to test the robustness of the results using snow cover and snow water equivalent datasets of the NOAA Snow Chart Climate Data Record (CDR; Estilow et al. 2015) and the NASA Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2; Gelaro et al. 2017).

The primary model output employed in this analysis is a 30-member ensemble of the Community Earth System Model, version 1 (CESM1; Hurrell et al. 2013) at a horizontal resolution of approximately 1° latitude and longitude, hereafter referred to as the CESM1 Large Ensemble (CESM-LE). Each of the 30 members of the CESM is forced by the same external forcing: the CMIP5-based (Taylor et al. 2012) historical radiative forcing for 1920–2005 and representative concentration pathway 8.5 (RCP8.5) radiative forcing for 2006–2100. The first ensemble member was initialized from a randomly selected date in a preindustrial (1850) control run of the CESM, and was then integrated forward from 1850 to 2100. Ensemble members 2–30 were all started on 1 January 1920 with slightly different initial conditions taken from ensemble member 1. A random round-off level difference was applied to the air temperature field for each ensemble member. The detail description of the simulations is given in Kay et al. (2015). An 1800-year preindustrial control integration of the CESM (PiCTL; Kay et al. 2015) is also used for the dynamical adjustment procedure. The model variables used in the present analysis include monthly mean SAT, SLP, sea ice concentration (SIC), snow cover fraction (SCF), snow water equivalent (SWE), shortwave radiation (SWR), longwave radiation (LWR), latent heat flux (LHF) and sensible heat flux (SHF).

The present analysis focuses on the period 1979–2019 in each observational dataset and model simulation. This time period was chosen because satellite-based observations of snow and sea ice became available since 1979. Boreal autumn, winter and spring refer to the average of September–November (SON), December–February (DJF), and March–May (MAM), respectively.

2.2 Methods

A dynamical adjustment method based on constructed circulation analogs proposed by Deser et al. (2016) is used to remove the circulation-induced component of changes. Dynamical adjustment technique has been employed to assess the contribution of atmospheric circulation variability to observed SAT trends over the Northern Hemisphere (Deser et al. 2016; Hu et al. 2019; Smoliak et al. 2015; Wallace et al. 2012), land–atmosphere interaction in the central United States (Merrifield et al. 2017), and attribution of the observed hydroclimate trends from the 1980s to the 2010s over the United States Southwest (Lehner et al. 2018). Constructed circulation analogs (CCAs) are employed to dynamically adjust monthly SAT in the CESM-LE and observations. The method is summarized briefly below and the reader is referred to Deser et al. (2016) for a detailed description of the method. Note that there is uncertainty in searching the analogues due to limited samples (Van Den Dool 1994).

For each month and year in the CESM-LE members, a set of closest SLP analogs are selected from the PiCTL in terms of Euclidean distance from the target SLP field over the domain of 20°–90° N and 60° W–180° E. Then, the analogs are randomly subsampled and their optimal linear combinations are computed to reconstruct the target SLP field. The same set of optimal linear weights is applied to the corresponding SAT fields to derive the dynamically induced component of SAT anomalies. This random subsampling and optimal reconstruction procedure is repeated 50 times. Finally, the 50 reconstructed SLP patterns and associated SAT anomalies are averaged to provide a “best estimate” of the dynamical-induced component of SAT.

The forced dynamical component is obtained by averaging the dynamical components in all CESM-LE members. The internal dynamical component in each ensemble member is obtained by subtracting the forced dynamical component from the total dynamical component. Thermodynamic components are obtained as residuals. Removing the internal component of dynamical-induced SAT trends from the total trends leaves the total forced component and internal thermodynamic component, referred to here as dynamically adjusted SAT trends.

In observations, the SLP pattern is constructed from the Twentieth Century Reanalysis (20CR) over the period 1901–2019 (excluding the year in question). Prior to computing the dynamically induced component, we apply a high-pass filter to the SAT data (removing a quadratic trend) to avoid potential influences from anthropogenically forced SAT following Deser et al. (2016). The internal dynamical component of observed SAT anomalies is obtained by conducting a separate dynamical adjustment procedure on the internal component of the observed SLP anomalies, obtained by removing a quadratic trend. Removing the CESM-LE ensemble-mean SLP anomalies from the observed total SLP anomalies as a way of detrending yields slightly different results (figures not shown).

In this study, the statistical significance of trends is evaluated based on a two-sided Student's t test. Regression analysis is used in the present study. The anomalies obtained by regression are estimated using the two-sided Student’s t test as well.

3 Dominant modes of dynamically adjusted SAT trends

The observed SAT trends during 1979–2019 display obvious differences among the three seasons. The SON SAT trends are large positive over Eastern Europe and central Siberia-Russian Far East (Fig. 1a). The former region is within the high pressure system and the latter region lies between low and high pressure systems. The SON SAT trends are relatively small over western Siberia. The DJF SAT trends are positive over western Europe and the coastal regions of central Siberia, but are large negative over central Eurasia (Fig. 1b). A region of large positive DJF SLP trends extends from the Arctic to Eurasia with the maximum magnitude of approximately 6 hPa/41 years (Fig. 1b). The MAM SAT trends are large positive almost everywhere with the largest warming (> 4 °C/41 years) over northern Siberia (Fig. 1c). The corresponding SLP trends feature a large cyclone over the Eurasian continent.

Fig. 1
figure 1

Observed total (a–c) and dynamically adjusted (d–f) SAT (color shading; °C/41 years) and SLP (contours; hPa/41 years) trends during 1979–2019 for boreal autumn (a, d), winter (b, e) and spring (c, f). CESM-LE ensemble-mean SAT and SLP trends for autumn (g), winter (h) and spring (i). The SLP interval is 1 hPa/41 years, with negative values dashed

Full size image

Removing the effects of internal atmospheric circulation variability through CCAs considerably reduces the areas of large warming and cooling and nearly eliminates the SLP trends (Fig. 1d–f). In addition, the dynamical adjustment improves the resemblance between the observed SAT trend patterns and those of externally forced component of the simulated SAT trends, in particular, in MAM (Fig. 1f versus Fig. 1i). The externally forced SAT trends, which are derived as the CESM ensemble-mean trends, tend to show large values (> 1.5 °C/41 years) over the coastal region of the Arctic Ocean and over the region around and to the north of the Caspian Sea (Fig. 1g–i). This warming pattern may be associated with reductions in sea ice and snow, which will be discussed later in Sect. 4. The ensemble-mean SLP trends are nearly zero throughout the Eurasia (Fig. 1g–i).

The simulated SAT trends during 1979–2019 vary widely in pattern and magnitude across the ensemble in all the three seasons (figures not shown). Dynamical adjustment greatly reduces the diversity of SAT trends within the CESM-LE (figures not shown). The remaining spread among the residual SAT trends is largely attributable to internal thermodynamic process as the CCAs method removes more than 90% of the variance of internal SLP trends across the individual ensemble members in all the three seasons (Fig. S1). The method error is small given the length of the control run, the repeated random subsampling and averaging of the optimal linear combination of analogs (Deser et al. 2016).

The spatial distribution of standard deviation of dynamically adjusted SAT trends in the three seasons is examined to understand what thermodynamic processes may contribute to the remaining spread among the residual SAT trends within the CESM-LE. Large variances in the dynamically adjusted SAT trends are observed over the coastal regions of the Arctic Ocean in all the three seasons (Fig. 2). In addition, the region north of the Caspian Sea displays large variance in winter and spring (Fig. 2b, c).

Fig. 2
figure 2

Standard deviation of dynamically adjusted SAT trends (°C/41 years) during 1979–2019 across the CESM-LE for boreal autumn (a), winter (b) and spring (c)

Full size image

The spatial distribution of dynamically adjusted SAT trends is further illustrated by the leading patterns of residual SAT trends. An empirical orthogonal function (EOF) analysis has been performed on the set of 30 dynamically adjusted SAT trend maps over the Eurasian continent (40°–70° N, 0°–140° E) in the three seasons. Note that the EOF analysis is applied in the ensemble domain rather than the traditional time domain and the ensemble-mean residual trends are removed before the EOF analysis (Deser et al. 2014; Hu et al. 2019).

The first EOF mode explains 29.3%, 32.3% and 32.5% of the total variance in autumn, winter and spring SAT trends across the ensemble, respectively. The EOF1 mode of SAT trends during autumn and spring is characterized by same-sign anomalies extending from northern Europe to the Russian Far East and to central Siberia, respectively (Fig. 3a, c), while same-sign anomalies are confined to western Siberia during winter (Fig. 3b). The second EOF mode explains 13.6%, 13.2% and 14.1% of the total variance during autumn, winter and spring, respectively. The spatial distribution of EOF2 mode in the three seasons is very similar and features a southwest-northeast dipole pattern (Fig. 4a–c). In the positive EOF2 phase, positive SAT trend anomalies are located in northern Siberia and negative SAT trend anomalies lie around and north of the Caspian Sea. The corresponding principal components (PCs) exhibit clear member-to-member variations (Figs. 3d–f, 4d–f). These two modes are well separated from each other and from the other modes based on the criterion of North et al. (1982). An analysis of SAT trends is conducted corresponding to both EOF1 and EOF2 modes, and the signal associated with the EOF2 mode is weaker and covers a smaller region. Thus, only results corresponding to the EOF1 mode are presented in the following.

Fig. 3
figure 3

Regressions of dynamically adjusted SAT trends ( °C/41 years) across the ensemble members upon the standardized PC1 corresponding to EOF1 of the residual SAT trends for boreal autumn (a), winter (b) and spring (c). d–f Are the corresponding standardized PC1. Stippling denotes SAT trends significant at the 95% confidence level based on the Student's t test

Full size image
Fig. 4
figure 4

The same as Fig. 3 except for EOF2

Full size image

4 Roles of sea ice and snow

This section investigates the contributions of sea ice and snow to the residual SAT trends. Widespread decline in sea ice extent (Stroeve and Notz 2018) and snow cover (Déry and Brown 2007) has occurred during the recent decades. Previous studies have demonstrated that sea ice plays an important role in the climate variability via the sea ice-albedo feedback and through modulating exchange of energy between the ocean and atmosphere (Gao et al. 2014; Honda et al. 2009; Liu et al. 2012; Wu and Wang 2020). Changes in snow may exert significant influences on regional and global climate via the snow-albedo feedback and snow-hydrological effect (Barnett et al. 1989; Thackeray et al. 2019; Yasunari et al. 1991). For the sea ice-albedo or snow-albedo effect, the change in sea ice or snow cover influences surface temperature through reflection of solar radiation back to the atmosphere due to high albedo. The snow-hydrological effect involves the demanding of heat for snow melting and moistening soil when it melts. Previous studies have demonstrated the impacts of the Arctic sea ice and Eurasian snow on the Eurasian temperature variability (Gao et al. 2014; Petoukhov and Semenov 2010; Wu et al. 2014; Ye et al. 2015). The effects of sea ice and snow tend to be large along the periphery of sea ice and snow covered region where the surface albedo change is large (Wu and Chen 2016; Wu and Wang 2020).

Here, climatological features of the Arctic sea ice and Eurasian snow cover are examined to provide a background for understanding the relationship between sea ice and snow changes and SAT anomalies. Figures 5a–c and 6a–c present climatological mean (contours) and standard deviation (shading) of SIC and SCF, respectively, in the three seasons based on the CESM-LE ensemble mean. Arctic sea ice and Eurasian snow display substantial seasonal variations in spatial extent. The sea ice periphery (areas with an average SIC < 70%) shows a large southward advance from SON to DJF (Fig. 5a, b). The largest standard deviation of SIC occurs along the sea ice periphery (Fig. 5a–c). The snow covered region moves southwestward from SON to DJF and retreats northeastward from DJF to MAM (Fig. 6a–c). The region of large standard deviation of SCF displays a seasonal move following the southern boundary of snow covered region (Fig. 6a–c). The climatological mean and standard deviation of observed SIC (Fig. 5d–f) and SCF (Fig. 6d–f) display similar seasonal move though the location of the periphery displays some differences between the model simulations and observations.

Fig. 5
figure 5

The CESM-LE ensemble-mean (a–c) and observed (d–f) climatology (contours) and standard deviation (shading) of year-to-year variation of SIC for boreal autumn (a, d), winter (b, e) and spring (c, f) based on the period 1979–2019. Unit is %

Full size image
Fig. 6
figure 6

The same as Fig. 5 except for SCF (%)

Full size image

In the following, the SIC trend anomalies corresponding to the leading EOF mode of dynamically adjusted SAT trends and associated surface heat flux anomalies are investigated to illustrate the effect of sea ice. Then, the SCF and SWE trend anomalies corresponding to the leading EOF mode of dynamically adjusted SAT trends are analyzed. After that, an analysis of surface heat flux anomalies is conducted to explore whether snow anomalies contribute to the SAT trends. Note that the dynamical adjustment removes the signals in SIC, SCF and SWE trends that are related to Eurasian SAT trends through atmospheric circulation changes.

Notable SIC trend anomalies are detected in the Barents-Kara Sea in the three seasons with differences in the location and magnitude of anomalies. Figure 7 shows the SIC trend anomalies associated with EOF1 of the residual SAT trends obtained by regression upon the PC1. Significant negative SIC anomalies are observed over the Barents-Kara Seas during autumn and spring (Fig. 7a, c), consistent with significant positive SAT trends in the above regions (Fig. 7g, i). The EOF1-related SIC trends during winter also show a decrease over the Barents-Kara Seas but less significant (Fig. 7b), corresponding to less significant SAT trends with smaller areal coverage (Fig. 7h). The CCAs have been applied to dynamically adjust the simulated SIC anomalies in the same way as for SAT. The spatial distribution of the thermodynamic component of SIC trends (Fig. 7d–f) resembles closely that of total SIC trends in all the three seasons (Fig. 7a–c). Note that the large signal of SIC trend anomalies tends to occur in the region of climatological mean sea ice periphery and large standard deviation in the model (Fig. 5a–c). The seasonal change in the location of large SIC trend anomalies also follows that of climatological mean and large standard deviation of SIC.

Fig. 7
figure 7

Regressions of total (a–c) and thermodynamic (d–f) SIC trends (%/41 years) and SAT trends (g–i °C/41 years) upon the normalized PC1 of residual SAT trends in the CESM-LE for boreal autumn (a, d, g), winter (b, e, h) and spring (c, f, i). Stippling denotes anomalies reaching the 95% confidence level

Full size image

The forced dynamical component of SIC trends during autumn, winter and spring from the CESM-LE has negative values over the Barents-Kara Seas, but much smaller compared to the total forced SIC trends (Fig. S2). However, the dynamical component may contribute a substantial part to the total SIC trends in individual realizations during all three seasons (figures not shown). The observed SIC trends are negative over most of the Barents-Kara Seas with small contribution from dynamics in boreal autumn, winter and spring (Fig. S3).

Surface heat flux trend anomalies are analyzed corresponding to EOF1 of the residual SAT trends to illustrate the effect of sea ice on the SAT trends. As noted above, the sea ice effects include modulation of surface albedo and heat conduction. The former affects upward SWR and LWR. The latter affects SHF and LHF. Less sea ice may lead to less upward SWR so that more SWR is absorbed by ocean surface, which raises ocean surface temperature. In turn, this results in more upward LWR, SHF and LHF, which contributes to SAT increase. Figure 8 displays anomalies of upward SWR, upward LWR, SHF and LHF obtained by regression upon the normalized PC1 of residual SAT trends. In all the three seasons, upward SWR is reduced while upward LWR, SHF and LHF increase in the Barents-Kara Sea. This confirms the effect of sea ice trends. In comparison, negative upward SWR anomalies are largest in spring and smaller in winter and located at higher latitudes in autumn. The seasonal difference is partly related to the latitudinal and seasonal change of incoming SWR and the change in the location of sea ice reduction region. The smallest upward SWR anomalies in winter are attributed to small incoming SWR and small sea ice trends. The positive upward LWR anomalies are larger and cover a broader area in autumn and spring and are smallest with a small areal extent in winter due to the small sea ice trends. The positive SHF and LHF anomalies are located at highest latitudes in autumn, which is due to the effect of surface wind speed increase from coastal region to open ocean (Figure not shown).

Fig. 8
figure 8

Regressions of surface heat flux (upward: positive) for boreal autumn (a, d, g, j), winter (b, e, h, k) and spring (c, f, i, l) upon the normalized PC1 of dynamically adjusted SAT trends in the CESM-LE. SWR (a–c), LWR (d–f), SHF (g–i) and LHF (j–l) anomalies are all measured in W/m2/41 years. Stippling denotes anomalies reaching the 95% confidence level

Full size image

SAT anomalies over the high latitudes of Eurasia are further examined based on regression upon the SIC trends averaged over the Barents-Kara Seas (70°–80° N, 30°–80° E). Pronounced negative SAT trends are observed over the northern coast of Eurasia in all the three seasons (Fig. 9a–c). In comparison, SAT anomalies are more confined to the coastline during winter than during autumn and spring. The above relationship is further illustrated by Fig. 9d–f that display scatter plots of the internal thermodynamic component of SAT trends averaged over northern Siberia (65°-80°N, 30°-120°E) against the internally-generated SIC trends averaged over the Barents-Kara Seas (70°-80°N, 30°-80°E). It is clear that the SAT trends over northern Siberia have a significant negative relation with the SIC trends over the Barents-Kara Seas across the 30 ensemble members during autumn, winter and spring (Fig. 9d–f). The correlation coefficient in autumn, winter and spring is − 0.74, − 0.40 and − 0.59, respectively, all of which are significant at the 99% confidence level according to Student's t test. The results indicate that the internally-generated SIC variability over the Barents-Kara Seas contributes partly to the uncertainties of SAT trends over northern Siberia.

Fig. 9
figure 9

Regression of dynamically adjusted SAT trends (°C/41 years) onto the normalized SIC trends averaged over the Barents-Kara Seas (70°–80°N, 30°–80°E) for boreal autumn (a), winter (b) and spring (c). Stippling denotes anomalies reaching the 95% confidence level. Scatterplots of the internal thermodynamic component of SAT trends (°C/41 years) averaged over northern Siberia (65°–80°N, 30°–120°E) against the internally-generated SIC trends (%/41 years) averaged over the Barents-Kara Seas (70°–80°N, 30°–80°E) in CESM 30 members for boreal autumn (d), winter (e) and spring (f)

Full size image

Obvious signals are detected in snow variations over the mid- and high-latitudes of Eurasia in the three seasons with notable seasonal variation in the location and magnitude. Figure 10 displays SCF anomalies associated with EOF1 of the dynamically adjusted SAT trends in the three seasons obtained by regression on the PC1. Negative SCF anomalies are located at higher latitudes and have a larger areal coverage but with a smaller magnitude during autumn (Fig. 10a), confined to lower latitudes with a larger magnitude during winter (Fig. 10b), and located over Eastern Europe and western Siberia during spring (Fig. 10c). The spatial distribution of the thermodynamic component of SCF trends (Fig. 10d–f) obtained using the CCAs resembles closely that of total SCF trends in all the three seasons (Fig. 10a–c). The distribution of SWE trend anomalies (Fig. S4) bears resemblance to that of SCF anomalies (Fig. 10). In comparison, negative SWE anomalies cover smaller areas than negative SCF anomalies. Significant negative SCF anomalies over Eastern Europe and west Siberia (Fig. 10a–c) correspond to positive SAT trends (Fig. 3a–c). The large signal of SCF trend anomalies tends to coincide with large standard deviation of SCF located at the southern boundary of climatological mean snow covered region in the model (Fig. 6a–c). The seasonal change in the location and magnitude of large SCF trend anomalies also follows that of large standard deviation of SCF.

Fig. 10
figure 10

Regressions of total (a–c) and thermodynamic (d–f) SCF trends (%/41 years) upon the normalized PC1 of dynamically adjusted SAT trends in the CESM-LE for boreal autumn (a, d), winter (b, e) and spring (c, f). Stippling denotes anomalies reaching the 95% confidence level

Full size image

The forced dynamical component of simulated SCF trends is small compared to the CESM-LE ensemble-mean SCF trends over most of the mid- and high-latitude Eurasia during all the three seasons (Fig. S5). The fraction of the dynamical component of SCF trends during autumn, winter and spring shows large values in individual members with scattered features (Figures not shown). The dynamical contribution of SCF and SWE trends is small based on the point-wise partial least-square regression (PLSR) used in Smoliak et al. (2015) over most of Eurasia (Figs. S6 and S7). Consistent results are obtained using CDR and MERRA-2 datasets (figures not shown).

The correspondence between snow and SAT anomalies suggests a connection between snow changes and dynamically adjusted SAT trends. Surface heat fluxes changes are analyzed to illustrate the cause and effect in this connection. As noted earlier, changes in snow may affect surface heat fluxes in two ways. One is modulation of upward SWR through alteration in surface albedo. The other is turbulent heat exchange between the land and air. Less snow cover is followed by less upward SWR and thus more gain of net SWR by land surface, which leads to increase in land surface temperature. This induces increase in upward LWR, resulting in SAT increase. During snow accumulation seasons (boreal autumn and winter in the Northern Hemisphere), less snow amount increases SHF and upward LWR from land surface to the air, leading to an increase in SAT. However, in deep snow region, this effect is weak as heat conductivity of snow is low. During snow melting seasons (boreal spring in the Northern Hemisphere), less snow amount increases the part of heat consumed for warming the land surface, followed by more SHF and upward LWR, and thus resulting in an increase in SAT. In addition, less snow reduces the soil moisture when snow melts, leading to a decrease in upward LHF. It should be noted that changes in SAT and air humidity induced by atmosphere processes may contribute to variations in SHF and LHF.

The effects of snow on SAT trends are illustrated by surface heat flux anomalies corresponding to EOF1 of the dynamically adjusted SAT trends. Figure 11 displays anomalies of upward SWR, upward LWR, SHF and LHF obtained by regression upon the normalized PC1 of residual SAT trends. Corresponding to EOF1, autumn upward SWR decreases over the northeastern European Plain and west Siberia (Fig. 11a) where SCF is reduced (Fig. 10a). During winter, a significant decrease in upward SWR around and north of the Caspian Sea corresponds to reduced SCF (Figs. 10b, 11b). Significant decreases in spring upward SWR are observed from the eastern European Plain to northern Siberia (Fig. 11c), and are related to a decrease in SCF (Fig. 10c). The above correspondence between upward SWR and SCF anomalies confirms the effect of snow albedo change. The spatial pattern of upward LWR anomalies (Fig. 11d–f) bears close resemblance to that of SAT trends during all the three seasons (Fig. 3a–c). Over most regions of the Eurasian continent, the change in SHF is weak and insignificant in all the three seasons (Fig. 11g–i). A significant increase in LHF is observed over central Siberia during autumn (Fig. 11j) and to the north of the Caspian Sea during winter and spring (Fig. 11k, l). These LHF changes cannot be explained by SCF decrease (Fig. 10a–c). The above results verify the contributions of snow changes to SAT trends through modulating upward SWR and upward LWR.

Fig. 11
figure 11

Regressions of surface heat flux (upward: positive) for boreal autumn (a, d, g, j), winter (b, e, h, k) and spring (c, f, i, l) upon the normalized PC1 of dynamically adjusted SAT trends in the CESM-LE. SWR (a–c), LWR (d–f), SHF (g–i) and LHF (j–l) anomalies are all measured in W/m2/41 years. Stippling denotes anomalies reaching the 95% confidence level

Full size image

5 Conclusions and discussions

The present study investigates the contributions of sea ice and snow to the multi-decadal trends of SAT during autumn, winter and spring over the mid- and high-latitudes of Eurasia based on a 30-member ensemble of simulations by CESM. A dynamical adjustment technique is used to remove the internal component of atmospheric circulation-induced SAT trends. During all the three seasons, dynamical adjustment greatly reduces the diversity of SAT trends among ensemble members and brings observed trends closer to the ensemble mean. The remaining differences in dynamically adjusted SAT trends are largely attributed to internal thermodynamic processes.

An EOF analysis is applied to obtain the dominant patterns of residual SAT trends. The first mode features a broad pattern of the same sign loading that extends eastward from northern Europe to the Russian Far East in autumn and to central Siberia in spring, respectively, whereas the same-sign loading is confined to western Siberia in winter. The second modes are very similar in the three seasons and display a dipole structure with same-sign loading over northern Siberia and opposite-sign loading around and north of the Caspian Sea.

Analysis shows that the internal component of SIC trends over the Barents-Kara Seas feeds back onto SAT locally and over adjacent northern Siberia via thermodynamic processes, contributing partly to the spread in the thermodynamic component of internal SAT trends within the CESM-LE. The SIC-induced SAT trends extend inland in autumn and spring and are confined to the coastal regions in winter. The SIC effects and their seasonal changes are confirmed by upward SWR and upward LWR anomalies.

Eurasian snow changes contribute to the spread in dynamically adjusted SAT trends around the periphery of snow covered region where the SCF anomalies are large mainly by modulating surface albedo. The region of snow effects displays a change with the season following the change of snow anomalies and associated upward SWR and upward LWR anomalies. The snow effect is identified over northeastern European Plain-western Siberia in autumn, over the region north of the Caspian Sea in winter, and over eastern Europe-northern Siberia in spring.

The present study assesses the roles of dynamics and thermodynamics in the SAT trends separately. However, pronounced feedbacks may exist between dynamic and thermodynamic effects (Fischer et al. 2007; Hoerling et al. 2014). Outputs from large ensemble simulations with models other than CESM need to be examined to address the above issue in the future.