Remote sensing-based quantification of spatial variation in canopy phenology of four dominant tree species in Europe

Abstract Spatial variation of phenology is a central feature of global change research. Satellite remote sensing is used for continental to global monitoring due to the limitations of long-term field observations of plant phenology. A threshold method was used to estimate the start of the season, length of the season, maximum normalized difference vegetation index (NDVI), and integral NDVI for selected tree species using remote sensing based NDVI data acquired by the VEGETATION instrument on board Satellite Pour l’Observation de la Terre (SPOT VGT NDVI). Afterward, the spatial patterns in the satellite-derived phenological metrics for four dominant tree species (i.e., beech, birch, pine, and spruce) across Europe were characterized. The results indicate that: (1) The SOS occurs 1.6–2.9 days later and the average LOS is 2.7–3 days shorter per 1 deg of latitude increase from south to north. (3) The SOS occurs 0.7–1.8 days later and the LOS was 0.6–2 days shorter per 100-m increase in altitude for the four species. (4) The SOS and LOS across Europe are well correlated with the mean annual air temperature (1°C correlates with a 4.5-day advance in the SOS and a 7-day extension in the LOS). Our research is the first one to characterize the spatial and temporal variations of phenology for different tree species across Europe using remote sensing.


Introduction
Phenology includes the study of annual rhythms of biological phenomena mainly in relation to climatic parameters. 1 When applied to plants, phenology includes such critical events as the beginning of the season, end of the season, or start of such major growth phases such as budding and flowering. The research on plant phenology has improved the understanding of the variations in lifecycle events in the past centuries 2-6 using site-and species-specific phenological measurements. 7,8 However, phenological data from many parts of the world are comparatively scarce. Moreover, ecological models relating to climate-change studies require phenological information at large spatial scales rather than inventory data that focus mostly on specific plant species and are mostly point observations. As a result, the use of remote sensing data for inferring the phenological characteristics of vegetation is increasingly regarded as a key to understanding large-area seasonal phenomena. 9 Among all remote-sensing indices, the normalized difference vegetation index (NDVI) is the index most often used in studying seasonal variations, in growing vegetation, and the shifts between vegetation covers. [10][11][12] The study of vegetation phenology using remote sensing has been called land surface phenology (LSP) by de Beurs and Henebry. 13 Over the past two decades, LSP has served as a useful tool for many purposes, including monitoring biospheric activity, 14 developing prognostic phenological models, 15 and deriving land-cover maps. 16 Myneni et al. 17 presented evidence from satellite data to show that the photosynthetic activity of terrestrial vegetation yielded an increase in plant growth associated with a lengthening of the active growing season. Tucker et al. 18 used NDVI data from NOAA satellites from 1982 to 1999 to estimate variations in photosynthetic activity and growing season length at latitudes above 35°N and achieved the same results as did Myneni et al. 19 Wang et al. 20 revealed patterns and trends of landscape dynamics, LSP, and ecosystem production along the Appalachian Mountains using the time series data from global inventory modeling and mapping studies and advanced very high resolution radiometer global production efficiency model datasets. Dunn and de Beurs 21 used the MODIS satellite 8-days 500 m Nadir BRDF Adjusted Reflectance product to quantify the LSP. They presented that the elevation and latitude exhibit the most significant influences on the timing of start of the season throughout the North American mountainous area.
Phenology as a major indicator of forest ecosystem status 22 is pivotal in forest practices and strategies. Forests cover 44% of Europe's land area and they continue to expand. However, as explained by Alley et al., 23 forest ecosystems seem to be especially prone to climate change, and European forests are among the most intensively used forests in the world. Against this background, analyzing the status of forest ecosystems and their dynamics builds a substantial basis for future conservation practices and strategies. Furthermore, with the Kyoto protocol, interest has been focused on trees for their potential role in mitigation of greenhouse gasses and carbon sequestration. 24 Therefore, studying phenological variations among various tree species at a continental scale is of great importance.
The primary objective of this study was to use remotely sensed phenological data derived from SPOT-VGT time-series NDVI profiles to characterize the spatial pattern of LSP for four dominant tree species across Europe. Specifically, our strategies and objectives were as follows:

Data source and data preprocessing
In this study, the SPOT-VGT S10 NDVI data (1 km × 1-km resolution) acquired from 1999 to 2005 across the study area were used. The high-temporal resolution of 10 days allowed for image selection according to the highest quality, least cloud cover, and the optimal phenological stage of vegetation cover. 26 Atmospheric corrections were routinely done using the SMAC model 27 for evaporation, ozone, and aerosol effects. Cloud masking was applied using the 8-bit status map provided by the SPOT-VGT website to eliminate the influence of clouds and atmospheric variability.
The presence of snow and the polar night during winter complicated the use of NDVI time series for modeling vegetation activity at high latitudes. As snow has a negative effect on the NDVI, the disappearance of snow at the end of winter causes the NDVI to rise. In addition, the NDVI values are missing due to the polar night in northern Europe. Hence, the NDVI values outside of the growing season needed to be calculated to cover the shortage of data due to these effects. 28 We first used the Pan-European Map of Biogeographical Regions from the European Environment Agency to mask the NDVI data. 29 Depending on the climatic characteristics, the ecoregions are separated into three categories: boreal, continental, and Mediterranean. Then, the winter NDVI for the boreal region was calculated using the method proposed by Beck et al. 28 In addition, an up-to-date tree species inventory map of Europe with a resolution of 1 km and a set of global climate layers with a spatial resolution of about 0.25°C were used in this research. 30 The overall accuracy of the classification amounts to 86.6% across Europe. 25 We selected four tree species widely used in the timber industry and with a widespread natural distribution across Europe. To reduce the influence of mixed pixels for a single tree species due to the 1-km resolution, we used 100-m resolution coordination of information on the environment (CORINE) land-cover data to refresh the tree species map and chose only the pure pixels (with a 100% of correspondence to tree species), as shown in the map in Fig. 1. There are around 30,000 pure pixels selected. The CORINE data are the only Pan-European land-cover maps based on remote sensing data with a spatial resolution finer than 100 m. The CORINE land cover 2000, with the abbreviated name CLC2000, was derived using Landsat TM and ETM+ imagery and visual interpretation. The CORINE land-cover maps 2000 were adopted in this research for refreshing the tree species map. The majority of the training data were derived via overlaying of co-registered coarse-resolution novel tree species maps and the interpreted high-resolution datasets from CORINE.
A set of global climate layers (climate grids) with a spatial resolution of about 1 km were acquired from a gridded observational dataset for precipitation and temperature measured daily in Europe (E-OBS) based on the information from the European Climate Assessment and Dataset spanning the period 1999-2005. 31 The temperature resolution of the E-OBS data is approximately 0.25°C. The data files contain gridded data for five elements (daily mean temperature, TG; daily minimum temperature, TN; daily maximum temperature, TX; daily precipitation total, RR; and daily averaged sea level pressure, PP).

Calculation of phenological metrics
A widely used time-series analysis program, TIMESAT, was used to calculate the LSP metrics based on the SPOT-VGT data. [32][33][34] TIMESAT uses three processing methods, i.e., adaptive Savitzky-Golay, Gaussian, and double logistics, for time-series smoothing. However, as explained by Beck et al., many of the algorithms used for vegetation monitoring are ill-suited to boreal regions, 28 whereas double logistics are the most suitable for describing vegetation dynamics at high latitudes. During the process of deriving phenological metrics, we eliminated NDVI spikes larger than two times the standard deviation of the median values of the closest neighbors in the time series and fitted the remaining upper envelope. 35 SOS is defined as the interpolated composite period when the NDVI has increased by 20% of the seasonal amplitude above the growing season minimum level. [34][35][36] The seasonality metrics are shown in To assess the spatial variability in phenological metrics among different tree species, we computed the date of SOS, the annual MaxNDVI, the LOS, and the integral NDVI for different tree species. For preliminary analyses, the region was divided into areas of 1 degree of latitude by 1 degree of longitude. Elevation zones with the intervals of 100 m were assigned using the UNEP-WCMC 37 data. For each species and in each spatial region, naturally established populations were sampled using the 100% species-pixel criterion explained above, as indicated by the tree species map. The relationships between the annual timing of a phenological phase and altitude, latitude, longitude, and climate factors were then determined using a linear regression. 38 3 Results

Phenological Metrics
To examine the spatial variability of the phenological metrics, four different NDVI metrics were mapped, as shown in Fig. 3. These metrics were averaged over 7 years to show the general gradient along the geographical transitions. As seen in Fig. 3, phenological parameters vary from east to west and north to south. Pixels that correspond to water surfaces, rock, deserts, or other surface features with insufficient NDVI dynamics, where phenological metrics could not be calculated, are shown in white. Figure 3(a) shows the dates for SOS as estimated from SPOT-VGT data between 35°N and 70°N and averaged over the years 1999-2005. The phenological patterns depend strongly on latitude in the mid-to high-latitude regions. The dates of the SOS advance northward from March to early June and eastward from March to early May. The growing season starts in most parts of central Europe between the 20th of March and the 15th of April. Note that a late date of the SOS after June was observed in most of the coastal areas next to the Mediterranean Sea because heat stress and seasonal precipitation patterns differ in these areas versus those in neighboring areas. In mountainous regions (Fig. 4), such as the Alps Mountains, great differences in the phenological metrics can be seen: the growing season starts much later than in nearby areas of lower altitude. The latest dates of the SOS, i.e., after the 1st of June, are observed in the Alps and correspond to similar dates in Norway and north of the Arctic Circle.
Almost the same gradient was found with regard to the LOS. The longest growing seasons, i.e., over 220 days, can be found in the coastal regions of southern Europe. In large parts of central Europe, except for the mountainous regions, the growing season lasts between 180 and 220 days. Shorter growing seasons, i.e., <170 days, are seen in zones of high altitude and in nearly all of Scandinavia (55°N-70°N). High-altitude zones in Scandinavia and the areas north of the Arctic Circle show growing seasons as short as 120 days (Fig. 4).
The mean MaxNDVI is lower in the northern hemisphere (boreal area) and the Mediterranean area than in the central continental areas, while the value of MaxNDVI over the central continental area is relatively consistent [Fig. 3(c)]. The pattern of mean integral NDVI is the combined function of MaxNDVI and LOS, 39 and the largest value of integral NDVI was found in Eastern Europe. The integral NDVI increases when moving from the Mediterranean area to the central continental area and decreases again when moving farther north from the central continent to the boreal areas. The integral NDVI also generally increases moving from Western to Eastern Europe [ Fig. 3(d)].

Phenological Characteristics of Different Tree Species
To evaluate whether the differences in phenology are caused by the tree species or by geographical factors, samples of the four main tree species were selected at representative locations for each species and tested by one-way analysis of variance in SPSS. Since beech and birch trees tend to show little overlap in their locations, two groups of analyses were carried out: one for beech, pine, and spruce and the other for birch, pine, and spruce. The results indicate that the variability introduced by species is not significant in the case of beech/pine/spruce, i.e., they exhibit the same phenological characteristics within a given region. However, in the case of birch, the dates of the SOS and LOS are significantly different from those of pine and spruce even at the same latitude (Table 1).

Spatial Variations for Four Tree Species Along the Geographical Gradient
The results of the regression analysis for different phenological phases are listed in Table 2. The Mediterranean areas are excluded, considering the fact that the vegetation cycle begins at end of the year. 40 The dependence of the SOS, LOS, maxNDVI, and integral NDVI on altitude, longitude, and latitude, averaged over the years 1999-2005, can be derived from the regression coefficients.
As we can see from the results, the dependence of phenological metrics on latitude is apparent except in the case of beech trees, whose metrics are more dependent on altitude. The start  The relationship between SOS and longitude is significant only in the case of birch: SOS occurred later moving from west to east by 0.7 day per 1 deg of longitude. The LOS was 0.7 day shorter, the maximum NDVI showed a steady gradient, and the integral NDVI decreased by 0.8 with each 1 deg increase in latitude.
In addition, the dependence of phenological metrics on altitude was also significant, especially for beech. The SOS occurs 1 day later, except for beech, 0.7 day later for birch, and 1.8 days for spruce with each 100-m increase in altitude. The LOS was 1.3 days shorter for beech, 0.6 day for birch, and 2 days for spruce with each 100-m increase in altitude. Any gradients in the maximum NDVI and integral NDVI were not significant, i.e., no single later or earlier phenological gradient was observed.

Relationship Between Satellite-Derived Phenology and Temperature
To assess the relationship between annual air temperature (Tmean) and phenological metrics, the linear correlations between them [Figs. 5(a)-5(d)] for different tree species were accessed. The Mediterranean area is excluded considering the different vegetation cycles. 40 The correlations suggest that the temperature is a more significant controlling factor in NDVI response than in precipitation. Figure 5 shows that all species show a negative correlation between the SOS and Tmean and a positive correlation between the LOS and Tmean, meaning that higher temperatures promote an earlier SOS. The regression between the SOS and Tmean indicates that a 1°C increase in mean air temperature is associated with an earlier start of the season by about 4.5 days. Similarly, a 1°C increase in mean air temperature tends to extend the length of the season by 7 days (Fig. 5). Compared to the SOS and LOS, the maxNDVI and integral NDVI show less correlation with the temperature (Fig. 5).
In order to investigate the relationships between the phenological metrics and temperature, regression analysis using phenological metrics with several temperature factors are carried out: temperatures at the month of SOS (T0), mean temperature of the following two months of SOS (T1, T2), mean temperature of the month of SOS and the following two months (T012), mean temperature of the preceding two months of SOS (T-1, T-2), and mean temperature of the month of SOS and the two preceding months (T0-1-2) ( Table 3). Table 3 shows that the phenological metrics are mainly influenced by the air temperature during the month of the start of the season. Also, during two months preceding the SOS, significant correlation coefficients between temperature and phenological metrics were found.

Discussion
In this paper, remote sensing data were used to monitor the phenological variations among and between four tree species at a continental scale in Europe. With the assumption that NDVI profiles of natural vegetation can be characterized by a main growing cycle, the evolution of vegetation phenology (SOS, LOS, MaxNDVI, and integral NDVI) was assessed using SPOT-VGT NDVI data. To obtain continental-scale NDVI data, eliminating the influence of spurious factors like snow or polar light was crucial, and the winter NDVI was first calculated for boreal regions. Then, a double logistics algorithm was chosen to estimate the phenological metrics for the whole area. Finally, the spatial variations seen in four dominant tree species (beech, birch, pine, and spruce) were evaluated, and climatic data were introduced as possible explanatory factors for those variations.

Phenologic Metrics Derived from Time Series NDVI Data
The maps of the phenological metrics (Fig. 3) are derived using 7 years of averaged data rather than just one year of data. This decision may have excluded the influence of external bias, such as extremely abnormal years with high temperatures or an unexpected NDVI error when using data obtained from remote sensing. Note that the main period of growth of vegetation in Mediterranean regions is in the winter and spring because of heat stress at other times and seasonal precipitation patterns. These interpretations coincide with our findings presented in several published studies. 41,42 The average length of the growing season in Europe [ Fig. 3(b)], matches well with the results modeled by Chmielewski and Rotzer. 43 At the same time, we also developed the same finding reported by Köble and Seufert 44 of a gradual west-toeast spread of spring phenophases across Europe in warmer years. However, among the four species of trees we studied, only birch displayed a clear west-to-east gradient ( Fig. 3 and Table 2).

Phenological Characteristics of Different Tree Species
The timing of the phenological events of different species on a large (continental) scale is an interesting phenomenon. Which factors dominate the phenological metrics may still be poorly understood. This study revealed large differences in phenological metrics between different tree species, in common with other studies. 45,46 Deciduous forests types can clearly be distinguished from evergreen forests, which are mostly present in the boreal zone and alpine areas. The results from Table 1 suggest that site-specific factors, tree species, and phenology display a complex set of interrelationships. Tree species may be locally adapted to their particular environments, but variations in species still account for part of the phenological variations, especially on a large scale, i.e., birch/pine and spruce exhibit widely differing phenological patterns within the same selected area. The dependences of phenological metrics on spatial factors varies with the tree species. Rötzer and Chmielewski 47 explained that the growing season starts after the 25th of April in Scandinavia and the Baltic Sea region, which was confirmed with our findings. The SOS observed between latitudes 39°N and 49°N is generally 12 days later compared to that observed in the latitude zone 49°N-55°N, due to the presence of the Alps which can lead to a significant differences between the plant phenology observed on the south side of this mountain range versus that observed nearby to the north. 48 We were unable to develop consistent estimates of the SOS in the Mediterranean area, an ecoregion with pronounced, regular wet and dry seasons (Fig. 3). The switch from dry to wet occurs around the end of the December to early January, which means that the SOS occurs much later here than in other areas. 49 The dependence of phenological metrics on altitude was significant, especially in the case of beech trees. The findings of König and Mayer, specifically a beginning of spring phenophases that was 1 to 3 days later with each 100-m increase in altitude across various parts of Europe, corresponds well with our results. 50

Relationship Between Phenology and Climate Factors
In our results, phenological metrics are correlated with the temperature variables. Thus, we can deduce that analyses of the phenological data reveal evidence of climate change, although the plant responses may sometimes be inhomogeneous due to the local microclimate conditions, natural variations, genetic differences, or other nonclimatic factors. 51 Similar results were obtained by Chmielewski and Rotzer. 43 Those authors concluded that the length of the season depends highly on the mean annual air temperature. Several authors have established some relationships between phenology and climatic conditions at continental or regional scales. [52][53][54] Our results suggested that the phenology corresponds well with changes in air temperature associated with the onset of spring, especially with the mean temperature during the month of the start of the season. This finding was also described by several other researchers. Chmielewski and Rötzer indicated that most phenological phases correlate significantly with the mean monthly temperatures of the month of onset and the two preceding months. 52 Tryjanowski et al. also described the manner in which the SOS mainly depends on the temperatures from February to May. 45

Conclusions
Remote sensing is a useful tool for generating consistent estimates of the phenological metrics over large areas. In general, the phenology estimated using remote sensing is coarser than that derived from direct observations of individual plant phenology, such as bud burst or first leaf; yet, this method can yield useful pixelated representations of the constituents of vegetation variations. 55 Our results highlight both the challenge and potential for using remote sensing as a tool in phenological analysis. No other technology aside from remote sensing offers consistent longterm monitoring, yet its validation poses obstacles due to limitations in using existing historical datasets. Establishing consistent plant phenological monitoring networks is therefore essential. We ultimately conclude that remote sensing is a valuable tool for estimating phenological metrics over large areas. Thus, our research may represent the first characterization of the spatial and temporal variations in the phenology of various tree species across Europe using remote sensing. An understanding of how phenological variations occur for different tree species is important to the tree breeder, entomologist, physiologist, ecologists, and so on. Moreover, ful comprehensions of tree species phenology provide guidance on carbon and species distribution modeling in the future.