Relationship between Spatiotemporal Variations of Climate, Snow Cover and Plant Phenology over the Alps—An Earth Observation-Based Analysis

Abstract

Alpine ecosystems are particularly sensitive to climate change, and therefore it is of significant interest to understand the relationships between phenology and its seasonal drivers in mountain areas. However, no alpine-wide assessment on the relationship between land surface phenology (LSP) patterns and its climatic drivers including snow exists. Here, an assessment of the influence of snow cover variations on vegetation phenology is presented, which is based on a 17-year time-series of MODIS data. From this data snow cover duration (SCD) and phenology metrics based on the Normalized Difference Vegetation Index (NDVI) have been extracted at 250 m resolution for the entire European Alps. The combined influence of additional climate drivers on phenology are shown on a regional scale for the Italian province of South Tyrol using reanalyzed climate data. The relationship between vegetation and snow metrics strongly depended on altitude. Temporal trends towards an earlier onset of vegetation growth, increasing monthly mean NDVI in spring and late summer, as well as shorter SCD were observed, but they were mostly non-significant and the magnitude of these tendencies differed by altitude. Significant negative correlations between monthly mean NDVI and SCD were observed for 15–55% of all vegetated pixels, especially from December to April and in altitudes from 1000–2000 m. On the regional scale of South Tyrol, the seasonality of NDVI and SCD achieved the highest share of correlating pixels above 1500 m, while at lower elevations mean temperature correlated best. Examining the combined effect of climate variables, for average altitude and exposition, SCD had the highest effect on NDVI, followed by mean temperature and radiation. The presented analysis allows to assess the spatiotemporal patterns of earth-observation based snow and vegetation metrics over the Alps, as well as to understand the relative importance of snow as phenological driver with respect to other climate variables.

1. Introduction

Phenology, the science of the timing of annual recurring biological events, has been studied for centuries. Within the field of vegetation phenology, these events comprise plant stages such as bud-burst, flowering, leaf unfolding, or leaf-fall. Monitoring and understanding plant phenology is important in the context of global change, because the timing and magnitude of phenological events are strongly sensitive to seasonal and inter-annual climate variability [1,2,3]. A key aspect of phenological studies is thus determining causal relationships between physiological plant phenomena and their seasonal and long-term drivers.

Phenological ground networks provide the longest sources of observations, which are often species specific, detailed, and highly accurate, but at the same time only conclusive for one organism at a certain location [4,5,6,7,8]. To derive spatially exhaustive information, satellite remote sensing data that trace seasonal changes in the spectral signature of vegetation photosynthetic activity have increasingly been used during the last decade [9,10]. Such remote sensing-based analyses are referred to as land surface phenology (LSP) [11,12,13,14,15,16,17,18]. LSP observations provide a spatially integrative view of continuous biophysical canopy properties at coarse scales instead of plant-specific phenological stages. Nevertheless, LSP metrics proofed to be comparable to these development stages, good measures of phenophases at the ecosystem level, and a suitable proxy indicator of climate variations [19,20].

Due to this sensitivity to climate variability, phenology has gained importance as an indicator for ecological responses to climate change [4,21,22,23,24,25]. Concurrently, phenology provides feedbacks to climate by influencing albedo, temperature and precipitation patterns [26,27,28,29]. A reliable representation of vegetation phenology in climate models is therefore necessary to model carbon, energy, and water cycles on a regional to global scale [27,30,31,32]. In this context, prognostic phenology models have been developed to simulate phenological events based on driving factors such as temperatures, photoperiod, precipitation, or soil moisture [33,34,35,36,37,38,39,40,41,42].

Mountain ecosystems are assumed particularly affected by climate change with various effects on ecosystems, cryosphere, and hydrological regimes [43,44,45,46,47,48,49,50]. Observed effects of climate warming on mountain phenology are e.g., longer growing seasons [19], the migration of plant species to higher altitudes [51,52,53,54] and the associated impacts on niches and endemic species [53,54,55,56,57,58]. These changes are likely to have considerable implications for mountain ecosystem compositing, functioning, and services. The mapping of patterns of alpine phenology and the understanding of the underlying processes on large spatial scales is therefore of importance. Further, according to the space-for-time substitution assumption, the understanding of altitudinal patterns of phenology is useful to estimate future phenological behavior [59].

However, few LSP studies have focused on alpine environments, and in consequence, the processes and changes of mountain phenology are insufficiently investigated. Only in recent years, mountains have begun attracting more attention within the LSP literature, with some studies located in the European Alps. Studer et al. [20] compared in situ observed spring phases with LSP metrics derived from Advanced Very High-Resolution Radiometer (AVHRR) data in Switzerland between 1982 and 2001 and showed their sensitivity to temperature. Also, Fontana et al. [60] focused on the Swiss Alps, analyzing grassland phenology from 2001 to 2005 using AVHRR, SPOT (Satellite Pour l’Observation de la Terre)) VEGETATION, and MODIS (Moderate Resolution Imaging Spectroradiometer) Terra data. The phenological patterns of alpine larch forests and grasslands in the Aosta Valley of northwestern Italy have been studied by Busetto et al. [61] and Colombo et al. [62,63] using MODIS time series for different periods between 2000 to 2009 in relation to climatic factors and elevation. Choler [64] presented a study on grassland phenology in the French Alps, in which plant responses to snow cover duration are analyzed. However, there are only few studies covering the entire Alps. As part of a continental study in which they analyzed AVHRR time series over the period of 1982–2001, Stöckli and Vidale [19] observed for the Alps inter-annual and seasonal variations as well as broad patterns related to topography, but they did not provide information on variabilities within the Alps. Also on an alpine-wide scale, Jolly et al. [65] and Reichstein et al. [66] analyzed the responses of low and high elevation phenology to the 2003 heat wave in the Alps.

The knowledge on LSP patterns and trends over the Alps is hence scattered and limited to specific areas and land cover types. In addition, the available regional studies rely on coarse remote sensing data of 1–8 km spatial resolution. Considering that mountains are heterogeneous landscapes with strongly varying altitudinal gradients and microclimatic conditions, this might reduce the reliability of these analyses [67]. Fisher et al. [10] report that small-scale topographical differences in the order of 50 m can result in a 1–2 week difference in the start of season (SOS), and Inouye and Wielgolaski [68] as well as Kulonen et al. [69] stress the relevance of microhabitat differences. This highlights the necessity of using the highest spatial resolution data available for spatially explicit analyses.

Moreover, no alpine-wide assessment on the relationship between LSP patterns and its climatic drivers exist to date. In this context, an often neglected but very important driver of alpine phenology is snow. Many studies [62,64,68,70,71,72] stress the relevance of snow cover and snowmelt for mountain phenology. The relationship between snow seasonality and phenology is however hardly assessed over large scales and on a per-pixel basis [72]. In a study on the Tibetan Plateau [73], mean Normalized Difference Vegetation Index (NDVI) and snow cover duration (SCD) metrics have been correlated for some individual and accumulated months over ten years, but not looking into the strength of the identified correlations. Xie et al. [74] correlated inter-annual changes of phenological metrics to snow accumulation and melt date, but they conducted this analysis only for the Swiss Alps.

This paper aims at closing these gaps by presenting 17-year time-series of snow cover and phenology for the entire European Alps, using the highest possible spatial resolution of the MODIS land surface reflectance data (250 m). The main objective is to show the potential of the joint analysis of earth observation-based vegetation productivity measures and snow metrics, aimed at answering the following questions:

• Which patterns show the temporal and spatial variabilities of snow and vegetation phenology in dependency of topography and land cover over the Alps?
• Can statistical relationships between vegetation phenology and snow cover be detected in dependency of the altitude? Are there time lags?
• What is the common seasonality between vegetation phenology, snow, and climate parameters? Which parameters are most important in which altitude? Are there time lags?
• Is there a combined influence of climate parameters and snow on phenology?

The above-mentioned data sets are used to answer the research questions (i) and (ii) at an unprecedented temporal and spatial resolution (up to 250 m) and extent (the entire Alps for the years 2000–2017). Since a comprehensive set of climate observations was only available for the Italian province of South Tyrol, the influence of different climate parameters [questions (iii) and (iv)] was tested for this area.

2. Materials and Methods

2.1. Study Area

The study area includes the Alpine range (43.0°–48.6°N/4.2°–17.1°E). The Alps as defined by the Alpine Convention (green shape in Figure 1) cover an area of 191,000 km², stretching across 1200 km and eight states: Austria, France, Germany, Italy, Liechtenstein, Monaco, Slovenia, and Switzerland. Its maximum width is 300 km, between Bavaria and Northern Italy. The Alps are the highest and most extensive mountain range that entirely lies in Europe. The height distribution decreases from west to east with the lowest altitudes being the Mediterranean Sea level and Mont Blanc being the highest peak of the Alps (4810 m).

The climate in the Alps has a very distinct spatial pattern, with temperature tracing altitude and showing a general increasing gradient from north to south and from east to west. Precipitation is highest along the outer chains and more abundant on the northern slopes and at higher altitudes. The rainiest areas are the Jura Mountains, Bernese Alps, Lepontine Alps, Bavarian Alps and the Julian Alps, where precipitation exceeds 2000 mm per year [75]. The zones characterized by low rainfall are those in correspondence of the great longitudinal valleys such as the upper Rhône valley, the Dora Baltea valley, the Valtellina Valley, the Venosta Valley, and the Engadin. Here, annual precipitation is generally below 800 mm per year [75,76]. On smaller scales, the climate is highly related to topography, as altitude and exposure to sunlight and wind influence the individual slopes.

Almost half of the Alpine area (43%) is covered by forests [77]. In the north and the south, deciduous trees dominate the lower slopes, while the upper areas are mostly covered by evergreen forests. Conifers abound in the dry and high altitudes and in the inland areas. Agricultural areas cover almost 40% of the alpine area, of which meadows and mountain pastures make up 18% [78]. Mostly occurring in the highest altitudes, about 10% of the alpine areas are glaciers and perpetual snow, sparsely vegetated or bare areas. Urban areas, which make up about 5% of the area, are the living environment for 14 million people [79] (Figure 2).

As a comprehensive set of climate observations was only available for the Italian province of South Tyrol, an area which is comparable to the entire Alpine range in terms of land cover, altitudinal range, and climate, the influence of different climate parameters was tested here. The Autonomous Province of Bolzano (South Tyrol) is the northernmost province of Italy (red outline in Figure 1). The province has an area of 7400 km². Its altitudinal range varies from 200 m in the southern Adige Valley, to about 4000 m in the Ortler region in the northeast. South Tyrol is characterized by densely populated and intensely cultivated valley floors. Higher altitudes between 1400 m and 2000 m are covered by forests, while pastures, dwarf shrubs, natural grasslands, rocks and glaciers are found above 2000 m. Climate in South Tyrol varies from temperate, humid climate in the valleys, through boreal climate in the forest belt to alpine climate above 2000 m [78]. It is, however, also variable between its different geographical regions. At the valley bottom of the Vinschgau in the west an average of 500 mm of precipitation per year is registered, in contrast to 1300–1500 mm per year at the neighboring higher altitudes [80]. In the Pustertal area in the northeast, precipitation reaches values of 1200 mm per year in the valley and more than 2000 mm per year in the mountains [78,81].

2.2. Data

2.2.1. Elevation and Land Cover Data

The resampled Shuttle Radar Topographic Mission (SRTM) digital elevation model (DEM) v4.1 [82] was used to derive information on altitude and exposition. Using aspect in quantitative analysis is hampered by its circular nature. Instead, aspect, slope and latitude were used to derive the heat load index (HLI), an ordinal index of potential direct incident radiation ranging from 0 to 1 [83], which is used in this study to assess the effect of topography.

For discriminating differently vegetated land surfaces, the CORINE (COoRdination of INformation on the Environment) land cover classification for the reference year 2012 was used to apply one consistent classification for the entire Alpine area (Figure 2). All non-vegetated areas and classes with less than 1% coverage of the overall area were excluded from the analyses. Namely, the ten CORINE classes “211 Non-irrigated arable land”, “231 Pastures”, “242 Complex cultivation patterns”, “243 Land principally occupied by agriculture, with significant areas of natural vegetation”, “311 Broad-leaved forest”, “312 Coniferous forest”, “313 Mixed forest”, “321 Natural grassland”, “324 Transitional woodland/shrub”, and “333 Sparsely vegetated areas” were analyzed.

2.2.2. Snow Data

Daily snow cover maps over the Alps were derived using Terra MODIS surface reflectance data (MOD09GA and MOD09GQ, collection 6), taking into account specific characteristics of mountain areas. The two main improvements of the algorithm compared to the standard MODIS product are the higher ground resolution (250 m) and a tailored topographic correction [84]. The algorithm uses the MODIS bands in the red (0.62–0.67 μm) and near infrared (0.84–0.88 μm) spectrum as well as the NDVI for the recognition of snow, while clouds are classified using bands in the green (0.55–0.57 μm) and shortwave infrared (1.63–1.65 μm) at 500 meters and 1 km resolution. The approach was validated through 148 in situ measurements and shows an accuracy of 82%–94% [85]. The daily snow cover maps were generated for February 2000 to June 2017.

The availability of this daily time series allows the methodology for SCD calculation to be kept as simple as possible [86]. The daily snow cover maps are accumulated and gaps due to clouds or missing data are linearly interpolated. To produce SCD maps that integrate different phases of the winter season and that are comparable with the different aggregated mean NDVI maps (see Section 2.2.4), 16-day SCD, monthly SCD, seasonal SCD, and yearly SCD maps were calculated by applying the algorithm to the data of these different periods. Additionally, the snow cover area (SCA) percentage in each map was computed for the areas of the Alps and aggregated to monthly means.

2.2.3. Vegetation Time Series

An NDVI time series was derived at 250 m resolution for all vegetated areas from the daily Terra MOD09GQ product collection 6 [87]. We chose not to use the BRDF-corrected MCD43A4 MODIS product because of the long composition intervals of 16 days and its reduced spatial resolution of 500 m. The MOD09 product was used in conjunction with the daily MOD09GA product which includes geolocation statistics (solar/sensor angles) and acquisition quality flags at 1 km spatial resolution. 4-day maximum value composites were generated for the years 2000 to 2017. Only those observations in the composites that fulfilled certain observation and quality criteria were used. Namely, pixels with a reflectance value being out of a physically meaningful range, with a cloud flag (either the cloud flagged pixels identified by the Eurac Research snow product (see Section 2.2.2) or, if not available, through the MOD09GA flags “cloudy”, “cirrus”, or “internal cloud algorithm”), or recorded under a local sun or sensor zenith angle larger than 75° were omitted. The selection of the 75° threshold was a compromise with the aim to reduce the level of noise but at the same time not rejecting too many observations to still enable the generation of 4-day composites. Snow covered pixels were masked based on the Eurac Research snow product in order to achieve consistency in the joint data set analysis. In case of a missing snow mask for an acquisition, the MOD09GA quality and internal snow flags were used for masking. Of the remaining, good quality pixel values for each 4-day period, the respective highest NDVI pixel value was chosen for the composite.

The influence of topography was tested by comparing NDVI values calculated based on topographically corrected reflectance values using a C normalization [88] to NDVI based on non-corrected reflectances. The difference was an overall reduced NDVI by 0.34% (which would translate in 0.002 NDVI value difference); thereof, the difference was 0.6% on north-facing slopes and 0.06% on south-facing slopes. The effect of the relief on the NDVI metric was thus assumed negligible. In order to reduce Bidirectional Reflectance Distribution Function (BRDF) effects due to different overpass times, only the Terra acquisitions were used [89]. Further BRDF corrections were not conducted in order not to exclude observation for which the BRDF modeling would fail.

Through the use of a 4-day composite instead of one of the more commonly used 8-, 16- or 30-days composites, we aim at reducing the temporal resolution mismatch between (assumed) changes in phenology of a few days and the usually used (bi)weekly to monthly observations. 8- or more day products might discard acceptable observations and reduce the information content of the input data. The probability of utilizing valid observations is hence higher with a 4-day product. As at the same time the quality criteria for data omission were high, this inevitably led to a higher number of data gaps in the 4-day data set. However, for pixels where there is only one valid observation during 8 days, the information preserved in the 4-day product is not less valid; it is just followed or preceded by a data gap. As the data are used for phenological model fitting (see below) in the next step, individual data gaps in an overall denser time series are no drawback. Also, in the NDVI averaging (see below) potential data gaps do not impair the result, while the inclusion of more valid observation will benefit the accuracy of the mean NDVI.

2.2.4. Phenology Modeling and Monthly Mean NDVIs

The 4-day NDVI composites are used to model the day-of-year (DOY) of the SOS using the TIMESAT software [90]. First, we removed negative NDVI values because they indicate the absence of green vegetation. Then, long data gaps in the winter months were filled as they impair the chosen model fitting (see below) following an approach proposed by Beck et al. [91], that was successfully applied by Zhang et al. [92], Colombo et al. [63] and Busetto et al. [61]. Winter gaps have been defined as data gaps between November 29 and February 1. To fill the gaps, we modeled the winter NDVI, which was assumed to be pixel-specific but constant over years. The gaps are filled on a pixel-wise basis with values of the same date (i.e., 4-day compositing period) from neighboring years or, if in none of the years an observation was available for a date, through a stepwise increase of the search window to the neighboring dates (Figure S1).

Similar to vegetation in high-latitude environments [91], modeling the annual vegetation growth of high-altitude environments faces special challenges such as rather short vegetation periods due to the above described long winter gaps, and steep increases and decreases of the NDVI signal in spring and autumn. Beck et al. [91] showed that a double logistic function is an appropriate method to describe NDVI in such biomes, if the winter NDVI is provided as one of the six model parameters. Hence, in conjunction with applying a median filter (Spike value = 0.5) for outlier removal, a double logistic fitting method (2–3 envelope iterations and adaptation strength 3–8, depending on land cover) was applied in the TIMESAT software. In a last step, a 0.5 amplitude threshold was selected to estimate SOS, in order to cover different characteristics of ground phenology [93]. As one aim of the study is the assessment of relative changes in SOS, we assumed that the choice of threshold will not influence the result as long as the same method is applied to all data sets. The maps have been filtered by removing pixels with a SOS before DOY 30 (30 January) and after DOY 212 (31 July).

In addition to the SOS metric, NDVI time series have the potential to track the temporal development of vegetation activity continuously throughout the vegetation period. Therefore, 16-day and monthly mean NDVI maps were calculated based on the 4-day composites in addition to the SOS, in order to trace the temporal development of vegetation over the year and to correlate it with the temporally highly variable snow cover data (see Section 2.3).

As stated by many authors working in northern latitudes [19,20,94,95,96,97,98,99], detecting vegetation green-up from remote sensing-based vegetation indices (VIs) is difficult in areas affected by snow, as snow acts as a confounding factor in LSP. Changes in vegetative phenophases, i.e., the green-up due to the emergence of leaves, have a similar influence on the reflectance in the red and near-infrared spectral domain as has the melting of snow, i.e., an increase in the NDVI value [100]. Due to the strong linkage of both processes with temperature increase, snowmelt and vegetation green-up can occur very closely in time, rendering the distinction of both processes a challenging task. Some studies have suggested the use of other VIs in order to overcome this issue [11,12,94,97,99,101]. These VIs rely however always on additional spectral information or input data, which is available in the MODIS data only at reduced spatial resolution. In order to account for the above mentioned high spatial heterogeneity of the alpine environment, it was therefore decided to rely on the NDVI time series and to apply the Eurac MODIS snow product for masking snow-influenced observations instead. Through the masking of snow-covered pixels, low NDVI values are removed from the time series. The remaining minimum NDVI of a pixel is hence related to the vegetation minimum (“winter NDVI”) and an increase in NDVI from this minimum is not related to a decrease in snow coverage, but reflects changes of the vegetation canopy. In combination with the above-mentioned gap-filling and thresholding approaches, this ensures that snow melt is not affecting SOS estimation.

2.2.5. Climate Data

Hourly climate data (temperature, precipitation, radiation, relative humidity, wind, air pressure) were available on a 2 km by 2 km grid for the period 2004–2013 in the region of South Tyrol. The Weather Research and Forecasting (WRF) model reanalysis data were provided by the meteorological service company CISMA (www.cisma.it). Besides, day length was derived using latitude and date. All values were aggregated to 16-day means or sums to reduce the influence of noise. Mean temperature was calculated as the average of minimum and maximum temperature. Outlier checks were performed based on which temperature data from June 2005 had to be discarded because of anomalous low values.

2.3. Statistical Analysis

2.3.1. Analysis on Altitudinal and Temporal Variability of Vegetation Metrics and Snow Cover

In order to assess and quantify alpine-wide spatiotemporal patterns and trends in mountain vegetation phenology (SOS and monthly mean NDVI) and snow cover (SCD and SCA), graphs and descriptive statistics (median and standard deviation) were derived from the MODIS data sets. Since different phenological patterns and processes are assumed in different biomes and under different topographic site conditions, individual analyses were done for a range of land cover types (see Section 2.2.1), aspects (i.e., HLI zones) and altitudinal zones in the Alps. The HLI data was split in four classes according to the data sets’ quartiles (which are at 0.62, 0.71, 0.82, and 1.0 HLI in the study region) and assessed in 100 m steps between 100 m and 3000 m of altitude. Changes in phenology and snow cover over the years 2001–2017 were tested through the calculation of linear trends.

2.3.2. Pixel-Wise Correlations between NDVI and SCD

As both data sets, i.e., the vegetation and snow maps show a high regional and inter-annual variability (see Section 3.1), we aim at analyzing the inter-annual influence of snow cover on vegetation development over the entire Alps in a next step. As exploratory analysis we conducted a correlation analysis on the monthly mean NDVI and monthly SCD data sets. The aggregation of both parameters to monthly metrics allowed tracking variabilities throughout the year while at the same time to reduce the influence of noise and data gaps. Pixel-wise Pearson correlations between the 17-year time series of NDVI and SCD were calculated for each month separately (“month-month” correlations) to assess the possible relationship of linearity between them, similarly to the analyses proposed by Grippa et al. [102], Wang et al. [73] and Zhou et al. [103].

Both data sets, i.e., the monthly mean NDVI and the monthly SCD rely partly or entirely on the NDVI information. Therefore, we aimed at identifying a potential influence of using NDVI in the SCD algorithm on the month-month correlations. We conducted an experiment using a simplified version of the snow algorithm that relies solely on an NDVI threshold. Correlations between synthetic data sets of NDVI and SCD that were generated using random samples of NDVI (uniformly distributed between −0.2 and 1) were computed repeatedly (n = 1000). We expected high correlation values even with random values of NDVI if there was an influence of SCD being calculated with NDVI only. The Pearson coefficients of these correlations were however all very low (R² < 0.117). From this we draw the conclusion that using a SCD function of NDVI does not influence the later correlations which integrate also other spectral information.

To identify possible time lags in the relationship, correlations between selected winter month (December–April, “different month combinations”) as well as longer winter periods (“periods–month combinations”) and the remaining months of the vegetation period (March–December) have been calculated similar to Wang et al. [73] and Gessner et al. [104] (Figure 3). While Wang et al. applied their method to a shorter as well as lower temporal and spatial resolution data set than the one used in this study, the main difference to their analysis is the true cross-correlation character in the sense that all reasonable combinations types were tested systematically.

Pixels with a SCD of “0”, pixels for which less than three observations were available within the 17-years vector, and pixels with a standard deviation of zero were excluded from the correlation analyses. The t-test was done to determine whether the correlations are significant at an accuracy level of 90% (p < 0.1). Since the sample size is small (maximum 17 years), and the analysis is only exploratory, we decided to use a p-value of 0.1 instead of the usual 0.05. The type and strength of correlation between NDVI and SCD were assessed for different land cover, altitude and HLI classes.

2.3.3. Common Seasonality and Time-Lags of NDVI, Snow and Climatic Drivers

To evaluate also the importance of other climatic parameters apart from snow for vegetation development in mountain areas, the intra-annual relationship between vegetation, snow and climatic drivers was assessed individually and in combination.

In a first step we tested how closely the seasonality of NDVI is related intra-annually to its individual climatic variables (SCD, mean temperature (tmean), radiation (rad), day length, precipitation) by calculating cross-correlations (Figure 4) between time series in each 2 km grid cell for 2004–2013 over South Tyrol, because for this area and time span, all variables were available (see Section 2.2.5). Since the climate data show a high variability also over short time periods, for all variables (climate, NDVI and SCD metrics) 16-days sums and means were used in the local analyses. Time lags in the relationships up to four data records, i.e., 2 months before and after, were tested. The time lag for the maximum correlation between the variables was identified by taking the 16-day-time-step with the highest absolute correlation. Finally, the results were assessed in dependency of altitude using 100 m altitude classes.

2.3.4. Combined Effects of Snow and Climatic Drivers on NDVI

To account for interactions, i.e., to evaluate the combined intra-annual effects of climate variability on vegetation phenology, NDVI anomalies were associated to snow and climate parameter variabilities in dependency of altitude and HLI, similar to Busetto et al. [61]. This was performed using linear regression models with two-level interactions of the climate variables with topographic variables (altitude and HLI). Interactions were then evaluated at three levels for easier comparison; these were low, medium and high altitudes that correspond to 700, 1500 and 2300 m, and low, medium and high HLI that correspond to 0.55, 0.75 and 0.95; both of which are approximately the 5, 50 and 95% quantiles in the study region.

Climate variables are expected to be highly correlated to each other, thus inducing collinearity issues in the regression modelling. This is mainly due to the seasonal cycle of temperature, snow, and radiation. To remove collinearity, all variables (including the response variable) were first deseasonalized for each grid cell using penalized cyclic splines (mgcv-package in R):

y_ij = f(doy_j) + ε_ij,

(1)

where y_ij is a variable (NDVI, SCD, …) at year i and 16-day group j, f is a cyclic function, and doy_j is the day-of-year. The residuals ε_ij are the deseasonalized values, hereafter denoted with prefix d.

Then for each 16-day group separately, dNDVI was regressed on climate variables, altitude, HLI, and two-way interactions between climate and altitude as well as climate and HLI. All available climate variables were used, except for precipitation, because the cross-correlation analysis before showed no influence. The regression model was run using 50% of the data as training and the other 50% as test set (pixels were randomly selected for each 16-day group). The model formula is as follows:

dNDVI_si = β₀ + β₁ dSCD_si + β₂ dTmean_si + β₃ dRad_si + β₄ WinterSCD_si + β₅ Altitude_s + β₆ HLI_s + Altitude_s (γ₁ dSCD_si + γ₂ dTmean_si + γ₃ dRad_si + γ₄ WinterSCD_si) + HLI_s (δ₁ dSCD_si + δ₂ dTmean_si + δ₃ dRad_si + δ₄ WinterSCD_si) + ε_si,

(2)

where d* denotes deseasonalized variables, s is a pixel index, i is year, β are the main effects, γ and δ are interactions with altitude and HLI, WinterSCD is the SCD of the previous winter (December-February; variable only included in March-November models), and ε_si ~ N(0, σ²) errors [105]. Model selection was not performed, in order to keep models comparable across the 16-day groups. Instead we chose a multiplicity adjustment, and p-values for the coefficients were adjusted for multiple testing (23 models). The total number of estimated parameters for each model was 13 (or 16 if with WinterSCD) and the total number of observations in the training set varied between 7491 and 28,157 depending on 16-day group (less data available in winter because of cloud cover), so overfitting was not considered an issue. This was further confirmed by the small differences between training and test data in the evaluation of model metrics.

3. Results

3.1. Temporal and Spatial Variabilities of Vegetation Phenology and Snow Metrics

3.1.1. Vegetation

SOS maps of alpine vegetation were derived for most years up to an elevation of 2700 m, while for the remaining high alpine areas the available data was mostly too scarce to perform a statistical analysis due to snow and cloud cover. Approaches adapted even further to the characteristics of these high altitudes vegetation signals and noise might be necessary to describe these biomes. Also, the incomplete MODIS time series of the year 2000 proofed to be unreliable to derive phenological metrics.

The derived yearly SOS maps show a high variability with SOS values ranging from 30 to 212 (median = DOY 109.5). Spatial variability is characterized by distinct small-scale patterns that are tracing altitude (see Figure 5). On average over all years, median SOS at different altitudes (100–2700 m) has a time lag of up to 57 days, with median SOS being delayed on average by 2.5 days per 100 m step (adjusted R² = 0.90, p-value < 0.001). This tendency is not valid for the very low altitude ranges from 200–800 m, where SOS is stable around DOY 106 or even a little delayed towards lower ranges (see Figure S2). The year-to-year variability of median SOS in different altitudes follows a bell shape with 16 to 25 days below 800 m and above 1600 m, and its maximum around 30 days from 1100 m to 1400 m.

To test whether slope and aspect related warming and cooling influences vegetation development, the SOS was averaged for the four different HLI and 29 elevation classes. In Figure S2, the SOS curves of the different HLI classes are displayed for each year and the entire altitudinal range, illustrating that the variability between HLI classes in our data sets is small (standard deviation σ = 2.59 days) compared to the intra-annual variability (σ = 7.19 days) as well as to the effects of altitude on SOS (σ = 18.3 days).

In contrast to HLI classes, strong differences between SOS in different land cover classes can be observed. While all classes follow the general trend of later SOS in higher altitudes, land cover classes that are under stronger human influence, such as arable land or complex cultivation patterns, deviate more from a linear relationship with altitude. Furthermore, the intra-annual variability differs between classes, with coniferous forests, pastures and agricultural areas showing a wider spread (Figure 6).

The derived SOS metrics are further more variable among years, with the median of all land cover SOS ranging from DOY 99 (9 April) to 119 (29 April). Years with an overall early SOS are 2014, 2017 and 2001, while 2006, 2008, 2010 and 2012 had a delayed SOS. This variance among years is however not stable over all altitude ranges. While in the years 2017, for example, the median SOS reached a minimum with respect to the previous 16 years, this advance is only distinct below approximately 1500 m, while its SOS is in an average range at higher altitudes (Figure 6). In these elevations, other years such as 2003 and 2015 show the earliest SOS values.

An overall weak negative temporal trend of 4.1 days earlier median SOS per 10 years was observed (adjusted R² = 0.11). Higher altitudinal ranges have a slightly weaker negative trend than the low and medium altitude ranges, while the strongest advance is detected for altitudinal ranges between 1100–1600 m (Figure S3).

Alpine-wide monthly mean NDVI maps were available from February 2000 to June 2017. The median NDVI maps show a strong negative relationship (adjusted R² = 0.77, p-value < 0.001) with altitude, decreasing on average by 0.02 median NDVI per 100 m step. This relationship is however highly variable between different vegetation classes and not valid for the low altitudes below 1000 m, where NDVI slightly increases with altitude (Figure S4).

The temporal variability of the monthly mean NDVI differs for different land cover types (Figure S4) as well as for different periods and elevations (Figure 7). Overall, the NDVI standard deviation is highest from 2000–2200 m (σ = 0.22) and from April to June (σ = 2.2−2.5). Averaged over all altitudes, NDVI shows a weak, non-significant positive trend for 2000–2017 in all months, meaning that NDVI is overall increasing. The only time series deviating from this trend are the month of June and July in 800 m altitude; February, June and July in 1200 m and 1600 m altitude; February and June in 2000 m altitude; and February and March in 2400 m altitude, all in which NDVI is lightly decreasing. While none of the trends are highly significant, the strongest positive trends are those of NDVI in March, April, and August at 400 m, of January at 2000 m, and of October at 2400 m altitude are, though, at accuracy levels below p = 0.2.

3.1.2. Snow

The snow parameters SCD and SCA show a high spatial variability in the Alps, mainly tracing altitude (Figure 8). Averaged over all the years, aspects and land cover types, SCD decreases by 10 day/100 m. This relationship does not vary much over different altitudinal ranges (σ = 0. 37 days). The difference in SCD between HLI classes (σ = 3.7 days) is much smaller than the differences between years (σ = 15.3 days) altitudinal classes (σ = 89.2 days) (Figure S5).

Also, the seasonal variability of SCD between the years is high with alpine-wide median SCD ranging from 37 days in 2011 to 94 days in 2010. The overall trend in median yearly SCD is a non-significant (R² = 0.01, p-value = 0.67) decrease of 4.3 days/10 years. While also none of the trends for the shorter time periods (individual months and different winter periods) show a significant trend (Figure 9), differences among seasons in the general tendency can be detected. While the SCD of October, November, December, and March is slightly decreasing, the SCD of January and February increases (R² values from 0.002–0.07; p-values from 0.87–0.3). Accordingly, also the combined data set of the months January-February is the only period showing an increase in SCD. At the same time, mean SCA of the month December and January varies in the range from 18.1% in 2016 to 58.9% in 2009, with an overall negative trend of 10% decrease/10 years (R² = 0.17, p = 0.13).

The overall decreasing trend in SCD varies additionally among altitudinal ranges. While annual SCD in elevation from 400–1400 m is hardly showing a decreasing trend, the reduction is especially strong in altitudes of 1800 m and above, with the strongest decrease of 15.3 days/10 years in 3000 m elevation and 9.8 days/10 years in 2000 m elevation (Figure 10).

3.2. Inter-Annual Relationship between Monthly SCD and Mean NDVI

The analysis of the inter-annual relationship between the monthly SCD and the average NDVI of the respective month shows strong significant negative correlations (r < −0.5, p < 0.1) for 5–37% of all vegetated pixels of the different month-month analyses for which enough pairwise observations were available, with a maximum in the winter month (December to April) (see for example Figure 11). Positively correlating pixels occur only at 0–2.5% of the entire area, reaching the highest values in September and October. An additional analysis of the inter-annual relationship between the monthly SCD and the average NDVI in dependency of altitude shows that, assessed for each elevation class individually, even up to 55% of all valid observation correlated significantly (a table with all month-month correlation results for all elevations is given in Figure S6). The significantly correlating areas are located especially in altitudes from 1000–2000 m, with a shift of large proportions of correlating pixels to higher elevations in the course of the year (Figure 12). Accordingly, the month with the overall highest percentages of correlating pixels are December to March, while in high altitudes also the spring and early summer month show high proportions, with the month of June correlating best in altitudes above 2400 m (Figure S6).

The correlations between different months achieve overall lower percentages of strongly negative correlating pixels (2–20% for the whole are) than the same months’ correlations, with a relatively high influence of December SCD compared to the other winter months (Figure 13). The amount of strongly negatively correlating pixels is generally decreasing from spring (Mar–May) to summer (Jun–Sep). Thereof, the December SCD correlations follow a different development, with peaks showing also in May and October. Positively correlating pixels are never found on more than 3% of the area.

The winter period–month correlations (green and blue lines in Figure 13) reached similarly shares of negatively correlated pixels (<30%), with maxima in March, April, and May, and in elevations below 2000 m (Figure 14). Generally higher amounts of correlating pixel were found when applying the three-month accumulation period, probably due to the higher variability introduced to the data set by the December data. The percentage of positively correlated pixels does not exceed the 5% threshold, even for the individual altitudinal ranges. Relative maxima of amounts of positively correlating pixels occur only at very low-lying areas below 400 m, i.e., areas strongly influenced by human activities, and in the months March–May in Altitudes above 2400 m.

With regard to the period correlations differentiated among land cover classes, the temporal development is overall the same with 21.1–40.1% of negatively correlating pixels in March, 15.1–28.7% of negatively correlating pixels in April, and further decreasing shares of the alpine area monthly mean NDVI being correlated to the winter SCD, with a minimum of on average only 2.7% of correlating pixels in September. However, differences between the land cover classes exist. Biomes consisting of rather herbaceous and gramineous species, i.e., annual plants, such as ‘pastures’ (on average 11.6% of NDVI-SCD negatively correlated pixels in all month), ‘complex cultivation patterns’ (10.1%) or ‘sparsely vegetated areas’ (9.6%) show a higher relationship with SCD than the forest classes (7.1–7.9%) which have more slowly growing, woody vegetation.

3.3. Intra-Annual Relationships: Common Seasonality of Vegetation Activity and Climate

The seasonality of NDVI in South Tyrol correlated best with the yearly seasonality of different climate variables depending on altitude (Figure 15). In the lowest altitudes up to 300 m, NDVI was correlated highest with radiation (purple) and temperature (orange), then until 700 m almost only to temperature. From 700 m to 2000 m the share of pixels where NDVI was correlated best to SCD (blue) increased continuously reaching 70%, while at the same time the share of pixels correlating best to temperature decreased to 25%. Across all altitudes from 300 m onwards, approximately 5% of pixels had the best correlation of NDVI to day length (green). While precipitation was also compared, its correlations with NDVI were always lower than the best of the other four variables, and so it is not present in the figure.

Of the climate variables, SCD was negatively correlated to NDVI, while the others were positively correlated (Figure S7). The correlation coefficients were high for day length, radiation and temperature (>0.7 up to 1000 m, and >0.5 for higher altitudes) as wells as for SCD (<−0.5 for altitudes higher than 1000 m). Concerning the time lags, SCD had its best correlations at no lag, temperature at no lag for the lowest altitudes and partly at 1 time step (16 days) for higher altitudes. On the other hand, the lag between NDVI and day length or radiation was 2–3 time steps (32–48 days) (Figure S7).

3.4. Combined Effects of Climate Variability on Vegetation Activity Variability

By removing the seasonality from NDVI and climate variables (SCD, mean temperature, radiation, and winter SCD), the combined influence of climate variables on NDVI was assessed. For average altitude and HLI (Figure S8 and black lines in Figure 16), SCD had the highest effect on NDVI, followed by Tmean, radiation and winter SCD, however, the influence of all variables varied throughout the year. SCD had the highest effects in early December (~0.019/d), smaller from January to early April (~0.014/d), decreasing further until August (0.008/d), and then increasing again until November (~0.014/d). Temperature had the highest effects in March and April, peaking early May (0.018/°C), followed by low effects in late May and June (<0.007/°C). In late July and August, temperature had again higher effects (~0.011/°C), followed by minor to non-significant effects throughout September to December. Radiation had the highest effects mid-May (0.0014/W/m²), lasting until mid-July, minor effects August and September, and again higher effects October through December. Winter SCD had negative effects from March until May (0.0049/10d) and minimal positive effects from June until November (~0.0007/10d), except for some dates in early August and mid-September.

The effects of climate variability on NDVI variability depended less on HLI (Figure S10) than on altitude (see colored lines in Figure 16 and Figure S9). The effect of SCD was stronger in lower than in higher altitudes (average difference 0.005/d), but the interaction non-significant in spring (March to May). Higher effects of mean temperature were observed in January and February for lower altitudes, while from mid-April to mid-July effects were higher for higher altitudes. The rest of the year, interactions were non-significant or inconclusive. Radiation effects were higher in lower altitudes throughout the year, except for a few dates, and few non-significant interactions. Effects of winter SCD were more positive at low altitudes for June, October, and early November (~0.003/10d), more positive at higher altitudes for July and late August (~0.002/10d), and negative at higher altitudes late-September until mid-November (~0.002/10d). Interactions with HLI resulted in higher effects of mean temperature for higher HLI for mid-February until late-April (average difference 0.005/°C), and more negative effects of winter SCD for low HLI in March. For the other climate variables and other dates, interactions with HLI were mostly non-significant (Figure S10).

Since only half of the data was used to fit the models, the other half of the data could be used for validation purposes (Figure S11). Climate variability interacted with altitude and HLI could explain on average 30% of NDVI variation from late-October until late-May, and 9% on average for June until early-October. Explained variation was the highest in December and March with values around 40%. The differences between the training and the validation R² were less 1.6%. Median absolute error (MAE) was around 0.08 from mid-November until mid-July, and around 0.06 otherwise. Differences in MAE between the training and validation set were less than 0.0017.

4. Discussion

4.1. Assessment of Spatiotemporal Patterns of Snow and Vegetation Phenology Variability

As expected, a strong dependency of vegetation and snow metrics on altitude is visible from the data, with a delay of SOS of 2.5 days and a shortening of SCD of 10 days per 100 m. The only areas deviating from this trend are situated in low elevations below 800 m, which are densely populated valley bottoms and therefore probably stronger influenced by human activities, as which might be often be influenced by atmospheric inversion, i.e., deviating from the overall temperature patterns. While it has to be considered that the vegetation types occur at varying frequencies in different elevation classes, it could be shown that vegetation in low lying areas and those under stronger human influence and more often consisting of annual plants (e.g., arable land, pastures and complex cultivation patterns) behaves distinctively different than more natural vegetation. The high variability of coniferous forest is an outlier in this respect and needs more investigation, e.g., by testing if this variance might be attributed to difficulties of snow cover mapping in forested areas. This effect is reduced in rather natural classes such as broadleaved forest, natural grasslands, and sparsely vegetated areas. Additionally, elevation-dependent and intra-annual differences between land cover classes can be observed.

Temporal developments of the analyzed metrics tend towards an earlier onset of vegetation growth of about 4 days/10 years, increasing monthly mean NDVI values in spring and late summer, as well as to shorter SCD of 4.3 days/10 years and a reduction in SCA. However, none of these metrics showed a statistically significant trend. This is most probably caused by the high intra-annual variability presented in this study as well as by the insufficient length of the MODIS observations time series to derive robust trends [106]. Hence, the usage of a longer time series might be necessary to derive robust statistical measures, a requirement that hardly can be met at the moment due to the concurrent need of a high spatial resolution.

Although none of the 17-year trends derived from the MODIS data is significant at the 95% accuracy level, they are in the range of changes observed by the IPCC [107] and other studies [8,14,108,109]. Differences in these tendencies occur in different elevations and periods. While the SCD of January and February increased, SCD of October, November, December and March decreased, pointing to a reduction of the snow cover period in the transition months between seasons as well as a concentration of snow fall to a shorter time period. The simultaneous decrease of SCA in the Alps in this period seems contradictory, but might indicate that also the ‘core winter’ is affected by warming. These findings ask however for more investigations on the distribution and amount of snow in the alpine area, an aspect which is out of scope of this study.

The strongest advance is detected for altitudinal ranges between 1100–1600 m, which could be an indication for the enhancement of warming rates with elevation, so-called elevation-dependent warming [50]. The fact that this tendency could not be observed for the highest altitude ranges could be accounted to snow and cloud induced missing data in these ranges (see below). The high intra-annual variability of vegetation metrics in elevation ranges from 1100 m–1500 m compared to other elevation ranges as well as the described change in “early SOS years” around that altitude further indicates an increased sensitivity of vegetation in this altitudinal range to different climatic changes and to intra-annual climate variability.

Analyzing the monthly mean NDVI time series in addition to SOS further revealed developments in vegetation activity throughout the vegetation period and in different altitudes. The high variance of mean NDVI in the months April to June is in accordance to the detected SOS dates, and enables a robust correlation analysis with its climatic drivers (see Section 3.2, Section 3.3 and Section 3.4). While the high variance of NDVI in elevations from 2000–2200 m seems contradictory to the variance of the SOS, it can be probably attributed to the quick and strong changes that alpine vegetation undergoes in these altitudes. The increase of NDVI metrics in low altitudes in March and April as well as in January and October in altitudes above 2000 m indicate spatiotemporal variable shifts of vegetation periods in mountain areas. In this regard the combined use of SOS and NDVI metrics increased the information content extracted from the MODIS time series.

No SOS estimates could be derived for the very high altitudinal ranges due to the high amount of observations affected by clouds and snow, which is a limitation for the monitoring of ecosystems which are potentially very susceptive to climate change. Adapted LSP methods that allow for robust SOS estimation also given more frequent data gaps would be needed for analyzing these areas. Overcoming this constraint would also be beneficial for LSP assessments in boreal and polar regions (see e.g., [94,96,110]). Also, the expected influence of exposition on snow and vegetation phenology is not traceable in the data set. The fact that the investigated metrics did not vary depending on HLI might indicate that the used spatial resolution of 250 m is still too coarse to capture habitat differences that are created by small-scale topographic variations, a shortcoming that was also detected in a study using similar data [111]. As discussed by Comola et al. [112], the effect of aspects become uncorrelated and orientation differences average out at larger scales.

4.2. Derivation of Statistical Relationships between Vegetation Phenology and Snow Cover

In this study, different statistical techniques were used to assess different aspects of alpine vegetation variability such as seasonality and anomalies as well as possible causal relationships. As a logical relationship, for most valid mean NDVI and SCD matches, a generally negative correlation was observed at the landscape scale, i.e., longer snow cover corresponds to lower NDVI. Furthermore, for years with a shorter SCD such as 2016 and 2017, earlier SOS values were observed, especially in lower altitudes, which is comparable with the results of similar, not earth observation-based studies (e.g., [113]).

In general, months and altitudinal ranges with higher inter-annual variability in snow cover show the highest negative relationship of NDVI to SCD. The correlation strength is decreasing from winter (Dec–Mar) to spring (Mar–May) to summer (Jun–Sep). Weak positive correlations between winter SCDs and spring and late summer mean NDVIs in higher altitudes could point to a lagged water storage effect.

The overall highest number of correlating pixels in the same months’ NDVI and SCD maps stress the direct influence snow cover can have on LSP, i.e., the high sensitivity of green-up to snow accumulation and snow melt. However, two confounding effects could interact. While snow has an inhibiting effect on vegetation growth [68,94], snow also exacerbates the seasonal dampening of the vegetation signal, making it difficult to separate the green-up due to the emergence of leaves from an NDVI increase due to the snow disappearance.

The strength of the SCD-NDVI relationship differs further with altitude. Correlating areas are located especially in altitudes from 1000–2000 m, with a general tendency of upward moving correlating pixels throughout the vegetation period. Similarly, differences could be detected among land cover classes. The direct influence of winter SCD on March NDVI was visible on about 30–40% percent of pasture and agricultural areas, but only on 20–25% of forest areas. This illustrates the ability of the data sets to trace the highly variable vegetation responses to snow as a climatic driver in the most sensitive altitudinal ranges. It was also shown, that using SCD data aggregated to two or three month achieved better results than using single month SCD, which is most probably related to the increased necessary variance in the data set that is achieved by a longer integration period.

Given the low number of years investigated, as well as statistical limitations given when investigating an explanatory variable that very often is close to zero, this share of explained variance on the landscape scale is high and in line with similar studies [73,102,103]. However, further remote sensing based phenological metrics such as end of season, length of season and the first and last snow day could be applied to asses in more detail also the phenological mechanisms in the senescence phase, a process that is overall less good understood and modeled [18,114,115].

4.3. Assessment of the Common Seasonality and Combined Influence of Climate Parameters and Snow on Phenology

To understand the relative importance of snow cover variations as phenological driver with respect to the other variables, a more detailed statistical analysis of vegetation metrics, snow cover and a range of climatic variables was performed on the limited area of South Tyrol. Therefore it was shown, that mean temperature alone explains vegetation seasonality across all altitudes, but better in lower areas. SCD alone, on the other hand outperforms in higher altitudes. The transition altitude between these two dominating effects is also the elevation which showed the highest year-to-year variability in vegetation metrics (see Section 4.1). It can therefore be concluded, that in mountainous areas, snow and temperature can be combined to produce equally good results in phenology modeling across all altitudes.

In the combined effects, SCD had the highest effect on NDVI. Overall, climate variability interacted with altitude and HLI could explain on average 30% of NDVI variation from late-October until late-May, and 9% on average for June until early-October. Explained variation was the highest in December and March with values around 40%. The combined effects of winter SCD were negative from March until May, more positive at low altitudes for June, October, and early November, more positive at higher altitudes for July and late August, and negative at higher altitudes late-September until mid-November. This shift of the peak of the effect might be an indication of a water storage effect; however, the effect was several magnitudes lower compared to the other variables. As for the variability of snow and vegetation metrics, the effects of climate variability on NDVI variability depended less on HLI than on altitude. As discussed above, this is likely a scale-dependent result, which is even reinforced with the resampling of all data to the 2 km spatial resolution [112]. Nevertheless, it seems that small variations in highly dynamic periods, such as the higher effects of mean temperature for higher HLI for mid-February until late-April and more negative effects of winter SCD for low HLI in March, might be reinforced by exposition.

Overall the high variance explained by SCD stresses the importance of snow as a phenological driver that might alter plant development even more strongly in the future. These findings enable also a better interpretation of the outcomes of the analysis on the whole European Alps executed considering snow cover variations as the only driver of phenology (Section 3.2), imparting them a high validity even though only snow as climatic driver was analyzed.

5. Conclusions

In this study, we jointly analyzed the temporal and spatial variability of snow and plant phenology on an alpine-wide scale at 250 m resolution. The presented regional analysis allows to assess the spatiotemporal patterns of earth-observation based snow and vegetation metrics over different land covers, altitudes and topographic features, as well as to understand the relative importance of snow cover variations as phenology driver with respect to the other variables and in different elevations. Both snow and vegetation parameters showed clear patterns related to topography, with a delay of SOS of about 2.5 days and an increment of 10 days in SCD per 100 m. No significant temporal trends of earlier SOS as well as of shorter SCD could be detected on the alpine-wide scale, probably due to the high regional variability and insufficient length of the MODIS data set. The statistical interrelationships of plant photosynthetic activity and its drivers revealed the high explanatory power of SCD, also in combination with other climatic drivers and especially in altitudes above 1000 m. However, the strength of the correlations and the combined climate influence varies throughout the year and for different altitudes. While the time series used as a data base in this study is not long enough to derive statistical robust climate-relevant trends, it nevertheless stresses the relevance of monitoring mountain areas, which are among the regions of the globe that are warming the most, but for which at the same time, monitoring the rate and patterns of warming is challenging [49].