## 3 Methods
### Empirical Mode Decomposition  
Empirical Mode Decomposition (EMD) is a signal processing method that aims to uncover the internal structure of data by decomposing a signal into its sub-components at different frequencies [14]. EMD has been used for various applications, such as analyzing time series data, identifying periodic or continuous structures in data, and decomposing signals in the frequency domain. These components are called Intrinsic Mode Functions (IMFs). Each IMF represents the local frequency modulation of the data and reflects the time-variance relationship in the signal. IMFs preserve the intrinsic characteristics of the data while determining their fundamental frequencies. EMD decomposes a signal into IMFs through a series of progressive steps. For a signal component to be considered as an IMF, it must satisfy two conditions: (i) the number of local maxima and minima and the number of zero crossings must be the same or differ at most by one and (ii) the arithmetic mean of the envelope defined by the local maxima and the envelope defined by the local minima must be zero. The sum of the IMF and residual mode components, determined through the decomposition and elimination processes of a signal by EMD, reconstructs the original signal. EMD calculation steps are defined as follows. Step 1: All extrema in the time series are identified. The local maxima and minima are then used for interpolation and the upper \(U_{x}(t)\) and lower \(L_{x}(t)\) envelopes are derived using cubic spline. Step 2: The mean envelope value \(m(t)\) and the detailed component \(d(t)\) are calculated. \[m(t)=\frac{U_{x}(t)-L_{x}(t)}{2} \tag{4}\]  
\[d(t)=x(t)-m(t) \tag{5}\]  
Step 3: If the resulting \(d(t)\) meets the stopping criteria, it becomes the first IMF. Otherwise, the procedure is repeated on \(d(t)\) until the stopping criteria are met, where the stopping criteria are defined as follows:  
\[\sum_{t=1}^{T}\frac{\left|d_{j}(t)-d_{j-1}(t)\right|^{2}}{d_{j-1}^{2}(t)}\leq SD \tag{6}\]  
where \(T\) represents the length of the time series and \(j\) denotes the number of iterative calculations. The default value of SD is usually set between 0.2 and 0.3. Thus, in this study, the default SD value was set to 0.2. Step 4: Continue iterating steps two to three until all IMFs and the detailed signal are obtained. Finally, the original time series \(f(t)\) can be decomposed as given below:  
\[f(t)=\sum_{i}^{N}IMF_{i}(t)+R_{N}(t) \tag{7}\]  
where \(R_{N}(t)\) is the final residual and can be considered as the global trend signal. The residual band obtained through the EMD method is commonly referred to as the long-term trend component [88], which identifies changes in the nonlinear trend present in the time series over time. In the literature, its application is noted in various fields such as precipitation [89], air temperature [90], extreme sea levels [91], and regional drought [92]. We decomposed the GPP time series into sub-series with different frequencies using the EMD method. This decomposition process utilized the EMD library [93; 94] in R. While the EMD algorithm performs a data-based decomposition, it still contains certain hyperparameters including the maximum number of siftings and the maximum number of IMFs. The maximum number of siftings refers to the maximum iterations the algorithm will perform during the extraction of each IMF, which influences in how much detail the algorithm will process each IMF and how long it will run. A larger maximum sifting number can lead to longer extraction processes for each IMF. This can create a significant computational cost, especially for large datasets. Excessive sifting may sometimes result in the over-sifting of the signal, causing the IMFs to be overly detailed and to fail to reflect the main signal's characteristics [95]. Depending on the specific application of the algorithm and the characteristics of the dataset, the optimal value of this hyperparameter can vary.