Elevation-dependent response of vegetation dynamics to climate change in a cold mountainous region

Our study period is from 1982 to 2015. The data we used include: (i) the long-term AVHRR (Advanced Very High Resolution Radiometer) GIMMS3 g (Global Inventory, Monitoring, and Modelling Studies) NDVI product (GIMMS3 g thereafter; Tucker et al 2005); (ii) the MODIS NDVI product (MOD13C2, Collection 6, Didan et al 2015); (iii) the observation-based gridded climate data (mean air temperature (T_mean ), precipitation (P) and downward longwave radiation (R_ld )) (He et al 2011); (iv) the snow cover extent (SCE) product of Hori et al (2017); (v) the GLEAM (Global Land-surface Evaporation Model) evapotranspiration product (Miralles et al 2011); (vi) the SRTM 1 km DEM data (http://srtm.csi.cgiar.org/); and (vii) the hydrologic data (runoff and pan evaporation (E₀ )). Table 1 summarizes the period of observations, spatiotemporal resolution, and data source for each data set. More detailed descriptions are provided in the supplementary material.

2.2.1. Evaluating spatiotemporal consistency between GIMMS3 g and MODIS NDVI datasets

2.2.2. Impact of climate change on montane vegetation dynamics

We examined two important vegetation dynamic metrics: (i) the maximum growing-season greenness (represented by the monthly maximum NDVI); and (ii) the beginning of growing season (BGS). Due to the sparse ground observations, assessment of BGS changes on the HYRB is heavily dependent on remote sensing data. The cloud- and snow-removed biweekly GIMMS3 g NDVI data were used to retrieve BGS. We used a six-degree polynomial function (equation (1)) to reconstruct daily NDVI time series (Piao et al 2006). Then, following White et al (2009), we employed a dynamic local threshold NDVI ratio of 0.5 based on the annual minimum and maximum values to determine the beginning of growing season (BGS) for each pixel (equation (2)).

$$\begin{equation}V{I_t} = {a_1} + {a_2}t + {a_3}{t^2} + {a_4}{t^3} + \ldots + {a_n}{t^n}\end{equation} \tag{ 1 }$$

$$\begin{equation}{\text{NDV}}{{\text{I}}_{ratio\left( t \right)}} = \frac{{NDV{I_t} - NDV{I_{min}}}}{{NDV{I_{max}} - NDV{I_{min}}}}\end{equation} \tag{ 2 }$$

where t is the day of year (DOY), parameters a₁, a₂, ..., a_n are the fitted coefficients (using least squares). NDVI_ratio(t) is the estimated ratio of NDVI change which varies from 0 to 1.0; NDVI_min and NDVI_max are the annual minimum and maximum NDVI.

Although the biweekly composited GIMMS3 g NDVI product is widely used to detect vegetation phenology at both regional and global scales (e.g. White et al 2009, Julien and Sobrino 2009, Jeong et al 2011), there is still a debate regarding the sufficient temporal resolution for retrieving vegetation phenology. For example, Zhang et al (2009) evaluated the impacts of compositing period on vegetation phenology detection based on the 16 d MODIS Enhanced Vegetation Index (EVI) data. They found that time series with a temporal resolution of 6–16 d were sufficient for accurate detection of several phenology events. Pouliot et al (2011) compared satellite-derived spring green-up onset date with in situ observations for deciduous forests and found that compositing periods ranging from 7 to 11 d provided the best estimates for leaf out. Coarser or finer compositing periods resulted in a greater error. Comparing the estimated BGS with in situ observations during 2003–2011, we found that the mean estimated BGS was DOY 142 ± 8 and the mean BGS estimated from eight agro-meteorological stations (dominated by alpine grassland) was DOY 136 ± 10. Due to the biweekly temporal resolution of the NDVI, such errors are expected. We note that the temporal sensitivity for accurately detecting vegetation phenology from satellite-based vegetation data is still poorly understood. We strongly suggest that much effort be directed towards comparing satellite-based vegetation phenology from multiple temporal resolutions with field observations in the future.

Some studies have found that vegetation responds differently to daytime (T_max ) and nighttime (T_min ) warming in the temperate steppe of northern China (Shen et al 2016, Xu et al 2018), eastern Europe, southern Siberia, western America (Xia et al 2014, Tan et al 2015), and eastern Canada (Rossi and Isabel 2017). Our results indicate that the inter-annual variations in growing season NDVI were more strongly associated with changes in T_max than those in T_min in the HYRB (figure S3). We estimated the sensitivity of NDVI to T_max (or T_mean ) based on multiple linear regressions with NDVI as the dependent variable and T_max (or T_mean ), P and SCE as independent variables. We found that the mean sensitivity of NDVI to T_mean was slightly higher compared to that of T_max (figure S4). Therefore, we did not consider the non-uniform (daytime/nighttime) warming because this effect is smaller than changes in T_mean . We used the Mann-Kendall (MK) test (two-tailed test) to determine whether pixel-level and basin-scale temporal trends in vegetation greenness (monthly maximum NDVI), BGS and climate factors were statistically significant (p < 0.05). Positive and negative trends indicate greening and browning in NDVI (or delaying and advancing in BGS), respectively. The magnitude of the trend was determined using linear regression. We evaluated the strength of relationships between the two vegetation metrics (greenness and BGS) and climate factors using partial correlation analysis (Liu et al 2016a). This method can explore correlations between NDVI/BGS and single climate variable while removing effects of the two remaining climate variables in this study. We defined spring as March through May, summer as June through August, autumn as September through November, and winter as December through the next February.

2.2.3. Pattern of air temperature change with elevation and response of vegetation dynamics

2.2.4. Hydrological implications of vegetation dynamic change

As vegetation dynamics are strongly coupled with the latent heat flux (Schlesinger and Jasechko 2014), temporal trends in annual actual evapotranspiration (ET_a ) will be consistent with those in annual integrated NDVI and BGS. That is, increases in NDVI and/or earlier BGS should lead to increases in ET_a , and vice versa. As reliable estimation of long-term ET_a remains a challenging problem in hydrologic research (Lettenmaier et al 2015, Purdy et al 2018), we estimate ET_a using multiple methods to reduce uncertainties. Because transpiration and canopy interception are positively related to leaf area index (LAI) and soil evaporation is negatively related to LAI, their contribution to ET_a changes should be different. To evaluate this, we investigated temporal trends in ET_a and its components. In addition to ET_a calculated from water-balance and by GLEAM, we also estimated sequential 5-year averages of ET_a using the generalized Budyko equation (equation (3)) (Choudhury 1999). This equation assumes that changes in water storage can be ignored at long-term annual scale (Roderick and Farquhar 2011).

$$\begin{equation}E{T_a} = \frac{{P{E_0}}}{{{{\left( {{P^n} + {E_0}^n} \right)}^n}}}\end{equation} \tag{ 3 }$$

where n is a catchment-specific parameter, which modifies the partitioning of P between ET_a and runoff. We used a recommended default value of n = 1.8 (Choudhury 1999). E₀ is catchment-averaged pan evaporation.

During the 34-year period (1982–2015), trends in NDVI were positive for 84%, 58% and 79% of the HYRB pixels in spring, summer and autumn, respectively; with 40%, 16% and 29% being statistically significant (p < 0.05; see figures 1(a)–(c)). The relative contribution of pixels with negative trends to total pixels was 16%, 42%, and 21% in spring, summer and autumn, respectively. The corresponding seasonal percentages with significant negative trends was 2%, 8% and 19%, respectively. Regional inter-annual NDVI in spring showed significant increasing trends at a rate of 0.006 per year (not statistically significant in summer and autumn) (figures 2(a)–(c)). Pinzon and Tucker (2014) analyzed uncertainties or errors of the GIMMS3 g NDVI data and found that errors resulted from seasonal variations, which contributed the largest proportion of the total error, depend upon NDVI values. That is, large errors correspond to high NDVI values. They estimated a total error of ±0.005 NDVI units. Due to relatively small NDVI values in the HYRB (figures 2(a)–(c)), minor uncertainties are expected. However, one should use the NDVI data with caution for climate change research.

Zoom In Zoom Out Reset image size

We examined linkages between vegetation growth and climate factors and found that NDVI was positively correlated with T_mean during the three non-winter seasons (figures 1(e)–(g)). There were positive partial correlations in 90%, 75% and 85% of the HYRB in spring, summer and autumn, respectively, with 54%, 25% and 28% of correlations being statistically significant. In spring, 61% of the pixels had positive partial correlations between vegetation greenness and P, while in summer and autumn, the partial correlations were negative in 65% and 59% of the pixels, respectively. Partial correlations between NDVI and SCE were dominantly negative in spring and ambiguous in summer and autumn. The results suggest that growing-season greening was primarily stimulated by climate warming across the HYRB. This finding is consistent with previous studies of the entire TP (Zhang et al 2013), the western U.S. (Cayan et al 2001), and central Europe (Kudernatsch et al 2008). We note that evidence from field warming experiments demonstrates that rising temperatures relax low-temperature limitations and enhance alpine grassland growth (Kudernatsch et al 2008, Wang et al 2011, 2012), which are consistent with our own.

We observed early BGS in the southeastern HYRB, whereas later BGS was mainly in the northwestern HYRB (figure S5(a)), which followed the spring gradient of temperature and precipitation. The growing season started in the middle of May in most areas while other areas exhibited earlier (beginning of April) or later (end of May) growing season (figure S5(b)). BGS trend analysis revealed that 53% had positive trends, of which 6% were statistically significant (figure 1(d)). Regionally, we found no significant trends in BGS over the entire study period. However, we identified two distinctly different trends between 1982–1998 and 1998–2015. During 1982–1998, BGS significantly advanced by 4 d/decade which is similar to the mean advancing rate of growing season in North America (4 d/decade) and Eurasia (3 d/decade) (Zhou et al 2001). In contrast, BGS had a small delay of 1 d/decade during 1982–1998 (figure 2(d)). Inter-annual variations in spring temperature show that the HYRB exhibited continuous warming between 1982 and 2015 (figure S6(a)) which allows us to suggest a non-linear response of spring phenology to temperature change. Climate warming has been widely recognized as advancing spring vegetation phenology, while both in situ observations and satellite records reveal that warming in autumn and winter decreases the fulfilment of chilling requirements and result in later spring green-up onset (Zhang et al 2007, Guo et al 2015). This was supported by weaker correlations between BGS and spring temperature during 1998–2015 than the period of 1982–1998 (figure S7).

We found negative partial correlations between BGS and spring T_mean at 71% of pixels; 9% of correlations were statistically significant (figure 1(h)). Partial correlations between BGS and spring P were generally positive. Our results indicate that both spring warming and wetting could promote advancing of spring phenology. This is mainly because higher temperature and precipitation in spring enhances photosynthesis and alleviates water stress. In addition to temperature and precipitation, BGS may also have been affected by reduced winter snow cover. Partial correlations between BGS and SCE were positive for 80% of the area, with 6% of the correlations being statistically significant. This finding highlights the important role of snow-related processes in alpine grassland spring vegetation activity.

We observed enhanced warming at elevations below 4300 m, but warming rates greatly declined above this elevation (figure 3(a)). Examining patterns of possible causes with elevation, we found spring P increased with elevation, while summer and autumn P shifted from negative trends to positive trends at about 4000 m and increased at higher elevations (figure 3(b)). SCE declined at lower elevations but increased above 4650 m (figure 3(c)). Generally, R_ld increased at each elevation which could be attributed to increased water vapor (figure 3(d)). The pattern of ΔR_ld with elevation was very similar to that of ΔT_mean . A correlation analysis showed that ΔT_mean was statistically significant and positively related to ΔR_ld but negatively related to ΔP and ΔSCE. We used a stepwise regression model to investigate their combined effects and to identify the primary driver of ΔT_mean . Inclusion of the three independent variables (i.e. ΔR_ld , ΔP and ΔSCE) accounted for 83%, 91%, 83% and 90% variance of ΔT_mean during spring, summer, autumn and the growing season, respectively (table 2). Among the three variables, ΔR_ld played a leading role followed by ΔSCE and ΔP, as the ΔR_ld univariate model explained most of the variance of ΔT_mean (table 2). It should be noted that the gridded precipitation data were developed using more than one dataset as background (including station observations, TRMM satellite precipitation, and GLDAS precipitation data) which would introduce uncertainties and biases in the data time series (Beck et al 2020).

Zoom In Zoom Out Reset image size

Most studies reported enhanced warming at higher elevations (table S1; Nogués-Bravo et al 2007, MRI 2015), however, we observed EDW only below 4300 m. Our results are generally consistent with Gao et al (2018) who analyzed changes in temperature with elevation over the TP using observed, reanalyzed and WRF (Weather Research and Forecasting) model simulated data. They found a clear signal of EDW below 5000 m, but reported that neither ERA-Interim reanalysis nor the WRF model indicated EDW above this elevation. They explored the physical process based on WRF output and found that snowpack showed an increasing melt under 5000 m resulting in lower albedo. However, increased snowfall in a warmer climate partially compensates for increased snowmelt and thus leads to a little change in albedo above 5000 m. Gao et al (2018) further concluded that EDW above 5000 m elevation over the TP is unlikely to occur in future decades based on their analysis of WRF climate projections. Our findings support their conclusion although with a somewhat different elevation threshold. As there is no strong consensus on EDW even within the same region, future studies should aim to accurately detect the signal and mechanisms through combined analysis of observations, model-simulated and satellite-retrieved datasets.

Examining changes in the two vegetation metrics with elevation, we found that spring and autumn ΔNDVI corresponded well with the concurrent pattern of ΔT_mean (figure 3(g)), suggesting a strong control of the elevation gradient of ΔT_mean on vegetation greening in two seasons. Summer NDVI exhibited browning trends below 4000 m and increased greening trends above this elevation. Because this pattern was consistent with the elevation gradient of summer ΔP, it can be reasonably inferred that ΔP was likely an important regulator of summer vegetation greenness across the elevation range. This result suggests a high sensitivity of vegetation growth to summer precipitation in semi-arid environments (Poulter et al 2013). Carlson et al (2017) investigated changes in peak NDVI (NDVI_max) with elevation in the Ecrins National Park (within the French Alps) and found that the highest proportional increases in NDVI_max occurred in the rocky habitats at high elevation. They attributed possible causes to glacier retreat, decrease in snow cover duration, and increased atmospheric nitrogen deposition at high elevations. Noting that they calculate relative values of NDVI_max trends rather than absolute values which is likely to change the elevation-dependence patterns, we analyzed relative values of NDVI trends with elevation as well. The result shows that the general patterns of the relative change in NDVI with elevation were consistent with absolute values (figure S8). Therefore, the elevation-dependent response of greening is closely related to patterns of local climate change. ΔBGS showed an evident elevation gradient with greater rates of advance occurring at lower elevations (figure 3(h)), which are positively related to ΔSCE and spring ΔT_mean .

Other studies have also observed elevation dependence of changes in greenness and BGS (e.g. Piao et al 2011, Li et al 2015), but most attributed the cause to EDW. Our findings reveal that seasonal changes in ΔT_mean, P, and SCE with elevation played different roles in the vegetation dynamic response. This is the first time that the impacts of P and SCE changes with elevation on vegetation activity have been evaluated in the region. The aspect is another important topographic factor that promotes variations in energy and water in complex terrain (Fan et al 2019). We analyzed changes in growing season NDVI and BGS with aspect and found that vegetation greening was slightly higher on the south, southwest and west slopes, while there were no obvious differences in spring and autumn (figures S9(a)–(c)). BGS advanced on the southeast and northwest slopes and delayed at other slopes (figure S9(d)). The elevation dependency of the vegetation dynamic changes has important implications for energy and mass balance of the high-mountain regions and associated hydrology and also for species that live in narrow habitats (MRI 2015). We assume that climate change is the main controller of elevation dependency of vegetation dynamics. Our analysis may overstate the influence of climate on vegetation growth if non-climatic factors that vary with elevation help explain spatial pattern of vegetation dynamics.

Another problem that mountain regions face globally is the lack of observations at high elevations. Taking the elevational distribution of precipitation observations as an example, the GPCC (Global Precipitation Climatology Centre) archive shows that precipitation observations for elevations above 2500 m are under-represented (above this elevation only 0.2% of the 48 767 GPCC precipitation stations are located (excluding Greenland and Antarctica)). In terms of network density per continent, Asia has the lowest density (0.02 stations per 100 km² or 5633 km² per station) and the value above 2500 m is even lower (0.004 stations per 100 km² or 26 692 km² per station) (Viviroli et al 2011). Environmental monitoring stations (at which both hydrometeorological and ecological observations are made) should be installed at higher-elevation and ecologically sensitive locations which exhibit rapid and detectable changes, rather than being located (as at present) in lower elevations, notwithstanding a range of pragmatic challenges associated with higher elevations, experiencing more extreme conditions and affordability.

Annual HYRB-averaged ET_a calculated from water balance (P − runoff) varied between 278.3 and 475.1 mm/year with a mean value of 371.0 mm/year during 1982–2015. The Budyko-estimated ET_a ranged from 405.7 to 452.8 mm/year with a mean value of 431.1 mm/year. The GLEAM ET_a rates exhibited smaller temporal variability during 1982–2015 with a mean value of 391.9 mm/year (figures 4(a)–(c)). The considerable discrepancy between the measured and GLEAM-estimated ET_a series is possibly due to its great underestimation of potential evaporation in comparison with pan evaporation (figure S10). With differences in magnitude among ET_a time series neglected, all three showed significant positive trends during the past 34 years, suggesting that vegetation greening and/or earlier growing season are imposing positive effects on evapotranspiration. This agrees with both global and regional analyses that vegetation greening and growing-season advancing are important drivers of long-term ET_a variations (Ukkola et al 2016, Zhang et al 2016, Liu et al 2016b). For example, Zeng et al (2018) quantified the potential effect of the satellite-observed global greening on the terrestrial hydrologic circle by using a coupled land-climate model and found that 8% increase in LAI resulted in 55% ± 25% increase in observed ET_a during 1982–2011.

Zoom In Zoom Out Reset image size

We also examined temporal variations of the ET_a components from GLEAM products, the results of which showed that annual transpiration and canopy interception increased by 0.83 and 0.04 mm/year whereas soil evaporation increased by 0.21 mm/year, which further supports our assumption (figures 4(d)–(f)). Because reliable estimation of long-term ET_a and its components remains a challenging problem in hydrologic research (Lettenmaier et al 2015, Purdy et al 2018), comprehensive intercomparison within and between output of land surface models and satellite-derived ET_a products provides robust support for our assumption which could motivate future research. Analyzing runoff variation in the study domain, we found that annual runoff from three hydrological stations in the upper, middle and lower HYRB (see figure S2(c) for specific location of each station) all showed decreasing trends with different magnitudes (figures 4(g)–(h)). Zhang et al (2018) recently reported that, in arid and semi-arid China (which includes the HYRB), increased LAI has resulted in an 8.5% runoff reduction in 1990 s and 11.7% reduction in 2000 s compared to the 1980 s baseline. Our findings along with the abovementioned study provide evidence at catchment/regional scale for potential drying environment and reducing water resource.

Our analysis is a first step in establishing the variation of ET_a and runoff as related to vegetation change in the study domain. Hydrologic models assume that vegetation dynamics (vegetation greenness, types, phenology and density) are static, whereas both our findings and those of other research show that elevation dependence of vegetation greenness and phenology changes is likely (Piao et al 2011, Trujillo et al 2012, Li et al 2015). Goulden and Bales (2014) found that both eddy covariance measurements and remotely sensed annual ET_a (calculated from MODIS NDVI data) in California's upper Kings River basin peaks at mid-elevation (about 2000 m) and declines at lower elevation due to water availability limitation and upper elevation due to temperature limitations. They explored sensitivity of river flow to vegetation expansion and found that continuous warming is expected to expand ET_a upslope due to vegetation redistribution and to increase basin ET_a and decrease runoff. It is also important to elucidate the runoff sensitivity to vegetation change at different elevation bands to fully evaluate the impact of climate change on mountain hydrology of the HYRB. However, it is not possible at present due to the sparsely distributed weather and hydrologic observations in the harsh environment.