## Methods
### Time Series Processing  
Following pre-processing it was necessary to smooth the data further before extracting seasonality metrics. A wide variety of smoothing techniques has been implemented in previous studies, however, there is a considerable lack of clarity in the literature regarding the most appropriate methods to employ [51]. The overall objective is to eliminate as much noise as possible, but to preserve variations in the seasonal curves in the data, e.g., by applying moving average filters or fitting functions [37]. Four different methods were applied here to allow an assessment of the impact of the chosen smoothing methods on the derived phenology estimates. The selected methods were Gaussian function, double sigmoidal function, fast Fourier transform (FFT) and Savitzky-Golay filter based on their successful implementation in previous studies (Table 1). These were applied to the MODIS VI time-series using the R-package Phenex 1.1-9. Phenex facilitates smoothing of time-series data with irregular time intervals on a single-year basis when assigning dates with no value NA [52]. In addition, the best index slope extraction (BISE) algorithm was applied to the NDVI\({}_{\text{prod}}\) data prior to this to further reduce the noise in the data [52; 53]. Default settings were kept in all cases. The four smoothing techniques were also implemented for the NEE flux data and ground measured PAR (daily mean of NEE and PAR values recorded in the satellite overpass time interval of 10:00 a.m. to 12.30 p.m. and normalized between 0 and 1). In Boreal latitudes, PAR exhibits a strong seasonal pattern due to seasonal changes in sun elevation, which could have an impact on VI estimates. Therefore, PAR was included to test if the sun elevation dynamics correlate with NEE estimates of phenology to evaluate if any potential relationships observed between NEE and the VIs occur as a result of this.  
\begin{table}
\begin{tabular}{c c c} \hline
**Smoothing Technique** & **Description** & **Implemented in (Modified Versions)** \\ \hline \multirow{3}{*}{Gaussian function} & Applies symmetrical or asymmetrical & \\  & Gaussian function that has & [PERSON] and [PERSON] [36] \\  & the best fit to a series of data & [PERSON] et al. [54] \\  & points (C implementation). & \\ \hline \multirow{3}{*}{Double sigmoidal function} & Applies a double sigmoidal & [PERSON] et al. [55] \\  & function to data values that has the & [PERSON] et al. [56] \\  & best fit to a series of data & [PERSON] et al. [57] \\  & points (C implementation). & [PERSON] et al. [11] \\  & [PERSON] and [PERSON] [14] \\ \hline \multirow{3}{*}{Fast Fourier transform} & Fast Fourier transform & \\  & (smoothing intensity parameter & \\  & default value = 3) which decomposes & [PERSON] et al. [58] \\  & noise-affected time series into simpler & [PERSON] et al. [44] \\  & periodic signals using sine and & \\  & cosine functions. & \\ \hline \multirow{6}{*}{Savitzky–Golay filter} & Smooth data values with a & \\  & Savitzky–Golay filter, which is based & [PERSON] et al. [59] \\ \cline{1-1}  & on local least squares polynomial & [PERSON] and [PERSON] [51] \\ \cline{1-1}  & approximation (default window size & [PERSON] et al. [9] \\ \cline{1-1}  & setting = 7, degree of fitting & [PERSON] et al. [10] \\ \cline{1-1}  & polynomial = 2 and smoothing & \\ \cline{1-1}  & repetition = 10). & \\ \hline \end{tabular}
\end{table}
Table 1: Selected smoothing methods, which are available in the Phenex package by [PERSON] and [PERSON] [52].