Spatiotemporal Variations in Grassland Vulnerability on the Qinghai-Tibet Plateau Based on a Comprehensive Framework

: Grasslands are globally important for providing essential ecosystem services and maintaining ecological security. Monitoring and assessing grassland vulnerability are critical for developing long-term grassland management policies and strategies. The grassland vulnerability on the Qinghai-Tibet Plateau (QTP) is considered high, but its spatial and temporal variations in response to human activities and climate change are not well understood. In this study, a comprehensive grassland vulnerability index (GVI), which includes natural factors (VNF), environmental disturbances (VED), and socioeconomic impacts (VSI), was developed by using the analytic hierarchy process (AHP), principal component analysis (PCA), and environmental vulnerability distance index (EVDI). Our results showed that the spatial distribution of GVI had obvious heterogeneity, decreasing from northwest to southeast; the regions with serious and extreme vulnerability were mainly concentrated in the north-western alpine steppe and desert steppe. From 2000 to 2018, GVI decreased from 0.61 in 2000 to 0.60 in 2010 and then to 0.59 in 2018, demonstrating a healthy tendency. The normalized difference vegetation index (NDVI), land desertiﬁcation, and population were the factors that had the most signiﬁcant impact on VNF, VED, and VSI, respectively. The global Moran’s I index of grassland vulnerability was greater than 0, with a signiﬁcant positive spatial correlation. The number of High-High and Low-Low units decreased, indicating that the High-High and Low-Low cluster regions tended to be discrete. Moreover, our results suggest that understanding the variations in grassland vulnerability on the QTP is important for regional sustainable development in the context of intensiﬁed climate change and human disturbances.


Introduction
The assessment of ecosystem vulnerability, which focuses on the potential consequences of disruptions, has become a key issue in global change ecology and sustainability studies [1][2][3]. Climate change presents numerous challenges to global ecosystems, including ecosystem degradation, biodiversity loss, and land desertification, all of which are expected to become more frequent and intense in the future, especially for the sensitive and vulnerable grassland ecosystem [4,5]. The grassland ecosystem provides feed for grazing livestock as well as a food supply and a source of income for over 800 million people [6]. It also provides numerous key ecosystem functions (e.g., absorbing the carbon dioxide (CO 2 ) stored underground, boosting the carbon sink capacity of ecosystems, and providing the rich animal and plant gene bank) [4,7,8]. In summary, the assessment of grassland vulnerability is vital and can be used as a key factor in possible ecological issues [9,10].

Study Area
The Qinghai-Tibet Plateau, covering about one-fourth of the total area of China (2.5 × 10 6 km 2 ), is called the Third Pole. With an average elevation of more than 4000 metres, it is geomorphologically the largest and highest plateau on Earth [31]. It lies from 26 • to 39 • north latitude and 73 • to 104 • east longitude, covering Tibet and the Xinjiang Uygur Autonomous Regions and the Qinghai, Gansu, Sichuan, and Yunnan Provinces. Under the influence of the Asian monsoon, the plateau is also the source of nine large rivers, supplying billions of people downstream with freshwater, adequate food, and other ecosystem services [35,36]. Meanwhile, because of the unique and complex interactions between the atmosphere, cryosphere, hydrology, and geology, the region is highly vulnerable to environmental changes [33]. In the last century, environmental conditions have changed dramatically. The air temperature on the QTP has warmed dramatically faster than the global average [35,37]. These changes may cause glacier retreat [38,39], snow cover change [33,40], and permafrost melting [41], and drive environmental changes at local, regional, and even global scales.
Grassland is the dominant vegetation type across the QTP, which covers approximately 15.3 × 10 5 km 2 and accounts for 59.15% of the total QTP area. It is one of the most important pasture regions in China and even in Asia [42]. It is also vital for livestock husbandry and environmental security [1]. From southeast to northwest, given the influence of plateau climate changes, the grassland types change from alpine meadow to alpine steppe to desert steppe ( Figure 1).

Study Area
The Qinghai-Tibet Plateau, covering about one-fourth of the total area of China (2.5 × 10 6 km 2 ), is called the Third Pole. With an average elevation of more than 4000 metres, it is geomorphologically the largest and highest plateau on Earth [31]. It lies from 26° to 39° north latitude and 73° to 104° east longitude, covering Tibet and the Xinjiang Uygur Autonomous Regions and the Qinghai, Gansu, Sichuan, and Yunnan Provinces. Under the influence of the Asian monsoon, the plateau is also the source of nine large rivers, supplying billions of people downstream with freshwater, adequate food, and other ecosystem services [35,36]. Meanwhile, because of the unique and complex interactions between the atmosphere, cryosphere, hydrology, and geology, the region is highly vulnerable to environmental changes [33]. In the last century, environmental conditions have changed dramatically. The air temperature on the QTP has warmed dramatically faster than the global average [35,37]. These changes may cause glacier retreat [38,39], snow cover change [33,40], and permafrost melting [41], and drive environmental changes at local, regional, and even global scales.
Grassland is the dominant vegetation type across the QTP, which covers approximately 15.3 × 10 5 km 2 and accounts for 59.15% of the total QTP area. It is one of the most important pasture regions in China and even in Asia [42]. It is also vital for livestock husbandry and environmental security [1]. From southeast to northwest, given the influence of plateau climate changes, the grassland types change from alpine meadow to alpine steppe to desert steppe ( Figure 1).

Framework for GVI Assessment
Ecosystem vulnerability assessment is a complicated procedure that focuses on the relationship between environmental conditions and human interference [15,28]. Based on comprehensive analyses of the vulnerability of grassland, we accepted a novel framework presented by Sun et al. [28], which divided vulnerability into three subsystems, namely, vulnerability by natural factors (VNF), vulnerability by environmental disturbances (VED), and vulnerability by socioeconomic impacts (VSI), and was successfully used in assessing forest vulnerability on the QTP. Based on the characteristics of grasslands and the novel framework, we selected 13 indicators, considering the operability and availability, and built a framework of VNF-VED-VSI (Table 1). Natural factors can directly reflect ecological and environmental changes, which are vital for vulnerability evaluation [28]. Temperature and precipitation can impact the biodiversity and quality of the environment. Differences in elevation and slope change the topographical features. Indicators, including annual average temperature, annual average precipitation, elevation, slope, soil organic matter content, and the normalized difference vegetation index (NDVI), were chosen to characterize vulnerability by natural factors.
The QTP is vulnerable to environmental disturbances and natural disasters. The occurrence of freeze-thaw erosion, land desertification, landslides, and grazing pressure causes serious and direct damage to grasslands. Therefore, it is an important foundation for assessing the risk of grassland vulnerability [28].
Vulnerability by socioeconomic impact indicators is selected to quantify anthropic pressures across the QTP because human activities significantly change the evolution of environmental characteristics [28,34]. With the population increase and rapid economic development, socioeconomic activities (e.g., grazing and land cover changes) modify the vulnerability and increase such uncertainties [30]. Population, gross domestic product (GDP), and actual livestock are selected. Population and GDP can directly reflect human disturbance and economic activities, and actual livestock can reflect the pressure of the livestock industry on the grassland ecosystem.

Data Collection and Approaches
In this paper, the data used mainly included remote sensing, socioeconomic, and meteorological data. Meteorological datasets were obtained from the National Earth System Science Data Centre, National Science and Technology Infrastructure of China (http://www.geodata.cn/ (accessed on 10 August 2021)) and consisted of the annual average temperature data and annual average precipitation data in 2000, 2010, and 2018, with a spatial resolution of 1000 m. A digital elevation model (DEM) and slope data were derived from the sharing platform Geospatial Data Cloud (http://www.gscloud.cn/ (accessed on 26 August 2021)), with a spatial resolution of 90 m. Soil organic matter content was provided by the Harmonized World Soil Database (HWSD), Food and Agriculture Organization of the United Nations (https://www.fao.org/ (accessed on 10 August 2021)), with a spatial resolution of 1000 m. Grazing pressure and actual livestock datasets for 2000, 2010, and 2018 were obtained from the Qinghai-Tibet Plateau ecosystem data integration platform (http://www.gscloud.cn/qtp/qtpecosys/ (accessed on 2 September 2021)) with a spatial resolution of 1000 m. GDP, population, and NDVI data in 2000, 2010, and 2018 were obtained from the Resources and Environmental Sciences and Data Centre (http: //www.resdc.cn/ (accessed on 21 August 2021)), with a spatial resolution of 1000 m (population and GDP in 2015 were adopted as assessments in 2018).
Multiple factors were selected for the freeze-thaw erosion vulnerability assessment, including the annual temperature range, annual average precipitation, slope, aspect, and vegetation coverage. Based on previous research, we built an evaluation index system and weighted assignment of different factors using the weighted average method to calculate the impact of freeze-thaw erosion on grassland [43]. In the present study, the vulnerability of landslides was calculated by analysing hazard factors (population, slope, elevation, distance to roads, and distance to streams). Following the method of previous studies with some additions, we quantified the index and assigned weight values, forming a map of landslide vulnerability with ArcGIS 10.2 [28,44]. The NDVI is an index that can accurately assess the status, change, and trends of desertification [45]. Therefore, we used NDVI to reflect the vulnerability of land desertification based on previous research [17,45].
In light of the relationship between grassland vulnerability and evaluation factors, indicators were classified as positive or negative indicators. When the grassland vulnerability index values increased, the positively correlated indicators increased and the negatively correlated indicators decreased [15]. To reduce the influence of differences among various factors, we used the maximum-minimum standardization method so that the indicator values ranged between 0 and 1 [13,46]. The indicators were standardized using Equations (1) and (2): Negative indicator : where X i represents the normalized value of i in the grid cell; is the actual value of the index; and X min and X max are the minimum and maximum values of the index in the study region. All indicators were pre-processed in ArcGIS 10.2, including merging, clipping, reprojection, and resampling. All indicators were converted uniformly to raster with a spatial resolution of 1000 m.

Data Analysis
Combining the AHP-PCA-EVDI, autocorrelation between index factors was effectively solved, and the accuracy of the evaluation results was apparently improved. Accounting for prior experience and expert knowledge, AHP constructed a judging matrix and assigned the weights subjectively [15]. To reduce the dimensionality of a dataset, PCA extracted indicators of significant difference, and the related indicators were objectively transformed into a set [47]. In this study, indicators of VNF, VED, and VSI were calculated using an entropy weight model, integrating AHP with PCA. In addition, the EVDI described the minimum vulnerability point as a determinate point; the distances between each of the other points and the minimum point were calculated in a multidimensional space [22]. As a result, the GVI was calculated using the distances from VNF, VED, and VSI to their minimum point. The EVDI allowed a reasonable evaluation procedure and made each subsystem effective. The comprehensive method united subjective and objective approaches, which made the results reasonable and convincing.
Step 1: Use the AHP to determine the weights of indictors of VNF, VED, and VSI. Classify the importance level of each factor, assign a value, and form a judging matrix according to the pairwise comparisons. Calculate the consistency ratio (CR) and the weight of each indicator (W 1j ) using Equations (3) and (4): where CI is the consistency index; λ max is the maximum eigenvalue of the judging matrix; and n is the order of the judging matrix and the number of judgement factors.
where CR is the consistency ratio; CI is the consistency index; and RI is the random index.
In this study, all of the values of CR in VNF, VED, and VSI are 0.000, suggesting that the judgement is acceptable (Tables 2-4).   Step 2: Calculate the weights for the PCA using ArcGIS 12.0 spatial analysis. Obtain the contribution rate and cumulative contribution rate of each principal component. Apply the Equations (5) and (6) to calculate the weight of indicators (W 2j ): where H j is the variance of each factor; W 2j is the PCA weight of factors; j is the number of indicators; and k is the number of principal components.
Step 3: Calculate the subjective and objective weights using AHP and PCA for VNF, VED, and VSI. Apply the principle of minimum information entropy to calculate the combined weight. Then, calculate the comprehensive weight of VNF, VED, and VSI using Equation (7): where W j is the comprehensive weight; W 1j is the subjective weight; W 2j is the objective weight; and j is the number of indicators.
Step 4: Based on the results of VNF, VED, and VSI, the EVDI was used to assess grassland vulnerability. The smaller the distance index is, the better the ecological conditions of grassland. Calculate the distance index using Equation (8): where GV I represents the grassland vulnerability index; V NF, VED, and VSI represent the vulnerability by natural factors, vulnerability by environmental disturbance factors, and vulnerability by socioeconomic impact factors, respectively; and V NF min , VED min , and VSI min represent the minimum VNF, VED, and VSI, respectively.

Classification Method
In ArcGIS 12.0, based on Natural Breaks (Jenks), we divided grassland vulnerability into five levels of VNF, VED, VSI, and GVI: normal, slight vulnerability, moderate vulnerability, serious vulnerability, and extreme vulnerability (Table S1). An increase in the level of vulnerability was considered a significantly increased vulnerability, while a decrease in the level of vulnerability was considered a significantly decreased vulnerability.

Correlation Analysis
The spatial autocorrelation indicator was useful for determining whether the values of attributes were strongly associated with the values of adjacent spatial attributes [48]. In this paper, Moran's I index (global spatial autocorrelation) and the local indicator of spatial association (LISA) were used to analyse the spatial correlation of the GVI.
(1) Global Moran's I The correlation of attribute values for the spatially adjacent unit was shown by Moran's I coefficient. When Moran's I had an absolute value close to 1, the unit's spatial autocorrelation was strong. Moreover, when Moran's I had an absolute value close to 0, the unit's spatial autocorrelation was weak. Spatial correlation can be calculated using Equation (9): where x i , x j is the attribute value of the location of i, j; n is the total number of grids; ω ij is the weight of the matrix; and x is the average value of the GVI.
(2) Local Moran's I The local spatial autocorrelation index analyses the effectiveness of Moran's I at each spatial unit, effectively representing local spatial connections. Equation (10) is applied: where x i , x j is the attribute value of the location of i, j; ω ij is the weight of the matrix; x is the average value of the GVI; and S is the sum of the spatial weight matrix.

Spatiotemporal Variations in VNF, VED, and VSI
The spatial and temporal variations in VNF are shown in Figure 2. Western and northern grasslands, which were dominated by serious and extreme vulnerability, were more vulnerable than eastern and southern grasslands. NDVI, annual average temperature, annual average precipitation, and elevation were the main influencing factors of VNF ( Table 2

Spatiotemporal Variations in Grassland Vulnerability
The results for GVI are shown in Figure 5. In general, the trend of GVI decreased from northwest to southeast. Areas with normal, slight vulnerability, and moderate vulnerability were mostly found in the eastern and southern QTP. As a result of flat terrain, low elevation, high vegetation cover, suitable climate, and fertile soil, the main grassland type was alpine meadow, which was relatively stable. Areas with serious and extreme vulnerability were mostly found in the western and northern QTP, including Qinghai Province and Tibet, especially near the mountains. The main grassland types were alpine steppe and desert steppe, which had low precipitation, steep terrain, and high elevation. Because of these characteristics, the ecological environment was extremely fragile.

Spatiotemporal Variations in Grassland Vulnerability
The results for GVI are shown in Figure 5. In general, the trend of GVI decreased from northwest to southeast. Areas with normal, slight vulnerability, and moderate vulnerability were mostly found in the eastern and southern QTP. As a result of flat terrain, low elevation, high vegetation cover, suitable climate, and fertile soil, the main grassland type was alpine meadow, which was relatively stable. Areas with serious and extreme vulnerability were mostly found in the western and northern QTP, including Qinghai Province and Tibet, especially near the mountains. The main grassland types were alpine steppe and desert steppe, which had low precipitation, steep terrain, and high elevation. Because of these characteristics, the ecological environment was extremely fragile.  Table  5). The heterogeneity of grassland vulnerability across the entire area was indicated by the spatial distribution of the GVI rate of change. Areas with significantly decreasing vulnerability, which were largely located in the eastern and southern grasslands, accounted for 13.75% of the total study area. A few places with considerably elevated vulnerability existed, and these regions accounted for 2.60% of the total research area.    Table 5). The heterogeneity of grassland vulnerability across the entire area was indicated by the spatial distribution of the GVI rate of change. Areas with significantly decreasing vulnerability, which were largely located in the eastern and southern grasslands, accounted for 13.75% of the total study area. A few places with considerably elevated vulnerability existed, and these regions accounted for 2.60% of the total research area.

Spatial Autocorrelation Characteristics of GVI
In ArcGIS 12.0, we created hexagons with a side length of 10 km. For continually partitioning a two-dimensional space, hexagons are regarded as more efficient than square grids, but rectangular or square grids are commonly used [49]. Hexagonal forms provide an isotropic neighbourhood, which allows for more accurate neighbourhood analysis [49]. The global Moran's I value of GVI in 2000, 2010, and 2018 was 0.565, 0.587, and 0.598, respectively. These findings suggested that grassland vulnerability featured clustering. In addition, the clustering was even more obvious from 2000 to 2018. The p values were 0.001, and the Z scores were greater than 2.58, passing a 99.9% confidence test, indicating that the spatial autocorrelation of grassland vulnerability was extremely significant (Table S2).
The correlation of adjacent regions was measured by local spatial autocorrelation, which was generally based on LISA. The grassland of QTP was divided into five regions: (1) High-High cluster (H-H) regions, in which both the region and its adjacent areas had high grassland vulnerability, mainly including serious and extreme vulnerability; (2) Low-Low cluster (L-L) regions, in which both the region and its adjacent areas had low grassland vulnerability, mainly including normal, slight vulnerability, and moderate vulnerability; (3) Low-High outlier (L-H) regions, in which the grassland vulnerability of the region was low but the grassland vulnerability of the adjacent region was high; (4) High-Low outlier (H-L) regions, in which the grassland vulnerability of the region was high but the grassland vulnerability of the adjacent region was low; and (5) Not significant regions, in which there is no significant difference between the region and its adjacent areas [50]. Positive spatial correlations existed in the H-H and L-L regions, whereas negative correlations existed in the H-L and L-H regions. From 2000 to 2018, the spatial autocorrelation produced similar results ( Figure 6). The H-H regions were mainly found in the western and northern QTP, mainly including alpine steppe and desert steppe. The number of units decreased from 1020 in 2000 to 845 in 2018, indicating that the spatial distribution tended to be discrete. L-L regions were mainly found in the eastern and southern QTP, mainly including alpine meadows. The number of units decreased from 809 in 2000 to 778 in 2018, indicating that the spatial distribution tended to be discrete. number of units decreased from 1020 in 2000 to 845 in 2018, indicating that the spatial distribution tended to be discrete. L-L regions were mainly found in the eastern and southern QTP, mainly including alpine meadows. The number of units decreased from 809 in 2000 to 778 in 2018, indicating that the spatial distribution tended to be discrete.

Spatiotemporal Variations in Grassland Vulnerability
The QTP is a hotspot for vulnerability studies. Jiang et al. [51], Xia et al. [34], and Jiang et al. [52] assessed the vulnerability of the QTP or the Tibetan Plateau. The conclusions obtained were similar to those presented in this paper, with a dramatic trend of increasing vulnerability from the east and south to the west and north. Although the region had high VSI with strong human disturbance, the south-eastern QTP, which is composed of alpine meadow, had the lowest grassland vulnerability and was the main L-L region. Warm and humid air is brought from the Pacific, promoting vegetation growth and productivity, and making the grassland ecosystem less vulnerable [34,36]. Furthermore, the north-western QTP, which consisted of alpine steppe and desert steppe, had the highest grassland vulnerability and was the main H-H region. Because of the high elevations and low temperatures, it became barren and sparsely vegetated, making it vulnerable to environmental disturbances such as freeze-thaw erosion and desertification.
According to the temporal patterns, the grassland vulnerability of the QTP continuously improved, and the average value of GVI decreased from 2000 to 2018, with similar results found by Jiang et al. [52] and Jiang et al. [51]. However, based on trend analyses, Xia et al. [34] discovered that there were more areas with significantly increased vulnerability than areas with significantly decreased vulnerability. More detailed analyses showed that, following the environmental protection measures being implemented, the increasing trend of ecological vulnerability was alleviated, and more regions experienced significantly decreased vulnerability, suggesting that the government's ecological restoration programs were effective [34,51]. During the study period, many natural programmes on the QTP were proposed: the Ecological Environment Protection and Comprehensive Controlling Plan of Qinghai Lake Basin in 2007, the Protection and Construction of Ecological Security Barrier Plan in Tibet  in 2009, the Regional Ecological Construction and Environmental Protection Plan of Qinghai-Tibetan Plateau (2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021)(2022)(2023)(2024)(2025)(2026)(2027)(2028)(2029)(2030) in 2011, and so on [34,51]. These programmes promoted the return of farmland to grassland, biodiversity conservation, and improvement of the ecological environment. The Prohibited Grazing Policy (PGP) was another ecological initiative that had an impact on grassland. The PGP was a major ecological service project started by the Chinese government in 2003 to recover highly degraded grasslands and environmentally fragile regions [20]. Many studies concluded that the PGP had a significant impact on environmental protection and resulted in significant ecological benefits [20,53]. Meanwhile, due to global warming, increased precipitation, temperature, and growing season length caused improvements in the alpine meadows and steppes in some regions of the QTP [33,54]. The start of the thermal growing season (STGS), net primary production (NPP), NDVI, and other indicators improved, indicating that the general health of grasslands was improving on the QTP [54][55][56].

Driving Factors of Grassland Vulnerability
The driving factors of grassland vulnerability can be traced back to natural factors, environmental disturbances, and socioeconomic impacts. Based on the weights obtained, we analysed the driving factors of grassland vulnerability (Tables 2-4). For VNF, among the selected indicators, NDVI had a strong influence, and the weight of NDVI increased from 2000 to 2018. In related studies, NDVI was a key indicator of vulnerability [34,51]. Green vegetation played an important role in water and soil conservation, biodiversity, and climate adjustment [34]. The increase in NDVI resulted in a more stable grassland ecosystem. Moreover, both temperature and precipitation were essential environmental elements that influenced vegetation dynamics. The weight of the annual average temperature decreased, while the weight of the annual average precipitation increased. Possibly due to climate change, temperature was not a limiting factor in grassland vulnerability, but precipitation remained a limiting factor.
Based on the background natural factors, land desertification was the most important environmental disturbance, followed by freeze-thaw erosion and landslides. In comparison, the significance of the grazing pressure was relatively low, with similar results found by Jiang et al. [51]. Grassland, especially alpine and desert steppes, was vulnerable to land desertification. Permafrost degradation has resulted in increased land desertification in previously permafrost-covered areas. The duration and annual mean area of frozen soil have decreased significantly as a result of global warming, making land desertification more severe [57]. Increased precipitation and runoff, combined with the steep mountain terrain, intensified the potential energy conditions, resulting in unstable rocks and dirt sliding down, causing collapses and landslides [44,49]. In addition, the intervention of human social and economic activities was an essential driving factor in grassland vulnerability change. The most essential socioeconomic impact indicator was the population, which reflected the frequency and intensity of human disturbance. Among the indicators, the weight of GDP increased significantly from 0.24 to 0.44 over the years examined, suggesting that the QTP grassland was experiencing rapid economic growth. The importance of actual livestock was relatively low, as a result of local governments restricting overgrazing to safeguard the QTP.
The GVI incorporated several contributing elements, including natural factors, environmental disturbances, and socioeconomic impacts, showing the spatiotemporal distributions of grassland vulnerability. Understanding the mechanisms of driving factors is critical to reducing vulnerability. However, calculating the index of grassland vulnerability to make judgements of vulnerability is insufficient. Long-term categorized monitoring and projections of driving factors are required to gain a thorough understanding of grassland vulnerability [28,51].

Suggestions and Uncertainty Analysis
The goal of the vulnerability assessment was to provide decision-makers with recommendations on how to improve environmental quality [22]. To select priority regions in the vast grassland of the QTP, indicators that convey information about the condition of the grassland ecosystems must be used [1]. In this paper, the grassland of QTP was divided into three types, namely, alpine meadow, alpine steppe, and desert steppe, which accounted for 60.22%, 33.82%, and 5.96% of the total area, respectively. Our results were similar to those of Li et al. [1], who conducted an alpine grassland vulnerability assessment of the QTP using the SPR model. Desert steppe was the most vulnerable (0.78), followed by alpine steppe (0.76) and alpine meadow (0.64). Desert steppe was the transitional region from grassland to desert, and it was the driest grassland. It was particularly sensitive to human activity and climate change, and land desertification was a strong possibility. Therefore, to slow the spread of deserts, national funds should be invested in combating desertification and expanding the construction of shelterbelts [22]. The alpine meadow and alpine steppe were suitable for grazing, while overgrazing had negative impacts on soil, even leading to grassland degradation. In severely degraded grasslands, rational management policies and strategies should be implemented [58,59]. To ensure the long-term sustainability of grassland, it was recommended that more ecological projects should be implemented and that the management of the project's later stages should be strengthened.
This study reflected the spatiotemporal variations in QTP grassland vulnerability during the last 20 years. By selecting relevant indicators based on the VNF-VED-VSI framework and using a specific type of ecosystem as the study area, a comprehensive picture of ecosystem vulnerability on a regional scale could be created. However, this study also had limitations. First, because the assessment methods for factors have not yet been harmonized, the accuracy of the results requires further testing. A simple and accurate evaluation framework, as well as a method for factors, must be developed. Second, due to the large area of the QTP region, this study applied a 1000 m resolution, which is prone to accuracy errors. As a result, in future research, a smaller scale, higher resolution, and longer time series to assess the vulnerability of key areas should be considered. Third, ecosystem processes are complex and difficult to fully identify, especially when the ecosystem is disrupted by human activities. In this study, the impact of human activities by VSI was reflected using the indicators of population, GDP, and actual livestock, but it was far from sufficient. In future research, more appropriate indicators and more accurate models should be used to reflect complex ecosystem processes.

Conclusions
The grassland ecosystem accounted for the most significant proportion of the QTP and played important roles in conserving the soil and water, maintaining the biodiversity, and controlling the wind and sand. In this study, a comprehensive grassland vulnerability evaluation framework was developed by integrating the grassland vulnerability index with AHP, PCA, and EVDI approaches. Based on this framework, our results showed that the vulnerability by natural factors and environmental disturbances decreased from 2000 to 2018, while the vulnerability by socioeconomic impacts increased on the QTP. The factors that had the greatest impact on VNF, VED, and VSI were NDVI, land desertification, and population, respectively. The spatial distribution of GVI had obvious heterogeneity, decreasing from northwest to southeast. The regions with serious and extreme vulnerability were mainly located in the north-western alpine steppe and desert steppe, while the regions with normal, slight, and moderate vulnerability were mainly concentrated in the south-eastern alpine meadow. The global Moran's I index of grassland vulnerability on the QTP was greater than 0, with a significant positive spatial correlation. The number of H-H and L-L units decreased from 2000 to 2018, indicating that the H-H and L-L Cluster regions tended to be discrete. In addition, the H-H and L-L aggregate distributions were highly coincident with the high and low value regions of GVI. These results suggested that the comprehensive GVI based on the AHP-PCA-EVDI approaches is useful for evaluating grassland vulnerability at the regional scale. Overall, our results indicated that understanding the variations in grassland vulnerability on the QTP is important for regional sustainable development in the context of intensified human activities and climate change.