## 2 Materials and Methods
### Estimation of the Growing Season Length from the NDVI  
For each pixel, the length of the growing season was defined as the time interval between the SOS and EOS dates. To estimate the SOS and EOS from the NDVI data, we adopted the modified Savitzky-Golay filter to smooth out mutational noise in the NDVI time series that is generated by unstable atmospheric states, such as clouds [39]. Because various methods can be used to extract phenology dates from the NDVI time series, the estimated result of each method may be different. Therefore, to prevent the phenology estimation method from affecting the results, we selected four widely used methods to determine the SOS and EOS for the northern hemisphere.  
For Method 1, we first used the double logistic method to fit the NDVI time series at a daily temporal resolution [24; 34; 40], and it can be described as:  
\[\mathrm{NDVI}_{(t)}=a+b(\frac{1}{1+e^{c(t-d)}}+\frac{1}{1+e^{(t-f)}}) \tag{1}\]  
where \(\mathrm{NDVI}_{(t)}\) is the fitted NDVI at day \(t\); \(a\) is the initial background NDVI value; \(a+b\) is the maximum NDVI value; and \(c\), \(d\), \(e\), and \(b\) are parameters of this function.  
The annual variation curve of the NDVI was divided into two sections by the NDVI maximum value: the ascending part (indicating the vegetation recovery from the nongrowing season) in about the first half of the year and the descending part (indicating the process of plant senescence from the growing peak) in about the second half of the year. For the ascending part of the NDVI, the SOS date was extracted as the local maximum in the first derivative of the fitted NDVI, and for the latter part, the EOS date was extracted as the local minimum in the first derivative of the fitted NDVI. In other words, the date when the NDVI increased (or decreased) faster was determined as the SOS (or EOS).  
In Method 2, the NDVI data were fitted with a double logistic function, and the second-order derivative of the fitted curve was then calculated. The two dates corresponding to the two local maximum points in the first half of the year were the SOS and the onset of maturity. The two dates corresponding to the two local maximum points in the second half of the year were the onset of senescence and the EOS [24].  
In Methods 3 and 4, we also used the different threshold values method to extract the phenology. First, the NDVI data were fitted by the double logistic function shown in Equation (1). The fitted NDVI was then normalized by using the following function:  
\[\mathrm{NDVI}_{ntio}=\frac{\left(\mathrm{NDVI}-\mathrm{NDVI}_{min}\right)}{ \left(\mathrm{NDVI}_{max}-\mathrm{NDVI}_{min}\right)} \tag{2}\]  
where NDVI represents the fitted daily NDVI data by the double logistic function; \(\mathrm{NDVI}_{min}\) represents the minimum NDVI for each year, and \(\mathrm{NDVI}_{max}\) represents the maximum NDVI for each year. In Method 3, the SOS was determined as the day when the \(\mathrm{NDVI}_{ratio}\) increased to 0.2 in the spring, and the EOS was determined as the day when the \(\mathrm{NDVI}_{ratio}\) decreased to 0.2 in the autumn. Method 4 was similar to Method 3, except that the threshold ratio value was set as 0.5. Finally, we used the average value of the four SOS and EOS extraction results as the final SOS and EOS applied in our study to reduce uncertainties caused by any single calculation algorithm.