Precipitation Drives the NDVI Distribution on the Tibetan Plateau While High Warming Rates May Intensify Its Ecological Droughts

Climate change has significantly affected the ecosystem of the Tibetan Plateau. There, temperature rises and altered precipitation patterns have led to notable changes in its vegetation growth processes and vegetation cover features. Yet current research still pays relatively little attention to the regional climatic determinants and response patterns of such vegetation dynamics. In this study, spatial patterns in the response of the normalized difference vegetation index (NDVI) to climate change and its dynamic characteristics during the growing season were examined for the Tibetan Plateau, by using a pixel-scale-based geographically weighted regression (GWR) based on the Global Inventory Modeling and Mapping Studies (GIMMS) NDVI data, as well as data for temperature and moisture indices collected at meteorological stations, for the period 1982–2015. The results show the following. Spatial nonstationary relationships, primarily positive, were found between the NDVI and climatic factors in the Tibetan Plateau. However, warming adversely affected vegetation growth and cover in some arid and semiarid regions of the northeast and west Tibetan Plateau. Additionally, precipitation played a dominant role in the NDVI of the Tibetan Plateau in the largest area (accounting for 39.7% of total area). This suggests that increased moisture conditions considerably facilitated vegetation growth and cover in these regions during the study period. Temperature mainly played a dominant role in the NDVI in some parts of the plateau sub-cold zone and some southeastern regions of the Tibetan Plateau. In particular, the minimum temperature was the dominant driver of NDVI over a larger area than any of the other temperature indices. Furthermore, spatial regressions between NDVI dynamics and climatic variability revealed that a faster warming rate in the arid and semiarid regions impeded vegetation growth through mechanisms such as drought intensification. Moisture variability was found to act as a key factor regulating the extent of vegetation cover on the south Tibetan Plateau.


Introduction
Hydrothermal conditions in the climatic environment are the primary nonbiological factors that determine vegetation characteristics (e.g., phenology, productivity, and distribution patterns of plants) and their dynamic variation [1,2]. The normalized difference vegetation index (NDVI) is a commonly used index to study vegetation growth and cover, and it has been extensively applied on various regional scales [3,4]. Time-series NDVI data are often used to analyze the characteristics of changes in vegetation and their relationships with climatic factors [5]. However, while able to promote vegetation growth and development as well as the physiological and biochemical attributes of vegetation, changes in climatic conditions may also adversely affect vegetation growth and cover and display a pronounced nonlinear pattern of action [6,7]. Thus, climate change will inexorably alter vegetation habitats and thereby influence its growth conditions and, to a significant extent, its key characteristics (e.g., geographical distribution and diversity) [8]. Hence, studying the relationship between climate and vegetation can provide a theoretical basis for coping with climate change and improving and optimizing ecological functions, making it an important area of current research on global change [9].
Due to its geographical conditions-high altitudes and middle-low latitudes-and special underlying surface features (e.g., permafrost), the Tibetan Plateau, also known as the "Third Pole" of the Earth, is a region quite sensitive to global change, one that significantly influences the climate in its surrounding regions or even at the global scale [10]. Since the 1980s, the Tibetan Plateau has experienced notable warming, at a rate approximately twice the global average [11], which has had an immense impact on its terrestrial carbon sink [12]. Studies in recent years have demonstrated that the vegetation greening trend on the Tibetan Plateau is driven by climate change [13,14]. By integrating multiple models, research has shown that the contribution of climate change to gross primary productivity has surpassed that of land use and carbon dioxide concentration on the Tibetan Plateau [15]; however, this negative warming effect upon vegetation productivity generally occurs in relatively arid regions. This is because warming intensifies the evapotranspiration, which then amplifies soil moisture deficiencies [16,17]. Thus, the response of vegetation to climate change may be even more complex under the dual stresses of frigidity and aridity [18]. Evidently, examining the impact of climate change on the dynamics and regional distribution of vegetation on the Tibetan Plateau has paramount ecological value and practical significance for deepening our understanding of global climate change and protecting the ecological environment.
Current research concerned with ecosystems of the Tibetan Plateau tends to mostly focus on the pattern of change in alpine grassland and its correlation with hydrothermal factors amidst climatic warming. Nonetheless, it remains necessary to advance research into distinguishing the regional climatic determinants and the response mechanisms of vegetation dynamics and their corresponding spatial pattern on a macroscopic scale [19]. Further, the superimposition of multiple climatic factors compounds the uncertainty in the investigation of their relationships with vegetation growth and cover [20]. In this context, based on the eco-geographical regionalization of the Tibetan Plateau, this study analyzed the spatial nonstationary relationships of NDVI vis-à-vis various climatic factors and their respective variability over a recent time period (at least 30 years), by using such methods as geographically weighted regression (GWR), with two objectives in mind. Firstly, to determine the role of climatic factors in controlling alpine vegetation and its general distribution pattern, and secondly, to examine the spatial pattern of the response of vegetation dynamics to climate change from 1982 to 2015, both under the overarching aim to continuing to strengthen research on the relationships between vegetation and climate on the Tibetan Plateau.

Study Region
Situated in southwest China, the Tibetan Plateau is a vast region spanning six provinces (Tibet, Qinghai, Sichuan, Yunnan, Gansu, and Xinjiang) which, at an average altitude of over 4000 m, is known as the "Roof of the World". The Tibetan Plateau is an important ecological security barrier in China. Climatically, the Tibetan Plateau is characterized by intense radiation, low temperatures, and wide daily temperature difference (14~17 • C) and annual temperature range (−5.6~17.6 • C). Annual average precipitation in this region is mostly below 400 mm, and moisture and thermal conditions transition from warm and wet to cold and dry from southeast to northwest. Alpine grassland is the predominant type of vegetation in this region, which accounts for one third of China's total grassland area [21]. Zonally, vegetation transitions from forests to meadows, to steppes, and eventually deserts, when going from east to west on the Tibetan Plateau. For the analysis, it was divided into 10 ecoregions ( Figure 1, Table 1) according to the eco-geographical regions established for China by Zheng [22].

Climatic Variables
Two types of climatic variables, namely temperature variables (mean temperature, T ave ; maximum temperature, T max ; and minimum temperature, T min ) and moisture variables (precipitation, P; and relative humidity, RH), were used in this study. These data are available from the China Meteorological Data Service Center. We collected monthly data sets produced by 135 meteorological stations on the Tibetan Plateau from 1982 to 2015. Then, the station climate data were interpolated using the Auspline v4.2 to 8-km raster data, to match the spatial resolution of the NDVI, with a view to facilitating subsequent analyses and computations. Finally, climate data for a randomly chosen period, collected at the meteorological stations not previously used (because of missing data for some months) in the interpolation procedure, were selected to compare with and validate the values in the corresponding interpolated grid cells. The correlation coefficients (Pearson's r) between the interpolated data for each of the three temperature indices and their measured values were >0.99, while those for P and RH were 0.86 and 0.89, respectively. Hence, the interpolation accuracy was relatively high and thus met the requirements for a robust study of their overall patterns on the Tibetan Plateau. After these treatments, the monthly temperature variables and relative humidity during the growing season (from May to September) were averaged to obtain interannual data, and the annual total precipitation was obtained by summing the monthly data of the growing season.

Normalized Difference Vegetation Index (NDVI)
The NDVI dataset (temporal resolution: 15 days; spatial resolution: 8 km) for the 1982-2015 period provided by the Global Inventory Monitoring and Modeling Studies (GIMMS) is used here to derive vegetation growth and cover indices. Characterized by relatively high accuracy and a long time series, the GIMMS NDVI dataset has been extensively used to study global and regional large-scale changes in vegetation [23][24][25].
Monthly NDVI values were then calculated using the maximum value composite (MVC) method. On this basis, monthly average NDVI values for the growing season (based on previous studies obtained for the Tibetan Plateau, the growing season was defined here as the period from May to September each year) were obtained. In addition, regions with annual average NDVI values <0.1 were removed, to eliminate the effects of non-vegetation factors [26]. It should be noted that non-vegetation area accounted for a large proportion in some arid regions (ER1D1, ER2D1, and ER2D2), which were not considered in this study. Subsequently, the spatial relationship between the monthly average NDVI and climate variables was analyzed by the methods mentioned in the following sections.

Spatial Autocorrelation
The precondition for using geographically weighted regression (GWR, see Section 2.6) is that geospatially correlated independent and dependent variables differ spatially, i.e., the regression coefficient varies with the spatial location of the independent variable [27]. The Moran's I index is a measure of global spatial autocorrelation that reflects the extent of the similarity between the attribute values of a spatial variable in adjacent spaces. The Moran's I is calculated this way: where W ij is the spatial connection matrix between points i and j; the n is the total number of spatial locations, and the x i and x j are the spatial attribute values of points i and j, respectively, while x is the average value of the variable in all the spaces. The value of I ranges from -1 to 1. When I > 0, it means the variable is positively spatially autocorrelated. Conversely, when I < 0, the variable is negatively spatially autocorrelated. When I = 0, however, the variable is not spatially correlated.
Usually, the standard Z-statistic is used to examine the significance of the Moran's I: where E I and var(I) are the mean and variance of the Moran's I, respectively. A positive and statistically significant Z-score indicates a positive spatial autocorrelation and that similar attribute values tend to cluster spatially. A negative and statistically significant Z-score indicates a negative spatial autocorrelation and that similar attribute values tend to disperse spatially. Accordingly, Z-score of 0 indicates that the attribute values of the spatial variable are distributed independently and randomly.

Trend Analysis
In this study, the interannual trends of change in the growing-season climate and vegetation indices were analyzed using the grid-scale-based ordinary least-squares (OLS) method [28]. The equation for trend analysis is given by: where θ slope is the slope of the linear regression, which conveys the trend and rate of change in the climate or vegetation index for the Tibetan Plateau, Y i is the value of the climate or vegetation index of the ith year, and n is the total number of years examined in this study (n = 34).

Geographically Weighted Regression (GWR)
GWR, originally proposed by Brunsdon et al. [29], is a method used to extend ordinary linear regression (e.g., OLS) when carrying out local spatial analyses. Using this method can clarify changes in the spatial relationships within a given region, for better understanding. The GWR's parameters are a function of spatial location, in this way: where y i , x ik , and ε i are, respectively, the dependent variable (NDVI/NDVI variation), independent variable (climate variable/variability), and random error at the spatial point i; the (µ i , v i ) is the spatial location of point i; the k is the number of independent variables; the β k is the regression coefficient at point i; and the β 0 is the intercept. A Gaussian model is generally used as the weight function for the GWR method, as follows: where ω ij is the weight of observation point j of point i; the d ij is the Euclidean distance between points i and j; and the b is the bandwidth, which is a function that describes how weight and distance are related. When d ij > b, the ω ij = 0; conversely, when d ij = 0, the ω ij = 1. Here, the corrected Akaike information criterion (AICc) was employed to evaluate model complexity and accuracy, as well as to determine an optimum b value. Lower AICc values indicate better-performing models for a given simulation. Generally, when the difference of AICc values between two models is greater than 3, the model with the lower AICc value has the optimum b value [30].
Here, GWR considered both the spatial distribution and the dynamics patterns to reveal the spatial variations in relationships between NDVI and climate variables. First, the annual average climatic variables and NDVI should be conducted with normalization by using the maximum-minimum standardized treatment [31], making the scopes from 0 to 1, and then the normalized climatic variables and NDVI were used as independent and dependent variables, respectively, to calculate the regression coefficients, which indicated the spatial relationship between NDVI and climatic factors. The purpose of normalized coefficients was to compare the results from climate variables in different dimensions, and this treatment did not affect the spatial relationship between climate and vegetation. Second, the climatic determinants for NDVI can be obtained by comparing the regression coefficients for different climate variables, combined with the NDVI and climate trends. If variations in both temperature and moisture factors coupled with relevant spatial relationship were inconsistent with the varying trend of NDVI (that is, θ slope_climate ·β k ·θ slope_NDVI < 0), nonhydrothermal factors (such as radiation, terrain and grazing) could be the dominant driver of vegetation growth and cover. At last, the interannual variability of climate variable (hereafter, climate variability) and NDVI variation during the study period were used in the GWR as independent and dependent variables, respectively, to evaluate the spatial correlation of climate change and vegetation dynamics.

Nonstationarity test of NDVI and Climate Change
First, the GWR coefficients of the mean and variability of each climatic factor were subjected to a spatial autocorrelation analysis, to determine whether their relationships of NDVI and its dynamics to climate change are spatially stationary (Table 2). These results showed that the autocorrelation indices of the GWR coefficients were all > 0.8. This suggested that the selected coefficients were highly spatially autocorrelated; that is, they are spatially nonstationary. Additionally, each coefficient had a Z-score well above 2.58, which is statistically significant at an alpha level of 0.01. Thus, GWR can comprehensively display the quantitative effects of the climatic factors on the NDVI and its dynamics in different regions.

Spatial Nonstationary Relationship Between NDVI and Climatic Factors
It can be seen from the interannual variability of climate variables during the growing season of the study period ( Table 3) that the Tibetan Plateau had experienced a significant warming with a rate of 0.42 • C/decade. In arid and semiarid regions, the warming rates were faster, and the T min varied more obviously. Precipitation also showed an increasing trend. Areas with a faster increase rate of precipitation appeared in the semiarid region of the plateau sub-cold zone, but the precipitation rate in the sub-humid and arid regions increased slowly. The RH trend was nearly constant throughout the Tibetan Plateau during the study period, and all ecoregions except the semiarid region of the plateau sub-cold zone showed a downward trend. Spatially, the three temperature variables showed increasing trends in almost all regions of the Tibetan Plateau, and the precipitation in most regions was also an upward trend (Figure 2b). The RH trend, however, was decreasing in the southeastern parts of the Tibetan Plateau (Figure 2c). Table 3. Climatic variables and trends in ecoregions of Tibetan Plateau during the growing season of 1982-2015.  The GWR results for the NDVI and temperature indices revealed a positive spatial correlation in most regions of the Tibetan Plateau (Figure 3). This indicated that warming, to a certain extent, facilitated vegetation growth on the Tibetan Plateau. A negative spatial correlation was found between the NDVI and each of the three temperature indices in some arid and semiarid regions (e.g., ER1C1, ER1C2, ER2C1, and ER2D3). Temperature rises accelerated the evapotranspiration of soil moisture and further increased the extent of droughts and, as a result, reduced the vegetation cover. This phenomenon was particularly prominent in region ER2C1 and the northeastern part of region ER1C1.

T ave T max T min P RH
The GWR results for the NDVI and moisture indices showed a positive spatial correlation between the NDVI and precipitation in most regions of the Tibetan Plateau (Figure 4a), so an increase in precipitation could increase the NDVI in these regions. A negative spatial correlation between the NDVI and RH was found in some southern and northwestern regions of the Tibetan Plateau (Figure 4b). An increase in RH might result in stomatal closure in vegetation and weaken its transpiration [32] and, as a result, inhibit its growth and development. This phenomenon was the most prominent in some parts of southern regions, namely, ER1B1, ER2AB1, and ER2C2.

Climate Determinants for NDVI in Different Ecoregions
The spatial domain of climatic determinants for NDVI is shown in Figure 5. In terms of the overall pattern, precipitation was predominantly involved in NDVI, accounting for 39.7% of the Tibetan Plateau's total area, a proportion larger than that of any other index (Table 4). In the period of more than 30 years examined in this study, the NDVI increased mainly in those regions where precipitation increased prominently (Figure 2a). Thus, the increase in precipitation played a vital role in fostering a gradual increase in vegetation cover. RH was found to have played a dominant role in the NDVI in the western part of region ER1B1, the eastern part of region ER2AB1, and some parts of region ER1C2. In addition, a negative spatial correlation was mostly found between RH and NDVI. Thus, at least spatially, vegetation growth in these regions was hindered by RH.  The regions of dominant effect of temperature on the NDVI were mainly distributed in some parts of the plateau sub-cold zone (ER1) and region ER2AB1 of the Tibetan Plateau. Among the three temperature indices, the effect of T min was the most pronounced, influencing NDVI over the largest area (7.1%) of the Tibetan Plateau (Table 4). This dominant role of T min was most prominent in the southwestern part of region ER1C1 and the northeastern part of region ER1C2. A positive correlation was found between the NDVI and each of the temperature indices in these regions; this suggested that, at least spatially, gradual temperature rises could foster greater vegetation cover.

Response Pattern of the NDVI Variation to Climate Variability
From the GWR results for climatic variability and NDVI dynamics, a negative spatial correlation between NDVI dynamics and the variability of each of the three temperature in-dices was mainly in the semiarid regions (e.g., ER1C1, ER1C2, and south ER2C1) ( Figure 6). The NDVI mostly increased in all but some southern parts of region ER2C1 (Figure 2a). In other words, the NDVI increased slowly in these regions where temperature increased rapidly. In the western and southern parts of region ER1C1, the NDVI was predominantly affected by T min (Figure 5). In most of the other regions, the NDVI was predominantly affected by moisture.
The NDVI mostly increased in the eastern part of region ER1B1 and region ER2D3 (Figure 2a). Spatially, NDVI variability was positively correlated with T max variability yet negatively correlated with the respective variability of T ave and T min in the eastern part of region ER1B1 ( Figure 6). Thus, the adverse effects of temperature variability on vegetation growth and cover were primarily due to T min variability. Spatially, in region ER2D3, NDVI variability was also found positively correlated with T min variability and negatively correlated with the variability of T ave as well as that of T max (Figure 6). This suggested that the adverse effects of temperature variability on vegetation growth and cover were primarily a result of T max variability. Both regions were characterized by increased precipitation and decreased RH (Figure 2b, c). Spatially, a negative correlation was found between NDVI variability and each of the moisture indices in the eastern part of region ER1B1 (Figure 7). In other words, a high rate of decrease in RH led to a high rate of increase in the NDVI. By contrast, a positive correlation was found between the NDVI variability and the variability of each of the moisture indices in region ER2D3 (Figure 7). Apparently then, climate change was able to alter vegetation growth and cover in semi-humid and semiarid regions via differing mechanisms.
The variability of each climatic factor and NDVI variability were found to display the same spatial relationship in the southern regions (ER2AB1 and ER2C2) of the Tibetan Plateau. In addition, the variability of each climatic factor tested except RH was positively correlated with NDVI variability in these regions (Figures 6 and 7). This suggested similarities in the mechanisms by which climate change impacts vegetation growth and cover in these regions. In these two regions, the NDVI was predominantly affected by precipitation ( Figure 5). In region ER2AB1, there was a decrease in precipitation (Figure 2b) but a relatively marked decrease in its NDVI (Figure 2a). In other words, a high rate of decrease in precipitation led to a high rate of decrease in the NDVI. An opposite pattern characterized region ER2C2: there, both its precipitation and NDVI have increased (Figure 2a and b). Put differently, a faster rise in precipitation levels fostered a high rate of increase in the NDVI. Arguably, changes in precipitation were the principal driving factor for differential vegetation growth and cover between these regions.

Discussion
Overall, the climate on the Tibetan Plateau is characterized by low temperatures, large temperature differences, intense radiation, and high precipitation in summer [33]. In addition, those effect regions of climatic determinants on NDVI on the Tibetan Plateau in this study are consistent with those found across China (at a higher scale and lower resolution) [34]. As temperature rises, the response of vegetation growth to temperature becomes increasingly prominent [35,36]. Before reaching the optimum temperature for photosynthesis, temperature rises will facilitate photosynthesis [37], which we found in some parts of the plateau sub-cold zone (ER1). Beyond this optimum temperature, warming will, on the one hand, improve vegetation respiration and accelerate nutrient consumption (e.g., the northeast Tibetan Plateau) and, on the other hand, accelerate moisture evapotranspiration and reduce the accumulation of organic matter content [38], as it did in the southeast Tibetan Plateau. Changes in moisture can, to a certain extent, modify vegetation growth and cover in the largest areas of the Tibetan Plateau during the study period. However, more moisture can also hinder vegetation growth and reduce vegetation cover by increasing cloud cover and RH [32], which could be found in some sub-humid regions (e.g., ER1B1). As many studies have shown, the response of NDVI on the Tibetan Plateau to climate change follows a nonlinear process and generates a composite effect [39][40][41]. Consequently, this necessitates an examination of the spatial distribution and temporal dynamic patterns of vegetation growth and cover through multifactor comprehensive analysis, and provides the theoretical reference for natural restoration of alpine ecosystem.
Regression analysis of climatic variability is crucial for helping to quantify the NDVI's response to climatic factors, since by conveying the dynamic responses of vegetation growth and cover to climate change, it can better explain the dynamic relationships that temperature and moisture have with vegetation. The intensity of vegetation growth depends on the rates at which photosynthesis and respiration respond to climatic factors [42]. A positive correlation of climatic variability arises when the increase in the photosynthetic rate surpasses that in the respiratory rate [43]. Specifically, an increase in T min variability can have accelerated the consumption of biomass by increasing the nighttime assimilation rate of vegetation. This is perhaps why T min became the principal factor inhibiting vegetation growth in some regions (e.g., some parts of ER1C1). An increase in T max variability, however, can have inhibited vegetation growth and reduced vegetation cover by strengthening the daytime respiration of vegetation and weakening its photosynthesis (e.g., ER2D3). In addition, a faster warming rate accelerates moisture loss from soil and vegetation and reduced regional vegetation cover [17], as it did in some southern parts of the Tibetan Plateau. Evidently, analyzing the variation in NDVI and climatic variability as objects of investigation can better capture and convey the relationship between vegetation dynamics and climate change, as well as accentuate the impact of climate change upon mechanisms of vegetation growth and cover on the Tibetan Plateau. It also provides new ideas for the further processing of remote sensing data.
While climate is the principal driving factor affecting both vegetation growth and cover on the Tibetan Plateau, the NDVI is also influenced by non-climatic factors [44], which might appear in some humid/sub-humid regions ( Figure 5; Table 4). Of the factors contributing to spatial heterogeneity of the NDVI on the Tibetan Plateau, only climatic factors were taken into consideration in this study. In reality, a multitude of other environmental factors, including terrain [45] and soil conditions [46], will also shape the spatial heterogeneity of vegetation growth and cover. In comparison with other regions, the Tibetan Plateau' vegetation has remained relatively intact, that is insignificantly affected by human activity [47]. The GWR used in this study is a local range-based spatial regression model. Human activity displays certain similarities in the neighborhood. Thus, this statistical method can be employed to uncover and establish the effects of climatic factors on the spatial heterogeneity of NDVI.

Conclusions
In this study, spatial nonstationary relationships were simulated between the NDVI and climatic factors in the Tibetan Plateau. Then, based on the regressions between NDVI dynamics and the respective variability of moisture and temperature indices, the leading climatic factors and their mechanisms of action were established for different regions of the Tibetan Plateau. Ample moisture and heat enabled a gradual increase in both vegetation growth and cover on the Tibetan Plateau. Spatially, temperature and precipitation are positively correlated with NDVI in most regions. However, adverse effects of warming upon vegetation growth and cover have begun to emerge in a concentrated manner in some arid and semiarid regions of northeast and west Tibetan Plateau. Precipitation plays a dominant role in the NDVI in nearly 40% of the Tibetan Plateau, an area larger than any of the other indices. Temperature mainly governs the NDVI in some parts of the plateau sub-cold zone (ER1) and some southeastern regions of the Tibetan Plateau. Among the temperature indices, T min has the greatest influence on NDVI over the largest area of the Tibetan Plateau (7.1%).
The increase in the warming rate mainly inhibits the NDVI in the semiarid regions of the Tibetan Plateau. Furthermore, different temperature indices are capable of impeding vegetation growth and reducing vegetation cover by means of intensifying respiration and reducing the moisture available. This explains the difference in the climatic determinants found among regions. Spatially, temperature and precipitation variability are both positively correlated with NDVI dynamics on the south Tibetan Plateau. This suggests similarities in the mechanisms of moisture and thermal conditions acting upon the NDVI in plateau temperature zone (ER2). Looking ahead, further adjustments to the combination of moisture and thermal conditions could substantially affect vegetation growth and cover in this region.