Seasonal Variation of Vegetation and Its Spatiotemporal Response to Climatic Factors in the Qilian Mountains, China

: The purpose of this study is to reveal the seasonal difference in vegetation variation and its seasonal response to climate factors in the Qilian Mountains (QM) under the background of global warming. Based on the MOD13 A2 normalized difference vegetation index (NDVI) data and meteorological data, this study analyzed the spatiotemporal dynamics and stability of vegetation in different seasons by using the mean value method, trend analysis and stability analysis method, and discussed their seasonal responses to climatic factors based on the correlation analysis method. The results show that the vegetation cover in the QM experienced a signiﬁcant upward trend in the past 21 years, but there were obvious spatial differences in vegetation change in different seasons. The growth rate of vegetation in summer was the fastest, and summer vegetation provided the most signiﬁcant contribution to the growing season vegetation. The order of vegetation stability in the QM among the seasons was growing season > summer > spring > autumn. The vegetation change was obviously affected by temperature in spring, while it was mainly controlled by precipitation in the growing season and summer. The response of vegetation to climatic factors was not signiﬁcant in autumn. Our results can provide important data support for ecological protection in the QM and socioeconomic development in the Hexi Corridor. and W.K.; J.Z.; and Formal and X.J.; Resources, W.K.; and and


Introduction
The relationship between global change and terrestrial ecosystems (GCTE) is a vital component of global change research [1,2]. As the most important and sensitive component of terrestrial ecosystems to climate change [3][4][5], vegetation connects the soil, hydrosphere and atmosphere and promotes material migration and energy exchange in each layer [6][7][8]. Meanwhile, vegetation plays an irreplaceable role in climate regulation, the terrestrial carbon cycle, the maintenance of biodiversity, the prevention of desertification and water conservation through related biophysical and chemical processes (i.e., photosynthesis, transpiration, etc.) [9][10][11][12]. Global and regional climate change often affects the growth of surface vegetation by changing its structure and function [13,14]. The growth status of surface vegetation can also have a positive feedback effect on corresponding climate change [2,15,16]. Vegetation dynamics and their responses to climate change have become popular topics in global change research [17,18].
Related studies show that from the early 1980s to the late 1990s, vegetation activity in the high latitudes of the Northern Hemisphere increased significantly due to global warming [19][20][21]. However, the latest studies found that in the decade after 2000, the but little attention has been given to the seasonal difference in vegetation variation, and there have been few studies on the seasonal variation characteristics of vegetation in this region. Relevant studies show that there is a certain deviation or one-sidedness in simply selecting vegetation cover changes annually or in a certain season to represent vegetation time series variations, which makes it difficult to fully reflect the difference in vegetation change [52][53][54]. Some studies have shown that seasonal variation in climate is one of the key controlling factors for the formation of seasonal vegetation variation characteristics [3,55]. In addition, most studies focus on the interannual scale when discussing the relationship between vegetation dynamics and climate factors and ignore the impact of climate factors in different seasons on different vegetation growth stages [38,56,57]. Although a few studies considered the seasonal response of vegetation to climate change, they often only focus on the simultaneous impact of climate change on vegetation and do not consider the effect of climate change on vegetation growth in adjacent seasons [58][59][60][61]. In view of the problems and deficiencies in previous studies, this research analyzed the spatial distribution characteristics and spatial-temporal variation of vegetation in the growing season, spring, summer and autumn during the past 21 years, and revealed the seasonal differences and trends of vegetation in the QM. The research results comprehensively revealed the details of the seasonal differences and changing trends of vegetation in the QM. At the same time, the responses and differences between vegetation and climate factors in different seasons were discussed from the seasonal scale.
Based on the above analysis, to better understand the spatiotemporal evolution and stability of vegetation on the seasonal scale, and to reveal the controlling effect of climate change on the seasonal rhythm of vegetation growth, the following objectives were established for this study: (1) analyze the interannual variation trends of NDVI in different seasons; (2) explore the spatiotemporal variation of NDVI in different seasons; (3) conduct a stability assessment of vegetation dynamics in different seasons; and (4) determine the response difference of NDVI to climatic factors on the seasonal scale. Based on the above objectives, we will finally answer the following questions: first, what changes have taken place in the vegetation of the QM on the seasonal scale since 2000? Second, what are the spatial differences in vegetation stability in different seasons over the past 21 years? Third, what is the relationship between climatic factors and vegetation growth in different seasons in this region? The results will put forward a scientific theoretical basis for restoring and protecting vegetation and reasonable measurements in the QM to tackle climate change.

Study Area
The QM is the functional area of ecological protection and the source of many inland rivers in the Hexi Corridor, known as the "Ice Source Reservoir". This region is located in the intersection zone of the Qinghai-Tibet Plateau, Mongolia-Xinjiang Plateau and Loess Plateau and crosses the border between western Gansu Province and northeastern Qinghai Province ( Figure 1). The QM plays an important role in the division of natural climatic zones in Northwest China [40,62]. The terrain of the QM decreases from northwest to southeast. This area includes a series of northwest to southeast mountains and valleys, most of which are between 3000 and 5000 m (Figure 1a). The climate of this region is a typical alpine continental climate, with complex natural conditions and significant spatial differences in hydrothermal distributions. The precipitation in the QM presented a step-like spatial distribution pattern with more in the east and less in the west, with a varied range of 28.59-655.19 mm (Figure 1c). The annual precipitation was between 400 and 656 mm in the central and eastern regions and between 200 and 400 mm in the central and western regions. Due to the low temperatures throughout the year, approximately 15% of the annual precipitation is in the form of snow. A large number of modern glaciers are distributed in the mountain region above 3500 m. The rivers are mainly supplied by glacier meltwater and mountain precipitation and radiate around the center of the QM. The annual mean temperature ranged from −17.48 • C to 7.80 • C, and the temperature was higher at the edge of the QM but lower in the central and western mountains (Figure 1d). The annual mean temperature in the area above 4000 m was below 0 • C, and the spatial distribution of the annual mean temperature was highly consistent with the topography (altitude) in the study area. Under the combined influence of the southeast monsoon and topography, the hydrothermal conditions are redistributed, which leads to the unique vertical zonal characteristics of vegetation distribution (Figure 1b). From low to high elevations, there are desert steppes, subalpine steppes, forests, subalpine shrub meadows, alpine meadows and sparse alpine meadows [47,63]. distributed in the mountain region above 3500 m. The rivers are mainly supplied by glacier meltwater and mountain precipitation and radiate around the center of the QM. The annual mean temperature ranged from −17.48 °C to 7.80 °C, and the temperature was higher at the edge of the QM but lower in the central and western mountains (Figure 1d). The annual mean temperature in the area above 4000 m was below 0 °C, and the spatial distribution of the annual mean temperature was highly consistent with the topography (altitude) in the study area. Under the combined influence of the southeast monsoon and topography, the hydrothermal conditions are redistributed, which leads to the unique vertical zonal characteristics of vegetation distribution (Figure 1b). From low to high elevations, there are desert steppes, subalpine steppes, forests, subalpine shrub meadows, alpine meadows and sparse alpine meadows [47,63].

Data Sources and Preprocessing
The MOD13A2 NDVI product was used as the data source in this study, which was obtained from the National Aeronautics and Space Administration (NASA) (https://ladsweb.modaps.eosdis.nasa.gov/search/, accessed on 1 February 2021). The time span of the data in this study is from 2000 to 2020, the time resolution is 16 days, and the spatial resolution is 1 km. According to the temporal resolution, 23 NDVI images can be obtained every year. The Modis Reprojection Tool (MRT) and ArcGIS software were used to preprocess NDVI data through mosaics, projection transformation and image clipping. The

Data Sources and Preprocessing
The MOD13A2 NDVI product was used as the data source in this study, which was obtained from the National Aeronautics and Space Administration (NASA) (https: //ladsweb.modaps.eosdis.nasa.gov/search/, accessed on 1 February 2021). The time span of the data in this study is from 2000 to 2020, the time resolution is 16 days, and the spatial resolution is 1 km. According to the temporal resolution, 23 NDVI images can be obtained every year. The Modis Reprojection Tool (MRT) and ArcGIS software were used to preprocess NDVI data through mosaics, projection transformation and image clipping. The projection of the original image was converted from sinusoidal to Albers equal-area conic. To suppress the interference of noise on the accuracy of NDVI data, the Savitzky-Golay (S-G) smoothing filter method was adopted to further denoise and smooth the NDVI images [37]. Removing the residual clouds or other adverse effects in the subpixels can make the data more genuinely reflect the seasonal variation in vegetation. To avoid the influence of non-vegetation areas (desert and insufficient vegetation coverage) on the calculation results of the vegetation change in the whole region, according to relevant studies and the actual situation of the research area [15,64,65], the areas where the multiyear mean value of the growing season NDVI (GSNDVI) was less than 0.1 were defined as the non-vegetation areas for mask removal. On this basis, the NDVI time-series datasets in different seasons were reconstructed to analyze the vegetation change trends and their response to climate in different seasons.
The meteorological data (monthly precipitation and temperature) from 2000 to 2019, with a spatial resolution of 0.0083333 • (approximately 1 km), were obtained from the National Earth System Science Data Center (NESSDC) (http://www.geodata.cn/, accessed on 5 January 2021). The accuracy of the dataset has been verified before use [66]. ArcGIS software was used to convert the data format and values, clip and determine the seasonal composition, and resample the original data to 1 km × 1 km. Time series datasets of precipitation and mean temperature in the annual period, growing season and other seasons were created to analyze the spatiotemporal relationship between climatic factors and vegetation change.
Digital elevation model (DEM) data came from the provincial data of the Resource and Environmental Science and Data Center (RESDC) (https://www.resdc.cn/, accessed on 28 July 2020). After collating and splicing, the 90 m DEM data of the QM was generated, and it was mainly used to create the distribution map of the elevation in the study area ( Figure 1a).
The vegetation type data were derived from the "Qilian Mountain area vegetation distribution data" of the National Cryosphere Desert Data Center (NCDDC) (http://www. ncdc.ac.cn/portal/, accessed on 16 December 2020). Its format is a polygon shapefile, and the scale is 1:1,000,000. The data include 11 primary vegetation types and 54 secondary vegetation types, reflecting the distribution and zonality of different dominant vegetation communities. In the QM, there are nine primary vegetation types (Figure 1b).

Determination of the Growing Season
The characteristics needed for plant phenological recognition should be determined according to the law of time evolution for many years, which can improve the robustness of phenological remote sensing recognition methods [67]. Therefore, based on the NDVI time-series data reconstructed by S-G filtering, the NDVI average values of each DOY were calculated by taking the corresponding NDVI values in all years [68]. The time points corresponding to the maximum and minimum values of change rates represent the start of the growing season and the end of the growing season, respectively, and the length of the time period between them indicates the length of the growing season ( Figure 2). Here, the maximum positive and negative change rates of the NDVI corresponded to the 113th and 289th days of the year, respectively. Therefore, we defined the 113th (late April)-289th (mid-October) day as the growing season in this research, and the mean value of the NDVI during this period was defined as the value of the GSNDVI.

Mean Value Method
To more comprehensively investigate the change characteristics of vegetation cover in different periods. The mean value method was adopted to calculate the average NDVI of the annual growing season and each season from 2000 to 2020, and the annual GSNDVI and NDVI of each season were obtained. Their calculation formula is as follows [40]: where NDVI s represents the mean NDVI in different seasons. m and n represent the GSNDVI and each season NDVI starting from the m th 16-day to the n th 16-day, respectively. NDVI ij denotes the NDVI value of the i th 16-day in the j th year (m ≤ i ≤ n, 1 ≤ j ≤ 21). K represents the number of periods of 16-day NDVI in different seasons. In this study, the time spans of different seasons are as follows: growing season (DOY113-DOY289), spring (DOY065-DOY145), summer (DOY161-DOY241) and autumn (DOY257-DOY321).

Mean Value Method
To more comprehensively investigate the change characteristics of vegetation cover in different periods. The mean value method was adopted to calculate the average NDVI of the annual growing season and each season from 2000 to 2020, and the annual GSNDVI and NDVI of each season were obtained. Their calculation formula is as follows [40]: where represents the mean NDVI in different seasons. m and n represent the GSNDVI and each season NDVI starting from the m th 16-day to the n th 16-day, respectively. NDVIij denotes the NDVI value of the i th 16-day in the j th year (m ≤ i ≤ n, 1 ≤ j ≤ 21). K represents the number of periods of 16-day NDVI in different seasons. In this study, the time spans of different seasons are as follows: growing season (DOY113-DOY289), spring (DOY065-DOY145), summer (DOY161-DOY241) and autumn (DOY257-DOY321).

Trend Analysis Method
The unary linear regression analysis method can better reflect the variation trend of NDVI with time from the pixel scale [69]. Compared with the difference method, this method eliminates the influence of occasional anomaly factors at the end of the study period and more effectively reflects the long-term trend of vegetation cover [70]. This method was adopted to simulate the interannual changing trends of the NDVI in the different seasons from 2000 to 2020. Its formula is [8]:

Trend Analysis Method
The unary linear regression analysis method can better reflect the variation trend of NDVI with time from the pixel scale [69]. Compared with the difference method, this method eliminates the influence of occasional anomaly factors at the end of the study period and more effectively reflects the long-term trend of vegetation cover [70]. This method was adopted to simulate the interannual changing trends of the NDVI in the different seasons from 2000 to 2020. Its formula is [8]: where Θ _slope is the slope of the NDVI trend regression in the different seasons, i represents the corresponding year (1-21), and n is the span of the research period. NDVI i represents the NDVI value in the different seasons for the i th year. Θ _slope > 0 indicates an increasing trend of NDVI in each season during the research period and vice versa. The F test was adopted to evaluate the significance of NDVI trends, and the significance level was set to p < 0.05.

Stability Analysis
The coefficient of variation (CV) is a statistic that measures the degree to which the observed values deviate from the normal [71]. In this study, the CV was adopted to reflect the relative fluctuation of the NDVI during the past 21 years. The higher the value is, the greater the disturbance intensity and the more unstable the vegetation. The smaller the value is, the more stable the vegetation is [72]. Its formula is: where NDVI σ is the standard deviation of the annual NDVI and NDVI is the average value of the NDVI for multiple years. Based on the relevant research and the actual situation [73,74], the stability of the NDVI variation in the growing season and different seasons was divided into six grades by using the geometric interval classification method ( Table 1).

Correlation Analysis Method
By analyzing the correlation between geographical elements, we can accurately explain the closeness of their relationship [75]. This study analyzed the spatial correlation between the NDVI and climatic factors in different seasons based on the correlational analysis method. Its formula is [76,77]: , the NDVI value in most areas of the QM was below 0.3 in spring, which was lower than the values of the other seasons. The NDVI value in summer was higher than that in the whole growing season, which was the most vigorous stage of vegetation growth in the QM. The NDVI value of autumn was between spring and summer, but the overall distribution of vegetation was similar to that in spring. Thus, the variation range and distribution of NDVI varied greatly in different seasons. The eastern part of the QM had more precipitation due to the southeast monsoon, while the precipitation in the central and western parts was lower, mainly due to the influence of the dry airflow from the northwest. This is the main reason why the vegetation coverage in the eastern part of the QM is higher than that in the western part.
GSNDVI was mainly between 0.1 and 0.2, which is the transition zone from grassland to desert, and the vegetation types were mainly alpine steppe and alpine desert. Most of the western part of the QM is mainly unused land, such as deserts, bare rock and glaciers, and their NDVI values were all less than 0.1 in the growing season. Therefore, these regions were masked as no vegetation areas in this study. Viewed from the spatial distribution of the NDVI in each season (Figure 3a-d), the NDVI value in most areas of the QM was below 0.3 in spring, which was lower than the values of the other seasons. The NDVI value in summer was higher than that in the whole growing season, which was the most vigorous stage of vegetation growth in the QM. The NDVI value of autumn was between spring and summer, but the overall distribution of vegetation was similar to that in spring. Thus, the variation range and distribution of NDVI varied greatly in different seasons. The eastern part of the QM had more precipitation due to the southeast monsoon, while the precipitation in the central and western parts was lower, mainly due to the influence of the dry airflow from the northwest. This is the main reason why the vegetation coverage in the eastern part of the QM is higher than that in the western part.  Through the comparative analysis of the interannual variation trend and seasonal difference in vegetation growth in the last 21 years (Figure 4), we found that the vegetation increased significantly (p < 0.01) in each season, in which the growth rate of the NDVI in summer was the fastest (0.032/10yr), and the minimum (0.334), maximum (0.426) and annual mean value (0.387) of NDVI were higher than those in the other seasons. The growth rate of the NDVI in spring was the lowest (0.015/10yr), and the minimum (0.149), maximum (0.194) and annual mean values (0.167) of the NDVI were all lower than those of the other seasons. The main reason for this result is that the study area is located at a high altitude, the temperature is lower in spring, and the vegetation returns to green later. Affected by the low temperature, the fluctuation range of vegetation was small, and its value was lower than that of other seasons. The study found that the interannual fluctuation trend of the NDVI in summer was similar to that in the growing season, while the trend in autumn was different from that in the growing season and other seasons. The maximum value of the NDVI in the other seasons occurred in 2018 except in autumn. This result indicates that the contribution of vegetation to the growing season in summer was significantly higher than that in autumn. In addition, by comparing the fluctuation range of the NDVI trend line in different seasons, it was found that the fluctuation range of the in autumn was different from that in the growing season and other seasons. The maximum value of the NDVI in the other seasons occurred in 2018 except in autumn. This result indicates that the contribution of vegetation to the growing season in summer was significantly higher than that in autumn. In addition, by comparing the fluctuation range of the NDVI trend line in different seasons, it was found that the fluctuation range of the NDVI in autumn was the largest, indicating that the vegetation cover in autumn was the most unstable.

Spatiotemporal Variation of the NDVI in the QM from 2000 to 2020
Although the interannual and seasonal NDVI showed an increasing trend from 2000 to 2020 in the QM, the dynamic change in the NDVI was significantly different in terms of its spatial distribution. Therefore, it was necessary to simulate and evaluate the spatial dynamic trend of the NDVI based on the changing trend of the interannual and seasonal NDVI on the time scale.

Spatiotemporal Variation of the NDVI in the QM from 2000 to 2020
Although the interannual and seasonal NDVI showed an increasing trend from 2000 to 2020 in the QM, the dynamic change in the NDVI was significantly different in terms of its spatial distribution. Therefore, it was necessary to simulate and evaluate the spatial dynamic trend of the NDVI based on the changing trend of the interannual and seasonal NDVI on the time scale.

Spatiotemporal Variation Trend of the GSNDVI and Its Significance
The spatial variation trend of the vegetation GSNDVI was fitted by the trend analysis method in the QM from 2000 to 2020. As seen from the spatial distribution map and statistical results of the GSNDVI variation trend ( Figure 5, Table 2), the overall vegetation in the QM showed an upward trend in the past 21 years, and the increased area accounted for 93.27% of the total vegetation cover region (TVCR), which was much larger than that of the vegetation reduction area. The area showing a significant increase (p < 0.05) in vegetation accounted for 75.06% of the TVCR, which was mainly distributed in the western and marginal areas of the central-eastern regions of the QM. The area of vegetation reduction accounted for 6.72% of the TVCR and was mainly concentrated in the central part of the QM and the surrounding area of Qinghai Lake, in which the significantly reduced (p < 0.05) vegetation area accounted for only 0.7% of the TVCR, dispersed in northern Guide County, northeastern Gonghe County and northeastern Tianjun County. In conclusion, the vegetation cover in most regions of the QM improved significantly from 2000 to 2020, only a few regions had vegetation degradation, and the significantly improved areas were much larger than the degraded areas. < 0.05) vegetation area accounted for only 0.7% of the TVCR, dispersed in northern Gui County, northeastern Gonghe County and northeastern Tianjun County. In conclusio the vegetation cover in most regions of the QM improved significantly from 2000 to 20 only a few regions had vegetation degradation, and the significantly improved areas we much larger than the degraded areas.  In addition to discussing the spatial variation trend of the GSNDVI, this study a compared the spatial variation trends of the vegetation NDVI in different seasons. It c be seen from the distribution map that their trends had some similarities (Figure 6abut the NDVI change rates in different seasons had obvious differences in spat

Spatiotemporal Variation Trends of the NDVI in Different Seasons and Their Significance
In addition to discussing the spatial variation trend of the GSNDVI, this study also compared the spatial variation trends of the vegetation NDVI in different seasons. It can be seen from the distribution map that their trends had some similarities (Figure 6a-c), but the NDVI change rates in different seasons had obvious differences in spatial distribution. According to the change trends and significance grades (Figure 6a-f, Table 3), there were significant (p < 0.05) degradation areas in the central QM in spring, which were mainly distributed in northwestern Gangcha County, northern Tianjun County and southwestern Qilian County. The significant degradation area accounted for 1.57% of the TVCR, which was larger than that in the growing season (0.70%), summer (0.72%) and autumn (0.61%). The spatial distribution of the vegetation change trend and significance in summer was consistent with that in the growing season. The area of vegetation increased significantly (p < 0.05) in autumn, accounting for 44.17% of the TVCR, which was significantly smaller than that in spring (62.95%) and summer (72.74%). This region mainly included the southeastern and northwestern marginal zones of the QM, while the change was not obvious in most of the central regions. The vegetation degradation area in autumn was mainly distributed in northern Qinghai Lake, and the significant (p < 0.05) degradation area was distributed primarily in southwestern Menyuan Hui Autonomous County, which is a famous rape and highland barley planting base, as well as a tourist scenic area, so it was disturbed by human activities. The above analysis concluded that the area of vegetation degradation in spring was the widest, the changing trend of vegetation in summer was essentially consistent with that in the growing season, and the area of vegetation with a significant increasing trend in autumn was the smallest.
in autumn was mainly distributed in northern Qinghai Lake, and the significant (p < degradation area was distributed primarily in southwestern Menyuan Hui Autonom County, which is a famous rape and highland barley planting base, as well as a to scenic area, so it was disturbed by human activities. The above analysis concluded the area of vegetation degradation in spring was the widest, the changing trend of v tation in summer was essentially consistent with that in the growing season, and the of vegetation with a significant increasing trend in autumn was the smallest.

Stability of the Vegetation in Different Seasons from 2000 to 2020
The coefficient of variation method was used to calculate the CV of NDVI in the growing season and other seasons to reflect the stability of vegetation in the QM. According to the statistical results of CV classification in different seasons (Table 1), the proportions of stable vegetation regions in the growing season, spring, summer and autumn were 76.27%, 60.76%, 68.68% and 46.37%, respectively. Among the stability grades in different seasons, the proportion of low stable areas was the highest. The stable area in autumn was the smallest, and more than half of the total area was unstable. This result indicated that the fluctuation of vegetation in autumn was the largest, and the stability was also the lowest, followed by that in spring. This conclusion is consistent with that obtained from Figure 4d in Section 3.2: the proportion of highly stable regions was the smallest, and the stable regions were mainly moderate and low stability. Based on the statistical results of the unstable regions (Table 1), the proportions of moderately and highly unstable regions in autumn were 12.70% and 23.10%, respectively, which were much higher than those in other seasons, followed by those in spring. The proportions of moderate and highly unstable areas in the growing season were the smallest, indicating that the vegetation stability was the highest in the growing season. This scenario was also the main reason why we used the GSNDVI to reflect the interannual variation in vegetation in this study. According to the results of vegetation stability in each season, we found that the vegetation is relatively stable in the growing season and summer, but unstable in spring and autumn.
From the spatial distribution map of the NDVI variation coefficient and stability grades in the QM (Figures S1 and 7), the vegetation variability was smaller in the eastern region and larger in the western region, and the vegetation stability in spring and autumn was lower than that in the growing season and summer. The vegetation stability in summer was basically consistent with that in the growing season. The moderate and highly stable regions were mainly distributed in the high vegetation cover regions in the eastern QM, and the moderate and highly unstable regions were sporadically distributed in the low vegetation cover areas that transition from grassland to desert in the northwest. In addition, the stability of the southern margin and northeastern corner of the QM was relatively low. The unstable area in autumn was significantly larger than that in spring, especially in the highly unstable areas. The moderate and highly unstable areas were mainly the low vegetation cover areas in the middle and western parts and most of the slope areas below the snow line in the northern part of Qinghai Lake. The moderate and highly stable areas were primarily located in eastern QM and the surrounding area of Qinghai Lake. Based on the spatial distribution and statistical results of vegetation stability in different seasons ( Figure 7, Table 1), it was concluded that the order of vegetation stability was the growing season > summer > spring > autumn, and the sensitivity of vegetation in autumn and spring was the highest. Therefore, attention should be given to the changes in vegetation in spring and autumn in the QM.

Interannual and Seasonal Variations in Climatic Factors
According to the variation trends of precipitation and temperature in the interannua period, the growing season and other seasons were discussed in terms of the vegetation covered region of the QM. As shown in Figure 8a, the annual precipitation showed a sig

Interannual and Seasonal Variations in Climatic Factors
According to the variation trends of precipitation and temperature in the interannual period, the growing season and other seasons were discussed in terms of the vegetationcovered region of the QM. As shown in Figure 8a, the annual precipitation showed a significant (p < 0.05) upward trend from 2000 to 2019; the upward trend was not obvious from 2000 to 2012 but showed a sharp upward (p < 0.01) trend from 2013 to 2019. This result was generally consistent with the changing trend of the vegetation. The annual mean temperature showed a slow upward trend, but the changing trend was not obvious, with a slight fluctuation from 2000 to 2005 and a significant fluctuation from 2006 to 2019 (Figure 8f). Viewed from the variation trends of precipitation in different seasons (Figure 8a-e), the trends in the annual, growing season and summer showed significant (p < 0.05) upward trends overall, and their fluctuation trends were similar. However, the changing trend of precipitation in spring was not significant. Contrary to other seasons, autumn had a downward trend in precipitation, but the trend was not significant. Viewed from the variation trends of temperature in different seasons (Figure 8f-j), the spring temperature increased significantly (p = 0.05) and at the fastest increasing rate (slope = 0.040). The fluctuation trend was generally consistent with that of the NDVI in spring, indicating that the spring vegetation may be obviously affected by the spring temperature. The average interannual temperature, growing season, summer and autumn temperatures all showed an upward trend, but the trends were not significant. The variation trend of the summer temperature was very flat (Figure 8i). The NDVI and precipitation analysis results showed that both vegetation and precipitation increased significantly in summer. Combined with Figure 4c, it is found that the fluctuation shape of NDVI was similar to that of summer precipitation, so it can be speculated that the increase in vegetation may be related to precipitation in summer.

Correlation between the GSNDVI and Climatic Factors Interannually and in th Growing Season
According to the correlation distribution map (Figure 9), the distributions of the relation between the GSNDVI and annual precipitation and growing season precipit were generally consistent (Figure 9a,c). Their positive correlation areas were much l

Correlation between the GSNDVI and Climatic Factors Interannually and in the Growing Season
According to the correlation distribution map (Figure 9), the distributions of the correlation between the GSNDVI and annual precipitation and growing season precipitation were generally consistent (Figure 9a,c). Their positive correlation areas were much larger than the negative correlation areas; the significant positive correlation (p < 0.05) areas accounted for 58.87% and 60.67% of the TVCR, respectively, which were primarily distributed in the central and western QM, around Qinghai Lake and northeast of Wushaoling (Figures 9 and S2). The vegetation cover types in these regions were mainly steppe and desert steppe, with high temperatures and intense evaporation in the growing season, so precipitation was the dominant climatic factor controlling vegetation growth in these regions and had a strong correlation with the GSNDVI. The area with a negative correlation between the GSNDVI and annual precipitation and growing season precipitation was small, accounting for only 8.41% and 7.39% of the TVCR, respectively. Their significant negative correlation areas were less than 0.5%. They were mainly sporadically distributed in the eastern agricultural and pastoral areas (Figures 9 and S2). These areas were affected by human activities, and the farmland was supplied primarily by rivers, which were less dependent on precipitation. The correlation between the GSNDVI and annual mean temperature and growing season temperature was mainly positive (Figure 9b,d). The area where the GSNDVI was significantly and positively correlated (p < 0.05) with the growing season temperature was greater than the annual mean temperature; their significant positive correlation areas accounted for 14.42% and 5.84% of the TVCR, respectively ( Figure S2), while the other regions were not significantly affected by temperature. It can be concluded that the response of the GSNDVI to growing season temperature was more sensitive than to the annual mean temperature; that is, the GSNDVI was more closely related to concurrent temperature than to annual mean temperature. To explore the differences in the relationship between vegetation and climate change in different seasons, the correlation coefficients between the NDVI and its corresponding

Spatial Response of the Seasonal NDVI to Seasonal Climatic Factors
To explore the differences in the relationship between vegetation and climate change in different seasons, the correlation coefficients between the NDVI and its corresponding climatic factors were calculated in each season ( Figure 10). From the spatial distribution of the correlation coefficients in different seasons, it can be seen that the response of the NDVI to climate change has significant seasonal differences in space. In spring, NDVI was significantly (p < 0.05) and positively correlated with precipitation in 7.18% of the TVCR, which was primarily distributed in the surrounding area of Xining, east of Gonghe County and the junction of Shandan County, Yongchang County and Sunan Yugur Autonomous County in the north (Figures 10 and S3). NDVI was significantly (p < 0.05) and negatively correlated with precipitation in 7.38% of the TVCR, which was mainly distributed in the western part of the study area and the Lenglongling region (Figures 10 and S3). Under the control of low temperature, vegetation in these regions returned to green relatively late. The area where NDVI was significantly and negatively correlated with precipitation may be caused by the time lag effect of the NDVI response to precipitation. The correlation between the spring NDVI and spring temperature was the highest during the year, and the significant positive correlation (p < 0.05) areas were mainly concentrated in most parts of the central and eastern regions, accounting for 29.47% of the TVCR (Figures 10 and S3).
Sustainability 2022, 14, x FOR PEER REVIEW 16 analysis, the overall vegetation still increased despite the autumn precipitation decre slightly, which was closely related to the implementation of environmental prote projects in the QM. The NDVI in autumn was significantly and positively correlated the autumn temperatures only in the eastern region of the Laji Mountains and Da River Valley ( Figure S3), indicating that the increasing temperature in these regions certain promotional effect on vegetation growth in autumn.

Time Lag Effect of the NDVI Response to Climatic Factors in Different Seasons
This study not only analyzed the correlation between NDVI and correspondin matic factors pixel by pixel, but also conducted a statistical analysis on the correl between seasonal NDVI and climatic factors in adjacent seasons, and discussed the lag effect of NDVI on the seasonal response of climatic factors ( Table 4). The results The correlation between the NDVI and precipitation in summer was the highest during the year, and the proportion of significantly (p < 0.05) and positively correlated areas was 56.53%, which was primarily distributed in the central and western regions of the QM and Wushaoling (Figures 10 and S3). These areas belong to arid and semi-arid regions, and the main vegetation cover types are desert steppe and alpine grassland. Under the condition of no significant change in temperature, precipitation becomes the dominant factor controlling vegetation growth. The correlation between the NDVI and temperature in summer was not obvious overall, but the NDVI in summer was negatively correlated with precipitation in some parts of the northwest (Figures 10 and S3). From the previous analysis, the overall vegetation still increased despite the autumn precipitation decreasing slightly, which was closely related to the implementation of environmental protection projects in the QM. The NDVI in autumn was significantly and positively correlated with the autumn temperatures only in the eastern region of the Laji Mountains and Datong River Valley ( Figure S3), indicating that the increasing temperature in these regions has a certain promotional effect on vegetation growth in autumn.

Time Lag Effect of the NDVI Response to Climatic Factors in Different Seasons
This study not only analyzed the correlation between NDVI and corresponding climatic factors pixel by pixel, but also conducted a statistical analysis on the correlation between seasonal NDVI and climatic factors in adjacent seasons, and discussed the time lag effect of NDVI on the seasonal response of climatic factors ( Table 4). The results indicated that there was a significant correlation (p < 0.01) between the GSNDVI and summer precipitation, with a correlation coefficient as high as 0.662, followed by spring, with a correlation coefficient of 0.440 (p = 0.05). By comparing the correlation between the GSNDVI and temperature in each season, it was found that the GSNDVI had the highest correlation with spring temperature with a correlation coefficient of 0.411, but it did not pass the significance test, indicating that the GSNDVI was primarily affected by summer precipitation. The spring NDVI had little relationship with spring precipitation but had a significant (p < 0.01) correlation with spring temperature, and the correlation coefficient was as high as 0.704, indicating that the growth of vegetation in spring was susceptible to temperature variations. The summer NDVI had a strong correlation with both summer and spring precipitation, and their correlation coefficients were 0.658 (p < 0.01) and 0.491 (p < 0.05), respectively, indicating that vegetation variation in summer was not only affected by precipitation in the same period but also had a significant time lag relationship with spring precipitation. The increased rainfall in spring and summer promoted the growth of vegetation in summer, especially in the arid and semi-arid desert and desert steppe in the central and western regions, where the growth of vegetation was directly limited by precipitation. The summer NDVI was not significantly correlated with either spring or summer temperatures, indicating that the NDVI was primarily affected by precipitation in summer. The correlations between the autumn NDVI and summer precipitation and temperature were higher than those in spring and autumn, indicating that the response of autumn vegetation variation to summer precipitation and temperature had a certain time lag effect. Although the vegetation variation in autumn was influenced by summer precipitation and temperature to some extent, the effect was not significant.

The Difference in Fluctuation in the Changing Trend of the NDVI
The results of this study showed that the vegetation in the QM has increased as a whole since 2000, and the improved area of vegetation is much larger than the degraded area, which is consistent with previous studies [40,61,[78][79][80], even though they adopted different data sources and research methods. However, the range of interannual NDVI fluctuation was quite different in various studies. In our research, the GSNDVI was used to represent the interannual variation in vegetation, and the fluctuation range of the NDVI was mainly between 0.27 and 0.35 (Figure 4a). Dai et al. [49] analyzed the spatiotemporal variation of vegetation in the QM from 1999 to 2007 by using the annual average NDVI. The results showed that the vegetation had an increasing trend since 1999, and the fluctuation range of the annual average NDVI was mainly distributed at approximately 0.15. The main reason for the difference in NDVI fluctuation ranges is that the non-vegetation region was not masked in Dai's study, which leads to the low annual average NDVI value. Even in the eastern region with high vegetation cover, the NDVI fluctuates only from 0.25 to 0.3. Xu et al. [60] used the annual maximum NDVI to study the spatiotemporal changes of vegetation and its influencing factors in the QM from 2000 to 2010. The results showed that the annual maximum NDVI increased at a rate of 0.024/10 yr, and the fluctuation of NDVI ranged from 0.48 to 0.57, which was closest to the results of our study (0.026/10 yr). However, in Xu's study, the fluctuation range of the interannual NDVI was significantly larger than that of Dai and our results. This was mainly because Xu's research not only masked the non-vegetation region but also adopted the annual maximum NDVI. Even in the western desert steppe, the NDVI ranged from 0.3 to 0.4. The disadvantage of this method is that it is likely to be disturbed by extreme months in calculating the maximum NDVI, thus inflating the result. Therefore, it is necessary not only to remove the influence of non-vegetation areas to avoid low results but also to avoid high results caused by extreme values. Based on the results of our study, we suggest that it is more reasonable to remove non-vegetation areas and use the average NDVI in the growing season to monitor the spatiotemporal change in vegetation.
One of the most obvious differences between our study and previous studies is the definition of the growing season. A large number of studies have defined May-September as the growing season [40,48,81,82], and others have used April to October as the growing season [44,59,83]. Most studies in the QM used May to September as the growing season [56,58,79,84]. These studies all defined the growing season on a monthly scale. Meanwhile, some studies used the methods of remote sensing phenology to conduct more detailed research on the growing season and refine the beginning and end of the growing season to the scale of days [85][86][87]. However, because it is difficult to obtain NDVI and its corresponding climate grid data covering the whole region on the day scale, few studies use NDVI daily data to synthesize the GSNDVI and seasonal NDVI based on the phenological calculation results. Therefore, based on remote sensing products, this study took 16 days as the starting point combined with the phenological calculation method to determine the growing season based on the 16-day time scale. The length of the growing season is also the main reason for the different fluctuation ranges and variation trends of the GSNDVI. In addition, based on 16-day NDVI data, this paper also determined the time span of other seasons which is the first attempt in this region.

Differences in the Seasonal Response of the NDVI to Climate Factors
There are few seasonal studies on vegetation change and its response to climatic factors in this region, and there are still differences between them. Wu et al. [88] showed that the area of vegetation improvement in the QM in summer was the largest in the last 13 years, which was consistent with the results of our study. This study indicated that the improvement rate of vegetation in spring was the most obvious, which was inconsistent with our results. There may be two reasons: first, the different monitoring periods lead to the difference in the interannual change rate of NDVI. The second possibility is that the study did not mask the non-vegetation area. Song et al. [89] found that the vegetation improved significantly in all seasons except spring, especially in summer, and the NDVI in the growing season and summer was significantly correlated with precipitation. Our research was consistent with this conclusion. Qiu et al. [61] and Wu et al. [88] carried out detailed studies on vegetation changes and their influencing factors in this region. They discussed the correlation between NDVI and climatic factors in different seasons. Through comparative analysis with our research, it is found that the spatial distributions of the correlation between seasonal NDVI and concurrent climatic factors are basically consistent in other seasons except autumn. Different from our research, these two studies only analyzed the simultaneous response of seasonal NDVI and climatic factors but did not consider the relationship with other seasonal climatic factors and lacked a comparison with the growing season. Based on the problems existing in previous studies, we specifically discussed the seasonal response of vegetation change to climatic factors.

Limitations of the Study
Vegetation growth is affected by many factors, including not only precipitation and temperature but also humidity, solar radiation, altitude, soil, hydrology, human activities and other factors. Considering the importance of impact factors and the availability of data, this study only selected the two most critical climatic factors (temperature and precipitation) for related research, focusing on the response of vegetation change to the two main control factors. In addition, because the focus of this study is to explore vegetation dynamics and their response to climate change on a seasonal scale, their response and time lag effect to climate change on monthly and 16-day scales were not considered.

Conclusions
Combined with remote sensing and climate reanalysis data, this study investigated the spatial distribution, variation trend and stability characteristics of vegetation in the QM from multiple time scales and revealed the seasonal responses of vegetation to climatic factors. The main conclusions are as follows: (1) The vegetation cover in the east was noticeably higher than that in the west. The interannual and seasonal variations in NDVI showed an increasing trend in the QM from 2000 to 2020. The significantly improved areas were much larger than the degraded areas. The growth rate of NDVI was the fastest in summer (0.032/10 yr) and the slowest in spring (0.015/10 yr). Vegetation in summer had the most prominent contribution to the vegetation variation in the growing season, while vegetation in autumn had little contribution. The decreasing range of vegetation greenness was the widest in spring.
(2) The vegetation was relatively stable in the growing season and summer but unstable in spring and autumn. The stable area was mainly distributed in the eastern high vegetation cover area. The unstable area was mainly distributed in the low vegetation cover area that forms the transition from grassland to desert in the northwest.
(3) The vegetation response to climatic factors has significant spatiotemporal heterogeneity. The effect of precipitation on vegetation was more significant in arid and semi-arid regions of the central and western parts. The spring NDVI was obviously affected by spring temperatures (p < 0.01), while the growing season and summer NDVI were most closely related to summer precipitation (p < 0.01) and were also affected by the lag in spring precipitation (p < 0.05). This result indicates that the climatic factors affecting vegetation change were also different in different seasons.

Data Availability Statement:
The data presented in this study are available from the corresponding author upon reasonable request.