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).

FIGURE 1

FIGURE 1. The background and climate change of the Three-River Headwaters Region (TRHR). (A) Location of the TRHR; (B) Land cover of the TRHR; (C) Population density in 2015; (D) Station-based annual mean air temperature from 1961 to 2018; (E) Station-based annual mean precipitation from 1961 to 2018; (F) Growing season (May to September) NDVI from 1981 to 2015; (G) Grazing meat production from 1995 to 2015.

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₂ fertilizing, currently contributes to the greening of vegetation (Figure 1F) in the water-temperature-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.

Materials and Methods

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 ${\text{NDVI}}^{\ast }$. 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.

Model

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. ${\text{NDVI}}^{\text{*}}$. 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:

$\mathit{P}\left(\mathit{Vg}\right|\mathit{W},\mathit{X},\mathit{z}, \mathit{\zeta }, \mathit{\theta })={\displaystyle \prod }_{\mathit{t}=1}^{\mathit{T}}\mathit{P}\left({\mathit{z}}_{\mathit{t}}\right|{\mathit{z}}_{\mathit{t}-1},{\mathit{X}}_{\mathit{t}}, \mathit{\zeta })\mathit{f}\left(\mathit{V}{\mathit{g}}_{\mathit{t}}\text{|}{\mathit{W}}_{\mathit{t}},{\mathit{z}}_{\mathit{t}}, \mathit{\theta }\right)(1)$

where $P\left(z_t\right|z_{t-1},X_t,\zeta )$ is the Markov transition probability modelled by a multinomial logistic link function at time t, $f\left(V{g}_t\text{|}{W}_t,z_t,\theta \right)$ is the emission distribution modelled via a 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 $\zeta$ and $\theta$ 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 large-scale spatial phenomena and long-range temporal effects. In this study, the candidate variables for $X$ included four harmonic terms ($\text{sin}\left(\frac{\pi }{6}\text{t}\right)$, $\text{sin}\left(\frac{\pi }{3}\text{t}\right)$, $\text{cos}\left(\frac{\pi }{6}\text{t}\right)$, and $\text{cos}\left(\frac{\pi }{3}\text{t}\right)$), 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. $\text{sin}\left(\frac{\pi }{6}\text{t}\right)$, $\text{sin}\left(\frac{\pi }{3}\text{t}\right)$, 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$, $\zeta$, and $\theta$. 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 ($\text{K}$) and the number of mixed Normal distributions ($\text{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 $\text{NDVI}$ values (${\text{NDVI}}_{\text{ref}}^{\text{*}}$) 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 (${\text{NDVI}}_{\text{obs}}^{\text{*}}$) for each county, and vice versa in S2.

If we denote ${\text{NDVI}}_{\text{diff}}^{\text{*}}$ as the difference between ${\text{NDVI}}_{\text{ref}}^{\text{*}}$ and ${\text{NDVI}}_{\text{obs}}^{\text{*}}$, then the relationship between ${\text{NDVI}}_{\text{ref}}^{\text{*}}$ and ${\text{NDVI}}_{\text{obs}}^{\text{*}}$ for S1 can be written as:

$\mathbf{NDV}{\mathbf{I}}_{\mathbf{obs}}^{\ast }=\mathbf{NDV}{\mathbf{I}}_{\mathbf{ref}}^{\ast }\mathrm{+}{\epsilon }_{\mathbf{LD}}+{\epsilon }_{\mathit{r}}(2)$

where ${\epsilon }_{\text{LD}}$ is the error induced by livestock impacts, ${\epsilon }_r$ is the model random error, and ${\text{NDVI}}_{\text{diff}}^{\text{*}}={\epsilon }_{\text{LD}}+{\epsilon }_r$.

For S2, the relationship between ${\text{NDVI}}_{\text{ref}}^{\text{*}}$ and ${\text{NDVI}}_{\text{obs}}^{\text{*}}$ was ${\text{NDVI}}_{\text{obs}}^{\text{*}}={\text{NDVI}}_{\text{ref}}^{\text{*}}-{\epsilon }_{\text{LD}}+{\epsilon }_r$. Here, ${\text{NDVI}}_{\text{diff}}^{\text{*}}={\epsilon }_{\text{LD}}-{\epsilon }_r$.

For both scenarios, we assumed that $\left|{\epsilon }_r\right|\ll \left|{\epsilon }_{\text{LD}}\right|$, thus ${\text{NDVI}}_{\text{diff}}^{\text{*}}$ correlated with the difference in livestock density $(\text{LD}-{\text{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:

${\mathit{S}}_{\mathit{n}}=\frac{\mathbf{NDV}{\mathbf{I}}_{\mathbf{obs}}^{\ast }-\mathbf{NDV}{\mathbf{I}}_{\mathbf{ref}}^{\ast }}{\mathbf{LD}-\mathbf{L}{\mathbf{D}}_0}(3)$

where $S_n$ is sensitivity, ${\text{NDVI}}_{\text{obs}}^{\text{*}}$ corresponds to the livestock density at model application county, $\text{LD}$, while ${\text{NDVI}}_{\text{ref}}^{\text{*}}$ corresponds to the model calibration county, ${\text{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.

Results and Discussion

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 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.

FIGURE 2

FIGURE 2. Attribution of grazing effects. (A) RMSE for S1; (B) RMSE for S2; (C) ${\text{NDVI}}_{\text{diff}}^{\text{*}}$ for S1; (D) ${\text{NDVI}}_{\text{diff}}^{\text{*}}$ for S2; (E) Median ${\text{NDVI}}_{\text{diff}}^{\text{*}}$; (F) Slope for the linear relationship between ${\text{NDVI}}_{\text{diff}}^{\text{*}}$ and $(\text{LD}-{\text{LD}}_0)$ (Related scatter points are presented in Supplementary Figure S7).

We attributed the difference in standardized NDVI (${\text{NDVI}}_{\text{diff}}^{\text{*}}$) 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, ${\text{NDVI}}_{\text{diff}}^{\text{*}}$ 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 ${\text{NDVI}}_{\text{diff}}^{\text{*}}$ 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 ${\text{NDVI}}_{\text{diff}}^{\text{*}}$, 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 ${\text{NDVI}}_{\text{diff}}^{\text{*}}$ was induced by the difference in livestock density between the model application counties and the model calibration counties, we assumed that the magnitude of ${\text{NDVI}}_{\text{diff}}^{\text{*}}$ correlated with the livestock density difference. Therefore, we examined the linear relationship between ${\text{NDVI}}_{\text{diff}}^{\text{*}}$ and ($\text{LD}-{\text{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 ${\text{NDVI}}_{\text{diff}}^{\text{*}}$, 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 ${\text{NDVI}}_{\text{ref}}^{\text{*}}$ 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.

FIGURE 3

FIGURE 3. Seasonal grazing effects on vegetation. (A) Seasonality of grazing effects on NDVI; (B) Seasonality of grazing effects on biomass. The reference NDVI (${\text{NDVI}}_{\text{ref}}$) was estimated by deducting grazing effects using the ${\text{NDVI}}_{\text{diff}}^{\text{*}}\sim (\text{LS}-{\text{LS}}_0)$ relationship from the observed NDVI (${\text{NDVI}}_{\text{obs}}$).

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 (${\text{NDVI}}_{\text{mon}}$) of TRHR for the period 1982–2015 in Figure 4A, where August (${\text{NDVI}}_{\text{aug}}$) and growing season (${\text{NDVI}}_{\text{gs}}$) were marked. We also used the monthly values to calculate and compare their trends. For ${\text{NDVI}}_{\text{gs}}$ and ${\text{NDVI}}_{\text{mon}}$, the growing trends were statistically significant; on the other hand, the trend of ${\text{NDVI}}_{\text{aug}}$ was milder and insignificant. We also calculated the inter-annual trends by using the annual means of ${\text{NDVI}}_{\text{gs}}$ and ${\text{NDVI}}_{\text{mon}}$, and found similar results (Supplementary 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 $\text{NDVI}$ for analysis.

FIGURE 4

FIGURE 4. Assessment of grazing pressure on vegetation. (A) Inter-annual trends of ${\text{NDVI}}_{\text{aug}}$ ${\text{NDVI}}_{\text{gs}}$, and ${\text{NDVI}}_{\text{mon}}$ for the TRHR during the period 1982–2015; (B) Spatial distribution of inter-annual trend for ${\text{NDVI}}_{\text{gs}}$ during the period 1982–2015; (C) Spatial distribution of inter-annual trend for ${\text{NDVI}}_{\text{aug}}$ during the period 1982–2015.

Underestimation of the grazing effects by using ${\text{NDVI}}_{\text{gs}}$ trend, instead of ${\text{NDVI}}_{\text{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 ${\text{NDVI}}_{\text{gs}}$, the percentage of pixels with significant greening and browning trends were 37.44% and 5.84%, and while for ${\text{NDVI}}_{\text{aug}}$, the numbers were 23.83% and 9.49%, respectively. The greening trend estimated from ${\text{NDVI}}_{\text{aug}}$ was not as prevailing as that from ${\text{NDVI}}_{\text{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 ${\text{NDVI}}_{\text{aug}}$, and the locations of these pixels matched well with the regions with high population density (Figure 1C). The results also cautioned about the locations where pressures exerted by livestock grazing and other human activities were higher.

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 NDVI-biomass 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 ${\text{NDVI}}_{\text{diff}}^{\text{*}}$ of two scenarios and five sets of applications for all counties in August. The trend of August ${\text{NDVI}}_{\text{diff}}^{\text{*}}$ was calculated by linear fitting of the median values in the boxplots. The whole trend for August ${\text{NDVI}}_{\text{diff}}^{\text{*}}$ 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 ${\text{NDVI}}_{\text{diff}}^{\text{*}}$, was mainly caused by abundant precipitation (Chen et al., 2020). However, for the divided sub-periods (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 ${\text{NDVI}}_{\text{diff}}^{\text{*}}$ 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.

FIGURE 5

FIGURE 5. Trend and sensitivity of grazing impacts on August NDVI. (A) Inter-annual trend of August ${\text{NDVI}}_{\text{diff}}^{\text{*}}$; (B) Percentage of negative $S_n$ distributed in August climate space (with 0.5$\mathrm{℃}$ temperature interval and 10 mm precipitation interval).

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 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.

Author Contributions

CC, TL, and GW conceived the initial idea, CC and TL performed the analysis, CC, BS, AS, JA, LZ, and GW contributed to the manuscript, TL and GW provided the funding.

Funding

This study was funded by the Major Project of Qinghai Provincial Department of Science and Technology (2021-SF-A6), the State Key Laboratory of Hydroscience and Engineering of Tsinghua University (2017-KY-04), and the National Key Research and Development Program of China (2016YFE0201900, 2017YFC0403600).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenvs.2021.797971/full#supplementary-material

References

Albon, S. D., and Langvatn, R. (1992). Plant Phenology and the Benefits of Migration in a Temperate Ungulate. Oikos 65, 502–513. doi:10.2307/3545568

CrossRef Full Text | Google Scholar

Center for International Earth Science Information Network - CIESIN - Columbia University (2018). Gridded Population of the World, Version 4 (GPWv4): Population Density, Revision 11. Palisades, NY: NASA Socioeconomic Data and Applications Center. doi:10.7927/H49C6VHW

CrossRef Full Text | Google Scholar

Chen, C., Li, T., Sivakumar, B., Li, J., and Wang, G. (2020). Attribution of Growing Season Vegetation Activity to Climate Change and Human Activities in the Three-River Headwaters Region, China. J. Hydroinform. 22, 186–204. doi:10.2166/hydro.2019.003

CrossRef Full Text | Google Scholar

Chen, H., Zhu, Q., Peng, C., Wu, N., Wang, Y., Fang, X., et al. (2013). The Impacts of Climate Change and Human Activities on Biogeochemical Cycles on the Qinghai-Tibetan Plateau. Glob. Change Biol. 19, 2940–2955. doi:10.1111/gcb.12277

PubMed Abstract | CrossRef Full Text | Google Scholar

de Jong, R., Schaepman, M. E., Furrer, R., De Bruin, S., and Verburg, P. H. (2013). Spatial Relationship between Climatologies and Changes in Global Vegetation Activity. Glob. Change Biol. 19, 1953–1964. doi:10.1111/gcb.12193

CrossRef Full Text | Google Scholar

Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., et al. (2011). The ERA-Interim Reanalysis: Configuration and Performance of the Data Assimilation System. Q.J.R. Meteorol. Soc. 137, 553–597. doi:10.1002/qj.828

CrossRef Full Text | Google Scholar

Friedl, M. A., Sulla-Menashe, B., Schneider, A., Ramankutty, N., Sibley, A., and Huang, X. (2010). Sulla-MenasheMODIS Collection 5 Global Land Cover: Algorithm Refinements and Characterization of New Datasets. Remote Sensing Environ. 114, 168–182. doi:10.1016/j.rse.2009.08.016

CrossRef Full Text | Google Scholar

Funk, C., Peterson, P., Landsfeld, M., Pedreros, D., Verdin, J., Shukla, S., et al. (2015). The Climate Hazards Infrared Precipitation with Stations-A New Environmental Record for Monitoring Extremes. Sci. Data 2, 150066. doi:10.1038/sdata.2015.66

PubMed Abstract | CrossRef Full Text | Google Scholar

Ge, W., Deng, L., Wang, F., and Han, J. (2021). Quantifying the Contributions of Human Activities and Climate Change to Vegetation Net Primary Productivity Dynamics in China from 2001 to 2016. Sci. Total Environ. 773, 145648. doi:10.1016/j.scitotenv.2021.145648

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, Z., Song, W., Deng, X., and Xu, X. (2018). Grassland Ecosystem Responses to Climate Change and Human Activities within the Three-River Headwaters Region of China. Sci. Rep. 8, 9079. doi:10.1038/s41598-018-27150-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansen, P. M., and Schjoerring, J. K. (2003). Reflectance Measurement of Canopy Biomass and Nitrogen Status in Wheat Crops Using Normalized Difference Vegetation Indices and Partial Least Squares Regression. Remote Sensing Environ. 86, 542–553. doi:10.1016/s0034-4257(03)00131-7

CrossRef Full Text | Google Scholar

Hebblewhite, M., Merrill, E., and McDermid, G. (2008). A Multi-Scale Test of the Forage Maturation Hypothesis in a Partially Migratory Ungulate Population. Ecol. Monogr. 78, 141–166. doi:10.1890/06-1708.1

CrossRef Full Text | Google Scholar

Holben, B. N. (1986). Characteristics of Maximum-Value Composite Images from Temporal AVHRR Data. Int. J. Remote Sensing 7, 1417–1434. doi:10.1080/01431168608948945

CrossRef Full Text | Google Scholar

Holling, C. S. (1973). Resilience and Stability of Ecological Systems. Annu. Rev. Ecol. Syst. 4, 1–23. doi:10.1146/annurev.es.04.110173.000245

CrossRef Full Text | Google Scholar

Holsclaw, T., Greene, A. M., Robertson, A. W., and Smyth, P. (2016). A Bayesian Hidden Markov Model of Daily Precipitation over South and East Asia. J. Hydrometeorol. 17, 3–25. doi:10.1175/Jhm-D-14-0142.1

CrossRef Full Text | Google Scholar

Holsclaw, T., Greene, A. M., Robertson, A. W., and Smyth, P. (2017). Bayesian Nonhomogeneous Markov Models via Pólya-Gamma Data Augmentation with Applications to Rainfall Modeling. Ann. Appl. Stat. 11, 393–426. doi:10.1214/16-aoas1009

CrossRef Full Text | Google Scholar

Hughes, J. P., and Guttorp, P. (1994). A Class of Stochastic Models for Relating Synoptic Atmospheric Patterns to Regional Hydrologic Phenomena. Water Resour. Res. 30, 1535–1546. doi:10.1029/93wr02983

CrossRef Full Text | Google Scholar

Hughes, J. P., Guttorp, P., and Charles, S. P. (1999). A Non‐homogeneous Hidden Markov Model for Precipitation Occurrence. J. R. Stat. Soc. Ser. C (Applied Statistics) 48, 15–30. doi:10.1111/1467-9876.00136

CrossRef Full Text | Google Scholar

Huntzinger, D. N., Schwalm, C., Michalak, A. M., Schaefer, K., King, A. W., Wei, Y., et al. (2013). The North American Carbon Program Multi-Scale Synthesis and Terrestrial Model Intercomparison Project - Part 1: Overview and Experimental Design. Geosci. Model. Dev. 6, 2121–2133. doi:10.5194/gmd-6-2121-2013

CrossRef Full Text | Google Scholar

Klein, J. A., Harte, J., and Zhao, X.-Q. (2007). Experimental Warming, Not Grazing, Decreases Rangeland Quality on the Tibetan Plateau. Ecol. Appl. 17, 541–557. doi:10.1890/05-0685

PubMed Abstract | CrossRef Full Text | Google Scholar

Kohli, M., Mijiddorj, T. N., Suryawanshi, K. R., Mishra, C., Boldgiv, B., and Sankaran, M. (2021). Grazing and Climate Change Have Site‐dependent Interactive Effects on Vegetation in Asian Montane Rangelands. J. Appl. Ecol. 58, 539–549. doi:10.1111/1365-2664.13781

CrossRef Full Text | Google Scholar

Li, C., de Jong, R., Schmid, B., Wulf, H., and Schaepman, M. E. (2019a). Spatial Variation of Human Influences on Grassland Biomass on the Qinghai-Tibetan Plateau. Sci. Total Environ. 665, 678–689. doi:10.1016/j.scitotenv.2019.01.321

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, G., Jiang, C., Cheng, T., and Bai, J. (2019b). Grazing Alters the Phenology of alpine Steppe by Changing the Surface Physical Environment on the Northeast Qinghai-Tibet Plateau, China. J. Environ. Manage. 248, 109257. doi:10.1016/j.jenvman.2019.07.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L., Zhang, Y., Liu, L., Wu, J., Li, S., Zhang, H., et al. (2018a). Current Challenges in Distinguishing Climatic and Anthropogenic Contributions to alpine Grassland Variation on the Tibetan Plateau. Ecol. Evol. 8, 5949–5963. doi:10.1002/ece3.4099

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, P., Peng, C., Wang, M., Luo, Y., Li, M., Zhang, K., et al. (2018b). Dynamics of Vegetation Autumn Phenology and its Response to Multiple Environmental Factors from 1982 to 2012 on Qinghai-Tibetan Plateau in China. Sci. Total Environ. 637-638, 855–864. doi:10.1016/j.scitotenv.2018.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, W., Li, X., Zhao, Y., Zheng, S., and Bai, Y. (2018c). Ecosystem Structure, Functioning and Stability under Climate Change and Grazing in Grasslands: Current Status and Future Prospects. Curr. Opin. Environ. Sustainability 33, 124–135. doi:10.1016/j.cosust.2018.05.008

CrossRef Full Text | Google Scholar

Mårell, A., Hofgaard, A., and Danell, K. (2006). Nutrient Dynamics of Reindeer Forage Species along Snowmelt Gradients at Different Ecological Scales. Basic Appl. Ecol. 7, 13–30. doi:10.1016/j.baae.2005.04.005

CrossRef Full Text | Google Scholar

McKenna, C. M., Maycock, A. C., Forster, P. M., Smith, C. J., and Tokarska, K. B. (2021). Stringent Mitigation Substantially Reduces Risk of Unprecedented Near-Term Warming Rates. Nat. Clim. Change 11, 126–131. doi:10.1038/s41558-020-00957-9

CrossRef Full Text | Google Scholar

Mysterud, A., Hessen, D. O., Mobæk, R., Martinsen, V., Mulder, J., and Austrheim, G. (2011). Plant Quality, Seasonality and Sheep Grazing in an alpine Ecosystem. Basic Appl. Ecol. 12, 195–206. doi:10.1016/j.baae.2011.03.002

CrossRef Full Text | Google Scholar

Na, L., Genxu, W., Yan, Y., Yongheng, G., and Guangsheng, L. (2011). Plant Production, and Carbon and Nitrogen Source Pools, Are Strongly Intensified by Experimental Warming in alpine Ecosystems in the Qinghai-Tibet Plateau. Soil Biol. Biochem. 43, 942–953. doi:10.1016/j.soilbio.2011.01.009

CrossRef Full Text | Google Scholar

Nemani, R. R., Keeling, C. D., Hashimoto, H., Jolly, W. M., Piper, S. C., Tucker, C. J., et al. (2003). Climate-driven Increases in Global Terrestrial Net Primary Production from 1982 to 1999. Science 300, 1560–1563. doi:10.1126/science.1082750

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Neill, B. C., and Nakicenovic, N. (2008). Learning from Global Emissions Scenarios. Environ. Res. Lett. 3, 045014. doi:10.1088/1748-9326/3/4/045014

CrossRef Full Text | Google Scholar

Pettitt, A. N. (1979). A Non-parametric Approach to the Change-point Problem. Appl. Stat. 28, 126–135. doi:10.2307/2346729

CrossRef Full Text | Google Scholar

Piao, S., Nan, H., Huntingford, C., Ciais, P., Friedlingstein, P., Sitch, S., et al. (2014). Evidence for a Weakening Relationship between Interannual Temperature Variability and Northern Vegetation Activity. Nat. Commun. 5, 5018. doi:10.1038/ncomms6018

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, H., Han, G., Lan, Z., Wan, H., Schönbach, P., Gierus, M., et al. (2016). Grazing Effects on Herbage Nutritive Values Depend on Precipitation and Growing Season in Inner Mongolian Grassland. Jpecol 9, 712–723. doi:10.1093/jpe/rtw011

CrossRef Full Text | Google Scholar

Seddon, A. W. R., Macias-Fauria, M., Long, P. R., Benz, D., and Willis, K. J. (2016). Sensitivity of Global Terrestrial Ecosystems to Climate Variability. J. Plant Ecol. 531, 229–232. doi:10.1038/nature16986

CrossRef Full Text | Google Scholar

Tang, L., Zhong, L., Xue, K., Wang, S., Xu, Z., Lin, Q., et al. (2019). Warming Counteracts Grazing Effects on the Functional Structure of the Soil Microbial Community in a Tibetan Grassland. Soil Biol. Biochem. 134, 113–121. doi:10.1016/j.soilbio.2019.02.018

CrossRef Full Text | Google Scholar

Tucker, C. J., Pinzon, J. E., Brown, M. E., Slayback, D. A., Pak, E. W., Mahoney, R., et al. (2005). An Extended AVHRR 8‐km NDVI Dataset Compatible with MODIS and SPOT Vegetation NDVI Data. Int. J. Remote Sensing 26, 4485–4498. doi:10.1080/01431160500168686

CrossRef Full Text | Google Scholar

Väisänen, M., Ylänne, H., Kaarlejärvi, E., Sjögersten, S., Olofsson, J., Crout, N., et al. (2014). Consequences of Warming on Tundra Carbon Balance Determined by Reindeer Grazing History. Nat. Clim. Change 4, 384–388. doi:10.1038/nclimate2147

CrossRef Full Text | Google Scholar

Viovy, N., and Saint, G. (1994). Hidden Markov Models Applied to Vegetation Dynamics Analysis Using Satellite Remote Sensing. IEEE Trans. Geosci. Remote Sensing 32, 906–917. doi:10.1109/36.298019

CrossRef Full Text | Google Scholar

Wang, S., Duan, J., Xu, G., Wang, Y., Zhang, Z., Rui, Y., et al. (2012). Effects of Warming and Grazing on Soil N Availability, Species Composition, and ANPP in an alpine Meadow. Ecology 93, 2365–2376. doi:10.1890/11-1408.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, D., Zhao, H., Zhang, J., Qi, Y., and Wang, X. (2020). Human Activities Alter Response of alpine Grasslands on Tibetan Plateau to Climate Change. J. Environ. Manage. 262, 110335. doi:10.1016/j.jenvman.2020.110335

CrossRef Full Text | Google Scholar

Wessels, K. J., Prince, S. D., Zambatis, N., Macfadyen, S., Frost, P. E., and Van Zyl, D. (2006). Relationship between Herbaceous Biomass and 1‐km2 Advanced Very High Resolution Radiometer (AVHRR) NDVI in Kruger National Park, South Africa. Int. J. Remote Sensing 27, 951–973. doi:10.1080/01431160500169098

CrossRef Full Text | Google Scholar

Xu, X., Luo, Y., Shi, Z., Zhou, X., and Li, D. (2014). Consistent Proportional Increments in Responses of Belowground Net Primary Productivity to Long-Term Warming and Clipping at Various Soil Depths in a Tallgrass Prairie. Oecologia 174, 1045–1054. doi:10.1007/s00442-013-2828-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Gao, Q., Dong, S., Liu, S., Wang, X., Su, X., et al. (2015). Effects of Grazing and Climate Warming on Plant Diversity, Productivity and Living State in the alpine Rangelands and Cultivated Grasslands of the Qinghai-Tibetan Plateau. Rangel. J. 37, 57–65. doi:10.1071/rj14080

CrossRef Full Text | Google Scholar

Zhang, Y., Zhang, C., Wang, Z., Chen, Y., Gang, C., An, R., et al. (2016). Vegetation Dynamics and its Driving Forces from Climate Change and Human Activities in the Three-River Source Region, China from 1982 to 2012. Sci. Total Environ. 563-564, 210–220. doi:10.1016/j.scitotenv.2016.03.223

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, M., and Running, S. W. (2010). Drought-Induced Reduction in Global Terrestrial Net Primary Production from 2000 through 2009. Science 329, 940–943. doi:10.1126/science.1192666

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, Y., Han, J., Huang, Y., Fassnacht, S. R., Xie, S., Lv, E., et al. (2018). Vegetation Response to Climate Conditions Based on NDVI Simulations Using Stepwise Cluster Analysis for the Three-River Headwaters Region of China. Ecol. Indicators 92, 18–29. doi:10.1016/j.ecolind.2017.06.040

CrossRef Full Text | Google Scholar

Zhu, Z., Piao, S., Myneni, R. B., Huang, M., Zeng, Z., Canadell, J. G., et al. (2016). Greening of the Earth and its Drivers. Nat. Clim. Change 6, 791–795. doi:10.1038/nclimate3004

CrossRef Full Text | Google Scholar