Climate change-induced greening on the Tibetan Plateau modulated by mountainous characteristics

Global terrestrial vegetation is greening, particularly in mountain areas, providing strong feedbacks to a series of ecosystem processes. This greening has been primarily attributed to climate change. However, the spatial variability and magnitude of such greening do not synchronize with those of climate change in mountain areas. By integrating two data sets of satellite-derived normalized difference vegetation index (NDVI) values, which are indicators of vegetation greenness, in the period 1982–2015 across the Tibetan Plateau (TP), we test the hypothesis that climate-change-induced greening is regulated by terrain, baseline climate and soil properties. We find a widespread greening trend over 91% of the TP vegetated areas, with an average greening rate (i.e. increase in NDVI) of 0.011 per decade. The linear mixed-effects model suggests that climate change alone can explain only 26% of the variation in the observed greening. Additionally, 58% of the variability can be explained by the combination of the mountainous characteristics of terrain, baseline climate and soil properties, and 32% of this variability was explained by terrain. Path analysis identified the interconnections of climate change, terrain, baseline climate and soil in determining greening. Our results demonstrate the important role of mountainous effects in greening in response to climate change.


Introduction
Understanding changes in vegetation greenness in mountainous areas is of great importance because these changes are associated with a series of ecosystem processes and services, such as carbon and water cycling (Gottfried et al 2012, Hagedorn et al 2019. Changes in vegetation greenness over time usually consist of an alternating sequence of greening and/or browning, which are defined as statistically significant increases and decreases in vegetation greenness over a period of several years (De Jong et al 2012, Chen et al 2019a. The effects of climate change, notably global warming or wetting, on vegetation greenness have been widely recognized (Zhu et al 2016, Du et al 2018, Park et al 2020. However, the effects of climate change on greening trends are inconsistent depending on the mountainous characteristics, such as topography, climate and soil properties, which may interact with climate change to regulate vegetation greenness dynamics (Pugnaire et al 2019, Xu et al 2020. Identifying the processes underlying the spatial variability of greening under climate change is challenging, particularly in mountainous areas.
The Tibetan Plateau (TP), which is called the 'Third Pole' of the earth with an average elevation of 4000 m (Qiu 2008), has responded more sensitively to climate change than other zones at the same latitude (Duan et al 2017, Teng et al 2018. Over the past decades, satellite-derived vegetation indices have indicated widespread change on the TP (Shen et al 2015a, Zhu et al 2016. The longterm changes in vegetation greenness over the TP are recognized to be driven primarily by climate change (Piao et al 2014, Zhu et al 2016, specifically warming (Qiu 2008) and wetting (Chen et al 2013). However, the trend of greening shows great variability on the TP due to the strong vertical zonality of the topography (Anderson et al 2020), baseline climate (Shen et al 2011), and soil conditions (Kato et al 2006). The terrain is a significant regulator of energy and water availability even at fine scales and is very complex on the TP; thus, the terrain may play a vital role in regulating vegetation dynamics (Bolton et al 2018, Gao et al 2019. The baseline climate plays a key role in the surface energy balance and in modulating the land-climate interactions (Forzieri et al 2017), and may therefore impact the widespread greening of the TP. In addition, the heterogeneity of soil nutrient availability and the abiotic soil conditions may lead to profound changes in plant species diversity and community structure (Pugnaire et al 2019), ultimately driving diverse greenness dynamics on the TP. The comprehensive impact of climate change, terrain, baseline climate and soil, which influence the permafrost and groundwater, the flow of rainwater into or off the soil, and the heat absorbed by soil, induce extreme local variability of the vegetation composition and structure on the TP Jin 2013, Mu et al 2016). Thus, investigating the greening pattern on the TP and the underlying drivers is important for sustainably managing the TP ecosystems and is also important for the land surface energy balance of the region (Zhang et al 2013. However, it is still challenging to understand the underlying mechanisms underpinning the spatial variability of greening on the TP.
Long-term remotely sensed measurements of the normalized difference vegetation index (NDVI) have been widely used to indicate changes in large-scale vegetation dynamics and greenness (Huete 2016, Myers-Smith et al 2020. Here, we investigate the trend of the growing season NDVI across the TP for the period 1982-2015 using remotely sensed data sets. We present a new method to spatially resample the NDVI record and thus provide a new growing season NDVI data set at a 1 km resolution. Combined with climate records from meteorological stations across the TP, we also disentangle the relative importance of, and interactions among, climate change variables, terrain attributes, baseline climatic characteristics, and soil properties in controlling the spatial variability of the greening trend.

Satellite-derived NDVI data sets
Two remote sensing data sets, AVHRR GIMMS NDVI3g (hereafter NDVI GIMMS ) during the period from 1982 to 2015 and SPOT-VGT NDVI (hereafter NDVI SPOT ) covering the period from 1999 to 2013, were used in this study. NDVI GIMMS is the longest time series NDVI data set and provides information on terrestrial vegetation changes with a spatial resolution of 1/12 • (∼8 km) . This data set was compiled by merging segments (data strips) during a half-month period using the maximum value composite (MVC) method. The NDVI SPOT product has been pre-processed by the VEGETA-TION Processing Centre at the Flemish Institute for Technological Research in Belgium. The spatial resolution of NDVI SPOT is 1 km, and the temporal resolution is 10 days (Maisongrande et al 2004). To ensure temporal overlap and equal treatment with NDVI GIMMS , the 10 day images were aggregated into monthly MVC.
In this study, we used an empirical orthogonal teleconnections (EOT) model (Appelhans et al 2015, Detsch et al 2016 to assess the dynamics of the NDVI during 1982-2015 at a 1 km spatial resolution (SI appendix, Processing of NDVI data sets (available online at stacks.iop.org/ERL/16/064064/ mmedia)). After that, the new NDVI data set (hereafter NDVI EOT ) was used in the following analysis. In this study, the EOT analysis was performed using the package remote in R 3.6.3 (R Core Team 2020). To minimize the impacts of snow and ice, we considered NDVI EOT within the plant growing season from May to September during 1982-2015 (Shen et al 2015b). The mean values of the monthly NDVI EOT (hereafter NDVI mean ) during the growing season over 1982-2015 were calculated based on the corresponding NDVI EOT values to represent the average of growing season NDVI. To reduce the noise from non-vegetation signals, grid cells with longterm NDVI mean values less than 0.1 over the period 1982-2015 were excluded in this study .

Estimation of long-term NDVI trend
To quantify the long-term trend of growing season vegetation growth, we performed least squares linear regression analysis using the annual growing season NDVI mean as the dependent variable and the year as the independent variable over the period 1982-2015. The slope of the regression (hereafter Tr mean ) was then defined as the linear trend of the vegetation growth over the growing season. Positive and negative values of Tr mean over time can be referred to as greening and browning, respectively. The spatial patterns of the temporal trends in each data set were calculated using least squares linear regression for each grid. All the data calculations were accomplished using the function lm() in R 3.6.3 (R Core Team 2020).

Identifying possible driving factors 2.3.1. Climate variables
There were 64 meteorological stations, which contained daily climate variables of precipitation, temperature and evaporation, provided by the China Meteorological Administration over the TP during the study period 1982-2015. Among these meteorological stations, four had a corresponding NDVI mean less than 0.1 during the study period, and were excluded from this study. Thus, 60 meteorological stations were used in the following analysis (figure 1).
In this study, ten baseline climate variables, which determined a reference point for climate change based on the average climate over a 33 year period  based on monthly precipitation, temperature and evaporation data sets, were calculated from daily records from 60 meteorological stations over the TP. They included the average annual total precipitation (P re ), average total precipitation in the growing season (P re_gs ), average annual mean temperature (T mp ), average mean temperature in the growing season (T mp_gs ), average annual minimum temperature (T mn ), average minimum temperature in the growing season (T mn_gs ), average annual maximum temperature (T mx ), average maximum temperature in the growing season (T mx_gs ), average annual total evaporation (E vp ), and average total evaporation in the growing season (E vp_gs ). The corresponding ten climate change variables, which represented the long-term trend of climate dynamics (Tr Pre , Tr Pre_gs , Tr Tmp , Tr Tmp_gs , Tr Tmn , Tr Tmn_gs , Tr Tmx , Tr Tmx_gs , Tr Evp , Tr Evp_gs ), were also calculated (table S1). All ten climate change variables were quantified by performing least squares linear regression analysis using climate data sets as the dependent variable and year as the independent variable over the period 1982-2015. The regression slopes of the linear regression models were identified as the trends of the corresponding climate change variables (table S1).
Since the ten baseline climate variables and ten climate change variables were highly correlative (figures S2(a) and (d)), we used principal component analysis (PCA) to summarize/reduce the informative features and control for autocorrelations between these variables. We performed PCA on all variables reflecting the baseline climate variables (n = 10) and climate change (n = 10) to analyze the variation in baseline climate (PCb) and climate change (PCc). Based on the PCA results, three PCbs (PCb1 to PCb3) that explained 97.25% of the baseline climate variance and six PCcs (PCc1 to PCc6) that explained 96.08% of the variance in the assessed climate change factors were used in the following analysis (figure S2).

Terrain
At each of the 60 meteorological stations, seven terrain attributes, including the digital elevation model (DEM), slope, aspect, curvature (Curv), topographic wetness index (TWI), multiresolution index of valley bottom flatness (MrVBF), and terrain ruggedness index (TRI), were extracted from the 3 arcsecond (approximately 90 m) grid Shuttle Radar Topography Mission DEM (SI appendix, Explanation of terrain attributes).

Soil attributes
Maps of bulk density (BD), total nitrogen (TOTN), soil content (sand, silt and clay), effective cation exchange capacity (ECEC), total exchangeable bases (TEB), and available water capacity (TAWC) were obtained from a 30 arcsecond raster database of the harmonized soil property values for the world (WISE30sec) (Batjes 2016). Soil maps of pH (Chen et al 2019b) and soil organic matter (SOM) (Liang et al 2019) were produced using soil samples, environmental variables, and digital soil mapping techniques. In this study, the ten soil attributes mentioned above were extracted at the 60 meteorological station locations for the following analysis.

Ecosystem types
The ecosystem types were identified according to the vegetation map of the TP that was obtained from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (www.resdc.cn). The vegetation map was produced from a digitized 1:1000 000 Vegetation Atlas of China (Chinese Academy of Sciences 2001). In the vegetation map, four vegetation types related to forest (coniferous forest, coniferous and broad-leaved mixed forest, broad-leaved forest, and shrub) were classified as woodland ecosystems; three vegetation types belonging to grassland (steppe, herbs and meadow) were classified as woodland ecosystems. Other vegetation types, including alpine vegetation, cultivated vegetation desert and others, were classified as other ecosystems (figure 1).

Disentangling the contributions of driving factors on the spatial variability of the greening trends
We explored the variability of the relationships between indicators (i.e. climate change, terrain, baseline climate and soil) and the response variable (Tr mean ) by using a linear mixed-effects model (LMM) (Qie et al 2017) at the 60 meteorological stations. The LMM contains two components: a fixed effect (the explanatory variables) and random effects. Hypotheses were made prior to assessing the contributions of climate change, terrain, baseline climate and soil to the greening variation. First, we assumed that climate change impacted greening. Second, we assumed that the impact of climate change on the spatial variability of the greening trend was modulated by terrain, baseline climate and soil attributes. Third, we assumed that terrain characteristics amplified the spatial variability of the greening.
We built three sets of LMMs to examine these hypotheses. The first set of LMMs used the six climate change variables (PCc1 to PCc6, figure S2) as explanatory variables and Tr mean as the response variable to describe the response variable of the greening trend. The second set of LMMs used the variables of climate change (PCc1 to PCc6), terrain, baseline climate (PCb1 to PCb3) and soil as explanatory variables. The third set of LMMs was used similarly to the second set of LMMs, but without the terrain attribute.
Ecosystem types were treated as a random effect, allowing each ecosystem type to have a separate intercept in all sets of LMMs. The predicted greening rates by the three sets of LMMs driven by the identified variables mentioned above were compared with the corresponding modeled greening rates. The LMM enabled the estimation of the assessed RMSE of an explanatory variable (i.e. climate change, terrain, baseline climate and soil in this study) compared with other variables that were considered. In this study, we used the RMSE value to evaluate the relative importance of a given variable in the LMM (SI appendix, Relative importance of variables), and also used the resulting coefficient of determination (R 2 ) to present the explanation rate of the model. The LMMs were used at the locations of the 60 meteorological stations using the package nlme in R 3.6.3 (R Core Team 2020).

Partial least squares-path modeling (PLS-PM)
We applied PLS-PM with four latent variables, namely, climate change, terrain, baseline climate and soil, to assess their direct/indirect effects on the variability of Tr mean . The latent variables were reflected by indicators. For example, for the latent variable 'climate change,' we considered six indicators: PCc1 to PCc6. In the PLS-PM, the loading of each latent variable was the key to estimating the variable scores and was calculated as the correlation between a latent variable and its indicators. An iterative algorithm was used to estimate the loadings until the convergence of the loadings was reached to maximize the explained variance of the dependent variables (both latent and indicator variables) (Luo et al 2019). The R 2 values of the endogenous latent variables indicated the amount of variance in the endogenous latent variable explained by its independent latent variables. The result of the path coefficient can be positive or negative. A negative path coefficient means that increasing the causal variable causes a decrease in the dependent variable if all other causal variables are held constant.
Using PLS-PM, we explored how the greening trend responded to climate change, terrain, baseline climate and soil, and we also investigated the interactions among these variables. We first assumed that climate change, terrain, baseline climate and soil directly impact the spatial variability of the greening trend. Second, we assumed that these variables have indirect impacts on the spatial variability of the greening trend through their interactions with each other. We also tested the assumptions based on the variables without terrain (climate change, baseline climate and soil). The PLS-PM was performed using the package plspm in R 3.6.3 (R Core Team 2020). Figure 2 shows the long-term annual total precipitation, average mean temperature and annual total evaporation as well as their corresponding trends during 1982-2015. During the past three decades, the TP has experienced a warming trend based on the comparison of time series annual average temperature on the whole TP (0.05 • C yr −1 , p < 0.0001), woodland ecosystem (0.04 • C yr −1 , p < 0.0001) and grassland ecosystem (0.07 • C yr −1 , p < 0.0001) (figure 2). Based on the records of the 60 meteorological stations, the total annual precipitation showed interannual variation with an insignificant trend during 1982-2015 (0.13 mm yr −1 , p = 0.8206). Woodland ecosystem, however, showed a significant decreasing trend of annual total precipitation over the study period (−1.67 mm yr −1 , p = 0.0415). The annual total evaporation showed a universal decreasing trend over the whole TP, woodland ecosystem, and grassland ecosystem (figure 2).

Greening on the TP
The growing season NDVI EOT of the time series grid cells at a 1 km resolution shows a mean value of NDVI mean of 0.343, with an increasing trend of 0.011 per decade (i.e. greening, p < 0.0001) across the whole TP during the period 1982-2015 ( figure 3(c)). In general, the NDVI increases from the very dry western TP to the wetter eastern TP. At the pixel level, NDVI mean shows great spatial variability, ranging from 0.10 to 0.83 ( figure 3(a)). A 91% proportion of the vegetated land over the TP exhibits greening (i.e. Tr mean > 0, p < 0.0001), and only 9% experiences browning (p < 0.0001), that is, a decrease in NDVI mean ( figure 3(b)). Woodland ecosystems, which grow at lower altitudes (1706 ∼ 4440 m), have a higher NDVI mean value (0.514) (figure 3(d)) and show more rapid and significant greening rates (0.013 per decade, p = 0.0003) than the average value for the whole TP. Grassland ecosystems, which account for 45% of the total TP (figure 1(c)), show a lower greening rate than the whole TP (0.010 per decade, p < 0.0001), with an NDVI mean value of 0.286 ( figure 3(e)).

Association of greening with climate change only
The LMM results suggested that climate change explains only 26% of the spatial variability of greening trends (R 2 = 0.26) ( figure 4(a)). The random intercepts for woodland and grassland ecosystems were 0.08 (×10 −10 ) and 0.04 (×10 −10 ), respectively (table S2). Among the indicators of climate change, PCc4, which had a high correlation with Tr Tmx_gs over 1982-2015 based on the PCA results (correlation = 0.55, figure S2(f)), showed significant correlation with the greening trend (coefficient = −3.23, table S2).

Greening modulated by terrain, baseline climate and soil properties
In addition to climate change, we introduce terrain, baseline climate and soil variables (table S1) into the LMM to test the hypothesis that these factors modulate the impact of climate change on the spatial variability of greening rates. We find that variables of terrain, baseline climate and soil can explain another 58% of the greening trend variability (R 2 = 0.26 for climate change alone versus R 2 = 0.84 for all variables together). Among these variables, terrain can help explain more than 32% of the variability (R 2 = 0.52 with all variables except terrain) ( figure 4(a)).
The results regarding the relative influence of the LMM, which includes all environmental variables, suggest that terrain is the most important variable in terms of explaining the spatial variation in greening trends (relative importance = 42%). Altitude has been found to be the most important indicator within terrain variables in driving the spatial variation in greening trends (relative importance = 14%) ( figure 4(b)). PCb1, which was mainly correlated with temperature ( figure S2(c)), showed the highest contribution in the baseline climate to explaining the spatial variance of greening trends (relative importance = 13%). Overall, the contribution of soil to explaining the variance of greening was comparable to that of other variables (relative importance = 28%). According to our results, the ECEC showed the highest relative importance value (8%) for the greening trend among the soil attributes ( figure 4(b)).

Interconnections between greening trend and its drivers
PLS-PM suggests that climate change is directly and negatively correlated with the spatial variability of greening trends over the whole TP (correlation coefficients = −0.40) ( figure 5(a)). However, the correlation between climate change and greening could be different among ecosystems, which is 0.57 in the woodland and −0.11 in the grassland (figures 5(b) and (c)). According to the loading score of each indicator to the corresponding variable, PCc1, which was mainly correlated with temperature ( figure S2(c)), was a more powerful indicator of climate change for the model over the whole TP and woodland (loading score = 0.76 and 0.89, respectively). However, PCc2, which was mainly correlated with evaporation ( figure  S2(c)), was a more powerful indicator of climate change for the grassland (loading score = −0.96) (table S3). Climate change also indirectly affects the greening trend through its impact on soils (correlation coefficients = 0.19) ( figure 5(a)). Baseline climate, which was mainly derived by temperature on PCb2 (loading scores = −0.77, figure S2(f), table S3), was shown to have a higher direct contribution than climate change in terms of explaining the variance of greening (correlation coefficient = −0.72) ( figure 5(a)).
Terrain, which is mainly reflected by the TRI and TWI (loading scores in the path analysis are −0.90 and 0.90, respectively; table S4), had direct and negative control on the variability of greenness (correlation coefficient = −0.24) ( figure 5(d)). This finding could be amplified over woodland ecosystems, with a direct correlation coefficient of −0.65 between terrain and greening ( figure 5(e)). The terrain also had an indirect relationship with the greening rate via all hypothesized pathways, including those involving climate change, baseline climate and soil (figure 5). The direct impact of climate change on greening could be strengthened over woodland and grassland when the terrain attributes were included in the path analysis (figures 5(e) and (f)).

Discussion
As expected, climate change has a significant direct effect on the spatial variability of greening trends. The finding of a high correlation between temperature and vegetation dynamics (table S2 and figure S2(f)) was in line with previous research showing that the vegetation phenology at high latitudes is primarily determined by temperature (Shen et al 2016). However, the temperature increase may have variable effects on greening among ecosystems. In this study, woodland ecosystem experienced a relatively rapid greening trend during the warming period over 1982-2015. The increasing temperature, which primarily contributed to the indicator of PCc1 in the climate change variable, promoted the greening over the woodland ecosystem ( figure 5(b), table S3, figure S2(f)), indicating that warmer temperature may lead to an extended growing season for temperate forests due to earlier greenup and delayed senescence (Hwang et al 2014). However, increasing temperature will have a negative impact on the greening over the grassland ecosystem. Warming over the TP induced an increased preseason precipitation and corresponding reduced precipitation during the growing season. The deficient sunshine intensity and duration of preseason, and the drought stress induced by dramatic increases in temperature and inadequate precipitation during the growing season, could counteract early spring carbon assimilation (Angert et al 2005) and increase vegetation mortality (Adams et al 2009). Accelerated leaf senescence associated with drought stress along with warming has been studied specifically in grasslands (Rivero et al 2007), for which chemical mechanisms are largely known.
Moreover, climate change not only favors the rapid growth of plants and controls the structure and function of ecosystems but also accelerates permafrost thawing throughout the TP (Chen et al 2013). In past decades, the mean annual permafrost temperatures at a depth of 6.0 m increased by 0.12 • C-0.67 • C (Wu and Zhang 2008). Such permafrost degradation caused rapid soil carbon loss through the efflux of CH 4 and CO 2 . In turn, the decomposition of soil organic carbon over these areas has accelerated with climate warming (Chen et al 2013). On the other hand, a changing climate has the potential to alter the composition of plant and soil communities and the interactions between them. The interactions among plants and their associated abiotic soils, which are termed plant-soil feedbacks, can lead to complex feedbacks that regulate plant community dynamics and ecosystem processes (Pugnaire et al 2019).
There is growing evidence that the impact of climate change on the trend of vegetation greenness can be amplified by terrain attributes (Liu et al 2013, Radula et al 2018. The influence of elevation-dependent warming on the vegetation dynamics in mountainous environments has also been documented, revealing that warming is more rapid at higher elevations (Mountain Research Initiative EDW Working Group 2015, Xu et al 2020). For instance, warming promotes seedling establishment and vegetation infilling (Du et al 2018) and induces the upward migration of the treeline on the TP. The terrain also causes localized changes in soil moisture and soil temperature. The steepness, shape and slope of an area affect the flow of rainwater into or off the soil as well as the amount of heat absorbed by the soil. Moreover, the structure and function of ecosystems vary significantly with climate driven by fine-scale variations in topography. It is usually assumed that the spatial variability of the greening trend over the mid-and high-latitude forest phenology is primarily determined by temperature, while seasonal rainfall is a more important factor for the interannual vegetation greenness dynamics in arid and semiarid grassland regions (figure 2) (Hwang et al 2014).
Baseline climate also has an impact on the greenness trend through its relationship with climate change, soil, and vegetation type (figure 5). At higher latitudes, plants slow or postpone developmental processes to mitigate the high risk of freezing injury at very low temperatures (Shen et al 2016). Frozen soil water under low temperatures could limit water absorption by alpine vegetation roots. The dominant enhancing effect of the increasing minimum temperature was likely caused by the reduction in low-temperature constraints on alpine plants, which promoted vegetation growth. In this study, the grassland with a higher average temperature (i.e. meadow) showed a higher greening rate than that in a colder environment (i.e. steppe) (figures 2 and 4). Moreover, the combined effects of changes in snow depth and cloud cover in response to the increase in CO 2 concentrations result in greater surface heat storage in high-mountain regions and therefore amplify the warming rate with elevation on the TP (Palazzi et al 2017). The implications for the combination of cloud radiation and snow albedo feedbacks and the increases in water vapor content in the atmosphere can thus affect the greening trend (Mountain Research Initiative EDW Working Group 2015, Yan et al 2016).

Conclusion
This study used multiple data sources and approaches to evaluate the greening trend of the TP vegetation and revealed the drivers of the variation in this greening. The findings of this study have several important implications. First, our results suggest an overall greening trend over the TP in recent decades. Second, the variability of greening in recent decades has been driven not only by climate change but also by terrain, baseline climate and soil. Finally, the interconnections of climate change, terrain, baseline climate and soil need to be considered in understanding and predicting the vegetation dynamics over the TP. The methods and logic developed here could aid understanding of the influence of climate change on terrestrial vegetation and toward improving the simulation of future vegetation dynamics.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).