Quantitative Effects of Climate Change on Vegetation Dynamics in Alpine Grassland of Qinghai-Tibet Plateau in a County

: Alpine grassland in the Qinghai-Tibet Plateau is known to be sensitive to climate change. To quantify the impacts of climate change on alpine ecosystems at small scale, Wulan County in Qinghai Province was taken as the research object, and the relationships between vegetation dynamics and climate changes and the direct and indirect effects of climate factors on vegetation dynamics were analyzed using the methods of ordinary least squares, Pearson correlation analysis and path analysis, on the basis of MOD13A3 data and meteorological data from 2001 to 2020. The results showed that the Normalized Difference Vegetation Index (NDVI) in the growing season of the county and 5 vegetation types showed similar ﬂuctuation processes and relationships with climate factors during the study period. The growing season NDVI (GSN) of shrubland and desert steppe respectively were the highest and lowest. The yearly mean values of GSN over the county ranged from 0.151 to 0.264, and increased extremely signiﬁcantly with years at a rate of 0.0035 yr − 1 . Spatially, GSN gradually decreased from northeast to southwest, and 97.2% of the county area showed an increasing trend in GSN. With years, growing season evaporation (GSE) decreased extremely signiﬁcantly at a rate of 29.6 mm yr − 1 , while growing season average relative humidity (GSARH) showed a signiﬁcant increasing trend at a rate of 0.16% yr − 1 . The correlations and effects of GSE, GSARH, and growing season precipitation (GSP) on vegetation dynamics were weakened in turn. GSE was the main direct effect factor, and the latter two were the indirect effect factors through GSE. The total contribution rates of GSE, GSARH and GSP to vegetation dynamics was about 78.0%. However, growing season average temperature (GSAT) had little effect on vegetation dynamics. This study provides information for understanding the characteristics of vegetation dynamics of alpine grassland in the Qinghai-Tibet Plateau at a small scale. 944.2–1697.4 mm with a mean value of 1183.2 mm, and the minimum and maximum value appeared in 2019 and 2001 respectively; growing season average relative humidity (GSARH) was between 42.1% and 50.0% with a mean value of 46.1%, and the minimum and maximum value appeared in 2001 and 2019, respectively. The general linear regression analysis of the 4 climate factors revealed that the changing processes of GSAT and GSP had insigniﬁcant increase or decrease trends; GSE decreased extremely signiﬁcantly with the year having a reduction rate of 29.6 mm yr − 1 ( R 2 = 0.780, p = 0), and the GSE in 2020 had decreased by 786 mm compared with 2001; GSARH showed a signiﬁcant increasing trend with an increase rate of 0.16% yr − 1 ( R 2 = 0.237, p = 0.030).


Introduction
Global climate change has an extremely significant impact on surface vegetation cover. At the same time, the changing characteristics of surface vegetation cover will have corresponding feedback on regional climate and global climate systems [1][2][3][4]. As a key component of the terrestrial ecosystem [5,6], grassland is one of the most common vegetation types, providing various ecosystem and social service functions [7][8][9]. With a total area of 2.38 million square kilometers of grasslands in Northern China (about 9.92% of the world's total grasslands), they are not only the important strategic resources and the basis for developing animal husbandry production in China, but also the important ecological defense line to maintain ecological balance [10,11]. Therefore, understanding the vegetation dynamics of grassland and exploring its response to climate factors are very important and necessary for the protection of the grassland ecosystem in China.
The Qinghai-Tibet Plateau, known as the "Third Pole" of the earth, is an ecological barrier for China and even Asia dominated by alpine grassland [12]. Alpine grassland covers more than 66% of the area of the Qinghai-Tibetan plateau with the highest altitude and the largest area in the world, and plays an important role in biodiversity maintenance, carbon fixation and regional climate regulation [13,14]. Nevertheless, alpine grassland in the Qinghai-Tibet Plateau has a simple composition of vegetation community and suffers adverse living conditions, such as high altitude, barren land, lack of available water resources and sparse vegetation distribution [15][16][17], thus the stability of the ecosystem in alpine grassland is poor, and the vegetation growth is more sensitive to regional climate change [18][19][20]. With the continuous progress of earth observation technology and observation means, the research contents and methods for vegetation on climate change are becoming more and more diverse [21][22][23]. The Normalized Difference Vegetation Index (NDVI), a widely used indicator of vegetation cover, is derived from remote sensing and has been shown to be a sufficiently stable indicator for monitoring the intra-and inter-annual variations of vegetation greenness [24][25][26][27]. And the growing season NDVI (GSN) is often used by scholars to examine the relationships between vegetation dynamics and climate changes for alpine grasslands in the Qinghai-Tibet Plateau, for the reason that GSN has more accurate expression significance for vegetation dynamics [24,26,[28][29][30]. Guo et al. [28] reported that a decrease in growing season precipitation mainly affects the vegetation ecosystems with low precipitation and worse vegetation condition, however, an increase in growing season temperature could improve vegetation growth conditions in some higher/rugged regions of the Yellow River, China, on the basis of GSN generated by the 1 km advanced very high-resolution radiometer (AVHRR) NDVI data during 1990-2000. Hu et al. [29] concluded that water condition and temperature were found to be the most important factors affecting the variation of GSN in both the semi-arid and the semi-humid areas in the three-river-source region using NOAA/AVHRR 10-day composite NDVI data during 1982-2000. Bai et al. [30] found that climate warming benefits alpine vegetation growth in Three-River Headwater Region since 2000 with GSN acquired from multiple high-resolution satellite data. Huang et al. [31] reported that the effects of different climatic factors on alpine vegetation growth are significantly different in the Tibetan Plateau with the increase of the study period, based on GSN calculated from the MOD09A1 data product of a moderate-resolution imaging spectroradiometer (MODIS).
Throughout the above studies, most of them focus on large-scale study areas, such as headwater of the Yellow River, Three-River Headwater Region and Qinghai-Tibet Plateau. Few studies, however, have reported the quantitative analysis of the relationship between vegetation dynamics and climate change in the alpine grassland of the Qinghai-Tibet Plateau at a small scale. The climate characteristics in large-scale areas are more complex and changeable than those in small-scale areas because the climates have significant spatial heterogeneity with lots of microclimates in large-scale areas. For the microclimate differences, it is difficult to accurately identify the driving effects of climate change on vegetation dynamics in large-scale study areas [30,32]. While in small-scale areas, the spatial difference of climate can be almost ignored. And compared with the studies for alpine grassland at large scales, small-scale studies could reveal the driving mechanism of climate change on vegetation dynamics more deeply [33]. In addition, previous studies mostly focused on the impact of temperature and precipitation on vegetation change in alpine grassland, while less on other climate factors. The impact of multiple climate factors and their direct and indirect effects on vegetation dynamics are more complex and cannot be ignored, which are still not clear enough and need to be further explored.
Based on the above analysis, Wulan County, located in the northeast of the Qinghai-Tibet Plateau, was selected as the typical study area to quantitatively detect the effects of climate change on vegetation dynamics in alpine grassland of the Qinghai-Tibet Plateau at a small-scale area. Wulan County is rich in alpine grassland, the area of which accounts for about 88% of the county area. It is a typical alpine grassland pasture with animal husbandry as the main economic pillar. Affected by the Qaidam Basin Desert, the grassland ecosystem of Wulan County is more fragile and sensitive to climate change [34]. Besides, a county is the smallest unit for the management, utilization, and protection of alpine grassland in China [35]. This work is necessary for the utilization and sustainable development of alpine grassland in Wulan County. At the same time, it can also provide reference for the research of alpine grassland protection in other small regions of the Qinghai-Tibet Plateau. In order to explore the effects of the interaction of multiple climate factors on the vegetation dynamics, the Path analysis was introduced into our study, which is widely used for variable analysis in medicine, economy, biology, agriculture and other fields with good identification effects [2,36,37]. It is a technique to determine to what extent the independent variable affects directly the dependent variable, and to what extent the independent variable affects the dependent variable indirectly through the other independent variable. Path analysis can provide a theoretical basis for scientific statistical decision-making, which is useful for quantitatively evaluating the driving effects of temperature, precipitation, evaporation, and average relative humidity on GSN, so as to explore the main driving factors of vegetation growth in alpine grassland of the Qinghai-Tibet Plateau at small scale. Overall, the objectives were to: (i) accurately describe the change features of vegetation in alpine grasslands from 2001 to 2020 in the study area; (ii) assess the variations of 4 climate factors in alpine grasslands from 2001 to 2020 in the study area; and (iii) quantitatively identify the relationship between climate factors and vegetation dynamics, and the direct and indirect effects of climate factors on vegetation dynamics in the study area.

Study Area
Wulan County (97 • 01 E-99 • 27 E, 36 • 19 N-37 • 20 N) governs four towns with a total population of about 40,000 and covers a total area of 1.29 × 104 km 2 with an average altitude of about 4000 m in Qinghai Province, China ( Figure 1). It is one of the main study areas of our team's National Key Research and Development Program of China (grant number 2016YFC0400301), with strong availability and reliability of basic data, providing a good foundation for this study. Wulan County has a typical dry continental climate characterized by a short frostless period, sufficient sunshine, a large temperature difference between day and night, as well as rain and heat in the same season. The average annual precipitation is 187.5 mm, the average annual temperature is 1.6-3.3 • C, and the average annual water surface evaporation is 2000-3000 mm. The area of grassland in Wulan County is about 112.8 × 10 4 ha, and the dominant vegetation types are shrubland (8.9%), steppe (18.4%), meadow (11.7%), and desert steppe (57.0%), while forest, moor, cultural vegetation, and other alpine vegetation covering low size are merged into 'rest vegetation' (4.0%) ( Figure 2a) [38]. We conducted the on-site investigations on the grassland in Wulan County in 2018 and 2020, respectively. Figure 2b,c shows photos taken during the investigations, showing the two vegetation types and that sheep are the main livestock in the pastoral area of Wulan County. Desert steppe covers the largest area in Wulan County containing the western desert areas with sparse vegetation. The permanent non-vegetation area, e.g., water bodies, glaciers, and snow, were excluded in this study.

Data Collection and Pre-Processing
Vegetation status over the study area is measured in terms of the Normalized Difference Vegetation Index (NDVI), which is a widely used proxy of vegetation canopy greenness. In this study, the moderate-resolution imaging spectroradiometer MODIS-derived monthly composite product (MOD13A3) with a spatial resolution of 1 km was selected to generate a long time series of NDVI dataset from 2001 to 2020, which was pre-processed and released by the Atmosphere Archive and Distribution System (LAADS) of the National Aeronautics and Space Administration (NASA) (https://ladsweb.modaps.eosdis.nasa.gov/, accessed on 9 August 2021).

Data Collection and Pre-Processing
Vegetation status over the study area is measured in terms of the Normalized Difference Vegetation Index (NDVI), which is a widely used proxy of vegetation canopy greenness. In this study, the moderate-resolution imaging spectroradiometer MODIS-derived monthly composite product (MOD13A3) with a spatial resolution of 1 km was selected to generate a long time series of NDVI dataset from 2001 to 2020, which was pre-processed and released by the Atmosphere Archive and Distribution System (LAADS) of the National Aeronautics and Space Administration (NASA) (https://ladsweb.modaps.eosdis.nasa.gov/, accessed on 9 August 2021).
Meteorological data from 2001 to 2020, including temperature, precipitation, evaporation, and average relative humidity, were collected from 3 meteorological stations, namely, Wulan Station, Chaka Station, and Dulan Station, and the data were available from China Meteorological Data Service Center (http://data.cma.cn/, accessed on 9 August 2021).

Data Collection and Pre-Processing
Vegetation status over the study area is measured in terms of the Normalized Difference Vegetation Index (NDVI), which is a widely used proxy of vegetation canopy greenness. In this study, the moderate-resolution imaging spectroradiometer MODIS-derived monthly composite product (MOD13A3) with a spatial resolution of 1 km was selected to generate a long time series of NDVI dataset from 2001 to 2020, which was pre-processed and released by the Atmosphere Archive and Distribution System (LAADS) of the National Aeronautics and Space Administration (NASA) (https://ladsweb.modaps.eosdis.nasa.gov/, accessed on 9 August 2021).
Meteorological data from 2001 to 2020, including temperature, precipitation, evaporation, and average relative humidity, were collected from 3 meteorological stations, namely, Wulan Station, Chaka Station, and Dulan Station, and the data were available from China Meteorological Data Service Center (http://data.cma.cn/, accessed on 9 August 2021). Meteorological data from 2001 to 2020, including temperature, precipitation, evaporation, and average relative humidity, were collected from 3 meteorological stations, namely, Wulan Station, Chaka Station, and Dulan Station, and the data were available from China Meteorological Data Service Center (http://data.cma.cn/, accessed on 9 August 2021).
As we focused on vegetation change among annual growing seasons, we defined the yearly period of April-October as the plant growing season. The NDVI data were extracted from raw hierarchy format (HDF) images and converted into TIFF format with a new projection coordinate of WGS_1984_UTM_Zone_47N by using HDF-EOS to GeoTIFF Conversion Tool (HEG). The monthly NDVI raster images of the study area in the growing season from 2001 to 2020 were obtained by using the "Mask extraction" tool of ArcGIS combing with Python. The annual NDVI images were composited based on the monthly NDVI raster data in the growing season with the method of maximum value composite (MVC) [39,40]. Subsequently, the growing season NDVI (GSN) data of the county and all vegetation types were generated by calculating the mean value of pixel NDVI for the corresponding areas [41]. Considering that the spatial scale of the study area is a county, the climate in Wulan County was regarded as a single and uniform microclimate in this study. Growing season average temperature (GSAT), precipitation (GSP), evaporation (GSE), and average relative humidity (GSARH) annually from 2001 to 2020, as well as multi-year mean values of average temperature (MAT), precipitation (MP), evaporation (ME), and average relative humidity (MARH), were calculated for further analysis based on the monthly meteorological data calculated by Voronoi diagram method [36].

Study Method
The ordinary least squares method and Pearson correlation analysis were employed to explore the temporal dynamics of growing season NDVI (GSN) and climate factors and detect the relationships between them [30,[42][43][44][45]. Then, we further evaluate the direct and indirect effects of climate variables on vegetation GSN with the method of path analysis. Path analysis was proposed by American geneticist Sewall Wright in 1921. By decomposing the apparent direct correlation, the correlation coefficient between causal variables can be divided into direct effect (direct path coefficient) and indirect effect (indirect path coefficient), so as to study the data structure of causality and analyze the direct and indirect importance of independent variables to dependent variables. The principles of path analysis are as follows.
For an interrelated system, if there is a linear relationship between a dependent variable and independent variables, the regression equation is: where y is the dependent variable, i.e., GSN.; x m is the independent variables, i.e., climate factors, and m is the number of independent variables. The effect of each variable on the dependent variable can be characterized by the path coefficient, which can be expressed as: where r x i x j is the simple correlation coefficient of the variables x i and x j , i, j ≤ m; a i is the direct path coefficient between the independent variable x i and dependent variable y, that is, the partial correlation coefficient after the standardization of x i and y, indicating the direct effect of x i on y; r x i x j × a j is the indirect path coefficient that indicates the indirect effect of x i on dependent variable y through x j ; ∑ i =j r x i x j × a j represents the total indirect effect of x i on the dependent variable y through other variables; r x i y is the simple correlation coefficient between the independent variable x i and dependent variable y; a i × r x i y is the contribution rate of independent variable x i on dependent variable y.

Vegetation Dynamics Analyzed by Growing Season NDVI
The yearly changes in growing season NDVI (GSN) of the county and each vegetation type showed similar trends of wavelike increase with consistent fluctuation processes from 2001 to 2020, and their GSN values ranking from high to low were shrubland (0.385 ± 0.052), rest vegetation (0.365 ± 0.045), meadow (0.327 ± 0.032), steppe (0.258 ± 0.042), county (0.222 ± 0.029) and desert steppe (0.153 ± 0.020) (Figure 3). GSN of the whole county ranged from 0.151 (2001) to 0.264 (2012) and demonstrated an extremely significant increasing trend with an increasing rate of 0.0035 yr −1 over the study period. Similarly, GSN of different vegetation types also had extremely significant increasing trends and their increasing rates were shown in Table 1. from 2001 to 2020, and their GSN values ranking from high to low were shrubland (0.385 ± 0.052), rest vegetation (0.365 ± 0.045), meadow (0.327 ± 0.032), steppe (0.258 ± 0.042), county (0.222 ± 0.029) and desert steppe (0.153 ± 0.020) (Figure 3). GSN of the whole county ranged from 0.151 (2001) to 0.264 (2012) and demonstrated an extremely significant increasing trend with an increasing rate of 0.0035 yr −1 over the study period. Similarly, GSN of different vegetation types also had extremely significant increasing trends and their increasing rates were shown in Table 1.  The multi-year mean values of GSN, ranging from 0.027 to 0.717, had obvious spatial heterogeneity and gradually decreased from northeast to southwest in the study area (Figure 4). The low values of GSN between 0 and 0.1 accounted for 17.0% of the total area and were mainly distributed in the western desert region; GSN between 0.1 and 0.2 had the highest proportion, about 40.0% of the total area, and the primary vegetation type was the desert steppe, which mainly distributed in the middle and southeast of the county; GSN over 0.2 were mainly distributed in the northeast of the county with a proportion of 42.6% of the total area. GSN changes in individual pixels over the study period were fitted by the method of ordinary least squares ( Figure 5). The fitting slope ranged from −0.016 to 0.032, and decreased gradually from northeast to southwest as a whole; the area with GSN  The multi-year mean values of GSN, ranging from 0.027 to 0.717, had obvious spatial heterogeneity and gradually decreased from northeast to southwest in the study area ( Figure 4). The low values of GSN between 0 and 0.1 accounted for 17.0% of the total area and were mainly distributed in the western desert region; GSN between 0.1 and 0.2 had the highest proportion, about 40.0% of the total area, and the primary vegetation type was the desert steppe, which mainly distributed in the middle and southeast of the county; GSN over 0.2 were mainly distributed in the northeast of the county with a proportion of 42.6% of the total area. GSN changes in individual pixels over the study period were fitted by the method of ordinary least squares ( Figure 5). The fitting slope ranged from −0.016 to 0.032, and decreased gradually from northeast to southwest as a whole; the area with GSN fitting slope greater than 0 accounted for 97.2% of the total, indicating that the vegetation covers of almost the whole county were in the process of improvement, even in the desert areas at the northwest of the county. The vegetation dynamics were divided into 6 changing patterns based on the fitting slope and significance ( Figure 6), namely, Extremely significant increase (ESI, Slope > 0 and p < 0.01), Significant increase (SI, Slope > 0 and 0.01 ≤ p < 0.05), Non-significant increase (NSI, Slope > 0 and p ≥ 0.05), Non-significant decrease (NDI, Slope < 0 and p ≥ 0.05), Significant decrease (SD, Slope < 0 and 0.01 ≤ p < 0.05), Extremely significant decrease (ESD, Slope < 0 and p < 0.01). Approximately 53.9% of the total area showed ESI and 22.5% was SI from 2001 to 2020 and their corresponding pixels were evenly distributed throughout the county; for different vegetation types, about 80% of shrubland, steppe, and desert steppe were ESI and SI, while meadow and rest vegetation were about 50% and 60%, respectively. nificant decrease (NDI, < 0 and ≥ 0.05 ), Significant decrease (SD, < 0 and 0.01 ≤ < 0.05 ), Extremely significant decrease (ESD, < 0 and < 0.01 ).
Approximately 53.9% of the total area showed ESI and 22.5% was SI from 2001 to 2020 and their corresponding pixels were evenly distributed throughout the county; for different vegetation types, about 80% of shrubland, steppe, and desert steppe were ESI and SI, while meadow and rest vegetation were about 50% and 60%, respectively.   nificant decrease (NDI, < 0 and ≥ 0.05 ), Significant decrease (SD, < 0 and 0.01 ≤ < 0.05 ), Extremely significant decrease (ESD, < 0 and < 0.01 ).
Approximately 53.9% of the total area showed ESI and 22.5% was SI from 2001 to 2020 and their corresponding pixels were evenly distributed throughout the county; for different vegetation types, about 80% of shrubland, steppe, and desert steppe were ESI and SI, while meadow and rest vegetation were about 50% and 60%, respectively.

Characteristics of Climatic Factors
Monthly variations of meteorological factors in the study area showed that temperature, precipitation and relative humidity reached the peak in July, while evaporation reached the largest in April; over the vegetation growing season, average temperature was 2.7 times of annual average temperature, precipitation accounted for 94.0% of annual precipitation, evaporation was 75.8% of annual evaporation, and average relative humidity was significantly higher than that in non-growing season (Figure 7).

Characteristics of Climatic Factors
Monthly variations of meteorological factors in the study area showed that temperature, precipitation and relative humidity reached the peak in July, while evaporation reached the largest in April; over the vegetation growing season, average temperature was 2.7 times of annual average temperature, precipitation accounted for 94.0% of annual precipitation, evaporation was 75.8% of annual evaporation, and average relative humidity was significantly higher than that in non-growing season (Figure 7). The yearly changes in climate factors over the growing season showed different variation processes (Figure 8). Growing season average temperature (GSAT) ranged from Figure 6. Changing patterns of the pixel growing season NDVI with years from 2001 to 2020. CT, SH, ST, ME, DS, and RV represent county, shrubland, steppe, meadow, desert steppe, and rest vegetation, respectively. ESI, SI, NSI, NDI, SD, and ESD represent Extremely significant increase, Significant increase, Non-significant increase, Non-significant decrease, Significant decrease, and Extremely significant decrease, respectively. Figure 6. Changing patterns of the pixel growing season NDVI with years from 2001 to 2020. CT, SH, ST, ME, DS, and RV represent county, shrubland, steppe, meadow, desert steppe, and rest vegetation, respectively. ESI, SI, NSI, NDI, SD, and ESD represent Extremely significant increase, Significant increase, Non-significant increase, Non-significant decrease, Significant decrease, and Extremely significant decrease, respectively.

Characteristics of Climatic Factors
Monthly variations of meteorological factors in the study area showed that temperature, precipitation and relative humidity reached the peak in July, while evaporation reached the largest in April; over the vegetation growing season, average temperature was 2.7 times of annual average temperature, precipitation accounted for 94.0% of annual precipitation, evaporation was 75.8% of annual evaporation, and average relative humidity was significantly higher than that in non-growing season (Figure 7).  The yearly changes in climate factors over the growing season showed different variation processes (Figure 8). Growing season average temperature (GSAT) ranged from 10.0 • C to 11.6 • C with a mean value of 10.5 • C, the minimum and maximum value appeared in 2020 and 2016, respectively; growing season precipitation (GSP) ranged from 118.9 to 299.6 mm with a mean value of 220.9 mm, and the minimum and maximum value appeared in 2001 and 2018, respectively; growing season evaporation (GSE) was in the ranges of 10.0 °C to 11.6 °C with a mean value of 10.5 °C, the minimum and maximum value appeared in 2020 and 2016, respectively; growing season precipitation (GSP) ranged from 118.9 to 299.6 mm with a mean value of 220.9 mm, and the minimum and maximum value appeared in 2001 and 2018, respectively; growing season evaporation (GSE) was in the ranges of 944.2-1697.4 mm with a mean value of 1183.2 mm, and the minimum and maximum value appeared in 2019 and 2001 respectively; growing season average relative humidity (GSARH) was between 42.1% and 50.0% with a mean value of 46.1%, and the minimum and maximum value appeared in 2001 and 2019, respectively. The general linear regression analysis of the 4 climate factors revealed that the changing processes of GSAT and GSP had insignificant increase or decrease trends; GSE decreased extremely significantly with the year having a reduction rate of 29.6 mm yr −1 ( = 0.780, = 0), and the GSE in 2020 had decreased by 786 mm compared with 2001; GSARH showed a significant increasing trend with an increase rate of 0.16% yr −1 ( = 0.237, = 0.030).

Relationships between GSN Change and Climate Factors
The Pearson's correlation analysis illustrated high correlations among GSP, GSE, and GSARH, while low correlations between GSAT and other 3 climate factors; furthermore, GSN of the county and 5 vegetation types exhibited evident relationships with GSP, GSE, and GSARH (Figure 9). GSP showed a significant negative correlation with GSE ( = −0.51, < 0.05), while a significant positive correlation with GSARH ( = 0.54, < 0.05); and there was an extremely significant negative correlation between GSE and GSARH ( = −0.58, < 0.01). This indicated that the climatic factors were mutually affected. As the above result that GSN of the county and different vegetation types had similar fluctuation process with years, GSN of the county and different vegetation types also had similar relationships with climate factors. GSN trends of the county and the 5 vegetation types had no obvious correlation with GSAT, extremely significant negative correlation with GSE (−0.89 < < −0.81, < 0.01), and extremely significant positive correlation with GSARH ( 0.71 < < 0.76, < 0.01 ). For GSP, GSN trends of the county, shrubland, steppe, meadow, and desert steppe exhibited an extremely significant positive correlation with it (0.59 < < 0.61, < 0.01), while that of rest vegetation was a significant positive correlation with it ( = 0.56, < 0.05).

Relationships between GSN Change and Climate Factors
The Pearson's correlation analysis illustrated high correlations among GSP, GSE, and GSARH, while low correlations between GSAT and other 3 climate factors; furthermore, GSN of the county and 5 vegetation types exhibited evident relationships with GSP, GSE, and GSARH ( Figure 9). GSP showed a significant negative correlation with GSE (r = −0.51, p < 0.05), while a significant positive correlation with GSARH (r = 0.54, p < 0.05); and there was an extremely significant negative correlation between GSE and GSARH (r = −0.58, p < 0.01). This indicated that the climatic factors were mutually affected. As the above result that GSN of the county and different vegetation types had similar fluctuation process with years, GSN of the county and different vegetation types also had similar relationships with climate factors. GSN trends of the county and the 5 vegetation types had no obvious correlation with GSAT, extremely significant negative correlation with GSE (−0.89 < r < −0.81, p < 0.01), and extremely significant positive correlation with GSARH (0.71 < r < 0.76, p < 0.01). For GSP, GSN trends of the county, shrubland, steppe, meadow, and desert steppe exhibited an extremely significant positive correlation with it (0.59 < r < 0.61, p < 0.01), while that of rest vegetation was a significant positive correlation with it (r = 0.56, p < 0.05).
According to the Path Analysis, the direct and indirect effects of climate factors on vegetation dynamics were revealed (Table 2). For the county and 5 vegetation types, GSE had the maximum absolute values of direct path coefficients on GSN, ranging from 0.559 to 0.663, which were significantly higher than its absolute values of indirect path coefficients through other climate factors; moreover, the contribution rates of GSE to vegetation dynamics were also the largest, ranging from 42.4% to 58.8%; it can be seen that GSE was the primary climate factor affecting vegetation dynamics, and its effect on desert steppe was the most prominent. GSARH had the sub-maximum absolute values of direct path coefficients (0.197-0.329) on GSN but the maximum absolute values of indirect path coefficients through GSE (0.303-0.382), and the contribution rates of GSARH ranged from 13.9% to 25.0% and ranked second, that is, GSARH was the second important climate factor affecting vegetation dynamics and its effects was mainly generated through GSE. The effect of GSP on GSN was mainly through the indirect effect of GSE, because the absolute values of its indirect path coefficients through GSE ranked second with the range of 0.268-0.338, however, the contribution rate of GSP on vegetation dynamics were only 3.1%-10.2%. GSAT had low absolute values of direct and indirect path coefficients as well as contribution rates, indicating that the direct and indirect effects of GSAT on vegetation dynamics both were weak.
Atmosphere 2022, 13, x FOR PEER REVIEW 10 of 16 Figure 9. Correlation analysis between growing season NDVI (GSN) and 4 climatic factors over growing season. CT, SH, ST, ME, DS and RV represent county, shrubland, steppe, meadow, desert steppe and rest vegetation, respectively. The symbols below the diagonal were the correlation coefficient, and the symbols above the diagonal were the corresponding significance. * and ** represented a significant relationship at = 0.05, = 0.01 level, respectively.
According to the Path Analysis, the direct and indirect effects of climate factors on vegetation dynamics were revealed ( Table 2). For the county and 5 vegetation types, GSE had the maximum absolute values of direct path coefficients on GSN, ranging from 0.559 to 0.663, which were significantly higher than its absolute values of indirect path coefficients through other climate factors; moreover, the contribution rates of GSE to vegetation dynamics were also the largest, ranging from 42.4% to 58.8%; it can be seen that GSE was the primary climate factor affecting vegetation dynamics, and its effect on desert steppe was the most prominent. GSARH had the sub-maximum absolute values of direct path coefficients (0.197-0.329) on GSN but the maximum absolute values of indirect path coefficients through GSE (0.303-0.382), and the contribution rates of GSARH ranged from 13.9% to 25.0% and ranked second, that is, GSARH was the second important climate factor affecting vegetation dynamics and its effects was mainly generated through GSE. The effect of GSP on GSN was mainly through the indirect effect of GSE, because the absolute values of its indirect path coefficients through GSE ranked second with the range of 0.268-0.338, however, the contribution rate of GSP on vegetation dynamics were only 3.1%- Figure 9. Correlation analysis between growing season NDVI (GSN) and 4 climatic factors over growing season. CT, SH, ST, ME, DS and RV represent county, shrubland, steppe, meadow, desert steppe and rest vegetation, respectively. The symbols below the diagonal were the correlation coefficient, and the symbols above the diagonal were the corresponding significance. * and ** represented a significant relationship at p = 0.05, p = 0.01 level, respectively.

Discussion
The Normalized Difference Vegetation Index (NDVI) in growing season is widely used to explore the vegetation dynamics of grasslands [24][25][26][27]46]. In this study, the annual mean values of the growing season NDVI (GSN) in Wulan County are obviously different from previous studies at large spatial scales. The multi-year mean values of GSN in Wulan County are significantly lower than the results obtained by Dai et al. [47] in Qinghai Province from 2000 to 2015 (0.37-0.42), Huang et al. [31] in Qinghai Plateau from 2001 to 2017 (0.48-0.51), and Zhang et al. [48] in northern Tibet (0.112-0.492). The mean reasons are that Wulan County is located at the eastern edge of Qaidam Basin Desert, where the soil texture is mostly sandy soil and the vegetation is sparse. Desert grassland, the dominant grassland type, with low GSN accounts for more than half of the area of the county, resulting in the obviously low vegetation GSN of the whole county. The values reported by Bai et al. [30] were in the range of 0.232-0.266 from 1982 to 2015 in the Three-River Headwater Region, close to the results of this study, but the minimum value was significantly higher than that in Wulan County. On the one hand, due to the difference in vegetation types, Wulan County is dominated by desert steppe with low GSN that governs the lowest level of GSN. On the other hand, the research period in Bai et al. [30] is longer and spans two centuries. The increase of population and livestock led to over reclamation and overgrazing, resulting in the expansion of grassland degradation and low GSN in the Three-River Headwater Region during 1982-2000, then gradually improved with the implementation of a series of ecological recovery programs after 2000 [49][50][51].
For temporal trends, the overall increasing trend of GSN in Wulan County is consistent with the above reports. However, the increasing rate of GSN in Wulan County was significantly higher than those for Qinghai province [47], Northern Tibet [48], and Qinghai-Tibet Plateau [31,52], but a little lower than that for Three-River Headwater Region [30]. The scale of the study area is related to the complexity of regional vegetation community composition, the determination of the duration of vegetation growing season, the difference of climate and environmental conditions, etc. Thus, it had a nonnegligible impact on the regional GSN level and its changing rate in alpine grassland. Of course, there were some other impacting factors, such as the study period, the type and resolution of remote sensing data, the data processing methods as well as the implementation of a series of ecological restoration projects. With regard to the spatial distribution of GSN in Wulan County, the spatial heterogeneity was mainly explained by the distribution of vegetation types. The shrubland, steppe, meadow and rest vegetation with high values of GSN were mostly distributed in the northeast and a small amount in the southeast and middle of the county, while the desert steppe with low values of GSN was mainly distributed in the Midwest of the county and gradually sparse with the decrease of latitude. In addition, 97.2% of GSN showed an increasing trend, and 76.4% exhibited a significantly increasing trend with even distribution throughout the county. This indicated that the improvement of vegetation coverage was obviously better than that in Qinghai Province [47] and Qinghai-Tibet Plateau [31,52].
The area of Wulan County accounts for only 0.52% and 1.85% of the Qinghai-Tibet Plateau and Qinghai province, respectively [53,54]. The climate in Wulan County was relatively uniform and unique and was inconsistent with those reported by several previous studies at larger spatial scales [30,53,[55][56][57]. In particular, the change of temperature is obviously in disagreement with the overall rising trend of temperature in the Qinghai-Tibet Plateau [58,59]. These results suggest that relatively small-scale areas may undergo significant changes in climate [30,60]. As for the effects of multiple climatic factors on vegetation growth, our results showed that growing season evaporation (GSE) and growing season average Relative Humidity (GSARH) could better explain the vegetation dynamics than growing season average temperature (GSAT) and growing season precipitation (GSP) at small scale. Vegetation dynamics had an extremely significant negative correlation with GSE, while had extremely significant positive correlations with GSARH and GSP. The path analysis indicated that GSE was the direct factor affecting the vegetation growth, while GSARH and GSP were the indirect factors mainly through GSE. Wulan County is located in the arid climate area in the northeast of the Qinghai-Tibet Plateau, the evaporation is much higher than the precipitation, and the vegetation growth mainly depends on the available water [51,61]. When the GSP and GSAT remain relatively stable, the continuous decline of GSE slows down the dissipation of precipitation and soil moisture and increases the average relative humidity, indicating that the climate in Wulan County is developing towards warm and humid [58,62]. It is supposed that decreasing GSE will increase alpine grassland GSN due to easing water constraints and improving water use efficiency on vegetation growth [2,63]. Describing it differently, the increase of GSARH and GSP could limit GSE and form advantageous climatic conditions for the restoration of alpine grassland. Moreover, the total contribution rates of GSE, GSARH and GSP to vegetation dynamics was only about 78.0%, indicating that the factors affecting vegetation dynamics had not been fully considered, such as other climate factors, land use, grazing, soil and water conservation, and so on, which need to be further strengthened in future studies.

Conclusions
Vegetation is the key index of alpine grassland in the Qinghai-Tibet Plateau. Its condition can reduce or aggravate climate disasters and determine the regional livestock carrying capacity. Taking Wulan County as the example, this study quantitatively analyzed the vegetation dynamics of alpine grassland on the Qinghai-Tibet Plateau and the relationships with multiple climate factors at small scale on the basis of Normalized Difference Vegetation Index (NDVI) in growing season. The results deepened our understanding of the driving mechanism of climate change on the vegetation dynamics of alpine grassland in the Qinghai-Tibet Plateau. (1) The annual growing season NDVI of the study area ranged between 0.151 to 0.264 from 2001 to 2020, and extremely significant wavelike increased at a rate of 0.0035 yr −1 . For the 5 vegetation types, the annual growing season NDVI ranking from high to low were shrubland, rest vegetation, meadow, steppe, and desert steppe, and showed similar increasing trends with the whole area at rates of 0.0058 yr −1 , 0.0049 yr −1 , 0.0034 yr −1 , 0.0049 yr −1 , and 0.0026 yr −1 , respectively. (2) Over the study area, growing season NDVI gradually decreased from northeast to southwest with ranging values of 0.027-0.717. The fitting slopes of growing season NDVI in individual pixels from 2001 to 2020 were between −0.016 and 0.032, and the positive slopes account for 97.2% of the county area, and growing season NDVI increased significantly in about 76.4% of the county area. (3) During the study period, growing season evaporation decreased extremely significantly at a rate of 29.6 mm yr −1 , growing season average relative humidity increased significantly at a rate of 0.16% yr −1 , while growing season average temperature and growing season precipitation were relatively stable. (4) For the whole study area and 5 vegetations types, growing season evaporation had an extremely significant negative correlation with growing season NDVI, and was the primary climate factor that had a direct effect on vegetation dynamics with a contribution rate of 42.4%-58.8%; growing season average relative humidity had an extremely significant positive correlation with growing season NDVI, and had an indirect effect on vegetation dynamics through growing season evaporation with a contribution rate of 13.9%-25.0%; growing season precipitation had extremely significant positive correlation with growing season NDVI, and also had an indirect effect on vegetation dynamics through growing season evaporation, but its contribution rates were less than 11.0%; growing season average temperature had no significant effect on vegetation dynamics.