Spatiotemporal variations of forest phenology in the Qinling Mountains and its response to a critical temperature of 10°C

Abstract. Based on the fifth version MOD09A1 product of NASA LP DAAC (Land Processes Distributed Active Center) from 2001 to 2016, an enhanced vegetation index (EVI) time series dataset in the Qinling Mountains was reconstructed; it was used to extract the plant phenological parameters of the region by employing the maximal slopes of the alteration method and threshold method. The results show that the Whittaker filter method was better than other methods for the reconstruction of shrubland, farmland, and forest in the Qinling Mountains. Based on the reconstructed EVI, the characteristics of vegetation coverage, the start of growth season (SOG), the end of growth season (EOG), and the length of growth season (LOG) in the Qinling Mountains were all analyzed. There was an increasing trend of vegetation coverage in most regions (about 89.93% of the monitored areas) over the Qinling Mountains in the past decades, and the average phenological distribution pattern in the Qinling Mountains was closely related to the local hydrothermal conditions. From the high altitude mountainous areas to the farming areas, the SOG was gradually postponed, the EOG was gradually earlier, and the LOG gradually shortened. Furthermore, the time series variation of SOG, EOG, and LOG from 2001 to 2016 in the Qinling Mountains was also studied. The variation showed that the SOG shifted earlier, which was more prominent in high-altitude areas, while, for some southern and northern foothills altitude below 500 m and a few areas of the eastern Funiu Mountains, the SOG was delayed. The EOG shifted later, which was more apparent in the northern Qinling foothills and mid- to low-altitude areas of the eastern ranges, and the LOG was extended. Finally, the studies of correlation analysis between the plant phenology and a temperature threshold of 10°C showed that global warming was the major factor affecting the phytophenology of the Qinling Mountains, and the effects were concentrated mainly in the 1000- to 2000-m zones of the southern and northern Qinling Mountains. Moreover, the northern Qinling foothills showed a greater response to accumulated temperature than the southern foothills, and the low-altitude area of the eastern Funiu Mountains exhibited the lowest correlation with the accumulated temperature.


Introduction
Phenology refers to periodic changes of organisms after their long-term adaptation to climate conditions, thereby developing corresponding growth and development cycles. Phenology is the most intuitive indicator for the response of ecosystems to climate change. By studying forest phenology, we can understand the influences of climate change on the species composition and physiology of trees and on forest productivity, which are of vital importance for the conservation of biodiversity. As such, phenology is a hot topic in geology research, particularly large-scale vegetation surveillance, which is a key field in phenology that addresses a crucial scientific question of ecological and environmental studies. 1,2 Common methods of phenology surveillance include ground observation, remote sensing, and model-based surveillance. Among them, ground observation has high time precision and is easy to operate for small-scale phenology studies, but it has difficulty revealing phenology conditions at a regional scale. Satellites remote sensing technology can remedy defects of ground observation. Atzberger and Eilers, 3 Beck et al., 4 Chen et al., 5 Sakamoto et al., 6 Yves and Sobrino, 7 Wang et al., 8 White et al., 9 and Zhang et al. 10 used satellites data to monitor the vegetation activity and phenology over different areas. Since the threshold method was reported to extract phytophenological parameters in 1990, 11 remote sensing-based phenology surveillance has been applied in a variety of fields, including biodiversity, land cover classification, estimation of crop yield, and biomass assessment. Correspondingly, phenology has become an effective tool for studying terrestrial ecosystems. Many studies have been devoted to investigating approaches for extracting phenological parameters based on remote sensing data, such as the normalized difference vegetation index (NDVI)-based threshold method and the maximal slope of alteration method that are indicative of vegetative types. Meanwhile, the relationship between the phenological information of remote sensing monitoring and the actual observation data of terrestrial stations has received increasing attention. It was reported that, through the use of ground observation data as well as NVDI datasets of Landsat and MODIS, different sources of remote sensing results could be used to describe the observation data of ground phenology in an area. 12 A study by Kross et al. 13 indicated that an NDVI dataset with a time resolution of 10 or 15 days was suitable for extracting phenological information from ground vegetation. Atkinson et al. 14 compared four models for smoothing satellite sensor time-series data to estimate vegetation phenology. It has been demonstrated that over the last several decades, the plant growth seasons in the northern hemisphere have gradually extended, which leads to an increase in carbon sequestration in highlatitude areas. 15 In China, many studies show similar trends in most regions using ground observation and satellite monitoring, Fang et al. 16 and Zheng et al. 17 discussed the effect of global warming on plant phenological changes, and the results show that the response of ahead (or delay) of phenophase to increasing (or decreasing) of temperature was nonlinear. Zhang, 18 Lu et al., 19 Li et al., 20 He et al., 21 Fang et al., 22 Zhong et al., 23 Wang et al., 24 and Xie et al. 25 studied the trend of plant phenology and its responses to climatic change over different areas in China.
The Qinling Mountains lie in central China and are at the ecological transition belt between the warm temperature zone and subtropical zone. The mountain range, which is the geological boundary between southern and northern China, has two very different definitions. In the broad definition, the Qinling Mountains join the Kunlun Mountains in the west, connect the Minshan Mountains in the north, and travel eastward through Gansu, southern Shaanxi, to reach the Funiu Mountains in Henan in the east; in the narrow definition, Qinling refers to mountains between the Wei River and Han River in Shaanxi Province. In recent years, progressive urbanization and increased human activity have aggravated vegetation degeneration in some areas. 26,27 Although the alpine vegetation that is present at altitudes above 2800 m is barely affected by any human activity, it is highly impacted by climate change. Thus far, considerable effort has been dedicated to examining the phytophenological responses in areas of different altitudes in the Qinling Mountains, [28][29][30] and vegetation degeneration, its causes, and spatiotemporal changes of phenology in the Qinling Mountains have been analyzed. However, although spatiotemporal changes of phenology have been examined in the Qinling Mountains, there has been no thorough investigation into the causes of phenological changes and the resultant responses to climate change, which is due to the short or overly macroscopic time series. Here, we utilized satellite data from 2001 to 2016 to extract the phenological changes of different altitudes in the Qinling Mountains. In addition, we employed a biological threshold temperature of 10°C to analyze the causes of the phenological changes. As such, the study provides crucial evidence to understand the phytophenological response features at different altitudes that are under the dual influence of human activity and climate change.

General Conditions of the Study Region
The Qinling Mountains in a broad sense include 43 counties in Gansu, Shaanxi, and Henan provinces. 31 It stretches over 700 km, and there are many interconnected mountains, including Xiaolong Mountain, Taibai Mountain, Xiong'er Mountain, Huashan Mountain, and the Funiu Mountains. Overall, the orientation of the Qinling Mountains matches well with the 0°C isotherm of January. With an annual average temperature of 12°C to 14°C and precipitation of 700 to 900 mm, the mountain range has a rich coverage of vegetation, including forest and grassland in the high-altitude areas and farmland in low-altitude valleys. The central hills boast extraordinary altitude differences that accommodate a spectrum of vegetation types, which are the farming area (altitude 200 to 700 m), deciduous broad-leaved forest belt (altitude 700 to 1300 m), deciduous oak forest belt (altitude 1300 to 2300 m), birch forest belt (altitude 2300 to 2800 m), coniferous forest belt (altitude 2800 to 3200 m), and alpine shrub-meadow belt (above altitude 3200 m). The coniferous forest belt, which represents the high-altitude limit of the forest, is dominant with Farges' fir (Abies fargesii) forest and Shaanxi larch (Larix potaninii var. chinensis) forest. In addition, there is a little distribution difference spectrum of vegetation types between the southern and northern Qinling hills, which is attributed to the relevant climate differences. 32 To improve validation accuracy, nine classical ecological sample plots were selected based on high-resolution land use data ( Fig. 1) and ground survey data to represent the various vegetation types of the forest, shrubland, and farmland (Table 1). Among these,  Ninshan, with an altitude of 2630 m, is a national field observation station of a coniferous forest of the Qinling ecological system (Table 1).

Satellite Data
The fifth version MOD09A1 product from 2001 to 2016 was provided by NASA LP DAAC (Land Processes Distributed Active Archive Center); it included surface albedo bands 1 to 7, day of year (DoY) data, an 8-day time resolution, and a 500-m spatial resolution. In this study, the MODO09A1 product was used to generate pixels with the best observation quality in an 8-day interval. After atmospheric correction, the Qinling Mountains were found to be located in the H26V5 and H27V5 tiles. From 2001 to 2016, a total of 736 MOD09A1 images were generated, and all were subject to SIN GRID projection. To remove the influences of nonforest areas, we used the MOD12Q1 land cover data of 2012 to reveal that forest and farmland accounted for over 95% of the total area. The linear tendency method was used for the trend analysis; when B > 0, an increase of time T was associated with an increase of X; the value of B indicated the rate (or tendency) of ascent or descent [Eq. (1)]: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 5 0 1 where B is the slope of linear tendency, x is parameter of phenology, t is the year of phenology, and n ¼ 15.
The coefficient of variation (C v ) was used to evaluate the time series stability; it indicated a degree of dispersion such that the lower the C v , the lower the fluctuation [Eq. (2)]: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 3 9 9 where C v is the coefficient of enhanced vegetation index (EVI) variation, x i is the EVI of the i year, and x is the average of EVI from 2001 to 2016. The correlation coefficient (R) reflected the ability of the reconstructed time series curve to retain the features of the corresponding original curve; the greater the value is, the higher the reconstruction fidelity [Eq. (3)]: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 2 9 4 The standard error of the regression estimate (regression mean squared error, RMSE) was employed to describe the differential level between the individual time series; the lower the value is, the greater the representatively of the fitted value [Eq. (4)]: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 2 1 0 where EVI pi is reconstruction series data, EVI oi is original series data, and n is the number of samples.

Calculation of the Vegetation Index
This study analyzed the MODIS MOD13Q1 product with a 16-day interval, and it showed that when the NDVI of the high-vegetation coverage areas in the Qinling Mountains approached 0.8, the fluctuation decreased, which indicated that the chlorophyll responses of the vegetation reached saturation. This result was consisted with the previous report. 33 Fortunately, the EVI avoids the saturation issue of high-vegetation coverage areas and reduces the background noises of the atmosphere and soil. In addition, the EVI uses a time resolution of 8 days, which is considerably better than the 16-day resolution. As such, we utilized interactive data language along with the MRT tool of MODLAND to write the batch program, whereby the reflectivity data and DoY bands in the MODIS tiles of H26V5 and H27V5 were pieced together and projected. Correspondingly, the data were converted to coordinate projections of latitude and longitude in the WGS-84 coordinate system. The EVI was calculated using Eq. (5): ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 6 4 0 in which ρ 1 , ρ 2 , and ρ 3 represent the reflection values of the infrared, red, and blue bands, respectively, after atmospheric correction from the MODIS sensor, G (¼2.5) represents the gain factor, L (¼1) represents the soil-adjusted index, and C 1 (¼6.0) and C 2 (¼7.5) describe the influences of the blue-band correction on the red band. To further improve time precision, the DoY values of pixels were interpolated into day values during reconstruction of the time series, thereby avoiding the fluctuation of crests and troughs around the 8-day intervals.

Reconstruction of the EVI Time Series
To minimize the influences of random factors on the time series, a time series reconstruction was needed. [34][35][36][37][38][39] Common algorithms for this purpose include the filter smoothing method (e.g., the Savitzky-Golay filter 40 and Whittaker filter), 30 the function fitting method [also known as asymmetric Gaussian (AG)], double logistic (DL) function-fitting, and the frequency-domain transformation method [e.g., harmonic analysis of time series (HANTS), and wavelet analysis]. Different methods have different requirement stringencies on the original time series and are therefore suitable for different environmental conditions. In the Qinling Mountains, where there are considerable levels of data noise, the Whittaker, DL, HANTS, and AG methods were chosen to reconstruct the time series. Afterward, the resulting fidelity and antijamming capability were analyzed for individual sample plots to determine the optimal approach for each ecological type.

Quantitative comparison of the reconstruction algorithms
Based on the correlation performance and RMSE values of the four methods that were used to reconstruct the time series of shrubland, farmland, and forest ( Fig. 2) from 2001 to 2016, the results showed that the Whittaker method exhibited a better correlation performance than the DL and AG methods, whereas the HANTS method showed a relatively poor outcome. Also, an examination of the RMSE values showed that the Whittaker method showed the lowest deviation between the reconstructed and original EVI values, indicating that this method generated the best fitting performance, better than DL and AG, and the HANTS methods displayed a relatively poor outcome. Figure 3 showed the comparison of the time series curves of original EVI data from 2014 to 2016 reconstruction data derived from the four methods. The results showed that, in the shrub-land, the Whittaker, DL, AG, and HANTS methods displayed similar levels of fitting, but, in the peak value, the Whittaker method showed clearly superior antijamming performance over the other methods, as it delivered better fitting outcomes at points with considerable noise levels. In the farmland sample plots, the HANTS, AG, and DL methods generated fitting results that were relatively high in the early growth phase, although in the middle and peak growth stages, all four methods produced consistent fitting results that matched the original EVI time series curves. Fig. 3 The reconstructed EVI curves of the sample plots of (a) shrubland, (b) farmland, and (c) forest.

Comparison of the time series curves of several reconstruction algorithms
All three filtering methods except the Whittaker displayed excessive fitting in the crest period of 2016 and produced a single peak, whereas Whittaker-based fitting generated a dual-peak curve that was consistent with the original data. In the forest sample plots, the four methods generated similar fitting results in the early growth phase and end of growth season (EOG), but the HANTS produced relatively high fitting results in the crest period of 2014 to 2015 due to noise interference, as did the AG and DL in 2016. By contrast, the Whittaker generated a curve that was consistent with the original EVI time series curve. In general, the Whittaker method displayed the most stable performance among the four methods and is suitable for use in the Qinling Mountains, which is consistent with the finding by Atkinson. 14

Extraction of Phenological Parameters
The methods to extract phenological parameters include the threshold method, 40 maximal slope of alteration, 41 and the logistics function method. 42,43 The Qinling Mountains have extraordinary altitude differences, application of the threshold method alone may not adequately reflect topological influences. To avoid this problem, the threshold method was also used with the maximal slop of alteration to extract the phenological parameters. First, the Whittaker filter was employed to reconstruct the daily means of EVI data from 2001 to 2016; then the differences between adjacent points were calculated to generate an EVI difference curve. Subsequently, the maximal change value of the upswing (or downswing) curve is taken as the EVI threshold of the SOG (or EOG). Finally, year-by-year sequences of daily values were compared with the corresponding threshold value to generate the SOG and EOG, which were subtracted to generate the length of growth season (LOG).

Accuracy Validation of the Phenological Extraction in the Qinling Mountains
A traditional and straightforward approach to validating the accuracy of the phenological extraction is to adopt the occurrence dates of a given phenological phenomenon as the validation data. In an area of limited size, there is a certain level of similarity in the plant phenophase of a particular vegetation type. The Qinling Botanic Garden of Xi'an is located in the northern foothills of the mountain range and is representative of the local phenology. We therefore acquired typical phenological data of woody plants in the Xi'an station during 1963 to 2008 and adopted the daily average of their peak leafing period as the SOG and their leaf color-change period as the EOG. It was discovered that the peak leafing period was predominantly concentrated in the first 10 days of April (the 95th to 100th days of each year), whereas all the color-change periods were in the first 10 days of November (the 300th to 310th days). With regard to multiyear trends, it was reported that the peak leafing period advanced earlier by 1.5 day∕year and the color-change period was postponed by 1 day∕year. 24,44 These findings overall were corroborated by the results derived from the satellite data.

Characteristic of Vegetation Cover Changes of the Qinling Mountains From 2001 to 2016
The vegetation cover changes and spatial distribution of vegetation variability (Fig. 4) of the Qinling Mountains from 2001 to 2016 (Fig. 4) showed that, during this period, except for a few zones that displayed a decreased trend, the mountain range overall manifested increased vegetation coverage; the areas showing increased (or decreased) vegetation coverage accounted for 89.93% (and 11.07%) of the total area, respectively. Among the areas showing increased coverage, 38.57% displayed significant differences (P < 0.05). Notably, the southeastern part of the Qinling Mountains, where the vegetation cover had originally been low, registered a more prominent ascending trend than other parts. This phenomenon could be explained by the unfavorable local conditions that initially prompted policies of returning farmland to forest and grassland and culminated in a remarkable rise in forest concentrations. We then interrogated the vegetation cover data by incorporating the altitude factor and discovered that areas with increased vegetation coverage were mainly the deciduous broadleaf forest belt, with an altitude of 800 to 1500 m. These forests had vegetation coefficients of variation between 0.02 and 0.6 [ Fig. 4(b)], indicating considerable spatial variation. Among these forests, those in high-altitude zones or valleys that were the least impacted by human activity showed the highest vegetation stability. On the other hand, the mixed farmland-forest vegetation in the shallow mountains of northern Qinling and the surrounding areas of the mountainous residential communities, which were highly influenced by human presence, exhibited apparent variations [ Fig. 4(b)]. Furthermore, the alpine coniferous forests also exhibited a relatively high coefficient of variation as they accommodated great concentrations of Chinese larches (Larix chinensis) and Farges' firs (Abies fargesii), which are highly sensitive to climate change.

Spatial Distribution of the Phenological Average Value Over Qinling Mountains from 2001 to 2016
The spatial distribution of the phenological average value over the Qinling Mountains from 2001 to 2016 is shown in Fig. 5. The data showed that from the alpine areas to farmland areas, the SOG was gradually postponed (decreased in value), the EOG was gradually advanced (decreased in value), and the LOG was gradually shortened. These trends echoed the topography of the Qinling Mountains. The shallow mountainous areas of the northern and southern Qinling foothills and farmland area of the Funiu Mountains had relatively early SOG dates, which began during the first 10 days of March. By contrast, the SOG of the alpine boreal zones was relatively late and began between the 120th and 135th days of the year. The EOG dates mostly appeared between the 270th to 310th days of the year. The EOG was relatively early in the alpine boreal zones but relatively late in the shallow mountainous areas (commonly after the 300th day of the year). The LOG was generally between 150 and 200 days∕year and was relatively long in the low-altitude areas (over 180 days) but was relatively short in the high-altitude forest areas (between 150 and 170 days). The results also showed that both the SOG and LOG displayed apparent changes in response to rising altitude. However, the areas that were under strong human influence, such as the shallow mountainous areas of the northern Qinling foothills and the surrounding areas of the mountainous residential communities, exhibited no clear correlation between the altitude and phenological changes.

Interannual Variation of the Phenological over Qinling Mountains
The spatial distribution of the multiyear means of the SOG, EOG, and LOG in the Qinling Mountains from 2001 to 2016 is shown in Fig. 6. The results show that the SOG dates in the Qinling Mountains advanced earlier, including 23.45% of the pixels having P values <0.01 (Fig. 7), and about 95% of the study region showed an advancement of 1 to 2 days∕year [ Fig. 6(a)], which was more apparent in the high-altitude areas. Nevertheless, an SOG delay was found in the southern and northern low-altitude foothills and the central area of the Funiu Mountain, which overall matched the scope under significant human interference. The EOG dates in the Qinling Mountains were delayed [ Fig. 6(b)], including 21.45% of the pixels having P values <0.01 (figure not shown here). An EOG delay was apparent (2 to 3 days∕year) in the northern foothills and eastern mid-to low-altitude areas in the Qinling Mountains [ Fig. 6(b)].  During the 2001 to 2016 period, the LOG overall displayed an extending trend in the forest areas of the Qinling Mountains [ Fig. 6(c)], including 70% of the areas that were prolonged by 1 to 3 days∕year; in areas with a prolonged LOG, 51% of the pixels had P values <0.01 (figure not shown here). These areas with extended LOGs were mainly distributed in the central and eastern mid-to low-altitude regions. However, ∼30% of the areas in the Qinling Mountains witnessed a shortened LOG (−1 to 0 days∕year), mainly in low-altitude areas that are vulnerable to the influence of human activity.

Correlation Between the Phytophenology of the Qinling Mountains and the
Threshold Temperature of 10°C Environmental and climate alterations are the major factors contributing to phytophenological changes. Hydrothermal variations were investigated in the Qinling Mountains, showing that the mountain range was experiencing rising temperatures (measured in the spring or as a multiyear average) [45][46][47][48][49] and reduced precipitation 50,51 and that these changes mainly took place in the southern and northern foothills. An annual accumulated temperature above 10°C is a critical temperature for the majority of plant growth; usually there was an accumulation of daily average temperature above 10°C, which was defined as an active accumulated temperature. In this study, the annual accumulated temperature above 10°C was analyzed to understand the real reason of the variation of phenological changes in the Qinling Mountains from 2001 to 2016. The active accumulated temperature above 10°C, annual precipitation, and sunshine duration from 2001 to 2016 at 55 stations across the Qinling Mountains was analyzed. Furthermore, the means of the SOG, EOG, and LOG within a 5-km diameter range of each station were extracted. A correlation analysis was then performed for the two groups of data. The results showed that active accumulated temperature exhibited a higher correlation with the SOG, EOG, and LOG than the precipitation and sunshine duration; the three climate parameters had an order of response to the phenological indices of ≥10°C accumulated temperature > precipitation > sunshine duration. This was consistent with a previous study, 52 and, collectively, the findings support that global warming is a main factor impacting the phytophenogical changes in the Qinling Mountains. Among the three phenological parameters, the SOG overall showed a negative correlation with the ≥10°C accumulated temperature [ Fig. 8(a)]. There were 13 stations approved by significant testing, mainly distributed in the counties of the northern foothills (e.g., Qianyang, Fengxiang, Qishan, Fufeng, and Huayin) and southern foothills (e.g., Mianxian, Chenggu, Yangxian, Shiquan, Zhen'an, Liuba, Danfeng, and Shangnan), where the altitude is 1000 to 2000 m and the annual mean of the ≥10°C accumulated temperature is 4000 to 4500°C. In these areas, the rising temperature clearly advanced the SOG. By contrast, the SOG was poorly correlated with the ≥10°C accumulated temperature in the alpine areas, where the annual accumulated temperature was below 4000°C, as well as in the areas that were under significant human influence, where the annual accumulated temperature was above 4800°C. In general, the rising temperature showed no apparent effect on the SOG in these areas.
Both the EOG [ Fig. 8(b)] and LOG [ Fig. 8(c)] overall exhibited a positive correlation with the ≥10°C accumulated temperature. The LOG showed a higher correlation with the accumulated temperature than the EOG, with 14 stations approved by significant testing. These stations matched well with the stations where the SOG was significantly correlated with the accumulated temperature. Together, these results indicated that the areas with a greater response to the ≥10°C accumulated temperature were distributed in the southern and northern foothills, where the altitude was between 1000 and 2000 m and the annual accumulated temperature was between 4000°C and 4500°C, and that the northern foothills displayed greater responsiveness than the southern foothills. In addition, the alpine areas (with an accumulated temperature of <4000°C) had complicated responses to the phenological changes, which might be attributed to the extraordinary sensitivity of Chinese larches (Larix chinensis) and Farges' firs (Abies fargesii) to temperature changes. The phenological indices in the alpine areas showed no apparent responses to the active accumulated temperature as the climate parameter grew slowly in the alpine areas. Last, stations in the Funiu Mountains in the east, where the annual accumulated temperature is over 4900°C, also displayed a poor correlation with the threshold temperature, indicating that this region was not significantly impacted by this temperature condition.

Summary and Conclusions
In this study, the C5 MOD09A1 product from 2001 to 2016 provided by NASA LPACC was used to rebuild the EVI time series. Four reconstruction methods were compared; the maximal slop of alteration method and the threshold method were employed to extract the phenological parameters in the Qinling Mountains.
The results showed that the Whittaker filter was better in the reconstruction of sample plots for shrubland, farmland, and forest than the DL, AG, and HANTS methods. Additionally, the reconstructed curves of the Whittaker filter displayed the lowest deviation from the original EVI data in all three vegetation types in the Qinling Mountains.
There was an increasing trend of vegetation coverage in most regions (about 89.93% of the monitored areas) from 2001 to 2016, and only 11.07% of the study region showed decreased trends. The coefficients of vegetation variation were between 0.02 and 0.6, and the forests that were less affected by human activity showed the highest vegetation stability. The mixed farmland-forest vegetation in the shallow hills of the northern Qinling Mountains and the surrounding areas of mountainous residential communities obviously exhibited alterations with the extraordinary impact of human activity.
The average variation of phenological parameters in the Qinling Mountains indicated that, as the altitude decreased (from the alpine areas to the farmland districts), the SOG was gradually delayed, the EOG was gradually advanced, and the LOG was gradually shortened. These trends were consistent with the local topographical and thermal conditions. Within the Qinling Mountains, the shallow mountainous areas and farmland area of the Funiu Mountains had relatively early SOG dates, which began in the first 10 days of March. By contrast, the SOG of the alpine boreal zones started between the 120th and 135th days of the year. The EOG dates typically started between the 270th and 310th days of the year, which was relatively early for the alpine boreal zones but relatively late in the shallow mountainous areas (typically after the 300th day of the year). The LOG typically lasted between 150 and 200 days∕year and was relatively long in the long-altitude areas (over 180 days) but was relatively short in the high-altitude forest areas (between 150 and 170 days).
Time series variation of phenological from 2001 to 2016 showed that the SOG in the Qinling Mountains was earlier, which was more obvious in high-altitude areas (about 1 to 2 days∕year), and, for altitude below 500 m, there was a little delay, mostly in the southern and northern foothills and the eastern Funiu Mountains. The EOG showed a delaying trend, which was clear (postponed by 2 to 3 days∕year) in the northern slope of Qinling and the eastern mid-tolow-altitude areas. The LOG overall showed a trend of extended duration, with 87.5% of the forest areas showing a change in the duration of the growth season by −1 to 2 days∕year.
The correlation analysis between the climate parameters (i.e., the active accumulated temperature, precipitation, and sunshine duration) and phenological indices (the SOG, EOG, and LOG) in 55 stations across the Qinling Mountains during the 2001 to 2016 period showed that, for the three climate parameters, the active accumulated temperature had the highest correlation with the phenological indices of the SOG, EOG, and LOG, while, for the precipitation and sunshine, duration was relatively low. Global warming was a key factor influencing the phytophenological changes in the Qinling Mountains, as the rising temperature advanced the SOG, delayed the EOG, and prolonged the LOG. This effect was more obvious in the southern and northern slopes of the Qinling Mountains with an altitude of 1000 to 2000 m. However, for the alpine areas, the effect of global warming to the phytophenological change was more complex, which might be related to the low increase of active accumulated temperature and the extraordinary sensitivity of Chinese larches (Larix chinensis) and Farges' firs (Abies fargesii). Stations in the Funiu Mountains in the east part of the Qinling Mountains displayed a poor correlation with the active accumulated temperature, indicating that this region was not significantly impacted by global warming.