## 2 Materials and Methods
### Method  
#### 2.3.1 Phenology Metrics Determination  
**GPP time series reconstruction**. Although the GLASS GPP product was widely validated with ground measurements and shown with high accuracy [7], there were still some outliers resulting from factors such as limitations of the GPP model and environmental variables. To eliminate various interferences during data collection and processing, and \begin{table}
\begin{tabular}{c c c c} \hline \hline
**Type Code** & **Landcover Type** & **Value** & **Major Type** \\ \hline WA & Water & 0 & - \\ \hline ENF & Evergreen needleleaf forest & 1 & \\ EBF & Evergreen broadleaf forest & 2 & \\ DNF & Deciduous needleleaf forest & 3 & forestlands \\ DBF & Deciduous broadleaf forest & 4 & \\ MF & Mixed forest & 5 & \\ \hline CS & Closed shrulands & 6 & \\ OS & Open shrulands & 7 & shrulands \\ \hline WS & Woody savannas & 8 & \\ SAV & Savannas & 9 & grasslands \\ GRV & Grasslands & 10 & \\ \hline PW & Permanent & 11 & \\ CRO & Croplands & 12 & \\ UAB & Urban and built-up & 13 & \\ OVM & Cropland/natural vegetation mosaic & 14 & - \\ SAI & Snow and ice & 15 & \\ BSV & Barren or sparsely vegetated & 16 & \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Reclassification of land cover types in boreal ecosystems. Figure 2: Main landcover types in boreal terrestrial ecosystem. (**a**) The original MCD12C1 landcover types. (**b**) The reclassified landcover types. ensure the high confidence of the complete cycle in vegetation development, we performed the following preprocessing steps. Firstly, the pixels on the original GPP time series dataset covered by bare land, snow, and ice all year round were excluded, and then the Savitzky-Golay (S-G) filter algorithm was used to smooth the time series and reduce the possible outliers. The S-G filter has been widely used for data smoothing and noise removal; it can remove outliers with significant errors and is not limited by spatial or temporal scales. **GPP time series curve fitting. ** To extract the phenology metrics, we used the double-sigmoidal (DS) model to fit the reconstructed GPP time series. The DS model is a double logistic function, in which the amplitudes of the left and right sides of the curve can be inconsistent, and the fitted curve is double S-shaped. The DS function could be expressed as Equation (1):  
\[y(t)=\alpha_{1}+\frac{\alpha_{2}}{1+\exp(-\partial_{1}(t-\beta_{1}))}-\frac{ \alpha_{3}}{1+\exp(-\partial_{2}(t-\beta_{2}))} \tag{1}\]  
where \(y(t)\) is the fitted GPP at a specific DOY (day of year, \(t\)). These parameters of \(\alpha_{1},~{}\alpha_{2},\alpha_{3},\partial_{1},\partial_{2},\beta_{1}, \beta_{2}\) represent different metrics of the fitted DS model. \(\alpha_{1}\) is the background (baseline) of the fitted GPP time series, \(\alpha_{2}\) and \(\alpha_{3}\) represent the amplitude of the green-up and leaf senescence process of the fitted annual GPP curve. The \(\partial_{1}\) and \(\partial_{2}\) represent the transitions curvature (normalized slope coefficients); \(\beta_{1}\) and \(\beta_{2}\) are the midpoints (in DOY) of these transitions for the green-up and leaf senescence processes, respectively. The functional model fitting can effectively reduce the influence of noise in the time series on the overall shape of the curve, providing a reasonable and effective time series for the phenology date extraction. In order to ensure the robustness of GPP annual curve fitting and phenology metric extraction, only the pixels that have high goodness of fit (\(R^{2}>0.8\)) were retained for the extraction in this study. **Vegetation phenology determination. ** There are methods proposed and successfully used for extracting vegetation phonology metrics, such as fixed threshold methods, dynamic threshold methods, and derivative extremum methods [20; 41]. In this study, we chose the derivative extremum method to determine the phenology metrics from the annual GPP time series [48].