## 2 Materials and Methods
### Preliminary Analysis and Optimisation of the Extraction Method
#### 2.1.3 Fitting Models and Extraction Methods  
In order to perform an accurate satellite analysis of phenological date trends over the period 2001-2021, different fitting models (FMs) and date extraction methods (EMs) were applied both to raw data of GPP and kNDVI to find the most performing FM\(\times\)EM approach. The methodology was then chosen taking into account the results of the comparison between SOS, POS and EOS extracted from GPP and kNDVI annual patterns.  
Consequently, raw data of GPP and kNDVI underwent a process for retrieving the phenological dates of SOS, POS and EOS. Specifically, using a work package for phenological analysis within the R work environment (_phenopix_, [17]), the annual patterns of GPP and kNDVI were subjected to the fitting operation and, subsequently, to SOS, POS and EOS extraction.  
\begin{table}
\begin{tabular}{c c c c c c c c} \hline ID & **Site** & **Country** & **Network** & **Lat** & **Lon** & **Mamsl** & **Years** \\ \hline AT-Neu & Neusitft [41] & Austria & FLUXNET & 47.1167 & 11.3175 & 970 & 2002–2012 \\ CH-Cha & Chamau [42] & Switzerland & FLUXNET & 47.2102 & 8.4104 & 393 & 2002–2008 \\ CH-Fru & Früebel [43] & Switzerland & FLUXNET & 47.1158 & 8.5378 & 982 & 2005–2014 \\ CZ-BKZ & Bily Kriz & Czech Republic & FLUXNET & 49.4944 & 18.5429 & 855 & 2006–2012 \\ DE-Gri & Grillenburg [44] & Germany & FLUXNET & 50.9500 & 13.5126 & 385 & 2004–2018 \\ DE-Rur & Rollesbroich [45] & Germany & FLUXNET & 50.6219 & 6.3041 & 514 & 2011–2018 \\ IT-Mal & Malga Arpaco & Italy & EFDC & 46.1140 & 11.70334 & 1662 & 2003–2004 \\ IT-Mbo & Monte Bondone [46] & Italy & FLUXNET & 46.0147 & 11.0458 & 1550 & 2004–2013 \\ IT-Tor & Torgnon [25] & Italy & FLUXNET & 45.8444 & 7.5781 & 2160 & 2009–2018 \\ \hline \end{tabular}
\end{table}
Table 1: List of grassland sites used to obtain daily GPP data. Mamsl represents meters above mean sea level (m), while EFDC represents the European Fluxes Database Cluster.  
The FM used in this comparison were 4: [PERSON] et al. [48] (ELM), [PERSON] et al. [49] (GU), [PERSON] et al. [50] (BEC) and [PERSON] et al. [51] (KLS). These methods fitted different double logistic curves to the raw data of GPP and kNDVI according to the aforementioned studies. Equations of ELM, GU, BEC and KLS [17] are shown below:  
\[f(t)=mn+(mx-mn)\cdot\left(\frac{1}{1+e^{(m_{3}^{2}-t)/m_{4}^{ \prime}}}-\frac{1}{1+e^{(m_{5}^{2}-t)/m_{6}^{\prime}}}\right) \tag{4}\] \[f(t)=y_{0}+\frac{a_{1}}{\left[1+e^{-(t-t_{01})/b_{1}}\right]c_{1} }-\frac{a_{2}}{\left[1+e^{-(t-t_{02})/b_{2}}\right]c_{2}}\] (5) \[f(t)=mn+(mx-mn)\cdot\left(\frac{1}{1+e^{(-sp\cdot(t-\cos))}}+ \frac{1}{1+e^{(-rau\cdot(t-\cos))}}\right)\] (6) \[f(t)=(a_{1}t+b_{1})+\left(a_{2}t^{2}+b_{2}t+c\right)\cdot\left( \frac{1}{\left[1+q_{1}e^{-h_{1}(t-n_{1})}\right]v_{1}}-\frac{1}{\left[1+q_{2}e ^{-h_{2}(t-n_{2})}\right]v_{2}}\right) \tag{7}\]  
The selected FMs optimise a different number of parameters and hence present diverse flexibility in fitting raw data. For a more comprehensive understanding of the methods and curve parameters, see the corresponding publications.  
The GPP and kNDVI curves obtained after FM application were then used to extract SOS, POS and EOS with four extraction methods [17]: three working on inflection points of the derivatives (_Klosterman, Gu, Derivatives_), and one that identifies the dates when a fixed threshold of the seasonal amplitude is reached (_Tresholds_). Specifically, _Klosterman_ is based on local extremes in the rate of change of curvature k [52], _Deritatives_ on local extremes in the first derivative, and \(Gu\) on a combination of local maxima in the first derivative [49]. Regarding the _Tresholds_ method, in this trial we tested different values for the fixed threshold: 10 (TRS\({}_{0.1}\)), 20 (TRS\({}_{0.2}\)), 30 (TRS\({}_{0.3}\)), 40 (TRS\({}_{0.4}\)) and 50% (TRS\({}_{0.5}\)) of the seasonal amplitude.