Introduction

Consistent estimates of plant phenology over continuous spatial and temporal domains are essential for quantifying climate change impacts on ecosystems (Walther et al. 2002; IPCC 2014) and for understanding the mechanistic basis of phenology (Pau et al. 2011) and interactions among involved biotic and abiotic factors (Richardson et al. 2013). Many studies based on direct human observations of vegetation have indicated strong phenological responses to climate during recent decades (Fu et al. 2014a; Menzel et al. 2006; Chmielewski and Rötzer 2001). For example, across Europe, the advancement of spring phenology matched the warming pattern during 1971–2000 (Menzel et al. 2006). Over western Europe, continued advancement of the start of the growing season (SOS) has been observed, including the two recent decades (Fu et al. 2014a).

Satellite observations are invaluable for complementing these point scale studies and investigating broad-scale phenology variations, allowing for regional and global studies of climate impact and sensitivities. However, several recent studies of spring phenology based on satellite data have produced inconsistent results. For example, Fu et al. (2014a), using the normalized difference vegetation index (NDVI) from the Advanced Very High-Resolution Radiometer (AVHRR) Global Inventory Modeling and Mapping Studies (GIMMS) dataset, showed a delayed SOS during 2000–2011 in western Central Europe, indicating a reversal of the 1982–1999 advancement. Furthermore, Wang et al. (2015b), who examined SOS over the northern hemisphere using NDVI datasets from different satellite sensors, found inconsistent SOS trends depending on satellite platforms. Other studies based on GIMMS NDVI suggested that the northern hemisphere SOS advancement had weakened or even reversed in the 2000s compared with its trend in the 1980s and 1990s (Jeong et al. 2011; Barichivich et al. 2013). These inconsistencies indicate considerable uncertainties in the remotely sensed phenology estimates, caused by, e.g., large inter-annual variations, short time periods, missing data, and inability of commonly used vegetation indices to correctly handle snow conditions and dense northern forest canopies (Huete et al. 2002; Delbart et al. 2005; Jönsson et al. 2010).

Air temperature is the main climatic factor regulating the onset of plant growth at boreal and temperate forests, including direct regulation by spring warming accumulation for bud burst and indirect regulation by winter chilling accumulation for bud rest break (Hänninen 2016). However, many studies have reported inconsistent spring phenology temperature sensitivities using different data sources, e.g., warming experiments (Wolkovich et al. 2012), ground phenological observations (Wolkovich et al. 2012; Menzel et al. 2006; Chmielewski and Rötzer 2001), and satellite observations of land surface greenness (Piao et al. 2006). Precipitation during winter and spring has also been found to affect the spring phenology in a complex manner at northern middle and high latitudes (Fu et al. 2015a, b; Piao et al. 2006; Cong et al. 2013; Yun et al. 2018). Environmental conditions during the previous year have a legacy effect on spring events of the next year: the warming-promoted summer growth and primary production will carry over into other seasons (Xia et al. 2014), whereas summer drought may negatively affect leafing in the next spring (Chmielewski et al. 2005). It was found that climatic factors in a complete annual cycle before spring had different importance and seasonal carry-over effects on leaf unfolding in Mediterranean ecosystems (Gordo and Sanz 2010) and on spring flowering in central England (Fitter et al. 1995).

In order to improve on recent satellite-based spring phenology estimates, this study uses a newly developed physically based vegetation index, the plant phenology index (PPI, Jin and Eklundh 2014). The PPI delineates vegetation canopy growth dynamics by addressing two major weaknesses of the NDVI: the sensitivity to snow and background variations (Huete et al. 2002; Jönsson et al. 2010), and the lack of sensitivity to leaf biomass variations in dense vegetation canopies (Fassnacht et al. 1994; Qi et al. 1994). PPI has been shown capable of disentangling remotely sensed plant phenology from snow seasonality (Jin et al. 2017), thereby proving its usefulness for spring phenology retrieval over the northern high latitudes (Karkauskaite et al. 2017).

To tackle the common problem of weak estimates from single-pixel analyses (Keenan et al. 2014), we use panel data analysis (Hsiao 2003) to estimate the phenology trend and phenology-climate sensitivities and correlations. We summarize our results over the study area using a meta-analysis on regression slopes (Becker and Wu 2007) instead of using a simple average method or discarding insignificant values (Garonna et al. 2014). The full annual cycle before the mean SOS dates of a pixel is divided into spring, winter, and summer, and the seasonal climate variables are generated and linked to SOS. With the newly proposed vegetation index PPI and advanced statistical analyses, the study aims to investigate recent trends in spring phenology and its sensitivities to temperature and precipitation in periods before the SOS, and to investigate the relative importance of seasonal climates on spring phenology over northern Europe (> 50° N) for the period 2000–2016.

Materials and methods

Study area

This study focuses on the northern European latitudes (> 50° N), covering a range of climate zones from subarctic in the north to temperate in the south, and a temperature gradient from − 2 to + 11 °C of mean annual temperature. The main land covers are evergreen needleleaf forest (ENF), deciduous forest (DBF), wetland (WET), and cropland and grassland (CRO) (Fig. 1).

Fig. 1
figure 1

Map of the study area shown as the major land cover types from CORINE Land Cover 2006 (EEA 2007) and the MODIS land cover type product MCD12Q1 (Friedl et al. 2002). WET, wetland; ENF, evergreen needleleaf forest; DBF, deciduous broadleaf forest, shrubland, and mixed forest; CRO, cropland and grassland

Full size image

Satellite data and phenology retrieval

The MODIS MCD43A4 NBAR dataset (Collection 5, Schaaf et al. 2002) at 500-m resolution and 8-day interval was used to compute the PPI for the period of 2000 to 2016 over the study area. The PPI was formulated using the difference between near-infrared reflectance (NIR) and red reflectance (red):

$$ PPI=-K\times \ln \frac{DV{I}_{\mathrm{max}}- DVI}{DV{I}_{\mathrm{max}}-0.09}, $$
(1)

where DVI is the difference vegetation index (DVI = NIR − red). DVImax represents the maximum DVI of a pixel. K is a gain factor depending on vegetation structure, diffuse fraction of solar radiation and sun zenith angle. (For details of PPI, see Jin and Eklundh 2014). PPI has been shown having a linear relationship with green leaf area index and strong correlation with gross primary productivity, and is a useful proxy for vegetation growth dynamics (Jin and Eklundh 2014). The raw PPI data were smoothed with double logistic fitting using the TIMESAT software (Jönsson and Eklundh 2004). Based on previous studies (Jin et al. 2017; Karkauskaite et al. 2017), the thresholds for SOS and EOS (end of growing season) were set at 20% of the mean 17-year amplitude of a pixel. The phenology metrics from 20% PPI threshold have been shown to be in good agreement with those derived from carbon flux towers and direct phenology observations (Jin et al. 2017). In this study, we only focused on SOS trends and climate sensitivities. The EOS was used to divide the climate variables in the full annual cycle before the average SOS into three periods: spring—the 3 months preceding the SOS; winter—from the previous EOS to the beginning of spring; and summer—the growing period from the SOS to the EOS. Since the three periods were determined from the 17 year’s average phenology of a pixel, the window period was the same for all years of the pixel and different among pixels. The phenology-based three-season division (no autumn) was determined from the partial correlation of the SOS to the climate of the individual month before the average SOS (Suppl. Material Fig. S1) and differs from a calendar-based four-season division. The direct effect of spring temperature and precipitation on SOS and the legacy effects of these climate variables (Gordo and Sanz 2010) during previous winter and summer were all investigated in this study.

Climate data and analysis

The gridded daily mean temperature data at 0.22° resolution from ENSEMBLES gridded observational dataset (E-OBS, version 16.0, Haylock et al. 2008) generated by the European Climate Assessment & Dataset project (https://www.ecad.eu/) were used in this study. The precipitation data in the E-OBS are incomplete over our study region and period; therefore, we instead used the NOAA CPC (Climate Prediction Center) unified gauge-based analysis of global daily precipitation from US National Center for Atmospheric Research (NCAR, ftp://ftp.cdc.noaa.gov/Datasets). We resampled the NCAR precipitation data from 0.5° to 0.22° resolution using bilinear interpolation to match the E-OBS grid size. We then calculated the mean temperature and total precipitation for each of the three periods for the years 1999 to 2016. Only the winter and summer periods were calculated for 1999, and only the spring period for 2016, so as to relate these climate variables to the SOS from 2000 to 2016 using panel data analysis. These climate variables had some temporal tendencies in the 17 years, but in general, there were no significant trends (Suppl. Material Fig. S2 and Table S1).

Satellite data analysis

We used the fixed-effect panel analysis technique (Hsiao 2003) to estimate SOS trends, climate sensitivities, and partial correlations between SOS and climate variables in the periods before the mean SOS. In panel analysis, data of many independent pixels within a predefined area are pooled together to increase the sample size and obtain stronger statistical inference. Panel data regression is not a simple regression on pooled data; instead, it is a multivariate regression on a common slope for all pixels and different intercepts for individual pixels. The pixel-specific intercept from the fixed effect accounts for heterogeneity within the panel that is not explained by the common slope.

Here we pooled 17-year SOS time series of all MODIS pixels (500-m pixel size) within a climate grid (0.22° pixel size, non-overlapped and covering about 3000 pixels on average) to estimate SOS trend, climate sensitivity, and SOS-climate partial correlation. The sensitivity was determined as the changes in SOS per unit change in climate factors from a linear model on the first-order differences of variables following Zhou et al. (2001), so as to avoid exaggerated statistical significant relations from spurious regression (Granger and Newbold 1974). The partial correlation was also estimated on the first-order differences of variables to avoid spurious correlation. The panels covering less than 20 valid MODIS pixels (non-vegetated area) were excluded from further processing. Ordinary least squares (OLS) regression was used to estimate the slope from the panel data, which gave the phenology trend (or sensitivity) of the panel. The significance level of the slope was estimated from heteroscedasticity and autocorrelation-corrected variance of the slope (Wooldridge 2010; Baltagi et al. 2011). The overall slope of the studied region was summarized from a meta-analysis of the slopes using the variance-weighted least squares method (Becker and Wu 2007). (For further information about the data analysis, see Suppl. Material).

Results

SOS trend

The panel analysis of satellite-derived SOS trends showed that 78.6% of the grids in the study area had significantly negative trends (p ≤ 0.05), 17.0% of the grids had significantly positive trends (p ≤ 0.05), and 4.4% of the grids had no significant trends (p > 0.05, Fig. 2). The standard error and the significance level of SOS trend are shown in Suppl. Material Fig. S3. The proportion of grids with significant trends is shown in Suppl. Material Table S2 for each land cover class. Overall, the SOS had advanced − 0.30 ± 0.01 day year−1 during 2000–2016, given from the meta-analysis of regression slope (Table 1). The ENF had experienced a relatively large SOS advancement (− 0.47 ± 0.02 day year−1) and the CRO a small SOS advancement (− 0.14 ± 0.02 day year−1) from 2000 to 2016.

Fig. 2
figure 2

Map and histogram of SOS trend (day year−1) during 2000–2016. Only significant trends (p ≤ 0.05) are shown, with red color (negative value) denoting the advancement of SOS and blue color (positive value) denoting the delay of SOS

Full size image
Table 1 SOS trend and sensitivities (mean ± SE) to temperature and precipitation of the preceding periods during 2000–2016
Full size table

SOS partial correlation with temperature and precipitation

The mean air temperature in spring had on average a negative correlation (− 0.38 ± 0.17) with the SOS during the years 2000–2016 (Fig. 3, Table 2). The winter temperature had a positive correlation (0.09 ± 0.20) with the SOS, with the partial correlation coefficients spanning from − 0.4 to + 0.5. The temperature during the summer growing season, in general, had a weak positive correlation with the SOS (0.05 ± 0.18). The spring precipitation had a positive correlation with the SOS (0.10 ± 0.17). The summer precipitation had a slight positive correlation with the SOS (0.01 ± 0.18). These partial correlation coefficients were all significant at p ≤ 0.05 (two-tailed t test, N > 2500). The winter precipitation had on average no correlation with the SOS (p > 0.05).

Fig. 3
figure 3

Violin plots of partial correlation of the SOS to a temperature and b precipitation during 2000–2016. The shaded areas of the violin plot denote distribution, the white dots denote the median values, and the bottom and top edges of the black bars denote the 25th and 75th percentiles respectively. Spring—the 3 months preceding SOS; winter—from the EOS to the start of spring; summer—from the SOS to the EOS

Full size image
Table 2 Partial correlation coefficients (mean ± std) of SOS to temperature and precipitation during three periods preceding SOS
Full size table

Among individual land cover types, the correlations between spring temperature and SOS did not differ much, with values from − 0.36 ± 0.19 for ENF to − 0.40 ± 0.16 for CRO (Table 2). Contrastingly, the correlations between winter temperature and SOS varied greatly among different land cover types. In ENF, the winter temperature positively correlated with SOS (0.18 ± 0.16), while in CRO, the winter temperature had on average a weak negative correlation with the SOS (− 0.03 ± 0.21). The average positive correlation between spring precipitation and SOS was observed across all land cover types, with the strongest positive correlation (0.15 ± 0.19) in WET.

SOS sensitivities to temperature and precipitation

The meta-analysis of panel data regression slopes showed that, across all land covers, the SOS was sensitive to spring temperature at a rate of − 2.47 ± 0.04 day °C−1 during the years 2000–2016 (Fig. 4 and Table 1). The SOS was positively sensitive to the winter temperature at a rate of 0.66 ± 0.04 day °C−1 and to the summer temperature at a rate of 0.56 ± 0.07 day °C−1. The SOS was sensitive to the spring precipitation at a positive rate of 0.18 ± 0.02 day cm−1, whereas the overall sensitivities of SOS to winter and summer precipitation were − 0.02 ± 0.01 and 0.02 ± 0.01 day cm−1 respectively (Table 1). However, the spatial heterogeneity is large. For example, the SOS of Swedish boreal forests responded to winter and summer precipitation differently to that of Finnish boreal forests (Fig. 5).

Fig. 4
figure 4

SOS sensitivities to a spring temperature (day °C−1) and b spring precipitation (day cm−1) during 2000–2016. Only significant sensitivities (p ≤ 0.05) are shown, with red color (negative value) denoting that the increase of the temperature or precipitation has led to advancement of SOS and blue color (positive value) denoting that the decrease of the temperature or precipitation has led to advancement of SOS

Full size image
Fig. 5
figure 5

SOS sensitivities to a winter temperature (day °C−1) and b winter precipitation (day cm−1), c summer temperature (day °C−1), and d summer precipitation (day cm−1) during 2000–2016. Only significant trends (p ≤ 0.05) are shown, with red color (negative value) denoting the advancement of SOS and blue color (positive value) denoting the delayed SOS

Full size image

Among the land cover types and periods, the SOS of WET had the strongest spring temperature sensitivity (− 3.21 ± 0.17 day °C−1) and the strongest spring precipitation sensitivity (0.50 ± 0.06 day cm−1). The SOS of each individual land cover type was positively sensitive to winter and summer temperature except for the CRO, which had negative sensitivities to both winter and summer temperature (Table 1). The SOS of CRO had an overall small spring precipitation sensitivity (0.09 ± 0.02 day cm−1), contrasting to other land cover’s spring precipitation sensitivity from 0.18 day cm−1 for DBF to 0.50 day cm−1 for WET.

The SOS had significantly negative spring temperature sensitivity over most of the study area (96.9% of all grids, p ≤ 0.05), and only 2.2% of the grids had positive sensitivity (p ≤ 0.05), which summed up to 99.1% of all grids that significantly responded to spring temperature (Fig. 4a and Suppl. Material Table S2). The SOS responded to winter and summer temperature heterogeneously in the spatial domain, and the sensitivity values switched sign from positive in most of the northern part of the study area to negative in the southern part (Fig. 5a, c). The spatially heterogeneous responses of SOS to precipitation were also seen in all periods (Fig. 4b and Fig. 5b, d).

Latitudinal patterns of SOS trend and sensitivities

The SOS trend and its climate sensitivities had no monotonic relationship to latitude (Fig. 6). Below 55° N, the advancement of SOS (negative trend) decreased from − 0.26 day year−1 at 50° N to about zero at 55° N with increasing latitude (Fig. 6a). Above 55° N, the SOS advancement increased up to − 0.62 day year−1 at about 65° N with increasing latitude. Above 65° N, the latitudinal mean SOS trend hovered about − 0.39 day year−1 along with the variation of latitude.

Fig. 6
figure 6

Latitudinal patterns of SOS trends and climate sensitivities in preceding spring, winter, and summer during 2000–2016. a SOS trend (day year−1), b temperature sensitivity (day °C−1), c precipitation sensitivity (day cm−1). Curves show variance-weighted latitudinal mean value. Shaded areas show one standard deviation of the mean

Full size image

The latitudinal mean sensitivity of SOS to spring temperature was around − 2.41 day °C−1 for the latitudes from 50 to 62° N (Fig. 6b). Above 62° N, the spring temperature sensitivity increased up to − 5.44 day °C−1 at around 70° N with a local peak of − 4.41 day °C−1 at 64° N. The winter temperature sensitivities were positive at latitudes 57–69° N, and negative at latitudes below 57 and above 69° N. The strongest positive winter temperature sensitivity was 2.35 day °C−1 at 67° N, and the strongest negative sensitivity was − 1.07 day °C−1 at around 54° N. The latitudinal pattern of summer temperature sensitivities was similar to that of winter temperature sensitivity except that the former was stronger with a larger latitudinal spread and was positive at latitudes above 57° N.

The SOS spring precipitation sensitivities were around zero at latitudes below 58° N (Fig. 6c). With latitude increasing from 58 to 70° N, the spring precipitation sensitivity increased from zero up to 0.80 day cm−1, with a local peak sensitivity of 1.15 day cm−1 at around 67° N. The SOS winter precipitation sensitivities hovered at zero along the whole latitudinal profile, with a noticeable positive sensitivity at latitudes above 69° N. The SOS summer precipitation sensitivities were around zero except for latitudes above 67° N.

Discussion

Spring phenology trend

We found that the spring phenology trends were spatially diverse and on average advanced at a rate of − 0.30 day year−1 at northern European latitudes for the period 2000–2016. However, the spatial heterogeneity was large, spanning from − 1.5 day year−1 to + 1.0 day year−1. Our trend estimation is comparable to the ground-observed phenology in Europe, e.g., − 0.25 day year−1 by Menzel et al. (2006) over Europe during 1971–2000 and − 0.38 day year−1 by Fu et al. (2014a) over West and Central Europe during 2000–2011. By contrast, other studies using NDVI have reported inconclusive or insignificant trends during the beginning of the new century (e.g., Karlsen et al. 2014; Zhang et al. 2014; Zeng et al. 2011) or a delayed SOS onset of spring phenology that reversed the advancement of preceding decades (Fu et al. 2014a). Many factors could lead to these discrepancies. Firstly, our study has overlapped but not identical spatial coverage, temporal span, and vegetation types with these studies. Secondly, different spectral indices may show different land surface and vegetation processes. The NDVI-based land surface phenology at northern latitudes was shown to relate to cryospheric seasonality, such as occurrence of snow (Jin et al. 2017; Jönsson et al. 2010) or ice melting (White et al. 2009), whereas the PPI-based phenology to relate to plant photosynthetic dynamics (Jin et al. 2017; Karkauskaite et al. 2017). Thirdly, different statistical methods could also lead to different results. It has been shown difficult to estimate significant phenology trends from short and noisy time series using pixel-wised satellite data (e.g., Karlsen et al. 2014; Zhang et al. 2014; Zeng et al. 2011), whereas a panel analysis can reach statistically significant estimates by pooling data together and increasing the degree of freedom (Keenan et al. 2014). Furthermore, most studies summarized the phenology trend from simple average or median values; the meta-data analysis used in our study considered data quality is also an advantage. The new estimation of satellite-based phenology trend shown here provides new information on plant phenology at the northern European latitudes and another option for remote sensing phenology studies.

It should be noted that we only considered the temporal autocorrelation of individual pixels and assumed independence of pixels in a panel. The potential spatial autocorrelation and contemporary spatial-time lag effects of neighboring pixels in the panel that violate the assumption could be tested and addressed further using more complex spatial dynamic panel data modeling (Lee and Yu 2011), which is out of the scope of this study. Caution should be taken when explaining the statistical significance to avoid confusing with biological significance.

Alongside the majority area (78.6% of the grids) experiencing an advancement of SOS during 2000–2016, we detected locally delayed SOS onsets (significant positive trends, 17.0% of the grids), mainly in the northern border area between Norway and Sweden, southern Norway, and Lithuania. Our estimation confirms the positive trend in parts of Fennoscandia reported by Høgda et al. (2013) using NDVI and local experience of birch budburst dates, which was suggested to relate to the increased winter precipitation at high-latitude mountains of the oceanic region (Høgda et al. 2013), because more heat accumulation was required for leafing with increasing precipitation (Fu et al. 2015a). In Lithuania, there were similar decadal positive trends for the years 1971–2000 shown based on ground observations (Kalvāne et al. 2009). The delayed spring events were also shown on winter rye in Poland in the 2000s (Blecharczyk et al. 2016). These spatially heterogeneous phenology trends imply that it is essential to use consistent and reliable wall-to-wall phenology estimates to comprehensively evaluate climate change impacts on ecosystems in order to create an unbiased estimate.

SOS correlation with temperature and precipitation

Our study showed strong negative partial correlations of SOS with spring temperature and positive correlations with spring precipitation. The earlier SOS onset with warmer temperature (negative correlation) is consistent with the warming-promoted advancement of bud break and leafing (Singh et al. 2017; Hänninen 2016). That more spring precipitation delays the growing onset (positive correlation) agrees with Fu et al. (2015a) using long-term observations from Europe and North America, and this explains that more heat accumulation was required for leafing with increasing precipitation. The spatial heterogeneity of the SOS relationship to spring precipitation reveals that the precipitation effects on the growth onset are complex and may entangle with soil moisture condition (Piao et al. 2006) and snowmelt date (Böttcher et al. 2014; Yun et al. 2018). Disentangling these effects on phenology is important for further understanding of the phenology process but out of the scope of this study.

The climate during the full annual cycle before spring was reported to relate to the phenology onset of different strength (Gordo and Sanz 2010; Fitter et al. 1995). Here, we found that the temperature and precipitation during winter and summer indeed affected the SOS, but in general weaker than the direct effect of spring climate on SOS. The winter and summer temperature had on average a positive correlation with the SOS, whereas the large spread (standard deviation) of these partial correlation coefficients indicated a large spatial heterogeneity of the SOS responses to climate variables during the different periods preceding the SOS. In addition to the spatial heterogeneity of climate, the spatial heterogeneity of SOS climate sensitivities explains the spatial heterogeneity of spring phenology trends at northern European latitudes. It should be noted that our method of dividing a year cycle into three periods may overlook possible opposite carry-over effects of individual months in a season (Gordo and Sanz 2010). Also, the forest classes in reality consist of several species with individual phenology responses, and our CRO class includes pastures and grasslands and is not one homogenous class. Furthermore, simply averaging correlations over a whole region will cancel out the positive correlation in one area with the negative correlation in another irrespective of data quality, and may result in no significant correlation over the region. In this regard, the meta-analysis of sensitivities in this study will generate more robust estimation over a whole region by considering individual data quality.

SOS temperature sensitivity

The SOS temperature sensitivity of − 2.47 day °C−1 in our estimation agrees well with the value of − 2.5 day °C−1 by Menzel et al. (2006) of leaf unfolding in Europe during 1971–2000, but is slightly larger than the values of − 2.3 day °C−1 by Fu et al. (2015b) of European tree leafing during 1999–2013 and − 2.21 day °C−1 by Wolkovich et al. (2012) from warming experiments of leafing dates mostly at northern middle and high latitudes. Other studies have reported considerably larger SOS temperature sensitivities. For example, Chmielewski and Rötzer (2001) showed an SOS advancement of 7 day °C−1 with a mean temperature of February to April in Europe based on ground phenological observations. Wolkovich et al. (2012) reported an apparent temperature sensitivity of − 6.42 day °C−1 from long-term observed leafing in the United States and Europe. Using the land surface greenness estimated from NDVI, Piao et al. (2006) reported a greenness onset sensitivity of − 7.5 day °C−1 to the mean temperature in March to early May in China during 1982–1999.

Apart from the fact that different study areas and time spans may lead to different temperature sensitivity estimation, other factors can also contribute to these inconsistencies. Firstly, land surface greenness onset estimated from NDVI does not necessarily match the plant canopy growth onset. Snowmelt concurs largely with the greenness onset over high northern latitudes (Delbart et al. 2005; Jönsson et al. 2010) but may differ in other areas. Secondly, an improper preseason window may bias the temperature sensitivity estimation over a large region. In Chmielewski and Rötzer (2001), to avoid circular reasoning, a common preseason of February to April was used for the SOS in March and May in Europe. This is not appropriate for all sites over the large region because the April temperature at a site has no influence on the SOS if it occurs in March. Thirdly, a linear model of SOS to temperature generates an estimation of apparent sensitivity (Wolkovich et al. 2012), which attributes all SOS changes to spring temperature regardless of other factors or deterministic trends in SOS (Zhou et al. 2001). Lastly, statistic bias by either considering only significant results or simply averaging of all results without considering data quality will generate improper summarization. In this study, we used the physically based vegetation index PPI to estimate vegetation spring growth onset, and related the SOS variations to all seasonal climate variable changes using panel analysis and partial correlation. A site-specific preseason window fixed for all 17 years of the site was used in a local regression model for sensitivity estimation to avoid both improper preseason bias and circular reasoning. We summarized the overall sensitivities of the area using a meta-analysis of regression slopes, which considered both the sample number and the precision of individual regression. The overall sensitivity estimation is more robust than a simple average or biased estimation by neglecting insignificant regression results (e.g., Garonna et al. 2014).

The warming winter would lead to a delayed SOS onset in most of the Fennoscandian region (Fig. 5a), which can be explained by the decreased chilling accumulation for woody plant dormancy break at high latitudes. Many studies have shown that warming winter resulted in delayed spring phenology onsets at northern latitudes (e.g., Fu et al. 2014a; Zhang et al. 2007; Murray et al. 1989) due to insufficient chilling accumulation required to break winter dormancy of trees. In croplands and grasslands at the southern part of the study region, the winter warming–promoted earlier spring growth recovery may evidence that no chilling accumulation was required to break winter dormancy for these plants. Here, we also showed that summer temperature was important for the coming spring phenology. One explanation is that the high temperature in later summer or autumn could reduce soil moisture availability later and delay the spring events like leafing (Chmielewski et al. 2005) and flowering (Fitter et al. 1995). Also, later summer temperature may affect bud dormancy initiation and consequently have legacy effects on spring bud break (Gordo and Sanz 2010). These unclear mechanisms merit further studies considering that more summer warming and extreme events are projected in future at European latitudes (IPCC 2014).

SOS precipitation sensitivity

A few studies have reported precipitation sensitivity of phenology over the northern latitudes. Piao et al. (2006) showed diverse precipitation influences on the growing season depending on vegetation type and growth conditions in China. Cong et al. (2013) found that precipitation may regulate the temperature-determined spring phenology to some extent over temperate China. In this study, we found that, among the diverse spatial patterns of spring precipitation sensitivity, the spring phenology advanced along with decreasing precipitation (positive sensitivity), opposite to the negative temperature sensitivity (increasing temperature induced SOS advancement). This may be explained by promoted growth onset due to increased solar radiation and temperature accompanying decreasing precipitation over Nordic boreal forests, in which plant growth is generally not water limited (Bergh et al. 1999) but light limited (Singh et al. 2017). Wetter winters were also found to delay the SOS onset (positive sensitivity) in Swedish boreal forests (Fig. 5b), which was related to the decreased heat flux and accumulated growing degree days by winter precipitation mostly in form of snow in boreal forests (Yun et al. 2018). We found that dry summers could delay the SOS onset (negative sensitivity) in Swedish boreal forests (Fig. 5d), which was probably related to dry summer–induced soil moisture deficit during the growth onset (Seyednasrollah et al. 2018), a similar legacy effect of warmer summer on the delayed SOS as mentioned before.

The linkage of precipitation and vegetation phenology depends on soil water potential, the timing of precipitation events, and plant types, rather than precipitation amount (Jolly and Running 2004; Peñuelas et al. 2004). We note here the importance of precipitation on growth onset, while the complicated spatial pattern of precipitation sensitivity revealed here calls for further studies. Species-specific phenological variations and possibly non-linear relationships between the SOS and climate variation may result in substantial uncertainties.

Latitudinal patterns

The latitudinal pattern of the SOS trend in our study shares some common properties with that in Garonna et al. (2014) of Europe during 1982–2012, such as the near-zero SOS trend at about 55° N and the decreasing SOS advancement with increasing latitude after 65° N. However, due to the mismatch of the study area and time span, the latitudinal pattern in our study does not exactly match that in Garonna et al. (2014) over the short study period, during which the temperature and precipitation had no overall significant trends in our study (Suppl. Material Fig. S2 and Table S1), and were largely affected by complicated atmospheric circulation patterns (de Beurs and Henebry 2008; de Beurs et al. 2018). The SOS trend value in our estimation was about three times the value reported by Garonna et al. (2014), more prominent at latitudes over 60° N. This is probably because the phenology estimation from different spectral indices represents different land surface and vegetation processes as mentioned before. Moreover, our meta-analysis of regression trend is advantageous and statistically sound compared to Garonna et al. (2014), who considered only significant trends and discarded insignificant trends in their analysis. The statistic bias report has been criticized by Menzel et al. (2006) and is still ignored by researchers.

The spring temperature sensitivity was low and stable at latitudes below 62° N. Larger sensitivities were seen at higher latitudes. This finding supports to some degree that high sensitivity was found in higher latitudes or colder area (Güsewell et al. 2017; Seyednasrollah et al. 2018), contradictory to the viewpoint that large sensitivity occurred at low latitudes (Wang et al. 2015a). The positive winter temperature sensitivities at latitudes 57–69° N were probably attributed to high chilling requirements for dormancy break in temperate and boreal forests towards higher latitudes (Fu et al. 2015b; Singh et al. 2017; Hänninen 2016). Apparently, the direct influence of spring precipitation on the SOS was stronger than the legacy effects of winter and summer precipitation, whereas the spatially complex and heterogeneous precipitation sensitivity could not be shown in the latitudinal summarization and was approximately zero. At latitudes above 69° N, the negative winter temperature sensitivities and the concurrent positive winter precipitation sensitivity showed that a wetter and colder winter would delay the spring onset in sub-Arctic areas, while the entangled effect of snowmelt on sub-Arctic spring phenology onset remains to be studied.

Conclusions

Using remote sensing and robust statistical methodology, we mapped significant phenology trends and its complex temperature and precipitation sensitivities over most parts of northern European latitudes (> 50° N) for the period 2000–2016. The SOS trend (− 0.30 day year−1) and temperature sensitivity (− 2.47 day °C−1 to spring warming) differ from those estimated from land surface greenness onset and apparent sensitivity estimation but are consistent with those from ground-observed phenology with careful analyses. The generated maps show that regional variations in phenology trends and their sensitivities to driving factors in different periods and areas may be strong. We conclude that the methodology is suitable for generating spatially continuous large-area estimates and climate sensitivity analyses. Our study brings a new perspective on climate influences on phenology and demonstrates a useful application of the newly developed vegetation index PPI and advanced statistics for improved assessment of phenology variations and its climate sensitivities over large areas.