Snow season variability in a boreal-Arctic transition area monitored by MODIS data

The duration and extent of snow cover is expected to change rapidly with climate change. Therefore, there is a need for improved monitoring of snow for the benefit of forecasting, impact assessments and the population at large. Remotely sensed techniques prove useful for remote areas where there are few field-based monitoring stations. This paper reports on a study of snow season using snow cover area fraction data from the two northernmost counties in Norway, Troms and Finnmark. The data are derived from the daily 500 m standard snow product (MOD10A1) from the NASA Terra MODerate Resolution Imaging Spectroradiometer (MODIS) sensor for the 2000–2010 period. This dataset has been processed with multi-temporal interpolation to eliminate clouds. The resulting cloud-free daily time series of snow cover fraction maps, have subsequently been used to derive the first and last snow-free day for the entire study area. In spring, the correlation between the first snow-free day mapped by MODIS data and snow data from 40 meteorological stations was highly significant (p < 0.05) for 36 of the stations, and with a of bias of less than 10 days for 34 of the stations. In autumn, 31 of the stations show highly significant (p < 0.05) correlation with MODIS data, and the bias was less than 10 days for 27 of the stations. However, in some areas and some years, the start and end of the snow season could not be detected due to long overcast periods. In spring 2002 and 2004 the first snow-free day was early, but arrived late in 2000, 2005 and 2008. In autumn 2009 snowfall arrived more than 7 days earlier in 50% of the study area as compared to the 2000–2010 average. MODIS-based snow season products will be applicable for a wide range of sectors including hydrology, nature-based industries, climate change studies and ecology. Therefore refinement and further testing of this method should be encouraged.


Introduction
The duration of the snow-free season controls the rhythm of ecosystems and feedback to the climate system. This is particularly the case in boreal and Arctic-alpine areas, characterized by a short and intensive snow-free season. Changes in the timing of snowmelt in spring may influence the timing of the growing season and thereby the growth, seed production and competition between plants, and is therefore the first indication of changes in plant diversity and vegetation cover (Wipf 2010, Wipf and Rixen 2010, Cooper et al 2011. Changes in snow season may influence the migration patterns and survival rate for birds, mammals, and insects at high northern latitudes (Kausrud et al 2008) including reindeer husbandry (Riseth et al 2011). Moreover, unseasonal snow thaw caused by midwinter warming events often have severe impacts on plants and animals (Bjerke et al 2014, Hansen et al 2014. Bokhorst et al (2016) provide a recent interdisciplinary review on advances in snow monitoring and modeling, and the impact of snow changes on ecosystems and society in Arctic regions.
Arctic and boreal temperatures are increasing due to positive feedbacks (e.g., albedo-temperature feedback) of increased greenness and shrubification (Xu et al 2013) and this is also linked to a shortening of the winter season. Hence, the exact timing of the onset and disappearance of snow is a critical parameter within a range of climatic and ecological studies. A series of methods have been developed to detect the start and end of the snow season. Hall et al (2002) used the normalized differential snow index (NDSI) based on Landsat and MODerate Resolution Imaging Spectroradiometer (MODIS) data, while Xu et al (2013) used daily records of freeze-thaw data obtained from passive-microwave remote sensing using the scanning multichannel microwave radiometer and special sensor microwave imager sensors with a spatial resolution of 25 km×25 km in order to monitor the snow-and ice-free season (photosynthetic active period). Liston and Hiemstra (2011) studied changes in the Pan-Arctic snow cover in the period 1979-2009 using a 10 km dataset based mainly on atmospheric reanalysis data, and found that the number of snow days per year decreased on average by 2.5 days per decade. This trend is supported by many similar studies referred to in SWIPA (2011). Dietz et al (2012) presented an inventory of snow cover characteristics based on MODIS Terra and Aqua for Europe for the period 2000-2011.
Within the study region (see section 2.1) there exists two operational services that provide regional snow cover data based on satellite and snow services. These are based on hydrological models provided by the hydrological authorities in Norway, Sweden and Finland (see e.g. senorge.no in Saloranta 2012). The satellite-based services are provided by SYKE (The Finnish Environment Institute) in Finland and Globesar in Norway. The SYKE snow map (Metsämäki et al 2012) is based on Terra MODIS data whereas the Globesar service is based on a multi-sensor snow map algorithm using Terra MODIS and Sentinel-1 data (Malnes et al 2005, Solberg et al 2010. This service was developed by Norut and the Norwegian Computer Center in Norway. Common to both these products is that they only provide data for the melting season (1 April-31 July) and focus on near-real time data for hydrological applications. Due to the near-real time approach, it is much harder to provide cloud gap-filled data sets than when data can be processed afterwards, as is the approach in this study. The main objective of this study is therefore to develop a method for detecting the first and last snow-free day geographically. To assess the accuracy of the method, we compare the outputs with the first and last snow-free days as measured at all meteorological stations within the study area. Finally, we discuss some applications of the dataset.

Study area
The study area covers the two northernmost counties in Norway, Finnmark and Troms, and nearby areas in Sweden, Finland, and Russia. It extends from approximate 68.1°N, 16.2°E to 70.8°N, 31.8°E (figure 1). A mountain range divides the area into two main climatic regions; (1) a western oceanic part, characterized by high precipitation rates and mild winters and (2) an eastern part characterized by drier and colder winters. The lowland in the western part is dominated by deciduous birch forest whereas in the southeastern part, lowland areas are covered by coniferous forest. However, the forest line is mostly below 500 m in the south and can reach sea level in the north, leaving large areas with arctic-alpine treeless tundra, within the mountain range and in the north (Moen 1999).
The vegetation in the study area is sparse. The dominating tree species along the coast are sparse birch (Betula pubescens) forests which are among the northernmost forests in the world (Wielgolaski and Sonesson 2001). The inner parts of Troms and central northern Finnmark county (Norway) host the world's northernmost Scots pine forests (Wielgolaski and Sonesson 2001). In the north-eastern part, arctic climate conditions (period: 1961-1990), prevail with an average temperature for the warmest summer month July of <10°C (Tveito et al 2000). The lowest monthly winter temperatures are found in the inland part of Finnmark (−16°C). Precipitation ranges from around 300 to 500 mm year (Tveito et al 2001) in Finnmark and Finland with rather low monthly winter precipitation (15-20 mm) in the continental parts of Finnmark and Finland (total winter precipitation of 87 mm). On the coast in Troms county, however, the monthly winter precipitation can be more than 100 mm. The snow cover is deepest around mid-March, with an average snow depth of 60-90 cm in eastern and northern Finland and at the coast of Troms County. On an annual basis, and especially during winter, inland Finnmark is the coldest and driest area amongst the Nordic countries (Tveito et al 2000).

MODIS data processing
This study uses the NASA MODIS/Terra Snow Cover Daily L3 Global 500 m Grid, Version 5 (MOD10A1) for the 2000-2010 period. The product provides snow cover fraction (SCF in the range 0 to 100%) as daily products (Hall et al 2002, Riggs et al 2006. SCF is calculated from the NDSI, from MODIS bands 4 (545-565 nm) and 6 (1628-1652 nm) using the formula, where NDSI band 4 b and 6 band 4 band 6 . 1 The product masks out areas that are covered by clouds, ocean and winter darkness. We converted the product to WGS 1984 UTM Zone 33 N and performed a refinement by applying multi-temporal interpolation techniques to obtain a cloud-free daily time series of snow maps for the study area. In total we have processed daily data sets from 24 February, 2000 (DOY 55) to 1 November, 2010 (DOY 304). With the exceptions of a few days (particularly in 2001, DOY 166-183) we have daily datasets for the whole period. For the period of polar night in the study area from 15-22 November to 20-30 January, we have no information and pixels corresponding to this period are flagged in the MODIS-products. We have, however, experienced that the products probably underestimate the snow cover fraction during the period before the onset of polar night (1 November-20 November) and during the period after the Sun returns (21 January-21 February). All snow maps have thus been 'padded' with 100% snow cover fraction in this period (21 November-20 February) to avoid confusion and create more likely results.
We have used a land mask for the region based on the N50 standard Norwegian map (www.kartverket. no) and maps of similar quality from Sweden and Finland to better delineate ocean, lakes, rivers and glaciers than can be done using the masks provided by NASA. This dataset has been further processed with multitemporal interpolation to eliminate clouds (SI: 1.1 Multitemporal interpolation method). For each day in the period we made an SCF map by mosaicking, reprojecting and resampling the MOD10A1 data to the fixed grid in UTM WGS z33N projection (see SI: figure S3).
Gap-filling techniques are frequently used in remote sensing data to mitigate lack of data due to clouds or other reasons (darkness, lack of coverage and instrument failure). Several authors (e.g. Høgda et al 2013) use the best observations in 7 or 8 day periods to find the most reliable data points. This strategy significantly alleviates the problems due to cloud cover significantly, but at the expense of a coarser temporal resolution. Our strategy maximizes the number of observations, but increases the risk for including noisy samples in the dataset. Such noise can for example occur in the vicinity of clouds. Since our study has a focus on detecting the first and the last snow-free day, we are confident that our approach reduces the sampling bias and increases the temporal accuracy for detection of the first and last snow free days.
Hall et al (2010) reviewed several cloud mitigation strategies for the MODIS snow products (see references therein) and also developed their own approach. Several approaches use both MODIS Terra and Aqua to provide an increased number of daily views, in combination with using temporal filtering/interpolations of observations from the days before. Spatial filtering has also been applied. Foppa and Seiz (2012) and Hüsler et al (2014) applied forward/backward gap-filling techniques for MODIS and AVHHR snow maps over the Alps, respectively. They used the average of the closest observation before and after the gap, instead of using a weighted mean (depending on the number of days before and after) as we do in this paper. The advantage of their approach is that it might be more tolerant to noise in the closest observation, in particular during stable snow conditions. The disadvantage is that the temporal information is smeared and the uncertainty in estimates of the transition time between snow-covered and snow-free conditions will increase.

Detection of the first and last snow free day
To select the best threshold for detecting the first and last snow-free day we use snow depth data from meteorological stations (see figure 1(b)), downloaded from www.eklima.no. Due to changes from manual to automatic measurements at some of the stations, and that the automatic stations often do not provide snow depth data, the number of stations with snow data were reduced from 40 in 2000 to 33 in 2010.
The mean SCF values were calculated for a 3×3 pixel area centered on each of the meteorological stations. The day when SCF passes below or above the threshold set at SCF=50% for an extended period (more than 10 days in the spring and more than 5 days in the fall) is set as the first and last snow-free day respectively. We also require that this can only happen in the period from DOY 75 (15 March) to DOY 230 (15 August) for the first snow-free day and in the period from DOY 264 (20 September) to DOY 310 (5 November) for the last snow-free day. These days are chosen to avoid problems due to low solar angles (in March and November) and confusion with random, short-lived snow falls during summer/early autumn. The detection rules are summarized in table 1.
In addition to the rules mentioned above, we also decided not to make a detection of first/last snow-free day for a given pixel if the age of the day for first/last snow free is above the threshold set at 15 days. It was found experimentally in the validation section (section 3.4) that a higher threshold led to increased bias/standard deviation in the first/last snow estimates. Likewise, a lower threshold led to fewer years (in particular autumns), where we were able to estimate the first/last snow.
The rules for detecting the first snow-free day and first snow-covered day in the meteorological dataset is very similar to the rules for satellite data, and was adapted from Vikhamar-Schuler and Hanssen-Bauer (2010). The only exception is that we now set the threshold to 2 cm snow depth. We have found that by systematically varying the threshold that 2 cm results in a minimum for the overall bias in the estimates for first and last snow cover in the two datasets.
These methods are illustrated for the meteorological station Kautokeino in figure 2 and S2. In figure 2, we show a comparison between SCF and snow depth at the Kautokeino meteorological station. Based on a total number of 43 cloud-free SCF detections during the year 2004, we were able to estimate the first and the last snow-free day with 1-2 days bias in comparison to the snow depth. The middle panel shows age, i.e. the number of days between the estimated SCF and an actual cloud-free measurement of SCF and this is used as an error estimate throughout this paper.

Development of daily snow maps
For each day in the period we made an SCF map by mosaicking, re-projecting and resampling the MOD10A1 data to the fixed grid in UTM WGS z33N projection (see figure S3(a)). In general, the maps have a high cloud fraction (see figure S3(a)). By applying the multi-temporal interpolation method discussed in 2.4, we obtained cloud-free snow maps for the study area. The method described above is now applied to the complete dataset consisting of 4015 (11×365) snow maps. For further information about development of daily snow maps see also SI, section 1.3.
2.5. Annual maps for first and last day of snow cover and the length of the snow-free season Based on the methods mentioned above and the multi-temporal MODIS dataset we are now able to provide maps for the first snow-free day in the spring and the last snow-free day in the autumn. In figure 3 we show an example for 2006. Each pixel in the maps is coded with the first/last day of snow. We can also calculate the pixel age, or uncertainty in this estimate of first/last snow cover. Finally, we have also calculated the number of snow-free days per snow season and the number of cloud-free days per year.

Snow season trends and variability
Snow seasons varied greatly between different years, both with regards to start, end and duration (figure 4). As expected, the alpine areas had much longer seasons than the lowlands. For instance the first snow-free day was delayed for the years 2000, 2005, 2007 and 2008. In contrast, the first snow-free day appeared to arrive  is shortest in the mountains and on the Varanger Peninsula (see figure 1(b)) and longest in the forested continental regions (south-east) and along the coast.
3.2. Average maps for first and last day of snow cover and the length of the snow free season After obtaining the first and last day of snow cover for each year, we can also find the average for the first and last day of snow cover. These estimates are shown in figure 5 together with the overall estimate for the length of the snow-free season. As we can observe from figure 5, the coastal parts and the river valleys as well as the more forested areas in Sweden and Finland experience a quite early start of the snow-free period (spring) as well as late start of the snow season in the autumn/winter. Hence, the snow-free period is longest in these parts of the study area.

Annual variability in the snow cover
We have applied the results from sections 3.3 and 3.4 to study the annual variability in the snow cover for the entire study area and the geographical variability per year. In figure 6, we show the total daily SCF per day for the entire 11 year time series. We also show the mean value for all the years. By comparing the year to year variation we can observe that the start of the snow-free periods  Table 2 shows the relationship between the first and last day of snow cover measured from MODIS data and from the meteorological stations. The number of years with data is less than 11 for most of the stations due to either lack of MODIS data (due to clouds/polar darkness) or lack of snow depth data from the meteorological stations. In spring, the correlation (r=0.79) between the MODIS-based and the meteorological measurements of the last day of snow cover was very high, where 31 of the 40 stations show significant values at the 1% level or better. Only four of the stations did not show significant correlation (p>0.05). The average standard deviation is six days. Only six of stations show a bias of 10 days or more, where the MODIS-based measurements were earlier than the meteorological measurements in these cases. The relationship between MODIS-based and meteorological measurements of the first day with snow in autumn was weaker (r=0.51) compared with the spring period. On average only 6 years of data were available for the stations. However, 27 of the stations show significant values at the 1% level or better, and the average standard deviation is 7 days. Three of the stations have a bias of 10 days or more and for all of these stations the MODIS-based measurements of the first day of snow were earlier than the meteorological measurements. Only 2 of 40 stations have a correlation of less than 0.5 for the first snow-free day. The standard deviations and bias were also higher for the end of the growing season than for the spring. The reason for this can be found in the number of years applied. On average, we use only 7 of the 11 available years to estimate the last snow-free day. The reason for this is mostly due to late first snowfall in the autumn, which yields a lower probability of detecting the transition. Also, some stations lack snow depth measurements for some of the years. Some of the stations are located in challenging terrain as the stations are often localized close to the lowest point in the 3×3 pixel area used for estimating SCF or in coastal areas with high cloud cover probability (e.g Gamvik, Ytre Holmbukt, Sørkjosen). The Sørkjosen station (r=0.25) also has a high standard deviation and an inspection of this site indicates that the spring of 2006 and 2009 were exceptions, leading to large errors. The snow depth record from 2006 is  Table 2 also shows that the correlation between the MODIS estimate and the snow depth estimate of the last snow-free day is significantly poorer (r=0.65).

Assessment of the MODIS based snow cover products
The average values indicate that the first snow-free day for the entire area is approximately 1 May and that the last snow-free day is around 15 October. Due to the bias in the distribution of meteorological stations, these values are applicable for the coastline.

Critical discussion on methods and results
We are confident that the methods presented here give reasonable results for the estimation of first and last snow-free day. The time sequences of cloud-free (multi-temporal) snow cover maps as well as the annual maps for first and last snow-free day are reasonable at the scale we are studying (500 m resolution). The quality of the MOD10A1 product is questionable (e.g. Metsämäki et al 2015). However, the overall absolute clear-sky accuracy of the well-studied 500 m resolution swath (MOD10_L2) and daily tile (MOD10A1) products was assessed to be >93%, but varies with land-cover type and snow condition (Hall and Riggs 2007). The cloud mask in the MOD10A1 product is probably too conservative most of the time, and we observed frequently that the cloud cover was overestimated. Since the MOD10A1-algorithm was optimized in order to perform on a global scale, this is  expected. Other algorithms have improved the performance for forested areas Rune Solberg 2003, Metsämäki et al 2012). We also observed occasionally falsely detected snow during summer (related to erroneous cloud classification) and falsely detected areas with low snow cover fraction during the winter period (related to poor light conditions and uncompensated topography).
The method we have developed for detecting first and last day of snow cover could also be applied on improved SCF-time series (e.g. Globesnow2-data, Metsämäki et al 2015), but we suspect that the results will only improve slightly since our method is robust against many of the existing problems in the MOD10A1 products. Problems with conservative cloud estimates are handled by the multi-temporal interpolation scheme. Effects associated with falsely detected snow during summer and too low SCF in mid-winter will seldom be a problem as we have developed detection rules that excludes such detections in the majority of cases.
The linear interpolation approach may also be questioned. It is a well-known fact that the snow depletion curve is seldom linear in studies where mid- Unmarked correlation values have significance at the 1% level, * significance at 5% level, ** insignificant result.
to large-scale hydropower basins with a large altitude gradient are studied (Kolberg and Gottschalk 2010). On the scale we study however, the snow depletion curve will in most cases approach a sigmoid function with a rather steep slope. The steeper the slope, the more important it becomes to have cloud free SCF estimates around the timing of the snow melt. This problem could be alleviated by including SCF estimates from synthetic aperture radars (SAR) in the multi-temporal time series. Malnes et al (2005) have developed methods to fuse SAR and optical time series, but the density of SAR time series is seldom as consistent as optical datasets. We experienced in our study (table 2) that the average correlation between the MODIS estimate and the snow depth estimate of the last snow-free day was significantly poorer (r=0.51) than the correlation for the first snow-free day (r=0.79). We do no attempt to extract trends from the current dataset. An 11 year dataset is too short, and it can easily be seen from e.g. Figure 4 that the individual years with alternating dominating regional patterns of early and late snow cover will dominate any trends at the current stage. Future studies will attempt to use the whole time series of MODIS data (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) and might provide some more clarity concerning observable trends.
The validation provided in table 2 gives overall trustworthy values for correlation and significance. We may thus state that our methods developed to assess the first and the last dates of snow cover are overall trustworthy. It should however, be emphasized that the distribution of meteorological stations used for the validation are strongly biased towards the coastline and low altitudes.
A weak point in our validation assessment is the method we use for removing outliers. Since MOD10A1 data gives little to no information following the onset of polar darkness, we are forced to remove all winters with late onset of snow. This yields a clear bias towards the years when the onset of snow is early, and as long as we do not have better sensors to capture late snowfalls this will be a clear limitation. However, we believe that the information provided is better than no information, and we thus choose to provide it.

Ecological interpretation
The averaged map (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) for the length of the snow-free season can be compared with phenological maps showing the extent of the growing season. The length of the snow-free season is a good proxy for the growing season  and the photosynthetic active period (Xu et al 2013). We find similarities between the map in figure 5 and the maps produced by ), Høgda et al (2013 and Xu et al (2013). However, the start and the end of the growing season cannot be compared exactly since vegetation activity is included in the phenological maps (Høgda et al 2013). In our approach, we detect the start and end of the snow season, hence there can be several days/weeks difference from when the snow-free period starts to when vegetation green-up starts (e.g. bud burst of downy birch). Similarly in autumn, there may be several days/weeks from when the vegetation becomes yellow/brown or leaves have fallen until the start of the snow season. However, it is clear that the map of the duration of the snow-free season (figure 5) has provided valuable information in order to detect the extent of the growing season as well as data for assessing and modeling the reindeer range conditions In addition to differences in the snow-free season due to land cover (vegetation type), we also notice that height gradient and continentality are major factors affecting the length of the snow-free season and this is consistent with . We find that the longest lasting snow-free areas are located along the coast and in the river valleys as well as in the areas around the Inari lake (Finland), while the shortest snow-free periods are found in the mountains and in the north-eastern parts of Finnmark (Varanger peninsula), and here we also find the sparsest vegetation .

Hydrological interpretation
Daily snow maps have been used in assimilation with lumped hydrological models (Engeset et al 2003) with limited success, mainly due to the problems with overcalibration against run-off. The authors emphasize the importance of timely snow maps in the melting season when snow reservoirs are uncertain. In particular flood warnings can be issued in such a context. When used in conjunction with available models, snow cover data based on MODIS measurements may improve the runoff forecast significantly when assimilated in distributed hydrological models for Snake River basin in USA (Andrealis and Lettenmaier 2006).

Climate interpretation and other application areas
The current study provides detailed maps for the length and variability of the snow season in the study area. Similar studies using MODIS have also provided maps at a coarser scale for the northern hemisphere (Chen et al 2015) and Europe (Dietz et al 2012). The current dataset is not sufficiently long to derive trends, but the spatial resolution is sufficiently high to show many regional details which can be related to the complex interaction between large scale weather patterns affecting the coastal climate in northern Norway, the topography and the continental sub-arctic climate zones in the interior of the land. The overall picture, where the length of the snow season is one of several factors that determine the living conditions for many plants, show remarkable resemblance to vegetation maps for the same region. Regions of similar snow phenology also have similar vegetation covers. This shows that the maps, when updated on a regular basis, could be used to determine where there are significant and long-lasting deviations from the average that also could indicate ongoing climate changes.

Conclusions
We have described a method for generating long-term, cloud-free daily time series of the snow cover fraction based on optical medium resolution sensors. The method is applied to the MOD10A1 dataset to provide an 11 year record for our study area in northern Norway. The time series of snow maps shows a convincing development of the snow cover that can be used for generating climatological records as well as in ecological monitoring and assessments. We have calculated annual maps for the first and last snow-free day for each year in the period. The maps give reasonable estimates which are supported by validation against data from 40 meteorological stations. We have also calculated the average first and last snow-free day for the 11 year period, and the average length of the snow free period.
The correlation between the first snow-free day mapped by MODIS data and by snow data from 40 meteorological stations was highly significant (p<0.05) for 36 of the stations, and with a bias of less than 10 days for 34 of the stations. The correlation for the last snow-free day for 31 of the stations was significant (p<0.05), and the bias was less than 10 days for 27 of the stations. The region of interest is located at high latitudes and is affected by low solar angles/ polar darkness during the first snowfall in autumn. These conditions are highly noticeable in the dataset, and cannot be totally removed by any algorithm.
When the first snowfall arrives later than onset of polar darkness (1-20 November), it will remain undetected until we acquire new observations in January-February. Our approach places an emphasis on maximizing the number of observations, and provides in our opinion, the best compromise between using only real observations and refraining from guesswork.
The map showing the average length of the snowfree period is consistent with the land cover in the region, and reflects differences due to continentality and height gradients. By inspecting the map manually and comparing it with local observations (meteorology) we find very convincing trends. Some of the regions with earliest snow-free conditions are well known in local municipalities. E.g. the southern part of Tromsø island becomes snow-free earlier than in the northern part. Similarly we also find early snowfree conditions in farmland areas.