## Materials and Methods
### Time-Series Smoothing Methods
#### 2.2.2 Spline Smoothing (SP)  
We adapted the Fortran code by [PERSON] [64] to achieve spline smoothing based on the description in Craven and Wahba [63]. The Whittaker smoother [66; 67] can be considered as a special case of spline smoother [68]. Spline smoothing is employed to find a spline function \(S_{p}(t)\) that minimizes the value of a criterion function \(C_{p}\) for a selected smoothing parameter \(p\geq 0\),  
\[C_{p}=\sum_{i=1}^{n}\biggl{\{}w_{i}\bigl{[}y_{i}-S_{p}(t_{i})\bigr{]}^{2}+v_{i} p\int|S_{p}^{\prime\prime}(t)|^{2}dt\biggr{\}} \tag{2}\]  
where \(t_{i},i=1,\ldots,n\) is a regularly spaced time grid and \(y(t_{i})\equiv y_{i}\) the corresponding NDVI values. Each point in the time-series is associated with a weight \(w(t_{i})\equiv w_{i}\). \(v\) is a vector containing an adaptive stiffness weight for each \(i\). The smoothing parameter \(p\) controls the shape of the spline, varying from an exactly interpolating spline (\(p=0\)) to a straight line (\(p\rightarrow\infty\)) (Figure 5). Instead of manually selecting suitable \(p\), an alternative automatic approach is to use generalized Cross-Validation (CV), which was described by [PERSON] and [PERSON] [63] and implemented in the Fortran code based on the improvement by [PERSON] [64]. In this study, we tested SP with CV (SPCV) and without CV, respectively. We also implemented CV in our code for the SG and LO methods to automatically choose the \(q\) parameter and to be able to compare the effect of CV on these methods with the effect on SP. The range of \(p\) in SP was limited to produce local fits, which means that we forced the performance of SP to be fairly close to SG or LO. In this case \(p\) was set to the range 0.1-10 with exponential intervals, i.e., 0.1, 0.25, 0.63, 1.58, 3.98 and 10 (Table 2). \begin{table}
\begin{tabular}{c c c c c c} \hline
**Smoothing** & **Upper-Envelope** & **Upper-Envelope** & **Smoothing** & **Stiffness** & **Total Number** \\
**Method** & **Iterations (\(u\))** & **Strength (\(a\))** & **Parameter (\(q\)/\(\log(p)\))** & **Strength (\(l\))** & **of Settings** \\ \hline SG & 1–3 & 1–10 & 2–7 & - & 126 \\ LO & 1–3 & 1–10 & 4–12 & - & 189 \\ SP & 1–3 & 1–10 & (\(-\)1)–1 & 0.1–0.9 & 630 \\ DL & 1–3 & 1–10 & - & - & 21 \\ AG & 1–3 & 1–10 & - & - & 21 \\ \hline \end{tabular} Note: \(u\): number of iterations of upper-envelope, \(a\): adaption strength of upper-envelope iterations (only valid for \(u>1\)), \(q\): half window size, \(p\): smoothing parameter, \(l\): adaption strength of spline stiffness.  
\end{table} 
Table 2: Number of parameters of each smoothing method. Altogether, 987 simulations were made for each test site. Differently from the upper-envelope iterations in the other smoothing methods, the upper-envelope iterations in SP and SPCV also updates each \(v_{i}\) based on the formula:  
\[v_{i}=1-l\frac{\|S_{p}^{q}(t_{i})\|}{\|S_{p}\|} \tag{3}\]  
where \(\|S_{p}\|\) is the amplitude of all \(S_{p}\), and \(\|S_{p}^{q}(t_{i})\|\) is the amplitude of \(S_{p}\) in the window from \(t_{i}-q\) to \(t_{i}+q\), in which \(q\) was set to 5, i.e., a quarter of one year's data; \(l\) is the stiffness adaptation strength, which is in the range from 0 to close to 1, where 0 will give no local adaptation and close to 1 will give very high local adaptation. The adaptation factor \(v_{i}\) was integrated into SP to locally adjust the curve shape. In the original code spline stiffness is governed by \(p\). A large value of \(p\) gives a stiff spline, and a low value of \(p\) gives a spline that is flexible. However, a large value of \(p\) may result in low absolute values of the derivative in the spring and the autumn, resulting in poor fits during these seasons. The idea with the adaptive method is to give each point \(i\) a stiffness by multiplying the global stiffness parameter \(p\) with an adaptation factor \(v_{i}\).