Variations of Vegetation Phenology Extracted from Remote Sensing Data over the Tibetan Plateau Hinterland during 2000–2014

How vegetation phenology responds to climate change is a key to the understanding of the mechanisms driving historic and future changes in regional terrestrial ecosystem productivity. Based on the 250-m and 8-day moderate resolution imaging spectroradiometer (MODIS) normalized difference vegetation index (NDVI) data for 2000–2014 in the Three-River Source Region (TRSR) of Qinghai Province, China, i.e., the hinterland of the Tibetan Plateau, we extracted relevant vegetation phenological information (e.g., start, end, and length of growing season) and analyzed the changes in the TRSR vegetation in response to climate change. The results reveal that, under the increasingly warm and humid climate, the start of vegetation growing season (SOS) advanced 1.03 day yr−1 while the end of vegetation growing season (EOS) exhibited no significant changes, which led to extended growing season length. It is found that the SOS was greatly affected by the preceding winter precipitation, with progressively enhanced precipitation facilitating an earlier SOS. Moreover, as the variations of SOS and its trend depended strongly on topography, we estimated the elevation break-points for SOS. The lower the elevations were, the earlier the SOS started. In the areas below 3095-m elevation, the SOS delay changed rapidly with increasing elevation; whereas above that, the SOS changes were relatively minor. The SOS trend had three elevation break-points at 2660, 3880, and 5240 m.


Introduction
Phenology generally concerns the seasonality of recurring biological events and the mechanisms controlling them, especially with regard to meteorological phenomena (Justice et al., 1985). In response to climate changes, phenology affects the carbon flux and storage in the terrestrial ecosystems through changes in the length of the growing season and feedback to climate changes via vegetation growth (Jeong et al., 2011;Cao et al., 2015;Fan et al., 2015;Wang C. Z. et al., 2015;Ding et al., 2016a). As one of the key sources of uncertainty in terrestrial biosphere models and even the marked uncertainty in projections of the future climate (Wu et al., 2017), the understanding of the seasonality and its mechanisms is currently a hot topic, especially against the background of global climate changes.
As the earth's third pole, the Tibetan Plateau (TP) is one of the most important regions in the world because of its unique location, high elevation (topography), and overall climate. Under global warming, the vegetation phenology of the TP has been investigated for its responses to experimental warming (Fu and Zhong, 2016;Zhu et al., 2016), spatiotemporal patterns, their climatic dependencies (Piao et al., 2011;Shen et al., 2011;Jin et al., 2013;Zhang et al., 2013;Shen et al., 2014;Wang H. S. et al., 2015;Ding et al., 2016a, b), and the methods for detecting the key phenological dates from remotely sensed time-series data (Cao et al., 2015).
Nowadays, there are various explanations for the trends of vegetation phenology and their mechanisms. Across most of the plateau, temperature is the main factor that affects the start of vegetation growing season (SOS), whereas precipitation and solar radiation are known to also play key roles in the SOS in a few areas . Specifically, in the relatively humid area of the plateau, the SOS is more sensitive to the preseason temperature than to the pre-season precipitation; whereas in the arid area, the SOS is more sensitive to the pre-season precipitation than to the pre-season temperature . However, due to the large extent of this region and its complex topography, climate variations, and diverse vegetation types, the mechanisms that control the TP vegetation phenology also vary with the elevation gradient. Thus, the regional average phenological change does not represent the specific local SOS very well, and the mechanisms that control the regional phenological changes at more localized scales are poorly understood.
As the source of the Yangtze River, Yellow River, and Lancang River, the "Three-River Source Region" (TRSR) is the hinterland of the TP where ecosystem changes have been attracting considerable attention (Liu et al., 2008). Due to the rapid global warming, intense human activities, and its special location, the ecosystem in TRSR is more fragile than that of the other parts of the TP. Over the past few decades, the grassland in TRSR has continued to deteriorate. However, after 2005, with the implementation of a series of ecological engineering projects, this trend of degradation has been alleviated to some extent (Shao et al., 2017). However, how the degradation impacts the phenology of alpine ecosystems and whether the projects have had an influence on the phenology are still open questions.
To address some of these issues, this paper presents an investigation using moderate resolution imaging spectroradiometer (MODIS) reflectance data (MOD09Q1) at a fine spatial resolution for a 15-yr period (2000-2014) to explore the variation features of vegetation phenology in TRSR under climate change, through using phenological information extracted from the remote sensing data. The ultimate goal of this study is to better understand the mechanisms controlling grassland productivity, to optimize grazing, and to facilitate the protective utilization of grassland resources over the TP.

Study area
The TRSR is located in the hinterland of the TP over 31°39′-36°12′N, 89°45′-102°23′E (Fig. 1a). This region has a typical plateau continental climate. Grassland accounts for 69% of the total area, followed by desert, which accounts for 16% ( Fig. 1b; based on the land cover data provided by the National Earth System Science Data Center, China). According to the interpolated meteorological grid data of precipitation and air temperature for 2000-2014 (Wang J. B. et al., 2017), the annual total precipitation (ATP) averages 470 ± 134 mm, and gradually decreases from southeast to northwest (Fig. 1c). The annual mean air temperature is −3 ± 3.7°C, and is significantly higher in the eastern and southeastern areas than in the western and northwestern areas, and significantly higher in the valleys than elsewhere (Fig. 1d).

Remote sensing data
The MOD09Q1 data were used to calculate the normalized difference vegetation index (NDVI) for 2000-2014 and conduct an phenological analysis. The MODIS NDVI dataset was obtained from NASA's Earth Observing System (http://reverb.echo.nasa.gov/). The spatial resolution of MOD09Q1 is 250 m, corresponding to a maximum synthetic reflectance for 8 days. First, we used the MODIS Reprojection Tool (MRT) of NASA for data mosaic, re-projection, and band information extraction. Then, we clipped the data and converted the data format in ArcGIS 10.2 from the Environmental Systems Research Institute, Redlands, California. Based on the extracted band information, we calculated NDVI = (NIR − R)/(NIR + R), where NIR and R are the near-infrared and red spectral reflectances, respectively. The smoothing and noise reduction of the NDVI data were conducted with the S-G filtering method of TIME-SAT3.2 (Jönsson and Eklundh, 2004), which was used to reconstruct the NDVI time series. We also calculated the annual total NDVI (ATN) and the summer mean NDVI for analysis.

Meteorological data
The air temperature and precipitation data used in this study were from 753 meteorological stations across China from the China Meteorological Science Data Sharing Service Network, including 18 stations in the TRSR and 345 stations nearby. We used the ANUSPLINE interpolation software program (Hutchinson, 1991(Hutchinson, , 1998Tan et al., 2016) for the spatial interpolation, and thereby obtained temperature and precipitation data with a spa- The studies have shown that temperature and precipitation interpolations obtained by using ANUS-PLINE can explain 90% and 77% of the spatial variance, respectively, i.e., such data can represent their spatial variations very well . Based on the 8-day interval interpolated meteorological data, we calculated the monthly, seasonal, and annual temperature and precipitation data for further analysis.

Phenological information extraction
We used long-term data reconstructed with S-G filtering and TIMESAT3.2 (Jönsson and Eklundh, 2004). We extracted 11 types of phenological parameters, denoted by a-k (Fig. 2a). Specifically, a and b denote the SOS and end of vegetation growing season (EOS), respectively; c is the length of vegetation growing season (LOS), which is the time span between the SOS and EOS; d is the base value; e is the middle time of the growing season; f is the maximum vegetation index during the growing season; g is the amplitude of the growing season; h is the NDVI integral during the growing season; h + i denotes the upper integral value; and j and k are the NDVI left and right derivatives, respectively. The start and end dates of the growing season were based on the NDVI amplitude; that is, the SOS was defined as the date when the increasing NDVI reached a certain percentage of the amplitude of that year, and the EOS as the date when the NDVI decreased to a certain percentage of the amplitude of that year (Fig. 2b). Applying this method to a large area by using multi-year data avoids the disadvantage of the threshold method, which needs to determine a reasonable threshold given the differences in NDVI among vegetation types across a large area. In addition, the NDVI variations in different years will not affect one another, so it is possible to obtain appropriate results. On the basis of a large number of pixel-based analyses, the thresholds were determined as 20% for SOS and 30% for EOS.

Analysis methods
The methods of linear regression and Pearson coefficients were used to analyze the influences of climate on the trend of SOS. In order to explore the influences of the climatic factors on the SOS change in a time series, the SOS data for the main alpine grassland types in the headwater areas of the three rivers were extracted. At the same time, the temperature and precipitation data for every 8-day interval from July of the previous year to July  of the current season, which were consistent with the spatial resolution of phenological information, were selected for correlation analysis with the SOS of the main alpine grassland types in the three different headwater areas. We averaged the fitted correlation coefficients for every five values.
To analyze the impact of topography on the relationship between climate variables and NDVI, mean values were calculated at elevation intervals of 50 m for ATP, AMT, annual summer mean NDVI (AMN), ATN, and SOS, as well as their trends from 2000 to 2014. AMN was averaged in summer (June-August) of 2000-2014. Moreover, to analyze the influences of climatic factors such as temperature and precipitation on the interannual changes of SOS in the TRSR, the SOS's trend was averaged within precipitation intervals of 20 mm, and its changes along the precipitation gradient in the study area were then probed. Based on the overall regional mean temperature ± 0.5 standard deviation (SD), i.e., an interval of 1 SD, the study area was classified into three temperature zones, i.e., colder/lower temperature (−6.8 to −4.2°C), prevailing/median temperature (−4.2 to −1.6°C), and warmer/higher temperature (−1.6 to 1.0°C) zones, to disclose how the trends of the SOS change with temperature and precipitation variations. Figure 3 shows the spatial distributions of the start and end dates, expressed as day of year (DOY), and the length of growing season averaged for the multi-year data from 2000 to 2014. The overall SOS was delayed from southeast to northwest (Fig. 3a). The growing season begins between April 30 and May 20 (120-140 DOY) in the eastern and central regions; while it is between May 20 and June 20 (140-170 DOY) in the northwest region. However, more than 84% of the region began to grow in the period of 120-160 DOY with a regional mean date of 147 ± 10 DOY (15-30 May), the mean end date was 285 ± 6.5 DOY (Fig. 3b), and the average LOS was 140 ± 14 days (Fig. 3c). In the trend analysis, no significant changes were found in the EOS. The changes of LOS were mainly affected by SOS, so we mainly focused on the changes of SOS in the following results.

Temporal and spatial patterns
This headwater area experienced warming-wetting climate changes during 2000-2014 (Fig. 4). The ATP showed a significant increasing trend at a rate of 73.21 mm decade −1 (R 2 = 0.383, p < 0.05), while the annual mean air temperature significantly increased at a rate of 0.86°C decade −1 (R 2 = 0.500, p < 0.05). Responding to the climate changes, the SOS advanced significantly at a mean rate of 1.03 day yr −1 (R 2 = 0.335, p < 0.05) in the TRSR during 2000-2014. The SOS of alpine meadow advanced 0.63 day yr −1 (R 2 = 0.254, p < 0.05) while the alpine steppe advanced 0.65 day yr −1 (R 2 = 0.526, p < 0.05; Fig. 5). The SOS of the alpine meadow occurred mostly between the end of May and beginning of June (140-155 DOY) during 2000-2007, and advanced significantly to mid-May during 2008-2009. The areas where the SOS advanced accounted for 80% of the region, and this advance was greater in the eastern and southern valley areas (Fig. 6). Figure 7 shows the headwater boundary and main grassland types (grassland type data were provided by   the National Earth System Science Data Center, China) in the TRSR. The main grassland types in the Yellow River and the Yangtze River sources were alpine meadow and alpine steppe, but in the Lancang River source, it was mainly alpine meadow.

Responses in phenology of the main grassland types to climate in TRSR
Analysis of the correlations between the precipitation, temperature, and SOS of the main grassland types revealed that the impacts of precipitation and temperature on the spring vegetation phenology were basically con-sistent in the three different headwater regions (Fig. 8). However, in the time series, the changes in precipitation and temperature affected the SOS differently. The precipitation during the preceding year, particularly in October, was significantly negatively correlated with the SOS, while the precipitation after November of the preceding year gradually transitioned from a negative correlation to a positive correlation. The precipitation approximately 30 days ahead of the SOS was slightly negatively correlated    790

Journal of Meteorological Research
Volume 34 with the SOS. Water required for the SOS in the TRSR is derived mainly from melting snow and ice, and the spring precipitation (i.e., the cumulative precipitation in the previous winter) plays an important role. Spring temperatures were negatively correlated with the SOS, especially in January and April (Fig. 8). Figure 9 shows the mean and 95% confidence interval of the SOS trend in each 20-mm precipitation interval in the whole TRSR and in three temperature zones. First, the variance of the SOS trend in the colder area was very different from those in the warmer or median-temperature areas. The whole-region interannual changes in the SOS (Fig. 9a) were consistent with those in the lower-temperature area (Fig. 9b). Under different precipitation amounts, the interannual changes in the SOS generally displayed a bimodal pattern. In the zone receiving 160-240 mm of precipitation, the SOS gradually advanced. In the 240-320-mm zone, the interannual advance of the SOS weakened drastically with the increase in precipitation. Meanwhile, in the 320-520-mm zone, the interannual advance of the SOS continuously increased. In the 520-660-mm zone, the interannual advance showed a single peak; particularly at 660-mm, the interannual advancement reached the maximum. With precipitation exceeding 660 mm, the interannual advance gradually weakened. The interannual advance of the vegetation growing season is more obvious in the areas with more precipitation. However, the appearance of the second peak indicates that in the areas with lower precipitation, precipitation was not the dominant factor affecting the advance of the SOS.

Interannual trends of SOS under different precipitation amounts
The pattern of changes in the interannual advance of the SOS with the changes in precipitation in the lowertemperature zone was generally consistent with the pattern in the entire region, which was significantly different from the changes in the median-and higher-temperature zones, although the changes in the latter two zones were consistent with each other (Fig. 9b). In the 360-660-mm zone, the interannual advance of the onset of the vegetation growing season was relatively stable in the median-and higher-temperature zones. In the 660-800-mm zone, the weakening of the interannual advance of the onset of the vegetation growing season occurred in all three temperature zones. However, this weakening of the advance of the SOS gradually en-  hanced with the increase in temperature.

Changes of NDVI, AMT, ATP, and SOS under different elevation levels
The ATN and the summer mean NDVI each showed a strong dependence on the elevation level, and in both cases, they were more controlled by the dependence of precipitation on the topography than that of air temperature over the TRSR (Fig. 10a). The AMT linearly decreased by a rate of 4.6°C km −1 with increasing elevation. But the precipitation showed a unimodal curve with the maximum around the 4000-m elevation. According to the segmented regression, three break-points were found, and the changes of precipitation could be classified into four elevation ranges: < 3093, 3093-3929, 3929-5001, and > 5001 m. In the region with elevations below 3093 m, ATP increased by 28.14 mm km −1 with rise in elevation, and by 265.85 mm in the 3093-3929-m zone, while it decreased by 197.06 mm km −1 in the 3929-5001-m zone, and showed an increase of 109.76 mm km −1 at elevations over 5001 m (Table 1).
Following the precipitation and temperature changes along the elevation gradient, both the ATN and the AMN showed almost the same patterns of change along the increasing elevation level (Fig. 10a). The changes of NDVI can be explained well by the climatic factors, and the regressed slopes and estimated elevation break-points can be found in Table 1.
The changes of the SOS were more interesting in the area with elevations around 3095 m, which probably represented ridges for the SOS trend (Fig. 10b). In the 2000-3095-m zone, the SOS was delayed 36.64 day km −1 increase in elevation (R 2 = 0.97, p < 0.05), whereas in the 3095-5800-m zone, the SOS advancement showed some fluctuations (Fig. 10b). The regressed slopes and the estimated elevation break-points can be found in Table 1.

Trends of AMT, ATP, and SOS under different elevation levels
The trends of AMT, ATP, and SOS also showed strong dependence on elevation in the study region (Fig.  11a), and we estimated their elevation break-points (Table 1). Although the AMT was still rising in the whole area, the rate of the rise varied depending on the elevation levels, with break-points at 3080, 3920, and 5030 m. Over the area with elevations lower than 3080 m, the climate warming trend was slow, with a speed of around 0.2°C decade −1 , but the speed showed an almost linear increase at the rate of 0.08°C decade −1 km −1 . Meanwhile, in the 3080-3920-m zone, the rate was 0.66°C decade −1 km −1 ; and in the 3920-5030-m zone, the rate was weakened to 0.08°C decade −1 km −1 . Above 5030 m, the rate was 0.34°C decade −1 km −1 . However, the trend of precipitation changes showed a unimodal pattern along the increasing elevation gradient, with the wetting trend increasing as elevation increased, reaching a maximum at elevation 3920 m, and then the wetting speed decreased with the increasing elevation.
The interannual trend of the SOS was averaged at 50-m elevation intervals, and the SOS advanced between the 2000-and 6000-m elevations, although it varied with different elevation levels (Figs. 11a, b)  In the 2660-3880-m zone, the relative delay shifted at a rate of 0.26 day yr −1 km −1 . In the 3880-5240-m zone, the SOS advanced greatly, at a rate of 0.69 day yr −1 km −1 . In the zone above 5240 m, the SOS was delayed greatly, at a rate of 1.72 day yr −1 km −1 (Table 1). In the high elevation/low temperature area, the annual advance trend of vegetation growth was more sensitive to the temperature increase than that in the low elevation/high temperature area, and the annual advancement in days of SOS was larger than that in the low elevation/high temperature area.

Climate influence on the SOS under different elevation levels
Elevation plays an important role in the redistribution of temperature and precipitation, which were both negatively correlated with the SOS, although their effects varied with elevation (Fig. 12). In the zone below 3000 m, precipitation was the dominant factor affecting the SOS, and the negative effects of temperature increased with increasing elevation. In the 3000-4400-m zone, precipitation remained the dominant factor affecting the start of the vegetation growing season. In the 4400-5500-m zone, temperature was the dominant factor. At higher elevations, the effects of both precipitation and temperat-ure gradually weakened.

Climatic factors affecting SOS
Temperature and precipitation are generally regarded as the main factors that affect vegetation phenological changes. The studies have shown that rising spring temperatures play a dominant role in advancing the SOS (Piao et al., 2011;Zhao et al., 2016), and the pre-season temperature and precipitation affect the vegetation growing season differently (Shen et al., 2011). In this study, more precipitation promoted an earlier SOS. Temperature and precipitation 30 days ahead of the SOS exerted the maximum effects. This pattern suggests that in relatively wet areas, the sensitivity of SOS to the pre-season temperature was greater than the sensitivity to the preseason precipitation; whereas in more arid areas, the sensitivity of SOS to the pre-season precipitation was higher than the sensitivity to the pre-season temperature. Overall, both temperature and precipitation affect the SOS (Shen et al., 2011); however, the factor that dominates varies depending on the location and period of time.

The elevation effect
Elevation plays an important role in the redistribution    of temperature and precipitation. In the TRSR, the amplitude of temperature increases intensify with increasing elevation. This study found that the sensitivity of the interannual advance of the SOS to the temperature increase was greater in the high-elevation/low-temperature areas than in the low-elevation/high-temperature areas, and the interannual advance of the growing season was also greater than in the low-elevation/high-temperature areas. This interannual change of the SOS with elevation was consistent with the overall change throughout the TP during 2001-2010 (Qiu et al., 2017), whereas the magnitude differed (Fig. 11b). Figure 11b shows that at a given elevation level, the interannual trend of the SOS in the TRSR advanced at a rate of approximately 1 day more than in the TP. This comparison means that the SOS is more sensitive in the TRSR than in the TP.

Uncertainty
Research based on Global Inventory Modeling and Mapping Studies (GIMMS)-NDVI showed that the SOS in the TP area was advancing at a rate of 0.88 day yr −1 during 1982-1999, whereas the advancing trend has been insignificant since 2000 and has even lagged in the southwestern TP and other regions (Piao et al., 2011;Ding et al., 2016b). However, research based on three different vegetation indexes [GIMMIS-NDVI, SPOT (Systeme Probatoire d'Observation de la Terre) NDVI, and MODIS NDVI], and phenological observations from field stations indicated that the SOS in the TP area is advancing significantly at a mean interannual rate of 1.04 day yr −1 (Zhang et al., 2013) The applications of different remote sensing datasets and different analytical methods to extract vegetation phenological information and analyze changes, in addition to the regional complexity, all lead to uncertainty in the results of various studies. We used the 250-m and 8day MODIS NDVI data to extract the vegetation phenological information for this study. This dataset has relatively high temporal and spatial resolution, and thus can yield more accurate depictions of changes in the vegetation phenology. We used spatially interpolated precipitation and temperature data, which may have affected the interpolations.
Regional plant ecosystems are complex and self-regulating. Therefore, the SOS cannot advance continuously as a result of factors such as climate warming because the plants themselves possess certain adjustment mechanisms for responding to their external environment. The delays of the SOS in 2009 and 2011 relative to past years, and the subsequent fluctuations in the advances, may be related to these self-regulating mechanisms in the plants.
Anthropogenic activities such as grazing will also impact the SOS. The SOS during 2000-2014 in the TRSR clearly advanced since 2005 (Fig. 5). A program called AUGUST 2020 Overall Planning for TRSR Ecological Protection and Construction was implemented in 2005, and a program of ecological protection and ecological construction project was launched. This program has greatly improved the TRSR ecological environment. A decrease in human activity may also have played a role in the advancement of the SOS, although this correlation requires further quantitative analysis.

Conclusions
We used the 250-m and 8-day MODIS NDVI dataset over 2000-2014 to analyze the spatiotemporal changes in the SOS in the TRSR. Due to regional warming and wetting, the SOS advanced significantly, and the growing season became increasingly longer. Winter precipitation, particularly in October-November of the previous year, greatly affected the advance of the SOS, and more winter precipitation promoted earlier growing season starts. The precipitation and temperature in 30 days before SOS had the largest effects on the start date. The relative contributions of temperature and precipitation on the start date of the vegetation growing season differed and depended on the elevation. Below the 4400-m elevation, precipitation played a dominant role, whereas above 4000 m, the dominant role was played by temperature. In contrast to the trend of the entire TP, at a given elevation level, the interannual trend of the SOS in the TRSR advanced at a rate of approximately 1 day more than in the TP, which suggests that the vegetation in the TRSR area is the most sensitive vegetation on the plateau.