Combined Effects of Warming and Grazing on Rangeland Vegetation on the Qinghai-Tibet Plateau

Climate warming has increased grassland productivity on the Qinghai-Tibet Plateau, while intensified grazing has brought increasing direct negative effects. To understand the effects of climate change and make sustainable management decisions, it is crucial to identify the combined effects. Here, we separate the grazing effects with a climate-driven probability model and elaborate scenario comparison, using the Normalized Difference Vegetation Index (NDVI) of the grassland on the Qinghai-Tibet Plateau. We show that grazing has positive effects on NDVI in the beginning and end of the growing season, and negative effects in the middle. Because of the positive effects, studies tend to underestimate and even ignore the grazing pressure under a warming climate. Moreover, the seasonality of grazing effects changes the NDVI-biomass relationship, influencing the assessment of climate change impacts. Therefore, the seasonality of grazing effects should be an important determinant in the response of grassland to warming in sustainability analysis.


INTRODUCTION
Grasslands provide food and habitats for humankind and animals (de Jong et al., 2013), but are facing increasing pressure under climate change, in addition to intensifying human activities and increasing food demands (Chen et al., 2013;Seddon et al., 2016;Zhu et al., 2016). A key concern, at the current time, is whether grassland systems are sustainable under these double threats. While grasslands are among the most important indicators of the impacts of climate change (Li et al., 2018a), only very little research has been done on the role of human activities on the grasslands in climate change models (Väisänen et al., 2014).
The Three-River Headwaters Region (TRHR) on the Qinghai-Tibet Plateau is a crucial rangeland ecosystem with great importance to water resources and ecological security of China ( Figure 1A) (Zhang et al., 2016;Han et al., 2018). Dubbed as China's Water Tower, the TRHR is the source of the Yellow River, the Yangtze River, and the Lancang River (the upper part of the Mekong River). The Plateau grassland, covering 91% of the land surface with a mean altitude of 4,500 m above mean sea level (Zheng et al., 2018), is mainly used as rangeland of nomadic Tibetans, providing grass for Tibetan yak and sheep, as well as wild animals, to graze ( Figures 1B,C).
Weather station data revealed that air temperature and precipitation had increased significantly from 1960 to 2018 in the TRHR, and the increasing trend had particularly accelerated since the 1990s ( Figures 1D,E). The warming trend in the TRHR (at a rate of 0.036°C/yr) had been much stronger than that at the global level (0.016°C/yr), with more than double the rate ( Figure 1D). Climate change, including warming and CO 2 fertilizing, currently contributes to the greening of vegetation ( Figure 1F) in the watertemperature-constrained TRHR (Nemani et al., 2003). On the other hand, the grazing intensity (represented by meat production throughout this paper) almost doubled from 1995 to 2015 ( Figure 1G), exerting great pressure on the grassland ecosystem. The sustainability of the rangeland as well as the heritage of this unique Tibetan nomadic culture are facing huge challenges. However, the effects of grazing and climate change combine together to impact vegetation [e.g., yielding an overall rising Normalized Difference Vegetation Index (NDVI)], thus making it difficult to analyze and assess the influences and impacts, not to mention the complex nature of the individual factors involved as well as their interactions.
The combined effects of climate change and grazing on vegetation productivity, phenology and ecosystem properties have been explored by some studies (Klein et al., 2007;Zhang et al., 2015;Li et al., 2018c;Kohli et al., 2021). On the one hand, grazing alters the response of vegetation to climate change by modulating the dependency of vegetation growth on temperature (Wei et al., 2020). On the other hand, warming offsets the grazing effects on vegetation through altering the vegetation living state, i.e. increasing the plant height and the aboveground biomass (Zhang et al., 2015;Tang et al., 2019). The complexity associated with climate change and grazing and, hence, their combined effects on vegetation has generated contrasting results (Li et al., 2011;Wang et al., 2012;Xu et al., 2014;Ge et al., 2021). For example, grazing could exert both positive and negative influences on biomass due to the complex relationships between livestock grazing and biomass (Li et al., 2019a). How climate change and grazing interactively affect the biomass, biodiversity and ecosystem sustainability needs further investigation (Li et al., 2018c).
We analyzed the combined effects of climate change and grazing with elaborate scenario comparison (O'Neill and Nakicenovic, 2008;Huntzinger et al., 2013;McKenna et al., 2021), using a climate-driven probability model for time series, i.e., the Non-Homogeneous Hidden Markov Model (NHMM) (Holsclaw et al., 2017). First, we used the NHMM to simulate the vegetation dynamics month by month at county scale (a total of 414 months and 22 counties, Figure 1B, Supplementary Table  S1), with the standardized climate inputs and NDVI. Then, we devised two scenarios to separate the grazing effects. We find that, besides the negative effects on NDVI in the middle of the growing season (May to September), grazing shows surprising but obvious positive effects in the beginning and end of the growing season. This seasonal pattern has important implications for better understanding of the grazing pressure on grasslands and for modifying the existing NDVI-biomass relationships, towards a more reliable indication of climate change.

Data
Climate data. In general, climatic factors impose complex and varying limitations on vegetation growth (Nemani et al., 2003). Therefore, monthly precipitation and temperature data were obtained from the Climate Hazards Group Infrared Precipitation with Stations (CHIRPS) (Funk et al., 2015) and ERA-Interim (Dee et al., 2011), respectively, to feed the NHMM. The data from CHIRPS have a spatial resolution of 0.05°for the period from 1981 to the time being, and the monthly-mean surface temperature at 2-m from ERA-interim has a resolution of 0.125°for the period from 1979 to the present. To evaluate the amplitude of climate change on the Qinghai-Tibet Plateau study area, the monthly gauging station data from 1960 to 2018 provided by the China Meteorological Administration were used. Moreover, global temperature anomaly data were obtained from the National Oceanic and Atmospheric Administration (NOAA) to make comparison between regional and global warming, and to demonstrate the relatively extraordinary effects of warming on rangeland vegetation on the Qinghai-Tibet Plateau. Since the non-stationary climate system might also be related to vegetation growth, the climate oscillation indexes from NOAA were used as candidate predictors in our model. Vegetation data. NDVI, a proxy of vegetation greenness and production, was acquired from the Global Inventory Monitoring and Modelling Studies (GIMMS) NDVI3g (Tucker et al., 2005). The biweekly datasets span from July 1981 to December 2015 with a spatial resolution of 0.083°. Monthly NDVI was composited from the biweekly data using the Maximum Value Composite method (Holben, 1986). The moderate-resolution imaging spectroradiometer (MODIS) land cover data (MCD12Q1) (Friedl et al., 2010), with 500 m resolution, were utilized to composite a binary pasture mask. The grassland pixels were classified as the pasture area. Apart from the non-grassland pixel, the surrounding eight pixels were masked out to make the classifying result more robust.
Anthropogenic data. Yak and sheep meat production data, for the period 1995-2015, were obtained from the records of the Qinghai Statistical Almanacs, considering the lack of information and reliability on the livestock population data. Here it is noted that during our study analysis, we have found a clear correlation between the meat production and the livestock population data, thus meat production could be used as a convincing surrogate for the grazing effect. The livestock density (LD) for each county was defined as meat production divided by pasture area. Annual time series were converted to monthly series based on linear interpolation. The Gridded Population of the World, Version 4 (Center for International Earth Science Information Network -CIESIN -Columbia University, 2018) was used to demonstrate population density of the study area.
The linear regression method was applied for the trend analysis, with Student's two-tailed t-test for significance test, and the Pettitt change point test was used for change point detection (Pettitt, 1979). Since livestock data were at the county scale, the gridded climate data and NDVI were aggregated to county scale after eliminating the non-pasture area. The assumption for the attribution of grazing effects was that the climate zone was homogeneous across the TRHR and, therefore, it was reasonable to calibrate the NHMM for one county and apply it to another. To meet this assumption, the climate data and NDVI were further processed with standardized normalization. The standardized NDVI was denoted as NDVI p . The period from July 1981 to December 2015 (a total of 414 months) was selected for the model construction, as GIMMS NDVI3g was available only for this time period. Viovy and Saint (1994) adopted the Hidden Markov Model (HMM) for modelling the vegetation dynamics, but it has the issues of poor representation in the temporal variability and spatial variability (Holsclaw et al., 2016). In our study, the NHMM (Hughes and Guttorp, 1994;Hughes et al., 1999) was applied for modelling the monthly vegetation time series, i.e. NDVI * . We chose the NHMM for its simplicity; further, it could not only simulate the complicated temporal variability in the response using hidden state variables, but also consider climate variables as the external factors. The model involved an underlying hidden, or not observed, stochastic process, during which the hidden state evolved according to a first-order Markov chain (Supplementary Figure S1). In vegetation modelling, the hidden state was usually interpreted as a weather state. The observed process, e.g. vegetation index, was conditional on the hidden states through the state-dependent emission distribution. The conditional likelihood for the model can be written as:

Model
where P(z t |z t−1 , X t , ζ) is the Markov transition probability modelled by a multinomial logistic link function at time t, f(Vg t |W t , z t , θ) is the emission distribution modelled via a Frontiers in Environmental Science | www.frontiersin.org December 2021 | Volume 9 | Article 797971 mixture of Normal distribution, z is the hidden state, X t are exogenous variables that regulate the transition of hidden Markov state, W t are exogenous variables that control the possibility of emission distribution, and ζ and θ are for the transition parameters and emission distribution parameters, respectively. In the NHMM, we introduced two types of exogenous variables (Holsclaw et al., 2016;Holsclaw et al., 2017). Time series X were the first exogenous covariate series, reflecting largescale spatial phenomena and long-range temporal effects. In this study, the candidate variables for X included four harmonic terms (sin( π 6 t), sin( π 3 t), cos( π 6 t), and cos( π 3 t)), representing seasonal variability, and five climate oscillation indexes (Pacific Decadal Oscillation, Arctic Oscillation, North Atlantic Oscillation, Southern Oscillation Index, and Pacific North American Oscillation) related to the study area, representing inter-annual variability. The other exogenous covariate time series, W, were introduced to allow the underlying distributional characteristics of vegetation to change for each hidden state. Here, monthly precipitation and temperature, 2-4 months accumulated precipitation and temperature, and 2-3 months lagged precipitation and temperature were chosen as the candidate variables for W. In total, we proposed nine candidate exogenous variables for X and 12 for W. The least absolute shrinkage and selection operator (LASSO) was performed to reduce the dimension of the candidate variables and to avoid overfitting. To ensure that the NHMM functioned well and was parsimonious (Supplementary Figure S2), five exogenous variables were retained, i.e. sin( π 6 t), sin( π 3 t), monthly precipitation, 2-months accumulated precipitation, and 2-months accumulated temperature.
The Markov chain Monte Carlo (MCMC) algorithm (Holsclaw et al., 2017) was applied to estimate the posterior full conditional distributions for z, ζ, and θ. Two metrics for quantifying the performance of the NHMM, i.e. Bayesian Information Criteria (BIC) for model calibration performance and Predictive Log Score (PLS) for model validation performance (Holsclaw et al., 2017), were utilized to determine the number of hidden states (K) and the number of mixed Normal distributions (nmix) in emission distribution. The model with the minimum BIC and maximum PLS was the preferable one. The BIC was estimated from the model calibrated with the first 31 years of data (372 monthly data), and the PLS was calculated from the model validated with the remaining 42-months data. That is to say, the first 372-months data were used for the model calibration, and the remaining 42-months data were used for the model validation. For parsimony and interpretability, we chose the model with seven states and two mixed Normal distributions, as BIC score and PLS score were found to be the best for this model (Supplementary Figure S3).

Scenario Setting
For the purpose of attribution of grazing effects, we devised two scenarios in opposite directions. Scenario 1 (S1) was to calibrate and verify the model in five counties that had the least grazing density, to get five sets of model parameters (i.e., five similar calibrated models), and then apply them to the remaining 17 counties (Supplementary Table S1; Supplementary Figure S4).
For each of the 17 application sites, we got five sets of simulation results. Scenario 2 (S2) was the opposite of S1, i.e. to calibrate the models in five counties that had the highest grazing density. For either of these scenarios, the grazing effects were estimated by comparison between the five-county and 17-county groups, and the results were cross-verified by the opposite directions for confirmation. The NHMMs always explicitly used local climate forcing as inputs for each county, and then yielded NDVI results for each county. This way, the results for the application counties would reflect the grazing impacts, as at the counties where the model was calibrated, implicitly through parameters. Therefore, we hypothesized that the difference in the results and model performance of the NHMMs would reflect the livestock density that differs between model calibration and application counties. In S1, the simulated standardized NDVI values (NDVI * ref ) in model application counties have the least impacts of grazing (because S1 models were calibrated for the counties that had the least grazing density), and should have been greater (if grazing decreased NDVI) than the observed standardized values (NDVI * obs ) for each county, and vice versa in S2. If we denote NDVI * diff as the difference between NDVI * ref and NDVI * obs , then the relationship between NDVI * ref and NDVI * obs for S1 can be written as: where ε LD is the error induced by livestock impacts, ε r is the model random error, and NDVI * diff ε LD + ε r . For S2, the relationship between NDVI * ref and NDVI * obs was NDVI * obs NDVI * ref − ε LD + ε r . Here, NDVI * diff ε LD − ε r . For both scenarios, we assumed that |ε r | ≪ |ε LD |, thus NDVI * diff correlated with the difference in livestock density (LD − LD 0 ) between the model application counties and the model calibration ones.

Sensitivity Analysis
Consistent with the definition of ecological resilience (Holling, 1973), vegetation sensitivity to grazing was defined as the ratio of change in vegetation to change in livestock density: where S n is sensitivity, NDVI * obs corresponds to the livestock density at model application county, LD, while NDVI * ref corresponds to the model calibration county, LD 0 .
It was assumed that vegetation sensitivity to grazing was a function of climate forcing. Therefore, vegetation sensitivity to grazing was evaluated at different temperature-precipitation ranges.

Attribution of Grazing Effects
Consistent with our hypothesis, the performance metric of the model used in this study, i.e. root-mean-square deviation (RMSE), correlated with livestock density. For S1, the performance of the five model applications became worse Frontiers in Environmental Science | www.frontiersin.org December 2021 | Volume 9 | Article 797971 when they were migrated to counties with increasing livestock density ( Figure 2A). Meanwhile, for the models calibrated in counties that had the highest grazing density (S2), model performance became worse in counties that had less grazing density ( Figure 2B). The results suggested that the livestock density was a strong controlling factor of NDVI. As the results from the two scenarios yielded reliable cross-verification, we have general confidence that the model could detect the signal of the impacts of grazing. We attributed the difference in standardized NDVI (NDVI * diff ) between the results from the simulations and observations to grazing effects; a negative value means a grazing-induced decrease in NDVI. As shown in Figures 2C,D, NDVI * diff was plotted as the box-plots, and each box-plot stood for one model application in 1 month representing 17 counties and 34 years. To our surprise, the direction of NDVI * diff changed with season for the five sets of applications as well as for the two scenarios. This, however, was not due to a systematic model error or phenology difference at different elevations ( Supplementary Figures S5 and  S6). For July, August, and September (the middle of the growing season), the grazing effects were generally negative across the counties from 1982 to 2015. However, grazing mostly showed positive effects on NDVI for May, June, and October, i.e., at the beginning and end of the growing season. In the non-growing season, grazing did not exhibit any apparent positive or negative effects, partly because of the small and unreliable NDVI values (Tucker et al., 2005). The median value of NDVI * diff , which was integrated from the two scenarios, clearly exhibited the seasonal pattern ( Figure 2E). Therefore, we can construe that grazing activities can have seasonal, including both negative and positive, influence on rangelands on the Qinghai-Tibet Plateau.
Since NDVI * diff was induced by the difference in livestock density between the model application counties and the model calibration counties, we assumed that the magnitude of NDVI * diff correlated with the livestock density difference. Therefore, we examined the linear relationship between NDVI * diff and (LD − LD 0 ) for each month using data points from all scenarios, years, and counties. The results ( Figure 2F, Supplementary Figure S7) showed that the slope (i.e., sensitivity) had a similar seasonal pattern as that of NDVI * diff , and the relationships were statistically significant for every month. Therefore, it may be inferred that, even though the direction of the effects of grazing differed seasonally, the magnitude of NDVI was sensitive to livestock density for each month.

Seasonal Grazing Effects on Vegetation
With these results, we quantified the seasonal grazing effects on rangeland NDVI for the TRHR, and only focused on the growing season. After recovering NDVI * ref to the same scale with the observation value, we found that, due to negative grazing effects, the peak value of NDVI decreased by 0.015 from an estimated 0.445 to the observed 0.430 in August ( Figure 3A). The NDVI decrease in July and September was less, with an amplitude of 0.005-0.010, respectively. In the early (May and June) and late (October) growing season, on the contrary, grazing seemed to benefit the grassland, as NDVI increased by 0.015 from an estimated 0.234 to the observed 0.249 in May. This observation is in agreement, to some extent, with the positive grazing effects reported by some past studies using controlled experiment (Klein et al., 2007), and also for other regions, such as the Inner Mongolian grassland (Ren et al., 2016) and the Norwegian tundra (Mysterud et al., 2011). For the early growing season, grazing could contribute to the removal of snow (Mårell et al., 2006) and dried grass, and stimulate young and protein-rich regrowth (Mysterud et al., 2011;Ren et al., 2016). For the late growing season, grazing could also melt snow, remove aged leaves, and delay the deterioration of plant biomass by keeping it in young phenological stages (Albon and Langvatn, 1992;Hebblewhite et al., 2008;Mysterud et al., 2011;Li et al., 2018b). Therefore, the effects of grazing, at least on NDVI, can be overall positive.
It is appropriate to note, at this point, that NDVI is an index of the earth surface's status of greenness (Nemani et al., 2003;Tucker et al., 2005). Grazing effects on increasing NDVI values in the early and late growing season are partly benefited from the exposure of chlorophyll-rich grass parts. This, however, does not guarantee the improvement of grassland productivity and ecology. An effective quantification of grassland productivity is biomass production. Here, we proposed a conceptual diagram to illustrate the seasonal processes of biomass production rate, grazing consumption rate, and residual biomass ( Figure 3B). During the growing season, the mean temperature was above zero and the amount of precipitation constituted most of a year. The biomass production rate was estimated from NDVI using the empirical relationship (Zhao and Running, 2010), and the grazing consumption rate was assumed to be constant, for the simplicity. The residual biomass was calculated from the accumulated difference between the biomass production rate and the grazing consumption rate. Within an annual cycle, the production rate of biomass only exceeded the grazing consumption rate during the middle of the growing season. Therefore, the residual biomass reached its maximum after maximum NDVI appeared in July or August without any coverage.

Implications of Seasonal Grazing Effects
From biomass point of view, the sustainability of rangeland is determined by considering whether the maximum residual biomass, reached in the late growing season, is sufficient for grazing until biomass production surplus occurs in the next growing season. Both climate and grazing factors influence the accumulation for the maximum residual biomass. The climate is related to the biomass production, and the grazing is linked with the biomass consumption. To represent the maximum residual biomass by NDVI, as well as the grazing effects on biomass, the best indicator is the NDVI in August. This is because August is the closest month before the maximum biomass appears in September and has the maximum in both NDVI value and grazing effects on NDVI.
In view of this, we plotted the monthly mean NDVI anomaly (NDVI mon ) of TRHR for the period 1982-2015 in Figure 4A, where August (NDVI aug ) and growing season (NDVI gs ) were marked. We also used the monthly values to calculate and compare their trends. For NDVI gs and NDVI mon , the growing trends were statistically significant; on the other hand, the trend of NDVI aug was milder and insignificant. We also calculated the inter-annual trends by using the annual means of NDVI gs and NDVI mon , and found similar results  Figure S8). The NDVI trend results not only provided further evidence for the negative grazing effects in August, but also cautioned about possible underestimation of the grazing effects if using whole-year or growing-season NDVI for analysis. Underestimation of the grazing effects by using NDVI gs trend, instead of NDVI aug trend, could also be observed in terms of spatial distribution. As shown in Figures 4B,C, the blue and red pixels indicate the greening and browning trends with significance (p < 0.05), respectively, and the white ones represent the trend without significance. For NDVI gs , the percentage of pixels with significant greening and browning trends were 37.44% and 5.84%, and while for NDVI aug , the numbers were 23.83% and 9.49%, respectively. The greening trend estimated from NDVI aug was not as prevailing as that from NDVI gs . Besides, a large number of pixels with browning trend were found in the southern parts of the TRHR (Circle I and Circle II in Figures 4B,C) by using NDVI aug , and the locations of these pixels matched well with the regions with high population density ( Figure 1C). The results also cautioned about the The positive grazing effects in the beginning and end of the growing season would offset the negative grazing effects on NDVI. However, they had little to do with biomass production and accumulation. Therefore, when the overall whole-year or growing-season trend was used to investigate the rangeland sustainability, it might have given an illusion that the climate change benefit overwhelmed the grazing impacts. This might not be true, however; the grassland condition index that represented the yearly biomass surplus best, i.e., August NDVI, did not have a significant increasing trend.
The difficulty in assessing the effects of grazing also emphasizes the need to better understand the climate change effects on vegetation. Grazing activities regulate the responses of the grassland to climate change, even if biomass is used instead of NDVI. Because remotely-sensed biomass is also derived from NDVI, the underlying relationship between NDVI (for greenness) and biomass (for production) is widely evaluated for different vegetation types and seasons (Hansen and Schjoerring, 2003;Wessels et al., 2006); however, grazing and other human activities can change the NDVI-biomass relationship by reducing biomass and shifting the phenological process (Li et al., 2019b), with inconsistent change in greenness. Unfortunately, this has seldom been considered in studies evaluating the large-scale effects of climate change on vegetation, which requires some modifications to consider the grazing effects. Particularly, since the seasonality in grazing effects changes the NDVIbiomass relationships, modifications to the existing relationships should be different for different seasons, and even for different months.
Finally, we examined the historical trend of grazing effects and also whether the rangeland ecosystem would be sustainable in the future. For Figure 5A, each boxplot represents August NDVI * diff of two scenarios and five sets of applications for all counties in August. The trend of August NDVI * diff was calculated by linear fitting of the median values in the boxplots. The whole trend for August NDVI * diff was insignificant from 1981 to 2015, but a change point was found in 2009. The sudden rise in NDVI in 2009 and 2010 ( Figure 1F), contributing to the change point to positive August NDVI * diff , was mainly caused by abundant precipitation (Chen et al., 2020). However, for the divided subperiods (i.e. before and after 2009), though the vegetation growth was under different climatic conditions (Supplementary Figure S9), the grazing impacts increased and drew August NDVI * diff back to negative. The results revealed that even though the degree of grazing pressure exerted on grassland was influenced by climate, the increasing trend of grazing impacts still existed. We also evaluated the response of grazing pressure on climate variables, using the percentage of negative S n in different temperature-precipitation ranges. Figure 5B was calculated from two scenarios and five sets of applications for all the counties and all the years, plotted by a 0.5°C temperature interval and 10 mm precipitation interval. We observed a remarkable pattern; the percentage of negative effects was the highest when both temperature and precipitation were lower than their historical means (dashed lines in Figure 5B). With increase in temperature and precipitation, the percentage of negative effects reduced. However, the overall percentage of negative effects still exceeded 50%, which revealed that grazing still exerted pressure on grassland under climatic conditions within this range.

CONCLUSION
This research revealed the seasonal grazing effects on rangeland NDVI, biomass, and its sustainability in the TRHR on the Qinghai-Tibet Plateau. It also revealed the possible underestimation of grazing effects when the growing-season or whole-year NDVI values were used in evaluation, as what has normally been done in most studies. The seasonality in grazing effects on NDVI also brought uncertainties in the application of the existing NDVI-biomass relationships to derive biomass data products, which have commonly been used for climate change impact assessment. Although climate change has been found to provide more favorable conditions for grassland on the Qinghai-Tibet Plateau, the grazing effects have been generally negative and their magnitude still growing. Therefore, the sustainability of the rangeland Frontiers in Environmental Science | www.frontiersin.org December 2021 | Volume 9 | Article 797971 8 ecosystem might degrade, provided grazing intensity continues to rise and warming temperature causes droughts in some regions by increasing the evapotranspiration demands. While there is a broad consensus among researchers that climate in the future will continue to change, evidence shows that the relationship between temperature and vegetation activity is weakening for the boreal region (Piao et al., 2014). Therefore, more attention should be paid to the combined effects of climate change and human activities on those fragile vegetation ecosystems. Further experimental and field observations will help to improve our statistical model and disclose the fundamental physics of the problem. This would constitute our future work.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.