Changing Spring Phenology Dates in the Three-Rivers Headwater Region of the Tibetan Plateau during 1960–2013

Shuang YU,Jiangjiang XIA,Zhongwei YAN,and Kun YANGKey Laboratory of Regional Climate-Environment for East Asia,Institute of Atmospheric Physics,

Chinese Academy of Sciences,Beijing 100029,China

2Key Laboratory of Tibetan Environment Changes and Land Surface Processes,Institute of Tibetan Plateau Research,

Chinese Academy of Sciences,Beijing 100107,China

3University of Chinese Academy of Sciences,Beijing 100049,China

Changing Spring Phenology Dates in the Three-Rivers Headwater Region of the Tibetan Plateau during 1960–2013

Shuang YU1,3,Jiangjiang XIA?1,Zhongwei YAN1,3,and Kun YANG2,31Key Laboratory of Regional Climate-Environment for East Asia,Institute of Atmospheric Physics,

Chinese Academy of Sciences,Beijing 100029,China

2Key Laboratory of Tibetan Environment Changes and Land Surface Processes,Institute of Tibetan Plateau Research,

Chinese Academy of Sciences,Beijing 100107,China

3University of Chinese Academy of Sciences,Beijing 100049,China

The variation of the vegetation growing season in the Three-Rivers Headwater Region of the Tibetan Plateau has recently become a controversial topic.One issue is that the estimated local trend in the start of the vegetation growing season(SOS)based on remote sensing data is easily affected by outliers because this data series is short.In this study,we determine that the spring minimum temperature is the most influential factor for SOS.The significant negative linear relationship between the two variables in the region is evaluated using Moderate Resolution Imaging Spectroradiometer–Normalized Difference Vegetation Index data for 2000–13.We then reconstruct the SOS time series based on the temperature data for 1960–2013.The regional mean SOS shows an advancing trend of 1.42 d(10 yr)-1during 1960–2013,with the SOS occurring on the 160th and 151st days in 1960 and 2013,respectively.The advancing trend enhances to 6.04 d(10 yr)-1during the past 14 years.The spatiotemporal variations of the reconstructed SOS data are similar to those deduced from remote sensing data during the past 14 years.The latter exhibit an even larger regional mean trend of SOS[7.98 d(10 yr-1)]during 2000–13.The Arctic Oscillation is found to have significantly influenced the changing SOS,especially for the eastern part of the region,during 2000–13.

start of growing season,normalized difference vegetation index,spring minimum temperature,Three-Rivers Headwater Region,Arctic Oscillation

1.Introduction

Vegetation plays a key role in the terrestrial ecosystem.For example,changes in the vegetation growing season influence crop yields(Chmielewski et al.,2004;Pengetal.,2004),water cycling(Peterson et al.,2002；Barnett et al.,2005)and carbon cycling(Piao et al.,2007,2011a),and hence provide feedback to climate change(Shen et al.,2015).Therefore,it is important to monitor the interannual and seasonal variations of vegetation for eco-environmental assessment(Chen et al.,2011a).The Three-Rivers Headwater Region(TRHR)is in the hinterland of the Tibetan Plateau,which fosters the Yangtze,Yellow and Lancang rivers with important ecological functions(Hu et al.,2011;Luo et al.,2014).Detection of shifts in the TRHR growing season has become an increasingly interesting topic from the perspective of climate change under global warming(Liu et al.,2004;Zeng et al.,2011).

In general,the annual variation of the vegetation growing season in the mid-to-high latitudes depends mainly on the date of the start of the vegetation growing season(SOS;Linderholm,2006),as opposed to the end of growing season(Piao et al.,2007;Piao et al.,2011a;Song et al.,2011).In recent years,studies of vegetation seasons in the TRHR have focused on the SOS dates(Song et al.,2011;Wang et al.,2013a;Liu et al.,2014b).These studies often used the phenological definition of growing seasons based on satellite remote sensing data of the normalized difference vegetation index(NDVI),owing to theirfi ne spatial resolution(Liu et al.,2004;Yu et al.,2010;Ding et al.,2013;Zhang et al.,2013a).However,the satellite data record is short for studying long-term trends of change.Many studies have investigated the relationship between SOS and other variables such as temperature(Root et al.,2003;Piao et al.,2006;Jeong et al.,2011),precipitation(Piao et al.,2011b;Shen et al.,2014),snow depth(Zhang et al.,2013b)and altitude(Shen et al.,2014).It has been evidenced that the negative correlation between the spring temperature and SOS is particularly significant in the TRHR(Piao et al.,2011a;Tian,2015).de Beurs and Henebry(2008)reported that the negative correlation between SOS in the Asian part of northern Eurasia and atmospheric circulation such as the Arctic Oscillation(AO)is significant.Whether and how the AO exerts influence on SOS in a particular region such as the TRHR remains an interesting and open topic.

It is notable that previous results regarding the trend of SOS during the last decade on the Tibetan Plateau remain controversial.For example,some researchers have reported a continuous advancement in SOS on the Tibetan Plateau during 1982–2011(Song et al.,2011;Zeng et al.,2011;Zhang et al.,2013a),whereas others have argued that there is no evidence to indicate such a conclusion(Shen et al.,2013).In particular,Shen(2011)reported an SOS-delaying trend from 1998 to 2003 and an advancement from 2003 to 2009,based on Système Pour l’Observation de la Terre Vegetation(SPOTVGT)NDVI data.

Three issues are relevant to this controversy.First,the quality of the Advanced Very High Resolution Radiometer(AVHRR)NDVI data(since 1982)is not as good as that of SPOT-VGT NDVI data(since 1998)and Moderate Resolution Imaging Spectroradiometer(MODIS)data(since 2000),as discussed by Zhang et al.(2013a),particularly for alpine steppe and alpine meadow regions during 2001–06.Second,although better-quality NDVI data from MODIS and SPOTVGT have been available since 2000 and 1998 respectively(Shen,2011;Song et al.,2011;Zhang et al.,2013a;Shen et al.,2014),the time series of these datasets are still too short for robust estimation of long-term trends.The estimate of a trend in vegetation activity can be easily biased by individual outliers(Zhang et al.,2013b).Having removed these outliers,Zhang et al.(2013b)reported that the advancing trend of SOS in the study region based on SPOT-VGT NDVI data during 1998–2011 is statistically significant.Third,the SOS dates defined using different methods can differ considerably,and hence can influence the estimate of a trend in a time series.In fact,even if based on the same AVHRR dataset,large uncertainty remains in estimating the phenological trend in a region using different methods.The trend can range from 0.4 to 1.9 d(10 yr)-1in China(Cong et al.,2013).

The main objective of the present study is to incorporate long-term series of local temperature records to investigate the long-term trend of SOS in the TRHR during 1960–2013 by using a novel method of time series analysis to define SOS,and understand the changes in SOS during the last decade as part of the long-term series of SOS.A preliminary analysis of the effect of the AO on the long-term series of SOS in the TRHR is made to understand the possible driving mechanism via large-scale atmospheric circulation.The data and methods used are described in section 2;the results are illustrated in section 3;and the conclusions and further discussion are given in section 4.

2.Data and methods

2.1.Data

The daily mean/maximum/minimum surface air temperature(Tmean/Tmax/Tmin)records from 28 National Standard Stations in the TRHR(31°39′–36°16′N,89°24′–102°23′E)during the period 1960–2013 are used in this study.The altitudes of all the stations are between 1800 m and 4500 m(Fig.1).This dataset has been homogenized and updated by Li and Yan(2009).

The NDVI data used in this study is the MODIS-NDVI product MODI3C1,which can be retrieved from the online Data Pool of the NASA Land Processes Distributed Active Archive Center at https://lpdaac.usgs.gov/dataset discovery/modis/modis products-table/mod13c1.This dataset has a spatial resolution of 0.05°×0.05°and a temporal resolution of 16 days for the period 2000–13.The average NDVI values within a box area of 1°latitude×1°longitude,centered at each temperature station,are applied.

The AO index is defined as the mean deviation from the average sea level pressure measured throughout the Northern Hemisphere at longitudes north of 20°N(Thompson and Wallace,1998).It is obtained from the NOAA Climate Prediction Center(www.cdc.noaa.gov/correlation/ao.data).The mean AO data in winter(December–January–February)and spring(March–April–May)during 1960–2013 are used.

2.2.Method

2.2.1.Determining the SOS date

Several methods can be used to determine the SOS date from the seasonal variation in NDVI for one year(Lloyd,1990;Reed et al.,1994;White et al.,1997;Yu et al.,2003;Zhang et al.,2003).In the present study,the SOS date is defined in a manner similar to that reported by Piao et al.(2006).First,the 14-year-average seasonal cycle of NDVI is calculated for the entire study area.Second,a threshold of NDVI is defined at the time when the rate of change in NDVI reaches its maximum in the mean annual cycle.For each year,the SOS date is defined as the first day when the value of NDVI exceeds this threshold.The rate of change in NDVI is calculated as

where t is time(temporal resolution of 16 days).Piao et al.(2006)analyzed biweekly NDVI time series data from January to September to determine the SOS dates for each year.We apply the ensemble empirical mode decomposition(EEMD)method to filter out the seasonal cycles from the NDVI series,rather than the polynomial function used by Piao et al.(2006).EEMD is a type of time-varying algorithm used for analyzing nonlinear and nonstationary time series such as those of climate(Huang et al.,1998;Huang and Wu,2008;Qian et al.,2010;Yan et al.,2011;Xia and Yan,2014).In this study,we apply the EEMD method to extract the seasonal component of NDVI datasets and smooth the NDVI data within 16-day windows.The specific steps are described in Qian et al.(2010).Figure 2 shows a betterfit of the multi-year average NDVI seasonal curve based on EEMD than that obtained by the six-degree polynomial function,which influences the determination of the SOS date. The procedure and algorithm of EEMD are summarized in Qian et al.(2010),Wu and Huang(2011),and Huang and Wu(2008).

Fig.1.Temperature observations recorded at 28 National Standard Stations.The 17 green points indicate stations satisfying the criterion for reconstruction outlined in section 2.2.2.The 8 yellow points indicate the stations used after further data processing outlined in section 3.2.2.The 3 red points indicate the stations eliminated in section 3.2.2.

Fig.2.Multi-year average NDVI seasonal cycle fitted using EEMD(red)and a six-degree polynomial function(green)for station 52868.The blue points indicate the original NDVI data.The horizontal black line indicates the threshold of NDVI at the time when the rate of NDVI change reaches its maximum in the mean annual cycle.The vertical lines show the SOS based on the EEMD fitting(red)and the six-degree polynomial function fitting(green)multi-year average NDVI seasonal cycle.The enlarged portion of the image shows that the SOS based on EEMD(red)is closer to the original NDVI data(blue points).

2.2.2.Establishing a linear relationship between spring temperature and SOS

In general,the pre-season temperature has a strong correlation with SOS(Yu et al.,2003;Piao et al.,2011a;Xia et al.,2013,2015).For each station,we choose the minimum value of the SOS date during 2000–13 as the base day.Similar to Yu et al.(2003),we use the term “pre-season”to refer to a period before the base day,and we analyze the correlation between SOS and the mean temperature of the preseason period.We set the length of the pre-season(N)from 1 to 100.Then,by using the least-squares method,we calculate the linear regression coefficient between SOS and the pre-season temperature with different N values for the period 2000–13 for each station.

The results show that the number of stations with significant correlation(P＜0.1)is largest when N is set to 40.We suggest that,for each station,if SOS and the pre-season temperature,with N from 30 to 50,show a significant linear relationship under the 0.1 significance level,as defined as the criterion for reconstruction,the station qualifies for reconstructing the SOS time series.

In addition,we define N with the largest correlation coeffi cient as the optimal pre-season length,which is used to reconstruct the SOS time series.Thus,the optimal pre-season length can differ among stations.

2.2.3.Reconstructing the 50-year series of SOS

The regression coefficient(a)of SOS onto the mean temperature during the optimal pre-season temperature is calcu-lated for each station for the period 2000–13 as

where D0is the mean SOS date depending on T0,the mean optimal pre-season temperature;δT is the temperature anomaly with respect to T0;and D is the SOS date corresponding to the temperature.On the basis of these regression coefficients and the optimal pre-season temperature records for 1960–2013,we can reconstruct the SOS time series for the same period for a given location near each station.

3.Results

3.1.Tminas the main factor influencing SOS

We calculate the correlation coefficient between SOS and the optimal pre-season temperature based on Tmin,Tmaxand Tmeanrecords.The Tmin-based optimal pre-season temperature index results in the most negative correlation with SOS for most of the stations(Fig.3a).Figure 3b shows that the linear correlation between regional mean SOS and the optimal pre-season minimum temperature is significant.Figure 3c illustrates that the regional mean spring minimum temperature during 1960–2013 has an increasing trend of 0.30°C(10 yr)-1.The warming trend during 2000–13 enhances up to 0.92°C(10 yr)-1and is significant at most stations.The warming trend during the last 14 years corresponds well with the rapid advancement of SOS during the same period.

Of all stations,17,or 61%(Fig.1),show a significant(P＜0.1)linear relationship between SOS and the optimal pre-season Tmin.Thus,for these 17 stations,the SOS time series are reconstructed on the basis of the optimal pre-season temperature and regression coefficient a(Table 1),using the method described in section 2.2.3.

Table 1.Linear regression between SOS and spring temperature at 17 stations.Triple,double and single asterisks indicate the correlation is significant at the 0.01,0.05 and 0.1 confidence level,respectively.

Fig.3.(a)Boxplot of the correlation coefficients between SOS and the pre-season temperature for 28 stations in the TRHR during 2000–13,based on Tmin,Tmaxand Tmeanseparately.The bottom and top of the boxes are the lower and upper quartiles,respectively;the bar near the middle of the box represents the median;and the solid points represent the outliers.The horizontal line indicates the significant criterion(P＜0.1)of-0.457;the degrees of freedom is 12.(b)Linear correlation between mean SOS and mean optimal pre-season Tminaveraged over 17 stations.(c)Regional mean spring Tminseries,marked with a linear trend during 2000–13(dotted line)and 1960–2013(solid lines).

3.2.Factors influencing the SOS–temperature relationship

It should be noted that 11 stations do not qualify for the aforementioned reconstruction analysis. Considering the limited number of stations in the TRHR,we utilize as many stations as possible to obtain relatively reasonable results.We explore the reason for the lack of an SOS–temperature relationship for these stations,and we examine the potential of utilizing as many data sources as possible.

3.2.1.Effect of outliers

Fig.4.The(a)Optimal pre-season temperature(Tmin;52943,for example)and(b)SOS(56144,for example).Linear trends are statistically significant(P＜0.1)after removing the data for 2002.

Zhang et al.(2013b)reported that the outliers of NDVI in 1998 and 2001 might considerably influence the SOS trend estimates.In addition,outliers may also influence the regression between the optimal pre-season temperature and SOS for the 11 unqualified stations.Figure 4 shows an example that the trends of both temperature and SOS are statistically significant at P＜0.1 after removing the year 2002.Having analyzed the correlation between the optimal pre-season temperature and SOS after removing each year during 2000–13,we determined that 4(52856,56064,56125 and 56146)of the 11 stations exhibit significantly improved correlation up to the criterion for reconstruction after removing the year 2002(Figs.5a–d).In contrast,the correlation is not enhanced after removing any one of the other years.This suggests possible bad records in either station temperature or remote sensing data in the year 2002.We therefore calculate the regression coefficients and reconstruct the SOS time series for these four stations back to the 1960s,according to the revised relation.To facilitate subsequent discussion,this method is referred to as Method 1.

3.2.2.Effect of local spring temperature data quality

Previous studies have reported that observed climate data may be fl awed because of artificial error,transference and changes to instrument types(Li and Yan,2009;Li and Yan,2010).In this study,we define each of the 11 unqualified stations as a base station,and we analyze the relationship between the SOS of the base station and the optimal pre-season temperature of each of the 28 stations.Taking station 52825 as an example(Fig.6),the correlation coefficients between the SOS at 52818 and the optimal pre-season temperature at adjacent stations,such as 52908,are higher and meet the criterion for reconstruction,whereas that between the SOS and local optimal pre-season temperature at the same station(52818)does not.This result could be due to temporally fl awed local observation records at this station.Therefore,we use the temperature data at the adjacent station,52908,which exhibits the highest correlation,to replace those at 52818 for calculating the regression coefficients and reconstructing the long-term SOS series during 1960–2013 for station 52818.We refer to this method as Method 2.Similarly,we apply this method to 7 out of the 11 unqualified stations to reconstruct the local long-term SOS series(Table 2).

Having applied the aforementioned two methods for excluding the influence of bad data and choosing the processing method with higher correlation for each station,we can add eight more stations for reconstructing the SOS series back to 1960,based on Method 2 for seven stations and Method 1 for one station(Table 2).The correlation between the SOS and pre-season Tmindata for the eight stations is improved and meets the criterion for reconstruction after further data processing(Fig.5a–h).Three stations(55299,56137 and 56144)still do not qualify,which implies other factors of influence beyond the scope of consideration in the present analysis.Nevertheless,the results of the 25 chosen stations(Tables 1 and 2)can form an effective base for an overview of the changes in SOS in the TRHR during 1960–2013.

3.3.Changes in SOS for the period 1960–2013

Fig.5.Improvement in the reconstruction-required relation for eight stations based on further data processing.Method 1 involves removing outliers in 2002 and Method 2 uses temperature data at a neighboring station.In each panel,the curves show the correlation coefficients between SOS and pre-season temperature,in which the length of pre-season changes from 1 to 100 days.The horizontal line indicates the(P＜0.1)significant criterion of the correlation coefficient.

Table 2.Linear regression between SOS and spring temperature at eight stations.Triple asterisks indicate the correlation is significant at the 0.01 confidence level,respectively.Method 1 and Method 2 are two solutions used to compensate for unreliable data.Method 1 removes outliers in 2002 and Method 2 uses the temperature data from a neighboring station.

Fig.6.Correlation coefficients between the SOS of base station 52825 and the pre-season Tminat 52825(red solid line),52908(red dotted line)and the other 26 stations(black dotted lines)in the TRHR during 2000–13.Horizontal red line indicates the significant criterion(P＜0.1)for the correlation coefficients.

Using the method described in section 2.2.3,the SOS time series for the period 1960–2013 are reconstructed for the 25 chosen stations.The geographical pattern of mean SOS dates based on the reconstructed data is consistent with that based on MODIS-NDVI during 2000–13(r=0.999,P＜0.05,where r is the spatial correlation coefficient;Figs.7a and c).Moreover,the correlation between the spatial distribution of the SOS trends based on MODIS-NDVI and the reconstruction is significant(r=0.628,P＜0.05).In particular,the stations in the central,southwest and northwest regions of the study area have higher spatial correlation,while those in the northeast and southeast regions are weaker(Figs.7b and d).These results reinforce the reasonability of the reconstructed SOS for most stations in the TRHR.

Tables 1 and 2 show the equations of linear regression between the pre-season Tminand SOS for the chosen 25 stations.Most of the stations show an advancing SOS trend for the period 1960–13.Figure 7e shows that the geographical pattern of mean SOS has generally earlier dates in the east in early April than those in the west by the end of May in the TRHR during 1960–2013.These results are consistent with the results of 2000–13(r=0.992,P ＜ 0.05).The spatial distribution of the SOS trends during the last half-century is also consistent with that during the last decade(r=0.669,P＜0.05),particularly in the central and northeast regions of the study area(r=0.878,P＜0.05).The SOS trend maxima in the study area are distributed mostly in the east,at 1.8–2.9 d(10 yr)-1,in contrast to those in the central and southwest regions at 0.3–1.3 d(10 yr)-1during 1960–2013(Fig.7f).

As shown in Fig.8,the regional mean SOS series during 1960–2013 indicate an advancing trend of 1.42 d(10 yr)-1(R2=0.357,P＜0.05,where R2is the coefficient of determination),with SOS occurring on approximately the 160th day of the year in 1960 and the 151st day in 2013.The reconstructed SOS series for the region is similar to that based on the remote sensing data for the last 14 years,although the latter appears to be more variable from year to year.The regional mean advancing trend of SOS in the TRHR during the last 14 years is quite large,supporting the conclusion of Zhang et al.(2013a).The regional mean reconstructed SOS series exhibits an advancing trend of 6.04 d(10 yr)-1during 2000–13,with the SOS occurring on approximately the 158th day of the year in 2000 and the 152nd day in 2013.The trend based on the remote sensing data is even larger,at 7.98 d(10 yr)-1.

3.4.Effects of the AO on the regional SOS trends

According to a sequential Mann–Kendall test(α =0.05),a significant jump in spring temperature in the TRHR is apparent in the late 1980s.This jump is similar in timing to the increase in the AO index for February–April.Negative correlation between the AO index for spring(especially for March)and regional SOS is found for the last 10 years,with the largest correlation coefficient(r=-0.260,P＞0.1)for March,especially for the eastern part of the region(r=-0.496,P＜ 0.1),during 2000–13(Fig.9a).The result suggests that the AO could exert its influence but is hardly a dominant factor for SOS in the study region.The influence of the AO exists mainly for winter temperature.The positive correlation between the winter AO index and regional spring surface air temperature is not significant,though the correlation with winter temperature is significant during 1960–2013(r=0.268,P＜0.1)(Fig.9b).When the AO index is in a high-value period,the mean surface air temperature in winter in the TRHR is also high.This preliminary relationship is useful for exploring the possibility of a link between SOS in the TRHR and large-scale atmospheric circulation.

Fig.7.Spatial distributions of average SOS dates and linear trends based on(a,b)MODIS-NDVI during 2000–13,(c,d)reconstruction during 2000–13 and(e,f)reconstruction during 1960–2013.Black circles indicate stations with non-significant advancing trends.Dots show the stations with significant trends.Units for SOS date:order number of the day in a year.

4.Conclusions and discussion

In the TRHR,the spring minimum temperature is the most influential factor for SOS,and the negative correlation is significant during the last decade.The increasing temperature during spring is key to enhancing the enzyme activity of vegetation and therefore accelerating this phenological process(Wang et al.,2010).Specifically,it leads to an earlier occurrence of photoinhibition,which takes place at low temperature,and hence causes earlier green-up of the vegetation(Wang et al.,2013b).The present result is useful for studying growing season and terrestrial ecosystem responses to climate change.In this study,we reconstruct SOS time series based on temperature data for 1960–2013 and investigate the spatiotemporal change of SOS in the TRHR during the last 50 years.The regional mean SOS shows an advancing trend of 1.42 d(10 yr)-1during 1960–2013,which enhances up to 6.04 d(10 yr)-1during the last 14 years.These results support the conclusion of Zhang et al.(2013a)and Song et al.(2011),whereas trends of delay(Yu et al.,2010)and reversal(Shen,2011)are not found during the last 10 years.The methods used in our study reduce the impact of individual outliers in short series of satellite data when calculating the trend of SOS,and the results are meaningful when attempting to clarify controversies in the SOS trend in the TRHR during the last decade.

Fig.9.(a)Linear correlation between mean SOS and AO index in March over most stations in the eastern part of the TRHR region during 2000–13.(b)Linear correlation between winter AO index and winter Tminof the TRHR region during 1960–2013.

Fig.8.Interannual variations in the reconstructed SOS for the entire study area from 1960 to 2013(black dotted curve).The black(red)line indicates the linear trend of the reconstructed SOS for the period 1960–2013(2000–13).Red points indicate the SOS dates calculated on the basis of MODIS-NDVI datasets.The red dotted line indicates the linear trend.

In addition,we attempt to identify possible links between SOS in the TRHR and large-scale atmospheric circulation via the AO.Preliminary results indicate a negative correlation between the AO index for spring(especially March)and SOS for the last 10 years,with the correlation particularly significant in the eastern part of the region.This builds upon the conclusion of de Beurs and Henebry(2008),albeit the relatively weak correlation implies more factors of influence are involved in this region.The possibility of teleconnection mechanisms involving multiple factors is deserving of further study.

In short,the present study establishes a linear relationship between pre-season temperature and SOS in the TRHR.The correlation is significant in most of the study region,and the reconstructed results reflect the long-term SOS trend in the TRHR during the last 50 years.Further research using accurate spring phenology models would still beneficial.It also remains interesting to study the synergistic effect of multiple factors on SOS,such as temperature(Root et al.,2003;Piao et al.,2006;Jeong et al.,2011),precipitation(Piao et al.,2011b;Shen et al.,2014)and the vegetation component(Chen et al.,2011b),as well as atmospheric circulation.The effects of climate variables on spring phenology might be discontinuous and synergistic(Kong et al.,2017)and are also deserving of further study.

Acknowledgements.This work was supported by the National Key Research and Development Program of China(Grant Nos.2016YFA0600400and2016YFA0602500).This work was also supported by the Open Research Fund of the Key Laboratory of Tibetan Environmental Changes and Land Surface Processes,Chinese Academy of Sciences,and the National Natural Science Foundation of China(Grant No.41405082).

Barnett,T.P.,J.C.Adam,and D.P.Lettenmaier,2005:Potential impacts of a warming climate on water availability in snow-dominated regions.Nature,438(7066),303–309,https://doi.org/10.1038/nature04141.

Chen,H.,Q.A.Zhu,N.Wu,Y.F.Wang,and C.-H.Peng,2011a:Delayed spring phenology on the Tibetan plateau may also be attributable to other factors than winter and spring warming.Proceedings of the National Academy of Sciences of the United States of America,108(19),E93,https://doi.org/10.1073/pnas.1100091108.

Chen,H.L.,Y.J.Liu,Z.X.Du,Z.Y.Liu,and C.H.Zou,2011b:The change of growing season of the vegetation in Huanghe-Huaihe-Haihe region and its responses to climate changes.Journal of Applied Meteorological Science,22(4),437–444,https://doi.org/10.3969/j.issn.1001-7313.2011.04.006.(in Chinese with English ABSTRACT)

Chmielewski,F.-M.,A.Müller,and E.Bruns,2004:Climate changes and trends in phenology of fruit trees and field crops in Germany,1961-2000.Agricultural and Forest Meteorology,121(1–2),69–78,https://doi.org/10.1016/S0168-1923(03)00161-8.

Cong,N.,T.Wang,H.J.Nan,Y.C.Ma,X.H.Wang,and R.B.Myneni,2013:Changes in satellite-derived spring vegetation green-up date and its linkage to climate in China from 1982 to 2010:A multi method analysis.Global Change Biology,19(3),881–891,https://doi.org/10.1111/gcb.12077.

deBeurs,K.M.,and G.M.Henebry,2008:Northernannular mode effects on the land surface phenologies of northern Eurasia.J.Climate,21(17),4257–4279,https://doi.org/10.1175/2008 JCLI2074.1.

Ding,M.J.,Y.L.Zhang,X.M.Sun,L.S.Liu,Z.F.Wang,and W.Q.Bai,2013:Spatiotemporal variation in alpine grassland phenology in the Qinghai-Tibetan plateau from 1999 to 2009.Chinese Science Bulletin,58(3),396–405,https://doi.org/10.1007/s11434-012-5407-5.

Hu,M.Q.,F.Mao,H.Sun,and Y.Y.Hou,2011:Study of normalized difference vegetation index variation and its correlation with climate factors in the three-river-source region.International Journal of Applied Earth Observation and Geoinformation,13(1),24–33,https://doi.org/10.1016/j.jag.2010.06.003.

Huang,N.E.,and Z.H.Wu,2008:A review on Hilbert-Huang transform:Method and its applications to geophysical studies.Rev.Geophys.,46(2),RG2006,https://doi.org/10.1029/2007RG000228.

Huang,N.E.,and Coauthors,1998:The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis.Proc.Roy.Soc.A,454(1971),903–995,https://doi.org/10.1098/rspa.1998.0193.

Jeong,S.-J.,C.-H.Ho,H.-J.Gim,and M.E.Brown,2011:Phenology shifts at start vs.end of growing season in temperate vegetation over the northern hemisphere for the period 1982-2008.Global Change Biology,17(7),2385–2399,https://doi.org/10.1111/j.1365-2486.2011.02397.x.

Kong,D.D.,Q.Zhang,W.L.Huang,and X.H.Gu,2017:Vegetation phenology change in Tibetan Plateau from 1982 to 2013 and its related meteorological factors.Acta Geographica Sinica,72(1),39–52,https://doi.org/10.11821/dlxb 201701004.(in Chinese with English ABSTRACT)

Li,Z.,and Z.-W.Yan,2009:Homogenized daily mean/maximum/minimum temperature series for China from 1960-2008.Atmospheric and Oceanic Science Letters,2(4),237–243,https://doi.org/10.1080/16742834.2009.11446802.

Li,Z.,and Z.W.Yan,2010:Application of multiple analysis of series for homogenization to Beijing daily temperature series(1960-2006).Adv.Atmos.Sci.,27(4),777–787,https://doi.org/10.1007/s00376–00909052-0.

Linderholm,H.W.,2006:Growing season changes in the last century.Agricultural and Forest Meteorology,137(1–2),1–14,https://doi.org/10.1016/j.agrformet.2006.03.006.

Liu,G.H.,Q.H.Tang,X.C.Liu,J.H.Dai,X.Z.Zhang,Q.S.Ge,and Y.Tang,2014b:Spatiotemporal analysis of ground-based woody plant leafing in response to temperature in temperate eastern China.International Journal of Biometeorology,58(7),1583–1592,https://doi.org/10.1007/s00484-013-0762-8.

Liu,L.M.,Y.T.Liang,H.Y.Ma,and J.Huang,2004:Relationship research between MODIS-NDVI and AVHRR-NDVI.Geomatics and Information Science of Wuhan University,29,307–310,https://doi.org/10.3321/j.issn:1671-8860.2004.04.006.(in Chinese with English ABSTRACT)

Liu,X.F.,X.F.Zhu,W.Q.Zhu,Y.Z.Pan,C.Zhang,and D.H.Zhang,2014a:Changes in spring phenology in the threerivers headwater region from 1999 to 2013.Remote Sensing,6(9),9130–9144,https://doi.org/10.3390/rs6099130.

Lloyd,D.,1990:A phenological classification of terrestrial vegetation cover using shortwave vegetation index imagery.Int.J.Remote Sens.,11(12),2269–2279,https://doi.org/10.1080/01431169008955174.

Luo,C.F.,J.Wang,M.L.Liu,and Z.J.Liu,2014:Analysis on the change of grassland coverage in the source region of three rivers during 2000-2012.IOP Conference Series:Earth and Environmental Science,Vol.17,No.1,p.012062.,IOP Publishing.

Peng,S.B.,and Coauthors,2004:Rice yields decline with higher night temperature from global warming.Proceedings of the National Academy of Sciences of the United States of America,101(27),9971–9975,https://doi.org/10.1073/pnas.0403720101.

Peterson,B.J.,R.M.Holmes,J.W.McClelland,C.J.V¨or¨osmarty,R.B.Lammers,A.I.Shiklomanov,I.A.Shiklomanov,and S.Rahmstorf,2002:Increasing river discharge to the arctic ocean.Science,298(5601),2171–2173,https://doi.org/10.1126/science.1077445.

Piao,S.L.,J.Y.Fang,L.M.Zhou,P.Ciais,and B.Zhu,2006:Variations in satellite-derived phenology in China’s temperate vegetation.Global Change Biology,12(4),672–685,https://doi.org/10.1111/j.1365-2486.2006.01123.x.

Piao,S.L.,P.Friedlingstein,P.Ciais,N.Viovy,and J.Demarty,2007:Growing season extension and its impact on terrestrial carbon cycle in the Northern Hemisphere over the past 2 decades.Global Biogeochemical Cycles,21(3),https://doi.org/10.1029/2006GB002888.

Piao,S.L.,M.D.Cui,A.P.Chen,X.H.Wang,P.Ciais,J.Liu,and Y.H.Tang,2011a:Altitude and temperature dependence of change in the spring vegetation green-up date from 1982 to 2006 in the Qinghai-Xizang plateau.Agricultural and Forest Meteorology,151(12),1599–1608,https://doi.org/10.1016/j.agrformet.2011.06.016.

Piao,S.L.,X.H.Wang,P.Ciais,B.Zhu,T.Wang,and J.Liu,2011b:Changes in satellite-derived vegetation growth trend in temperate and boreal Eurasia from 1982 to 2006.Global Change Biology,17(10),3228–3239,https://doi.org/10.1111/j.1365-2486.2011.02419.x.

Qian,C.,Z.H.Wu,X.B.Fu,and T.J.Zhou,2010:On multitimescale variability of temperature in China in modulated annual cycle reference frame.Adv.Atmos.Sci.,27,1169–1182,https://doi.org/10.1007/s00376–009-9121–4.

Reed,B.C.,J.F.Brown,D.VanderZee,T.R.Loveland,J.W.Merchant,and D.O.Ohlen,1994:Measuring phenological variability from satellite imagery.Journal of Vegetation Science,5(5),703–714,https://doi.org/10.2307/3235884.

Root,T.L.,J.T.Price,K.R.Hall,S.H.Schneider,C.Rosenzweig,and J.A.Pounds,2003:Fingerprints of global warming on wild animals and plants.Nature,421(6918),57–60,https://doi.org/10.1038/nature01333.

Shen,M.G.,2011:Spring phenology was not consistently related to winter warming on the Tibetan Plateau.Proceedings of the National Academy of Sciences of the United States of America,108(19),E91–E92,https://doi.org/10.1073/pnas.1018390108.

Shen,M.G.,Z.Z.Sun,S.P.Wang,G.X.Zhang,W.D.Kong,A.P.Chen,and S.L.Piao,2013:No evidence of continuously advanced green-up dates in the Tibetan Plateau over the last decade.Proceedings of the National Academy of Sciences of the United States of America,110(26),E2329,https://doi.org/10.1073/pnas.1304625110.

Shen,M.G.,G.X.Zhang,N.Cong,S.P.Wang,W.D.Kong,and S.L.Piao,2014:Increasing altitudinal gradient of spring vegetation phenology during the last decade on the Qinghai-Tibetan Plateau.Agricultural and Forest Meteorology,189-190,71–80,https://doi.org/10.1016/j.agrformet.2014.01.003.

Shen,M.G.,and Coauthors,2015:Evaporative cooling over the Tibetan Plateau induced by vegetation growth.Proceedings of the National Academy of Sciences of the United States of America,112(30),9299–9304,https://doi.org/10.1073/pnas.1504418112.

Song,C.-Q.,S.-C.You,L.-H.Ke,G.-H.Liu,and X.-K.Zhong,2011:Spatio-temporal variation of vegetation phenology in the Northern Tibetan Plateau as detected by MODIS remote sensing.Chinese Journal of Plant Ecology,35(8),853–863,https://doi.org/10.3724/SP.J.1258.2011.00853.(in Chinese with English ABSTRACT)

Thompson,D.W.J.,and J.M.Wallace,1998:The arctic oscillation signature in the wintertime geopotential height and temperature fields.Geophys.Res.Lett.,25(9),1297–1300,https://doi.org/10.1029/98GL00950.

Tian,L.Q.,2015:The research of green-up date in the Qinghai-Tibetan plateau based on remote sensing technology.M.S.thesis,Northwest A&F University.(in Chinese)

Wang,L.X.,H.L.Chen,Q.Li,and Y.D.Wu,2010:Research advances in plant phenology and climate.Acta Ecologica Sinica,30(2),447–454.(in Chinese with English ABSTRACT)

Wang,T.,S.S.Peng,X.Lin,and J.F.Chang,2013a:Declining snow cover may affect spring phenological trend on the Tibetan Plateau.Proceedings of the National Academy of Sciences of the United States of America,110,E2854–E2855,https://doi.org/10.1073/pnas.1306157110.

Wang,Y.T.,L.Q.Yu,F.G.Wang,N.Wang,and J.Liu,2013b:Effects of low-temperature stress on chlorophyll fluorescence parameters of alfalfa in returning green period.Chinese Journal of Grassland,35(1),29–34,https://doi.org/10.3969/j.issn.1673-5021.2013.01.005.(in Chinese with English ABSTRACT)

White,M.A.,P.F.Thornton,and S.W.Running,1997:A continental phenology model for monitoring vegetation responses to interannual climatic variability.Global Biogeochemical Cycles,11(2),217–234,https://doi.org/10.1029/97GB00330.

Wu,Z.H.,and N.E.Huang,2011:Ensemble empirical mode decomposition:a noise-assisted data analysis method.Advances in Adaptive Data Analysis,1(1),https://doi.org/10.1142/S1793536909000047.

Xia,J.J.,and Z.W.Yan,2014:Changes in the local growing season in eastern China during 1909-2012.SOLA,10(1),163–166,https://doi.org/10.2151/sola.2014-034.

Xia,J.J.,Z.W.Yan,and P.L.Wu,2013:Multidecadal variability in local growing season during 1901-2009.Climate Dyn.,41(2),295–305,https://doi.org/10.1007/s00382-012-1438-5.

Xia,J.J.,Z.W.Yan,G.S.Jia,H.Q.Zeng,P.D.Jones,and W.Z.Zhou,2015:Projections of the advance in the start of the growing season during the 21st century based on CMIP5 simulations.Adv.Atmos.Sci.,32(6),831–838,https://doi.org/10.1007/s00376-014-4125-0.

Yan,Z.W.,J.J.Xia,C.Qian,and W.Zhou,2011:Changes in seasonal cycle and extremes in China during the period 1960-2008.Adv.Atmos.Sci.,28,269–283,https://doi.org/10.1007/s00376-010-0006-3.

Yu,F.F.,K.P.Price,J.Ellis,and P.J.Shi,2003:Response of seasonal vegetation development to climatic variations in eastern central Asia.Remote Sensing of Environment,87(1),42–54,https://doi.org/10.1016/S0034-4257(03)00144-5.

Yu,H.Y.,E.Luedeling,and J.C.Xu,2010:Winter and spring warming result in delayed spring phenology on the Tibetan Plateau.Proceedings of the National Academy of Sciences of the United States of America,107(51),22 151–22 156,https://doi.org/10.1073/pnas.1012490107.

Zeng,H.Q.,G.S.Jia,and H.Epstein,2011:Recent changes in phenology over the northern high latitudes detected from multi-satellite data.Environmental Research Letters,6,045508,https://doi.org/10.1088/1748–9326/6/4/045508.

Zhang,G.L.,Y.J.Zhang,J.W.Dong,and X.M.Xiao,2013a:Green-up dates in the Tibetan Plateau have continuously advanced from 1982 to 2011.Proceedings of the National Academy of Sciences of the United States of America,110(11),4309–4314,https://doi.org/10.1073/pnas.1210423110.

Zhang,G.L.,J.W.Dong,Y.J.Zhang,and X.M.Xiao,2013b:Reply to Shen et al.:No evidence to show nongrowing season NDVI affects spring phenology trend in the Tibetan Plateau over the last decade.Proceedings of the National Academy of Sciences of the United States of America,110(26),E2330–E2331,https://doi.org/10.1073/pnas.1305593110.

Zhang,X.Y.,M.A.Friedl,C.B.Schaaf,A.H.Strahler,J.C.F.Hodges,F.Gao,B.C.Reed,and A.Huete,2003:Monitoring vegetation phenology using MODIS.Remote Sensing of Environment,84(3),471–475,https://doi.org/10.1016/S0034-4257(02)00135-9.

1 December 2016;revised 9 June 2017;accepted 14 August 2017)

:Yu,S.,J.J.Xia,Z.W.Yan,and K.Yang,2018:Changing spring phenology dates in the Three-Rivers Headwater Region of the Tibetan Plateau during 1960–2013.Adv.Atmos.Sci.,35(1),116–126,https://doi.org/10.1007/s00376-017-6296-y.

?Corresponding author:Jiangjiang XIA

Email:xiajj@tea.ac.cn

?Institute of Atmospheric Physics/Chinese Academy of Sciences,and Science Press and Springer-Verlag GmbH Germany 2018