HTML
--> --> -->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 their fine 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\`eme Pour l'Observation de la Terre Vegetation (SPOT-VGT) 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 SPOT-VGT 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)-1 in 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.1. Data
The daily mean/maximum/minimum surface air temperature (\(T_\rm mean/T_\max/T_\min\)) 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).
Figure1. 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.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.2. Method
2.2.1. Determining the SOS dateSeveral 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 \begin{equation} {\rm NDVI}_{{\rm ratio}}(t)=[{\rm NDVI}(t+1)-{\rm NDVI}(t)]/{\rm NDVI}(t) , \ \ (1)\end{equation}
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 better fit 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).
Figure2. 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 pre-season 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 coefficient 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 calculated for each station for the period 2000-13 as \begin{equation} D=a\delta T+D_{0}(T_{0}) , \end{equation} where D0 is 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.1. \(T_\min\) as the main factor influencing SOS
We calculate the correlation coefficient between SOS and the optimal pre-season temperature based on \(T_\min\), Tmax and T mean records. The \(T_\min\)-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)-1 and 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 \(T_\min\). 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.
Figure3. (a) Boxplot of the correlation coefficients between SOS and the pre-season temperature for 28 stations in the TRHR during 2000-13, based on \(T_\min\), Tmax and T mean separately. 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 \(T_\min\) averaged over 17 stations. (c) Regional mean spring \(T_\min\) series, marked with a linear trend during 2000-13 (dotted line) and 1960-2013 (solid lines).2
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
(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.
Figure4. The (a) Optimal pre-season temperature (\(T_\min\); 52943, for example) and (b) SOS (56144, for example). Linear trends are statistically significant (P<0.1) after removing the data for 2002.
Figure5. 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.3.2.2. Effect of local spring temperature data quality
Previous studies have reported that observed climate data may be flawed 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 flawed 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 \(T_\min\) data 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.
2
3.3. Changes in SOS for the period 1960-2013
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.
Figure6. Correlation coefficients between the SOS of base station 52825 and the pre-season \(T_\min\) at 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.
Figure7. 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.Tables 1 and 2 show the equations of linear regression between the pre-season \(T_\min\) and 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)-1 during 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 R2 is 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)-1 during 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.
2
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.
Figure8. 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.
Figure9. (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 \(T_\min\) of the TRHR region during 1960-2013.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.
