Is Eurasian snow cover in October a reliable statistical predictor for the wintertime climate on the Iberian Peninsula?

In this study, the recently found lead–lag relationship between Eurasian snow cover increase in October and wintertime precipitation totals on the Iberian Peninsula is re‐visited and generalized to a broad range of atmospheric variables on the synoptic and local scale. To this aim, a robust (resistant to outliers) method for calculating the index value for Eurasian snow cover increase in October is proposed. This ‘Robust Snow Advance Index’ (RSAI) is positively correlated with the wintertime (DJF) frequency of cyclonic and westerly flow circulation types over the Iberian Peninsula, while the corresponding relationship with anticyclonic and easterly flow types is negative. For both cases, an explained variance of approximately 60% indicates a strong and highly significant statistical link on the synoptic scale.


Introduction
Wintertime precipitation totals on the Iberian Peninsula were recently found to be statistically related to Eurasian snow cover increase during the previous October (Brands et al., 2012). A possible dynamical pathway for this formerly unknown lead-lag relationship has been identified by observational and idealized generalized circulation model studies, linking Eurasian snow cover in fall to the Northern Hemisphere extratropical circulation during the following winter (Cohen and Entekhabi, 1999;Saito et al., 2001;Gong et al., 2003;Cohen et al., 2007;Fletcher et al., 2007Fletcher et al., , 2009Smith et al., 2011;Mote and Kutney, 2008;Peings et al., 2012;Orsolini et al., 2013), the latter commonly described by the Arctic Oscillation (Kutzbach, 1970;Thompson and Wallace, 1998). Following the conceptual model by Cohen et al. (2007), a positive snow-cover anomaly in October leads to the early appearance of a strong Siberian cold high, to large amplitudes in the Rossby-wave train and to an upward wave activity flux that weakens the stratospheric polar vortex (Polvani and Waugh, 2004). Owing to the relatively long de-correlation time of the latter (Baldwin et al., 2003), this weakening/warming persists for several months and propagates downward to the troposphere (Baldwin and Dunkerton, 1999), favouring a negative tropospheric AO during the following winter months. As the North Atlantic Oscillation (NAO) (Walker and Bliss, 1932) can be interpreted as a regional manifestation of the AO (Thompson and Wallace, 1998), this remote snow forcing is also expected to favour anomalous climate conditions over the North Atlantic and adjacent land areas, such as the Iberian Peninsula (Zorita et al., 1992;Hurrell, 1995;Rodo et al. 1997;Rodriguez-Puebla et al., 2001;Goodess and Jones, 2002;Lorenzo et al., 2008).
This study is dedicated to the regional manifestation of this hemispheric-wide teleconnection for the Decemberto-February (DJF) mean climate on the Iberian Peninsula. To this aim, a robust method for calculating the index value of October Eurasian Snow cover increase (Cohen and Jones, 2011) is proposed. This 'Robust Snow Advance Index' (RSAI) is shown to be significantly associated with the DJF circulation characteristics over the Iberian Peninsula, which in turn control the concurrent

Data and methods
Two types of predictand data covering the Iberian Peninsula are used in this study: (1) large-scale circulation data for calculating weather type frequencies and (2) in situ station data.
The large scale circulation is represented by daily instantaneous mean sea-level pressure (MSLP) fields at 12 UTC from the ERA-Interim reanalysis (Dee et al., 2011), which were downloaded from ECMWF's public server (http://data-portal.ecmwf.int/ data/d/interim_daily/). In situ station data were provided by the European Climate Assessment and Dataset Project (ECA&D, http://eca.knmi.nl/dailydata/ predefinedseries.php) (Tank et al., 2002;Klok and Tank, 2009), documenting daily precipitation amount (in mm), sun hours, cloud cover (in octas), wind speed (daily mean value in m/s) and diurnal temperature range (DTR), the latter obtained by subtracting the daily minimum from the daily maximum temperature and assuming a missing value in case this difference is negative. A time series is excluded if the percentage of missing values exceeds the 5% threshold. Finally, DJF averages are calculated upon the daily values, not taking into account the 29th of February in leap years.
Following the study of Cohen and Jones (2011), snow cover increase over mid-latitudinal Eurasia (25-60 • N and 0-180 • E) is calculated for each October between 1997 and 2011 (n = 15), using daily satellite retrievals from the Interactive Multisensor Snow and Ice Mapping System (Ramsay, 1998) obtained at ftp://sidads.colorado. edu/pub/DATASETS/NOAA/G02156/24km/. For a given October, the snow cover extension over the abovementioned spatial domain is calculated for each of the 31 days, yielding a sample of 31 km 2 values. The index value describing snow cover advance is then defined as the regression coefficient (i.e. the slope) of the linear regression line fitted to this sample. A visual inspection of the daily snow cover time series revealed the presence of large and rapid snow cover increases, which are especially prominent in October 2011 (see last two days in Figure 1(b)). As this snow cover 'surges' are outliers from a statistical point of view, and because ordinary least-squares regression is known to be sensitive to outliers, a robust linear regression method for calculating snow cover increase is proposed as an improvement of the original definition of the 'Snow Advance Index' (SAI) (Cohen and Jones, 2011). This method gives less weight to outlier data points when fitting the regression line and essentially removes outlier-related uncertainty (Street et al., 1988). The slope/regression coefficient obtained from robust linear regression is then defined as the RSAI and the 15 index values for October 1997-2011 are standardized to have zero mean and unit variance. As shown in Figure 1 for the case of October 2011 (panel b) as compared to October 2009 (panel a), the modified index differs considerably from the original one if outliers are present in the underlying data. Figure 1(c) shows the comparison between both indices for the 15 October months between 1997 (when satellite-sensed snow cover data on daily timescale became available) and 2011, exhibiting large differences for October 2011. For a detailed description of the applied robust linear regression method, see Appendix.
To compute discrete weather types from continuous daily MSLP patterns, the automated Lamb weather typing (LWT) approach is applied (Jenkinson and Collison, 1977;Jones et al., 1993). The LWT-specifications described by Trigo and DaCamara (2000) have been adopted, using daily MSLP data at 12 UTC from ERA-Interim covering all DJF-days between 1997DJF-days between /1998DJF-days between and 2011DJF-days between /2012. The reanalysis data are linearly interpolated to the 16 grid-boxes shown in Figure 3(a), forming a 'cross' centred over the Iberian Peninsula. Following the study of Trigo and DaCamara (2000), we opt for classifying all days, that is, work with 26 classes instead of the original 27 classes. Composite maps showing the temporal mean of the MSLP values corresponding to a given weather type (i.e. the conditional mean) have been calculated to assure that the LWTs are physically meaningful. These composite maps are similar to those obtained by Trigo and DaCamara (2000) and the 14 WTs considered in this study are displayed in Figure 3.
In contrast to other automated classification techniques like 'Self Organizing Maps' (Hewitson and Crane, 2002;Gutierrez et al., 2005) or 'k -means clustering' (Gutierrez et al., 2004), LWT is a rule-based classification scheme where the classes are pre-defined based on meteorological expert knowledge (Lamb, 1972). This is convenient for the present type of study, as applying stochastic classification schemes would lead to slight differences in the obtained frequencies of weather types, caused by the fact that some days would be assigned to different classes in different realizations. This, in turn, would inhibit a proper estimate of statistical significance when correlating weather type frequencies against another variable (Hewitson and Crane, 2002), that is the RSAI in this study.
To reveal the statistical relationship between the RSAI and the target variables on the Iberian Peninsula, the Pearson correlation coefficient (hereafter 'Pearson correlation' or 'r') is used. Owing to the short sample size (n = 15), which is imposed by the availability of daily snow cover data (Ramsay, 1998), outlier-presence can falsify the results of the correlation analysis. For instance, the DJFseason 2009/2010 was characterized by an extremely negative phase of the AO and NAO (Cohen et al., 2010), associated with largely anomalous values for the concurrent climate conditions on the Iberian Peninsula (Vicente-Serrano et al., 2011). To check the robustness of the results to this 'outlier-winter', the Spearman rank correlation coefficient (hereafter 'Spearman correlation' or 'rs') is used in addition (Wilks, 2006). To test if the DJF-mean target variables can be hindcast from the October RSAI, simple linear regression is applied in a one-year-out cross-validation framework (Michaelsen, 1987), using the RSAI as the only predictor variable. Note that the predictor-predictand pairs withheld in each step of the cross-validation are not truly independent as all predictor-predictand pairs (i.e. the whole available time series) have been used for searching the statistical relationship between October Eurasian snow cover and the DJF-mean climate on the Iberian Peninsula (DelSole and Shukla, 2009). Consequently, our results obtained from cross-validation (see Section 4) might suffer from artificial skill, that is, might not reflect the skill the statistical 'forecasting' method would obtain in real/future forecast situations (see also Section 7). Therefore, when referring to statistical 'predictions' obtained from cross-validation, we will use the term 'hindcast' instead of 'forecasts'.
Note also that using robust instead of ordinary regression as a statistical forecasting method leads to virtually identical results. Therefore, the results obtained from the simpler method (i.e. ordinary regression) will be shown in Section 4. To exclude the possibility that the results of the one-year-out cross-validation could be biased by linear trends, the predictor/predictand samples used to obtain the forecast equation (sample size: n − 1) are linearly de-trended and centred to have zero mean. To eliminate a further potential source of dependency/artificial skill, the trend and mean removal is repeated in each step of the cross-validation (von Storch and Zwiers, 1999), that is, the i th forecast is obtained from predictor/predictand samples having zero linear trend and zero mean in any case. Note that the trend and mean values obtained in the i th step of the cross-validation are also removed from the i th withheld predictor and predictand values, respectively.
To assess the skill obtained from cross-validation, the hindcast time series are compared to their observed counterparts by using the Pearson correlation coefficient. The corresponding results will be referred to as 'hindcast correlations' (Folland et al., 2012) hereafter. Local statistical significance is assessed with a two-sided t-test (H 0 = zero correlation) and global significance is tested by repeating the cross-validation procedure 2000 times. In each repetition, the RSAI time series is randomly re-ordered following the ranks of a random sample of integers drawn from the standard uniform distribution. Then, the percentage of significant local hindcast skill (α local = 0.05) obtained by chance is calculated and saved. The 99th percentile of the resulting sample of 2000 percentage values is the critical value above which global significance (α global ) is assumed at a test-level of 1%.
Apart from forecasting with the October RSAI, the one-year-out cross-validation approach is additionally applied for the October SAI, that is, for the original Snow Advance Index described by Cohen and Jones (2011), as well as for the October monthly mean NAO and AO indices (both based on Principal Component Analysis) as defined by Hurrell et al. (2003) (http://climatedataguide. ucar.edu/es/guidance/hurrell-north-atlantic-oscillationnao-index-pc-based) and the Climate Prediction Center (www.cpc.ncep.noaa.gov/products/precip/CWlink/daily_ ao_index/ao_index.html), respectively.
In addition to the Pearson correlation, the root mean squared error skill score (rmsess) is applied for assessing the out-of-sample skill (Jolliffe and Stephenson, 2003): where rmse is the root-mean squared error of the time series predicted by the statistical forecasting method described above and rmse ref is the rmse obtained from always predicting the climatological mean, which is zero in any case since anomalies are forecasted. Thus, rmsess gives the percentage with which the purely climatological forecast is outperformed by predicting from October Eurasian snow cover increase.

Relevance of serial correlation
There are two potential reasons why positive serial correlation could adversely affect the results of this study. First, and if not accounted for (see Equation (2)), positive serial correlation leads to too-many type 1 errors because of an artificial lowering of the correlation coefficient's p-value, arising from the fact that the number of temporally independent data pairs is lower than the sample size (Trenberth, 1984;Kristjansson et al., 2002). Therefore, the p-value is calculated upon the effective sample size (n*): where n is the sample size and l 1 (l 2 ) is the lag1 autocorrelation coefficient of time series 1 (2) (Bretherton et al., 1999;Beranova and Huth, 2007). Note that the time series are assumed to follow a first-order autoregressive process and that the effective number of degrees of freedom for the two-sided t-test is n* − 2. Second, positive serial correlation questions the applicability of the one-year-out cross-validation approach which assumes zero serial correlation for the predictor/predictand time series (Michaelsen, 1987). To assess the degree to which our forecast skill estimates are affected by serial correlation, Figure 2 displays the spatial distribution of the lag1 autocorrelation coefficients (r − lag1) obtained from the 61-66 samples of each predictand variable (note that the number of station time series slightly varies from one variable to another). As can be seen from figure, the median (bar of the boxplot) ranges between ±0.1 and the interquartile range (box of the boxplot) between ±0.2 for any of the applied variables, that is, the samples are approximately centred around 0. Owing to the limited sample size, critical values for significantly positive r − lag1 (note that the t-test should be one-sided, as only positive values of r − lag1 decrease the effective sample size) would be high even for a local test level of 10%. Therefore, r − lag1 > +0.25 was defined as an alternative threshold above which serial correlation would have a measurable effect on cross-validation estimates of forecast skill, as was originally proposed by Michaelsen (1987). The percentage of time series exceeding this threshold (which will hereafter referred to as 'problematic') can be obtained from Table 1 for each predictand variable under study. For all predictand variables except wind speed, at most 5% of the applied time series suffer from a problematic serial correlation.
To additionally assess if the areal fraction of problematic serial correlations is globally significant, the following Monte-Carlo technique is applied separately for each of the five predictand variables. First, the temporal sequence of all time series corresponding to a given predictand variable is randomly re-ordered following the ranks of a random sample of integers drawn from the standard uniform distribution. As this re-ordering is identical for all time series of a given predictand variable, the spatial autocorrelation of the field is maintained whereas the serial correlation is eliminated. Then, the areal fraction for which r − lag1 > +0.25 by chance is calculated and saved. After repeating the whole procedure 2000 times, the 90th (95th) percentile of the resulting 2000 areal fractions for which r − lag1 > +0.25 by chance is assumed as the critical value above which the fraction of local r − lag1 > +0.25 obtained from the correct time  series (i.e. having the correct temporal order) is globally significant at a test-level of 10% (5%). Even for a global test-level of 10% (α global = 0.10), in which case the global H 0 ('the observed fraction of r − lag1 > +0.25 arises from chance') is easier to reject than for assuming α global = 0.05, the H 0 cannot be rejected for any single predictand variable (see Table 1). On the basis of these results, we conclude that the hindcast skill estimates obtained from one-year-out crossvalidation (see Section 4) are generally not seriously affected by serial correlation. Figure 3 shows the composite maps of the 14 weather types relevant for this study. For the ease of understanding, the panels of the figure are ordered to follow the cardinal directions, that is, westerly flow types are displayed on the left hand side and easterly ones are shown on the right hand side. 'CNW' is the acronym for 'cyclonic northwest', 'ANE' for 'anticyclonic northeast', etc. Above each panel, two numbers are displayed. The first refers to the number of the weather type and the second to its relative frequency (in %) over the whole period under study (DJF 1997(DJF /1998(DJF to 2011(DJF /2012. For an adequate visualization of the direction of the geostrophic flow (hereafter: 'flow'), 10 isobars are displayed in each panel. The pressure gradient can be derived from the colour shading.

Results
With a relative frequency of 55.8%, the wintertime circulation over the Iberian Peninsula is dominated by anticyclonic and easterly flow conditions, while cyclonic and westerly flow conditions occur on 22.5% of the days. The spatial patterns of the rare hybrid weather types (4, 5, 6 and 9, 10, 11) are similar to their more frequently occurring purely directional flow counterparts (1, 2, 3 and 12, 13, 14).
In order to obtain an adequate sample size for each of the 15 winter seasons, and similar to the approach applied by Trigo and DaCamara (2000), the DJFfrequencies/counts of the cyclonic and westerly flow types on the one hand and those of the anticyclonic and easterly flow types on the other are aggregated to two groups with opposite vorticity and flow characteristics. The corresponding two composite maps are shown in Figure 4(a) and (b). With a standard deviation of 12 and 14 d from a mean of 20 and 50 d respectively, the yearly DJF-counts of both groups are characterized by a large inter-annual variability, especially in case of the cyclonic and westerly flow group.
The second row of Figure 4 indicates that this interannual variability is statistically related to Eurasian snow cover increase in October as represented by the RSAI, yielding a fraction of explained variance of approximately 60% for both of the above mentioned groups. With a Pearson/Spearman correlation of +0.76/+0.83, this relationship is positive for the frequency of cyclonic and westerly flow types (see Figure 4  inverse relationship of −0.80/−0.86 is found for the anticyclonic and easterly flow types (see Figure 4(d)). Note that the linear trend of the frequency time series has been removed before computing these correlations and that the anomalies of the de-trended time series are displayed in Figure 4(c) and (d). The RSAI is displayed in standardized anomalies. As both the Pearson and Spearman correlations are significant at a test level of 1%, and as the results for the non-detrended time series (see parentheses in the headers of Figure 4(c) and (d)) are similar, the strong statistical relationships are (1) very unlikely to be a product of chance, (2) insensitive to outlier values and (3) insensitive to possible trends in the underlying data.
To link synoptic and local scale predictability, Figure 5 shows the mean values, conditioned to the two defined groups, for precipitation amount, DTR, sun hours, cloud cover and wind speed at the available weather stations. In the left column, the mean of the local values pertaining to cyclonic and westerly flow days is displayed at each site, while in the right column the corresponding results for the anticyclonic and easterly flow days are shown. The mean surface climate associated to cyclonic and westerly flow types is generally characterized by wetter conditions, a reduced DTR, less sun hours, cloudier skies and windier conditions than it is the case for the anticyclonic and easterly flow types.
As it has been shown that Eurasian snow cover increase is a significant predictor for the DJF-circulation dynamics over the Iberian Peninsula, which in turn control the concurrent mean conditions of various climate variables on the local scale, the latter will be directly hindcast from the October RSAI in the next step of the study. Figure 6 shows the hindcast correlations obtained from one-year-out cross-validation, which are only shown in case they are significant at a test-level of 5% (the critical value is not constant since it depends on the effective sample size n * defined in Equation (2)). Note that spurious hindcast correlations are marked by a black cross and that the mean and maximum of the significant values is shown in the upper left corner of each panel. With hindcast correlations of up to 0.92, 0.89, 0.88, 0.80 and 0.79 for precipitation amount, DTR, sun hours, cloud cover and wind speed respectively (see left column of Figure 6), the skill of the proposed statistical forecasting method is significant over a large fraction of the study area, this fraction being smallest for the case of wind speed. The climatological (no-skill) hindcast is outperformed by up to 60, 55, 52, 39 and 37% for the five above mentioned variables (see right column of Figure 6), demonstrating that the skill is robust to applying an alternative measure. Note that these results are obtained from detrending the predictor and predictand samples in each step of the cross validation. Using the original/non-detrended timeseries (not shown) yields similar skill levels, indicating that the results are not sensitive to possible trends in the underlying data. Note also that the Pearson correlation between the RSAI and the target variables, i.e. the within sample relationships (not shown) are systematically stronger than the hindcast correlations obtained from cross-validation. The sign of the correlation between the RSAI and the target variables is spatially homogeneous and is shown in the lower left corner of each panel.
Using the RSAI instead of the outlier sensitive SAI (Cohen and Jones, 2011) Figure 5); (c) and (d) Relationship between the Robust Snow Advance Index (RSAI) for October and the DJF-frequency of the above mentioned weather types (in days); the time series for these DJF-counts are detrended and centred to have zero-mean. Also shown are the Pearson (r) and Spearman (rs) correlation coefficients. Correlation coefficients for the original/non-detrended predictand time series are shown in parentheses.
hindcast skill obtained from cross-validation for all applied variables. Table 2 documents this by comparing the areal fraction of locally significant (α local = 0.05) hindcast correlations obtained from the RSAI to the respective areal fraction obtained from the SAI (see columns 2 and 3). For any of the five target variables, the 99th percentile of the 2000 fractions obtained from the randomly re-shuffled RSAI is much lower than the fraction obtained from the 'correct' RSAI (i.e. having the correct temporal order). Hence, field significance (α local = 0.05, α global = 0.01) is given in any case. Moreover, the areal fractions obtained from using the Octobermean NAO (or AO) index as single predictor instead of the RSAI are comparatively low (see columns 4 and 5). This indicates that the hindcast skill stemming from the NAO (or AO) anomaly in October, which potentially could persist throughout the following winter months, is negligible.

Discussion and Conclusions
This study has provided further statistical evidence for the existence of a formerly unknown lead-lag relationship between Eurasian snow cover increase in October and the winter climate on the Iberian Peninsula (Brands et al., 2012). It has been found that an anomalously high increase of Eurasian snow cover in October favours an above normal frequency of cyclonic and westerly flow weather situations over the Iberian Peninsula during the following December-to-February season, whereas the frequency of anticyclonic and easterly flow spatial mean = 4,8 mm spatial mean = 0,7 mm spatial mean = 3,2 hours spatial mean = 5,5 hours spatial mean = 5,4 octas spatial mean = 3,6 octas spatial mean = 3,6 m/s spatial mean = 2,7 m/s spatial mean = 7,5ºC spatial mean = 9,4ºC   situations is below normal. With an explained variance of approximately 60% for both groups of circulation types, this statistical relationship is strong and highly significant (α = 0.01). At the local-scale, this favours below-normal DJF-mean conditions for diurnal temperature range and sun hours, while the corresponding values for cloud cover, wind speed and precipitation amount are above-normal.
On the basis of these results, it has been additionally shown that the above mentioned variables can be skillfully hindcast using simple linear regression in a oneyear-out cross-validation framework. Locally significant hindcast correlations of up to 0.92, 0.89, 0.88, 0.80 and 0.79 where found for precipitation amount, diurnal temperature range, sun hours, cloud cover and wind speed respectively, the corresponding skill patterns being globally significant in any case (α local = 0.05, α global = 0.01). Applying robust linear regression instead of ordinary linear regression (Cohen and Jones, 2011) for calculating October snow cover increase was found to improve the statistical relationship. However, owing to the limited sample size of the applied snow-indices (n = 15 imposed by data availability), we cannot judge the significance of this improvement yet.
The conducted tests for local and global significance and the consistency of the results for a broad range of atmospheric variables on the local and synoptic scale indicate that the described teleconnection is very unlikely to be by chance and that the question posed in the title can be affirmed. Again, the limited size poses some restrictions on this conclusion. First, it was not possible to test the validity of the teleconnection / skill of the statistical forecasting scheme for a large independent time period. Actually, the strength of the statistical link between Eurasian snow cover in October and the DJF-mean AO is known to be non-stationary (Peings et al., 2013) and a similar behaviour would be expected for the teleconnection suggested here. In this context, it is also important to note that the data withheld in each step of the one-yearout cross-validation is not a surrogate of truly independent/future data since, prior to cross-validation, all data pairs had been used for detecting the teleconnection, as is commonly done in this type of studies (see DelSole and Shukla, 2009 and references therein). Consequently, the proposed teleconnection should be re-tested in the future, when larger samples of independent predictor-predictand pairs become available (Labitzke et al., 2006). Finally, it is recommended to assess the atmospheric precursors of October Eurasian snow cover increase in order to challenge the causal relationship suggested in this study.
Our results are expected to be of value for the purpose of statistical seasonal prediction and its applications (Brands, 2013). One evident message is to include the Robust Snow Advance Index as an additional informative predictor of multiple linear or non-linear predictions schemes (Tangang et al., 1997;Rodriguez-Fonseca and de Castro, 2013;Hertig and Jacobeit, 2011;Lorenzo et al., 2011;Folland et al., 2012). Our results are also expected to be of interest for the numerical climate modelling community. First, general circulation models (GCMs) run in seasonal prediction mode are known to have little skill in predicting the boreal winter climate (Doblas-Reyes et al., 2009;Frias et al., 2010;Kim et al., 2012;Doblas-Reyes et al., 2013). Second, transient GCM simulations run over climatic periods of the historical past (Taylor et al., 2012) are known to overestimate the boreal winter westerlies over the North Atlantic Zappa et al., 2013). These GCM-errors might be attributed to poor snow-atmosphere/troposphere-stratosphere coupling and improvements in these fields may consequently help to improve the models (Hardiman et al., 2008;Charlton-Perez et al., 2013).
On the other hand, the purely statistical relationships described in this study are incomplete without assessing the physical background of the teleconnection with idealized numerical model studies. In this context, it is important to note that, using reanalysis data, only approximately 20% of the variance of the boreal winter AO can be explained by downward propagation from the stratosphere (Baldwin et al., 2003). This is in disagreement with the much larger fraction of variance of winter climate anomalies on the Iberian Peninsula that can be explained by October Eurasian snow cover increase (e.g. 60% for the case of weather type frequencies). To put it in another way, if one-way downward propagation accounts for only 20% of the variance of the hemispheric-wide circulation in boreal winter (as described by the AO), how is it possible that two-way troposphere-stratosphere coupling accounts for 60% of variance of the regional winter climate on the Iberian Peninsula? Consequently, alternative teleconnection pathways involving other slowly varying components of the climate system should be taken into account, the ocean being one promising candidate (Barnett et al., 1989;Henderson et al., 2013).
S.B. would like to thank the 'Consejo Superior de Investigaciones Científicas' for financial support. The authors acknowledge the ERA-Interim (http://www.ecmwf.int/ research/era/do/get/era-interim), ECA&D (http://eca. knmi.nl/) and E-OBS (http://ensembles-eu.metoffice. com/) datasets. They are grateful to Dr. Colin Harpham (Climate Research Unit) for providing programming details on the automated Lamb weather typing approach and would like to thank Dr. Nieves Lorenzo (University of Vigo), acting as a reviewer, and one anonymous reviewer for their helpful comments on the former version of this manuscript. Finally, the authors would like to thank Dr. Michael Riemer (University of Mainz) for a critical comment on the interpretation of our results.

Appendix
Consider the regression equation: where y i is the daily snow cover extension for the day x i , β is the regression coefficient, σ is the error scale parameter and e i is the error assumed to be independent and identically distributed. Then, iteratively re-weighted least-squares regression is used as follows (Street et al., 2002): 1. Obtain an initial estimate β for β, as well as the residuals r i by performing a least squares regression on y i = βx i 2. Obtain an estimate σ for σ , where σ = MAD r /0.6745, and MAD r is the median absolute deviation of the residuals from their median. 3. Re-calculate the residuals r i = y i − x i β /σ 4. For any R i < π, a weighting term w i = sin(r i )/R i is defined, where R i = r i / 1.339 * (MAD r /0.6745) * √ (1 − h) and h is the leverage obtained from a least-squares fit. Note that this weighting function was published by Andrews (1974) and that alternative functions (Huber, 2013) yielded virtually identical results. 5. Update the estimate β as well as the residuals r i by performing a least squares regression with the weights w i . 6. Repeat steps (2)-(5) until MAD r is minimized.
The RSAI is then defined as the regression coefficient β obtained from this iterative optimization procedure, performed with the robustfit.m function of the programming environment Matlab. Finally, the 15 RSAI index values obtained for each October are z -transformed to yield standardized anomalies.