The Role of Climate and Land Use in the Changes in Surface Albedo Prior to Snow Melt and the Timing of Melt Season of Seasonal Snow in Northern Land Areas of 40 ◦ N–80 ◦ N during 1982–2015

: The rapid warming of the Northern Hemisphere high latitudes and the observed changes in boreal forest areas affect the global surface albedo and climate. This study looks at the trends in the timing of the snow melt season as well as the albedo levels before and after the melt season in Northern Hemisphere land areas between 40 ◦ N and 80 ◦ N over the years 1982 to 2015. The analysis is based on optical satellite data from the Advanced Very High Resolution Radiometer (AVHRR). The results show that the changes in surface albedo already begin before the start of the melt season. These albedo changes are signiﬁcant (the mean of absolute change is 4.4 albedo percentage units per 34 years). The largest absolute changes in pre-melt-season albedo are concentrated in areas of the boreal forest, while the pre-melt albedo of tundra remains unchanged. Trends in melt season timing are consistent over large areas. The mean of absolute change of start date of melt season is 11.2 days per 34 years, 10.6 days for end date of melt season and 14.8 days for length of melt season. The changes result in longer and shorter melt seasons, as well as changed timing of the melt, depending on the area. The albedo levels preceding the onset of melt and start of the melt season correlate with climatic parameters (air temperature, precipitation, wind speed). The changes in albedo are more closely linked to changes in vegetation, whereas the changes in melt season timing are linked to changes in climate.


Introduction
Climate change over the Northern Hemisphere high latitudes and boreal forest zone has affected the snow and vegetation cover and, thus, the surface albedo [1][2][3][4][5][6][7][8][9][10][11][12]. Previous studies have shown that the duration and timing of the snow melt season has changed differently in different areas [13][14][15][16][17]. The snow cover extent has decreased especially significantly in the spring [3,5,7], and the surface albedo during the melt season months has also decreased, largely due to a decline in the area covered by snow [10]. These changes affect the local and global energy budgets [6]. The changing climate also changes the vegetation. The increased size of vegetation decreases wintertime albedo by covering the land surface and casting larger shadows on snow-covered surfaces. This is particularly relevant in the late winter. Changes in vegetation affect the local climate and the scattering properties of the forest (with larger shadows and increased multiple scattering due to increased forest height or density) [11]. The changes in the vegetation are also linked to changes in the spatial coverage of permafrost. Table 1. GlobCover2009 classes found to be most common within one surface albedo (SAL) resolution unit in the study area and the number of occurrences of each class as the most common land-use class in one resolution unit of melt season data. Water bodies 115239 220

LUC Class Label Number of Occurrence
Permanent snow and ice ERA-Interim data was 0.25°. The role of land use in the trends in melt season albedo and timing is assessed using data from GlobCover2009 [47]. The data was coarsened to the same resolution as the melt season data (0.25°) by choosing the most common land-use class within the melt season grid cell. Figure 1 shows the GlobCover data at CLARA-A2 SAL resolution. The GlobCover land-use classes present in the study area are listed in Table 1.  Table 1. GlobCover2009 classes found to be most common within one surface albedo (SAL) resolution unit in the study area and the number of occurrences of each class as the most common land-use class in one resolution unit of melt season data.

Methods
The melt season parameters were determined by fitting sigmoids to CLARA-A2 SAL pentad (5 day) mean albedo values using non-linear regression [48]. For each grid cell and year the 5-day mean albedo values from the end of January until the end of August were used for the sigmoid fitting. Figure 2 shows two examples of sigmoid fitting. To include all changes in albedo during the melt season, the dates of snow melt at the onset was taken to be the date at which the sigmoid reached 99% of its variation range (i.e., a change of 1% (relative) from the pre-melt albedo level). Likewise the end of the snow melt season was defined to be the date at which the sigmoid reached 1% of its variation range. These thresholds were chosen in order to include the whole dynamic change in surface albedo during melt season. The length of the melt season was then the difference between the start and end date of melt. The albedo values corresponding to the dates of the onset and end of melt were used as the representative albedo values preceding and following the melt season.

Methods
The melt season parameters were determined by fitting sigmoids to CLARA-A2 SAL pentad (5 day) mean albedo values using non-linear regression [48]. For each grid cell and year the 5-day mean albedo values from the end of January until the end of August were used for the sigmoid fitting. Figure 2 shows two examples of sigmoid fitting. To include all changes in albedo during the melt season, the dates of snow melt at the onset was taken to be the date at which the sigmoid reached 99% of its variation range (i.e., a change of 1% (relative) from the pre-melt albedo level). Likewise the end of the snow melt season was defined to be the date at which the sigmoid reached 1% of its variation range. These thresholds were chosen in order to include the whole dynamic change in surface albedo during melt season. The length of the melt season was then the difference between the start and end date of melt. The albedo values corresponding to the dates of the onset and end of melt were used as the representative albedo values preceding and following the melt season.  In some cases the snow melt onset could not be determined, because the melt had already started before the first cloud-free albedo pentad of the year was available for the grid cell in question. The final analyses included only the grid cells for which (1) both the snow melt onset and end days were retrieved successfully; (2) the albedo difference between the start and end date of melt was larger than 5% absolute albedo units; and (3) data meeting these two criteria were available for at least 10 years. The mean values of R 2 and RMSE of the final sigmoid regressions were 0.989 and 5.55 (albedo percentage), respectively. The corresponding median values were 0.993 and 4.75. In the analyses only the grid cells for which R 2 > 0.95 and RMSE < 20 (for the sigmoid fitting) were included. This lead to discarding about 2% of the data. The final dataset consisted of 2.46 million grid cell level melt seasons. Figure 3 shows the number of years with successful melt season retrievals per resolution unit.
The effect of random error in the albedo data on the derivation of the melt season timing was In some cases the snow melt onset could not be determined, because the melt had already started before the first cloud-free albedo pentad of the year was available for the grid cell in question. The final analyses included only the grid cells for which (1) both the snow melt onset and end days were retrieved successfully; (2) the albedo difference between the start and end date of melt was larger than 5% absolute albedo units; and (3) data meeting these two criteria were available for at least 10 years. The mean values of R 2 and RMSE of the final sigmoid regressions were 0.989 and 5.55 (albedo percentage), respectively. The corresponding median values were 0.993 and 4.75. In the analyses only the grid cells for which R 2 > 0.95 and RMSE < 20 (for the sigmoid fitting) were included. This lead to discarding about 2% of the data. The final dataset consisted of 2.46 million grid cell level melt seasons. Figure 3 shows the number of years with successful melt season retrievals per resolution unit.
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 20 The trends for the melt season parameters (start and end times of melt, the length of the melt season and albedo levels before and after the melt season) for each grid cell over the 34 years were detected using linear regression. The trends were determined using rolling 5-year means. Only 5year means with data from at least 3 years were included in the trend fitting, and the fitting was carried out only for grid cells that had at least 20 mean values during the 34 years. This gave 72092 grid cell level estimates of trends in pre melt season albedo. From these 30% had R 2 larger than 0.5. The climatic dependencies of the melt season parameters were studied using ERA-Interim reanalysis data [46] for air temperature at 2 m, precipitation and wind speed prior to melt. The time interval used in the analysis was 14 days prior to melt. Correlation analysis between ERA-Interim climate data and surface albedo or day of onset of melt were carried out only for the grid cells for which data were available for at least 20 years. The analysis was performed by looking at the linear correlation coefficient between the melt season parameter in question and the climatic parameters.
The GlobCover2009 land-use data set ( Figure 1) was coarsened to the same resolution as the melt season data (CLARA-A2 resolution of 0.25°). The area of the boreal forest was determined by looking at GlobCover classes 70, 90 and 100. The tundra areas were determined as the grid cells with landuse class 140, 150 and 200 that are north of 55°N in North America and 65°N in Europe.

Albedo Before and After the Melt Season
The change in albedo of Northern Hemisphere land areas between 40°N and 80°N before the melt season shows large spatial variations ( Figure 4). The changes are concentrated in large areas with homogeneous change characteristics. These areas are listed in Table 2. The areas are shown on a map together with place names in Figure 5. For the observations with clear trends (R 2 > 0.5, 30% of all observations) the pre-melt albedo decreased by 2.4 absolute albedo percentage units on average over the whole study area over the 34 years of record. The statistically reliable trends in pre-melt season albedo are concentrated in the region of the boreal forest zone. This can be seen in Figure 4b, which shows the land-use classes for the resolution units with albedo trends with R 2 > 0.5. The trends in the boreal forest zone show decreasing pre-melt albedo in the southern half of the Central Siberian Plain and Scandinavia, and increasing pre-melt albedo in the north of Mongolia and China. In North America the direction of changes varies over short distances, whereas in Eurasia there are larger areas with similar direction of change in albedo.
Most of the tundra areas show no significant change in albedo prior to melt. This can be seen in Figure 4c, which shows the R 2 values for all albedo trend retrievals. For all the areas with an R 2 value, the retrieval of melt season parameters was successful, but the data show no reliable trends (having R 2 values lower than 0.5). The slope values of these excluded trends are typically close to 0. In many areas the annual variability of the pre-melt albedo values is so large that even 34 years is not long enough to determine a small trend in the albedo value. In southern Eurasian tundra the pre-melt albedo shows weak negative trends, but no trends in the higher latitudes. The Kola Peninsula and northern Finnish Lapland show strong negative trends (−0.29 albedo percentage units per year). The effect of random error in the albedo data on the derivation of the melt season timing was estimated by using 6 different sigmoids with different levels of albedo prior to melt. After constructing artificial data around these sigmoids, the data was modified by introducing relative random error of +12.5% to −12.5% to it. This error is larger than the level of typical variation of albedo values. These data with random error were then used to produce sigmoids and to extract the melt season parameters. The analysis was repeated for 100 different cases for each of the six chosen pre-melt levels with random error. The effect of relative random error in the albedo data on the start date of melt was 1.3 days (standard deviation) and 1.1 days for the end dates of melt.
The trends for the melt season parameters (start and end times of melt, the length of the melt season and albedo levels before and after the melt season) for each grid cell over the 34 years were detected using linear regression. The trends were determined using rolling 5-year means. Only 5-year means with data from at least 3 years were included in the trend fitting, and the fitting was carried out only for grid cells that had at least 20 mean values during the 34 years. This gave 72092 grid cell level estimates of trends in pre melt season albedo. From these 30% had R 2 larger than 0.5.
The climatic dependencies of the melt season parameters were studied using ERA-Interim re-analysis data [46] for air temperature at 2 m, precipitation and wind speed prior to melt. The time interval used in the analysis was 14 days prior to melt. Correlation analysis between ERA-Interim climate data and surface albedo or day of onset of melt were carried out only for the grid cells for which data were available for at least 20 years. The analysis was performed by looking at the linear correlation coefficient between the melt season parameter in question and the climatic parameters.
The GlobCover2009 land-use data set ( Figure 1) was coarsened to the same resolution as the melt season data (CLARA-A2 resolution of 0.25 • ). The area of the boreal forest was determined by looking at GlobCover classes 70, 90 and 100. The tundra areas were determined as the grid cells with land-use class 140, 150 and 200 that are north of 55 • N in North America and 65 • N in Europe.

Albedo Before and After the Melt Season
The change in albedo of Northern Hemisphere land areas between 40 • N and 80 • N before the melt season shows large spatial variations ( Figure 4). The changes are concentrated in large areas with homogeneous change characteristics. These areas are listed in Table 2. The areas are shown on a map together with place names in Figure 5. For the observations with clear trends (R 2 > 0.5, 30% of all observations) the pre-melt albedo decreased by 2.4 absolute albedo percentage units on average over the whole study area over the 34 years of record. The statistically reliable trends in pre-melt season albedo are concentrated in the region of the boreal forest zone. This can be seen in Figure 4b, which shows the land-use classes for the resolution units with albedo trends with R 2 > 0.5. The trends in the boreal forest zone show decreasing pre-melt albedo in the southern half of the Central Siberian Plain and Scandinavia, and increasing pre-melt albedo in the north of Mongolia and China. In North America the direction of changes varies over short distances, whereas in Eurasia there are larger areas with similar direction of change in albedo.
causing more shadowing of the surface and increased multiple scattering. The darkening of the southern tundra prior to melt could be explained by the reported shrubification of tundra [49].
The role of climate change in altering the pre-melt season albedo was studied using the linear fitting between the melt season data and the ERA-Interim reanalysis data [45] on air temperature, precipitation, wind speed and the number of days with maximum temperature above 0 °C, −4 °C and −10 °C for 14 days prior to onset of melt. Figure 7 shows the 34-year trends in air temperature, accumulated precipitation and wind speed. Changes in all these climatic parameters contributed to the changes in the pre-melt albedo (mean R 2 for the whole area being 0.64 and the 80th percentile being 0.79) (Figure 4d). In the area around the borders of China, Mongolia and Russia ( Figure 5) the climatic parameters explained almost all of the albedo change. The mean air temperature was the dominant influence (mean R 2 = 0.51 for the whole area). It was the largest explanatory factor in particular in Yablonovyy and Verhoyansk Mountain Ranges, Northern West Siberian Plain, Kola Peninsula, Baffin Island and Central Siberian Plain.    The mean trends (R 2 > 0.5) for all the melt season and climatic parameters in areas with observed consistently homogenous change. The trends for snowfall and rain were not included in the table due to the fact that the trends for the pre-defined regions presented in the table were so weak that they were less than the yearly variation of the parameter.  Table 2.
Most of the tundra areas show no significant change in albedo prior to melt. This can be seen in Figure 4c, which shows the R 2 values for all albedo trend retrievals. For all the areas with an R 2 value, the retrieval of melt season parameters was successful, but the data show no reliable trends (having R 2 values lower than 0.5). The slope values of these excluded trends are typically close to 0. In many areas the annual variability of the pre-melt albedo values is so large that even 34 years is not long enough to determine a small trend in the albedo value. In southern Eurasian tundra the pre-melt albedo shows weak negative trends, but no trends in the higher latitudes. The Kola Peninsula and northern Finnish Lapland show strong negative trends (−0.29 albedo percentage units per year).
The role of vegetation in the observed pre-melt albedo changes can be estimated by looking at the albedo levels right after the snow has melted, before the vegetation has started greening. Figure 6 shows the trends for albedo after melt as well as the corresponding R 2 values for the trend fitting. The trends in the level of albedo after melt are much weaker than those in the pre-melt season albedo before melt season, which can be expected since the differences in the albedo between different biomes are much smaller than the changes in the snow cover. The post-melt albedo of the northern Eurasian tundra decreased over the study period. In more southerly areas of tundra, however, there are no clear trends. In the boreal forest zone the trends are towards lower albedo values in most of the area, except for the area west of the River Ob in Russia. One potential reason for the darkening of the boreal forest zone after the melt season could be the increased size of trees and denser forests, causing more shadowing of the surface and increased multiple scattering. The darkening of the southern tundra prior to melt could be explained by the reported shrubification of tundra [49].
The role of climate change in altering the pre-melt season albedo was studied using the linear fitting between the melt season data and the ERA-Interim reanalysis data [45] on air temperature, precipitation, wind speed and the number of days with maximum temperature above 0 • C, −4 • C and −10 • C for 14 days prior to onset of melt. Figure 7 shows the 34-year trends in air temperature, accumulated precipitation and wind speed. Changes in all these climatic parameters contributed to the changes in the pre-melt albedo (mean R 2 for the whole area being 0.64 and the 80th percentile being 0.79) (Figure 4d). In the area around the borders of China, Mongolia and Russia ( Figure 5) the climatic parameters explained almost all of the albedo change. The mean air temperature was the dominant influence (mean R 2 = 0.51 for the whole area). It was the largest explanatory factor in particular in Yablonovyy and Verhoyansk Mountain Ranges, Northern West Siberian Plain, Kola Peninsula, Baffin Island and Central Siberian Plain.

Melt Season Timing
In addition to the changes in albedo, the timing of the melt season has also changed (Figures 8-10 and Table 2). The changes are, as with albedo, significant and spatially consistent but they also vary within the study area. The observations with high values for coefficient of determination (R 2 > 0.5) are concentrated in large distinct areas. The changes were in general towards longer melt seasons and earlier onset of melt. The mean start date of melt season in the pixels for which melt season data are available for the whole 34 years, became 6.1 days earlier over the 34 years. Similarly, the melt ended on average 5.2 days earlier and the melt season, therefore, became 1 day longer on average. The majority of the observations showed no clear reliable trends (Figure 11), but in many areas the changes were significant (Figures 8-10). In many areas the inter-annual variation in the start and/or end dates of the melt season were so large that it was not possible to detect a statistically significant trend. In Eurasia all the parameters showed changes over large homogenous areas, but in North America the trends are typically more localized and variable.

Melt Season Timing
In addition to the changes in albedo, the timing of the melt season has also changed (Figures 8-10 and Table 2). The changes are, as with albedo, significant and spatially consistent but they also vary within the study area. The observations with high values for coefficient of determination (R 2 > 0.5) are concentrated in large distinct areas. The changes were in general towards longer melt seasons and earlier onset of melt. The mean start date of melt season in the pixels for which melt season data are available for the whole 34 years, became 6.1 days earlier over the 34 years. Similarly, the melt ended on average 5.2 days earlier and the melt season, therefore, became 1 day longer on average. The majority of the observations showed no clear reliable trends (Figure 11), but in many areas the changes were significant (Figures 8-10). In many areas the inter-annual variation in the start and/or end dates of the melt season were so large that it was not possible to detect a statistically significant trend. In Eurasia all the parameters showed changes over large homogenous areas, but in North America the trends are typically more localized and variable.
in the Southern West Siberian Plain and the Southern Byrranga Mountains. Changes in wind conditions affect snow surface scattering by affecting sublimation, mechanical metamorphism of the surface crystals, and distribution of the snow, thus affecting the albedo and depth of snow, and the length of the melt season. The distribution of impurities, such as litter from vegetation on the snow cover, can change due to wind conditions. Impurities increase both the absorption of solar energy into the snow pack and the melt rate. Precipitation (together with air temperature) correlates with the start of melt particularly in the Yablonovyy Range and south of Taymur Peninsula. The trends in the 5-year rolling mean of climatic parameters in individual grid cells show similar spatial patterns as the melt season parameters. The trends in these areas are summarized in Table 2.     The distribution of all the computed trends for start and end day of melt can be seen in Figure 11. The majority of the trends had R 2 values lower than 0.5. and the rate of change of these were typically close to zero but slightly towards earlier onset and end of melt.
The trends in melt season timing showed no significant dependency on the land use. In the Central Siberian Plain the melt started and ended earlier (Table 2), resulting in earlier melt seasons across the whole region regardless of the vegetation type (both boreal and tundra). The decreasing trends for mean air temperature before the melt onset in the Central Siberian Plain (Figure 7a, Table 2) show similar spatial patterns as the melt season timing parameters. This can be explained by the fact that the air temperatures prior to melt are not derived from the same time of year, but change together with the start date of melt. With earlier onset of melt, the air temperatures are also derived from an earlier period. In the mid-winter the air temperature is more heavily influenced by the lack of heating from the Sun, whereas later in the spring other climatic factors start to affect the air temperature more significantly. Earlier onset of melt can be associated with colder air temperatures prior to melt onset and more rapid change in the air temperature from cold mid-winter values to melting conditions.
In the area around the borders of Russia, Mongolia and China the melt starts earlier and ends later (Table 2), resulting in longer melt seasons. This is also the case for the Canadian Rocky Mountains. In North America, the northern parts of Labrador Peninsula, which are tundra, also show trends towards a shorter melt season and earlier end of melt.
Using the climatic data from ERA-Interim, three parameters (mean air temperature, mean wind speed and accumulated precipitation) are required to explain the changes in the start date of melt, giving a mean R 2 value of 0.65 for the whole study area (Figure 8). In some regions, the mean wind speed and accumulated precipitation (for 14 days prior to melt) are strongly correlated with the starting time of the melt, while the mean air temperature is not (Figure 12), thus supporting the multivariate explanation. For example, wind speed affects the start of melt more than air temperature in the Southern West Siberian Plain and the Southern Byrranga Mountains. Changes in wind conditions affect snow surface scattering by affecting sublimation, mechanical metamorphism of the surface crystals, and distribution of the snow, thus affecting the albedo and depth of snow, and the length of the melt season. The distribution of impurities, such as litter from vegetation on the snow cover, can change due to wind conditions. Impurities increase both the absorption of solar energy into the snow pack and the melt rate. Precipitation (together with air temperature) correlates with the start of melt particularly in the Yablonovyy Range and south of Taymur Peninsula. The trends in the 5-year rolling mean of climatic parameters in individual grid cells show similar spatial patterns as the melt season parameters. The trends in these areas are summarized in Table 2.

Discussion
The magnitude of changes in albedo prior to melt season seem to be linked strongly to vegetation characteristics, whereas the changes in melt season timing seem to be more strongly linked to climatic factors. The influence of vegetation on surface albedo prior to melt can be understood by considering the difference between the albedo of vegetation and the albedo of snow. Even a small increase in the vegetation cover can alter the surface albedo by several absolute albedo percentage units [49,50], whereas changes in climatic factors prior to melt affect the snow depth and/or the surface crystal structure, and have smaller (but still significant) impact on the surface albedo. Changes in the number of days with snow on the trees are an obvious exception to this, putting emphasis on the influence of climatic parameters on surface albedo. Increased precipitation in vegetated areas would mean increased albedo while, in open areas, such as tundra, more precipitation would more strongly affect the depth and surface crystal structure of the snow pack, with a weaker effect on the surface albedo.
The start date of melt is linked to the air temperature. Increasing temperatures would result in earlier melt onset. However, in this study the focus was on the temperature prior to melt. That means that the time period in question changes from year to year. An increase in air temperatures prior to melt onset would mean that at the start of melt the snow pack is already warmer and the surface crystals may have been affected by melt and sublimation. Conversely, lower air temperatures, such as in the area around the borders of Russia, Mongolia and China, would result in a colder snowpack, and, in the case of a thin snow cover, resulting from low rates of precipitation, colder ground. This would mean slower melt. Both the end date of melt and the length of the melt season are affected by the air temperature, snow depth, snowpack characteristics, and the ground temperature. These are all regulated by climatic factors. For changes in vegetation to cause significant changes in the melt season timing, such as observed in this study, the vegetation would need to change considerably, for example from bare ground to tree cover. At larger spatial scales, such changes take more than 34 years.
The relationship between vegetation and climate is not straightforward. Existing studies of changes in vegetation, permafrost and impurities in snow [26,[51][52][53][54][55] reveal similar spatial patterns of consistent change as do the melt season parameters, but the effect is not the same in all areas. None of the aforementioned factors alone are able to explain the changes over the whole study area. In fact, they can have opposite effects on the albedo in different regions.
The inconsistency in the effects of vegetation on winter time albedo can be partly explained by the different response of vegetation to the changing climate. According to Xu et al. [26], in Eurasia the normalized difference vegetation index (NDVI) is positively correlated with warming temperature, whereas in North America the effect varies in different regions and patterns of greening and browning in North America are fragmented [26]. It is noteworthy that while the NDVI data describe growing season conditions and are thus not directly translatable to winter conditions, they should be considered an indicator of changes that may also be visible in winter, such as growth in shrub size and coverage. The observed changes in vegetation cannot be directly translated to changes in surface albedo. According to Myers-Smith et al. [56] shrub growth is responsive to different drivers in different regions. Sturm et al. [50] found that if shrubs protrude above the snow and cover 10% of the surface, the albedo will decrease by 30%. The spatial coverage of continuous permafrost shows similar spatial patterns with many of the trends for the melt season parameters. For example, at the southern edge of Central Siberian Plain and in Labrador Peninsula continuous permafrost ends in the same area as the trends for melt season change. The difference in trends for melt season parameters can be due to changes in permafrost coverage causing changes in land use, such as vegetation or formation of melt ponds, or by different response of vegetation and snow cover to climatic changes in areas with and without permafrost.

Conclusions
The surface albedo before the onset of melt shows clear but spatially varying trends. The clearest trends are observed in the boreal forest zone, and can be as large as 9 absolute albedo percentage units over the 34-year-long study period. This is about 10-15% of the albedo of clean snow cover with no vegetation protruding above the snow surface. At grid cell level, both the albedo prior to melt onset and the start date of the melt season are responsive to changes in air temperature, wind speed and precipitation amount, with air temperature being the most significant driving factor. At most latitudes the mean albedo before the onset of melt has decreased over the 34 years.
The timing of the melt season shows strong rates of change in localized areas. These changes are better explained by climatic factors than land-use changes. Areas with consistent homogenous changes are larger in Eurasia than in North America. In the Central Siberian Plain the melt season takes place earlier than before. In the Canadian Rocky Mountains and the area around the borders of Russia, China and Mongolia, the melt season starts earlier, ends later, and lasts longer than before.
All in all, both the timing of the melt season and the albedo prior to the melt season changed in large areas between 1982 and 2015. In most areas, both the start and end dates of the melt season have advanced, and the albedo prior to melt onset has decreased, indicating a darkening of winter snow surfaces.