Detecting InterAnnual Variations in the Phenology of Evergreen Conifers Using Long-Term MODIS Vegetation Index Time Series

Long-term observations of vegetation phenology can be used to monitor the response of terrestrial ecosystems to climate change. Satellite remote sensing provides the most efficient means to observe phenological events through time series analysis of vegetation indices such as the Normalized Difference Vegetation Index (NDVI). This study investigates the potential of a Photochemical Reflectance Index (PRI), which has been linked to vegetation light use efficiency, to improve the accuracy of MODIS-based estimates of phenology in an evergreen conifer forest. Timings of the start and end of the growing season (SGS and EGS) were derived from a 13-year-long time series of PRI and NDVI based on a MAIAC (multi-angle implementation of atmospheric correction) processed MODIS dataset and standard MODIS NDVI product data. The derived dates were validated with phenology estimates from ground-based flux tower measurements of ecosystem productivity. Significant correlations were found between the MAIAC time series and ground-estimated SGS (R2 = 0.36–0.8), which is remarkable since previous studies have found it difficult to observe inter-annual phenological variations in evergreen vegetation from satellite data. The considerably noisier NDVI product could not accurately predict SGS, and EGS could not be derived successfully from any of the time series. While the strongest relationship overall was found between SGS derived from the ground data and PRI, MAIAC NDVI exhibited high correlations with SGS more consistently (R2 > 0.6 in all cases). The results suggest that PRI can serve as an effective indicator of spring seasonal transitions, however, additional work is necessary to confirm the relationships observed and to further explore the usefulness of MODIS PRI for detecting phenology.


Introduction
Developing robust methods for monitoring terrestrial photosynthetic activity is key to improving our understanding of vegetation dynamics and the carbon cycle in the context of climate change.Temporal shifts in vegetation phenological events, such as green-up and senescence, are effective indicators of climate change [1] and can have a significant impact on the annual carbon uptake by forest ecosystems [2,3].Phenological changes can be assessed from various data sources, including field observations and gas exchange measurements from flux tower sites.At regional to global scales, satellite remote sensing provides the most efficient means to observe spatial and temporal variations in vegetation productivity and phenology.Phenological parameters can be extracted from per-pixel time series of spectral vegetation indices (VIs) that have been closely linked with vegetation productivity [4,5].Long-term studies using VI time series have revealed widespread greening occurring in various regions and ecosystems [6][7][8].Ongoing efforts are being made to establish the accuracy with which satellite imagery can be used to predict timings of seasonal transitions by validating the estimates with ground-based observations [7,[9][10][11].
Remote sensing time series studies of phenology are commonly based on the Normalized Difference Vegetation Index (NDVI) [12].The NDVI relies on the contrast between absorption of red light by leaf pigments and strong scattering of near-infrared radiation by foliage.As a consequence, changes in NDVI are linked to changes in various vegetation biophysical properties including leaf area, green biomass and the fraction of absorbed photosynthetically active radiation (fAPAR) [13].While studies utilizing NDVI and similar broad band vegetation indices have provided valuable insights into the response of vegetation to climate change, such indices are only partially linked to photosynthetic activity.Seasonal and periodic down-regulation of photosynthesis can occur while vegetation continues to absorb radiation with very little change in leaf area and greenness and therefore causes no corresponding changes in NDVI.Thus, the NDVI fails to capture the efficiency with which vegetation converts the absorbed light into biomass and therefore is unable to detect all photosynthetic variation [14].Vegetation light use efficiency (LUE) is influenced by many environmental and physiological factors, so it varies temporally and across different biomes, species and plant functional types [15,16].Taking the LUE component into account may have the potential to improve the accuracy of satellite-based estimates of vegetation productivity and phenology.This is of particular importance in coniferous forests and other evergreen-dominated ecosystems as these exhibit less temporal variability in NDVI compared to deciduous vegetation [14].
A promising approach for assessing variations in LUE involves the use of narrow spectral bands to detect changes in plant pigments.Pioneering studies by Gamon et al. [17,18] linked changes in xanthophyll pigment cycle associated with reduced LUE and plant stress with changes in reflectance at 531 nm, which enabled the development of a Photochemical Reflectance Index (PRI).The PRI combines a narrow band centered at 531 nm and a normalizing reference band unaffected by LUE.Xanthophyll cycle pigments provide a means of dissipating excess energy when light absorbed by plants exceeds the photosynthetic capacity to avoid damage to the photosynthetic system, which influences the PRI signal on a diurnal timescale.It was later demonstrated by Porcar-Castell et al. [19] that on seasonal scales the main physiological mechanism controlling the relationship between PRI and LUE is not the xanthophyll cycle, but rather dynamics in the relative pigments pools of chlorophyll and carotenoids.This was further confirmed by Gamon et al. [20] and Wong and Gamon [14], who observed high correlations between PRI responses and adjustments in chlorophyll:carotenoid ratios.
The PRI was developed from leaf-scale measurements, however, an increasing number of studies have demonstrated the sensitivity of PRI to photosynthetic activity at the canopy and ecosystem scale across different vegetation types [21][22][23].Nichol et al. [24] found that the PRI served as a useful indicator of LUE specifically during spring transitions in boreal forest canopies.The Moderate Resolution Imaging Spectrometer (MODIS) satellite sensor has commonly been applied in PRI studies [25][26][27][28][29][30][31][32].Although MODIS is not specifically designed for the purpose, it has several narrow ocean bands (10 nm), including a band centered at 531 nm, and a high temporal resolution (1-2 days revisit time globally), which makes the sensor suitable for time series studies of photosynthetic activity and LUE.
MODIS-based studies use different versions of the PRI (i.e., different reference bands) since the sensor does not include a band centered at 570 nm, identified as the most useful reference wavelength at the leaf scale by Gamon et al. [17].The optimal wavelength bands vary from site to site [27,30,33], which highlights some well-known limitations of using PRI as well as other VIs.When applied at larger scales, the VI signal is affected by various confounding effects, which are not related to the canopy itself [34].These are the result of bidirectional reflectance distribution function effects caused by varying viewing and illumination geometry [22], canopy structure and background effects (e.g., soil color and shadows) [30] and changing atmospheric conditions [34].These include compositing methods such as the CV-MVC algorithm used to generate global MODIS NDVI products [35] and temporal smoothing of the data points by applying filters or fitting functions [36][37][38].A new study by Middleton et al. [39] shows that these confounding factors can be largely ameliorated by use of standardized geometry and processing approaches.Snow and cloud cover is particularly problematic when observing phenology in boreal regions as the noise arising as a result tends to decrease the VI signal significantly and consequently reduce the temporal resolution of the VI time series.A wide variety of strategies have been developed to minimize the impacts of noise and improve the stability of the time series Recent work by Wong and Gamon [14] showed that changes in PRI and pigment pool sizes were closely timed with spring photosynthetic activation in conifers measured by gas exchange on the ground.The study was carried out at the leaf and stand scales and the authors encourage further work using MODIS PRI to demonstrate these relationships at larger scales and at multiple sites.The application of MODIS in PRI studies has been hindered by the fact that the narrow bands used are designed for ocean color detection and are therefore not included in standard processing over land regions.Rahman et al. [32] and Garbulsky et al. [33] highlight this lack of processed and quality controlled MODIS ocean band data over terrestrial regions as a major hindrance to long-term studies of PRI.However, recently NASA has processed a large dataset, which consists of 16 years of gridded daily surface reflectances in band 1-12 (including the 531 nm band) for North America, South America and Europe.The dataset is atmospherically corrected using a multi-angle implementation of atmospheric correction (MAIAC), which uses time-series and image-based processing to derive aerosol optical thickness and surface reflectance [40,41].In comparison with a standard atmospheric correction (6S), Hilker et al. [42] found that MAIAC-processed data markedly improved the MODIS signal in a PRI-based study of vegetation LUE.The MAIAC dataset thus presents new opportunities for testing the applicability of PRI as an optical indicator of phenological events and exploring long-term trends in vegetation productivity in evergreen ecosystems.
The overall aim of this study was to investigate the utility of MODIS-derived spectral indices for inferring photosynthetic activity and detecting seasonal transitions in a southern boreal coniferous forest.The specific objectives were: 1.
To derive phenological parameters from a 13-year NDVI and PRI time series based on MODIS MAIAC-corrected data and validate the observations with estimates from in-situ ecosystem gas exchange measurements.2.
To compare the performance of MAIAC data in predicting phenology against the standard MODIS NDVI product.

3.
To assess the relative utility of NDVI versus PRI as optical indicators of seasonal variations in vegetation photosynthetic activity.

Study Site
The study site is located at the Station for Measuring Forest Ecosystem-Atmosphere Relation (SMEAR II) in Hyytiälä, southern Finland (61 • 51 40"N, 24 • 17 13"E) (Figure 1).The surrounding terrain is dominated by Scots pine (Pinus sylvestris) planted in 1962 [43].The forest area is homogenous with less than 10% of the canopy made up of other species within 200 m of the flux tower [44].
The flux tower at the site is equipped with instrumentation to measure several gas concentrations, atmosphere-biosphere exchange of energy and matter and meteorological parameters since 1996.Detailed panoramic images of the forest can be accessed via the University of Helsinki SMEAR web portal (https://www.atm.helsinki.fi/SMEAR/index.php/smear-ii).
Remote Sens. 2017, 9, 49 4 of 21 homogenous with less than 10% of the canopy made up of other species within 200 m of the flux tower [44].The flux tower at the site is equipped with instrumentation to measure several gas concentrations, atmosphere-biosphere exchange of energy and matter and meteorological parameters since 1996.Detailed panoramic images of the forest can be accessed via the University of Helsinki SMEAR web portal (https://www.atm.helsinki.fi/SMEAR/index.php/smear-ii).

Flux Tower Data
The net ecosystem exchange (NEE) was measured using the eddy covariance (EC) technique [45].Shortly, the measurement system, located 23.3 m above the ground, included an ultrasonic anemometer (Solent Research 1012R2, Gill Instruments Ltd., Lymington, Hampshire, UK) for measuring three wind speed components and temperature and a closed-path infrared gas analyzer (LI-6262, LI-COR Biosciences, Lincoln, NE, USA), which measured CO2 and H2O concentrations.The vertical CO2 fluxes were calculated over 30 min averaging periods by block averaging approach using the EddyUH software [46].The fluxes were corrected for high frequency response underestimation by using the co-spectral transfer functions calculated according to Mammarella et al. [47] together with in-situ parameterization of the co-spectral model derived from ensemble mean of measured temperature co-spectra [47].The CO2 flux was then corrected for storage change to obtain NEE.The 30 min averaged NEE was partitioned into total ecosystem respiration (TER) and gross primary productivity (GPP) as described in Kolari et al. [48].The air temperature and the incoming photosynthetically active radiation (PAR) were also measured at the site above canopy.

Satellite Data
Two different datasets covering the period 2002-2014 from the MODIS sensor onboard the Terra and Aqua satellites were utilized.A MAIAC pre-processed MODIS time series dataset was provided by NASA's Goddard Space Flight Center with 1-2 daily observations from July 2002 to June 2014.The data consist of gridded surface reflectances in band 1-12 for a 3 × 3 pixel window centered at the position of the flux tower at a spatial resolution of 1 km × 1 km.The MAIAC algorithm utilizes time-

Flux Tower Data
The net ecosystem exchange (NEE) was measured using the eddy covariance (EC) technique [45].Shortly, the measurement system, located 23.3 m above the ground, included an ultrasonic anemometer (Solent Research 1012R2, Gill Instruments Ltd., Lymington, Hampshire, UK) for measuring three wind speed components and temperature and a closed-path infrared gas analyzer (LI-6262, LI-COR Biosciences, Lincoln, NE, USA), which measured CO 2 and H 2 O concentrations.The vertical CO 2 fluxes were calculated over 30 min averaging periods by block averaging approach using the EddyUH software [46].The fluxes were corrected for high frequency response underestimation by using the co-spectral transfer functions calculated according to Mammarella et al. [47] together with in-situ parameterization of the co-spectral model derived from ensemble mean of measured temperature co-spectra [47].The CO 2 flux was then corrected for storage change to obtain NEE.The 30 min averaged NEE was partitioned into total ecosystem respiration (TER) and gross primary productivity (GPP) as described in Kolari et al. [48].The air temperature and the incoming photosynthetically active radiation (PAR) were also measured at the site above canopy.

Satellite Data
Two different datasets covering the period 2002-2014 from the MODIS sensor onboard the Terra and Aqua satellites were utilized.A MAIAC pre-processed MODIS time series dataset was provided by NASA's Goddard Space Flight Center with 1-2 daily observations from July 2002 to June 2014.The data consist of gridded surface reflectances in band 1-12 for a 3 × 3 pixel window centered at the position of the flux tower at a spatial resolution of 1 km × 1 km.The MAIAC algorithm utilizes time-series processing of MODIS measurements to take into account different view angles and simultaneously build a dynamic land-water-snow classification and cloud mask to aid in the correction [40,41].MAIAC data for most of the globe are available through NASA GSFC supercomputer and can be obtained from Dr. Lyapustin by request.
In addition, MOD13Q1 250 m spatial resolution NDVI products were acquired for the 3 × 3 pixel window centered on the site.The NDVI product is calculated as where R is the reflectance at the given wavelengths in nm (near infrared and red reflectance).The product dataset is comprised of 16-day composite NDVI values (23 observations per year), which were built using the Constrained View Maximum Value Composite algorithm (CV-MVC).The algorithm accounts for negatively biased noise resulting from weather conditions by selecting the maximum NDVI value for each pixel over a set period and is designed to favor observations closest to nadir view [35].
The size of the footprint of the flux tower varies depending on micrometeorological conditions, but approaches that of a pixel from the MODIS sensor (250 m-1000 m depending on wavebands).In accordance with previous studies [10,44], the single pixel centered on the tower was assessed to be representative of the tower footprint.

Data Pre-Processing and Computation of VIs
Three different VI time series were constructed based on the MODIS satellite data: NDVI maiac and sPRI ref from the MAIAC dataset and NDVI prod from the MOD13Q1 product.The data were initially screened to exclude pixels with erroneous values due to presence of snow, clouds and cloud shadows.Despite the CV-MVC algorithm applied to the MOD13Q1 product, the NDVI prod values still had a considerable level of noise.Quality assurance data provided with the MOD13Q1 product were used to exclude observations that were flagged as lower quality based mainly on clouds, snow and aerosol quantity..The MAIAC dataset was provided pre-screened with approximately 88% of the observations recorded as missing data values and thus the remaining non-erroneous 12% of the data could be used in the analysis.
NDVI maiac was calculated like NDVI prod using Equation (1).Due to the lack of clarity regarding the optimal MODIS reference bands for the PRI, a number of different versions of the index using MODIS band 11 (546-556 nm) were calculated based on bands that have been used in previous studies [26,27,30].These included bands 1 (620-670 nm), 4 (545-565 nm), 10 (483-493 nm) and 12 (546-556 nm), modifying the original PRI formula by Gamon et al. [17]: where R ref indicates the reference band reflectance (i.e., reflectance in bands 1, 4, 10 or 12).To obtain only positive PRI values, which was a requirement in a later stage of the analysis, the PRIs were scaled using the following equation [32]: To determine the most useful reference band for the PRI ref , the MODIS indices were linked temporally with flux tower measurements of GPP and LUE.The satellite overpass times were matched with the nearest time stamp in the flux data.In cases where the time difference exceeded 10 min, the average of the two closest measurements was calculated.The statistical relationships between GPP/LUE and the differently configured PRIs and the NDVI were then explored through regression analysis to evaluate their suitability as indicators of photosynthetic activity.The PRI with the strongest relationship with GPP and LUE (highest R 2 ) was carried forward in the analysis.
Several authors suggest setting a baseline for VI values for the winter months in boreal regions due to persistence of snow and low sun zenith angles which introduces a high degree of noise [38,50].A pixel-specific baseline VI value, assumed to be constant across years, can be used because boreal vegetation goes into a state of stable dormancy during the winter [38].The winter VI values were estimated based on approximate lowest usable values over the full time series within a time window in each year proposed by Beck et al. [38] in the autumn following senescence, but prior to persistent snow cover (mid-October to mid-November).The baseline for NDVI and PRI was set at 0.4 and 0.3 respectively in the months December to February.Any values below the winter baseline VI were replaced with these values.
The final pre-processing steps involved calculating the mean VI on days with multiple overpasses and excluding satellite overpasses occurring before 10:00 a.m. to reduce effects of solar zenith angle variations.Days with no observations were assigned N/A values.

Time Series Processing
Following pre-processing it was necessary to smooth the data further before extracting seasonality metrics.A wide variety of smoothing techniques has been implemented in previous studies, however, there is a considerable lack of clarity in the literature regarding the most appropriate methods to employ [51].The overall objective is to eliminate as much noise as possible, but to preserve variations in the seasonal curves in the data, e.g., by applying moving average filters or fitting functions [37].Four different methods were applied here to allow an assessment of the impact of the chosen smoothing methods on the derived phenology estimates.The selected methods were Gaussian function, double sigmoidal function, fast Fourier transform (FFT) and Savitzky-Golay filter based on their successful implementation in previous studies (Table 1).These were applied to the MODIS VI time-series using the R-package Phenex 1.1-9.Phenex facilitates smoothing of time-series data with irregular time intervals on a single-year basis when assigning dates with no value NA [52].In addition, the best index slope extraction (BISE) algorithm was applied to the NDVI prod data prior to this to further reduce the noise in the data [52,53].Default settings were kept in all cases.The four smoothing techniques were also implemented for the NEE flux data and ground measured PAR (daily mean of NEE and PAR values recorded in the satellite overpass time interval of 10:00 a.m. to 12.30 p.m. and normalized between 0 and 1).In Boreal latitudes, PAR exhibits a strong seasonal pattern due to seasonal changes in sun elevation, which could have an impact on VI estimates.Therefore, PAR was included to test if the sun elevation dynamics correlate with NEE estimates of phenology to evaluate if any potential relationships observed between NEE and the VIs occur as a result of this.

Extraction of Phenology Metrics
A variety of phenology metrics can be derived from time-series VI data.Of interest here was the timing of the start and the end of the growing season (SGS and EGS respectively) given as day of the year (doy).Definitions of SGS and EGS vary significantly across studies, which introduces a degree of subjectivity in the extraction of the parameters [46].The approach taken here was to test two general extraction methods for the flux tower data (only the second method was used to extract phenology metrics for PAR): 1.
Determining the metrics based on the raw flux data according to Tanja et al. [60] and Thum et al. [61].SGS and EGS were defined as the first and the last date at which the daily mean NEE first fall below 20% of the maximum summer uptake of CO 2 defined as the lowest one percentile of NEE values.

2.
Determining the metrics based on the smoothed NEE and PAR data [9] implemented in the "Phenex" package.To ensure consistency, a relative threshold set at 20% of NEE/PAR curve amplitude was used to determine SGS and EGS dates.
The level of agreement of the results generated from the two different methods could then provide an indication of the reliability of the results.Method (2) was applied to the smoothed MODIS NDVI and PRI data with a threshold of 20% of the VI curve amplitude to extract the satellite-based phenology metrics.

Evaluation of Metrics
A correlation analysis was performed, comparing the phenology metrics derived from MODIS VIs and in-situ.Based on the assumption that the metrics estimated from the NEE data can function as ground truth data, the uncertainty associated with MODIS-derived metrics was quantified and assessed by calculating the Spearman rank correlation coefficient (r), the coefficient of determination (R 2 ), the Root Mean Squared Error (RMSE) and bias.

Time Series Smoothing
The NEE and VI time series with fitted curves over the full study period from 2002 until 2014 are displayed in Figure 2.Although four different PRIs were calculated, only the scaled PRI using band 1 (sPRI 1 ) is displayed as this version showed to be most suitable for estimating photosynthetic activity at this site (Table 2), as also demonstrated in other northern forests [39].All three VI time series displayed distinct seasonal variations with values peaking in July-August.NDVI maiac and sPRI 1 showed highly consistent seasonal variation in terms of the shape and magnitude of the smoothed curves regardless of the method applied.Despite the application of the BISE algorithm to the product data, a considerable amount of noise was still present in the NDVI prod time series and consequently there are more pronounced differences among the four smoothed curves, particularly for FFT and Savitzky-Golay (Figure 2d).Although daily means rather than 30-min intervals were used to produce the smoothed NEE curves, the NEE time series (Figure 2a) has a considerably higher temporal resolution compared with the satellite data, allowing more precise definition of the seasonality metrics.

Ground-Derived Metrics
SGS and EGS derived from both the raw and smoothed NEE data as well as the mean of the predicted dates are displayed in Figures 3 and 4. SGS was well defined for each year as the values were relatively consistent based on all methods except the Gaussian smoothing method, which produced slightly different results.SGS dates were mostly between Julian dates 70 and 110; however, in 2014, SGS were predicted to occur earlier than all other years, around Day 60-80.NEE-derived SGS tended to follow the timings of air temperatures exceeding a couple of degrees above 0 °C.An example of this is given in Figure 5, which displays SGS dates predicted from raw NEE values and seasonal variations in air temperatures (5-day moving average) for 2013 and 2014.EGS was less clearly defined from the ground data.The different methods applied produced values that were highly inconsistent relative to the inter-annual variability in the data (Figure 4).No overall long-term trends in the SGS and EGS metrics could be detected in the flux data.

Ground-Derived Metrics
SGS and EGS derived from both the raw and smoothed NEE data as well as the mean of the predicted dates are displayed in Figures 3 and 4. SGS was well defined for each year as the values were relatively consistent based on all methods except the Gaussian smoothing method, which produced slightly different results.SGS dates were mostly between Julian dates 70 and 110; however, in 2014, SGS were predicted to occur earlier than all other years, around Day 60-80.NEE-derived SGS tended to follow the timings of air temperatures exceeding a couple of degrees above 0 • C.An example of this is given in Figure 5, which displays SGS dates predicted from raw NEE values and seasonal variations in air temperatures (5-day moving average) for 2013 and 2014.EGS was less clearly defined from the ground data.The different methods applied produced values that were highly inconsistent relative to the inter-annual variability in the data (Figure 4).No overall long-term trends in the SGS and EGS metrics could be detected in the flux data.
SGS and EGS estimates based on ground-measured PAR varied considerably among the different methods relative to the low inter-annual variability and the dates were generally predicted to occur earlier in the year than those predicted by NEE.

Satellite-Derived Metrics
SGS dates derived from the MODIS VIs are graphed in Figure 3. NDVImaiac and sPRI1 exhibited a relatively high degree of consistency regardless of the methods used.Only in 2002 did the estimates diverge considerably.Since there are no data points in the MAIAC VI dataset before July 2002, SGS for this year was excluded in the following validation of the results.EGS in 2014 was likewise excluded due to lack of data after June 2014.In comparison the MAIAC VIs, NDVIprod produced less reliable results as indicated by the significant variation in the SGS estimates, which in some years diverged highly from the mean ground-estimated SGS.The Savitzky-Golay smoothing method produced unrealistic results when applied to the NDVIprod data (specifically in 2003, 2004, 2012 and 2014).For all three VI time series, EGS estimates varied considerably depending on the methods applied (Figure 4).

VI Predictability of Photosynthetic Activity and LUE
When testing the correlation between temporally linked ground observations and MODIS data, it was found that none of the MAIAC VIs strongly correlated with LUE (R 2 < 0.17 in all cases) (Figure 6) and, with sPRI4 and sPRI12, no significant relationships were found (Table 2).However, the correlation analysis yielded significant results for GPP and all the MAIAC VIs.NDVImaiac and sPRI1 plotted in Figure 7 both showed relatively high R 2 values (R 2 = 0.47 and R 2 = 0.63).Out of the four varieties of PRI tested, sPRI1 had the strongest relationship with both GPP and LUE and thus this version of the index was used in the time series analysis.SGS and EGS estimates based on ground-measured PAR varied considerably among the different methods relative to the low inter-annual variability and the dates were generally predicted to occur earlier in the year than those predicted by NEE.

Satellite-Derived Metrics
SGS dates derived from the MODIS VIs are graphed in Figure 3. NDVI maiac and sPRI 1 exhibited a relatively high degree of consistency regardless of the methods used.Only in 2002 did the estimates diverge considerably.Since there are no data points in the MAIAC VI dataset before July 2002, SGS for this year was excluded in the following validation of the results.EGS in 2014 was likewise excluded due to lack of data after June 2014.In comparison the MAIAC VIs, NDVI prod produced less reliable results as indicated by the significant variation in the SGS estimates, which in some years diverged highly from the mean ground-estimated SGS.The Savitzky-Golay smoothing method produced unrealistic results when applied to the NDVI prod data (specifically in 2003, 2004, 2012 and 2014).For all three VI time series, EGS estimates varied considerably depending on the methods applied (Figure 4).

VI Predictability of Photosynthetic Activity and LUE
When testing the correlation between temporally linked ground observations and MODIS data, it was found that none of the MAIAC VIs strongly correlated with LUE (R 2 < 0.17 in all cases) (Figure 6) and, with sPRI 4 and sPRI 12 , no significant relationships were found (Table 2).However, the correlation analysis yielded significant results for GPP and all the MAIAC VIs.NDVI maiac and sPRI 1 plotted in Figure 7 both showed relatively high R 2 values (R 2 = 0.47 and R 2 = 0.63).Out of the four varieties of PRI tested, sPRI 1 had the strongest relationship with both GPP and LUE and thus this version of the index was used in the time series analysis.

Validation and Inter-Comparison of MODIS-Based Metrics
mean ground-estimated SGS and EGS values estimated from the various techniques were used as validation of the MODIS data, although the ground-estimated EGS dates were considered less reliable due to the high degree of variability in the results.Significant relationships were found between all MODIS VI-predicted SGS dates and ground-derived SGS, except for the FFT and Savitzky-Golay smoothed NDVI prod time series (Table 3).The MAIAC-based metrics yielded considerable higher R 2 values in comparison with NDVI prod .SGS dates estimated using NDVI maiac all strongly correlated with NEE SGS with R 2 values between 0.6 and 0.79 (Table 3).SGS estimates based on the sPRI 1 time series smoothed using a double sigmoidal function exhibited the highest R 2 value (R 2 = 0.8) and had a relatively low RMSE (RMSE = 9.8 days).All except one of the MODIS VIs predicted the SGS earlier in the year than the ground data (the bias values are all negative) (Table 3).Conversely, all MODIS EGS dates were estimated to be later in time than that predicted from the ground-measured NEE.PAR did not correlate significantly with the NEE phenology estimates using any of the four methods.Likewise, no significant correlations were found between any of the satellite and ground-predicted EGS dates (Table 4).

Discussion
The study set out to achieve three objectives: to estimate phenological parameters from MODIS VI time series, to compare the utility of MAIAC-processed data against standard product data and to assess the predictability of seasonal transitions and vegetation productivity from MODIS PRI.Numerous studies have employed similar methods linking NDVI time series with flux tower measurements [9,44]; however, there is currently very limited research published on the use of satellite-derived PRI for estimating phenology.

Performances of MAIAC VIs and NDVI prod
Overall, relatively high correlations were found between the MAIAC VI time series and the ground estimated SGS (R 2 = 0.36-0.80).Few studies have found that MODIS VIs provide good estimates of phenology at a single site in evergreen conifer forests.Böttcher et al. [10] found significant correlations between NDVI derived from daily MODIS reflectances and CO 2 fluxes, however, Balzarolo et al. [9] and Liu et al. [62] concluded that the inter-annual variations in the timings of the seasonal transitions in evergreen forest ecosystems were too low relative to the error in the MODIS NDVI time series to accurately predict phenological parameters.The noise introduced by varying viewing and illumination conditions, atmospheric conditions and snow cover becomes too significant to detect the phenological variations.A significant challenge lies in the tendency for photosynthetic activation of boreal coniferous forests to occur in the spring while there is still snow on the ground, affecting the VI observations and thus their utility as predictors of phenology.Furthermore, frequent cloud cover in boreal regions can cause large temporal gaps in the data.Some of these issues may be cause of the weak correlations found between the SGS estimated from NDVI prod and the ground data.The processing of the MOD13Q1 product using the CV-MVC algorithm attempts to correct for some of the effects mentioned, however, the data still exhibited a high degree of noise in comparison with the MAIAC data.The product data are also limited by their coarse temporal resolution.Since each sample is taken at the maximum value within a 16-day interval, the time gap between two samples can be up to 31 days, which is significant in this context as spring reactivation of photosynthetic activity occurs rapidly in evergreen conifers [14].The significant relationships observed between SGS derived from the MAIAC data and in situ are thus remarkable and they indicate that the MAIAC aerosol-surface algorithm works effectively.This is also evident from the fitted curves in Figure 2, which exhibit steady seasonal curves and low degree of noise.The algorithm takes into account many of the factors that introduce error in the VI data by deriving average bidirectional reflectance factors [41], retrieving aerosol information [40] and by removing cloud-contaminated pixels through time series processing [63].Although the application of the MAIAC algorithm to the MODIS data had effectively removed 88% of the observations in the study period (i.e., observations were recorded as missing data), the MAIAC dataset had a markedly better temporal resolution during the growing season compared to the 16-day composite NDVI product.Consequently, interpolation between data points in the NDVI maiac and sPRI 1 time series was less susceptible to biases.
Despite the significant relationships observed, the SGS dates were generally predicted to be earlier in the year for all three MODIS time series compared to SGS estimated from the ground data.This was also found in previous MODIS-based phenology studies by Garrity et al. [57], Hmimina et al. [11] and Balzarolo et al. [9].The reduced number of useful observations during the winter months tends to extend into the beginning of the green-up phase.Hmimina et al. [11] found that gaps at the transitional stage in an NDVI time series smoothed with a sigmoidal function, shifted green-up estimates towards earlier in the year.The MAIAC dataset had gaps in the data at the beginning of the transition phase in all years (between the winter baseline and the first non-negative MODIS observations), which may be the cause of the early predictions.Furthermore, defining the phenological metrics is somewhat subjective and attempts to describe a gradual temporal transition as a single date.Therefore, some time lag between the estimates is not surprising considering that the same methods have been applied to two types of data based on fundamentally different observation techniques.The 20% threshold was selected according to Tanja et al. [60] and Thum et al. [61], who found this to produce similar predictions to those of temperature data, however, relative thresholds used in the reviewed literature vary between 10% [62] and 50% [9].Derivative methods to identify the points of maximum increase on the seasonal curve have also been applied in previous studies (e.g., Tan et al. [54]), which is particularly beneficial when applied to regional and global scale studies and ecosystems characterized by multiple growth cycles [55].
EGS could not be successfully estimated using any of the MODIS time series.Previous studies have similarly found that EGS could not be determined with as high accuracy as SGS [11,62].The factors controlling the transition of vegetation into winter dormancy are not as well understood as for the onset of greenness during spring [64].However, the decline in the VI signal taking place during senescence tend to occur at a slower rate and is less pronounced to the spring transition [62], which may explain why it was difficult to establish a specific EGS date from both the ground data and the remote sensing data.Furthermore, an important factor that causes GPP to decline in the autumn at high latitudes is the significant decrease in the amount of available PAR.The importance of LUE in vegetation productivity is thus decreased during this season as there is insufficient PAR to drive photosynthesis.

PRI as an Indicator of Photosynthetic Activity and Phenology
The initial testing of the ability of the different PRIs to describe variations in photosynthetic activity and LUE revealed mixed results.While previous studies have found significant relationships between MODIS PRI and LUE [26,28], none of the PRIs correlated well with LUE in this case (Figure 6).A number of factors could contribute in obscuring the relationships observed between satellite observations and flux tower measurements.At the ecosystem scale, variations in canopy structure (e.g., leaf orientation), background effects and discrepancies in the spatial footprint of flux measuring station and the MODIS sensor can affect the PRI signal significantly [33].The estimation of seasonal phenological dates from optical remote sensing observations of VIs in boreal regions may also be attributed to the high number of missing observations due to cloud cover, and the occurrence of snow at the same time as the vegetation is actually active.
An important source of uncertainty also lies in the calculation of LUE, which was based on ground-measured GPP and PAR.While a number of different definitions of LUE have been expressed in the literature, the original model of LUE defined by Monteith [49] utilizes the amount of incident photosynthetically active light absorbed by the canopy (APAR).However, lack of data on the fraction of radiation absorbed by green vegetation meant that APAR could not be included.Gitelson and Gamon [65] investigated the implications of using different definitions of LUE.They found that the definition of LUE based on either PAR, total absorbed PAR and PAR absorbed by photosynthetically active radiation had an impact on the estimations of seasonal vegetation productivity as well as the temporal behavior of LUE.Using APAR, rather than simply incident radiation could thus potentially produce different PRI-LUE relationships than observed here.
While no strong correlation was found between PRI and LUE, sPRI 1 , which uses the broad red band as a reference band, yielded the highest R 2 .Results from studies by Goerner et al. [29] and Drolet et al. [26] indicate that red wavelengths (MODIS bands 1, 12 and 13) generally tend to be useful as reference bands for estimating vegetation LUE from MODIS PRI.sPRI 1 was also the only version of PRI, which correlated strongly with ground-measured GPP.It is noteworthy that sPRI 1 had a stronger relationship with GPP (R 2 = 0.63) compared with NDVI (R 2 = 0.47), which is more commonly used to observe vegetation productivity.Further other studies have found that NDVI and similar VIs that are related to leaf biomass do not capture the seasonal variation in LUE and photosynthetic activity as well as PRI in evergreen forests [66][67][68].
The strong relationships observed between SGS derived from NEE and sPRI 1 (when applying Gaussian and double logistic functions) suggest that PRI could serve as an effective indicator of seasonal transitions.This is in agreement with the findings of Wong and Gamon [14] who tested the use of PRI for estimating spring phenology in conifers at the leaf and canopy scale.Wong and Gamon [14] found that the seasonal increase in PRI was nearly coincident with spring activation of photosynthetic activity with only a couple of days lag.The authors could conclude that the main physiological process accounting for the response of PRI observed was the seasonal change in pigment pool sizes.While the highest R 2 in this study was found for sPRI 1 , the NDVI maiac predictions of SGS exhibited more consistently high correlations with the ground estimates (R 2 > 0.6 in all cases).The results can thus not confirm that PRI is more useful for predicting phenological events in conifers compared to NDVI.Nevertheless, it is significant that PRI shows potential for detecting physiological transitions in evergreen conifer forests as the seasonal pattern of photosynthetic activity is mainly driven by air temperature [43,69], which makes boreal vegetation highly susceptible to climate change.

Uncertainties
Selecting the most suitable data smoothing and phenology extraction methods is not straightforward due to the variety of recommendations expressed in literature regarding the best performing techniques [37,38,50].The reliability of the methods can be difficult to assess and Ekhlund and Jönnson [51] emphasize that smoothing processes do not always result in an accurate depiction of seasonal variations in photosynthetic activity, particularly when the data are very noisy.For instance, the Savitzky-Golay filter did not produce reliable results for the NDVI prod time series, which is consistent with previous studies that have demonstrated the sensitivity of this method to noise [36,50].Applying four different methods provided some insight into the dependability of the phenology estimates.A comparative study by White et al. [5] tested 10 different methods for estimating spring phenology and found that different techniques produced significantly different results.Likewise, in this study, the choice of filtering or fitting technique did indeed influence the estimates and none of the methods applied were consistently superior at this site.However, the less noisy time series (NEE, NDVI maiac and sPRI 1 ) were less impacted by the choice of smoothing method (Figures 3 and 4) and the long study period of 13 years provided enough data points to observe the relationships between SGS predicted by the MAIAC time series and the NEE data with a fair degree of certainty as indicated by the correlation coefficients (Table 3).
The ability of NDVI maiac and sPRI 1 to estimate timings of spring seasonal transitions at the study site is very promising, however, the results should be interpreted with caution.It is important to note that the seasonal signal identified from the VI time series could occur as a results of other factors than photosynthetic activity [5].At boreal forest sites, the seasonal cycle in vegetation activity is highly consistent, driven by temperature and to some extent by available PAR.It can therefore be observed in several co-correlating variables, including snow cover [61] and any variations in the timings of the growing seasons identified from satellite-derived VIs could indicate temperature changes that result in green-up and senescence, but do not necessarily directly track increases or decreases in photosynthetic activity.The impact of sun elevation was tested by examining the correlation between the annual estimates of phenology based on NEE and PAR.The fact that no significant relationships were found between these demonstrates that the NEE-VI estimates are not an artefict of inter-annual variations in PAR.Other measures should be taken to distinguish the effects in snow cover on the VI signal.For instance, Delbart et al. [70] developed a snow detection algorithm specifically for forest ecosystems, which could potentially be applied to the data.
To gain a better understanding of the factors that complicate comparison of the satellite data and eddy covariance flux data, in situ radiation measurements have been incorporated in previous studies [9,44].Obtaining ground-measured PRI and NDVI would help bridge the gap between the two different types of observations and thus provide a more direct assessment of the accuracy of the MODIS-estimated timings of seasonal transitions.

Conclusions
By linking ground measurements of NEE with MODIS data, it was demonstrated that VI time series have significant potential for predicting the onset of the growing season in evergreen conifers.The work presented in this paper is based on potentially the longest PRI time series in satellite-based vegetation studies to date, made possible by the MAIAC-processed MODIS data.While no long-term greening trends were identified at the study site over 13 years, the MAIAC-derived NDVI and PRI allowed detection of inter-annual variations in SGS dates, which could not be achieved using the standard MODIS NDVI product.The range of results generated from each dataset highlights some inherent issues with the lack of reliable, standardized methods for deriving phenology from satellite-derived VIs and flux data.The MODIS data could not describe the timings of senescence accurately since the magnitude of in the estimated EGS dates was of the same order as the magnitude of inter-annual variation, even in the ground data.
The results of the study are encouraging and provide support for additional work to confirm the relationships observed and to further explore the usefulness of MODIS PRI for estimating phenological events.It would be interesting to see how well PRI performs in a wider range of vegetation types and ecosystems, even if the value of the index particularly lies in the potential for improving estimates in evergreen forests.Future work could involve processing of spatially continuous MAIAC time series data to study regional patterns of phenology.Pixel-based analysis of time series of images is for instance facilitated by TIMESAT software developed by Jönsson and Ekhlund [36].The program is currently limited to data with even time intervals, but efforts are being made to develop new ways to incorporate irregular time series as well as extending the smoothing techniques into the spatial domain [51].
The study demonstrates that there are significant advantages in utilizing MAIAC-processed data for detecting phenology, in terms of both the improved noise reduction and the possibility of studying long-term PRI time-series.The use of satellite-derived PRI for detecting phenology could potentially enable systematic monitoring of changing patterns in the seasonality of evergreen vegetation and thus reduce significant uncertainties in global carbon budgets associated with boreal regions.

Figure 2 .
Figure 2. Times series of flux tower measured scaled Net Ecosystem Exchange (NEE) (a) and MODIS vegetation indices smoothed according to four different methods: Application of Gaussian functions, double sigmoidal functions, fast Fourier transform (FFT) and a Savitzky-Golay filter.The MODIS time series are: MAIAC-processed scaled PRI (b); and NDVI data (c); and MODIS standard NDVI product MOD13Q1 (d) and (e) scaled PAR.

Figure 2 .
Figure 2. Times series of flux tower measured scaled Net Ecosystem Exchange (NEE) (a) and MODIS vegetation indices smoothed according to four different methods: Application of Gaussian functions, double sigmoidal functions, fast Fourier transform (FFT) and a Savitzky-Golay filter.The MODIS time series are: MAIAC-processed scaled PRI (b); and NDVI data (c); and MODIS standard NDVI product MOD13Q1 (d) and (e) scaled PAR.

Figure 3 .
Figure 3.Estimated start of growing season (SGS) based on daily Net Ecosystem Productivity (NEE), the three different MODIS vegetation index time series (NDVImaiac, sPRI1 and NDVIprod) and groundmeasured photosynthetically active radiation (PAR).

Figure 3 .
Figure 3.Estimated start of growing season (SGS) based on daily Net Ecosystem Productivity (NEE), the three different MODIS vegetation index time series (NDVI maiac , sPRI 1 and NDVI prod ) and ground-measured photosynthetically active radiation (PAR).

Figure 4 .
Figure 4.Estimated end of growing season (EGS) based on raw and smoothed daily Net Ecosystem Productivity (NEE) data, the three different MODIS vegetation index time series (NDVImaiac, sPRI1 and NDVIprod) and ground-measured photosynthetically active radiation (PAR).

Figure 4 .
Figure 4.Estimated end of growing season (EGS) based on raw and smoothed daily Net Ecosystem Productivity (NEE) data, the three different MODIS vegetation index time series (NDVI maiac , sPRI 1 and NDVI prod ) and ground-measured photosynthetically active radiation (PAR).

Figure 5 .
Figure 5. Daily mean air temperature measurements (five-day moving averages) at the Hyytiälä Forestry Field Station for 2013-2014.The start of the growing season (SGS) dates predicted from the NEE data occur almost immediately after air temperatures exceed 2-3 °C.

Figure 5 .
Figure 5. Daily mean air temperature measurements (five-day moving averages) at the Hyytiälä Forestry Field Station for 2013-2014.The start of the growing season (SGS) dates predicted from the NEE data occur almost immediately after air temperatures exceed 2-3 • C.

Figure 6 .
Figure 6.The relationships between the four different PRIs and flux tower measured light use efficiency (LUE).The red lines are a simple linear regression lines.

Figure 7 .
Figure 7.The relationship between NDVImaiac and sPRI1 and flux tower measured gross primary productivity (GPP).The red lines are simple linear regression lines.

Figure 6 . 21 Figure 6 .
Figure 6.The relationships between the four different PRIs and flux tower measured light use efficiency (LUE).The red lines are a simple linear regression lines.

Figure 7 .
Figure 7.The relationship between NDVImaiac and sPRI1 and flux tower measured gross primary productivity (GPP).The red lines are simple linear regression lines.

Figure 7 .
Figure 7.The relationship between NDVI maiac and sPRI 1 and flux tower measured gross primary productivity (GPP).The red lines are simple linear regression lines.

Table 2 .
Statistics of correlation analysis between GPP/LUE and NDVI and four different PRIs.The highest R 2 values observed are highlighted in bold and insignificant relationships are marked in grey.R 2 : Coefficient of determination, r: correlation coefficient (rho from Spearman Rank's test significant at p < 0.05).Sample size = 766.

Table 2 .
Statistics of correlation analysis between GPP/LUE and NDVI and four different PRIs.The highest R 2 values observed are highlighted in bold and insignificant relationships are marked in grey.R 2: Coefficient of determination, r: correlation coefficient (rho from Spearman Rank's test significant at p < 0.05).Sample size = 766.

Table 3 .
Statistics of correlation between the start of the growing season (SGS) derived from ground measured NEE and SGS derived from MODIS VI time series and PAR.R 2 : Coefficient of determination; r: rho from Spearman Rank's correlation; RMSE: Root mean squared error.R 2 > 0.7 are highlighted in bold.Insignificant relationships are marked in grey (p-value > 0.05).

Table 4 .
Statistics of correlation between the end of the growing season (EGS) derived from ground measured NEE and EGS derived from MODIS VI time series and PAR.R 2 : Coefficient of determination; r: rho from Spearman Rank's correlation; RMSE: Root mean squared error.All tested relationships were insignificant (p-value > 0.05).