Impacts of Snow Cover on Vegetation Phenology in the Arctic from Satellite Data

Remarkable shifts in key parameters of vegetation phenology have been detected from satellite and field observations in the past three decades over various regions and biomes, especially in northern high latitudes (Menzel and Fabian, 1999; Chmielewski and Rotzer, 2001; Ho et al., 2006; de Beurs and Henebry, 2008; Jia et al., 2009; White et al., 2009; Zeng et al., 2011; Jeong et al., 2012). Such changes in vegetation seasonality may have important consequences for regional and global carbon budgets and climate feedback. Most phenological studies based on satellite data have focused on the phenological changes that occurred before 1999 (Myneni et al., 1997; Zeng et al., 2011), but different trends in climate change might have taken place since 1999 (Buermann et al., 2007; Cane, 2010). Moreover, both an earlier start and a later termination resulted in a significant lengthening of the growing season during the 1980s and 1990s, but since 2000 different trends in phenological shifts have been gradually reported (Jeong et al., 2011; Piao et al., 2011). Thus, given that phenology is a preferred indictor of climate change, it is essential to examine those phenological shifts that might have occurred during the most recent decade (2000-10).

Temperature shifts have widely been reported as a critical factor driving changing vegetation phenology at global and regional scales (Menzel et al., 2003; Gordo et al., 2009; Zeng et al., 2011; Jeong et al., 2013), but knowledge about the response of vegetation phenology to the dynamics of snow cover is still limited. Snow cover has changed significantly in response to warming in the Northern Hemisphere (NH) over past decades (Brown et al., 2000; IPCC, 2007), with a continuous decrease in the length of the snow season in the north (Raisanen, 2008). The spatial extent of snow cover, especially in the spring, has decreased significantly during recent decades in NH land areas, as seen from satellite and field observations (Brown et al., 2000; Easterling et al., 2000). In addition, it has been reported that snow cover duration decreased by 1.8±1.7 days per decade in tundra ecosystems of northern high latitude land areas from 1988 to 2002 (Smith et al., 2004). The shifts in the extent and timing of terrestrial snow cover may be important factors resulting in changes to vegetation phenology for several reasons: first, snow cover is closely related to surface temperature (Jonas et al., 2008) and soil moisture (Manabe and Wetherald, 1987) in the early growing season of northern vegetation, and therefore strongly affects onset timing (Menzel et al., 2006). Second, variations in snow cover in high latitudes have a more important role in energy feedback via albedo, surface heating and water budgets (Serreze et al., 2002; Yang et al., 2003), and snow-vegetation feedback is considered a key mechanism behind snow-driven phenology change (Groisman et al., 1994b). Therefore, understanding how vegetation phenology responds to variations in snow cover is a critical topic for projecting future ecosystem dynamics.

Decreases in snowpack during the last 30-50 years and the associated advance in "start of season" (SOS) by 1-4 weeks have been reported over the western United States and Canada (Stewart et al., 2005), and the early thaw of snow cover has been attributed to a longer length of season (LOS) ( ila-Jimàez and Coulson, 2011). In high latitude areas, snow cover duration largely determines the length of the growing season for vegetation (Oberbauer et al., 1998). However, there still exist uncertainties and gaps in our understanding regarding the relationships between changes in snow cover and vegetation phenology, for instance in terms of how the timing of snow cover influences vegetation seasonality in northern high latitudes. To gain better insight into snow-phenology interactions over northern high latitudes, we selected the Yamal Peninsula in Arctic Russia as a case study area for which recent field observation and satellite data could be used for analysis. The Yamal region is considered an Arctic tundra "hot spot" because of the rapid changes in climate and to the ecosystem that have occurred there (Forbes et al., 2009; Walker et al., 2009; Forbes et al., 2010).

The main objectives of our study were: (1) to quantify general trends and spatial heterogeneity of vegetation phenology and snow cover in the most recent decade (2000-10); (2) to determine possible shifts in vegetation phenology and the timing of snow cover during the period; (3) to assess the relationships between the timing of snow cover and phenology among different tundra vegetation types and latitudinal gradients in the region.

Circumpolar Arctic Vegetation Map (CAVM) (Walker et al., 2005) datasets were used to categorize vegetation classes in the case study areas as follows: graminoid tundra (GT), prostrate-shrub tundra (PT), erect-shrub tundra (ET), wetlands, and carbonate mountain complex (CMC). About 31.34% of the study region was covered by GT, 4.32% by PT, 45.5% by ET, 15.2% by Wetlands, and 3.64% by CMC (Table 1). GT and ET covered nearly 80% of the areas in the study region. Fractional land cover was spatially summarized at 0.5° latitude intervals along the Yamal Peninsula based on CAVM data at a scale of approximately 1:7 500 000 (Fig. 1). More details about the study area and each site can be found in the 2007 and 2008 Yamal expedition reports and some other previous studies (Walker et al., 2005; Walker et al., 2009; Goetz et al., 2011).

Figure 1.  Fractional land cover of vegetation with 0.5° latitude intervals along the latitudinal gradients.

The Normalized Difference Vegetation Index (NDVI) (Tucker, 1979) calculated by near-infrared and red visible reflectance is an indicator of vegetation greenness (Myneni et al., 1997) and is widely used to monitor vegetation phenology in the north (Jia et al., 2006; Cleland et al., 2007; Jia et al., 2009). NDVI data spanning from 2000 to 2010 at a spatial resolution of 500 m and over eight-day intervals were derived from a pair of moderate resolution imaging spectroradiometer (MODIS) products. The first of these two products (MOD09A1) provides data in bands 1-7 and NDVI is calculated from the red band ( B₁: 0.62-0.67 m) and the NIR band ( B₂: 0.841-0.876 m) using the following equation: NDVI=(B₂-B₁)/(B₂+B₁). (1)

Pixels with a maximum NDVI in a year less than 0.09 were masked off to minimize the bare and sparsely vegetated regions in our study area. The HANTS (Harmonic Analysis of NDVI Time Series) algorithm was applied to the NDVI data to reduce the possible noise existing in the MODIS-NDVI profile because of cloud, aerosol and snow cover. TIMESAT was used to explore and extract the phenology parameters, including SOS and "end of season" (EOS) over the region during the period 2000-10. There is not much difference between "double logistic" and "Savitsky-Golay" as effective smoothing methods for NDVI time series (Jonsson et al., 2010), and so the "double logistic" method was chosen in TIMESAT to smooth the NDVI time series in our study. SOS and EOS were defined as the periods with the fastest NDVI increase in spring and decrease in autumn, and therefore the maximum and minimum values of variability of NDVI were set as the thresholds (Piao et al., 2006; Jeong et al., 2011).

The second of the two MODIS products adopted in this study (MOD10A2; MODIS snow cover eight-day L3 global 500 m grid) was used for the determination of maximum snow cover extent (Hall and Riggs, 2007). The Normalized Difference Snow Index (NDSI) has been widely used to determine effective snow cover mapping (Hall et al., 2002; Nagler et al., 2008), which is based on MODIS bands 4 ( G₁: 0.545-0.565 m) and 6 ( S₆: 1.628-1.652 m) as follows: NDSI=(G₁-S₆)/(G₁=S₆). (2)

When the green reflectance, NIR and NDSI are greater than the thresholds of 0.1, 0.11 and 0.4, respectively, a pixel is indicated as snow. MODIS snow cover data were widely selected to define "last day of snow cover" (LSC) in the spring and "first day of snow cover" (FSC) in the fall by defining 20% or 30% as the snow cover extent percentage thresholds (Reed et al., 2009; de Beurs et al., 2010). It was possible to apply this algorithm to calculate LSC and FSC for the whole study area. To obtain accurate LSC and FSC in each pixel, the LSC date was determined as when each pixel was identified to be snow-free for at least three continuous days during the period 1 April to 31 August (Brown et al., 2010). FSC was determined when the snow cover fraction increased rapidly from 0% to 100%. Cloud-covered pixels were considered as the main limitation for the usage of MODIS snow cover data, and thus we applied a cloud-removal algorithm(Gafurov and Bardossy, 2009)to remove the cloud pixels.

Least-squares linear regressions were applied in statistical analyses to examine the trends of phenology during a given period. One-way analysis of variance (ANOVA) was used to test for significance in correlations between vegetation phenology and snow cover timing. Correlations of |r|>0.6 were highlighted in our spatial datasets.

Figure 2.  Fig. 2.11-yr averaged (2000-10) spatial distributions of climatology of vegetation for (a) SOS and (b) LSC, and its linear trends at 0.5° latitude intervals along the latitudinal gradients based on MODIS data during the period 2000-10.

Figure 3.  Fig. 3. 11-yr averaged (2000-10) spatial distributions of climatology of vegetation for (a) EOS and (b) FSC, and its linear trends at 0.5° latitude intervals along the latitudinal gradients based on MODIS data during the period 2000-10.

Figure 4.  Fig. 4. Spatial patterns of (a) SOS and (b) LSC linear trends, and their area-averaged values at 0.5° latitude intervals along the latitudinal gradients from 2000-10.

Figure 5.  Spatial patterns of (a) EOS and (b) FSC linear trends, and their area-averaged values at 0.5° latitude intervals along the latitudinal gradients from 2000-10.

The decadal mean spatial distributions of the dates of phenology and snow cover during 2000-10 were analyzed (Figs. 2 and 3). Both SOS and LSC ranged widely from approximately day 120 to day 180 in the year (hereafter "day of year"; DOY) (Fig. 2). Vegetation turned green and snow melted much earlier in the south, with SOS and LSC respectively shifting to almost 4.1 days (r²=0.982, p<0.0001) later and 3 days (r²=0.984, p<0.0001) later with every one degree northward change in latitude. SOS was several days earlier (less than ten days) than LSC in the south, with the two dates becoming almost the same from around 71°N where the dominant vegetation turned from ET to GT (Fig. 2). The average EOS and FSC dates ranged from approximately 265 to 285 DOY in the last decade. A general trend of 1.9 days earlier for EOS (r²=0.938, p<0.0001) and 1.2 days earlier for FSC (r²=0.48, p<0.0042) was found with every one degree northward change in latitude (Fig. 3). EOS was several days later than FSC across almost the entire region, except at latitudes of 71°-72.5°N. However, spatial patterns of phenology and the timing of snow cover were heterogeneous even at the same latitude mainly due to diverse regional features. For example, comparing the same latitude, EOS and FSC were earlier in the region where CMC dominated.

To examine the interannual variations in phenology related to shifts in the timing of snow cover during the last decade, the spatial distributions of the linear trends of phenology and snow cover period were calculated and summarized (Figs. 4 and 5). Across the region during the study period (2000-10), the number of days by which SOS and LSC shifted ranged from an advance of about 2 days to a delay of more than 3 days, depending on latitude (Fig. 4). There were clearly different patterns of changes in SOS and LSC. For instance, the region from 70°N to 73.5°N experienced greater variations in LSC than variations in SOS. SOS changed by 1.0 to 2.0 days per year in the region from 66_°N to 69.5°N, indicating a stronger delay to SOS in the region. Along the lines of latitude from south to the north, there was an obvious delayed SOS in the south and advanced SOS in the north during the last decade, with the turning point being at about 71°N north. Almost the entire region experienced a delayed LSC, apart from between 72°N and 74°N. The value of changes in SOS (r²=0.894, p<0.0001) and LSC (r²=0.692, p<0.0001) increased significantly with decreasing latitude during 2000-10 along the latitudinal gradient.

Changes in EOS and FSC were from less than -2 days per year to more than 2 days per year in the last decade across the region (Fig. 5). EOS in most areas showed a strong agreement with FSC, except the most northern part of the region (72.5°-73.5°N). There was a significant relationship between latitude and EOS (r²=0.681, p<0.0001), while a weaker relationship was found between latitude and FSC (r²=0.302, p<0.05). Both EOS and FSC in the region between 65.5°N to 68°N, and between 70°N and 73°N advanced significantly, while a different trend existed in the other parts of the regions.

In order to understand pre-growing season shifts in snow cover related to phenological variations, correlation coefficients were calculated from 2000 to 2010 (Fig. 6). Large areas in south and central parts showed strong correlation coefficients (r² > 0.6) between SOS and LSC. However, there were no clear latitudinal gradients in the spatial distributions of strong correlation coefficients between EOS and FSC across most of the region. To discover any meaningful relationships between the timing of snow cover and the associated phenology responses among vegetation types, we examined the response of SOS to LSC and the response of EOS to FSC among the five dominant vegetation types (Fig. 7 and Table 1). There were significant relationships between SOS and LSC in the areas where GT (r²=0.78,p<0.001), ET (r²=0.82, p<0.0001), wetlands (r²=0.84,p<0.0001), and CMC (r²=0.93,p<0.0001) dominated. However, the relationship was slightly weaker in areas corresponding to high fractional cover of PT (r²=0.44,p<0.05). The correlation coefficient (k) for SOS responding to snow melting was a little more than one in the areas where ET, wetlands and CMC dominated, indicating a slightly delayed SOS (Table 2). LSC was greater than SOS over the entire region except in PT-dominated areas. LSC and SOS associated with the region where ET and CMC dominated showed a significant advanced date, compared to those in areas where GT and PT dominated. SOS of these various land cover features responded to snow cover differently. In the region where ET, wetlands and CMC dominated, k was larger than one, with a strong and significant relationship between SOS and LSC. LSC happened earlier than SOS for less than five days (Table 2). The response of EOS to snow cover was almost the same in each vegetation type, with k coming out as almost one, a strong and significant relationship. Both SOS and LSC in areas dominated by ET, wetlands and CMC were significantly earlier (95% confident level) than for other types. The five vegetation types showed no significant differences in FSC. The difference (LSCþÿþÿSOS) was larger in the region where wetlands and CMC dominated, and it showed that LSC was 4 days earlier than SOS. EOS was later than FSC for almost all types in the areas, while quite a different trend was found in areas dominated by CMC. EOS was earer in the region where CMC dominated.

Figure 6.  Spatial distributions of statistical correlations between phenology and the timing of snow cover in the region from 2000-10 from MODIS products. Here, r > 0.6 was used to show a relatively strong positive correlation. (a) SOS vs. LSC; (b) EOS vs. FSC.

Key vegetation phenology parameters and related snow cover variables were evaluated from 2000-10 over the Yamal Peninsula (Figs. 2 and 3) based on satellite sensor (MODIS) data. Vegetation can be detected by such means after snow melt, and by doing so the snow-free period was found to be closely related to the length of the growing season. SOS was several days earlier than LSC from 65.5°N to 70.0°N, which meant vegetation had turned green before snow cover disappeared in these areas. Farther north along the latitudinal gradient, from 70.0°N to 73.5°N, SOS was almost the same as LSC. Thus, we used this method to determine that SOS is less affected by snow cover in the south of the region studied (65.5°-70°N), where the vegetation is dominated by ET, with a much higher NDVI and level of vegetation production (Fig. 1). EOS was found to have a similar date as FSC in most study areas, apart from southern parts of the region (Fig. 2). NDVI did not show a sharp decrease in the fall, and snowfall could have influenced the vegetation indices, which identified a "sudden change" in the time series of NDVI. That is, as a consequence of snowfall in the fall, vegetation may have been suddenly covered by snow, and thus EOS would have been detected by the sensors. However, snowfall had a weaker influence on determining EOS in the south of the region, where taller and denser vegetation types, which are less sensitive to occasional snowfall, are dominant. We found delayed EOS in some areas, compared to FSC. Thus, FSC can be used to predict EOS in those areas with low vegetation density and high sensitivity to snow cover. In addition, MOD10A2 maps the maximum snow cover extent during a period of 8 days, meaning some areas of snow cover would not have been identified as snow if snow fell at the end of an 8-day period, and also if the pixel was identified as cloud for the following compositing period (Hall et al., 2002). This could also have been the reason for differences between phenology and the timing of snow cover in some parts of the region.

Comparing the same latitude, SOS and LSC occurred later in the region where CMC dominated, because snow melts earlier in open areas than the surrounding terrain (Brown et al., 2003). Another possible explanation for later snowmelt relates to the higher elevation and associated lower temperatures of this mountainous area, as compared to other areas at the same latitude. Therefore, the snow pack remains for a longer period of time and a delayed SOS subsequently takes place. There was no evident relationship between FSC and latitude, implying a much more complex control over the distribution of FSC.

Recent studies have suggested an earlier LSC (Brown, 2000; Dye, 2002) and SOS (e.g. Jia et al., 2006) across many regions of the NH. However, a delayed LSC and SOS were found in the middle and south of the region in the present study. An earlier SOS may have been stalled due to a different tendency of temperature during the first decade of this century (Jeong et al., 2011; Piao et al., 2011). Additionally, in our previous study, temperature and anthropogenic factors were used to examine the possible reason for phenology in the region. The results indicated that temperature is a key factor controlling phenological spatial gradients and variability, while anthropogenic factors (reindeer husbandry and resource exploitation) might explain the delayed SOS in southern Yamal.

Reducing uncertainties in determining the annual timing of snow thaw and onset could be helpful for effectively detecting changes in climate (Ye, 2001) and ecosystems in the north (Running et al., 1999). Variations in snow cover are mainly driven by warmer air temperatures in the NH (Brown et al., 2010; Brown and Robinson, 2011; Huang et al., 2011). Meanwhile, snow cover variations play an important feedback role in climate change over the NH (Hall and Qu, 2006; Baringer et al., 2010; Flanner et al., 2011) through a snow-albedo mechanism (Saunders et al., 2003), affecting terrestrial ecosystems by modifying soil and air temperatures, soil water recharge, and the availability of photosynthetically-active radiation to vegetation (Gerland et al., 2000). A number of studies have suggested that snow may be an important controlling factor over vegetation seasonality and spatial distributions in northern biomes (Peterson and Peter-son, 2001; Litaor et al., 2008; Wipf and Rixen, 2010). In addition, variations in spring snow cover are considered valuable indicators of climate change due to its sensitivity to temperature, as well as variations in temperatures (Karl et al., 1993; Groisman et al., 1994a). The annual cycle of snow cover, including LSC in spring and FSC in fall, can be an essential factor for vegetation phenology through its influence on the timing and strength of surface energy and hydrological budgets (Pielke et al., 2000; Arora and Boer, 2001; Serreze et al., 2002; Yang et al., 2003; Zhang, 2005; Kreyling, 2010), especially in cold ecosystems (Chapin et al., 2005). In our study, we found strong correlations between phenology and the timing of snow cover over the case study area. The positive feedback from reduced snow cover was thought to be one of the reasons for enhanced spring warming (Groisman et al., 1994a), and the warmer temperature could cause earlier SOS. Meanwhile, delayed LSC was found in the middle and south of the study region, which may increase the surface albedo because the change in the state of the land cover from snow to bare soil or short grass takes occurs later. Subsequently, a decrease in temperature is found, which potentially contributes to the delayed SOS. Consequently, changes in both the extent and timing of snow cover could alter vegetation phenology. LSC may be controlled by more complex factors than FSC, and we also found greater temporal variability in LSC than FSC (Reed et al., 2009).

Responses of vegetation phenology to the latitudinal gradients of the timing of snow cover were also analyzed. Both SOS and EOS showed significant relationships with latitude, and there were significant relationship as well between SOS and LSC in most regions, implying LSC might be an important factor influencing latitudinal distributions of SOS.

Shifts in snow cover and their feedback to local and regional climate are considered particularly influential in the high latitude regions above 50°N, and there might be large differences among vegetation types through albedo and surface heating ((Betts and Ball,1997; Eugster et al., 2000). To gain a better understanding of the effects of snow cover on vegetation phenology, the spatial heterogeneity of phenology in response to variations in the snow cover dates was investigated along spatial gradients and among different land cover types. Vegetation phenology responding to shifts in the timing of snow cover was different among land surface types (Table 2), indicating influences may vary among tundra types because of different sensitivities to seasonal changes in snow cover. In winter, low-lying plants are buried by snow and relate to higher albedo values (0.85), while tall plants with exposed canopies relate to albedo values of 0.60 (Sturm et al., 2005). When spring arrives and snow begins to melt, areas dominated by low-lying plants can absorb more solar radiation during the snow-cover period, resulting in earlier melting and SOS. Owing to significant warming, taller shrubs and boreal forest may be partly replaced by short tundra vegetation in the southern Arctic (Sturm et al., 2001), possibly making these regions more sensitive to changes in snow-cover duration.

During the time of LSC and FSC, when snow depth is very low and cannot cover an entire pixel, the pixel will be defined as snow-covered or snow-free solely by the different reflectance values in the visible and NIR bands. Moreover, MODIS data have their own limitations in terms of mapping shallow snow depths (<4 cm) (Wang et al., 2008). In addition, land cover type is considered to have a greater impact than snow depth in the algorithm used by the MODIS snow-cover product when the snow depth is very shallow (Huang et al., 2011). Thus, we must consider the effect of land cover type on the accuracy of MODIS snow cover mapping as a factor. Different responses of the vegetation types to snow cover may decide the spatial heterogeneity of phenology. The lagged effect of snow cover changes on phenology depends primarily on vegetation type. Fractional land cover possibly contributed to changes in species cover and diversity of the vegetation, and these various land cover features respond to snow cover differently. Responses of different vegetation types to changes in snow cover varied along the latitudinal gradients. The largest change in SOS was found in the south (65.5°-70°N), where ET mainly dominates, while variations in EOS did not show a clear latitudinal distribution. Compared to other types, a CMC areas were shown to have a late EOS and early FSC, thus demonstrating a negative trend (FSC-EOS) and indicating vegetation turned brown earlier than snow began to fall.

Temperature is used widely as a main factor to explain phenological shifts over many biomes throughout the world. Consideration of the timing of snow cover may help to explain vegetation phenology in the northern high latitudes. Here, we demonstrated that the distribution of vegetation phenology and the timing of snow cover are closely linked. Shifts in LSC were found to be an important indicator of changes in temperature and showed strong agreement with SOS. Delayed LSC could contribute to delayed SOS through an albedo-temperature relationship. EOS was found to be determined mainly by "sudden" snowfall, and the warming associated with changes in snow cover is likely the dominant control for vegetation phenology shifts in the northern high latitudes. The effect of anthropogenic factors on vegetation phenology cannot be ignored. Therefore, future work should consider how to effectively remove the impact of anthropogenic factors and analyze the quantitative relationships between climate change and vegetation phenology. Here, linkages between the timing of snow cover and vegetation phenology were investigated by analyzing MODIS data and vegetation classification. Because of the extreme logistical problems involved in taking ground measurements, such as the high cost and manpower required, the MODIS snow cover product is critical, but it does possess limitations when mapping shallow snow depths (<4 cm). Therefore, despite the constraints, in situ observations in the region are crucial for validating remotely sensed estimations of snow cover.