Deriving Snow Cover Metrics for Alaska from MODIS

1 Southwest Alaska Network Inventory and Monitoring Program, National Park Service, 240 W 5th Ave., Anchorage, AK 99501, USA; E-Mails: amy_e_miller@nps.gov (A.E.M.); peter_kirchner@nps.gov (P.K.); tammy_wilson@nps.gov (T.L.W.) 2 Geographic Information Network of Alaska, University of Alaska Fairbanks, Box 75-7275 GINA WRRB 111, 909 Koyukuk Dr., Fairbanks, AK 99775, USA; E-Mail: jiang@gina.alaska.edu


Introduction
There is a pressing need for better understanding of snow cover and snow season dynamics at northern high latitudes due to the region's tendency toward greater temperature variability (i.e., Arctic amplification) [1] and its importance to global climate [2].Warmer winter temperatures and reduced cold-season variability have been observed in the mid-and high latitudes of the Northern Hemisphere in recent decades [3,4].Surface air temperature anomalies over the Arctic Ocean have shown a pronounced warming in fall and winter [2] in agreement with model simulations that project stronger winter warming in the Arctic over the 21st century [5] and large decreases in winter temperature variation due to changes in meridional heat advection, snow cover and sea ice extent [2,4,6].
Winter precipitation patterns in the mid-and high northern latitudes can be highly variable and long-term trends difficult to detect.Challenges to identifying trends in snow data include sparsely distributed observations, a lack of long-term records and data discontinuities, low snowfall amounts, redistribution of snow by wind, and repeated snow on-off events [7].In Alaska, there is significant reliance upon modeled climate data products such as the Parameter Elevation Relationships on Independent Slopes Model (PRISM) [8] and sparse in situ observations such as Snow Telemetry (SNOTEL) stations for assessing precipitation due to the limited availability of precipitation data across the state.Snow cover monitoring is one way that we can assess spatial and temporal distributions of precipitation regionally.
Snow cover can be monitored using in situ measurements, modeling studies, and remote sensing applications.Ground-based, in situ monitoring networks constitute a long-term reliable data set for specific locations.These in situ data provide accurate measurements of depth, density, snow water equivalent (SWE) and meteorological conditions and support temporal analyses at point scales [9][10][11], but are of limited use for broad spatial analyses due to their sparse distribution.Snow cover characteristics and timing can also be modeled [12]; however, because of the lack of data on initial conditions, the accuracy of modeling results is often low.Remotely sensed data from satellites provide global spatial coverage and frequent temporal coverage of snow cover but do not measure depth, density or SWE.Therefore, a combination of these methods is required to characterize snow season dynamics at regionally relevant scales [13].The Moderate Resolution Imaging Spectroradiometer (MODIS) sensor provides important optical satellite remote sensing datasets for monitoring snow cover [14] and MODIS snow products are readily available at 500 m spatial resolution for the years 2000 to present.
Three daily MODIS snow cover products, based on different algorithms are globally available.The MODIS Terra Snow Cover Daily L3 Global 500 m Grid MOD10A1 [15] from the National Snow and Ice Data Center (NSIDC) is perhaps the most well-known.The NSIDC also provides the MODIS snow cover MYD10A1 product from the Aqua satellite [16].Both these products have snow cover, snow albedo, fractional snow cover, and quality assessment data.However, the Aqua satellite has a shorter record (July 2002 to present) and sustained a band 6 detector failure just after launch [17].A third product, MODIS Snow Covered-Area and Grain size (MODSCAG) generates fractional snow covered area and snow grain size from MODIS surface reflectance data and is available through the NASA Jet Propulsion Laboratory [18].
A major challenge for optical satellite remote sensing of snow comes from cloud cover.Polar darkness poses an additional challenge in high latitudes.Several methods have been developed to reduce data gaps in MODIS snow products, thereby improving the accuracy of classification of snow.These methods include the use of a combination of Terra and Aqua MODIS data, spatial and temporal filtering algorithms, and the compositing of maximum snow cover extent over multiple days (e.g., eight-day composite climate-modeling grid) [19].
Terra and Aqua satellites cross the equator at local time 10:30 AM and 1:30 PM, respectively, and the three-hour difference between the two satellites may allow MODIS sensors to view different sky conditions (i.e., cloudy or clear) in the same 24-hour period.A combination of Terra and Aqua MODIS data has been found to reduce the uncertainties in snow cover resulting from cloud cover in locations with significant differences between morning and afternoon cloud cover [20][21][22].A second approach to cloud reduction is the use of spatial filtering.A spatial filter based on the state of orthogonal neighboring pixels surrounding a cloudy pixel can be used to reclassify cloudy pixels, often with low levels of error [23].More complicated spatial filtering methods have been used to consider topographic effects such as an eight-neighboring pixel and snow transition elevation method [24].Further dividing scenes into different zones based on topographic character such as basin orientation, major aspect, and land-cover has been used to derive a snowline for each zone [23].Resultant snowlines have then been used to reclassify the higher-elevation cloudy pixels in each zone.Temporal filtering can likewise be applied at varying levels of complexity.Simple temporal filtering using one-day-forward and one-day-backward filtering of cloudy pixels has been used [14] as well as a more complicated temporal filtering based on snow cycle that reclassifies cloudy pixels by their temporal position in the snow accumulation-melt season [23].Compositing maximum snow cover extent can largely mitigate cloudy pixels [25], but temporal resolution is sacrificed in the process.
The aim of this study was to derive a spatially and temporally detailed set of seasonal snow cover metrics for Alaska, portions of western Canada and the Russian Far East for 2001-2013 from MODIS MOD10A1 distributed data products; including snow cover duration, start of snow season, end of snow season (snowmelt) dates, and inter-seasonal variability.In doing so, we applied spatial, temporal and snow-cycle filters to reduce data gaps in the MOD10A1 snow product caused by clouds and polar darkness, and we validated our results against data from three types of in situ weather stations.

Study Area
The study area, located approximately between 53°25'N to 73°15'N and 122°30'W to 173°35'E, includes Alaska and portions of western Canada and the Russian Far East (Figure 1).The total land area is approximately 3.76 × 10 6 km 2 .Elevations across the study area range from 0 m to 6168 m above sea level (a.s.l.), with 80% of the area located at or below 1000 m a.s.l., and less than 1% located above 3000 m a.s.l.
In Alaska, which comprises 1.48 million km 2 or 40% of the total land area in our study area, glaciers, permanent snowfields, and/or barren areas (13.8%) characterize most areas above 1500 m a.s.l.[26,27].Two major mountain ranges, the Brooks Range and the Alaska-Aleutian Range, divide the state into five climate zones [28]: (1) a maritime zone that includes southeastern Alaska, the south coast, and southwestern islands; (2) a maritime-continental zone that includes the western portions of Bristol Bay and west-central zones; (3) a transition zone between the maritime and continental zones in the southern portion of the Copper River basin and Cook Inlet, and the northern extremes of the south coast of Alaska; (4) a continental zone consisting of the northern Copper River basin and the interior Yukon River basin; and (5) an arctic zone, located north of the Brooks Range.
Mean annual temperatures range from approximately 4.5 °C in the maritime areas of southeast Alaska to −12.2 °C along the arctic slope [8].The central and eastern continental interior is characterized by the greatest temperature extremes, with mean maximum and minimum temperatures of 25 °C and −34 °C, respectively.Winter storm tracks originate in the North Pacific, bringing precipitation into southeast Alaska and coastal British Columbia, and to a lesser extent southcentral and southwest Alaska.Precipitation falls primarily as snow in most areas of the state, with mean annual precipitation ranging from 5000 mm in the southeastern panhandle and Aleutian Islands to 150 mm on the North Slope [28].Maximum precipitation occurs in most areas between July and October, with the highest variance in precipitation occurring in the coastal areas of southeast, southcentral and western Alaska [28,29].Winter precipitation is likewise greatest at higher elevations along the Gulf of Alaska, in southeast and southcentral portions of the state [8].The lowest precipitation occurs in the Upper Yukon River basin, but variance there is also significant [30].Dominant vegetation types in Alaska include forest (26.1% of land area) up to 1000 m a.s.l. in the eastern interior, and shrub/scrub (24.7%) and wet and dry tundra (26.6%) up to 1500 m a.s.l.[27].Polar darkness affects the region north of ~62°N but true polar night (i.e., continuous darkness) occurs only north of the Arctic Circle (~66.5°N).At the extreme, on the shortest day of the year, approximately 50% of the study region is affected to some degree by polar darkness, limiting the potential for optical satellite remote sensing of snow cover [22].
Seasonal snow cover classes [31,32] were used to group in situ weather station locations, described in Section 2.1.2,into geographic subsets to examine the effects of broad climate patterns and cloud cover on MODIS-derived snow metrics.These cover classes have been mapped at 0.5° spatial resolution using a binary system of three climate variables that affect snow: winter wind, precipitation, and air temperature [31].Across the study area, seasonal snow cover is predominately classified as tundra (48.4%) or taiga (33.6%) class north of the Alaska-Aleutian Range (~60-62°N) and as maritime (11.9%) or, to a lesser extent, ice (2.6%), prairie (2.5%) and alpine (1.0%) south of the Alaska-Aleutian Range (Figure 1).Tundra and taiga snow cover classes are found where winter temperatures are coldest.Snow in the tundra class is more susceptible to sublimation and redistribution by wind and snowpack is shallower than in the forested taiga [31].In contrast, maritime and alpine snow cover classes are often characterized by deep snowpack and melt features [31] resulting from warmer winter temperatures.Prairie snow is typically thin, except where it has drifted, and is found in windy areas where temperature is transitional between warmer maritime and colder continental environments [31].For the purpose of examining effects of cloud cover on MODIS-derived snow metrics at in situ station weather station locations, we grouped stations into Tundra/Taiga (northern) and Maritime/Alpine (southern-coastal) sites.The two combined classes differ in cloudiness due to elevation, physiography, and proximity to the southern coast and the predominant track of winter storms, with the Maritime/Alpine sites bordering the Gulf of Alaska having warmer, wetter conditions, and higher cloud cover.The Maritime/Alpine class also includes in situ station locations in the ice and prairie cover classes, which together comprise a small number (14%) of all in situ stations.Regional snow cover characteristics are summarized by snow cover class [32] and ecoregion [33], the latter of which is based on physiographic, climatic, and vegetative characteristics.[32].In situ weather station locations are grouped as Tundra/Taiga and Maritime/Alpine to differentiate between colder, drier conditions to the north and warmer, wetter conditions along the southern coast, respectively.

Methods
We calculated 12 snow season metrics across thirteen years (2001-2013) to characterize the spatial and temporal distribution of snow cover at a daily scale using the MODIS Terra MOD10A1 data product.We accomplished this through a three-step process to fill data gaps, primarily caused by clouds and extended hours of winter darkness.We assessed the accuracy of these metrics using seasonally-paired dates of snow onset and snow melt from 244 in situ locations throughout the area, and we present and describe these temporal variables in a spatial context.

MODIS Snow Product
We used Version 5 MODIS Terra Snow Cover Daily L3 Global 500 m Grid data (MOD10A1) [19] from the NSIDC to calculate our snow metrics.MODIS Aqua data were also evaluated for compositing with Terra data but were ultimately not used due to the limited difference between morning and afternoon cloud cover in our study area [34].MODSCAG distributed snow cover products were not available at the beginning of our study in 2011 [35].The Terra MOD10A1 data set contains snow cover, snow albedo, fractional snow cover, and quality assessment data in compressed Hierarchical Data Format-Earth Observing System (HDF-EOS) format with corresponding metadata.The snow cover data are based on a snow mapping algorithm that employs a Normalized Difference Snow Index (NDSI) using at-satellite reflectances in bands 4 (0.545-0.565 µm) and 6 (1.628-1652 µm): and other criteria tests [15,17].The data are distributed as 1200 km by 1200 km tiles of 500 m resolution data gridded in a sinusoidal map projection.We mosaicked and re-projected 28 daily MOD10A1 files into the Alaska Albers Conical Equal Area projection (a tailored and centered conic equal area projection) using the North American Datum of 1983, and then converted them into four stacks of single-band GeoTIFF files for calculation of snow metrics for each snow year.We did not filter canopy effects in forested regions, as the application of temporal and cloud filters have been found to reduce the majority of misclassification errors even in areas of high canopy cover [36].To capture the full seasonality of high latitude snow cover, we use "snow year" defined as August 1 of the previous year (day-of-snow-year 213) to July 31 of the snow year (day-of-snow-year 577) [37,38].

In Situ Snow Onset and Snow Melt Dates
We obtained seasonally-paired dates of snow onset and snow melt between 1 August 2000 and 31 July 2013 from 36 Snow Telemetry (SNOTEL) stations, 16 National Weather Service (NWS) first-order (1 st Order) observing stations, and 192 locations from the Global Historical Climatology Network (GHCN) where snow depth is measured daily (Figure 1).SNOTEL dates were manually derived by Natural Resources Conservation Service staff from snow pillow data and represent the first and last dates of persistent seasonal snowpack [39].1st Order observing station dates represent the first and last dates where any snow depth (>0 cm) was measured [40,41].The GHCN [42] includes snow depth observations from the NWS Cooperative Observer Program (COOP), Automated Surface Observing System (ASOS) airport weather stations, Climate Reference Network (CRN) stations and some volunteer observations from the Community Collaborative Rain Hail and Snow (CoCoRaHS) network.The GHCN data reported here represent the first and last dates that an observer measured at least one-half inch (≥1.27 cm) of snow depth at the time of observation.

Filtering Algorithm
We used a three-step method to fill data gaps in the MOD10A1 product for our study area that are primarily caused by clouds or extended hours of darkness.This three-step method is adapted from techniques described by [14] and [23] and involves spatial filtering, temporal filtering, and a forward/backward filtering approach that is applied to a snow cycle composed of three segments.Each step, described in the following subsections, was performed sequentially.These methods were accomplished using an algorithm [43] created in the IDL programming language with ENVI subroutines [44].

Spatial Filtering
The first step was to apply spatial filtering to daily snow cover data.We reclassified cloud pixels as snow if three of four orthogonal neighboring pixels were snow or as no-snow if the neighboring pixels were snow-free.After spatial filtering, the daily data were then stacked to produce four files per year, one for each of the following: snow cover, albedo, quality, and snow fraction.Occasionally, source data for an entire day or more were not available for our study area.Pixel values for these missing data were not interpolated in the stacked files.

Temporal Filtering
The second step was temporal filtering of each pixel in the stacked snow cover file.Temporal filtering was first used to identify whether pixels were correctly classified as land or water, and then used to identify if pixels were snow-covered or snow-free.If a pixel was classified as ocean or inland water more than 10 times in the snow cover time-series for one snow year, then the pixel was classified as such for the entire snow year; otherwise it was classified as land.The snow flag in the snow cover time-series was used to determine whether cloudy days should be classified as snow or no-snow.Temporal filtering was used to infer the condition of cloud-covered pixels based on values derived one day forward (t + 1) and one day backward (t -1) from a pixel on day t [14].This approach is based on the assumption that if the surface condition (snow or no-snow) did not change from the day before (t -1) to the day after the cloudy day (t +1), then the surface condition likely remained constant on the cloudy day (t), as well.Thus, if the corresponding pixel of both the previous and following day were classified as snow the intervening cloudy pixel was also reclassified as snow [23].Likewise, a cloudy pixel was assigned a no-snow value if the corresponding pixel of both the previous (t -1) and the following (t + 1) days were no-snow.In all other situations, cloudy pixels remained classified as cloud.Pixels classified as snow in the source data were not reclassified.

Snow Cycle Filtering
In the third and final step, we segmented each snow year into a snow accumulation, snow cover, and snow melt period, using methods adapted from [23].Our approach used a time series of accumulation and melt patterns over the course of the snow season, but was not confined to geographic snow zones.Using two snow season definitions adapted from [37], we defined the full snow season (FSS) as the interval between the first snow day (FSD) and the last snow day (LSD) of snow cover, and the continuous snow season (CSS) as the longest interval of the season for which a pixel was snow covered.The CSS is a period for which a pixel was snow covered for at least 14 days [37], with intervening snow-free periods limited to two days or less.If more than one uninterrupted 14-day snow covered event occurred, we selected the longest one of these as our CSS.Unlike [37], we did not constrain the timing of the CSS.
Fractional snowcover (FSC) was not considered when determining FSD and LSD.However, thresholds for FSC (≥0.5) and albedo (≥0.3) were used for determining CSS est.start day and CSS est.end day to divide the time series into three segments.We used FSC ≥ 0.5 because this value was the basis for the binary MODIS snow mapping algorithm for MOD10A1, prior to the addition of FSC in Collection 5 [17], and a commonly used method [23,37].We used albedo ≥ 0.30 [45] to better capture the snow to ice transition, during periods of melting and refreezing, in addition to late season firn and ice.
For snow cycle filtering, we assumed that snow accumulation (segment 1) spanned the period from the first day of the snow year (August 1) to the estimated start of the continuous snow season (CSS est.start day).The snow cover segment (segment 2) was defined as the period from the CSS est.start day to the estimated end day of the continuous snow season (CSS est.end day), and snow melt (segment 3) was the period from the CSS est.end day to the last day of the snow year (July 31).
We calculated CSS est.start day and CSS est.end day for each pixel by identifying all snow days in which the fractional snow cover was ≥0.5 and albedo was ≥0.3 and we forward checked each prospective CSS est.start day, beginning with the first day in the time series.If that day and the days following it comprised a duration ≥14 days classified as snow, cloud, fill, night, or missing data, it was retained as CSS est.start day.If not, we evaluated the next in the time series until we found one that met the above criteria of fractional snow cover, albedo and duration.If none of the prospective CSS est.start day met the continuous snow season criteria, then we used 31 December as the CSS est.start day.A similar method was used to estimate CSS est.end day, except that each CSS est.end day was backward checked until 1 January, which is defined as the CSS est.end day if no other days satisfied the backward check.
After obtaining both CSS est.start day and CSS est.end day, we applied different filtering methods to each snow year segment.Pixels classified as fill, night or missing data were reclassified as cloud prior to filtering.After reclassification, we used a two-step cloud filtering approach to remove the cloud pixels.First, we applied a backward filtering method to each segment, beginning with the last day of the snow calendar year.We reclassified cloud days as snow if a day of the snowmelt segment (segment 3), t, was snow, and the preceding consecutive days (t -1, t -2,…, t -n) were cloud.The same filtering and reclassification approach was applied to the snow cover segment (segment 2), working from the CSS est.end day backward in time.For the snow accumulation segment (segment 1), we reclassified cloud days as no-snow if a day of the segment (CSS est.start day), t, was no-snow, and the preceding consecutive days (t -1, t -2,…, t -n) were cloud.
Second, we applied a forward filtering method to each segment, starting with the first day of the snow year.Beginning with the snow accumulation segment (segment 1), we reclassified cloud days as snow if a day of the segment, t, was snow, and the following consecutive days (t + 1, t + 2,…, t + n) were cloud.We applied the same filtering and reclassification approach to the snow cover segment (segment 2), working from the first day of the segment (CSS est.start day) forward in time.For the snowmelt segment (segment 3), we reclassified cloud days as no-snow if a day of the segment (CSS est.end day), t, was no-snow, and the following consecutive days (t + 1, t + 2,…, t + n) were cloud.

Reclassification of Permanent Snow and Ice
Following the reduction of pixels classified as cloud through backward and forward filtering and reclassification as snow or no snow, we interpreted all land pixels with an absence of no-snow days as representing glaciers or permanent snowfields.Glacier and permanent snowfield pixels were then reclassified as snow for the entire duration of the snow year.

Definitions and Calculations
Twelve snow metrics were calculated using two principle snow season definitions (FSS and CSS) adapted from [37], as outlined in Section 2.2.3.The metrics were compiled as individual GeoTIFF files, one for each snow year, with date or number of days reported for each 500 m pixel in one of 12 bands (Table 1).The first snow day (FSD) and last snow day (LSD) represent the first and last occurrences within a snow year where the value for the MODIS snow cover field is snow.The length of the FSS was calculated as the difference in days between the FSD and LSD.The CSS was constrained to days between the FSD and LSD where there were ≥14 consecutive days of snow with no more than two intervening, consecutive days of no-snow.A CSS segment could start or end with cloud cover.In this case, the midpoint between the first or last cloud day and the following or previous snow day was used to estimate the start or end of the CSS.For example, if the last day of a CSS segment was classified as cloud, we stepped backwards to the last previous snow day in that segment and selected the midpoint between the two days, cloud and snow, as the CSS end date.
In cases of intermittent snow cover, a pixel might yield more than one CSS segment in a snow year.If so, the number of CSS segments (CSN) was documented and we report the first and last day of the longest CSS segment (LCFD and LCLD, respectively).The length of the CSS was estimated as the difference in days between the LCFD and LCLD.Finally, we recorded the sum of all days classified as snow, snow-free or cloud (CD) after filtering, and the total number of days within the CSS.

Accuracy Assessment
In situ snow onset and melt dates from weather stations were converted to a day-of-snow-year for direct comparison to snow metrics.Day-of-snow-year corresponds to the familiar day of year (i.e., 1 to 365), but spans two calendar years (e.g., the first day of the second year is day 1 + 365 = 366).Because our snow year is defined as August 1 to July 31 [37,38], possible values for day-of-snow-year are 213 to 577 (1 August is day of year 213; 31 July is day of year 212 + 365 = 577).Gridded metrics results were masked so that pixels classified as other than land (e.g., lake, ocean) more than 10 times in a snow year were ignored.This was done because of misalignment of the land/water mask (coastline) and snow/cloud misclassification errors that occur at the edges of a snow-covered area and water (e.g., lake/ocean) in the MOD10A1 product [46].Because of the combined uncertainty in the location of in situ data and each MODIS pixel and the problem of comparing point scale data to area scale data, we used the median value in a 4-pixel neighborhood for accuracy assessment [25,47].
We used bias and root mean square error (RMSE) from paired snow onset and melt dates from in situ weather stations and those derived from MODIS over the entire time series (2001-2013) to assess the accuracy of the MODIS-derived metrics: start of FSS (FSD, nobs = 1870); start of CSS (LCFD, nobs = 1856); end of FSS (LSD; nobs = 1881); and end of CSS (LCLD, nobs = 1853).Bias and RMSE metrics were computed by comparing MODIS-derived dates with in situ dates recorded at three types of weather stations (SNOTEL, 1st Order, and GHCN).RMSE and bias are frequently used measures of the difference between values predicted by a model and the values observed from the environment that is being modeled.We defined the error as the in situ date minus the MODIS-derived date and therefore report RMSE in the same units as the input data, i.e., as number of days [48].We used boxplots, and scatterplots of the predicted MODIS-derived snow season dates plotted as a function of the in situ measured dates with the least squares best-fit of the dates to visualize overall bias.We then used a general linear model (GL) to evaluate sources of variation in bias cause by station type (SNOTEL, 1st Order, and GHCN), combined snow class-Tundra/Taiga and Maritime/Alpine; (Section 1.1)-and number of cloudy days (CD; post-filtering cloud cover) using the difference between the in situ date and the MODIS derived date (bias) as the response variable.After determining that the three-way interactions were not significant, we fit the model with all possible two-way interactions between the two factors (station type, snow class) and one covariate (cloud cover) for each of the snow metrics using maximum likelihood (GLM function) in the program R [49].Our objective was to calculate separate means and slope parameters for cloud cover for each weather station type and snow class category, so further reductions of the full model were not tested.

Effectiveness of Cloud-Filling Procedures
Spatial, temporal and snow-cycle filtering were conducted for the entire data set.Results from a subset for the 2012 snow year (1 August 2011 to 31 July 2012) are presented in Figure 2 and Table 2.
The relative area percentage of classes for snow, snow-free, cloud and no data (e.g., polar night) for the study area change significantly from the unfiltered scenes (Figure 2A) to the final filtered scenes (Figure 2D).However, changes in the intermediary spatial and temporal filtering, shown in Figure 2B,C, are much smaller.In the subset for the 2012 snow year, the percentage of pixels in the study region classified as cloud was reduced from 27.7% to 3.1% (−24.6%)following spatial and temporal filtering.At the same time, the percentage of pixels classified as snow increased more than three-fold, from 8.0% to 28.3%, as the number of cloud pixels decreased (Table 2).Pixels classified as fill, night (polar darkness), or missing data were reduced from 57.7% to 51.5%.These results are comparable to other studies that used similar filtering procedures.A study from the Himalaya found that snow cycle filtering was the most effective filtering process in reducing cloud pixels (−25.2%) and that temporal filtering was more effective than spatial filtering [23].Our results are consistent with those of [23] although they reported a proportionally greater reduction in cloud pixels during these latter filtering steps than we achieved in our study.A study from the European Alps [50] also found that temporal filtering was more effective than spatial filtering in reducing the number of cloud pixels.We found that the effectiveness of filtering can vary significantly across space, in agreement with other studies [14,50].After filtering, the number of cloudy days (CD) per pixel per snow year averaged 18 days across the study area.The mean number of cloudy days/year for the combined the Tundra/Taiga and Maritime/Alpine snow classes was 15.8 and 25.5 days, respectively.Cloud filtering was most effective in continental locations north of the Alaska-Aleutian Range and least effective in mountainous terrain and along the southern coastal margin, where cloud cover is highest.There are two additional approaches to reduce cloud cover in MOD10A data that have been applied in other regions but that were not applied to these data.The first method, a zonal snow line approach to cloud cover as described in [20] and [23], was not used because the method is computationally impractical for the size of the region (3.76 × 10 6 km 2 ).A recent study from the European Alps [50] that applied the snow line estimation method reported only a small decrease in median cloud cover, from 37.2% to 35.7%, whereas the temporal filter they applied reduced median cloud cover from 35.7% to 0.1%.As with our study, the increase in the number of clear sky pixels following spatial and temporal filtering resulted in a marked increase in snow-covered pixels in the final snow product [50].The second cloud-reduction method, the compositing of Terra MOD10A, and Aqua MYD10A data (as described in Section 2.1.1)yields little improvement to cloud filtering because there is generally little difference in cloud cover between morning and afternoon at the higher latitudes and often slightly greater cloud cover in the afternoon during the fall and winter [34].

Accuracy Assessment
The difference between in situ and MODIS-derived snow onset dates for the FSD of the full snow season were positive, indicating that the MODIS-derived start of full snow season dates were early (positive bias) relative to the weather station data (Figures 3 left; Table 3).This bias was true for all weather station types, and was more pronounced for locations within the Maritime/Alpine snow classes than in the Tundra/Taiga snow classes.The overall RMSE for FSD was less than 26.0 days at half of all the locations that were investigated (nstations = 221).The onset date for the continuous snow season (LCFD) was late (negative bias) relative to observations from all weather stations, with a greater effect evident in the combine Maritime/Alpine snow (Figures 3 right).RMSE for snow onset dates for the LCFD of the continuous snow season was <21.0 days at 50% of all the locations that were investigated (nstations = 220) and < 11.9 days at half of the SNOTEL stations.MODIS-derived dates were biased early relative to in situ observations for the last day (LSD) of the full snow season for 1st Order weather stations in both snow classes (Figures 4C and 5C), and biased late for the GHCN stations in the Maritime/Alpine snow class (Figures 4E and 5C).The LSD metric performed well in reference to SNOTEL sites, where no significant bias was observed (Figures 4A and 5C).The difference between MODIS-derived snow melt dates for the LSD of the full snow season and in situ snow melt dates was <21.7 days at 50% of all the locations that were investigated (nstations = 233) and <15.4 days at half of the SNOTEL stations.The end of the continuous snow season (LCLD) was biased early relative to the weather station observations for all weather station types in both snow classes (Figures 4B,D,F and 5D).RMSE for metrics snow melt dates for the LCLD of the CSS was <22.5 days at 50% of all the locations that were investigated (nstations = 220) and <11.4 days for half of the SNOTEL stations.Both station type and snow class showed significant interactions with cloud cover (Table S1), indicating that cloud cover affected bias differently in most metric/station type/snow class combinations.In general, increasing cloud cover tended to exacerbate both positive and negative bias at the start and end of the snow season (Table S1) for metric/station type/snow class combinations were bias was large.The effects of cloud cover were minor and mixed for metric/station type/snow class combinations where bias was small (Table S1).The MOD10A1 snow mapping product, the basis of our algorithm for snow metrics, has been shown to have an overall accuracy of about 93% with the main sources of error being land-cover type, snow/cloud misclassification, and thin or patchy snow cover [17].However, the accuracy of MODIS generally only considers cloud-free conditions.Figure 5 shows bias for MODIS-derived snow onset and melt dates by in situ station type with respect to the combined seasonal snow classifications [31,32] discussed in Sections 1.1 and 2.3.2.Bias is greater for stations in the combined Maritime/Alpine snow classes.These stations are primarily located along Gulf of Alaska coast, in areas of steep topography, and/or at developed locations (Figure 1).Bias is less for stations in the combined Tundra/Taiga snow classes.Overall the GL model demonstrated that clouds contributed to bias regardless of whether the MODIS-derived dates were earlier or later than the in situ observations.Although adequate for explaining coarse sources of potential bias, there remains quite a bit of unexplained residual variance in our models.Additional sources of error could include elevation and tree cover [36,51] and scale mismatch between our observations (point estimates) and the large MODIS pixels.In mountainous terrain, which is most prevalent in the combined Maritime/Alpine snow class, there could be considerable variation in the start/end of the snow season within single pixels, contributing to apparent bias.Some in situ snow onset and melt dates may not be well suited for comparison with MODIS derived metrics because of the surrounding topography, land cover, or urban environment.Different thresholds have been presented to infer snow cover from snow depth [20].We addressed this by using a variety of snow depth thresholds (>0 cm, ≥1.3 cm) to approximate the fractional snow cover of 0.5 required by our filtering algorithm for the CSS.However, it is unlikely that stations located in developed locations (e.g., Anchorage Hillside SNOTEL, Juneau Downtown Cooperative Observer Program station in GHCN) accurately reflect snow cover conditions in the adjacent area as sampled by MODIS.Other sources of uncertainty include: polar darkness [22], snow cover in closed-canopy forest [51], steep topography and persistent cloud cover over predominately snow-covered or glaciated areas-all of which are characteristic of our study area.Some of the error we observed (e.g., early FSS onset dates) can be attributed directly to the source data (MOD10A1) in which pixels are sometimes anomalously classified as snow on one date (e.g., August 1) followed by weeks where that pixel is classified as no-snow or cloud [52].It is likely that most of the differences between our metric dates and in situ dates are due to error propagation from the MOD10A1 product and the subsequent amplification of these errors during filtering and re-classification of cloud pixels to snow or no-snow to yield metrics.These combined uncertainties can result in large errors in regions with extended polar darkness or persistent cloud coverage.
MODIS-derived dates for snow melt from this study have recently been compared to snow-free dates derived from Landsat during 1985-2011 for northwest Alaska.Snow melt dates for the end of the continuous snow season (LCLD) were 3.4 days earlier, on average, than snow-free dates derived from Landsat and mean absolute difference between the two data sets was 4.2 days [53].

Full Snow Season
Snow cover metrics with 12 attributes (Table 1) were derived from MOD10A1 snow cover using the methods described in Sections 2.2-2.3.1.The metrics are available as a gridded product through a Web Coverage Service (WCS) for use in Geographic Information Systems at [54] and through Hypertext Transfer Protocol (HTTP) from [55].FSS metrics (FSD, LSD, FLSDR, SD; Table 1) reflect the timing and duration of the period between the first snow fall and last snow melt.The FSS was primarily dependent on elevation, latitude and proximity to coastal areas.Areas of shortest snow season duration (FLSDR) were evident in coastal areas of western and southwestern Alaska, classified predominately as Maritime snow cover [32], and to a lesser degree, in the low elevation river basins in interior Alaska, classified largely as taiga.Mean snow season duration (FLSDR; 2001-2013) for these low-lying areas ranged from 179 to 206 days in the ecoregions of southwest Alaska and the Aleutian Islands (e.g., Kodiak Island, Alaska Peninsula; Figure 6), to 240-248 days in the interior river basins (e.g., Tanana-Kuskokwim Lowlands, Yukon River Lowlands, Copper River Basin).Areas of longest snow season duration (249-311 days) were found in mountainous areas, including the Brooks Range, Alaska Range, Wrangell Mountains, and the Chugach-St.Elias Mountains.Mean FSS snow onset (FSD; Figure 7) and snow melt (LSD; Figure 8) dates also varied with elevation and geographic location.Mean FSD occurred earliest at the highest elevations, followed by the interior uplands and river basins (Figure 7).Coastal areas in southwestern (Aleutian Islands, Bristol Bay Lowlands, Kodiak Island, Alaska Peninsula, and Ahklun Mountains) and western Alaska (Bering Sea Islands and Yukon-Kuskokwim Delta) experienced the latest snow onset, approximately 15 to 55 days later.Patterns in the timing of snow melt (Figure 8) were the inverse of those for snow onset but showed a greater dependence on elevation, particularly in western Alaska (Nulato Hills and the Seward Peninsula) and the Akhlun Mountains.[33].Ecoregion data courtesy of U.S. Geological Survey.

Continuous Snow Season
Over the course of the FSS, snow cover can melt, sublimate, or be redistributed by wind.Three CSS metrics (LCFD, LCLD, LCDR) reflect the timing and duration of the longest period within the full snow season for which a pixel is snow covered (≥14 consecutive days of snow with no more than two intervening, consecutive days of no-snow).The remaining CSS metrics reflect the number of CSS segments (CSN) and the total number of days that fall within these segments (TCD).Collectively, CSS metrics provide a way to monitor intraseasonal variability in snow cover.
The maximum number of CSS segments per pixel during any year (of 2001-2013) ranged from 0 to 8 (Figure 9).The majority (58%) of the study area experienced only one CSS during the years studied; however 31% of the area experienced two CSS segments and 10% experienced ≥ 3 CSS segments, indicating that snow on-snow off events were more common in some areas than others.When summarized by ecoregion [33], mean maximum frequency of CSS segments was highest in coastal areas of southwest and southcentral Alaska (e.g., Gulf of Alaska Coast, Bristol Bay, Kodiak Island, Alaska Peninsula and Aleutian Islands; approximately 2-3 segments per year).Conversely, the most consistently snow-covered areas, characterized by 1 CSS segment per year, occurred at the northernmost latitudes, regardless of elevation (e.g., Beaufort Coastal Plain, Brooks Foot Hills, Kobuk Ridges, Seward Peninsula, Davidson Mountains and Yukon-Old Crow Basin; Figure 9).Although our results show that roughly 40% of the study area experienced intermittent snow cover over the course of a year, intermittent snow cover has largely been ignored or assumed to be unimportant in other studies seeking to characterize snow cover seasonality [22,50].Decreases in winter precipitation over northern and interior Alaska during the last half century, and increases over the southern and western coasts [29], have the potential to lead to more frequent instances of intermittent snow cover across the state.
Interpretation of the CSS metrics is based on the assumption that onset and melt dates for the longest continuous segment within the FSS are representative of the majority of the snow-covered season.For the limited areas within the study area that experience ≥3 CSS segments, this may not be a valid assumption, because multiple CSS segments could be of comparable length, and thus the longest CSS could substantially underestimate the actual snow-covered season.Across the study area, CSS onset dates for the longest continuous segment (mean LCFD; 2001-2013) lagged behind snow onset dates for the FSS (mean FSD) by 5-48 days.The greatest difference between mean FSD and LCFD dates were found for ecoregions bordering the Gulf of Alaska and Bristol Bay (Boundary Ranges, Copper River Basin, Cook Inlet Basin, Aleutian Islands, Alaska Peninsula, Bristol Bay Lowlands, Kodiak Island, Gulf of Alaska Coast, and the Alexander Archipelago), where LCFD lagged behind mean FSD by 22-48 days.We found little difference between CSS and FSS onset dates for ecoregions in the northernmost latitudes (e.g., Beaufort Coastal Plain, Brooks Foothills) and for those at higher (e.g., Chugach-St.Elias Mountains) where the mean CSS onset dates were 5-14 days later than mean FSD.Therefore, the delay in onset of the continuous snow season also appeared to be dependent on latitude and elevation and, to a lesser degree, proximity to the Bering Sea and Bristol Bay.
CSS melt dates for the longest continuous segment (mean LCLD) were earlier than for the FSS (mean LSD) by 2-51 days for 2001-2013.In general, areas that had later CSS onset dates also had earlier CSS melt dates.Exceptions included the Yukon-Tanana Uplands, Tanana-Kuskokwim Lowlands, and Ray Mountains, where CSS melt dates preceded FSS melt dates by 14-22 days.It is possible that earlier CSS melt dates for these regions result from error propagation and subsequent amplification of these errors during filtering as discussed in Section 3.2; however, the accuracy of MODIS-derived metrics in the seasonal snow classes that characterize these regions (Taiga and Tundra) are among the highest observed.

Conclusions
The snow cover metrics derived for Alaska and adjacent areas include dates of snow onset and melt for the full snow season and the continuous snow season.These metrics are presented at a daily scale at 500 m resolution for an area of 3.76 × 106 km 2 .An accuracy assessment using in situ snow onset and melt dates from 244 weather stations where snow depth was monitored showed that MODIS-derived dates generally precede those measured at in situ weather station locations.Modeled mean bias varied from −12.2 to 33.9 days and was dependent on weather station type and snow class.This bias was lowest (−2.6 to 11 days) for SNOTEL stations in Tundra/Taiga snow classes.Differences between MODIS-derived snow metrics and in situ data are <41 days for snow onset and <34 days for snow melt at the 75th percentile for the full snow season.Increasing cloudiness tended to exacerbate both early and late bias.
The gridded spatial data produced by this study has been made publically available for use in Geographic Information Systems through WCS or HTTP.Our analysis of these data show a dependence of snow cover duration on latitude, distance from open (ice-free) ocean and elevation.We also found that snow cover was intermittent in 41% of the study area, primarily in lower latitude areas located along the northern margin of the Gulf of Alaska.
These MODIS-derived metrics provide a useful tool for documenting inter-annual variability and long-term trends in snow cover characteristics.Two other satellite data sets, the Advanced Very High Resolution Radiometer (AVHRR) and the Visible Infrared Imaging Radiometer Suite (VIIRS) may potentially extend these metrics to climatologically relevant timespans.AVHRR data could be used to extend these results back to 1984, at 1.1 km spatial resolution and data from VIIRS could potentially extend these results beyond the MODIS mission, at 750 m spatial resolution.

Figure 1 .
Figure 1.Study area showing locations of in situ snow onset and melt dates for Alaska with respect to seasonal snow class[32].In situ weather station locations are grouped as Tundra/Taiga and Maritime/Alpine to differentiate between colder, drier conditions to the north and warmer, wetter conditions along the southern coast, respectively.

Figure 2 .
Figure 2. Stepwise results from the spatial, temporal, and snow cycle cloud removal process for 2012 snow year.Filtering was applied to snow, snow free, cloud, and non-valid pixels (A) for unfiltered scenes; (B) after spatial filtering; (C) after temporal filtering; and (D) after snow-cycle filtering.Note missing data on May 11, 2012 from data outage due to spacecraft anomaly.

Figure 3 .
Figure 3. MODIS derived snow onset dates for first snow day (FSD) of the full snow season (FSS) and first day of the longest segment (LCFD) of the continuous snow season (CSS) versus in situ dates for 2001-2013.(A) FSD and SNOTEL stations; (B) LCFD and SNOTEL stations; (C) FSD at 1 st order stations; (D) LCFD and 1 st Order stations; (E) FSD and GHCN stations; and (F) LCFD and GHCN stations.Axes are day-of-snow-year, dashed line represents 1:1 agreement and solid lines show correlation.

Figure 4 .
Figure 4. MODIS derived snow melt dates for last snow day (LSD) of the full snow season (FSS) and last day of the longest segment (LCLD) of the continuous snow season (CSS) versus in situ dates for 2001-2013.(A) LSD and SNOTEL stations; (B) LCLD and SNOTEL stations; (C) LSD and 1st order stations; (D) LCLD and 1st Order stations; (E) LSD and GHCN stations; and (F) LCLD and GHCN stations.Axes are day-of-snow-year, dashed line represents 1:1 agreement and solid lines show correlation.

Figure 5 .
Figure 5. Boxplots summarizing bias (observed-predicted) for four MODIS-derived snow season metrics: (A) FSD, (B) LCFD, (C) LSD and (D) LCLD, compared to three in situ weather station types.The boxes depict the median, 1st and 3rd quartiles; the whiskers depict the 10th and 90th percentiles; and outliers are depicted by the open circles.Colors correspond to combined snow class type [32].

Figure 7 .
Figure 7. Mean first snow date (FSD) for the full snow season during 2001-2013 snow years with respect to ecoregion[33].Ecoregion data courtesy of U.S. Geological Survey.

Figure 8 .
Figure 8. Mean last snow date (LSD) for the full snow season during 2001-2013 snow years with respect to ecoregion [33].Ecoregion data courtesy of U.S. Geological Survey.

Figure 9 .
Figure 9. Maximum number of continuous snow season segments (CSN) in any winter during 2001-2013 snow years with respect to ecoregion [33].Ecoregion data courtesy of U.S. Geological Survey.

Table 1 .
Snow metrics, definitions, and band numbers.

Table 2 .
Quantitative assessment of classification after each type of filtering shown in Figure2.

Table 3 .
RMSE, multivariate mean bias and standard error (SE) from the general linear model, and Pearson's correlation coefficients (r) between metrics and in situ dates from weather stations for snow onset and melt during 2001-2013.