2.1 Snow products

2.1.2 ERA5-Land

The ERA5-Land reanalysis has a replay of the land component of the ERA5 climate reanalysis, forced by meteorological fields from ERA5. It has been produced with the ERA5-Land model H-TESSEL (version IFS CY45R1) without coupling the atmospheric model and without directly assimilating observations to make the simulations computationally affordable (Muñoz Sabater, 2019). This allowed the implementation of some improvements for land surface applications such as a finer spatial resolution of around 9 km. The snow model is the same as in ERA5, but observations are not directly assimilated. Neither in situ snow depth measurements nor the IMS product is directly assimilated by ERA5-Land. ERA5-Land is still influenced indirectly by the snow observations (and observations of other variables) assimilated by ERA5 because ERA5 atmospheric variables are used to control the simulated ERA5-Land fields, which is known as ERA5 atmospheric forcing.

Similar to ERA5, ERA5-Land spans 1950 to the present. The main ERA5-Land dataset goes from 1981 to present, and the ERA5-Land back extension covers 1950–1980. Compared to ERA5, ERA5-Land provides SD as a diagnostic parameter in the CDS, besides SWE and ρ_s. However, we did not use the diagnostic SD, and we calculated the daily SD from the hourly SWE and ρ_s to keep consistency with the method applied in ERA5.

Table 1Description of the snow products used in the study.

Download Print Version | Download XLSX

2.2 In situ snow measurements

In situ daily snow observations were obtained from the Global Historical Climatology Network (GHCN) managed by NOAA's National Centers for Environmental Information (NCEI), the All-Russia Research Institute of Hydrometeorological Information-World Data Centre (RIHMI-WDC), and the Environment and Climate Change Canada (ECCC) network.

The GHCN is an integrated database with more than 100 000 stations across the globe providing daily measurements of land variables since 1981 (Menne et al., 2012). The stations are divided into four main groups based on the data source: (i) US collection, the largest one; (ii) international collection, obtained through personal contacts; (iii) government data exchange, data collected through official GCOS or bilateral agreements; and (iv) the global summary of the day from SYNOP reports.

The RIHMI-WDC contains 620 stations measuring snow since 1882 (Bulygina et al., 2011, 2009). They measure both the snow depth at the station and the snow cover fraction in the surrounding region, which is estimated visually every morning. Snow course surveys are also available every 5 or 10 d depending on the season but are not used in the present study. The dataset includes an automatic quality control procedure that flags potentially erroneous snow depth measurements. All values flagged by this procedure were discarded.

The ECCC network is the primary in situ network for monitoring snow cover trends in Canada (Brown et al., 2021). The ECCC network is composed of manual and automated stations. Manual stations based on ruler measurements started in 1883, but global coverage of Canada was not achieved until the mid-1950s. The manual network peaked around the 1980s with more than 1600 stations measuring daily snow depth, but it has been declining since the 1990s due to the closure of stations and curtailment of manual SD-observing programs. Automated stations with Campbell Scientific SR50 or SR50 A sonic snow depth sensors have been replacing manual ones since the 1990s, but these efforts have not compensated for the high number of manual stations closed. Only manual stations fulfilled the selection procedure for this study. As suggested by Brown et al. (2021), gap-filled SD values within 14 d of a measurement were used in this study.

All stations measuring daily SD from each network between 1950 and 2020 were used. Over Canada, this period was reduced to 1955–2015 due to the low number of available ECCC stations before 1955 and after 2015. Only stations with more than 10 snow days per year and at least 90 % of valid years in the study period were exploited. Missing values are frequent in snow measurements due to the practice of only recording days with snow presence (Pirazzini et al., 2018). Filling them systematically with zeros would introduce a negative bias in the measurements since some missing values could be truly missing. To avoid this, only years with less than 5 % of missing days were used. In some stations, missing values had already been filled with zeros. These cases were identified by flagging years without snow in stations with more than 40 snow days per year. Flagged years were removed after visually inspecting their time series.

Based on the previous methodology, a reference dataset of 527 stations (228 RIHMI, 242 GHCN, 57 ECCC) was created (Fig. 1). Currently, 313 of these stations are assimilated into ERA5. These stations were kept for the stability analysis, since their addition to ERA5 may explain some of the discontinuities observed, but were removed from the accuracy analysis to guarantee the independence of the validation set. The reference dataset was manually divided into the following spatial regions based on the snow patterns and the performance of the snow products: ECCC N (northern Canada, 5), ECCC E (eastern Canada, 31 stations), ECCC W (western Canada, 21 stations), GHCN AK (Alaska, 12 stations), GHCN USA-W (western USA, 26 stations), GHCN USA-E (eastern USA, 116 stations), GHCN EU (central Europe, 43 stations), GHCN CH (Switzerland, 7 stations), GHCN NO (Norway, 38 stations), RIHMI EU (European Russia, 49 stations), RIHMI Ural (Ural region, 59 stations), RIHMI Siberia (46 stations), RIHMI S (southern Siberia, 41 stations), and RIHMI E (eastern Russia, 33 stations).

Figure 1Spatial distribution of the stations used in the study.

2.3 Spatial representativeness of in situ snow observations

The spatial representativeness of in situ observations is critical to conduct point-to-pixel validations, particularly when evaluating coarse products such as global reanalyses. The extent to which point observations are representative of their larger surrounding depends on the geophysical variable and the characteristics of the surrounding terrain, among other factors. The spatial representativeness of in situ observations was assessed based on the method proposed by Schwarz et al. (2017) for downward solar radiation measurements. This method uses a high-resolution product to evaluate the variability of a geophysical variable within the pixels of the coarser product being validated. For that, the high-resolution pixel collocated with the station is compared against the mean of the high-resolution pixels contained by the coarser pixel. The method includes (i) a correlation analysis and (ii) an estimation of the spatial sampling error (SSE), which arises when estimating a variable over a large area (i.e., the coarse pixel) from a point observation.

The NOAA's IMS high-resolution 1 km product was used, which has been providing daily binary snow cover since 2014 (Chiu et al., 2020). The coarse products evaluated are ERA5 ($\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$) and ERA5-Land ($\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$). The spatial representativeness was evaluated using all daily IMS maps of snow cover during 2015. For each station, the daily snow cover from the IMS pixel collocated with the station (${\text{SC}}_{\mathrm{d}}^{\text{station}}$) and the mean of IMS pixels contained by coarser ERA5/ERA5-Land pixels (${\text{SC}}_{\mathrm{d}}^{\text{area}}$) were extracted. The coefficient of determination (R²) between both variables was calculated. The spatial sampling error was estimated as the mean absolute deviation (MAD) of daily SC:

This metric slightly differs from that proposed by Schwarz et al. (2017), who used a more conservative 95 percentile instead of the mean. Additionally, a new metric was calculated and referred to as the spatial sampling bias (SSB) to quantify the systematic error introduced:

Both SSE and SSB were originally dimensionless because SC is a binary variable. However, both metrics were multiplied by 365 to analyze the results in terms of annual snow cover duration [d yr⁻¹] which is easier to interpret. Stations with (${\text{SSE}}_{\text{ERA5}}>\mathrm{10}\mathrm{d}{\mathrm{yr}}^{-\mathrm{1}}$) or (${\text{SSE}}_{\text{ERA5Land}}>\mathrm{10}\mathrm{d}{\mathrm{yr}}^{-\mathrm{1}}$) were flagged as low representative and subsequently removed from the validation after visually inspecting their SCD map around the station (Sect. 3.1).

2.4 Validation of snow products

The SD and SCD estimations from the snow products were validated against the reference ground measurements. The daily SD of ERA5 and ERA5-Land was directly compared against the daily SD measurements. For snow cover duration, daily SC had to be calculated first from daily SD values. The NOAA CDR was also added to this second part of the validation. Besides low spatially representative stations, stations falling within pixels masked as sea/ocean in the different products (ERA5 – 18 stations, ERA5-Land – 13 stations, NOAA CDR – 18 stations) were also removed for the validation. Note that some sites met the spatial representativeness criteria but fell within a pixel masked as sea by the product.

The conversion from snow depth into snow cover is one of the most sensitive aspects when validating snow products. In a previous study performed with the same group of stations by JAXA (Hori et al., 2017), a threshold of 2.5 cm was used to calculate SC based on a local minimum found in the station measurements. The snow cover fraction in the surrounding of the station is visually assessed at the RIHMI network (Bulygina et al., 2011). We used these measurements to analyze the correlation between SD at the station and the surrounding SCF (Fig. 2). For a surrounding SCF of 0.5 (50 %), the mean and median SD at the station was 3.95 and 1 cm, respectively. The 2.5 cm threshold used by JAXA falls just in the middle of these values, so we decided to keep the same value in our study. Additionally, we performed a sensitivity analysis to evaluate the impact of the threshold used in the snow cover validation metrics (Fig. A2). We analyzed how the SCD bias changes when varying the threshold from 0 to 10 cm by intervals of 2.5 cm. The analysis was done in 2015 to include the IMS 1 km product.

Figure 2Relationship between the snow cover fraction (SCF) around the station and the snow depth (SD) at the station in all RIHMI sites. Solid and dashed lines represent the mean and median value, respectively. The gray shaded region shows the 5th and 95th percentiles. The red dot represents the threshold of 2.5 cm selected to convert SD into SC.

Download

The threshold used to convert SD into SC varies between snow products. Besides, some products provide SCF and others the binary SC. The latter is needed to calculate the SCD. For ERA5 and ERA5-Land, all pixels with SCF>50 % (SD>5 cm) were considered snow-covered. The NOAA CDR already provides the binary SC but at a weekly resolution. The annual SCD was calculated by considering that if a week was flagged as snow-covered, all the days within that week were snow-covered.

2.5 Analysis of snow cover trends in the Northern Hemisphere

The trends in SD and SCD were analyzed using the in situ observations due to the artificial discontinuities and trends found in the snow products. Stations flagged as not spatially representative were kept in this part of the study. The methodology to calculate SCD from SD was the same as that used for the validation. For each variable, decadal trends and annual anomalies for the period 1955–2015 were analyzed. Compared to the stability analysis, the study period was reduced due to the low number of Canadian stations before 1955 and after 2015. The significance of the trends was evaluated with the Mann–Kendall test (Mann, 1945; Kendall, 1975). Note that the density of stations was too low for a complete analysis of NH snow cover trends. Even in regions with good coverage, the heterogeneous density of the stations as well as their different spatial representativeness also prevents the calculation of spatially representative trends. However, our main goal was to estimate the trend magnitude to evaluate the significance of the artificial trends and discontinuities introduced by each product.

The trends in the total NH SCE were analyzed using the three snow products, taking into account the stability issues detected during the validation. The SCE trends and anomalies were calculated over the period when the three datasets were available simultaneously (1972–2020). In the NOAA CDR, the total NH SCE was calculated by summing the area of all snow-covered pixels. In ERA5 and ERA5-Land, snow-covered pixels were summed taking into account the fraction of the pixel covered by snow.

3.1 Spatial representativeness of in situ snow measurements

The spatial representativeness criteria were met by 387 of the 527 snow stations (Table 2 and Fig. 3). The stations removed were primarily located in coastal and mountain regions. On the coast, stations overestimate the mean snow cover over the coarse reanalysis pixel because they are located over land, while the pixel covers both land and sea (see Fig. 3c). In mountainous regions, the spatial representativeness of the stations decreases due to the large elevation gradients (Fig. 3b). The SSE and SSB were very similar in all the stations, which indicates that most of the error introduced by spatial sampling is systematic. Some of these stations showed sampling errors above 100 d yr⁻¹ (>50 %). Norway was the region most affected by the station removal, with 87 % of its stations discarded due to the combination of an irregular coast surrounded by high mountains. Most of the stations were removed because they did not meet the SSE threshold for the coarser ERA5 grid. However, of the 104 stations that were removed, 17 passed the threshold for ERA5 grid ($\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$) but not for ERA5-Land ($\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$). This was the case in sites with high spatial variability of snow cover, where the mean SC over the coarse pixel agreed with the station value just by chance. Implementing the spatial representativeness test at different spatial resolutions, or taking into account the spatial variability of the geophysical variable, is therefore critical to identify and remove these cases.

Table 2Spatial representativeness metrics (R², SSE – spatial sampling error, SSB – spatial sampling bias) per region in the group of stations selected for the validation after discarding low spatially representative sites.

Download Print Version | Download XLSX

Figure 3(a) Spatial sampling error (SSE) of in situ snow measurements with respect to ERA5 grid. Stations in gray have an IMS snow cover equal to zero. Triangles indicate stations assimilated into the ERA5 model. Spatial distribution of annual snow cover duration (SCD) from the IMS 1 km product in two of the stations flagged as low spatially representative: (b) coastal station and (c) mountain station. Red lines represent the ERA5 grid. The red cross shows the station location. Terrain map tiles by Stamen Design, under CC-BY 3.0. Data by © OpenStreetMap, under ODbL.

The resulting group of 387 stations used for the validation has an average R²>0.91 and SSE<4.01 d in all regions. Non-representative stations were removed for both accuracy and stability assessment, despite the effects of low spatial representativeness on stability being smaller because SSE and SSB are usually constant in time. The spatial representativeness of the stations was not analyzed for the NOAA CDR grid. Despite the version used in this study having a similar resolution to ERA5, the resolution of the original input data (historical NOAA weekly charts) is much coarser (∼190.5 km) and not appropriate for point-to-pixel comparisons. Therefore, this product was not used in the accuracy assessment and was only kept in the stability evaluation as a reference, i.e., to discard potential artifacts/trends in the stations when evaluating the temporal evolution of the bias.

3.2 Temporal stability of the products

The analysis of the temporal stability reveals different step discontinuities and trends in the annual bias of both SD and SCD (Figs. 4 and 5) for the three products evaluated. Additional information about the seasonal stability of the bias is available in Figs. A3–A8.

Figure 4Temporal evolution of the annual bias (product – station) in snow depth (SD) per network. Vertical lines show the years when the potential discontinuities in each product occur.

Download

Figure 5Temporal evolution of the annual bias (product – station) in snow cover duration (SCD) per network. Vertical lines show the years when the potential discontinuities/trends in each product occur/start.

Download

The annual bias of ERA5 significantly decreases in time for both SD and SCD, presenting three negative step discontinuities in 1977–1980, 1991–1992, and 2003–2004, as well as a negative trend between 1980–1991. From 1950 to 2020, the interquartile range (IQR) of the annual bias in SD decreases from [3.5,11.1 cm] to $[-\mathrm{0.2},\mathrm{0.4}\mathrm{cm}]$ at RIHMI, from [0.3,1.9 cm] to $[-\mathrm{0.2},\mathrm{0.2}\mathrm{cm}]$ at GHCN, and from [0.1,5.3 cm] to $[-\mathrm{0.3},\mathrm{0.3}\mathrm{cm}]$ at ECCC. The greater initial bias and greater decrease at RIHMI stations are explained by the longer snow season and deeper snowpack over Russia (Bulygina et al., 2011, 2009). The decrease of the bias is more evident in DJF and MAM seasons (Figs. A4 and A5). The magnitude of the bias in SCD is more similar between both networks, with the bias IQR decreasing from [5,25 d] to $[-\mathrm{11},-\mathrm{2}\mathrm{d}]$ at RIHMI, from [2,26 d] to $[-\mathrm{13},-\mathrm{1}\mathrm{d}]$ at GHCN, and from $[-\mathrm{5},\mathrm{15}\mathrm{d}]$ to $[-\mathrm{15},-\mathrm{3}\mathrm{d}]$ at ECCC. The negative change of the bias in both SD and SCD is driven by three stepwise discontinuities. They appear in both networks except for the 1992 discontinuity, which is only observed at RIHMI stations. The hypothesis of an artificial discontinuity in RIHMI in situ measurements was discarded, because the 1992 step discontinuity was not observed in the other products, ERA5-Land, and NOAA CDR. Instead, the three discontinuities are more likely caused by the assimilation of new observations into ERA5. In both 1978–80 and 1992 discontinuities, not only the median bias but also the bias variability, a measure of ERA5 random error, was reduced.

The ERA5-Land bias also decreases temporally for both SD and SCD but in a more gradual way without showing any stepwise discontinuity. The absence of discontinuities is explained by the lack of direct data assimilation in the ERA5-Land model. ERA5-Land is still indirectly influenced by the observations assimilated in ERA5 because ERA5 is used to control the simulated ERA5-Land fields (atmospheric forcing). Therefore, the gradual negative trends in ERA5-Land could be indirectly caused by the three stepwise discontinuities observed in ERA5. Despite being more stable than ERA5, ERA5-Land always exhibits a positive bias but larger bias variability in both SD and SCD. Both the magnitude and variability of ERA5-Land bias are comparable to that of ERA5 before 1980 when no data were being assimilated, whereas ERA5 has clearly been outperforming ERA5-Land since 1992. Despite having a finer spatial resolution and being tailored to land surface applications, the quality of ERA5-Land snow estimates is constrained by the lack of data assimilation.

The NOAA CDR showed a positive overestimation in SCD and a large bias variability. Both issues were somewhat expected. The positive bias could be explained by changes in the method used to produce the snow charts, which since 1999 have considered a pixel to be snow-covered when only 42 % of the IMS pixels within the pixel were snow-covered. On the other hand, the large bias variability could also be related to the coarse resolution of the original product (∼190.5 km), making it inappropriate for point-to-pixel comparisons. This was already stated in the product manual, which recommends using this CDR only for SCE studies over large regions (Estilow et al., 2015). Anyway, as discussed in Sect. 3.1, we decided to use the NOAA CDR in the stability assessment because the goal was not to make point SCD estimations but to include a satellite product in the comparison that helps to discard artificial trends and discontinuities in the in situ measurements. In this sense, the NOAA CDR shows an overall good temporal stability in spring and summer, but a positive trend in the bias has been observed since 1990 in fall (Fig. A6) and winter (Fig. A7). The positive trend in fall has been previously reported as problematic in several studies (Brown et al., 2017; Brown and Derksen, 2013; Derksen, 2014; Hori et al., 2017) and it is further investigated in Sect. 3.2.2.

3.2.1 ERA5 stepwise discontinuities

The magnitude of each ERA5 discontinuity is estimated by calculating the difference in the bias between the 4 years after and before the discontinuity occurred (Fig. 7 and Table A1).

The 1977–1980 discontinuity is the most important overall and could be explained by the assimilation of the first satellite products and in situ observations into ERA5 (Fig. 6). As reported in the ERA5 documentation web page (ECMWF, 2021), significant improvements were observed in the ERA5 forecast skill after 1978 over regions with scarce conventional observations. Carrying out simulations prior to the satellite era is the main challenge of the ERA5 back extension. The annual Δbias is larger in SD ($\text{RIHMI}=-\mathrm{19.4}$ %, $\text{GHCN}=-\mathrm{24}$ %, $\text{ECCC}=-\mathrm{49.8}$ %) than in SCD ($\text{RIHMI}=-\mathrm{2.5}$ %, $\text{GHCN}=-\mathrm{7.5}$ %, $\text{ECCC}=-\mathrm{8.2}$ %). The magnitude of the discontinuity is correlated with the snow depth, having more impact at RIHMI stations and particularly at mountain regions such as the Alps, southern Russia, or Norway. Similarly, the most affected seasons (ranked from largest to smallest) are those with more snow: DJF, MAM, and SON. A positive ERA5 bias in SD of similar magnitude was also observed in SD (Orsolini et al., 2019) and SWE (Bian et al., 2019) over the Tibetan Plateau, a region where neither in situ observations nor the IMS product is assimilated. Orsolini et al. (2019) suggested that the most likely cause was an excessive snowfall precipitation over the Tibetan Plateau, discarding other effects such as snow sublimation due to blowing snow or the SCF threshold. In this study, ERA5 also shows a large positive bias in periods with low data assimilation (before 1980). Besides, ERA5-Land, which does not directly assimilate observations, also shows a predominantly positive bias in SD. As suggested by Orsolini et al. (2019), the most likely cause of the snow depth overestimation in both ERA5 and ERA5-Land could be a precipitation bias, which is only corrected in ERA5 by the assimilation of snow depth observations after 1979.

Figure 6Temporal evolution of the in situ snow stations and observations assimilated by ERA5: (a) all stations assimilated by ERA5 and (b) stations assimilated by ERA5 within our validation set.

Download

Figure 7Change in the ERA5 bias in snow depth (SD) and snow cover duration (SCD) during 1977–1980, 1991–1992, and 2003–2004 discontinuities. The 4 years before and after the discontinuity are compared ($\mathrm{\Delta }\text{bias}={\text{bias}}_{\text{after}}-{\text{bias}}_{\text{before}}$). (a) Seasonal change in the bias per network. (b) Annual change in the bias per station.

The 1991–1992 discontinuity presents a similar seasonal pattern in Δbias_SD than the one in 1977–1980, having more impact (sorted from largest to smallest) in DJF (−52.8 %), MAM (−33.1 %), and SON (−8.1 %), respectively. The main difference is that the 1991–1992 discontinuity is only observed over Eurasia. This step is most likely caused by the assimilation of new in situ observations in Russia and China starting from 1992 (Fig. 6). The assimilation of these observations further corrects the large positive bias exhibited by most Russian stations ($\mathrm{\Delta }{\text{bias}}_{\text{SD}}=-\mathrm{18.2}$ %, $\mathrm{\Delta }{\text{bias}}_{\text{SD}}=-\mathrm{6.3}$ %), falling within a similar range to that observed over Europe and North America. Similar to the 1980 discontinuity, this discontinuity not only reduced the median bias but also its variability.

The 2003–2004 discontinuity, already reported by Mortimer et al. (2020), is caused by the assimilation of the satellite-based IMS product. Compared with the previous ones, this discontinuity has a greater impact on snow onset and snowmelt detection than on snow depth. The change of the bias is larger in SCD than in SD, and spatially, the discontinuity is larger at GHCN stations ($\mathrm{\Delta }{\text{bias}}_{\text{SD}}=-\mathrm{17.2}$ %, $\mathrm{\Delta }{\text{bias}}_{\text{SCD}}=-\mathrm{16.1}$ %) than at RIHMI ones ($\mathrm{\Delta }{\text{bias}}_{\text{SD}}=-\mathrm{0.6}$ %, $\mathrm{\Delta }{\text{bias}}_{\text{SCD}}=-\mathrm{1.5}$ %). This could be explained by how the IMS product is integrated into the ERA5 model. Snow-covered pixels are assimilated as 5 cm of snow depth, explaining the relatively low impact in snow depth and the large improvements in the detection of the start/end of the snow season, which are more evident in SON, MAM, and regions with a short snow season. Particularly large changes are also observed in coastal or mountain regions (Rocky Mountains, southern Siberia), which could be related to the benefits of assimilating a product with a finer resolution over these regions where assimilated snow observations are also scarce. Note that the IMS product was only assimilated below 1500 m, so large improvements observed in mountain regions mostly occur at stations located in the valleys.

Figure 6 shows that the observations assimilated into ERA5 also increased significantly in 2001–2002 (more European stations), 2011–2015 (new stations in eastern Europe and the Scandinavian Peninsula), and 2019–2020 (new Kazakh stations and more observations in China and Europe). However, discontinuities related to these changes were observed in neither the SD/SCD bias time series nor the global SCE time series.

Table A1 summarizes the number of stations where ERA5 shows acceptable stability according to the GCOS requirements. This metric was only calculated for SD (stability = 10 mm) because no explicit requirement is made for SCD. The three discontinuities introduced a Δbias_SD above this threshold in most of the stations, particularly if analyzing the winter trends alone. Only 10.6 %, 12.2 %, and 58.2 % of RIHMI stations; 41.1 %, 69.7 %, and 50.0 % of GHCN stations; and 0 %, 67.7 %, and 69.7 % of ECCC stations were below the stability threshold in winter during the 1977–1980, 1991–1992, and 2003–2004 discontinuities, respectively. Note that these values would be even larger if looking only at snow-covered days in regions such as the USA or Europe where snow does not last the full winter season.

3.2.2 NOAA CDR bias trend in fall and winter

The NOAA CDR exhibits a positive trend in the SCD bias during SON and DJF (Fig. 8 and Table A2). The trend steadily starts from 1990–1995, almost extending to present day. The positive trend in the bias explains the artificial positive trend in the snow cover extent observed in Fig. 11 during SON, which was already documented by different authors (Hori et al., 2017; Brown and Robinson, 2011; Brown and Derksen, 2013). Brown and Robinson (2011) reported that the positive October SCE trend was an artifact of the NOAA CDR since this positive trend was opposed to several independent snow products and to in situ measurements (Peng et al., 2013). They suggested that the most likely cause could be the increase of satellite data ingested by the NOAA CDR analysts, as well as the increase in the temporal and spatial resolution of these products, which led to a more accurate snow onset detection. In the same line, Mudryk et al. (2017) found that the NOAA CDR trends in October and November are non-physical and not consistent with other datasets. Our study corroborates the existence of a significant positive trend in the bias around 5–10 d per decade in many Eurasian stations and in some stations of northern USA. Additionally, our study reveals that the same trend in the bias appears in DJF in some European (10 d per decade) and eastern USA stations (7.2 d per decade), regions where snow onset takes place later than in Russia. Figure 5 also corroborates the fact that there is no step discontinuity related to the transition between the two methodologies in 1999, though the positive trend in the bias may have been aggravated after 1999 due to the improved resolution and the increasing number of satellite products ingested by the IMS product.

Figure 8Decadal trend of the annual bias in seasonal snow cover duration (SCD) of the NOAA CDR from 1992 to 2015 (a) per network and (b) per station. Only significant trends (p<0.05, Mann–Kendall) are shown. The MAM map was excluded due to the lack of an artificial trend globally during that season.

Brown and Derksen (2013) suggested that the opposite effect during the spring season could be expected but was not observed. Theoretically, improved detection of snow melting could lead to a stronger spring trend, introducing an artificial negative trend in the CDR. In this line, Derksen (2014) reported a tendency of the NOAA CDR to map less snow in spring since 2007 than the multi-dataset composed by the NOAA CDR, MERRA, and ERA-Interim. Mudryk et al. (2017) also found that the NOAA CDR has a spring trend stronger than other datasets. We analyzed this issue by evaluating the snow cover duration trends in spring. Negative trends in the spring bias only appear in some Russian stations (Fig. 8a). However, the number of stations showing significant trends in spring is smaller, and the magnitude of these trends is lower than those in fall and winter. Despite the fact that this issue could exist in some regions, the impact at a global scale is negligible (Fig. A8).

3.3 Spatial accuracy of the reanalyses after the last ERA5 discontinuity

Figure 9 shows the performance of ERA5 and ERA5-Land after the last ERA5 discontinuity (2005–2020). Only spatially representative stations that are not assimilated by ERA5 are used (152 of 527). The ERA5 estimations are mostly unbiased for SD, with the annual IQR bias within $[-\mathrm{0.1},\mathrm{0.1}\mathrm{cm}]$ in most regions. Large positive biases only remain over the mountains (Rocky Mountains, southern Russia ranges), which could be related to not assimilating IMS snow above 1500 m. On the contrary, ERA5-Land constantly overestimates SD in most regions, with a bias IQR of [0.8,2.9 cm]. Despite its finer resolution, ERA5-Land quality also degrades in the mountains. Regarding the absolute error, ERA5 shows an RMSE below 1.5 cm in most stations that increases up to 12 cm in mountain stations. On average, 82.6 % of daily ERA5 snow depth values meet the GCOS accuracy requirements, while this number decreases to 10.5 % for ERA5-Land.

Figure 9Bias (product – station) in snow depth (SD) and snow cover duration (SCD) after the last ERA5 discontinuity (2005–2020). Stations assimilated by ERA5 have been excluded.

In SCD, ERA5 presents a constant underestimation (IQR) of around $[-\mathrm{9.4},-\mathrm{5.5}\mathrm{d}]$, while ERA5-Land keeps overestimating [2.4,11.2 d]. As mentioned above, the SCD bias strongly depends on the threshold used to convert SD to SC. Both ERA5 and ERA5-Land use a threshold (5 cm) larger than the one applied to the stations (2.5 cm). This could explain why ERA5 has a negative SCD bias despite having an unbiased snow depth. Indeed, when the ERA5 threshold is applied to the stations (Fig. A2), the ERA5 SCD bias is close to zero in the three networks. We could be tempted to use the same threshold for in situ and satellite data. However, the thresholds applied by the products need to be validated as well, and we can only do it by deriving independent thresholds for the station measurements. In this study, we have used RIHMI visual observations of snow cover in the station, but other data sources such as high-resolution satellite imagery could be used.

We investigated this issue further with a sensitivity analysis that evaluates how the SCD bias changes with different snow depth to snow cover thresholds (Fig. A2). The magnitude of the SCD bias is similar between networks, suggesting a good consistency between their measuring protocols. However, the magnitude of the SCD bias strongly varies between products. When a threshold of 2.5 cm is used, the mean SCD bias varies as follows: 24.8 d (NOAA CDR), 14.3 d (IMS), 8.0 d (ERA5-Land), and −6.7 d (ERA5). These differences are the result of the different thresholds applied by the products, as well as their different snow depth biases (in the case of reanalysis). Orsolini et al. (2019) already pointed out that the different thresholds applied by reanalysis datasets were one of the main limitations for inter-comparing them. The sensitivity analysis also shows that changing the threshold by 1 cm leads to changes in the annual SCD bias of around 2–3 d. These changes are constant between products but vary between networks (ECCC=2.8–4.3 d cm⁻¹, GHCN=1.8–2.1 d cm⁻¹, RIHMI=2.6–3.2 d cm⁻¹) due to the different snow conditions in each region. Stations with many SD values near the threshold will be more affected by a change in the threshold.

3.4 Snow cover trends in the Northern Hemisphere

Linear decadal trends in SD and SCD were calculated annually and seasonally over the period 1955–2015 using data from the ground stations (Fig. 10 and Table 3). The temporal representativeness of the linear trends was further analyzed by plotting the temporal evolution of the anomalies per spatial region (Figs. A9 and A10).

Figure 10(a) Annual and (b) seasonal decadal trends in snow depth (SD) and snow cover duration (SCD) from 1955 to 2015 based on in situ measurements. Only statistically significant trends (p value<0.05, Mann–Kendall) are shown.

Table 3Annual and seasonal decadal trends (median with its 95 % CI) in SD [cm per decade] and SCD [d per decade] from 1955 to 2015. Only statistically significant trends (p value<0.05, Mann–Kendall) are shown. N depicts the number of stations with significant trends in each region. The median was not calculated in regions that have fewer than five stations with significant trends.

Download Print Version | Download XLSX

Figure 11Annual (a) and seasonal (b) NH snow cover extent (SCE) [10⁶ km²].

Download

The SD trends show large spatial variability. Significantly negative trends are observed over Europe (Norway – −0.9 cm per decade, central Europe – −0.1 cm per decade) driven by a strong decrease in winter SD, particularly between 1980–1990. On the contrary, significantly positive SD trends are observed in most of Russia, specifically over the Ural region (+0.9 cm per decade), Siberia (+1.3 cm per decade), and the Sea of Okhotsk (+1.7 cm per decade), driven by a strong increase of both winter and spring snow depth. These trends agree with those reported by Brown et al. (2017) for the Russian Arctic over the 1966–2014 period (SDmax, +0.7 cm per decade). However, as mentioned by Brown et al. (2017), a tipping point is observed around 2000 that reverses the SD increase during the latest years in some Russian regions (e.g., European Russia, Ural region). Negative trends in SD are also observed in eastern USA and in most of Canada (−0.9 to −1.6 cm per decade)

The SCD trends are more spatially homogeneous. A predominantly negative SCD trend of around $[-\mathrm{2},-\mathrm{4}$ d per decade] is observed globally, driven by a strong negative trend during the melting season. The largest reductions in annual SCD are observed over Europe (Norway – −6 d per decade, central Europe – −2.9 d per decade, European Russia – −3.8 d per decade). In Russia, most regions experience a decrease in annual SCD despite their positive SD trends (Siberia – −2.3 d per decade, southern Siberia – −2.2 d per decade). Only a few stations in the Ural region and the Sea of Okhotsk show a longer snow season during the last 70 years. When the trends are recalculated for the period 1981–2020, an acceleration of the SCD decrease is observed as well as an increasing importance of the later snow onset in the annual SCD trends. Again, few stations in eastern USA show significant trends. The low number of significant trends compared to that reported by Knowles (2015) could be explained by a recent recovery in winter and spring SCD starting around 2000 (Fig. A10). Still, the few significant trends observed in the USA are predominantly negative with some exceptions around the Great Lakes that Knowles (2015) attributed to an increased precipitation pattern. The SCD trends are also consistently negative in most of Canada (−1.5 to −5.3 d per decade), driven by negative DJF trends in coastal regions and negative MAM trends in inland and polar regions.

The large spatial variability in SD trends is explained by the non-linear interactions between temperature and precipitation (Brown and Robinson, 2011). At high latitudes, increasing temperatures lead to increasing precipitation due to a moister climate (Thackeray et al., 2019), but snowfall depends on the precipitation phase as well. In relatively warmer climates and maritime regions (e.g., central Europe, Scandinavia, European Russia), negative SD trends could be related to a shift in the form of precipitation towards a rainfall-dominated winter (Luomaranta et al., 2019). On the contrary, in colder and drier climates such as Siberia, snow accumulation is limited by moisture availability (Kunkel et al., 2016), so the positive SD trends may be due to warmer and moister weather (Bulygina et al., 2009) and/or to more extreme snow events, which are more likely to occur below the freezing point (Kunkel et al., 2016). Despite these heterogeneous SD trends, SCD trends are consistently negative globally. The SCD reductions from 1955 to 2015 are mainly driven by an earlier melt that is strongly correlated with the increasing spring temperatures amplified by the snow-albedo feedback (Brown et al., 2017; Luomaranta et al., 2019; Bulygina et al., 2009; Matiu et al., 2021). In regions such as Europe, both SD and SCD are decreasing, with the trend towards shallow snow depth amplifying the shorter snow season. In Russia, spring SCD is also decreasing despite the positive trends in SD. This means that the spring melt driven by warming temperatures overrides any increase in snow accumulation during winter.

Larger variability has been reported for SCD trends during the snow onset season (Brown et al., 2017). However, as also recently suggested by Mudryk et al. (2020), this study evidences the increasing importance of negative SON trends in regions such as Europe, Russia, and the Rocky Mountains, where they have a higher impact than spring trends during recent years.

The Northern Hemisphere presents an average annual SCE of 23.9×10⁶ km² (NOAA CDR) over the 1972–2020 period (the common period between the three products). The three products show an annual decrease in NH SCE (Fig. 11), though the SCE trends should be interpreted cautiously, taking into account the trends and discontinuities discussed in Sect. 3.2. The NOAA CDR is the product typically used for assessing the NH SCE trends. It shows the smallest trend (1972–2020) in annual SCE overall ($-\mathrm{0.15}\times {\mathrm{10}}^{\mathrm{6}}$ km² per decade, −0.63 % per decade), which is driven by a significant decrease in MAM ($-\mathrm{0.61}\times {\mathrm{10}}^{\mathrm{6}}$ km² per decade, −2.13 % per decade) and JJA ($-\mathrm{0.71}\times {\mathrm{10}}^{\mathrm{6}}$ km² per decade, −14.2 % per decade). These seasons are when most snow melts, and again, these reductions are strongly related to the increasing temperature and the snow-albedo amplification. The small decrease in annual SCE compared to that in MAM and JJA is explained by the SCE positive trends of +0.62 and $+\mathrm{0.19}\times {\mathrm{10}}^{\mathrm{6}}$ km² per decade in SON and DJF, respectively. This is due to the artificial positive trend in SON and DJF SCD described in Sect. 3.2.2. In this sense, Derksen (2014) estimated that the artificial trend in October SCE for this product could be around 1×10⁶ km² per decade, which would revert the sign of the SON trend. Other snow satellite products such as JAXA's GHRM5C (Hori et al., 2017) have also reported negative trends of −0.94 and $-\mathrm{0.39}\times {\mathrm{10}}^{\mathrm{6}}$ km² per decade for SON and DJF, respectively, during 1980–2020. Negative trends have been observed as well in SON for the group of stations used in this study. All of this corroborates the underestimation of the snow cover retreat in the fall and winter seasons by the NOAA CDR product.

The ERA5-Land reanalysis has better stability than ERA5, showing a small negative trend caused by the ERA5 atmospheric forcing. This is somewhat corroborated in terms of SCE, with ERA5-Land showing just a slightly smaller trend in MAM than NOAA CDR (−0.54 vs. $-\mathrm{0.61}\times {\mathrm{10}}^{\mathrm{6}}$ km² per decade). On the contrary, ERA5 strongly overestimates the SCE decrease throughout all seasons, showing the largest negative trend in annual SCE ($-\mathrm{1.07}\times {\mathrm{10}}^{\mathrm{6}}$ km² per decade). Among the three ERA5 discontinuities detected, the assimilation in 2004 of IMS snow products has the largest impact overall, leading to a large step discontinuity of around −13 % (annual SCE) and −30 % (SON SCE) in just 1 year (Fig. 11b). As discussed above, the large impact of IMS on the onset and melting period is explained by how this product is assimilated into the model. Overall, ERA5 should be avoided to analyze the NH SCE trends before 2004. The ERA5-Land reanalysis has better stability but still overestimates the actual snow cover retreat.