*2.1. Study Area*

The QLMA lies at the intersection of the Tibetan, Mongolian and Loess plateaus (35.84◦–39.97◦N, 93.61◦–103.90◦E), with Qinghai Province in the south and Gansu Province of China in the north (Figure 2). The average elevation of the QLMA is over 3000 m, higher in the center and lower in the surroundings. This area belongs to the mid-latitude high elevation region. Most of the QLMA is in the temperate semiarid zone of the highlands [43], where solar radiation is strong. The average annual precipitation is 300–500 mm, more in the east than in the west. These make QLMA a sensitive area for climate change.

As an important ecological barrier in northwestern China, the QLMA is home to headwaters of rivers fed by snow and glacial meltwater. Grassland vegetation, desert vegetation, shrub and alpine vegetation cover more than 90% of the QLMA. Additional vegetation types include coniferous forest, broadleaved forest and swamp vegetation [44,45]. Because of its rich gene pool of alpine species, the QLMA is considered a key area and a priority area for biodiversity conservation in China [46]. A pilot Qilian Mountain National Park in the QLMA was established in 2017 to protect forests, grasslands, wetlands, surface water resources and glaciers [47]. Collectively, the QLMA is an ideal laboratory for studying climate change and ecosystems in cold and arid regions.

**Figure 2.** Map of the (**a**) land cover types, (**b**) elevation, (**c**) annual mean temperature, and (**d**) annual mean precipitation of the study area.

#### *2.2. Data Sources and Pre-Processing*

The normalized difference snow index (NDSI) was used to create snow pixels with a resolution of 500 m provided by MO(Y)D10A1 V6. The Terra(O) and Aqua(Y) satellites provide a separate daily NDSI from 1 September 2002 to 1 June 2020 [48]. 'NDSI\_Snow\_Cover\_Class' and 'NDSI\_Snow\_Cover\_Algorithm\_Flags\_QA' are the bands we used to eliminate invalid pixels and 'NDSI\_Snow\_Cover' is the band to extract NDSI. Data were from NASA's NSIDC DAAC at CIRES, with a large percentage of cloud pixels [49–51]. Some preprocessing was required to eliminate cloud obscuration [52,53]. First, pixels with a non-zero value of the 'NDSI\_Snow\_Cover\_Class' and 'NDSI\_Snow\_Cover\_Algorithm\_Flags\_QA' band were masked out, as these pixels represent invalid values such as missing data and clouds. Second, daily MOD10A1 and MYD10A1 data were combined using daily maximum values to reduce the number of cloud pixels. Finally, the max of the cloud-free pixels within a three-day time window surrounding the cloud pixels was taken as a replacement value. Experiments proved that these operations can effectively reduce the proportion of cloud pixels and improve the accuracy of snow products [51,54]. These preprocesses were completed using the Google Earth Engine.

The MOD09A1 V6 products provide the surface spectral reflectance of Terra MODIS bands 1–7 at 500 m resolution [55] containing seven bands that have been corrected for atmospheric conditions such as gasses and aerosols. It was used to calculate NDPI and thus estimate land surface phenology from 2003 to 2021. Data were provided by NASA LP DAAC at the USGS EROS Center.

The QLMA vegetation distribution dataset (vegetation pattern data (1:1,000,000) in the Qilian Mountains) was obtained from the National Cryosphere Desert Data Center (NCDC, http://www.ncdc.ac.cn/, accessed on 6 September 2021) and the digital elevation model (DEM) from 30 m ASTER GDEM was downloaded from the Geospatial Data Cloud (http://www.gscloud.cn/, accessed on 13 January 2022). We used these data to distinguish the six types of vegetation and every 100 m in elevation of QLMA. Figure 3 shows the technical workflows.

**Figure 3.** Overview of technical workflows.

#### *2.3. Calculation of Snow Cover Seasonality*

Value equal to 40 in the preprocessed 'NDSI\_Snow\_Cover' layer was chosen as a threshold for identifying snow pixels [56]. Pixels recognized as snow pixels were reclassified as 1, otherwise 0. Previous studies have shown that the snowpack on the Tibetan Plateau increases sharply in September and decreases sharply in the beginning of May [57]. Therefore, 1st September to 1st June of the following year is considered to be the potential snow season to exclude outliers that appear in the summer. A new time series was constructed: 1st September was set to be the first day of year, and 31st May of the following year was set to be the last. The metrics of the snow season were calculated based on methods in [58], where *FSD* and *LSD* are measured in day of year (DOY) and *SSL* and *SCD* are measured in days. The formulas are as follows:

$$FSD = SDOY\_{\text{min}} \tag{1}$$

$$LSD = SDO\mathcal{Y}\_{\text{max}} \tag{2}$$

$$SSL = LSD - FSD + 1,\tag{3}$$

$$\text{SCD} = \sum\_{i=1}^{n} S\_{i\prime} \tag{4}$$

where *SDOY* represents the distance between the date when a pixel is identified as a snow pixel and previous 1st September, *n* denotes the total number of days in the potential snow season, and *Si* denotes the state of snow cover for any given day within *SDOY*; *Si* equals 1 if there is snow and 0 if not.

#### *2.4. Calculation of Land Surface Phenology*

Shortwave infrared reflectance was combined with near-infrared and red reflectance to calculate the *NDPI*. It is an index for extracting accurate surface phenology to achieve high contrast between vegetation and background [59]. It was proven to be an accurate vegetation index for estimating land surface phenology in QLMA [42]. The formula is as follows:

$$NDPI = \frac{\rho NIR - \left[a \times \rho RED + (1 - a) \times \rho SWIR\right]}{\rho NIR + \left[a \times \rho RED + (1 - a) \times \rho SWIR\right]},\tag{5}$$

where *ρNIR* represents the near-infrared band, *ρRED* represents the red band and *ρSW IR* represents the shortwave infrared band. In MOD09A1, *ρRED* corresponds to band 1, *ρNIR* to band 2 and *ρSW IR* to band 6. *α* is a constant value for a given sensor and is taken as 0.74 for MODIS products [59].

In this study, we first masked the pixels of cloud shadow, snow and cloud, then calculated *NDPI* values for the filtered images. Then we reconstructed the time series

to unify the number of images into 72 in each year by calculating the average value. Finally, we smoothed the curves using the Savitzky–Golay (SG) filter. The window size was set to 5 and the number of polynomials was set to 3. Several methods have been developed to detect land surface phenology based on vegetation indices [60–63]. A simple and effective dynamic threshold method was used to extract land surface phenology from the reconstructed *NDPI* time series (Figure 4) [64]. Pixels with an intra-year *NDPI* change of less than 0.1 were excluded, and these were assumed to be areas without significant seasonality, as the *NDPI* of vegetation with distinct phenological stages generally increases from very small values (close to 0) to above 0.4 in study area. The formula is as follows:

$$thd(t) = \frac{NDPI(t) - NDPI\_{\text{min}}}{NDPI\_{\text{max}} - NDPI\_{\text{min}}},\tag{6}$$

where *NDPI*(*t*) represents the *NDPI* value at the calendar year date sequence *t*; *NDPImin* represents the maximum value of *NDPI* in a year; *NDPImin* represents the minimum value of the time series vegetation curve on the left and right part of the curve (SOS corresponds to the left half of the curve and EOS corresponds to the right) in a year bounded by *NDPImax*; and *thd*(*t*) represents the percentage corresponding to *NDPI*(*t*) after stretching it in time to a range of 0–1. Critical thresholds of 30% and 70% were selected [34].

**Figure 4.** Schematic illustration for extracting LSP metrics using dynamics threshold method. The black curve represents the interannual NDPI for a given pixel. Red lines represent SOS, blue lines represent EOS and green lines represent LOS.
