## Materials and Methods
### 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.