How Does Spring Phenology Respond to Climate Change in Ecologically Fragile Grassland? A Case Study from the Northeast Qinghai-Tibet Plateau

Abstract

Vegetation phenology is an important indicator of global climate change, and the response of grassland phenology to climate change is particularly sensitive in ecologically fragile areas. To enhance the ecological security of the Tibetan Plateau, it is crucial to determine the relationship between fluctuations in the start of the growing season (SOS) and the response to environmental factors. We investigated the trends of the intra-annual (ten-day) and interannual spatiotemporal dynamics of the SOS on the Northeast Qinghai-Tibet Plateau (NQTP) from 2000–2020 with MOD09GA data. We identified the response relationships with environmental factors (climate, terrain) using the maximum value composite method and the Savitzky–Golay filtering and dynamic threshold method. The SOS was concentrated from the 110th to 150th days; the average annual SOS was on the 128th day, with a spatial pattern of “early in the east and late in the west”. The overall trend of the SOS was advanced (45.48%); the regions with the advanced trend were mainly distributed in the eastern part of the NQTP. The regions with a delayed SOS were mainly concentrated in the higher-altitude regions in the southwest (38.31%). The temperature, precipitation and SOS exhibited a reverse fluctuation trend around the midpoint of 2010. Precipitation affected the SOS earlier than temperature. When temperature became a limitation of the SOS, precipitation had a more significant regulatory effect on the SOS. The SOS and aspect, slope and altitude were distributed in axisymmetric, pyramidal and inverted pyramidal shapes, respectively. The SOS on shaded slopes was earlier and more intensive than that on sunny slopes. With increasing slope, the area of the SOS decreased, and it occurred later. The SOS area was largest at 4500–5000 m and decreased at lower and higher altitude intervals. The SOS occurred later as altitude increased.

1. Introduction

Human activities are the main catalyst for current global warming [1]. The massive emissions of greenhouse gases since the industrial revolution have led to varying degrees of land surface and ocean surface warming, with the rate of warming in the last 20 decades being twice that of the 20th century [2]. Vegetation ecosystems play an extremely important role as a bridge for material and energy exchange between the land and the atmosphere [3]; however, they are significantly changing due to the warming climate [4]. It has been shown that climate warming advances spring phenology in forests and grasslands in cold-temperate regions [5]. Vegetation phenology refers to the cyclical growth of vegetation over the course of seasons [6], where vegetation responds to temperature-driven climate change by regulating the length of the growing season. Spring phenology is the beginning of vegetation growth in spring, also known as the start of the growing season (SOS), which is a key variable in determining the length of the growing season. Additionally, the SOS is particularly sensitive to climate change and is regarded as the best “indicator” of global warming [7]. The advancement or postponement of SOS will indirectly change the productivity and carbon cycle of ecosystems, which in turn will cause changes in the regional climate system [8]. Therefore, it is important to clarify the relationship between current climate change and SOS and its impact on vegetation ecosystems.

Warming already has a significant impact on the SOS in the Northern Hemisphere [9], with many of the global mid- and high-latitude regions experiencing an earlier SOS due to warmer temperatures [10]. From 1982 to 2011, the SOS in the United States, China and Europe advanced by an average of 0.9, 4.2 and 4.7 days per decade, respectively [9], and the trend of advancement was more pronounced on the Qinghai-Tibet Plateau (QTP) than in other regions at the same latitude. The rapid climate change on the QTP has led to obvious spatial variability in the distribution of SOS, which is mainly manifested by the obvious delay of SOS in the high-altitude areas in the southwest and the trend of advancement in most other regions, while the interannual variation of SOS also shows an increasingly complex trend with increasing altitude, and the region has become a research hotspot for scholars in China and abroad. Zhang et al. [11] found that the SOS on the Tibetan Plateau continued to advance at a rate of 1.04 days per year from 1982 to 2011 and that the advancement trend was strongly associated mainly with rising winter and spring temperatures. Based on Systeme Probatoire d’Observation de la Terre (SPOT)/normalized difference vegetation index (NDVI) data, Ding et al. [12] found that the spatial and temporal distribution of alpine grassland phenology from 1999 to 2009 was closely related to hydrothermal conditions, and the SOS was delayed with the gradual deterioration of hydrothermal conditions. Wang et al. [13] used Global Inventory Modeling and Mapping Studies (GIMMS)/NDVI and SPOT/NDVI data to analyze vegetation phenology trends on the eastern edge of the Tibetan Plateau and found that the SOS in the region showed an early trend from 1982 to 1998, while a delayed trend was observed from 1999 to 2018. In a small-scale regional study, An et al. [14], based on Landsat 7 ETM+ (Enhanced Thematic Mapper Plus) and Landsat 8 OLI (Operational Land Imager) data from 2013 to 2015 along a 2600 km transect in the center of the Tibetan Plateau, investigated how vegetation changes in response to differences in slope aspect and elevation. The results showed that a 100 m increase in elevation delayed the SOS by 1.52 days and that the variation in slope temperature due to solar radiation was the main factor for the difference in grass SOS between different slopes. Li et al. [15] used a phenology camera to study the vegetation phenology of four alpine grasslands with different grazing intensities from 2015 to 2017, and the results revealed differences in the SOS of alpine grasslands with different grazing levels. The QTP is a geographic unit with a complex structure and presents diverse phenological patterns. In previous studies, GIMMS/NDVI products (15 days, 8 km) and STOP/NDVI data were widely used for SOS surveys on the QTP [16], and these data have been advantageous in large-scale studies; however, remote sensing data with low spatial and temporal resolutions are not well representative of regional SOS changes. Some studies have shown that GIMMS data from 2001 to 2006 have serious data quality problems in the Tibetan Plateau region, while Moderate Resolution Imaging Spectroradiometer (MODIS) data can maintain better data quality [11]. On the other hand, although Landsat data and phenology cameras can provide remote sensing data with high resolution, they still cannot be applied to long-term and large-scale phenology monitoring because of their small monitoring area and poor spatial and temporal continuity. The response of upland grassland SOS to climate change is extremely sensitive, so remote sensing data with high temporal and spatial resolution are more advantageous for studying SOS trends.

Qinghai Province is located in the hinterland of the Tibetan Plateau, which is an important ecological barrier of the “Third Pole”. Rich grassland resources play an important role in the development of animal husbandry, biodiversity and prevention of soil erosion. However, grassland ecosystems have undergone significant changes due to the unique natural environment and the drastic disturbance of human activities. In recent years, many ecological restoration projects carried out by the Chinese government on the Tibetan Plateau have achieved good results, and relevant data show a close link between the restoration of grassland ecosystems and the advancement of SOS caused by climate change [17]. However, the focus of many studies is on the Tibetan Plateau, but the Tibetan Plateau is large in area, with considerable topographic variability, diverse vegetation types, strong spatial heterogeneity and SOS-sensitive response relationships to climate change; thus, studies of the entire Tibetan Plateau are not representative of the weather-driven patterns in Qinghai Province [18]. Therefore, in this study, based on MOD09GA satellite images from 2000 to 2020 and data on climate (temperature, precipitation) and topography (elevation, slope and slope aspect) for the same period, the dynamic threshold method proposed by Jönsson and Eklundh [19] was used to extract vegetation SOS based on the TIMESAT 3.0 platform, and the SOS of Qinghai Province was extracted and analyzed at intra-annual (ten-day) and interannual scales. The analysis was conducted to explore the dynamic fluctuations and driving relationships of SOS over the past 21 years to clarify the SOS regulation mechanism of grasslands in Qinghai Province. The results provide a theoretical basis for the targeted formulation of grassland resource conservation policies and ecological research.

2. Materials and Methods

2.1. Study Area

Qinghai Province is an important part of the QTP, with a total area of 69.66 × 10⁴ km², making it the fourth-largest province in China, and it is located in the interior of Northwest China (31°9′–39°19′ N, 89°35′–103°04′ E). The province has the largest saltwater lake, “Qinghai Lake”, and is also one of the densest areas of highland lakes in China, where the Yangtze, Yellow and Lancang Rivers originate; thus, the region is known as the “Water Tower of Asia”. There are 2 cities and 6 states under its jurisdiction, and its capital is Xining (Figure 1). The average altitude of the province is 3000 m, with a high altitude in the west and a low altitude in the east, and the high mountains crisscross the province and are widely distributed. The Qilian Mountains-Altun Mountain System and the East Kunlun Mountain System constitute an important physical geographic boundary, dividing Qinghai Province into the Eastern Agricultural Region (EAR; the area of cultivated land in this area accounts for 0.60% of the total area of the province, which is relatively small, and the crop growing season is fixed, so the impact of farming activities on the dynamic changes in the SOS can be ignored), the region around Qinghai Lake and Qilian Mountain (ALQM), the Qaidam Basin (QB), and the Southern Qinghai Plateau (SQP). Qinghai Province has a continental arid and semiarid plateau climate, which is extremely unevenly distributed and affected by the terrain. The summer is short, and the annual temperature in the whole province is generally low, with an average annual temperature of −5.6 to 9.0 °C. The average annual temperatures in the EAR and QB are 3 to 9 °C and 2 to 5 °C, respectively, which represent warm areas and subwarm areas in the province, respectively. The SQP and ALQM, with higher altitudes, have an average annual temperature of −4 to −6 °C or less, representing the cold regions of the province. There are significant regional differences in precipitation distribution, and the average annual precipitation is 17.6–764.4 mm. The precipitation in the southeastern part of the province is high, and it decreases from southeast to northwest. For example, the annual precipitation in the QB in the northwest is only 25–50 mm, making it one of the driest areas in China.

Qinghai Province is rich in grassland resources, with a natural grassland area of 4.19 × 10⁵ km². It is one of the five largest pastoral areas in China, and animal husbandry is one of the pillar industries in Qinghai Province. The grassland growing season is usually from early May to late October. A total of 3.86 × 10⁵ km² of natural grassland can be grazed and produced, and this grassland is mainly distributed in the southeast of the SQP, ALQM and QB. Alpine grassland, alpine meadow and temperate grassland are the main grassland types in Qinghai Province. In the process of experiencing global warming, the overall abundance of grassland has increased, but to some extent, the biomass of Cyperaceae and miscellaneous grasses has decreased.

2.2. Data and Preprocessing

2.2.1. Remote Sensing Data

We used the National Aeronautics and Space Administration (NASA, https://search.earthdata.nasa.gov/search (accessed on 14 March 2021)) daily surface reflectance product MOD09GA, which has a 500 m resolution and rigorous radiation correction, atmospheric correction and geometric correction. The time series was from 24 February 2000 to 31 December 2020. The tiles were H25V05 and H26V05, and there were a total of 15,080 scenes during the 21 years. The red and near-infrared bands were extracted after image mosaicking, projection and format conversion using MODIS Reprojection Tools (MRT) software.

2.2.2. Meteorological Data

The China surface climate dataset (V3.0) was selected from the China Meteorological Data Network (https://data.cma.cn/ (accessed on 3 April 2021)), and temperature and precipitation data from 109 meteorological stations in and around Qinghai Province (Figure 1B) were extracted from 10 April to 20 May of each year from 2000 to 2020 according to the boundary of the study area. The average value of precipitation accumulation in the preseason (1 March to 10 April) over the last 20 years was 11.36 mm. Therefore, the effect of preseason precipitation on the SOS was not considered. Based on ANUSPLINE software, the meteorological data were interpolated using the local thin disk smooth spline method with digital elevation model (DEM) data as covariates, and finally, raster meteorological data with a spatial resolution of 500 m were obtained. We read the interpolation results and defined the projection as WGS84 in ArcGIS 10.5 software. After the mask, we obtained the ten-day average temperature and cumulative precipitation raster data of Qinghai Province from mid-April to mid-May of 2000–2020 (the first ten days, middle ten days and last ten days of the month corresponded to the 1st to 10th, 11th to 20th, and remaining days of the month, respectively).

2.2.3. Terrain Data

We selected 3 topographic factors—altitude, slope and aspect—graded with the actual conditions of the study area (

Table 1. Altitude, slope and aspect direction grading table.

Category	I	II	III	IV	V	VI	VII	VIII
Altitude (m)	<3000	3000–3500	3500–4000	4000–4500	4500–5000	>5000	Null	Null
Slope (°)	0–2	2–5	5–8	8–12	12–16	16–20	20–35	35–47
Aspect	North	Northeast	East	Southeast	South	Southwest	West	Northwest

) and then extracted data (from the Resources and Environment Science Center of the Chinese Academy of Sciences, http://www.resdc.cn (accessed on 26 April 2021), the resolution was 500 m), which were based on the Chinese DEM.

2.3. Methods

2.3.1. Vegetation Information Inversion

The normalized difference vegetation index (NDVI) was used to characterize the vegetation information and was calculated as follows:

$$\mathrm{NDVI}=\frac{{\rho }_{\mathit{nir}}-{\rho }_{\mathit{red}}}{{\rho }_{\mathit{nir}}+{\rho }_{\mathit{red}}}$$

(1)

where ρ_nir and ρ_red are the near-infrared and red wavelengths, respectively, and the NDVI value range is [−1, 1].

2.3.2. Time Series Reconstruction

We used the maximum value composite (MVC) pixel by pixel to calculate the ten-day NDVI using the TIMESAT 3.3 software package to smooth the NDVI based on MATLAB 2019 to eliminate abnormal fluctuations. The Fourier transform method [20], harmonic analysis method [21] and Savitzky–Golay (S-G) filter method [22] are commonly used methods for fitting the reconstruction of time series; among them, the S-G filter method can effectively eliminate the outliers that deviate from the vegetation growth trajectory in the smoothing process and has a better fitting effect in the QTP region [23,24]. The NDVI time-series data after S-G filtering (Equation (2)) effectively reduced the anomalous noise and more realistically reflected the vegetation growth in Qinghai Province (Figure 2; the green dots are the SOS for each year).

$$Y_j=\frac{{\displaystyle \sum _{i=-m}^{i=m}{C}_i{Y}_{j+i}}}{N}$$

(2)

where Y_j is the time series data at the end of the fit, Y_j+i is the original data, C_i is the filter coefficient and N is the moving window size (2m + 1).

2.3.3. Extraction of Vegetation Phenology Information

The main methods of extracting vegetation phenology information are the dynamic threshold method [25], maximum slope method [26] and sliding average method [27]. Combined with the actual condition of the study area, the dynamic threshold method is more suitable to extract SOS information from areas with large internal differences [28]. Based on the suggestion of Jönsson and Eklundh [19] and the research experience of Guan et al. [29] on the QTP, after several experiments, we set the threshold value of the SOS to 0.2 (that is, the NDVI of the growing season reaches 20% of the annual amplitude) and the sliding window size to 4. In addition, we excluded the areas where the NDVI < 0.05 to reduce the disturbance of the unvegetated area in the central part of the QB [30]. Since the damage to the original remote sensing image in 2003 caused the failed extraction of vegetation regreening information, the average value of the regreening period from 2000 to 2006 was used instead of the value of the regreening period in 2003 to reduce the influence of the data anomaly in that year on the experimental trend (Day 132, Figure 2).

2.3.4. Analysis of Spatial and Temporal Trends in the SOS

Theil–Sen median trend analysis (Equation (3)) is an important method used to evaluate the trend of long time series data [31,32] and is widely used in the field of vegetation growth [33]. The Mann–Kendall test [34,35] (Equation (4)) is a nonparametric test method that does not require the data to obey certain distribution rules to minimize the impact of data errors [36].

$$\beta =Median\left(\frac{X_j-X_i}{j-i}\right)$$

(3)

where β is the vegetation growth trend and X_i and X_j denote the values of the i-th and j-th terms in the time series (2000 ≤ i < j ≤ 2020), respectively. When β > 0, the SOS has a delayed trend, and when β < 0, the SOS has a decreasing trend.

$$S={\displaystyle {\sum }_{i=1}^{n-1}{\displaystyle {\sum }_{j=i+1}^{n}sgn}}\left(x_j-x_i\right)$$

(4)

$$sgn\left(x_j-x_i\right)=\{\begin{array}{l}+1, x_j-x_i>0\\ 0, x_j-x_i=0\\ -1, x_j-x_i<0\end{array}$$

(5)

where S is the test statistic; X_j and X_i are the data values corresponding to the j-th and i-th moments, respectively; n is the total number of time series; and sgn is the sign of the function. When n > 10, the S statistic can be calculated from the Z-values as follows:

$$Z_s=\{\begin{array}{l}\frac{s-1}{\sqrt{var(s)}},ifs>0\\ 0, ifs=0\\ \frac{s-1}{\sqrt{Var(s)}},ifs<0\end{array}$$

(6)

Among them:

$$Var(s)=\frac{n(n-1)(2n+5)-{\displaystyle {\sum }_{i=1}^m{t}_i(t_i-1)(2t_i+5)}}{18}$$

(7)

where Z_s is the standardized statistic, V_ar(s) is the variance, m denotes the number of recurring sets in the series and t_i is the number of repeated data in the i-th repeated data set. A bilateral test is used at a given significance level α. |Z_s| > Z_1-α/2 indicates a significant trend of data change in the time series.

2.3.5. SOS Change Trend Statistics

We combined the Theil–Sen median (Sen_sos) trend analysis with the Mann–Kendall (MK) test to reflect the trend of the SOS in Qinghai Province from 2000 to 2020. According to the actual situation of the study area, we took the standards of −0.1 and 0.1 to divide the Sensos. The results were divided into significant (|Z value| ≥ 1.96) and nonsignificant (|Z value| ≤ 1.96) when the confidence level of the MK significance test was set to 95% and then superimposed with the results of the Sensos trend analysis. The results were divided into five change types (

Table 2. Statistics of the SOS trend.

Sensos	|Z Value|	Trend of SOS
≥0.1	≥1.96	Significant delaying trend
≥0.1	≤1.96	Slight delaying trend
−0.1 ≤ Sensos ≤ 0.1	≤1.96	Stability
≤−0.1	≤1.96	Slight advancing trend
≤−0.1	≥1.96	Significant advancing trend

) to analyze the trend of vegetation SOS in Qinghai Province.

2.3.6. Correlation Analysis of SOS and Environmental Factors

We used correlation analysis to reveal the relationship between vegetation SOS and temperature and precipitation. When there is a relationship between these variables, the partial correlation coefficient (Equation (9)) can indicate the closeness of the relationship; t-test was used to test for significance (Equation (10)).

(a)

Correlation analysis

$$R_{xy}=\frac{{\displaystyle \sum _{i=1}^n\left[\left(x_i-x_p\right)\left(y_i-y_p\right)\right]}}{\sqrt{{\displaystyle \sum _{i=1}^n{\left(x_i-x_p\right)}^2{\displaystyle \sum _{i=1}^n{\left(y_i-y_p\right)}^2}}}}$$

(8)

where R_xy is the correlation coefficient of variables x and y, the value range is [−1, 1], and the closer the absolute value is to 1, the higher the correlation is; x_i and y_i denote the values of variables x and y in year i, respectively; x_p and y_p denote the multiyear averages of x and y, respectively; and n is the number of sample years, 1 ≤ n ≤ 21.

• (b) Partial correlation analysis

$$r_{xy·\mathrm{z}}=\frac{r_{xy}-r_{xz}{r}_{yz}}{\sqrt{\left(1-r_{xz}{}^2\right)\left(1-r_{yz}{}^2\right)}}$$

(9)

• (c) t-test

$$t=\frac{r}{\sqrt{1-r_{xy\cdot z}{}^2}}\times \sqrt{n-m-1}$$

(10)

where r_xy·z is the partial correlation coefficient between x and y assuming a constant variable z, and so on; m is the number of independent variables. Due to the large sample size of the data (2,499,626 pixels) and the fact that the data used (MODIS satellite images and rasterized meteorological data) were not real data, p < 0.1 was taken for the significance test.

3. Results

3.1. Factors Influencing the SOS Fluctuation

3.1.1. Climatic Factors

The results of the precipitation and SOS bias correlation showed (Figure 3A–E) that there were maximum positive and negative correlations in late April and early May, respectively. The positive correlation between the precipitation and SOS in late April was as high as 77.09%, which indicates that precipitation plays a positive role in the SOS in late April, while the central-eastern part of the EAR showed a significant negative correlation, and the rest of the regions, especially the southwestern part of the ALQM and the central-western part of the SQP with higher elevation, showed a significant positive correlation. In early May, the area of negative correlation was 75.18%, mainly concentrated in the central and southeastern regions of Qinghai Province. The significance test (Figure 4A–D) is similar to the partial correlation results, with the largest significant (p < 0.1) area in late April, when the number of pixels was 95.53%, 42.49% and 70.05% higher than in mid-April, early May and mid-May, respectively. The spatial distribution of the significant regions (p < 0.05) was generally consistent with the positively correlated regions in late April.

The time of the positive and negative bias correlation between temperature and SOS was exactly opposite that of precipitation (Figure 3F–J). In late April, 51.15% of the areas showed a negative bias correlation with temperature, and these areas were mainly distributed at the junction of the ALQM, QB and SQP, while most of the other regions showed a slight bias correlation with the temperature. In early May, 82.93% of the regions had a positive bias correlation, and the significantly biased correlated regions were mainly located in the northwestern part of the ALQM, the southern part of the QB and the northwestern part of the SQP, which are all high-altitude regions. The regions with a slightly biased correlation were mainly distributed in the northeastern part of the EAR and SQP, and the regions with a negatively biased correlation accounted for less than 20%; these areas were mainly distributed in the central part of the SQP and the northeastern corner of the QB. Interestingly, the significance (Figure 4E–H) is also somewhat consistent with the correlation in terms of temporal and spatial distribution, i.e., the area of significance is also the largest in early May (when the number of pixels was 49.45%, 114.41% and 97.29% higher than in mid-April, late April and mid-May, respectively) and somewhat consistent with the area of positive correlation in terms of spatial distribution. Overall, temperature and precipitation were the decisive factors affecting vegetation SOS in Qinghai Province; however, precipitation affected vegetation SOS approximately earlier than temperatures.

According to the bias correlation results between vegetation SOS and climate factors, the most significant influences of temperature and precipitation on SOS in Qinghai Province were in late April and early May, respectively. Therefore, based on the results from previous studies, we selected the temperature, precipitation, and annual average SOS data from late April to early May from 2000 to 2020 for time series analysis.

A time series analysis of the temperature, precipitation, and SOS for late April and early May over a 21-year period was performed (Figure 5), and we found that there were large fluctuations in temperature, precipitation and SOS during the study period. The years with maximum and minimum precipitation were 2001 (10.043 mm) and 2003 (1.404 mm), respectively, and the temperature was highest in 2004 (14.84 °C) and lowest in 2020 (−2.63 °C). Precipitation showed an increasing trend of 0.66 mm/year, while temperature showed a decreasing trend of 0.21 °C/year. Although the SOS experienced a large delay trend from 2015 to 2019, overall, the SOS experienced an advancing trend at a rate of 0.28 days/year.

Through further analysis, we found that the temperature, precipitation and SOS were all centered around 2010 and showed a reverse trend of fluctuation before and after. The temperature showed an increasing trend of 0.22 °C/year before 2010 but began to decline sharply with a cooling trend of 0.56 °C/year after 2010; similarly, precipitation and SOS showed a decreasing/advancing trend before 2010 (0.337 mm/year, 1.19 days/year) and a trend of increase/delay (0.339 mm/year, 0.24 days/year) after 2010. From the consistency of the three fluctuation trends, the SOS and precipitation maintained a relatively synchronous fluctuation trend, the overall relationship was closer than that of temperature, and the closeness of this trend had obvious periodicity. From 2000 to 2009, the annual average temperature was 6.33 °C and showed an increasing trend each year, and the synchronous fluctuation trend of SOS and precipitation was not obvious. After 2010, the annual average temperature was 4.15 °C and showed a fluctuating and decreasing trend. In this case, the SOS and precipitation showed a more consistent fluctuation. The results showed that when the temperature was not the limiting factor of the SOS, the time of the SOS was regulated by both variables. When the temperature was low, the advance or delay of the SOS had a greater correlation with precipitation. When the precipitation increased, the SOS was delayed; when the precipitation decreased, the SOS was advanced; and when the temperature was lower, the regulation of precipitation on SOS time was more obvious.

The partial correlation results between climate factors and SOS indicated that both temperature and precipitation were determinants of SOS in Qinghai Province, and precipitation had an earlier and greater effect than temperature on SOS. The time series analysis of temperature, precipitation and SOS showed the fluctuating relationship between different factors, and a main conclusion was that SOS was regulated by both temperature and precipitation when the temperature was higher, while precipitation maintained a more synchronized relationship with SOS when the temperature was lower.

3.1.2. Terrain Factors

The terrain of Qinghai Province is high; there is considerable variability in altitude, and the altitude increases from southwest to northeast, with mountain peaks and extremely complex terrain. Different external agents lead to a differential distribution of vegetation cover, which also determines the spatial and temporal differences in the SOS at different slope aspects, slope gradients and altitudes.

From the 21-year average SOS data, it can be found that the distribution of SOS varies greatly in different slope directions (Figure 6A) and basically shows an axisymmetric distribution, with SOS on shady slopes (east- and north-facing slopes) occurring earlier than on sunny slopes (south- and west-facing slopes) (Figure 6A), which was basically axisymmetric. The SOS areas were mainly dominated by the north (17.27%) and south (14.13%) slopes, followed by the northeast (14.13%) and southwest (12.78%) slopes, and the east (9.97%), southeast (9.91), west (9.43%) and northwest (10.02%) slopes all accounted for approximately 10%. Among them, the north slopes reached the peak of the SOS on the 130th (4.23%) and 140th (4.42%) days, which was earlier than that of the south slopes; however, the peak SOS on south slopes was concentrated on the 140th day, accounting for the largest proportion of all slope aspects (4.44%) in the same period. Additionally, there was still 4.31% of the area in the SOS from the 150th to 180th days. On the northeast slopes, the first SOS peak (2.28%) appeared on the 120th day, which gradually increased to the maximum on the 130th day (3.50%) and 140th day (3.55%) and then decreased (2.23%) on the 150th day. The SOS on the southwest slopes were also slightly later than on the northeast slopes, and the first SOS peak appeared on the 130th day (2.94%). The area of the SOS increased to the maximum on the 140th day (3.61%) and then decreased rapidly (2.50%); the time of the SOS was slightly late and concentrated. The SOS on the east and southeast slopes was also slightly earlier than that on the west and northwest slopes, and both reached a peak on the 140th day. Overall, the east slopes and the north slopes of the shady slopes in Qinghai Province had an earlier SOS, reaching the first peak on the 120th day and the maximum on the 140th day. The sunny slopes were dominated by the south and west slopes, which had the first small SOS peak on the 130th day and the largest peak on the 140th day; furthermore, the area of SOS after the 150th day was still higher than that of the east and north slopes.

The overall relationship between the SOS and slope is “pyramidal” (Figure 6B), and the area of the SOS is mainly in the area with less slope. The percentage of SOS area from 0° to 2° was 32.11%, followed by that from 2° to 5° (26.55%) and that from 5° to 8° (14.40%), while that from 8° to 47° accounted for only 25.94%, and the SOS area decreased with increasing slope. Among them, there were three identical peaks of SOS from 0°–2° and 2°–5°; that is, the SOS increased on the 130th day, the maximum SOS increase rate occurred on day 140, and the SOS decreased rapidly after the 150th day. When the slope was >5°, the peak SOS was on the 140th day. Overall, the larger the slope was, the smaller the SOS area was and the more concentrated the SOS time was.

Contrary to the slope, the relationship between the SOS and altitude has an inverted pyramidal shape (Figure 6C), and the SOS area is mainly distributed in the high-altitude area due to the alpine climate. The SOS area of 4500–5000 m was the largest (34.07%), followed by that of 4000–4500 m (27.75%) and that of 3500–4000 m (14.46%), while the SOS areas of >5000 m and <3500 m accounted for only 23.72%. In the regions with altitudes higher than 5000 m, the SOS characteristics were not obvious. However, in the four intervals in which the altitude was lower than 4500 m, the peak SOS occurred on the 130th or 140th day, the proportion of the SOS on the 140th day gradually increased with increasing altitude, and the first peak period of the SOS on the 150th day appeared after it rose to 4500–5000 m, accounting for 7.18%. The altitude of Qinghai Province is high overall, with 22.93% and 29.11% of the area at 4000–4500 m and 4500–5000 m, respectively, and similarly, the area of SOS accounted for the largest proportion, followed by the area below 3000 m, with an area of 16.20%. The main area was located in the middle of the QB, which has basically no vegetation-covered area, with only 8.24% of the SOS area. However, the peak of the SOS was on the 130th day, slightly earlier than that of the other altitudinal gradients, likely because the EAR’s farming activities cause its SOS to be slightly advanced. In general, the SOS was positively related to altitude, and with increasing altitude, the SOS was gradually delayed, with the most obvious delay effect at 4500–5000 m.

The difference in SOS between different topographic factors mainly depends on the different hydrothermal conditions. The SOS on shaded slopes is earlier than that on sunny slopes because of more adequate moisture conditions and lower evapotranspiration due to less solar radiation. A steeper slope will cause more soil erosion, while an area with a more gradual slope is more suitable for vegetation growth and development, so the area of SOS is largest in the gradually sloping area. The average elevation of the study area is high, which is the main reason why the overall relationship between SOS and elevation is inverted pyramidal, and it can be clearly seen that the peak of SOS is gradually delayed from low to high elevation.

Both climatic and topographic factors influence the SOS, but climatic factors play a more decisive role. In highland areas, elevation, slope gradient and slope aspect are important factors in determining SOS, but with the continuous change in topographic factors, the water and heat conditions in the environment also change, which causes the differences in SOS among different elevations, slope gradients and slope aspects; thus, temperature and precipitation play a decisive role, and the difference in SOS caused by different topographic factors is a phenomenon that can be expressed by climate regulation.

3.2. SOS Time Series Fluctuations and Spatial Distribution

The SOS series in Qinghai Province from 2000 to 2020 showed a fluctuating trend of “delayed–advanced–delayed” but advanced overall (Figure 7). Except for the obvious delay of the SOS before 2004 and after 2015, there was no drastic fluctuation in other years. The average annual SOS time was the 128th day; the earliest and latest SOS times were the 117th day (2019) and the 142nd day (2016), respectively; and the maximum amplitude was 25 days. The SOS was concentrated from the 110th to 150th day (90.65%), that is, from mid-April to late May, and mid-May and early May were the two largest SOS peaks, accounting for 26.45% and 24.65%, respectively.

The bar chart in Figure 7 shows a gradual increase in the degree of fluctuation of SOS in later years. The percentages of SOS occurring on the 110th to 120th days in 2009 and 2019 were 10.22% and 7.59% higher, respectively, than the average percentage over 21 years, which indicates that the SOS was significantly earlier in those years. The percentage of SOS occurring on the 110th to 120th days increased year by year, while the percentage on the 140th to 150th days showed a decreasing trend, which indicates that SOS tended to be earlier year by year during the study period.

The spatial distribution of the average SOS from 2000 to 2020 was uneven (Figure 8A), with an overall pattern of “early in the east and late in the west”. In the EAR, QB and the southeastern and western parts of the SQP, the SOS was concentrated from the 110th to 130th days, making them regions with an earlier SOS. The region with late SOS extended from the northern part of the ALQM to the southwestern corner of the SQP in a “V”-shaped zonation. The SOS of the EAR was the earliest, and the area was the largest (43.53%), as was the case in most of the regions in the southeastern SQP (24.35%), which reached the first SOS peak on the 130th day (Figure 8B). In most areas of the ALQM (34.16%), QB (18.37%) and SQP (24.88%), the peak SOS was on the 140th day.

3.3. SOS Fluctuation Trend

From 2000 to 2020, the SOS of Qinghai Province showed an overall advancing trend (45.48%). A total of 41.91% and 3.56% of the areas showed slight and significant advancing trends, respectively (Figure 9B). The regions with large advancing trends were mainly distributed in the ALQM, EAR, and southeastern part of the SQP (Figure 9A), accounting for 62.11%, 55.32% and 42.96% of their respective regions. Among them, the regions with significant advancing trends were distributed in the Hehuang Valley, which is at the junction of ALQM and SQP, accounting for 3.56% of the total area of Qinghai Province.

In Qinghai Province, 16.21% of the area showed a stable trend, and these areas were widely distributed across the whole province. Among them, except for QB, which accounted for the largest proportion (22.73%), the other three regions accounted for approximately 14.00% of their respective regional areas. The regions with a delayed SOS accounted for 38.31% of the total area of the province, and the regions with a slight delay and a significant delay accounted for 35.83% and 2.49%, respectively. Among them, the regions with a slight delay were mainly distributed in the southwest of the SQP and the southeast of the QB, accounting for 39.45% and 35.36% of their respective regions, while ALQM and EAR accounted for less than 30%. The regions with a significant delay were mainly concentrated in the high-altitude regions of the QB and SQP, accounting for 2.49% of the total area of the province. Overall, the spatial variation in the SOS change trend was significant, and all regional SOS trends were dominated by slightly advancing trends, supplemented by slightly delayed trends. Spatially, the SOS was gradually delayed from the northeast to southwest, the regions with an advancing trend were mainly concentrated in the low-altitude areas in the southeast of Qinghai Province, and the regions with a delayed SOS were mainly in the high-altitude areas of the SQP.

4. Discussion

4.1. The Spatiotemporal Trends in the SOS

Although the SOS in Qinghai Province was delayed before 2004 and after 2016, the overall trend of the SOS from 2000–2020 was an advance of 0.28 days/year, which was similar to the conclusion of Chang et al. [37] on the QTP, who found a 0.25 day/year advance. In recent years, the trend of SOS advance on the QTP has been confirmed in several studies [38]. Similarly, Shen et al. [39] found that the SOS on the QTP was delayed in 1998–2003 and then advanced in 2003–2009, which was more consistent with the results of this study. Moreover, vegetation phenology is a complex ecological phenomenon that is a joint result of the interaction between the regional environment and vegetation. It has become a consensus that climate warming will lead to an earlier SOS [40], but the latest research shows that vegetation becomes less sensitive to warming driven by long-term climate warming [41]. This result suggests that a continuously warming environment does not lead to an early SOS to maintain an increasing trend because the stimulation of a low temperature is needed for vegetation rejuvenation (vernalization). However, the warming in winter and spring causes vegetation to take a longer time to reach a low temperature [42], which may be one of the reasons for the different SOS trends before and after 2016. In the study area, the SOS was mainly concentrated from the middle of April to late May, and the peak of the SOS was the earliest in the EAR (early May), while those in the rest of the regions were in mid-May. Guo et al. [43] conducted phenological observations in the Xilingol grasslands of northern China and found that the SOS was concentrated from early April to mid-May overall, which was 10 days earlier than the results of this study. The Xilingol league has a continental monsoon climate, and the different hydrothermal conditions in the study area were the main reason for the differences in the SOS.

The overall spatial distribution of the SOS in Qinghai Province showed a pattern of “early in the east and late in the west”, with the early region concentrated in the southeast and the late region distributed in a “V”-shaped belt from the north to the southwest corner; this pattern reflected the geomorphological pattern from the northwest by west to southeast by east in Qinghai Province. The northwestern Qilian Mountains–Altun Mountain System and the middle East Kunlun Mountain System divide Qinghai Province into three regions; specifically, the QB is tightly surrounded by the two mountain systems, precipitation in the basin is scarce, wind is strong, and vegetation growth conditions are relatively poor [44], delaying the SOS. From the ALQM to the EAR, the altitude gradually decreases, the hydrothermal conditions improve, and the SOS gradually advances. South of the East Kunlun Mountain system is the vast SQP, and the high-altitude region in the central and western parts of the plateau becomes a hindering condition for vegetation SOS [45]; thus, the SOS in this region was earlier in the southeast and later in the northwest.

In most of the alpine regions, increased temperature had a significant advancing effect on the SOS of vegetation [46], which was the same as the overall advancing trend of the SOS in Qinghai Province found in the present study; however, there was a large area with a delayed trend in the southwestern high-altitude region. This phenomenon could be explained by the research of Dorji [47], who, after exploring the effects of temperature and precipitation on the semiarid grassland of the QTP, found that the SOS of plants with shallow root systems would be delayed due to the drought caused by the increase in temperature. This finding indicated that the SOS of vegetation in the arid area was limited by the soil surface water to a certain extent. Thus, the soil moisture becomes a limiting factor for the SOS in the high-altitude region of southwestern Qinghai Province.

4.2. SOS Response to Climatic Factors

The response of vegetation SOS to temperature and precipitation varies among habitats [48]. In this study, the correlation between air temperature and precipitation and SOS was determined using a fixed time period (mid-April to mid-May). SOS was most positively correlated with precipitation and temperature in late April and early May, respectively, and precipitation affected vegetation SOS to a greater extent than temperature. The phenological observation experiment of Zhu et al. [49] on alpine meadow grasslands in Nagqu, Tibet, also showed that soil water stress caused by temperature warming would delay the SOS. It was found that the delay or advancement of most SOSs in alpine meadows was highly correlated with the concentrated precipitation before the SOS. Similarly, the results of Meng et al. [50] showed that the positive effect of climate warming on vegetation growth on the QTP decreased each year, leading to a more sensitive demand of vegetation for precipitation, which indicated that precipitation was a key factor affecting the SOS under arid and semiarid climatic conditions.

Interestingly, we found that the most positive correlation between precipitation and SOS occurred at exactly the same time as the maximum negative correlation between temperature and SOS, suggesting that increased preseason temperature hinders the positive effect of precipitation on the SOS. Shen et al. [51] found that in most regions of the QTP, especially in the southwest and northeast, increased preseason precipitation would promote an advanced SOS; however, increased preseason temperature delayed the SOS. However, in the regions with more adequate moisture, the response of the SOS to temperature was more sensitive, and an increase in precipitation did not advance the SOS; instead, a lack of light intensity or a shortening of daylight hours indirectly reduced surface temperature and thus delayed the SOS [52]. This result suggested that after precipitation provided the basic conditions for the SOS of grasslands in Qinghai Province, the increase in temperature would promote vegetation growth and development, while the precipitation at the SOS would lead to a decrease in surface temperature, thus inhibiting vegetation growth. Thus, the temperature became the dominant influencing factor at this time, which could also largely explain why precipitation and SOS showed more consistent fluctuations due to the decrease in preseason temperature after 2010.

In arid and semiarid areas, precipitation before the SOS can offset the delayed SOS phenomenon brought by warming, and warming and sufficient precipitation can better promote an advanced SOS [53]. In summary, both temperature and precipitation were determinants of the SOS of the vegetation in Qinghai Province, but precipitation played a more important role in regulating the SOS in Qinghai Province.

The bias correlation analysis shows that the correlations between temperature, precipitation and SOS are not high. This low correlation occurs because, on the one hand, the study area is large, so low-resolution (500 m) MODIS images are used, and clear cloud-free data are rare in this area (convective cloud activity is frequent on the Tibetan Plateau) [54]. On the other hand, the complex topography of the study area, the small number of meteorological stations, and the large and uneven distribution of stations limit the simulation of the real meteorological environment to a certain extent, resulting in the low correlations between temperature and precipitation and SOS. However, the significance test also showed that temperature and precipitation had the largest areas of significance in early May and late April, respectively, and largely overlapped with the areas of greater correlation in terms of geographic distribution, indicating that our findings still reveal significant spatial differences and trends in the influence of climatic factors (temperature and precipitation) on the SOS of grassland vegetation in the study area.

4.3. SOS Response to Topographic Factors

The SOS in Qinghai Province is influenced by terrain factors and is later than that of other ecosystems at the same latitude due to the high-altitude and low-temperature conditions [55]. In addition to climatic factors, the SOS in small areas is influenced by the combination of slope, aspect and altitude to some extent [14]. We found that the SOSs of the shady slopes, mainly on the eastern and northern slopes, were earlier than those of the sunny slopes, mainly on the southern and western slopes. Similarly, Chen et al. [56] found that the SOS of grasslands on the northern slope was earlier than that on the southern slope and that there was a higher concentration of SOS on the northern slope than on the southern slope on the QTP, especially in the context of climate “warming and drying”. Although the effective radiation and heat were higher on the sunny slopes than on the shady slopes, the soils were drier on the sunny slopes due to the faster evaporation rate, while the shady slopes had more adequate moisture and a lower evaporation rate [57]; thus, the shady slopes were more favorable for vegetation growth.

Slope is one of the main terrain factors that causes spatial differences in soil nutrients [58]; the steeper the slope is, the faster the soil nutrient loss is, and then the worse the soil texture and the lower the vegetation cover are [59]. We found that the SOS in Qinghai Province was mainly distributed in the region with slopes <5°, i.e., the steeper the slope was, the smaller the area of SOS was. Additionally, with the increase in slope, the SOS was more concentrated, and the soil conditions were more suitable for vegetation growth in the region with a gentle slope. The steeper slopes were mainly distributed in the high-altitude region, so the steeper the slope was, the smaller the area of the SOS was and the shorter the optimal time for vegetation growth was. Therefore, the overall relationship between the SOS and slope was pyramidal.

Altitude is one of the key factors that determines the growth of vegetation [17]. On the one hand, altitude causes changes in hydrothermal conditions [60]; on the other hand, it affects biomass and species richness [61]. The grasslands in Qinghai Province are mainly distributed in the areas where the altitude is higher than 3500 m, and the higher the altitude is, the worse the vegetation growth conditions. At this time, soil microbial activity lays the foundation for improving the ability of grasslands to cope with the severe environment in high-altitude areas [62]; therefore, the area of SOS was the largest in the range of 4500–5000 m, showing an inverted pyramid shape. The SOS in the low-altitude area most commonly occurred on the 130th day, and with increasing altitude, low temperature became the main limiting factor of SOS in the high-altitude area [63]; thus, when the altitude was more than 4000 m, the proportion of SOS occurrences on the 140th day increased gradually.

Remote sensing tools have been widely used for vegetation phenology monitoring [64], and the method provides data support to realize long time series and large-scale phenology studies. The advancement of the SOS due to global warming has been confirmed by numerous studies, and predicting the response of the SOS to climate change has become a hot topic globally [65]. The response of the SOS of grassland vegetation to climate factors is complex. Some studies have shown that precipitation in temperate grasslands of North America has no significant effect on the flowering date of vegetation [66]; however, other experiments have shown that reduced water supply advances SOS and flowering date to some extent. Moreover, in semiarid grasslands, precipitation largely controls the direction of phenological changes [67]; however, in meadow grasslands, precipitation is dominated by temperature, which indicates the variability among different geographical regions. To better represent the relationship between climate change and SOS, researchers have also used modeling and other means to estimate phenology [68]. However, due to the complex relationships between phenology and meteorological variables, there are still controversies among the results of various studies, and further research in this area is needed.

5. Conclusions

From 2000 to 2020, the SOS in Qinghai Province was concentrated between the 110th and 150th days, with an average SOS on the 128th day, and the overall SOS in the study area showed an advancing trend. The spatial distribution of the average SOS from 2000–2021 shows that the early zone is mainly concentrated in the southeastern low-elevation region and the late region is mainly in the high-altitude region in a “V”-shaped belt distribution, showing an overall pattern of “early in the east and late in the west”.

The temperature, precipitation and SOS all showed inverse fluctuation trends before and after the midpoint of 2010, and the overall trend of temperature decreased, while the fluctuation trends of precipitation and SOS were more consistent. Both showed an increasing/advancing trend, and precipitation affected the SOS earlier than temperature. When temperature became the limiting factor of the SOS, precipitation had a greater effect on the SOS.

The overall SOS on shady slopes (east and north slopes) is earlier than that on sunny slopes (south and west slopes), and the SOS area decreases as the slope continues to increase when the SOS is delayed. The SOS area was largest in the 4500–5000 m interval and decreased at lower and higher altitudes. The SOS occurred later as altitude increased.