## 2 Materials and Methods
### Data Processing  
#### 2.3.1 Trend Analysis and Turning Point Year Detection  
The trends of NPP during 2001-2014 were calculated at both pixel level and regional level. Due to autocorrelation among the inter-annual time series data, a robust non-parametric Mann-Kendall (M-K) trend analysis [46] was applied. This method did not require the independence and normality of the time series data [47], which has been widely used in trend analysis [19]. Previous studies have reported that the M-K test statistic \(Z\) was approximately normally distributed when the sample size was \(n\geq 8\). A positive or a negative \(Z\) value indicated an increasing or a decreasing trend respectively, which were all monotonic [46]. The formulas for the M-K method are described in detail in Reference [46]. In addition, trends of the NPP were tested at a significance level of \(\alpha=0.05\).  
The Theil-Sen median slope estimator was applied to estimate the rate of change of NPP, which was more appropriate for assessing the rate of change in short or noisy time series [48]. The Theil-Sen median slope was computed as (Equation (9)):  
\[\beta_{i}\ =\ \text{Median}(\frac{x_{j}\ -\ x_{k}}{j\ -\ k})\ \text{for}\ i\ =\ 1,\ldots,\ N \tag{9}\]  
where \(x_{j}\) and \(x_{k}\) are the data values at times \(j\) and \(k\) (\(j>k\)), respectively.  
To detect the timing and magnitude of NPP changes, a linear regression model and a piecewise linear regression model (Equation (10)) were used. This method has been widely used [49; 41] as it can detect potential turning points (TPs) in a trend of time-series data. In addition, TP was limited to the years 2004 to 2010 to avoid obtaining a too-short segment before or after the TP [50]. The maximum number of TPs was specified as one given that too many TPs would make the result more complex and thus create more uncertainty in the understanding of NPP trends [23]. The models used were:  
\[Y=\left\{\begin{array}{ll}\beta_{1}\times t+\beta_{0}+\varepsilon&(t\leq \alpha)\\ \beta_{1}\times t+\beta_{2}\times(t-\alpha)+\beta_{0}+\varepsilon&(t>\alpha) \end{array}\right. \tag{10}\]  
where \(Y\) is the NPP; \(t\) is the year; \(\alpha\) is the turning point of the NPP time series; \(\beta_{0}\) is the intercept; \(\beta_{1}\) is the magnitude of the NPP trend before the TP; \((\beta_{1}+\beta_{2})\) is the magnitude of the NPP trend after the TP; \(\alpha\) is the year of the TP; and \(\varepsilon\) is the residual random error.  
The two models were then fitted to the NPP time-series data using the least-squares method, and the Akaike Information Criterion (AIC) [51] was used to determine whether the TPs were significant as it provides a means for model selection [41; 51]. The AIC values of these two models were calculated as (Equation (11)):  
\[AIC\ =\ n\ \times\log\ (\frac{RSS}{n})\ +\ 2k\ +\ \frac{2k(k\ +\ 1)}{n\ -\ k\ -\ 1} \tag{11}\]  
where \(RSS\) is the residual sum of squares for the estimated model; \(k\) is the number of parameters; and \(n\) is the sample size.  
Finally, \(\Delta\)AIC was defined as the difference of the AIC2 of the piecewise linear regression model (model 2) and the AIC1 of linear regression model (model 1); if the \(\Delta\)AIC was less than \(-2\), then model 2 was significantly preferred [8].