*2.2. Datasets*

NDVI is the most commonly applied vegetation index to characterize vegetation greenness and is strongly correlated with vegetation photosynthetic activity [32]. Climate change is the main factor affecting the change of vegetation greenness, and this change can be reflected by the spectral information of NDVI images. In this study, we used the NDVI datasets generated from NOAA/AVHRR series satellite images by the NASA MODIS13A2 group (Table 1). We used this dataset to extract the LSP dates for QMs from 2001 to 2019. We also excluded areas of bare soil/sparse vegetation with an annual average NDVI of less than 0.1 [33].


**Table 1.** Datasets and sources in the study area.

Meteorological data, including daily mean air temperature, daily total precipitation, and daily mean shortwave radiation and soil data, including daily mean soil temperature and moisture in 0–100 cm soil layer, from 2001 to 2019 were used in this study. These gridded data were derived from the ERA5-Land hourly data (Table 1). Moreover, we transformed the hourly climate data (24 hourly data were averaged for temperature, solar radiation, soil temperature and humidity, 24 hourly data were summed for precipitation) to daily-scale temporal resolution and resampled meteorological data to the same resolution as MODIS13A2 data. A time lag of 30 days before SOS and 60 days before EOS is defined as preseason.

The 300 m spatial resolution Climate Change Initiative Land Cover (CCI-LC) maps from 2001 to 2019 were available from the European Space Agency (ESA) (Table 1). CCI-LC discriminates 22 classes of land cover. In this study, we resampled these maps to 1 km and analyzed only pixels of unchanged vegetation types containing evergreen needle leaved forest (ENF), evergreen broadleaved forest (EBF), deciduous broadleaved forest (DBF), mixed forest (MF), shrubland (SL), grassland (GL), and cropland (CL).

#### *2.3. Retrieval of Phenology Metrics from NDVI Time Series Data*

The premise of quantitatively analyzing the phenology changes is to derive several key phenology metrics: SOS, EOS, LOS, MN, and MD (Figure 2). In this study, we firstly stacked the NDVI images from 2001 to 2019 in chronological order and smoothed the NDVI time series with a Savitzky-Golay (SG) filter for each pixel per year. The SG filter was chosen because it can best preserve the temporal vegetation dynamics and minimize atmospheric contamination and has also been integrated into the processing of the MODIS phenology product [34]. The smoothed data was used further for extracting phenology metrics of different vegetation types by detecting the inflection point (i.e., date) when the NDVI time series begins to ascend or descend for the specific year. This is the derivative method which the phenology metrics were extracted for each pixel per year, whereby the

maximum value of NDVI*ratio* corresponds to the greatest change of the smoothed NDVI time series [35]. Equation (1) is given as

$$\text{NDVI}\_{ratio(t)} = \frac{\text{NDVI}\_{t+1} - \text{NDVI}\_{t}}{\text{NDVI}\_{t}} \tag{1}$$

where NDVI*ratio*(*t*) is the calculated relative changing rate of NDVI at time *t* and NDVI*t* is the NDVI value at time *t*. Occurrence dates were obtained using these smoothed NDVI time series. SOS and EOS dates were determined as the day with the maximum and minimum NDVI*ratio*. LOS was determined to be the difference between EOS and SOS. MN was defined as the peak of vegetation growth, i.e., the NDVI value corresponding to NDVI*ratio* closest to zero. MD date was the middle date between the EOS date and the SOS date. The description of phenology metrics correlations is shown in Figure 2.

**Figure 2.** The description of phenology metrics correlations extracted using the NDVI time series datasets.

#### *2.4. Method and Statistical Analysis*
