Abstract
Over the past century, vegetation change has been reported at global, national, and regional scales, accompanied by significant climate change and intensified human activities. Among the regions is the rangeland of the Three-River Headwaters Region (TRHR) in China. However, which factor dominates in causing vegetation change in this region is still under considerable debate, and how would the grasslands adapt to the changing environment is largely unknown. To address these issues, we attribute growing season vegetation activity to climate change and human activities, investigate the interactions among different driving variables, and explore the dynamic relationship between vegetation activity and the driving variables. We perform Mann–Kendall trend analysis, Pearson correlation analysis, and partial correlation analysis. The results indicate that the dominant factor for vegetation growth, during the period 1995–2014, was temperature for the southeastern and southern parts of the TRHR, precipitation for the western part, and solar radiation for the northeastern part. The regulation effects of temperature on precipitation and cloud cover contributed to vegetation growth, while grazing activity and population activity offset the positive contribution of climate change. The dynamic relationship between vegetation activity and the driving variables reflected the acclimatization and adaption processes of vegetation, which needs further investigation.
Notation
- NDVI
Normalized Difference Vegetation Index
- NDVIGS
Growing season NDVI
- GP
Growing season precipitation
- GT
Growing season temperature
- GC
Growing season cloud cover
- LP
Annual livestock population
- POP
Annual human population
- RNDVI-GP
Partial correlation coefficient between NDVIGS and GP
- RNDVI-GT
Partial correlation coefficient between NDVIGS and GT
- RNDVI-GC
Partial correlation coefficient between NDVIGS and GC
- RNDVI-LP
Partial correlation coefficient between NDVIGS and LP
- RNDVI-POP
Partial correlation coefficient between NDVIGS and POP
INTRODUCTION
In the terrestrial biosphere, vegetation plays a vital role in providing food and habitats for humankind and animals (de Jong et al. 2013). Vegetation is one of the most important components of the ecosystem, as it bridges water, energy, and carbon transfer between the land and the atmosphere (Chen et al. 2013; Mao et al. 2016; Gao et al. 2017). Among the indicators of vegetation activity, the Normalized Difference Vegetation Index (NDVI) is a remote-sensed proxy for vegetation cover and photosynthetic capacity, and has been commonly used in vegetation-change studies (Peng et al. 2013; Wu et al. 2015). Using this index, vegetation change, including greening trend and browning trend, has been reported in global, national, and regional scales, generally based on simple linear models (Liu et al. 2014; Mao et al. 2016) or seasonal models (de Jong et al. 2011; Forkel et al. 2013). The study of vegetation dynamics under climate change has the potential to provide a theoretical basis for ecosystem-based adaptation and thus has become a major topic of global change research (Chen et al. 2013; Mao et al. 2016; Gao et al. 2017).
Vegetation activity is driven by both climatic factors and anthropogenic factors. In general, climatic factors include, among others, precipitation, temperature, and solar radiation, which entangle and interact with each other, imposing complex and varying limitations on vegetation activity (Nemani et al. 2003). For instance, Karnieli et al. (2010) reported that the relationship between NDVI and land surface temperature is regulated by water and solar radiation. Apart from climatic factors, vegetation change may be induced by anthropogenic drivers, such as land use change, fertilization, irrigation, urbanization, and grazing, among others (Diaz et al. 2007; Neigh et al. 2008; Zhu et al. 2016). Recently, attempts have been made to disentangle the contributions of multiple climatic and anthropogenic drivers from each other. For instance, Jiang et al. (2017) applied a residual analysis trend method to distinguish between the effects of precipitation, temperature, and the excessive exploitation of water resource, oil, and natural gas extraction on vegetation dynamics in Central Asia. Liu et al. (2015) analyzed the correlations of global NDVI trends with climatic and human-induced factors using the Human Influence Index. For the rangeland ecosystem, grazing effect is one of the most important anthropogenic drivers and has been assessed by Diaz et al. (2007) and Walker et al. (2009). However, which factor dominates in causing vegetation change in rangeland ecosystems is still under considerable debate (Archer 2004). Some believe that overgrazing and other anthropogenic impacts are the primary causes of grassland change, while others believe that climate change is the main driver (Dong et al. 2010). Debate on this is particularly relevant and important for large semi-arid rangeland regions, such as the Qinghai-Tibet Plateau in China (Wang et al. 2015; Li et al. 2018). Li et al. (2018) reported that four views exist for alpine grassland changes: climatic driver, non-climatic driver, the combination of climatic and anthropogenic drivers, or the alternation of climatic and anthropogenic drivers. They concluded that spatiotemporal heterogeneity of human activity intensity should be taken into account in the models for attribution. We address this debate and discuss the spatiotemporal heterogeneity of anthropogenic drivers in the present study.
In recent decades, the Qinghai-Tibet Plateau has undergone significant climate warming, with the magnitude of warming reported for the region to be larger than that for the Northern Hemisphere and for the entire globe (Kuang & Jiao 2016). The increase in precipitation has not been as pronounced as the increase in temperature (Xie et al. 2010). Nevertheless, it has shown more complicated spatially heterogeneous pattern (Li et al. 2014), which has resulted in different impacts on vegetation growth. Generally, grassland activity is extremely sensitive to temperature and precipitation changes, and also to non-climatic factors, including grazing, fires, nitrogen deposition, and rising CO2 levels (Zelikova et al. 2014; Li et al. 2018). With aggravating global climate change and increasing human activities, quantifying the effects of different driving variables on vegetation dynamics has become an important issue for undertaking countermeasures and management policies (Zhang et al. 2016). The gap in knowledge lies in how to detect and assess the impacts of climate change and human activities on the vegetation dynamics (Gao et al. 2016). To our knowledge, only a few studies have assessed the impacts of climate change on vegetation activities and the influencing mechanisms (Reyer et al. 2013; Gao et al. 2017). The present study attempts to bridge the gap in the current studies.
Over the past century, numerous methods have been developed to detect and assess the influence of variables/components in dynamic systems. In terms of vegetation ecology, the application of the partial correlation analysis to study the effects of climatic factors on vegetation activity has been attempted by a few studies, with some promising outcomes. Piao et al. (2014) argued that the weakening relationship between inter-annual temperature variability and northern sphere vegetation activity is induced by the effect of drought trends and increasing occurrence of extreme hot days. Cong et al. (2017) found that the impact of warming on vegetation activity in the Tibetan Plateau is more positive during periods with more precipitation. Peng et al. (2013) and Yang et al. (2017) have tried to explain the asymmetric responses of vegetation to daytime and night-time warming. In general, these studies have mainly focused on the effects of only climatic variables on vegetation dynamics. Only a few studies have discussed the interactions among different variables on vegetation dynamics (Zhang et al. 2013b; Cong et al. 2017), especially the interactions between climatic variables and anthropogenic variables. This provides the motivation for the present study.
In this study, we attempt to properly attribute the growing season vegetation activity to both climatic factors (precipitation, temperature, and cloud cover) and human activities (grazing effect and population pressure) in the Three-River Headwaters Region (TRHR), the hinderland of the Qinghai-Tibet Plateau. For this purpose, we perform Pearson correlation analysis and partial correlation analysis. First, we conduct trend analysis for the growing season vegetation, growing season climatic variables, and anthropogenic variables. To do this, we use NDVI data from the GIMMS3 g NDVI dataset, precipitation data from the Climate Hazards Group Infrared Precipitation with Stations (CHIRPS), temperature data from ERA-Interim, cloud cover data from the Climate Research Unit (CRU, version TS 4.01), and livestock production and human population data from the Qinghai Statistical Almanacs. Then, we perform Pearson correlation analysis and partial correlation analysis to detect the dominant factor(s) governing the vegetation dynamics. Here, we endeavour to figure out the interactions among the climatic variables as well as those between climatic variables and anthropogenic variables. Finally, we investigate the dynamic relationship between vegetation activity and the driving factors, and describe the mechanism that might be responsible for this dynamic relationship.
MATERIALS AND METHODS
Study area
The TRHR (Figure 1) is the source region for three major rivers in China, namely the Yangtze River, the Yellow River, and the Lancang River. It covers an area of 35,000 km2 and features a fluctuating terrain, dense river networks, extensive snowy mountains, and criss-crossing glaciers (Zhang et al. 2016). The region has a plateau continental climate and belongs to climate zones varying from semi-arid to semi-humid to humid areas from the northwest to the southeast. Along with Tibet, the region is regarded as the Asian Water Tower (Immerzeel et al. 2010), with altitude ranging from 2,000 m to nearly 6,800 m. Mainly covered with alpine meadow and alpine steppe, the region plays a key role in the ecological security of China due to its sensitive and fragile ecological environment (Zhang et al. 2016).
The TRHR in China. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2019.003.
The TRHR in China. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/hydro.2019.003.
Over the past half a century, the TRHR has experienced significant climate change (Li et al. 2010; Cuo et al. 2013; Wang et al. 2018). Overall, climate change in the TRHR has supposedly altered the biogeochemical cycles, hydrological regime, and vegetation activity profoundly (Chen et al. 2013; Mao et al. 2016; Shi et al. 2016). Apart from climate change, human activities in the TRHR have intensified during the past few decades. On the one hand, urbanization, land use change, overgrazing, and population growth have resulted in rangeland degradation in the TRHR (Harris 2010; Pan et al. 2017). On the other hand, ecological protection and restoration programmes have significantly improved the grassland quality and reduced grazing pressure, as revealed by Pan et al. (2017) and Zhang et al. (2017). The complexity associated with climate change and human activities and, hence, on their effects on vegetation has led to a key question: what are the major causes of vegetation change in the TRHR (Wang et al. 2015)?
Data
Normalized Difference Vegetation Index
The NDVI data used in this study is the NDVI3 g version from the Global Inventory Monitoring and Modelling Studies (GIMMS) group, derived from the NOAA/AVHRR land dataset. The dataset spans from July 1981 to December 2015 at a spatial resolution of 0.083° and at a temporal resolution of 15 days (Tucker et al. 2005). Monthly NDVI is composited from the fortnightly data using the Maximum Value Composite (MVC) method to remove noises from clouds and aerosols (Holben 1986). Areas with sparse or no vegetation cover (NDVI < 0) are masked out.
Climate data
The climate data used in this study are precipitation, temperature, and cloud cover. They are obtained from CHIRPS, ERA-Interim, and CRU TS 4.01, respectively. Monthly CHIRPS product, developed by the Climate Hazards Group, has a spatial resolution of 0.05° and available from 1981 to real time. The ERA-Interim dataset is a global atmosphere reanalysis dataset available from 1979 to real time. We make use of monthly-mean surface temperature at 2 m with a resolution of 0.125°. The CRU TS 4.01 dataset is interpolated from over 4,000 meteorological stations with spatial autocorrelation functions and has a spatial resolution of 0.5° (Harris et al. 2014). Besides, we use the gridded precipitation and temperature datasets from the China Meteorological Administration (CMA) to test the robustness of our analyses. The CMA gridded datasets have a spatial resolution of 0.5° and are interpolated from data obtained at the ground-level gauging stations.
Livestock and population
Livestock data (beef and mutton production) and human population data from 1995 to 2014 are acquired from the records of the Qinghai Statistical Almanacs. The datasets are grouped by the county scale.
In this study, growing season NDVI (NDVIGS), growing season precipitation (GP), growing season temperature (GT), growing season cloud cover (GC), annual livestock production (LP), and annual human population (POP) are used in the analysis. The growing season is defined as the 5-month period from May to September, and the growing season value is taken as the average of the five monthly values. NDVIGS, GP, GT, and GC for the TRHR and each county are averaged from pixels within the TRHR and each county, respectively. It should be noted that growing season vegetation is assumed to be related to annual LP and POP, due to data availability. This assumption is reasonable, for the growing season grass biomass constitutes most of the total annual biomass supporting livestock and population in the study area.
We select different periods for different analysis. In the trend analysis, NDVIGS, GP, GT, and GC are selected from 1982 to 2014 based on the temporal coverage of NDVIGS, to reflect the temporal variations as long as possible, while LP and POP are selected from 1995 to 2014 due to data availability. In the Pearson correlation analysis and partial correlation analysis, all the variables are selected from 1995 to 2014, since the anthropogenic variables are available only over this period.
A panel presentation of the overall datasets is presented in Figure 2. Spatially, the distribution of NDVIGS, GP, and GC roughly follows an increasing pattern from the northwestern part to the southeastern part, and the distribution of GT, LP, and POP exhibits an increasing pattern from the western part to the eastern part. Temporally, all the six variables show an upward trend from 1982 to 2014.
Spatial and temporal distribution of vegetation and the driving variables. (1) Spatial distribution of (a) NDVIGS, (c) GP, (e) GT, (g) GC, (i) LP, and (k) POP. NDVIGS, GP, GT, and GC are averaged from 1982 to 2014, and LP and POP are averaged from 1995 to 2014. (2) Temporal variation of (b) NDVIGS, (d) GP, (f) GT, (h) GC, (j) LP, and (l) POP in the TRHR.
Spatial and temporal distribution of vegetation and the driving variables. (1) Spatial distribution of (a) NDVIGS, (c) GP, (e) GT, (g) GC, (i) LP, and (k) POP. NDVIGS, GP, GT, and GC are averaged from 1982 to 2014, and LP and POP are averaged from 1995 to 2014. (2) Temporal variation of (b) NDVIGS, (d) GP, (f) GT, (h) GC, (j) LP, and (l) POP in the TRHR.
Methods
Trend analysis
The non-parametric Mann–Kendall test, which is widely used in meteorology, hydrology and ecology studies (Liuzzo et al. 2015; Mao et al. 2016; Shi et al. 2017), is adopted for trend detection in vegetation and the driving variables. The assumption of the Mann–Kendall test is that the data are serially independent (Liuzzo et al. 2015). Therefore, the trend-free pre-whitening method is applied to the data to eliminate the influence of autocorrelation before the application of Mann–Kendall test (Yue & Wang 2002). Considering that the coefficients obtained for lag-2 and lag-3 autocorrelation are small, only the lag-1 autocorrelation signal is removed from the original signal.
Pearson correlation analysis



Partial correlation analysis
To investigate the net effect of one independent variable on the dependent variable, there are two ways to control the effects of other variables: experimental and statistical (Garrett & Woodworth 1971). The experimental approach may reduce drastically the size of samples by selecting the appropriate observations. The statistical approach can eliminate the effects of other variables statistically without losing any samples. In the present study, the statistical approach is adopted. In particular, the partial correlation analysis is performed, as such an analysis is widely used in vegetation studies (Peng et al. 2013; Shen et al. 2015a; Wu et al. 2015; Cong et al. 2017). The function of the partial correlation analysis is to find the correlation between two variables by holding the other variables constant.
In this study, we conduct the partial correlation analysis between NDVIGS and only any one of the six driving variables at a time, controlling the other four variables. We try to disentangle the coordination effects between different variables in performing the partial correlation analysis.
RESULTS
Spatiotemporal change in vegetation
First, we apply the trend analysis to investigate the regional vegetation change in the TRHR. Overall, NDVIGS in the TRHR experienced, during 1982–2014, a slightly greening trend (0.36%/decade, p = 0.16) (Figure 2(b)). This result is in agreement with that reported by some earlier studies (Liu et al. 2014; Pan et al. 2017).
Next, we analyze the trend in NDVIGS data at grid and county scales in the TRHR for the same period (from 1982 to 2014). The distribution of greening and browning trend of vegetation is shown in Figure 3. As seen, the grids with greening trend outnumbered those with browning trend, and about 42.2% of grids in total experienced a significant vegetation change from 1982 to 2014 (p < 0.05). The areas with a browning trend, accounting for 24.1% of the total number of grids, were mainly located in the southeastern and southern parts of the TRHR, but only 19.3% of the browning grids show significant negative trend (p < 0.05, Figure 3(a)). The areas with a greening trend, accounting for 75.9% of the total number of grids, were distributed in the western and northeastern parts of the TRHR. Here, 49.4% of the greening grids exhibited significant positive trend (p < 0.05).
Distribution of trend in vegetation and the driving variables in the TRHR. (a) NDVIGS trend at the county scale; (b) NDVIGS trend at the grid scale; (c) GP trend at the county scale; (d) GP trend at the grid scale; (e) GT trend at the county scale; (f) GT trend at the grid scale; (g) GC trend at the county scale; (h) GC trend at the grid scale; (i) LP trend at the county scale; (j) POP trend at the county scale. The dots indicate significance at p< 0.05.
Distribution of trend in vegetation and the driving variables in the TRHR. (a) NDVIGS trend at the county scale; (b) NDVIGS trend at the grid scale; (c) GP trend at the county scale; (d) GP trend at the grid scale; (e) GT trend at the county scale; (f) GT trend at the grid scale; (g) GC trend at the county scale; (h) GC trend at the grid scale; (i) LP trend at the county scale; (j) POP trend at the county scale. The dots indicate significance at p< 0.05.
The trend of NDVIGS at the county scale is consistent with that at the grid scale (Figure 3(b)), with browning trend located in the southeastern and southern parts and greening trend located in the northeastern and western parts. The browning trend ranged from −0.08%/decade in the Tongren County to −0.22%/decade in the Gadee County, while the greening trend ranged from 0.09%/decade in the Baima County to 0.81%/decade in the Jianzha County. Overall, the trend of NDVIGS in three out of 22 counties in the TRHR was significant, and all of them showed a greening trend.
Spatiotemporal change in climatic and anthropogenic variables
We perform trend analysis for climatic and anthropogenic variables, and investigate the spatial and temporal characteristics of the trend in the TRHR. During the past three decades, climatic and anthropogenic variables in the TRHR have experienced a significant change. From 1982 to 2014, GP increased with a change rate of 0.78 mm/yr (p < 0.001), GT increased with a rate of 0.043 deg C/yr (p < 0.001), and GC increased significantly (p < 0.001) with a change rate of 0.11%/yr (Figure 2(d), 2(f), and 2(h)). Moreover, a significant increasing trend (p < 0.001) was found in LP and POP from 1995 to 2014 (Figure 2(j) and 2(l)).
With reference to trend in the grid scale and the county scale, both GP and GT in most grids and counties experienced an upward trend from 1982 to 2014 (Figure 3(c)–3(f)). The areas with a significant increase in GP were located in the northeastern and western parts of the TRHR, and the maximum rate was as large as 1.45 mm/yr (Figure 3(c) and 3(d)). It was found that 89.4% grids of the whole region, or 19 out of 22 counties, showed a significant change in GP. The rate of change in GT was found to be uniformly positive in all grids and counties (p < 0.05), ranging from 0.019 to 0.067 °C/yr (Figure 3(e) and 3(f)). The trend in the northwestern and southeastern parts is among the strongest. The spatial distribution of change in GC showed remarkable patterns. The significant increasing trend was mainly in the western part of the TRHR, while the decreasing trend was distributed in the southern part (Figure 3(g) and 3(h)).
From 1995 to 2014, LP and POP went up rapidly in most counties (Figure 3(i) and 3(j)). The rate of change in LP ranged from −0.01 to 0.08 ton/km2/yr, with the maximum rate observed in the northeastern part of the TRHR. The rate of change in LP was relatively smaller in the western part. The distribution of change in POP was similar to that of LP.
Pearson correlation between vegetation activity and the driving variables
To examine the relationship between vegetation activity and the driving variables, we perform the Pearson correlation analysis. Since the anthropogenic variables are available only at the county scale, we first check the relationship between vegetation activity and climatic variables based on the grid scale. For the period 1995–2014, we aggregate the original vegetation (NDVIGS) and climatic data (i.e. GP, GT, and GC) of different spatial resolutions into a single constant resolution (0.5°) and then conduct the Pearson correlation analysis between them, i.e. NDVIGS and GP, NDVIGS and GT, NDVIGS and GC. We use the Pearson correlation coefficient as an indicator of the relationship.
Figure 4 presents the relationship, indicated by the Pearson correlation coefficient, of NDVIGS with the three climatic variables GP, GT, and GC. During the period 1995–2014, the correlation between NDVIGS and GP was positive in over 75% of the total number of grids, but the positive relationship was significant in only about 25% of the grids. The negative relationship was mainly observed in the southeastern part of the TRHR. The relationship between NDVIGS and GT was also mainly positive, but negative relationship was found in the northeastern part of the TRHR. The relationship between NDVIGS and GC was relatively weak, with the median of the Pearson correlation coefficients only slightly above 0. The positive relationship was mainly distributed in the western part of the TRHR. The Pearson correlation coefficients seem to suggest that precipitation and temperature were the main constraints for vegetation growth for the TRHR. The impacts of precipitation on vegetation growth was complicated, as the range was huge for the coefficients.
Relationship of NDVIGS with climatic variables at the grid scale. The boxplots represent 0%, 25%, 50%, 75% and 100% quantiles of Pearson correlation coefficients for all grids, and the dashed line stands for significant level at p = 0.05.
Relationship of NDVIGS with climatic variables at the grid scale. The boxplots represent 0%, 25%, 50%, 75% and 100% quantiles of Pearson correlation coefficients for all grids, and the dashed line stands for significant level at p = 0.05.
With respect to the relationship between vegetation activity and climatic variables as well as anthropogenic variables at the county scale, we perform a similar analysis as the above for the period of 1995–2014, as data for the anthropogenic variables are available only for this period. The results are generally in agreement with those obtained for the grid scale. Overall, NDVIGS was positively correlated with GP and GT (Figure 5(a) and 5(b)), and was negatively correlated with GT only in the northeastern part of the TRHR. The response of vegetation to temperature for the northeastern part was different from that of the other parts, as two clusters could be identified. This might be due to the situation that the temperature is relatively high in this region (Figure 2(e)) while water availability is limited (Figure 2(c)) (see the Discussion section below for further details). The relationship between NDVIGS and GC was roughly positive (Figure 5(c)). Besides, growing-season vegetation coverage was higher with more LP and POP across the whole region (Figure 5(d) and 5(e)). As discussed previously, the counties with the highest livestock density (Figure 2(i)) and population density (Figure 2(k)) were located in the eastern parts of the TRHR, and this spatial variability might lead to slightly different vegetation behaviour to LP and POP between the eastern and western parts.
Relationship of NDVIGS with the driving variables at the county scale. (a) NDVIGS and GP; (b) NDVIGS and GT; (c) NDVIGS and GC; (d) NDVIGS and LP; (e) NDVIGS and POP; and (f) locations of counties.
Relationship of NDVIGS with the driving variables at the county scale. (a) NDVIGS and GP; (b) NDVIGS and GT; (c) NDVIGS and GC; (d) NDVIGS and LP; (e) NDVIGS and POP; and (f) locations of counties.
In terms of spatial heterogeneity, the results also indicated different responses of vegetation to the driving variables. For the period 1995–2014, NDVIGS was negatively correlated with GP in the southeastern part of the TRHR, positively correlated with GC and POP, respectively, in the western and northeastern parts. Except for the northeastern part, the relationship between NDVIGS and GT was positive across the entire region, with a particularly significant relationship in the southeastern and southern parts of the TRHR (p < 0.05). The correlation between NDVIGS and LP was positive in 16 out of 22 counties, indicating that better vegetation coverage could have supported more livestock, as is normally the case.
Discriminating climatic and anthropogenic drivers with partial correlation analysis
To further distinguish between climatic drivers and anthropogenic drivers of vegetation change, we first detrend all the variables to focus on correlations in interannual variability over the period 1995–2014. We then calculate the partial correlation coefficient between the vegetation and climatic/anthropogenic variables, i.e. between NDVIGS and GP (RNDVI-GP), NDVIGS and GT (RNDVI-GT), NDVIGS and GC (RNDVI-GC), NDVIGS and LP (RNDVI-LP), and NDVIGS and POP (RNDVI-POP). The results show (not presented here) that, over the period 1995–2014, vegetation activity was enhanced by climate warming (RNDVI-GT = 0.43, p = 0.10) in the TRHR, whereas it was degraded by increasing grazing pressure (RNDVI-LP = −0.28, p = 0.29), with the former having a more dominant effect. Precipitation (RNDVI-GP = 0.05, p = 0.84), cloud cover (RNDVI-GC = 0.18, p = 0.51), and population (RNDVI-POP = 0.09, p = 0.75) had no significant impacts on the greening trend of vegetation for the TRHR during 1995–2014. These observations are also in good agreement with those reported by earlier studies (Piao et al. 2014; Cong et al. 2017; Chen et al. 2018).
We conduct a similar analysis at the county scale for the period 1995–2014. Figure 6 presents the partial correlation coefficients between NDVIGS and the driving factors. The partial correlation coefficients between NDVIGS and GP, controlling the effects of GT, GC, LP, and POP, showed a remarkable spatial pattern (Figure 6(a)). The results suggest that NDVIGS in the northeastern and western parts of the TRHR was positively correlated with GP over the period 1995–2014, but the link was weak (0.1< p < 0.2). However, NDVIGS in the southeastern and southern parts of the TRHR was negatively correlated with GP, but again the correlation was not significant (p > 0.2). The relationship between NDVIGS and GT also showed a spatial pattern (Figure 6(b)). A significant positive partial correlation between NDVIGS and GT existed in the southeastern and southern parts of the TRHR, while a weak negative correlation between NDVIGS and GT existed in the northeastern part. The partial correlation coefficient between NDVIGS and GC was negative in the southeastern part of the TRHR, while it was positive in the other parts (Figure 6(c)). These results are consistent with those reported in a previous study (Liu et al. 2014).
Partial correlation coefficients of NDVIGS to the driving factors. (a) RNDVI-GP; (b) RNDVI-GT; (c) RNDVI-GC; (d) RNDVI-LP; and (e) RNDVI-POP. The labels on the colour bar, R = ±0.62, R = ±0.50, R = ±0.42, and R = ±0.34, correspond to the 1%, 5%, 10%, and 20% significance levels.
Partial correlation coefficients of NDVIGS to the driving factors. (a) RNDVI-GP; (b) RNDVI-GT; (c) RNDVI-GC; (d) RNDVI-LP; and (e) RNDVI-POP. The labels on the colour bar, R = ±0.62, R = ±0.50, R = ±0.42, and R = ±0.34, correspond to the 1%, 5%, 10%, and 20% significance levels.
The effects of anthropogenic factors on vegetation variation during 1995–2014 in the TRHR were complicated. NDVIGS in the western part of the TRHR was negatively correlated with LP (Figure 6(d)), indicating that the counties therein faced an increasing grazing pressure. The eastern part of the TRHR exhibited positive RNDVI-LP values. The link between vegetation activity and grazing effects was weak (p > 0.1) for most counties. The relationship between NDVIGS and POP was positive in most areas, except for the southern part (Figure 6(e)).
To test the robustness of our analysis, we substitute CHIRPS and ERA-Interim datasets with CMA-gridded datasets and perform the same analysis. We find that the partial correlation coefficients between NDVIGS and the driving factors show spatial patterns similar to those observed with our earlier analysis (Supplementary Figures S1 and S2, available with the online Supplementary Material). Since the definition of ‘growing season’ (May–September) might also have an impact on the results, we perform the analysis also considering a slightly extended period of growing season, i.e. from April to October. The results indicate some differences in the magnitude of partial correlation coefficients, but nevertheless similar distributions (Supplementary Figure S3).
Since it is reasonable to consider that the climatic factor with the biggest absolute value of partial correlation coefficient is the primary climatic driving factor of vegetation change (Peng et al. 2015), we now extend the definition of the primary climatic driving factor to primary driving factor (since the anthropogenic factors have only less influence) and compare the absolute values of RNDVI-GP, RNDVI-GT, RNDVI-GC, RNDVI-LP, and RNDVI-POP for each county. The results indicate that the primary driving factor during 1995–2014 was temperature in most counties of the TRHR, including in the southeastern and southern parts. For counties where GP was relatively low in the western part, the primary driving factor was precipitation. Besides, the results also reveal that vegetation activity in northeastern part was mainly driven by cloud cover (solar radiation).
DISCUSSION
Trend analysis of different periods
Due to limitations in data availability for the different variables, we considered the period 1982–2014 for analyzing the trends in NDVIGS, GP, GT, and GC, but only the period 1995–2014 for detecting the trend in LP and POP. Therefore, as one would normally expect for complex dynamic systems, such as the one under consideration in this study, the trend obtained for the shorter period (1995–2014) can be a little different from the one for the longer period (1982–2014). The results obtained for the two periods reflect this to a certain extent. Although the distribution of the trend of NDVIGS for the shorter period was similar to that of the longer period, the magnitude of the trend was found to be larger. This difference could be attributed to the fact that the enhancement in vegetation activity became stronger after 1995 in the TRHR. The GP across the region was relatively stable from 1982 to 1995, but the trend accelerated after that. This observation is consistent with the one reported in a previous study that the precipitable water in the atmosphere exhibited an increase since the 1990s (Zhang et al. 2013a). Past studies have reported that the TRHR has been exhibiting a generally uniform warming trend since 1979 (Kuang & Jiao 2016). The results from our study certainly reflect this as well. The biggest difference in the results obtained based on the analysis of different periods is that GT in the northeastern part of the TRHR experienced a slightly decreasing trend over the period 1995–2014, while the area went through an overall increasing trend over the period 1982–2015. Besides, GC over the period 1982–1995 showed a significant increasing trend, while it experienced a slightly decreasing trend during 1995–2014. Therefore, the overall GC trend was positive from 1982 to 2014.
Interactions of climatic and anthropogenic variables on vegetation activity
Climatic variables and anthropogenic factors might interact with each other to impact vegetation activity, as has been revealed by several previous studies (Piao et al. 2014; Cong et al. 2017; Pan et al. 2017). In this study, we examined the effects of interactions among the climatic variables as well as those between climatic variables and anthropogenic variables. The results reveal possible reasons as to why NDVIGS responded to the driving factors in a heterogeneous manner across the TRHR.
Vegetation in the western and northeastern parts responded to GP more sensitively, while it responded negatively in the southeastern part (Figure 6(a)). This heterogeneity could have been mainly due to temperature difference. As temperature is suggested as the main constraint on vegetation activity (Kato et al. 2006), warm weather could contribute to the effect of precipitation on vegetation. This is supported by the positive relationship between RNDVI-GP and GT obtained in our analysis (Figure 7(b)). The heterogeneity could also be explained by the difference in water availability. In arid zones, water could be the main factor that put a limitation on vegetation growth, as a higher metabolic rate requires more water (Cong et al. 2017). The positive relationship might be explained by the limitation of water in the western and northeastern parts of the TRHR, as these parts belong to semi-arid climate zones (Figure 2(c)). In humid zones, where precipitation exceeds the amount the vegetation needs, more water means higher relative humidity and less solar radiation (Ukkola et al. 2016). This puts a negative impact on vegetation growth and explains why the relationship between NDVIGS to GP was negative in the southeastern part of the TRHR. This mechanism is supported by the negative relationship between RNDVI-GP and GP (Figure 7(a)).
Partial correlation coefficients of NDVIGS with GP, GT, GC, LP, and POP for the different counties in the TRHR. GP, GT, GC, LP, and POP are calculated as the mean of 20 years' (1995–2014) growing season temperature, growing season precipitation, growing season cloud cover, livestock production, and population within a county. Each point represents a county. (a) RNDVI-GP to GP; (b) RNDVI-GP to GT; (c) RNDVI-GP to POP; (d) RNDVI-GT to GT; (e) RNDVI-GT to LP; (f) RNDVI-GT to POP; (g) RNDVI-GC to GT; (h) RNDVI-GC to LP; and (i) RNDVI-LP to GP. The sub-figures are selected from Supplementary Figure S4 based on the significant level.
Partial correlation coefficients of NDVIGS with GP, GT, GC, LP, and POP for the different counties in the TRHR. GP, GT, GC, LP, and POP are calculated as the mean of 20 years' (1995–2014) growing season temperature, growing season precipitation, growing season cloud cover, livestock production, and population within a county. Each point represents a county. (a) RNDVI-GP to GP; (b) RNDVI-GP to GT; (c) RNDVI-GP to POP; (d) RNDVI-GT to GT; (e) RNDVI-GT to LP; (f) RNDVI-GT to POP; (g) RNDVI-GC to GT; (h) RNDVI-GC to LP; and (i) RNDVI-LP to GP. The sub-figures are selected from Supplementary Figure S4 based on the significant level.
RNDVI-GT was significantly and negatively correlated with GT (r = −0.60, p = 0.003, Figure 7(d)), while weakly correlated with GP (r = 0.21, p = 0.355, Supplementary Figure S4) and GC (r = 0.25, p = 0.265, Supplementary Figure S4). Past studies have reported that climate warming contributed to the greening trend in the Tibetan Plateau and in the TRHR (Piao et al. 2015; Zhu et al. 2016). Such findings support the positive RNDVI-GT value in most areas of the TRHR. When the temperature reaches above the range of vegetation acclimation conditions, the relationship between NDVIGS and GT might become negative because of extreme hot days (Piao et al. 2014; Seneviratne et al. 2014) or drought (Zhang et al. 2015). It could be further found that RNDVI-GT significantly correlated negatively with both LP (r = −0.57, p = 0.006, Figure 7(e)) and POP (r = −0.65, p = 0.001, Figure 7(f)). These results suggest that both grazing pressure and human activity offset the contribution of climate warming on vegetation growth in the TRHR, and in general.
The results also indicate that the effect of GC on vegetation was regulated by GT, as a positive relationship between RNDVI-GC and GT was found (r = 0.24, p = 0.274, Figure 7(g)). Grassland in the TRHR normally favours cloudy and warm weather during the growing season. The coordination of climatic variables and anthropogenic variables might drive vegetation towards a favourable direction. The wet and cloudy weather facilitates the grassland to flourish, thus making it possible for the grassland to raise more livestock. This mechanism is supported by the significant positive correlation observed between RNDVI-GC and LP (r = 0.63, p = 0.002, Figure 7(h)) and the positive correlation between RNDVI-LP and GP (r = 0.34, p = 0.125, Figure 7(i)). Moreover, vegetation was enhanced by warm weather in most areas of the TRHR, and population pressure on vegetation was relieved by warming climate and increasing precipitation. This is confirmed by the positive correlation between RNDVI-POP and GT (r = 0.23, p = 0.305, Supplementary Figure S4) and the positive correlation between RNDVI-GP and POP (r = 0.49, p = 0.021, Figure 7(c)).
Both the Pearson correlation analysis and the partial correlation analysis support that temperature was the dominant variable for vegetation growth in most areas of the TRHR in recent decades and that vegetation has become better over the past three decades with warming climate and increasing precipitation. However, vegetation might also be driven in an opposite direction with variations in climatic and anthropogenic variables. Increasing precipitation and cloud cover (Figure 3) could have made a positive contribution to vegetation growth in the TRHR (Figure 7(c), 7(h), and 7(i)), but the relationship between vegetation activity and temperature was weakening with an increase in temperature (Figure 7(d)). This means that the role of climate warming has become less dominant, or even negative. Moreover, the continuously increasing grazing and population have posed pressure on vegetation activity, which might offset the contribution of climate warming further.
Dynamic relationship between vegetation activity and the driving factors
We assess the dynamic relationships between vegetation activity and the driving variables by performing the partial correlation analysis with an 11-year moving window. Figure 8 presents the results from this analysis. The results reveal that RNDVI-GP, RNDVI-GT, RNDVI-GC, RNDVI-LP, and RNDVI-POP changed with time over the past two decades in the TRHR. As seen, RNDVI-GT exhibited a slightly downward trend, indicating the weakening relationship between vegetation dynamics and GT in the TRHR. This observation has also been made by a previous study (Piao et al. 2014). However, RNDVI-GP exhibited an upward trend, which means that GP became more dominant in the vegetation growth. With the warming of the TRHR, evapotranspiration increased and soil moisture decreased, and precipitation became more important in providing water for vegetation growth. RNDVI-LP experienced a significant upward trend, indicating that vegetation became better under warming climate and increasing precipitation, thus making it possible for the grassland to raise more livestock, as discussed earlier. There is no clear trend in RNDVI-GC and RNDVI-POP.
Changes in the partial correlation coefficients between NDVIGS and GP, NDVIGS and GT, NDVIGS and GC, NDVIGS and LP, and NDVIGS and POP with an 11-year moving window. A value on the x-axis corresponds to the centre of the 11-year moving window.
Changes in the partial correlation coefficients between NDVIGS and GP, NDVIGS and GT, NDVIGS and GC, NDVIGS and LP, and NDVIGS and POP with an 11-year moving window. A value on the x-axis corresponds to the centre of the 11-year moving window.
At the county scale, the relationships between vegetation and its driving variables also went through significant changes. An increasing trend in RNDVI-GP was observed in the southeastern part of the TRHR, whereas a decreasing trend was found in the northeastern part (Figure 9). An increasing trend in RNDVI-GT and RNDVI-LP was mainly observed in the northeastern part and western part, respectively. These temporal changes in the response of vegetation activity to climatic variables have also been reported in other studies, based on different periods and study areas (Zhang et al. 2013b; Piao et al. 2014; Cong et al. 2017). Although these results are only based on 20 years of data, which is too short to capture the long-term trend, they still suggest the need to exercise caution in using the results from one period to constrain the response of vegetation to its driving variables in another period (Piao et al. 2014).
Changes in the partial correlation coefficients by applying 11-year moving windows at the county scale. (a) Change in RNDVI-GP; (b) change in RNDVI-GT; (c) change in RNDVI-GC; (d) change in RNDVI-LP; and (e) change in RNDVI-POP. The change means the Pearson correlation coefficient between partial correlation coefficient and time. The labels on the colour bar, r = ±0.77, r = ±0.63, r = ±0.55, and r = ±0.44, correspond to the 1%, 5%, 10%, and 20% significance levels. The dots indicate significance at p < 0.05.
Changes in the partial correlation coefficients by applying 11-year moving windows at the county scale. (a) Change in RNDVI-GP; (b) change in RNDVI-GT; (c) change in RNDVI-GC; (d) change in RNDVI-LP; and (e) change in RNDVI-POP. The change means the Pearson correlation coefficient between partial correlation coefficient and time. The labels on the colour bar, r = ±0.77, r = ±0.63, r = ±0.55, and r = ±0.44, correspond to the 1%, 5%, 10%, and 20% significance levels. The dots indicate significance at p < 0.05.
Multiple reasons may have contributed to the variability of the response of vegetation to the driving variables in the TRHR, including climatic and non-climatic considerations. It is now clear that both climate change and anthropogenic change were responsible for this variability, and still other effects, such as frost damage during spring, fire disturbances (Huxman et al. 2004), land use and land management changes (Epstein et al. 2013), may have also played important roles.
Our results also seem to suggest that the direction of temporal changes in the relationship (between vegetation and the driving variables) was to offset the restrictions that climate and human activity put on vegetation. More specifically, the change in RNDVI-GP was negatively correlated with RNDVI-GP (r = −0.37, p = 0.093, Supplementary Figure S5a), the change in RNDVI-GT was negatively correlated with RNDVI-GT (r = −0.66, p = 0.001, Supplementary Figure S5b), and the change in RNDVI-GC was negatively correlated with RNDVI-GC (r = −0.19, p = 0.399, Supplementary Figure S5c). That is to say, the magnitude of the response of vegetation to changes was to alleviate the pressure that these changes put on vegetation. This function might be fulfilled by adjusting respiration, photosynthesis, or evaporation (Shen et al. 2015b). Under intensified climate change and human activities, species must acclimate, adapt, move, or die (Corlett & Westcott 2013). The mechanism of this acclimatization and adaption needs further investigation.
CONCLUSIONS
In this study, we analyzed the changes in the growing season vegetation, growing season climate change, and annual human activities in the TRHR. We attributed the growing season vegetation activity to climate change and human activities by applying the Pearson correlation analysis and the partial correlation analysis. We also investigated the interactions among the climatic variables as well as between climatic variables and anthropogenic variables, and explored the dynamic relationship between vegetation activity and the driving factors. The results indicated that vegetation in the TRHR experienced a slightly greening trend during 1982–2014, accompanied by significant changes in the climatic and anthropogenic variables. The southeastern and southern parts of TRHR exhibited a browning trend, while the western and northeastern parts exhibited a greening trend.
During the period 1995–2014, the dominant factor for vegetation dynamics was temperature for the southeastern and southern parts of the TRHR, precipitation for the western part, and solar radiation for the northeastern part. The regulation effects of temperature on precipitation and cloud cover contributed to vegetation growth, while grazing activity and population activity offset the contribution of warming climate. Although vegetation has become better over the past three decades with warming climate and increasing precipitation, it might be driven in the opposite direction under climate change and intensified human activities. The dynamic relationship between vegetation activity and the driving variables reflected the acclimatization and adaption processes of vegetation to climate change and human activities. This mechanism, however, needs further investigation, which we will look into in our future studies.
ACKNOWLEDGEMENTS
This paper was financially supported by the Research and Development Special Fund for Public Welfare Industry of the Ministry of Water Research in China (No. 201501028) and the National Key Research and Development Program of China (2017YFC0403600 and 2016YFE0201900). The first author gratefully acknowledges the Tsinghua Scholarship (2018015) for the financial support for a 6-month overseas graduate study visit at the University of New South Wales.
SUPPLEMENTARY MATERIAL
The Supplementary Material for this paper is available online at https://dx.doi.org/10.2166/hydro.2019.003.