The spatial-temporal distributions of controlling factors on vegetation growth in Tibet Autonomous Region, Southwestern China

Due to cold and arid climate of Tibet Autonomous Region, vegetation growth is considered to be controlled by both moisture availability and warmth. In order to reveal the patterns of regional climate change and the mechanisms of climate-vegetation interactions, long term (1982–2013) datasets of climate variables and vegetation activities were collected from Climatic Research Unit (CRU) and Global Inventory Monitoring and Modeling System (GIMMS). Principal regression analysis and (partial) correlation analysis were conducted to reveal the contributions of controlling factors on vegetation growth. Study results showed that (1) Annual mean air temperature (TMP) had increased by 0.38 °C per decade (P = 0.00) and annual precipitation (PRE) had increased by 17.25 mm per decade (P = 0.15). A significant change point around the year 1997/1998 was detected by Mann-Whitney-Pettit test, coinciding with the occurrence of El Niño event. (2) Normalized Difference Vegetation Index (NDVI) had an insignificant positive trend. Spatially, pixels of high NDVI values, great NDVI trends and high inter-annual deviations are distributed in the densely vegetated eastern part. Principal regression analysis revealed that, alpine grassland (northern and western part) is mostly controlled by temperature, steppe meadow (middle and southern part) is mostly controlled by precipitation, and shrub/mixed needle leaved and broad leaved forest (eastern part) is mostly controlled by cloud coverage. (3) Partial correlation analyses showed that regions with high sensitivity to precipitation nearly overlapped with regions of high sensitivity to minimum temperature. And the high importance of cold index (CDI, accumulated negative difference between TMP and 5 °C) revealed in this study implied the effects of regional glacial melting and permafrost degradation. We concluded that the regional climate change can be characterized as warming and wetting. Different regions and vegetation types in Tibet Autonomous Region demonstrated different driving climate factors and climate-vegetation relationships.


Introduction
The Tibetan Plateau is located in southwestern China and considered as the 'Third Pole' of the Earth. It provides various ecosystem services for human beings, including climate regulation, water resources retention, biodiversity conservation and regional/national eco-security protection. With high sensitivity to climate Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
warming and an average altitude exceeding 4000 m a.s.l., TP is also a hot spot for studies of climate-vegetation interactions (Kim et al 2012, Huang et al 2016, Cong et al 2017. In recent years, a growing number of studies on the regional climate-vegetation relationship provided controversial results. For example, studies revealed that water availability is considered as the main limiting factor for plant growth (Yang et al 2010, Kulkarni 2012, Zeppel et al 2014 and phenological changes (Shen et al 2015, Wang et al 2015 over the Tibetan Plateau. Whereas, other studies highlighted temperature as an important constraint due to the cold environments , Shi et al 2019. In addition, more studies have recognized the role of asymmetric warming, which is characterized as higher trend of diurnal minimum temperature than maximum (You et al 2016), played in controlling the plant phenology and vegetation dynamics .
Besides the controversial conclusions on controlling factors for vegetation activity, the spatial pattern of vegetation growth responding to climate change remains unclear. Different vegetation types (Zhu et al 2017) and regions (Wang et al 2015) have different pattern of climate-vegetation interactions, resulting in a considerable spatial heterogeneity of vegetation activity (Huang et al 2016).
Therefore, sensitivity analyses of vegetation dynamics to climatic variables are essential to understanding the relationship between the regional climate change and vegetation dynamics. Moreover, understanding the ecological effects of climate change on the Tibetan Plateau is of great significance to improve environmental management and to promote regional sustainable development .
With datasets collection of climate variables and vegetation activities in the Tibet Autonomous Region, this study tries to address the following questions: (1) what are the temporal patterns of changes in climate factors and vegetation activity in this region?
(2) what are the ecological effects of climate change (especially asymmetric warming) on vegetation dynamics? (3) what is the spatial distribution of the vegetation sensitivity to different climate factors?
2. Materials and methods 2.1. Study site Tibet Autonomous Region (26°50′-36°53′N; 78°25′-99°06′E) is located in the south-western part of Tibetan Plateau. With average annual temperature of less than 5°C and long sunshine duration of more than 3000 h per year, this region is well known for its harsh environment and fragile ecosystems. Frost occurs more than 6 months of the year and permafrost occurs over extensive parts of the region. The climate of the Tibetan Plateau is alternatively controlled by the Indian summer monsoon in the summer and westerlies in the winter. The total amount of annual precipitation varies from over 1,000 mm in the southeast to less than 50 mm in the northwest.
This region consists of Qiangtang Semi Arid Zone, South Tibetan Semi Arid Zone, Eastern Tibetan Semi Humid Zone and Eastern Himalayas South Slope Humid Area (Academician of Chinese Academy of Sciences 2001). The vegetation distribution follows the spatial distribution of humidity: Alpine Steppe grows in the western and northern part, Meadow Grassland grows in the middle and southern part, scrubs and mixed needle/broad leaved forest grows in the eastern part (Yang et al 2010 (Figure 1). The major vegetation type of this region is alpine grassland, served as the dominant food supply for animal husbandry of Tibetan communities. Alpine steppe consists of cold and xerophytic herbaceous plants such as Stipa purpurea and Carex moorcroftii, mixed with alpine forbs (e.g., Polygonum viviparum). Whereas, meadow grassland is dominated by perennial grasses such as Kobresia pygmaea, K. humilis and K. tibetica.

Data collection 2.2.1. NDVI data
The Normalized Difference Vegetation Index (NDVI), derived from the contrast between the red and nearinfrared reflection of solar radiation, is commonly used to reflect the growth conditions of the vegetation activity. The latest version of the global vegetation inventory modelling and mapping studies (GIMMS) NDVI3g dataset was download form https://nex.nasa.gov/nex/projects/1349/. GIMMS NDVI3g datasets have been normalized to account for issues such as sensor calibration loss, orbital drift, and atmospheric effects such as volcanic eruptions (Pinzon and Tucker 2014). We followed He et al (2017) to do the cleaning (eliminating the Flag values from 3 to 7), then re-projected this NDVI3g data onto a geographic grid with WGS 1984 spheroid. We collected 32 years  of NDVI images with a temporal interval of 15 days and a spatial resolution of 1/12°.

Climate data
The Climatic Research Unit (CRU: University of East Anglia Climatic Research Unit, 2008) provides high resolution gridded datasets with global coverage. The data sets (CRU TS4.01, https://crudata.uea.ac.uk/cru/data) provide time series for a range of climate parameters (Harris et al 2014). In the CRU processing scheme, several station-based data sources were harmonized to obtain most reliable estimates (table 1). After spatial interpolation and superposition on the reference climatologies, these anomalies constitute the final 0.5°grids global climate dataset (Mitchell and Jones 2005). We extracted the monthly data of climate variables from the year 1982 to 2013 and cropped the dataset by the regional boundaries. Then, the original CRU datasets were disaggregated using a bilinear interpolation function to match the GIMMS NDVI3g datasets spatially and temporally.
The threshold of 5°C is considered as an important bio-climatic indicator in ecological studies especially for cold environment (Dong et al 2012). The continental vegetation distribution and species spatial distribution limits are supposed to be constrained by the accumulated temperature above/below this threshold (Kira 1991). In this study, Warm index (WMI) and cold index (CDI) (Kira 1945) were prepared for each pixel using equations (1) and (2), which counts the annual sum of positive and negative differences between monthly means and 5°C. As the vegetation activity is mostly influenced by coldness before the growing season, the CDI in this study is counting the negative temperature difference (TMP-5°C) between two consecutive summers.   (Kira, (1945)) CDI Cold index of Kira, accumulate negative difference between TMP and 5°C°C/month (Kira, (1945)) NDVI Normalized Difference Vegetation Index -GIMMS ndvi3g

Trend analysis
We analyzed the trends of climatic variables and vegetation activities with the Mann-Kendall test (Mann 1945, Hamed 2008. As serial autocorrelation could influence the trend significance detected by Mann-Kendall test, a pre-whitening procedure (MK-TFPW) was conducted to reduce the influence of autocorrelation on the significance of Mann-Kendall test results (Yue and Wang 2002). Then, the Mann-Kendall test was applied to prewhitened temperature series. Positive tau values of Mann-Kendall test indicated an increasing trend, whereas negative tau represented a decreasing trend. To detect the change points in these time series, the Mann-Whitney-Pettit test was conducted (Pettitt 1979). The slope of each variable was computed by simple linear regression against the year. Since a linear regression trend fails to represent the slope when abnormal anomalies are included in the time series, then the true slope (change per unit time) can be estimated by using a simple nonparametric procedure developed by Sen (1968), which is robust against outliers and has the ability to reject anomalies without affecting the slope. The Sen slope estimator is the median of the slopes computed for n values observed at all pairwise time steps for a total of n×(n-1)/2 slopes.

Correlation analysis
Correlation analyses were conducted to reveal the relationship between climatic variables and vegetation activities. In consideration of the colinearity among the climate variables, the partial correlation coefficient between the NDVI and each climatic factor was calculated, with the other climatic factor acting as control variable. The partial correlation of time series of x i and x j given x k is: Where r ij|k is the partial correlation coefficient between x i and x j controlling x k , r ij is the correlation coefficient between x i and x j , r ik is the correlation coefficient between x i and x k , r jk is the correlation coefficient between x j and x k .

Principal regression analysis
To reveal the spatial pattern of climate-vegetation interaction and eliminate the co-linearity among climate variables, we used principal components regression (PCR) to identify the vegetation sensitivities to the correlated climate variables. In each pixel, all time series of climate variables and NDVI were transformed to z-score anomalies using long-term mean and standard deviation. For the data frame of climate variables, principal component analysis was conducted to extract the major principals and the correlated climate variables. Then the loadings of each variables were multiplied by the coefficients derived from regression analysis between standardized NDVI and the major principals. Finally, the relative importance of each climate variables in driving the changes in NDVI can be revealed and then plotted. All of the calculations were conducted in R (R Core Team 2014). Figure 2 and  2). By contrast, PRE had a moderate positive rate of 17.80 mm·10 yr −1 from 1982 to 1997, followed by great decreasing rate of −47.94 mm·10 yr −1 from 1998 to 2013. The time series of CLD is characterized by significant positive trend but insignificant change point, which is different from the temporal patterns of PRE and temperature regime. Table 3 shows that, NDVI had an insignificant Pearson correlation coefficient but a significant Kendall correlation coefficient with PRE, indicating the nonlinear relationship between PRE and NDVI. Significant Pearson correlations between temperature regime (TMP, TMX, CDI) and NDVI suggests the complicated controlling factors of NDVI inter annual variation.

Correlation and principal component analyses
Partial Pearson correlation analysis shows that the correlation coefficient (R NDVI·TMP, PRE ) between NDVI and TMP (controlling PRE) is higher than R NDVI·PRE, TMP (correlation coefficient between NDVI and PRE when TMP is controlled), suggesting the higher NDVI sensitivity to TMP (table 4). Whereas, Kendall and Spearman methods revealed a contrary result, indicating the complicated climate-vegetation relationship. Table 4 also revealed that NDVI has high sensitivities to TMX and CDI than TMN and WMI respectively, indicating the higher importance of TMX and CDI in driving NDVI inter annual variation.
Spatially, more pixels display the significant positive partial correlation between NDVI and TMP (21.62%, controlled PRE), compared with 7.35% of pixels that display the significant positive partial correlation between NDVI and PRE (controlled TMP). A proportion of 14.49% pixels displays positive partial correlation between  NDVI and CDI (controlled WMI) whereas, only 3.92% of pixels displays positive partial correlation between NDVI and WMI (controlled CDI).

Principal regression analysis
Principal component analysis revealed that, the 1st principal component is highly correlated with TMP, TMN and TMX, representing the temperature component (table 5). The 2th principal component is highly correlated with PRE, representing the humidity component. And the 3th principal component is highly correlated with CLD, representing the radiation component. The first 6 principals explained more than 99.9% of variance.
Using the principal regression coefficients as the importance, CDI has the highest importance in controlling NDVI, CLD and PRE ranks the second and the third importance respectively.

Spatial distribution of climate-vegetation relationship
Spatially, high values of NDVI and NDVI standard deviation are distributed in the eastern part covered by densely vegetation (figures 3(a), (b)). And pixels of significant positive NDVI trends are mostly located in eastern part ( figure 3(c)). The negative NDVI trends are mostly located along the middle reaches of Yalu Tsangpo River. Using TMP, PRE and CLD as the independent variables, principal regression analysis revealed that the alpine steppe of western part is mostly sensitive to temperature, the meadow grassland of middle part is mostly sensitive to humidity and the scrubs/mixed needle leaved and broad leaved forest in eastern part is relatively sensitive to radiation ( figure 3(d)).
In order to reveal the spatial distribution of relative importance that controls the vegetation growth, partial correlation analysis at pixel scale was conducted. TMP has higher importance than PRE in northern Tibet Autonomous Region (figure 4(a)). Whereas, PRE has higher importance than TMP in South Tibetan Semi Arid zone ( figure 4(b)). The impact of TMX is higher than TMN in most of the Tibet Autonomous Region, whereas, the impact of TMN is more important than TMX in Yalu Tsangpo River overlapped with the regions with high sensitivity to PRE. Region with high partial correlation between NDVI and WMI (controlled CDI) is located in northern part and forest distributed southeastern edge (figure 4(e)). Whereas, regions with high sensitivity to CDI (figure 4(f)) nearly overlaps with the regions with high sensitivity to TMP ( figure 4(a)).

Characteristics of climate change and vegetation growth
In this study, air temperature showed a strong trend of 0.38°C per decade, which was higher than the global average level (Pachauri et al 2014). Combined with the positive trend in PRE, the regional climate change can be characterized as warming and wetting (Kuang and Jiao 2016). The greater warming trends in CDI and TMN than WMI and TMX represent the widely reported asymmetric warming Chen 2000, Xia et al 2014). As diurnal minimum temperature is considered to be influenced by the surface long wave radiation (You et al 2017), the increased cloud coverage is accountable for the increased TMN and the decreased diurnal temperature range (You et al 2016).
As a result of regional climate change especially for the change of thermal growing season length (Dong et al 2012), NDVI has a positive trend over the past decades, which is in line with the previous reports of regional vegetation greening (You et al 2016) and surface albedo decreasing (Tian et al 2014). Spatially, the positive trends of NDVI are mostly distributed in eastern part of the region, which is coincided with the previous studies of regional phenology changes (Ding et al 2013) and vegetation growth .
This study also reveals the scattered pixels with negative NDVI trend (mostly distributed along the Yalunzangbu river), implying the ecological effects of land cover change, grassland degradation, urbanization, deforestation and desertification (Cui andGraf 2009, Shen et al 2012). Therefore, the ecological service might have reduced in eastern region (Tang et al 2018), despite the regional climate warming and wetting over the past decades.

Driving factors for vegetation activity
The asymmetric warming and its impact on vegetation activity has been reported by Zu et al (2018). The limited impacts of TMX and WMI on NDVI variation suggests the insignificant role of warmth in controlling the vegetation growth, which is contradicted with the warmth-limited regions such as boreal biomes (Peng et al 2013). The high importance of TMN on controlling vegetation growth revealed in this study complies with the previous report of correlation between TMN and vegetation greenness over the Tibetan Plateau (Shen et al 2016). The insignificant regional NDVI trend and complicated driving factors can be attributed to the spatial distribution of vegetation types (Zhong et al 2010). As the Tibet Autonomous Region consists of the alpine steppe in the northern and western part, the meadow grass in the middle and eastern part, and the shrub/mixed needle leaved and broad leaved forest in the eastern part, principal regression analysis revealed the 3 different patterns of controlling factors corresponding to the spatial distribution of the 3 vegetation types. Previous studies had also reported the different ecological consequences of climate change determined by different vegetation types (Sun et al 2016. In this region, observed PRE is only a fraction of moisture availability for vegetation growth. The glacier melting, degradation of permafrost were neglected in this study. Previous study showed that, winter warming and greater TMN slope were considered as the indicators of permafrost degradation and water variability (Liu et al 2011). As a result, region with significant partial correlation between NDVI and TMN (controlled TMX) overlaps with the mass loss in southeastern Tibet and along the Himalayas

Conclusion
This study reveals the regional climate warming and wetting. The significant change point around the year of 1997/1998 coincides with the occurrence of global El Niño event and the starting of weakening India Monsoon. The regional NDVI had an insignificant positive trend with an insignificant change point, indicating the complicated climate-vegetation relationship. The asymmetric warming and the high importance of CDI and TMN possibly reflect the ecological effects of glacier melting and permafrost degradation.
Spatial distribution of controlling factors represents the spatial distribution of vegetation types. Alpine steppe is sensitive to temperature, meadow grassland is sensitive to humidity, and mixed needle leaved and broad leaved forest in eastern part is relatively influenced by radiation. The scattered negative NDVI trend implies the ecological effects of land cover changes, grassland degradation, urbanization, deforestation and desertification resulting from the increased human activity. Further studies are needed to reveal the impact of the asymmetric warming on the regional carbon/water flux and carbon/water balance.