Evaluation of the Plant Phenology Index (PPI), NDVI and EVI for Start-of-Season Trend Analysis of the Northern Hemisphere Boreal Zone

Abstract

Satellite remote sensing of plant phenology provides an important indicator of climate change. However, start of the growing season (SOS) estimates in Northern Hemisphere boreal forest areas are known to be challenged by the presence of seasonal snow cover and limited seasonality in the greenness signal for evergreen needleleaf forests, which can both bias and impede trend estimates of SOS. The newly developed Plant Phenology Index (PPI) was specifically designed to overcome both problems. Here we use Moderate Resolution Imaging Spectroradiometer (MODIS) data (2000–2014) to analyze the ability of PPI for estimating start of season (SOS) in boreal regions of the Northern Hemisphere, in comparison to two other widely applied indices for SOS retrieval: the Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI). Satellite-based SOS is evaluated against gross primary production (GPP)-retrieved SOS derived from a network of flux tower observations in boreal areas (a total of 81 site-years analyzed). Spatiotemporal relationships between SOS derived from PPI, EVI and NDVI are furthermore studied for different boreal land cover types and regions. The overall correlation between SOS derived from VIs and ground measurements was rather low, but PPI performed significantly better (r = 0.50, p < 0.01) than EVI and NDVI which both showed a very poor correlation (r = 0.11, p = 0. 16 and r = 0.08, p = 0.24). PPI, EVI and NDVI overall produce similar trends in SOS for the Northern Hemisphere showing an advance in SOS towards earlier dates (0.28, 0.23 and 0.26 days/year), but a pronounced difference in trend estimates between PPI and EVI/NDVI is observed for different land cover types. Deciduous needleleaf forest is characterized by the largest advance in SOS when considering all indices, yet PPI showed less dramatic changes as compared to EVI/NDVI (0.47 days/year as compared to 0.62 and 0.74). PPI SOS trends were found to be higher for deciduous broadleaf forests and savannas (0.54 and 0.56 days/year). Taken together, the findings of this study suggest improved performance of PPI over NDVI and EVI in retrieval of SOS in boreal regions and precautions must be taken when interpreting spatio-temporal patterns of SOS from the latter two indices.

1. Introduction

Growing consensus on severity of the ecological impacts posed by recent global warming and its amplification in the high latitudes of the Northern Hemisphere (NH) increases the importance of understanding the interaction between the climatic changes and terrestrial ecosystems [1,2,3,4,5]. Plant cyclical biological events (plant phenology) are strong indicators of the seasonality of the environment, revealing the implications of rising temperatures on vegetation functioning [6,7,8,9]. Depending on the vegetation cover and climate zone, plant phenology responds differently to climatic factors, such as air temperature, precipitation and photoperiod [10]. However, altered temperature conditions have been observed to have a dominating influence on vegetative activity in the Northern Hemisphere, which makes changes in plant phenology a significant indicator of biological responses to increasing temperatures [1,11,12,13,14].

Not only is plant phenology crucial for understanding the impacts of climate change, but phenology is also controlling various processes in the terrestrial ecosystems, such as net primary productivity and seasonal exchanges of CO₂, water and energy between the land surface and the atmosphere [15,16,17,18]. It has been observed that increasing temperatures in the NH are associated with advanced onset of growing season, which may possibly result in longer growing season and consequently alteration of the annual carbon (C) cycle [19,20]. Boreal forests—the dominating biome of the NH, are representing 32% (272 ± 23 PgC) of the global forest C stock, hence playing a crucial role in the global C cycle, either as a sink or source [21,22]. Investigating phenology of boreal forests is therefore important for understanding the impacts of climate change and associated changes in global C balance [23].

Satellite-based observations are becoming essential for boreal region plant phenology retrieval, since it allows for remotely capturing phenological variations both at small and large scales, while ground observations of the Northern mid and high latitude plant phenology are limited and represent local conditions [7,16,23,24,25,26]. Existing satellite-based research in plant phenology over the boreal zone of the NH has been mainly carried out by employing spectral vegetation indices (VI), calculated from the reflectance acquired from optical remote sensing using satellites [27,28]. The most widely applied VI that combines visible and near infrared light reflected by vegetation is the Normalized Difference Vegetation Index (NDVI) [29,30]. NDVI is considered to a proxy of plant phenology, since it observes vegetation canopy greenness at the same time as minimizing clouds, changing sun angles, topography and atmospheric contamination effects. However, it has been reported to be influenced by canopy background, atmospheric conditions and reduced sensitivity to variations in high biomass regions [31,32]. In order to improve upon these drawbacks, an Enhanced Vegetation Index (EVI) was developed by incorporating soil adjustment and atmospheric resistance factors that minimize soil brightness influences, increase sensitivity to variations in dense canopies and further reduce atmospheric influences [33].

However, existing studies based on NDVI and EVI reveal that existing methods for retrieval of boreal region plant phenology often involve uncertainties [34,35,36,37], leading to inconsistency among findings of satellite-based boreal forest phenology studies in relation to ground observations that have shown a continuous advance in the onset of the growing season both in North America and Eurasia [6,38,39]. The sources of uncertainty are mainly attributed to insufficient methods for overcoming the presence of snow that significantly affects VI signal [16,40,41,42,43], as well as difficulties in capturing subtle seasonal variation in canopy greenness of evergreen forests [15,44,45,46,47].

A physically-based Plant Phenology Index (PPI) was recently developed [36] to improve the satellite-based observations of vegetation phenology for the boreal regions. The index is derived from radiative transfer equations and is calculated from red and near-infrared (NIR) reflectance. The fundamental property of the PPI is its linearity to green leaf area index (LAI) and, unlike NDVI and EVI, PPI is reported to be well suited for phenology retrieval in boreal forests, due to its ability to minimize the snow influence on vegetation signal and increase the sensitivity to minor seasonal variation in dense canopy greenness of the boreal forests [36]. To date, PPI has been applied in a single study that investigated plant phenology trends over Northern Europe for the period of 2000–2014 [48], showing that trends derived from PPI were significantly correlated to phenology ground observations and flux tower estimated gross primary productivity (GPP).

Start of growing season (SOS) of the deciduous broadleaf and deciduous needleleaf forests corresponds to budburst and leaf emergence, whereas SOS of evergreen needleleaf forests is defined as the point in time when photosynthetic activity starts to increase, since evergreen forests retain green foliage all year round [15]. Research objectives of this study focus on comparing trend estimates of boreal zone SOS as derived from PPI in relation to the two widely applied indices: NDVI and EVI. Satellite-based SOS was evaluated against SOS derived from in situ measurements of GPP to test the performance of the VIs. The analysis focuses on the start of photosynthetic activity in spring SOS since onset of growth is known to be more distinct as compared to autumn phenology, which is often difficult to monitor and reveals scattered results [44,49,50].

2. Study Area

The study area is covering the boreal forest biome of North America and Eurasia (45°–75°N), delineated based on the Terrestrial Ecoregions of the World (TEOW) dataset [51] (Figure 1A). TEOW has been defined based on existing global maps of flora and fauna distribution [51]. The dataset was modified, revised and published by The Nature Conservancy (TNC) in 2009. It is considered to be a comprehensive map that provides a useful framework for conducting both global and regional scale studies [52].

Most of the boreal region is covered by needleleaf forests with a large distribution of evergreen needleleaf forests across the entire study area and deciduous needleleaf forests in the East Siberian Taiga (Figure 1B). The areas at the high latitudes of the study area are mainly covered by open shrublands and woody savannas. The southern latitudes are dominated by deciduous broadleaf and mixed forests. Even though tree cover is rather dense in the boreal region, a large part of terrain is covered with wetlands and lakes. The largest areas, classified as permanent wetlands, are found in the West Siberian Taiga and Hudson Plains of Canada. Moreover, since boreal forests are extending throughout both continuous and discontinuous permafrost zone, the melting active layer (uppermost layer of the permafrost) in the summer period yields swampy, poorly-drained soils across a substantial part of the boreal biome [53].

3. Materials and Methods

3.1. Satellite Data

The MODIS Nadir bidirectional reflectance distribution function (BRDF)—Adjusted Reflectance (NBAR) product was used for calculating VIs [54]. The bands used for calculation of VIs include near infrared (NIR)—band 2 (841–876 nm), red—band 1 (620–670 nm), blue—band 3 (459–479 nm), and, additionally, a zenith angle at local solar noon from MCD43C4 is used for calculation of PPI. MODIS NBAR data is widely used for vegetation analysis and is particularly suitable for PPI, which is a non-ratio based spectral index that requires standardized reflectance corrected for sun-sensor geometry [36]. We used two products of different spatial resolution for the analysis. Due to the large area studied, it was chosen to use the MCD43C4 (collection 5) product that has a 0.05° spatial resolution projected to a latitude/longitude Climate Modeling Grid (CMG). For the evaluation of VIs against flux tower data, we used the MCD43A4 (collection 5) with a 500 m spatial resolution for a better match with the spatial footprint of the flux towers. The average of a window of 3 × 3 pixels encompassing the flux site was used. Both products are characterized by a temporal resolution of eight days from continuous 16-day observations [54].

In order to account for possible influence from forest loss due to forest fires/logging on phenology trend assessment, the high-resolution dataset of global forest cover change by [55] was used in this analysis to mask out the areas classified as disturbed.

3.2. Flux Tower Data

Photosynthetic activity expressed through GPP was used in order to evaluate and validate the satellite-based SOS estimates in this study. Daily GPP were acquired from a SUBSET Data Product for the FLUXNET 2015 (fluxnet.fluxdata.org, 2016). FLUXNET is a global network of micrometeorological flux measurement sites measuring net CO₂ exchange using the eddy covariance method. The FLUXNET 2015 dataset is a new version that offers improved measurements in comparison to the previous versions—FLUXNET Marconi Dataset (2000) and the FLUXNET LaThuile Dataset (2007). The improvements consist of refined data quality, new methods for uncertainty quantification and use of reanalysis data to fill long gaps of micrometeorological variable records (http://fluxnet.fluxdata.org, 2016). In this study, GPP estimates from 14 sites with a total of 81 site-years across the study area were used (Figure 1B; Appendix A

Table A1. Location, data availability and land cover class according the International Geosphere Biosphere Project (IGBP) of the different Fluxnet sites.

Site ID: Site Name	Temporal Data Availability	Location (Lat; Long)	IGBP Vegetation Cover
FI-Hyy: Hyytiala	2000–2014	61.8475;	Evergreen needleleaf (ENF)
		24.2950
RU-SkP: Spasskaya Pad larch	2012–2014	62.2550;	Deciduous needleleaf (DNF)
		129.1680
RU-Che: Cherskii	2002–2005	68.6130;	Permanent wetlands (WET)
		161.3414
CA-Qfo: Quebec-Eastern Boreal, Mature Black Spruce	2003–2010	49.6925;	Evergreen needleleaf (ENF)
		−74.3421
CA-Gro: Ontario-Groundhog River, Boreal mixed wood Forest	2003–2014	48.2167;	Mixed Forest (MF)
		−82.1556
CA-NS1: UCI-1850 burn site	2002–2005	55.8792;	Evergreen needleleaf (ENF)
		−98.4839
CA-NS3:UCI-1964 burn site	2001–2005	55.9117;	Evergreen needleleaf (ENF)
		−98.3822
CA-NS4: UCI-1964 burn site wet	2002–2005	55.9117;	Evergreen needleleaf (ENF)
		−98.3822
CA-NS5: UCI-1964 burn site	2001–2005	55.8631;	Evergreen needleleaf (ENF)
		−98.4850
CA-NS6: UCI-1998 burn site	2001–2005	55.9167;	Open shrublands (OSH)
		−98.9644
CA-NS7:UCI-1998 burn site	2002–2005	56.6358;	Open shrublands (OSH)
		−99.9483
CA-SF1: Saskatchewan-Western Boreal	2001–2006	54.4850;	Evergreen needleleaf (ENF)
		−105.8176
CA-SF2: Saskatchewan-Western Boreal	2001–2005	54.2539;	Evergreen needleleaf (ENF)
		−105.8775
CA-SF3: Saskatchewan-Western Boreal, forest burned in 1998	2001–2006	54.0916;	Open shrublands (OSH)
		−106.0053

). For information of time span of the different sites, we refer to Table A1.

3.3. Vegetation Indices

PPI is a physically-based spectral index for retrieving plant phenology and is formulated to have a nearly linear relationship to green leaf area index (LAI) given a fixed soil effect [36]. The expression of PPI employs red and NIR reflectance and is based on modified Beer’s law that describes the relationship between canopy reflectance and leaf area index (LAI) [36]. PPI (unit: m² m⁻²) is defined as:

$$PPI = -K \times ln\frac{\left(M – DVI\right)}{\left(M – DV{I}_s\right)},$$

(1)

where $DVI$ (Difference Vegetation Index) is the difference between NIR and red reflectance; $M$ is a site specific canopy maximum $DVI$; $DV{I}_s$ is the $DVI$ of the soil. $K$ is the gain factor that is formed as a function of sun zenith angle $\mathsf{\theta }$, a geometric function of leaf angular distribution and instantaneous diffuse fraction of solar radiation [36,56]. $K$ is formulated as:

$$\mathrm{K}=\frac{0.25\cos \left(\mathsf{\theta }\right)}{\left(1-{\mathrm{d}}_{\mathrm{c}}\right)\mathrm{G}+{\mathrm{d}}_{\mathrm{c}}\cos \left(\mathsf{\theta }\right)}\frac{1+\mathrm{M}}{1-\mathrm{M}},$$

(2)

where $\mathsf{\theta }$ is a sun zenith angle; ${\mathrm{d}}_{\mathrm{c}}$ is an instantaneous diffuse fraction of solar radiation for clear sky and standard atmosphere, when the sun is at zenith angle $\mathsf{\theta }$; $\mathrm{G}$ is a geometric function of leaf angular distribution; and M is a site-specific maximum $DVI$.

In this study, the inputs in the $PPI$ formula were calculated as follows:

• M—a maximum difference between NIR and red MODIS NBAR was calculated for each pixel, for the time period 2000–2014.
• DVI value was set to 0.09, which is an empirically defined value [36]. DVIs can cause negative values if canopy DVI is lower than $DV{I}_s$, which can occur during snowy conditions or for an ecosystem with very bright soil. The empirical value of $DV{I}_s$ was still considered to be applicable in this study, since negative $PPI$ estimates will not alter linearity between PPI and LAI [36].
• zenith angle (θ) value for each pixel location was acquired from “Zenith Angle at local solar noon” data layer in the MCD43C4 products.
• d_c has been calculated for standard atmosphere from the equation formulated by [36] following [57,58]:

  $$d_c = 0.0336 + 0.0477/cos\left(\theta \right).$$

  (3)
• G value was set to 0.5 based on [59] who found that for a spherical needle (leaf) orientation G is 0.5, regardless of directions of the solar beam (different view directions). According to [36], G 0.5 is also valid for a flat leaf, when assuming a spherical (uniform) leaf angle distribution and may therefore be used in the calculation of PPI.

For a more complete description of the physical basis of PPI, we refer to [36].

NDVI, developed by [30], is a normalized ratio of the NIR and red bands. NDVI is defined as:

$$\mathrm{NDVI}=\frac{{\mathsf{\rho }}_{\mathrm{NIR}}-{\mathsf{\rho }}_{\mathrm{red}}}{{\mathsf{\rho }}_{\mathrm{NIR}}+{\mathsf{\rho }}_{\mathrm{red}}},$$

(4)

where ρ_NIR and ρ_red are the surface bidirectional reflectance for their respective MODIS bands. Enhanced Vegetation Index (EVI) is an index developed to improve the vegetation signal by incorporating the use of a soil adjustment factor, which ensures that the vegetation lines in the red-NIR space converts to a common point and a blue band that reduces the atmospheric influences [33]. EVI is defined as:

$$\mathrm{EVI}=\mathrm{G}\frac{{\mathsf{\rho }}_{\mathrm{NIR}}-{\mathsf{\rho }}_{\mathrm{red}}}{{\mathsf{\rho }}_{\mathrm{NIR}}+{\mathrm{C}}_1\times {\mathsf{\rho }}_{\mathrm{red}}-{\mathrm{C}}_2\times {\mathsf{\rho }}_{\mathrm{blue}}+\mathrm{L}},$$

(5)

where ρ_blue is surface reflectance for the blue band, L is the canopy background adjustment factor, and C₁ and C₂ are the coefficients of the aerosol resistance term, which uses the blue band in diminishing the atmospheric aerosol scattering in the red band [33]. The coefficients employed in the EVI algorithm are empirically determined, L = 1, C₁ = 6, C₂ = 7.5, and G (gain factor) = 2.5 [33].

3.4. Start of Season (SOS) Retrieval

Start of season (SOS) for VI and GPP data were estimated using the TIMESAT software [60]. By this tool, gap-filling and smoothing of original raw time-series was conducted using an adaptive Savitzky–Golay filtering technique in the SOS retrievals for the VI-GPP SOS comparison, whereas a double logistic smoothing was used for the SOS estimates across the NH. Adaptive Savitzky–Golay filtering uses local polynomial functions and can be referred to a moving average method that smooths the data by replacing the noisy values by the average of the values in the defined window. This method fits close to the raw time-series and should hence only be applied when raw data is relatively clean. When looking at local sites, we could empirically control each time-series and see that the smoothing technique performs well. The double logistic method smooth curves of the seasonality based on intervals around maxima and minima in the time-series, resulting in less sensitivity to noise. This method was reported to perform best on high-latitude vegetation time-series that are often noise contaminated caused by long periods of snow cover and persistent cloud cover [7,60,61]. The following parameters were used within the fitting: Savitzky–Golay: seasonal parameter = 1, number of envelope iterations = 2, adaptation strength = 3, Savitzky–Golay window size = 4; Double logistic: seasonal parameter = 1, number of envelope iterations = 2, adaptation strength = 2. To remove remaining outliers from the VI time series, a spike method of median filter was employed for both methods. A value was defined as an outlier if it diverges from the median of a moving window by more than the spike value multiplied with the standard deviation [62] (standard deviation was set to the default value of 2). A function force to minimum was used to prevent SOS assessment to be governed by an increase in VI values due to snowmelt instead of an increase in vegetation activity. This function enables disregarding the low, vegetation-unrelated winter values that otherwise affect the retrieval of phenological metrics [63]. We used the quality assessment dataset to define the influence (weight) of each pixel in the function fitting to decrease the influence from low quality pixels in capturing SOS. We used the quality assessment (QA) dataset “percent snow” to define the influence (weight) of each pixel in the function fitting to decrease the snow influence in capturing SOS. Pixels being partly/fully snow covered were assigned a weight of 0.1, missing values that were interpolated were given the weight of 0.5, and all of the other pixels were given full weight (1).

Within TIMESAT, SOS can be estimated in two ways [62]; the first method indicates SOS where the fitted curve reaches a % threshold of the seasonal amplitude, whereas the second method indicates SOS where a fitted curve reaches a specific absolute threshold value [60]. Using a % threshold allows automatic adjustment according to local conditions, and this method was chosen for the analysis across the NH. The analysis across the NH involves diverse vegetation and land cover types that does not allow for a unifying threshold value for determining SOS. The % threshold was set to 20% for all VIs [63,64,65]. For the evaluation of VI SOS against GPP SOS data, site-specific absolute threshold values were selected to allow for optimal evaluation of the performance of each of the VIs. When analysing local sites, we could empirically choose an optimum absolute threshold value that allowed for extracting the exact point in time when the phenology signal starts to increase from its base level.

3.5. Statistical Analysis

Trends were estimated by using a non-parametric median trend Theil–Sen (TS) operator [66,67] that calculates the median of the slopes between all n(n − 1)/2 pair wise combinations over time and expresses the trend as a rate of change per time unit, in this study per year. The Theil–Sen median trend was chosen as it is based on nonparametric statistics and is particularly effective for the estimation of trends in short and/or noisy time-series [68,69].

The significance of the time series trends was calculated by the non-parametric Mann–Kendall (MK) significance test. The MK significance test is commonly used as a trend test for the TS median slope operator [68] and produces outputs of z-scores, which allows for the assessment of both the significance and direction of the trend. A positive slope (z ≥ 1.96) represents a significant increase (α = 0.05) in VI for the period 2000–2014 and negative slopes (z ≤ −1.96) indicate a significant decrease (α = 0.05) over time.

4. Results

4.1. Evaluation of VI-Based SOS against in Situ GPP-Based SOS

SOS values of the VIs for the pixel location corresponding to FLUX towers were plotted against SOS estimates from in situ based GPP flux measurements to evaluate which of the three vegetation index SOS estimates that represent ground conditions most accurately (Figure 2). From the data that was available in this study, it can be seen that the relationship between satellite and ground observations is rather low. PPI SOS is better correlated with the GPP (r = 0.50; p < 0.01) phenology as compared to EVI (r = 0.11; p = 0. 16) and NDVI (r = 0.08; p = 0.24), indicating a more accurate capturing of SOS by PPI. NDVI shows SOS dates that start much earlier than SOS from GPP, PPI and EVI. Both EVI and PPI SOS dates are on average relatively close to GPP-based SOS (average ± std. dev. GPP SOS = 128 ± 21; PPI SOS = 126 ± 18; EVI 116 ± 25; and NDVI = 99 ± 21). Separating the data set into different land cover classes did not improve the relationships (Figure S1 and Table S1 in Supplementary Materials).

Seasonality of the VIs (estimated from the average of pixels used in VI phenology evaluation against GPP phenology) shows that NDVI seasonality differs significantly from the PPI seasonality by the earlier increase in spring and later decrease in autumn (Figure 2D). The average EVI curve follows the PPI curve closer than the NDVI curve with regards to spring increase and autumn decrease, but with less pronounced changes during the transition periods between growing and non growing seasons. Especially for NDVI, an abrupt shift in values can be observed (e.g., year 2000, 2005, 2007–2008, 2010–2011 and 2013), coinciding with changes in snow cover (checked by using the MODIS/Terra Snow Cover 8-Day product, MOD10A2).

4.2. Spatial Patterns of VI SOS Estimates

Averaged SOS of PPI (2000–2014) show a clear latitudinal gradient with a later onset towards higher northern latitudes (Figure 3A). PPI SOS on average starts around 5 May in the southern boreal forests. PPI SOS at the mid-latitudes shows an average SOS around 25 May, whereas the latest SOS is observed around 10–15 June at the high latitudes (East Siberian Taiga, the Eastern and Western Taiga Shield, Hudson Plains). Relative PPI SOS differences (ΔJulian days) to NDVI (NDVI_SOS–PPI_SOS) and EVI SOS (EVI_SOS–PPI_SOS) (Figure 3B,C) show that PPI generally yields later SOS dates as compared to EVI and NDVI, with EVI SOS dates being closer to PPI SOS as compared to NDVI SOS. Wetlands, croplands and mixed forest show the largest differences with NDVI and EVI SOS preceding PPI SOS by 20–42 and 8–34 days, respectively (

Table 1. Average ± one standard deviation Plant Phenology Index (PPI) start of season (SOS) dates for different land cover types (Figure 1) and the relative difference (in Julian days) in SOS to the Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI).

IGBP 1 Land Cover Class	PPI (Date)	NDVI (ΔDays)	EVI (ΔDays)
Evergreen Needleleaf forest	24 May ± 10	−16 ± 13	−7 ± 8
Deciduous Needleleaf forest	26 May ± 5	−15 ± 5	−13 ± 5
Deciduous Broadleaf forest	1 May ± 10	−16 ± 10	5 ± 9
Mixed forest	22 May ± 7	−20 ± 13	−8 ± 7
Open shrublands	2 June ± 8	9 ± 9	11 ± 9
Woody savannas	3 June ± 8	4 ± 10	5 ± 6
Savannas	4 June ± 8	−15 ± 10	−2 ± 7
Grasslands	26 May ± 15	10 ± 19	9 ± 16
Permanent wetlands	29 May ± 10	−42 ± 11	−34 ± 13
Croplands	27 May ± 11	−21 ± 8	−14 ± 6
Cropland/Natural vegetation	20 May ± 10	−20 ± 10	−12 ± 7
All	29 May ± 9	−13 ± 11	−6 ± 9

¹ International Geosphere-Biosphere Programme.

). For the two classes of needleleaf forest (evergreen and deciduous) covering vast areas of the study region (Figure 1B), a clear difference between PPI SOS and NDVI/EVI SOS is also observed (7–16 days difference). For areas of grassland and open shrublands, PPI SOS, however, precedes NDVI and EVI SOS with similar difference values to NDVI and EVI ranging from 5–11 days.

4.3. Trends in VI SOS Estimates

The results of VI SOS trends will be shown as two maps (TS trends and TS trends of significant changes) as relatively low amount of significant trends based on per-pixel trend analysis are observed caused by a relatively short period of analysis (2000–2014) and coarse resolution of the earth observation (EO) time-series [70].

4.3.1. PPI

The mean PPI SOS trend over the entire boreal region of the NH is showing an advance with an average of −0.28 (days/year) with a ±0.53 (days/year) standard deviation (

Table 2. Average ± one standard deviation trends in vegetation index (VI) start of season (SOS) dates (days/year) for different land cover types and boreal regions (Figure 1A,B).

IGBP Land Cover Classes	PPI	EVI	NDVI
Average for all land cover classes	−0.28 ± 0.53	−0.23 ± 0.72	−0.26 ± 0.65
Evergreen Needleleaf forest	−0.25 ± 0.58	−0.06 ± 0.92	0.04 ± 0.81
Deciduous Needleleaf forest	−0.47 ± 0.37	−0.62 ± 0.36	−0.74 ± 0.36
Deciduous Broadleaf forest	−0.54 ± 0.31	−0.34 ± 0.22	−0.18 ± 0.25
Mixed forest	−0.06 ± 0.45	0.02 ± 0.78	−0.18 ± 0.67
Open shrublands	−0.37 ± 0.51	−0.39 ± 0.57	−0.44 ± 0.50
Woody Savannas	−0.37 ± 0.55	−0.22 ± 0.71	−0.20 ± 0.57
Savannas	−0.56 ± 0.54	−0.37 ± 0.63	−0.41 ± 0.67
Grasslands	−0.12 ± 0.70	−0.12 ± 0.79	0.03 ± 0.84
Wetlands	−0.17 ± 0.51	−0.27 ± 0.70	−0.20 ± 0.61
Croplands	0.16 ± 0.63	0.24 ± 0.46	0.32 ± 0.46
Cropland/Natural vegetation	0.05 ± 0.49	0.14 ± 0.39	0.16 ± 0.50
Boreal regions
East Siberian Taiga	−0.35 ± 0.52	−0.43 ± 0.55	−0.56 ± 0.55
Iceland Boreal Birch Forests And Alpine Tundra	−0.23 ± 0.66	0.34 ± 0.97	0.50 ± 1.00
Kamchatka-Kurile Meadows And Sparse Forests	−0.27 ± 0.31	−0.46 ± 0.38	−0.47 ± 0.44
Kamchatka-Kurile Taiga	−0.33 ± 0.28	−0.32 ± 0.30	−0.30 ± 0.35
Northeast Siberian Taiga	−0.35 ± 0.42	−0.36 ± 0.49	−0.42 ± 0.41
Okhotsk-Manchurian Taiga	−0.11 ± 0.41	−0.14 ± 0.46	−0.15 ± 0.53
Sakhalin Island Taiga	−0.12 ± 0.32	−0.16 ± 0.47	−0.10 ± 0.41
Scandinavian and Russian Taiga (Russian part)	−0.30 ± 0.39	−0.04 ± 0.92	−0.10 ± 0.55
Trans-Baikal Conifer Forests	0.36 ± 0.43	0.18 ± 0.42	0.11 ± 0.70
Ural Montane Forests And Tundra	−0.22 ± 0.35	−0.11 ± 0.67	−0.25 ± 0.52
West Siberian Taiga	0.13 ± 0.48	0.03 ± 0.80	−0.15 ± 0.67
Alaska Range	−0.32 ± 0.55	−0.24 ± 0.72	−0.13 ± 0.68
Boreal Cordillera	−0.21 ± 0.54	−0.17 ± 0.83	−0.10 ± 0.67
Boreal Plains	−0.03 ± 0.50	0.02 ± 0.58	0.16 ± 0.73
Boreal Shield	−0.26 ± 0.63	0.02 ± 0.80	0.04 ± 0.87
Cook Inlet Basin	−0.09 ± 0.35	−0.03 ± 0.48	0.01 ± 0.51
Eastern Taiga Shield	−0.45 ± 0.53	−0.40 ± 0.49	−0.32 ± 0.52
Hudson Plains	−0.26 ± 0.50	0.06 ± 0.49	0.21 ± 0.47
Interior Alaska Taiga	−0.29 ± 0.51	−0.26 ± 0.70	−0.33 ± 0.61
Taiga Plains	−0.41 ± 0.43	−0.22 ± 0.70	−0.25 ± 0.68
Western Taiga Shield	−0.50 ± 0.47	−0.44 ± 0.77	−0.42 ± 0.60
Yukon Plateau And Flats	−0.32 ± 0.62	−0.33 ± 0.75	−0.43 ± 0.64
Scandinavian and Russian Taiga (Scandinavian part)	−0.53 ± 0.49	−0.25 ± 0.89	−0.27 ± 0.53

, Figure 4A,B). The strongest advance in SOS is observed for deciduous needleleaf and broadleaf forests (−0.47–−0.54) and savannas (−0.56), whereas croplands show a moderate delay in SOS (0.16) (Table 2). The advance in SOS is particularly pronounced in the Canadian boreal forests of the Eastern and Western Taiga shields (−0.45–−0.50) and Taiga plains (−0.45). In addition, the East Siberian Taiga and the Scandinavian and Russian Taiga show a significant advance (−0.35 and −0.30 days/year respectively). The advance in SOS over Alaska (Yukon Plateau, Alaska Range and the Interior Alaska Taiga) is on the same order of magnitude, with rates of −0.29–−0.32 (days/year). In contrast, the SOS trends over the West Siberian Taiga have shown a delay of 0.13 (days/year).

4.3.2. EVI and NDVI

The mean SOS trend for the NH boreal region is −0.23 (days/year) for EVI with a ±0.72 (days/year) standard deviation (Table 2, Figure 4C,D) and −0.26 (days/year) for NDVI with a ±0.65 (days/year) standard deviation (Table 2, Figure 4E,F) (Table 2). The result is close to PPI; however, trends in individual boreal regions or/and land cover classes show considerable diversity (Table 2) and overall the spatial variability in trends is smaller for PPI as compared to EVI and NDVI. NDVI showed an opposite trend as compared to PPI and EVI for evergreen needleleaf forests and grasslands and spatial variability in both EVI and NDVI is considerably higher than for PPI. EVI had a slightly positive trend for mixed forests, disagreeing with PPI and NDVI that showed an advance. The best agreement was observed in the East Siberian Taiga, where all three indices showed a strong significant advance in SOS and spatial variability of trends on the same order of magnitude. EVI was also in agreement with PPI in West Siberian Taiga, whereas NDVI showed an opposite direction in trends. Both EVI and NDVI showed delayed SOS over Iceland, whereas PPI showed an advance. The same patterns were observed over Boreal Plains in Canada, mainly caused by a strong delay of SOS over croplands that affected the mean trend.

5. Discussion

5.1. VI SOS Detection and Evaluation against GPP-SOS

PPI-based phenology retrieval for SOS assessment has shown an improved performance over NDVI and EVI with regards to several aspects. First of all, PPI was able to detect a more precise timing of SOS when evaluated against ground observed GPP data (Figure 2A–C). We found a surprisingly low and non-significant relation between EO based EVI/NDVI SOS and GPP based SOS (r-values of 0.11 and 0.08, respectively; 81 site-years of flux measurements used from 14 sites). In addition, a large difference in the timing of spring phenology between the NDVI-based SOS dates and in situ observations was reported in [39] for the period 2000–2011 in Western Central Europe. It is clear that a systematic difference between VIs is observed (Figure 2D), which is likely to be related to the different sensitivity to changes in snow cover and vegetation activity [36]. Our results comply with the tests of snow influence on VIs as implemented in [36,71]. The tests, as well as examination of VI relation to Normalized Difference Snow Index (NDSI), revealed a pronounced NDVI and EVI association to the changes in snow cover (NDVI was reported to be more sensitive to partly snow-covered surfaces, whereas EVI is more sensitive to high snow depth) [36,71]. The low r-values reported here indicate very limited to no sensitivity of EVI and NDVI to variability in ground observed SOS for areas characterized by winter snow. This result severely challenges the utility of NDVI/EVI for phenology assessment and change studies in high-latitude regions as such VI SOS estimates are likely to coincide with the disappearance of snow as also discussed in [11,35,42,50] unless careful data screening and preprocessing are done to avoid the impact from snow [28].

A possible source of uncertainty in the evaluation of SOS retrievals from the different VIs could be related to the threshold method as used here. Another commonly applied method (the derivative method) is to fit functions to the VI signal and to retrieve the SOS as the point of a specific growth rate (usually the maximum growth rate) of this function [72,73]. Wu et al. [72] tested different algorithms using both threshold and derivative methods on NDVI for estimating SOS against SOS derived from global FLUXNET sites and also reported moderate relationships regardless of the methods used. Nagai et al. [50] concluded that the derivative method is misleading for areas with snow cover on the forest floor as there is then two points in time with an increase in the phenology signal, one coinciding with the snowmelt and one with the leaf burst. In addition, in deciduous forests, there are two different leaf bursts during spring time with a corresponding increase in the phenology signal: the first for the ground vegetation and the second for the tree canopy. This may add interannual variability and thereby cause uncertainty in SOS estimates when using the derivative method.

EVI was found to perform moderately better than NDVI when compared against ground observed GPP SOS, similar to the results reported by [74]. PPI-SOS on the other hand showed a highly significant agreement with GPP-SOS, yet the correlation coefficient obtained was modest (r = 0.5). Possible sources of uncertainty impacting the relationship between EO VI SOS and ground observed SOS related to the footprint of satellite and ground measurements must be considered. The MODIS MCD43A4 500 m spatial resolution product was used (3 × 3 pixels) for the comparison with flux measurements for a better match with the spatial footprint of the towers whereas the MCD43C4 product of a 0.05° spatial resolution was used for the large scale analysis. However, using coarse resolution MODIS data when compared to flux tower measurements yielded higher r-values for PPI and NDVI (PPI = 0.56, EVI = 0.14, NDVI = 0.22). This shows that differences in the point to pixel scale mismatch between the two EO data sets used (500 m and 0.05° spatial resolution) and flux observations does not play a major role in explaining the low correlations obtained in the EO SOS evaluation. Uncertainty could also be related to that GPP SOS is a photosynthesis based estimate of phenology, whereas the vegetation indices are LAI-related [75]. For deciduous broad-leaved forests, it has been suggested that an earlier increase in the EVI signal in relation to GPP was caused by leafs burst earlier than the photosynthetic capacity [75]. For evergreen forests, needles or leaves are always present and it is instead temperature that drives the SOS of GPP [76]. Another source of uncertainty could be the GPP based estimates of SOS. Eddy covariance towers measures net ecosystem exchange of CO₂, a flux that is comprised of two different exchange processes: GPP and ecosystem respiration. At the beginning of the growing season, ecosystem respiration is the dominant flux process, and to disentangle the GPP signal can be difficult [77], and only when CO₂ uptake by the vegetation (GPP) starts to be the dominant flux process can the partitioning between the two processes be less problematic. The point in time when partitioning between GPP and ecosystem respiration can differ between years and sites, and thereby add uncertainty in the SOS estimates. A final plausible explanation could be that, as opposed to photosynthesis (GPP measured from flux towers), greenness indices are not reduced down to zero in evergreen needleleaf forest in dormant periods. Therefore, NDVI/EVI for such biomes do not provide an unambiguous estimate on when actual photosynthesis starts or ends as discussed in [43], in which used sun-induced chlorophyll fluorescence (SIF) from GOME-2 was used as a direct proxy of photosynthetic activity to monitor seasonality.

5.2. VI SOS Spatial Patterns

Average PPI-derived SOS occurs in the beginning and middle of May at southern latitudes but is delayed towards the beginning of June with an increase in latitudes (Figure 3a). EVI and especially NDVI have showed up to 20–40 days earlier SOS than PPI-derived SOS (Figure 3B,C). These differences in timing are likely to be caused by EVI and NDVI responding to variations in the snow season, rather than vegetation phenology. The greatest differences have been observed in the southern latitudes of the study area. The same pattern was observed in [48] where PPI phenology was compared to EVI phenology from MCD12Q2 land surface phenology product and showed that EVI detected SOS as much as one month earlier than PPI. This pattern was also explained by a larger temporal lag between the canopy growth onset and the snow melting time in the south than in the north [48]. The high latitude boreal zones are located within regions of subarctic and cold continental climate, which is defined by long and severe winters and short summers, limited by 50–100 frost-free days [78]. The short summer periods result in strong vegetation sensitivity to air temperature, leading to an almost immediate start of the growing season after snow melt, in order to ripen seeds successfully [78,79,80]. Hence, an overlap of the timing of snowmelt and the onset of canopy growth at the high latitudes is likely to cause the convergence of PPI, NDVI and EVI SOS estimates (Figure 3).

5.3. VI SOS Temporal Patterns

Temporal trends in observed SOS dates agree for the different VIs at the overall level (Figure 4) despite the significantly better performance of PPI in SOS detection as compared to EVI and PPI. However, the analysis for different land cover types and boreal regions reveals distinct differences between trends of VIs. Evergreen needleleaf forest is characterized by a markedly negative trend in PPI (0.25 days/year), suggesting a change towards earlier SOS dates, whereas NDVI shows the opposite (weak) trend towards later SOS dates (0.04 days/year) (Table 2). Some of the major boreal regions do also differ; e.g., the Boreal Shield, Boreal Plains and Hudson Plains do all show negative PPI trends towards earlier SOS dates (−0.26, −0.03 and −0.26 days/year), whereas NDVI shows positive trends (0.04, 0.16 and 0.21 days/year). These patterns possibly indicate a later snowmelt, but an earlier start of the growing season. For the West Siberian Taiga, the pattern is reversed, with PPI showing a trend towards later SOS (0.13 days/year) and NDVI shows an advance in SOS (−0.15 days/year). The delay over the West Siberian Taiga is likely to be associated with its geographical setting being the world’s largest continuous lowland, with more than a half of the area <100 m above sea level and densely covered by wetlands and floodplains [81]. The delay in SOS can therefore be caused by high sensitivity to altered hydrological conditions of the wetlands as was also reported for the East Siberian Taiga [82]. The delay in SOS observed over the Boreal Plains in Canada for areas dominated by croplands is likely related to changes in cropping practices. Trends for areas under substantial anthropogenic influence are generally expected to be decoupled from the patterns of trends of natural forested ecosystems dominated by changes in climate.

The observed patterns of trends (Figure 4) showed good agreement with existing studies [7,35]; however, we find fewer areas of significant trends in the current study. This is likely due to a short time-period of investigation (as compared to Global Inventory Modeling and mapping Studies (GIMMS) based studies starting in 1982), reduced statistical significance for phenology trends in recent decades [24], outliers in the noisy VI time-series of the high latitudes, and per-pixel based trend estimation method [48,70,83]. A better estimation of significant trends could have been achieved by applying panel analysis [84] instead of pixel-based trend analysis, which is reported to give non or minor significant trends in the high latitudes [85,86]. To avoid areas of abrupt phenology changes caused by fires or deforestation, a dataset of forest loss during 2000–2014 [55] was used to exclude such areas. However, it takes several decades for regrowth of a forest, and it is likely that pixels affected by forest fire or logging activities prior to 2000 are still present in the analysis. Forest fires have a negative effect on satellite-based phenology retrieval, since spectral properties of burnt scar-vegetation signal persist up to several years and may cause observations diverging from the overall trend of the landscape [85].

The overall lower spatial variability in trends for PPI as compared to EVI and NDVI (Table 2) suggest that PPI is less prone to many of the abovementioned confounding factors when determining SOS as compared to the conventional indices. The standard deviations of the trends reported (Table 2), however, also indicate that trends in SOS (days/year) for many International Geosphere-Biosphere Programme (IGBP) land cover classes and boreal regions are highly varying, and, in many cases, the standard deviations are higher than the trends reported. High latitude ecosystems are highly responsive to the interannual variation in the drivers of phenology, and even if the results presented here overall support the earlier onset of growth in the Northern Hemisphere boreal zone, the spatial variability reveals considerable differences for the period of analysis. This accentuates the need for developing improved EO-based metrics for accurate assessments of SOS before drawing generalized conclusions about changes in SOS at high latitudes.

6. Conclusions

This study analyzes the robustness of the satellite derived Plant Phenology Index (PPI) for estimating start of season (SOS) in boreal forests of the Northern Hemisphere (2000–2014), in relation to two other widely applied indices NDVI and EVI. These satellite-based observations were evaluated against GPP-retrieved SOS derived from a network of Flux tower observations. Results revealed that PPI performs significantly better in boreal regions as compared to NDVI and EVI, where very poor relations to ground observations were found. The main challenge related to SOS retrieval by using NDVI or EVI is their sensitivity to changes in snow-cover rather than vegetation canopy changes, and lower response to variations in vegetation growth dynamics. It was shown that NDVI and EVI based SOS occurred up to 20–40 days earlier than PPI in the southern part of the boreal region driven by an increased temporal decoupling between the date of snow disappearance and the date of onset of photosynthetic activity. At higher latitudes, the convergence of snow disappearance and onset of photosynthetic activity results in a larger similarity between PPI and NDVI/EVI SOS.

The SOS trends observed in this study have shown an overall advance in PPI SOS at the rate of 0.28 (days/year) as compared to −0.23 for EVI and −0.26 for NDVI, but with considerably larger differences between PPI and EVI/NDVI for specific land cover types and boreal regions. The most pronounced patterns of significant trends of advancing SOS when considering all three indices were observed for the land cover class of deciduous needleleaf forest, where PPI, however, revealed less dramatic changes as compared to EVI/NDVI (0.47 days/year as compared to 0.62 and 0.74). On the other hand, advancing PPI SOS trends were found to be even higher for deciduous broadleaf forests and savannas (0.54 and 0.56 days/year) not equally reflected in EVI (0.34/0.37) and NDVI (0.18/0.41). The largest patterns of significant PPI trends towards earlier SOS were observed in the Eastern and Western Taiga Shield in Canada (0.45 and 0.50 days/year) followed by the East and Northeast Siberian Taiga in Russia (both 0.35 days/year).

Taken together, the findings of this study support the suggested improved performance of PPI over NDVI and EVI in SOS retrieval for areas of seasonal snow cover. The considerable difference in the spatial patterns of SOS trends between PPI and NDVI/EVI over the 15-year study period suggests that precautions must be taken when using and interpreting NDVI/EVI derived SOS results from northern latitudes and improved accuracy in assessing spatio-temporal patterns of SOS can be gained by applying PPI.