Heterogeneous spring phenology shifts affected by climate: supportive evidence from two remotely sensed vegetation indices

The Northern Hemisphere spring greenup (SG) has advanced between 0–12 days per decade since early 1980s as inferred from multiple satellite time series. The wide range of SG shifts is mainly due to the fact that these studies cover different periods and regions, and using different satellite records. Assessing the spatial heterogeneity of SG trends associated with different satellites is essential for robustly interpreting phenological dynamics and their responses to climate. We investigated the heterogeneity of the SG trends and their responses to climate variability with two satellite products (1) Terra Moderate Resolution Imaging Spectroradiometer (MODIS) and (2) Advanced Very High Resolution Radiometer (AVHRR) over the period 2001–2013. Both MODIS and AVHRR agreed in showing the spatial distribution of mean SG, and SG advancement in northern Canada, the eastern United States, and Russia, and SG delay in western North America, parts of Baltic Europe, and East Asia. However, we identified contrasting MODIS and AVHRR SG trends in the northern high latitudes. Our analyses of correlations between SG and preseason climate drivers indicated that temperature dominated the interannual variability of SG. Preseason, the period preceding SG and highly correlated with the timing of SG has experienced much stronger warming than the spring season. MODIS and AVHRR indicated consistent temperature sensitivity of SG across biomes, even though the MODIS inferred SG is better correlated and more sensitive to temperature across biomes as compared to AVHRR. The sensitivities of SG to temperature across biomes is stable but with a slight increase over 2001–2013, in comparison with that over 1988–2000. The increased SG-temperature sensitivity is associated with increased precipitation during the spring season, which regulated the sensitivity of SG to spring temperature.


Introduction
Vegetation phenology plays an important role in regulating energy, water, and trace-gas exchanges between land and atmosphere (Galvagno et al 2013, Richardson et al 2013. As the satellite-based vegetation indices are extending to longer time spans, many studies have applied vegetation indices to investigate vegetation spring greenup across regional and global scales (table 1). The Normalized Difference Vegetation Index (NDVI) computed from radiances collected by the Advanced Very High Resolution Radiometer (AVHRR) has been applied in numerous vegetation phenology studies for their long time span. Because of the improved sensor calibration and atmospheric corrections, and higher spatial resolution, the latest NDVI products retrieved from Terra Moderate Resolution Imaging Spectroradiometer (MODIS) and Système Pour l'Observation de la Terre (SPOT) VEGETATION mission (1 km)(e.g., Durpaire et al 1995) are considered to be improved over AVHRR (Zhang et al 2003).
The SG has been advanced in the range of 0 to 12 days per decade in the Northern Hemisphere according to NDVI observations (table 1). The wide range of SG shift is from studies covering different periods and regions, and using different satellite records. According to AVHRR NDVI, mean SG over 45-75°N advanced 6.2 days per decade (days decade −1 ) over 1982-1991, while during 1992-1999, the advancement rate decreased to 2.4 days decade −1 . The Northern Hemisphere SG advanced about 3 days decade −1 during 1982-1999, however, the advancement was negligible during 2000-2008 (Jeong et al 2011). The magnitude of SG advancement differed across regions, e.g., North America was usually reported to have a higher advancement than Eurasia during 1981-1999(Zhou et al 2001), 1985-2000(de Beurs and Henebry 2005 based on AVHRR NDVI, and 2000-2010 based on MODIS NDVI (Zeng et al 2011).
As some studies covered similar periods and regional extents, the wide range of SG shift is likely to be contributed by different NDVI datasets and methods that have been used to derive phenology metrics. Different methods, when applied to the same NDVI products over the same period, can lead to consistent sign of SG trends across regions and vegetation types (Cong et al 2013). The differences in SG shift derived by different datasets are much larger than the differences led by different methods (Shen et al 2014b). The advancement of SG is about 10.4 days decade −1 from 2001 to 2012 in Tibetan Plateau, based on combined AVHRR NDVI over 1982-2000and SPOT-VGT NDVI over 2001(Zhang et al 2013. However, this magnitude of SG shift is very different from the SG trend over 2000-2011 inferred from single AVHRR NDVI, which inferred insignificant SG trend in the Tibetan Plateau (Ding et al 2016).
Discrepancies in SG inferred from different NDVI products were also found in northern mid-to highlatitudes. MODIS and SPOT-VGT agreed in inferring SG magnitudes and trends during their overlapping period in western Arctic Russia where, however, AVHRR SG showed a very different SG after 2000s in both magnitude and sign . The magnitude of SG trend over 2002-2012 inferred by MODIS NDVI (−5.9 days decade −1 ) is almost three times higher than that over the period 1982-2002 inferred using the AVHRR NDVI (−1.9 days decade −1 ) in the Northern Hemisphere . However, another study revealed an negligible SG advancing trend in Northern Hemisphere over 2000-2008 (−0.3 days decade −1 ) in relative to 1980-1999 using AVHRR NDVI time series (−2.9 days decade −1 ) (Jeong et al 2011).
In this study, we (1) analyze the heterogeneous SG trend distribution across regions and biomes and quantify how much of the heterogeneity is associated with NDVI products and (2) explore climate controls on the spatial heterogeneity of SG variation. SGs and their respective sensitivities to climate drivers over the period 2001-2013 are compared between AVHRR and MODIS. The preseason, a period preceding SG during which the climate factors determine the timing of SG and spring (March, April and May) climate and how these affect the temporal and spatial heterogeneities of SG were analyzed by using an independent climate reanalysis dataset.

Study area and biomes
We limited the analyses to north of 30°N because vegetation phenology in this region is mainly controlled by the annual cycle of temperature (Linderholm 2006, Fu et al 2014, Güsewell et al 2017, and may be regulated by water availability (Peñuelas et al 2004) and photoperiod (Singh et al 2017). In order to analyze SG and its climate influences across biomes, we masked the satellite-based SG results with MODIS IGBP classification of land cover types (MCD12Q1) with spatial resolution of 0.5°×0.5° (Channan et al 2014). The SG responses to climate drivers were analyzed for 9 dominant biomes in our studied region. The description of the 9 biomes is given in supplementary S1 is available online at stacks.iop.org/ERC/1/091004/mmedia.

Vegetation indices and SG determination
We used AVHRR NDVI time series (GIMMS NDVI3g) extended from July 1981 to December 2013 with spatial and temporal resolution of 1/12°and bimonthly, respectively (Pinzon and Tucker 2014). The 16-day MODIS NDVI composites (MOD13C1, collection 6) with 0.05°spatial resolution and 16-day temporal resolution were used.
Both AVHRR and MODIS NDVI data were resampled to 0.5°×0.5°with the mean value of the original high resolution NDVI within a 0.5°×0.5°pixel in order to match the spatial resolution of climate reanalysis. We screened the pixels with annual maximum NDVI <0.1 to exclude the non-vegetated pixels. For GIMMS NDVI3g, an average seasonal profile or spline interpolation was applied to fill the pixels identified with snow or ice cover (Pinzon and Tucker 2014), which were applied to MOD13C1 pixels flagged with snow/ice. In the dormant season, the filled values are approximate to zero and usually smoothed by the piecewise logistic method described below.
The regression methods that have been used to reconstruct NDVI time series and derive SG include the Savitzky-Golay fitting method, spline smoothing, asymmetric Gaussian functions, logistic functions, and harmonic analysis of times series. In order to exclude the uncertainty from different reconstruction methods and NDVI products, we used piecewise logistic method which uses least-square fitting applied to the first half of the growing season (Zhang et al 2003). We assessed the piecewise logistic method by comparing its inferred SG with that inferred by Savitzky-Golay and asymmetric Gaussian in different biomes (supplementary S2).
We first applied piecewise logistic method to the vegetation growth stage. We then defined the date of SG (D SG ) as the Julian date with the largest change in vegetation growth which is defined as the local maximum curvature of the fitting curve. We calculated the average of D SG (D SG ) in each pixel from 2001 to 2013. For pixels with multiple growth cycles in a year, this piecewise logistic method is applied to the first growth cycle.
We used CRU-NCEP v6 dataset to calculate the preseason period and climate separately for temperature and precipitation (Shen et al 2014a). In each pixel, we defined the preseason climate (T m and P c ) over the period preceding D SG from 15 to 120 days with an increment of 3 days and calculated T m and P c during the respective preseason periods. T m and P c were then detrended over 2001-2013. We expect that the relative change in precipitation is more relevant than the absolute changes in determining phenology, so we used the relative change in cumulative precipitation to represent the percentage change in precipitation (%) rather than the absolute cumulative precipitation change in millimeter (mm). For each period preceding D SG for a given pixel, we calculated the Pearson's correlation coefficients (PCC) between D SG and T m (and P c ). We further analyzed partial correlation between D SG and T m (and P c ) to consider the covariation of both temperature and precipitation.
We screened the data by removing the pixels with a positive interannual correlation between preseason temperature and D SG and preseason precipitation and D , SG respectively. The period with the greatest negative correlation between D SG and T m (and P c ) is defined as preseason P T (and P P ). Superscripts A and M represents variable derived from AVHRR and MODIS, respectively (e.g. D SG M and D SG A are D SG derived from MODIS and AVHRR, respectively. As the preseasons are inferred from the negative correlation between D SG and climate preceding D , SG the difference in AVHRR and MODIS inferred D SG could lead to different preseason length and climate. We used the spring temperature (T s ) and precipitation (P s ) to calculate the sensitivities of D SG to climate over 2001-2013 by applying linear regressions between D SG and T s (and P s ).

Heterogeneous SG trend
We inferred consistent spatial patterns of mean D SG (r=0.83, p<0.01) between MODIS and AVHRR. In the north of 30°N, about half of the pixels had an inter-annual correlation above 0.5 (p<0.1) ( figure 1(a)). The SG shift has shown strong spatial heterogeneity, according to both MODIS and AVHRR NDVI. Northern Canada, Eastern United States, and Russia appeared relatively consistent SG advancement, opposing D SG delay in western North America, parts of Baltic Europe and East Asia (figures 1(b) and (c)). In Eurasia, where D SG advanced, D SG M advancement was greater than D . SG The sign of D SG anomalies was consistent between MODIS and AVHRR over 30-60°N (figures 2(a), (c) and (e)). However, north of 60°N, D SG anomalies between MODIS and AVHRR were substantially different (figure 2(b)), which is mainly due to the opposite sign of D SG anomalies between MODIS and AVHRR in North America (figure 2(d)). In Eurasia, both MODIS and AVHRR indicated D SG anomalies with the same sign. However, the magnitude of D SG A anomalies was much smaller than that of D SG M (figure 2(f)). A large transition in the D SG A anomalies occurred around 2000 in North America where D SG A over 2001-2013 is 5-6 days later than that over 1982-2000.

Preseason climate regulating SG
Both partial correlation and Pearson correlation between D SG and preseason climate were calculated. The spatial patterns of partial correlation and Pearson correlation are similar, while the partial correlations are slightly lower than Pearson correlations for both temperature and precipitation. Here we used the Pearson correlation to analyze relationship between D SG and preseason climate. The correlation between D SG and preseason climate indicated a different coupling between D SG and preseason climate for MODIS and AVHRR (figure 3). The mean correlation between D and SG M T m is 0.57±0.2 during its corresponding preseason ( ) P T M in the Northern Hemisphere, which is much higher than the  (b)). The preseason length for AVHRR (62±38 days) and MODIS (41±31 days) that we inferred from the correlation between T m and D SG differed due to the differences between D SG G and D .
SG M The fraction of the northern mid-to high-latitude land surface correlated with preseason precipitation is less than that correlated with temperature for both AVHRR and MODIS. MODIS and AVHRR showed consistent spatial patterns in correlations between P c and D SG and between P c and preseason length (figures 3(c) and (d) The spatial pattern of preseason temperature is consistent with the spring temperature trend. However, the preseason temperature change is more intensified than spring temperature change (figure 4(c)).
The preseason precipitation trend is negligible and the spatial distribution of the trend is more uneven as compared to the temperature trend (figures 4(d) and (e)). The spatial pattern of P P M and P P A precipitation trends are less consistent (r=0.40, p<0.01), as compared to that of temperature trends. The preseason wetted in the mid to eastern United States, Western Canada, Northern Norway, and Northwestern Russia. The maximum wetting trend is about 7 mm yr −1 . The preseason only dried in southeastern United States and scattered in Eurasia. The largest preseason drying trend is about 4 mm yr −1 . The spatial pattern of preseason precipitation trends is consistent with the spring precipitation trend. However, the preseason precipitation change is milder than spring precipitation change (figure 4(f)).

SG shift in response to spring climate
We evaluated the SG and its response to climate based on D SG sensitivity to spring temperature and precipitation, to avoid the impact of varied preseason length and climate correlated with MODIS and AVHRR inferred SG. MODIS and AVHRR showed varied D SG sensitivity to spring temperature and precipitation. Both MODIS and AVHRR records had dominant negative correlation between D SG and spring temperature, i.e., the warmer spring, the earlier D SG (figures 5(a) and (b)). D SG M is more sensitive to spring temperature than D .

SG
A The mean sensitivity of D SG M to temperature is about -3.59 days°C −1 (per°C warming in spring), while mean sensitivity of D SG A to temperature is about −2.59 days°C −1 . The D SG is not significantly relied on the spring precipitation, as inferred by the correlation between D SG and preseason precipitation (figures 5(c) and (d)).
Due to the lower correlation between D SG and precipitation, and less extent that shows D SG sensitivities to precipitation, we only analyzed biome-scale sensitivity of D SG to T s (figure 6). The inter-biome variation of D SG to spring temperature is consistent between MODIS and AVHRR (r=0.89, p<0.001), although D SG M is more sensitive to T s than that of D SG A in all the biomes. The sensitivities of D SG to T s is less different between MODIS and AVHRR in the mid and low latitude biomes. The differences in D SG to T s sensitivity are especially significant in northern biomes, e.g., sensitivity of D SG M to T s in open shrublands and northern grasslands are more than 50% higher than sensitivity of D SG A to T s in these biomes. As the AVHRR NDVI product extended to the early 1980s, we performed a comparison of D SG A to T s sensitivity between the periods 1988-2000 and 2001-2013, with the same length of time (13 years). The interbiome sensitivity of D SG A to T s is stable (r=0.98, p<0.0001) during these two periods ( figure 6). The sensitivity of D SG A to T s was slightly increased for all the biomes over 2001-2013, as compared to that over 1988-2000. The sensitivity increase was the highest in the deciduous needleleaf forest (42.3%) and permanent wetland (31.5%), but relatively stable in deciduous broadleaf forests and evergreen needle leaf forest (<5%).

Heterogeneous SG trend
We analyzed spring greenup and how it was affected by preseason climate over the period 2001-2013 by using MODIS and AVHRR NDVI products. Whether it is MODIS or AVHRR-inferred SG, SG showed a strong spatial heterogeneity. Both MODIS and AVHRR showed large areas of delayed SG in North America and advanced SG in Eurasia. The main differences between the SG trend inferred by MODIS and AVHRR were in high latitude regions, i.e. Alaska, Russia Far East, and circumpolar Siberia. The delayed AVHRR SG in North America and Eurasia contrasted the strong advanced SG in these region over 1982-1999 with AVHRR NDVI (Jeong et al 2011), which

SG sensitivities to climate
The MODIS and AVHRR SG to preseason climate correlation showed that temperature is the dominant factor regulating spring phenology. The warmer the preseason, the earlier the spring greenup. The correlation between MODIS SG and preseason temperature is higher, indicating that the coupling between MODIS SG and climate variability is stronger. Kern et al (2016) analyzed the correlation between temperature and August NDVI Figure 4. The preseason temperature trend (°C yr −1 ) calculated from CRUNCEP correlated to spring greenup date inferred from AVHRR (a) and MODIS (b) NDVI and precipitation trend (mm yr -1 ) calculated from CRUNCEP correlated to spring greenup date inferred from AVHRR (c) and MODIS (d) NDVI. The shaded regions indicate that the trend is significant (p<0.1). Figure 5. Sensitivity of spring greenup to spring temperature (days°C −1 ) for AVHRR (a) and MODIS (b) and sensitivity of spring greenup to spring precipitation (days % −1 of precipitation increases) for AVHRR (c) and MODIS (d).  60°N (GLN), from temperate grassland in the south (GLS) due to their expected differences in climate and controls on phenology. The error bars represent standard deviation of the SG sensitivities to spring temperature for each biome. The number of pixels for calculating biome-scale SG sensitivities to spring temperature are provided in supplementary table S1. anomalies and found that the correlation between MODIS NDVI and temperature was higher than that between AVHRR NDVI and temperature in central Europe. The correlation between SG and temperature is stronger than precipitation, which is consistent with our previous SG to climate sensitivity over 1982-2005 by AVHRR NDVI and models (Xu et al 2018) and as reported by Shen et al (2015).
The sensitivity of SG to spring temperature across biomes is consistent between MODIS and AVHRR. Both inferred strong SG-temperature sensitivities in the forest biomes. The northern permanent wetlands and evergreen needleleaf forests are relatively more sensitive to spring temperature changes compared to the northern grasslands and open shrubland in circumpolar regions. The across-biome sensitivity is compatible with the finding with AVHRR NDVI  that vegetation with earlier growing season showed greater sensitivities to temperature (Shen et al 2014a). In contrast, within tundra biomes, phenology is more sensitive to temperature in colder, higher latitudes than in warmer regions, according to the site-level observation (Prevéy et al 2017).
The varied sensitivity of SG to temperature for boreal biomes is usually linked to multiple factors. Under warmer climate conditions, the temperature sensitivity of spring phenology tends to decline because (1) In spite of the continued warming in spring, winter chilling deficiency can delay SG (Yu et al 2010), (2) shorter photoperiods limit leaf development with earlier SG (Chmielewski andGötz 2016, Garonna et al 2018), (3) greenup may respond to temperature nonlinearly and be saturated at high temperature (Caffarra et al 2011), and (4) under warmer conditions, the thermal forcing is shifted or reduced, reducing SG sensitivity to temperature (Friedl et al 2014, Güsewell et al 2017. Water supply is an important regulator of the sensitivity of vegetation phenology and growth to temperature . Particularly in arid and semi-arid biomes, vegetation phenology is more sensitive to interannual variations in preseason precipitation in more arid areas and more sensitive to preseason temperature in wetter than drier areas (Shen et al 2015b). The water stress can limit phenology response to the preseason daily maximum temperature (Shen et al 2016). The spring droughts can ever delay the spring phenology in grassland (Kang et al 2018).
In our results, wetter springs during 2001-2013 than during 1988-2000 can explain the enhanced SG response to preseason temperature. Across all biomes, spring precipitation increased by 11 mm during 2001-2013 compared to 1988-2000. The largest enhanced sensitivity occurred in deciduous needleleaf forests, with spring precipitation increasing from 63 mm over 1988-2000 to 175 mm over 2001-2013. Despite greater SG sensitivities to spring temperature over 2001-2013, our results did not imply a greater SG advance than that over 1988-2000 as proposed by Shen et al (2015a), which is probably due to data quality issue raised by (Zhang et al 2013) and the large transition before and after 2000 found in AVHRR SG.

Differences in SG as derived by MODIS and AVHRR NDVI
The different MODIS-and AVHRR-inferred SG trends are mainly in high latitudes. The lack of ground calibration, the presence of evergreen vegetation types, and snow make the NDVI retrieval and phenological determination complex and uncertain in high latitudes. The amount of evergreen vegetation foliage in high latitudes remains relatively stable through seasons, so that the seasonal variation features of leaf activity can be masked by noise inherent to satellite-based radiance measurements (Gamon et al 2016). The snow cover affects the greenness determination on pixels in northern regions (Moulin et al 1997) and led to gaps of NDVI during the dormant season. Therefore, the NDVI time series cannot be adequately filled during the snow melting and vegetation greening transition season (Zhou et al 2015).
The overlapping periods of snowmelt and greenup lead to complexity in greenup determination. In the Yamal Peninsula, the spring greenup date was almost the same as the snow-end date , so that snow cover affects the accuracy of vegetation greenup time identification. In high latitudes, the start of the growing season is often depends on snowmelt rather than temperature (Semenchuk et al 2016). Although snow cover can affect the determination of SG, the causes of the difference in MODIS and AVHRR inferred SG might stem from other factors. In some locations of Canada and Sweden, MODIS NDVI was inconsistent with AVHRR NDVI, even if the pixels affected by snow cover were removed. For example, MODIS NDVI was lower than AVHRR NDVI in the dormant season (Fensholt and Proud 2012), which could lead to a slower MODIS-inferred transition from dormant to growing season.
The effects of interactions between vegetation and snow on NDVI retrievals require further investigation. MODIS NDVI inferred that the average SG advanced 0.96 days year −1 over 52-75°N between 2001 and 2013 in our results. The magnitude of SG advancement is greater than the MODIS snow-end date advancement of 0.37 days year −1 in this region between 2001 and 2014 (Chen et al 2015). The lag of snow phenology advancement implies that snow complexity in determining cold-region SG is still present in a warmer climate. To reduce the snow effect on SG determination, the normalized difference water index method (Delbart et al 2005(Delbart et al , 2006, plant phenology index method (Jin et al 2017), normalized difference vegetation index-normalized difference infrared index phase-space method (Thompson et al 2015) have been proposed to improve the NDVI-based phenological metrics.
The different SG trends after 2000 reported in other studies and our findings are more likely to be related to differences in NDVI products. For example, MODIS and AVHRR sensor channels have different spectrum ranges. MODIS and AVHRR were designed with different spectral ranges for both the red and near-infrared bands, which can lead to differences in NDVI. The retrieved NDVI is also a mixture of different vegetation species with diverse phenologies, bare soil, and even water bodies dependent on the spatial resolution (Helman 2018). Both AVHRR and MODIS NDVI were generated using daily surface reflectance products to a similar composite interval. However, the MOD13C1 applied the constrained-view angle-maximum value composite while NDVI3g applied maximum value composite. The maximum value composite cannot completely remove atmospheric effects (Pinzon and Tucker 2014) and the different composite technique can cause differences in the same interval (Gallo et al 2004). The large transition of GIMMS SG anomalies around 2000 may be associated with the sensor transition from AVHRR/2 to AVHRR/3 (Pinzon and Tucker 2014). These sensor transition caused data issues potentially influence the interpretation of the SG trend and its sensitivity to climate drivers.

Conclusions
We analyzed the heterogeneity of spring green up timing and its response to preseason and spring climate over 2001-2013 using MODIS and AVHRR NDVI. Both NDVI products have heterogeneous SG trends and sensitivities to spring temperature. The differences in MODIS-and AVHRR-inferred spring greenup trends are mainly in northern high latitudes (>60°N). The differences are primarily in opposite anomalies in the time of spring greenup in North America. Temperature is the dominant climate driver of for MODIS and AVHRR spring greenup and MODIS inferred SG has a better correlation with and more sensitive to temperature. The sensitivity of greenup to temperature has increased slightly over time for all biomes related to wetting spring. Vegetation index products are useful to quantify plant dynamics. However, substantial uncertainties still exist in these products, especially over high latitudes. We therefore recommend continued assessment of data reliability through multi-sensor comparisons and ground-based observations in high latitudes.