Geographic and Climatic Attributions of Autumn Land Surface Phenology Spatial Patterns in the Temperate Deciduous Broadleaf Forest of China

Autumn vegetation phenology plays a critical role in identifying the end of the growing season and its response to climate change. Using the six vegetation indices retrieved from moderate resolution imaging spectroradiometer data, we extracted an end date of the growing season (EOS) in the temperate deciduous broadleaf forest (TDBF) area of China. Then, we validated EOS with the ground-observed leaf fall date (LF) of dominant tree species at 27 sites and selected the best vegetation index. Moreover, we analyzed the spatial pattern of EOS based on the best vegetation index and its dependency on geo-location indicators and seasonal temperature/precipitation. Results show that the plant senescence reflectance index-based EOS agrees most closely with LF. Multi-year averaged EOS display latitudinal, longitudinal and altitudinal gradients. The altitudinal sensitivity of EOS became weaker from 2000 to 2012. Temperature-based spatial phenology modeling indicated that a 1 K spatial shift in seasonal mean temperature can cause a spatial shift of 2.4–3.6 days in EOS. The models explain between 54% and 73% of the variance in the EOS timing. However, the influence of seasonal precipitation on spatial variations of EOS was much weaker. Thus, spatial temperature variation controls the spatial patterns of EOS in TDBF of China, and future temperature increase might lead to more uniform autumn phenology across elevations.


Introduction
Tracking the vegetation growing season is crucial for assessing an ecosystem response to climate change [1] and its impacts on carbon, water and energy exchange between vegetation and the atmosphere [2][3][4][5]. Some studies show that autumn phenology may have a greater control than spring phenology over the growing season length [6] and net ecosystem production (NEP) [7]. However, current research is focused mainly on understanding the relationship between autumn phenology and climate drivers [8,9], which have shown the importance of autumn temperatures in determining the timing of leaf senescence [10,11]. On the other hand, the spatial variation of autumn phenology and its climatic controls have received far less attention. Earlier attempts for identifying environmental controls of phenological spatial patterns can be traced to Hopkins' work. He proposed a "Bioclimatic Law" to estimate the offset in onset of spring as a function of latitude, longitude and elevation in the eastern U.S. [12]. Subsequently, multiple regression equations among average phenological dates and geo-location factors were constructed to estimate phenological gradients, along with latitude, longitude and elevation in different regions [13][14][15][16]. As geo-location factors are not climatic elements, they could not explain the essential environmental causes of spatial difference of phenological occurrence dates or examine the spatial responses of phenological occurrence dates to climatic difference. To reveal climate drivers of phenological spatial patterns, Chen and Xu [5] analyzed relationships between the leaf fall end date of Ulmus pumila trees and the preseason temperature across geographic locations in China's temperate zone, and found that the spatial differences of preseason temperature control the spatial differences of Ulmus pumila leaf fall date. Whether this kind of spatial association based on an individual species is applicable to the ecosystem level is still unclear. Thus, a further examination of the spatial pattern of autumn phenology and its attributions at broad spatial scales over continuous vegetated land areas, such as through the use of remote sensing, is needed. Since spatial patterns of phenology reflect an adaptation of plant life cycle stages to different climates [17], modeling the spatial pattern and its interannual variation in autumn phenology will facilitate a better understanding of phenological responses to recent climate change.
Remotely sensed vegetation indices (VIs) have been commonly used to extract the start (SOS) and end (EOS) of the growing season [18][19][20]. However, few studies have assessed the reliability of different vegetation indices in extracting EOS. For deciduous broadleaf forests, some comparison studies indicated that the Normalized Difference Vegetation Index (NDVI) does not perform better than other VIs [21][22][23][24]. Meanwhile, a new vegetation index, the Plant Senescence Reflectance Index (PSRI), was proposed to determine the stage of leaf senescence and fruit ripening, as it is sensitive to carotenoid retention or accumulation in senescing leaves and ripening fruit [25]. The PSRI value increases during leaf senescence because of the increase in the ratio of carotenoids to chlorophylls. Spectral measurements in northern European boreal forests and moorland throughout the growing season showed that PSRI has similar sensitivity to vegetation dynamics with NDVI in spring, but higher sensitivity than NDVI in autumn [26,27]. Recently, Ren et al. firstly determined PSRI-derived SOS and EOS using moderate resolution imaging spectroradiometer data in the Inner Mongolian grassland, and found a similar performance of PSRI with NDVI [28]. More applications of PSRI should be made to evaluate its ability for identifying EOS in different vegetation types.
In this study, we first compared the spatial pattern of EOS dates retrieved from different VIs with a spatial pattern of ground-observed multispecies mean leaf fall dates at station/pixel scales, and then selected the best VI based on its performance in reflecting the ground reference. Then, we analyzed the relationships between EOS date and geo-location factors (such as latitude, longitude and altitude). Moreover, we constructed spatial phenology models to identify climatic controls of EOS spatial patterns in a multi-year average and for individual years.

Study Area
As the reliability of land surface phenological detection depends highly on selected vegetation types and their seasonality [29][30][31][32], our study focused on the temperate deciduous broadleaf forest (TDBF) with a distinct seasonal nature in northeastern China [33]. The existing TDBF is distributed mainly in mountainous areas, including the middle part of the Da Hinggan Ling Mountains, the Xiao Hinggan Ling Mountains, the Changbai Shan Mountains, the Yan Shan Mountains, the Taihang Shan Mountains and the Qin Ling Mountains from north to south with elevations ranging mainly from 200 to 2000 meters above sea level ( Figure 1). The vegetation formations and sub-formations in the TDBF mainly contain Acer, Ulmus, Quercus, Robinia, Populus, Betula and Salix forest or woodland. The climate type in the TDBF is a temperate monsoonal climate with a hot and rainy summer and a cold and dry winter. The annual mean temperature increases from −8.9 • C in the north to 14.9 • C in the south, while the annual mean total precipitation ranges from 234.7 mm in the northern Yan Shan Mountains to 1456.5 mm in the Qin Ling Mountains. Because of the remarkable spatial heterogeneity of thermal-moisture conditions across the study areas, the plant phenology displays significant spatial variations, which makes the TDBF a suitable biogeographical region for exploring spatial patterns of vegetation phenology and its climate drivers.

Phenological and Meteorological Data
Ground-observed phenological data were acquired from the China Meteorological Administration (CMA). The CMA phenological observation network is the largest phenological observation system in China, and came into operation in 1980. At present, 446 stations regularly undertake phenological observations across the country, in which 27 stations are located within the TDBF area ( Figure 1). According to the availability of phenological data and the representativeness of species, we selected the end of leaf fall stage of the following widely distributed species to validate satellite-retrieved leaf fall stage: Populus berolinensis, Populus simonii, Populus X canadensis, Populus tomentosa, Robinia pseudoacacia, Sophora japonica, Salix matsudana, Salix babylonica, and Ulmus pumila. The geographical coordinates and elevation of phenological stations and selected plant species at each station are shown in Table S1.
Moderate resolution imaging spectroradiometer (MODIS) data have been available since 2000, however, phenological data beyond 2012 were not yet available at the time of this study. Thus, we utilized a ground-observed 13 year continuous record of the leaf fall end date of each species and MODIS data from 2000 to 2012. For each species, at least three individual plants were selected as fixed observation objects at a given station. The leaf fall end date was identified as the day when more than 90% of the leaves fell off in more than half of the observed trees [34]. We calculated the average leaf fall end date of all indicator species at each station to represent roughly the leaf fall end date at plant community scales, which should be more comparable with the remotely sensed EOS than the leaf fall end date of a single species.
Daily temperature and precipitation data from 2000 to 2012 were acquired from the China Meteorological Forcing dataset (CMFD). This dataset is a fusion of data from multiple sources, including the Global Land Data Assimilation System (GLDAS) forcing data, Princeton forcing data, Global Energy and Water Cycle Experiment-Surface Radiation Budget (GEWEX-SRB) downward shortwave radiation, the Tropical Rainfall Measuring Mission (TRMM) satellite precipitation analysis data (3B42) and GLDAS precipitation and CMA station data [35]. This data product is available at a spatial resolution of 0.1 degrees and a temporal resolution of 3 h with seven basic climatic factors (e.g., near-surface air temperature, precipitation, maximum wind speed, etc.). In this study, daily mean air temperature was calculated by the average of the eight temperature values in each day, while daily precipitation was calculated by an arithmetic sum of the eight values of accumulated precipitation every 3 h within a day.

MODIS Data and Processing
We used Moderate Resolution Imaging Spectroradiometer (MODIS) Surface Reflectance 8-Day L3 Global 500 m product (MOD09A1, V5) distributed by the Land Processes Distributed Active Archive Center (LP DAAC). Surface reflectance was computed from seven MODIS Level 1B bands (band 1 to 7 centered at 648 nm, 858 nm, 470 nm, 555 nm, 1240 nm, 1640 nm and 2130 nm, respectively) after removal of the atmospheric effect [36], in which the surface reflectance data from bands 1-5 were used to calculate the six types of vegetation indices, including the Normalized Difference Vegetation Index (NDVI), Land Surface Water Index (LSWI), Enhanced Vegetation Index (EVI), Wide Dynamic Range Vegetation Index (WDRVI), Green-Red Ratio Index (GRVI) and Plant Senescence Reflectance Index (PSRI). Detailed descriptions of these indices are shown in Table 1. Table 1. Descriptions of vegetation indices used in this study. [25] MOD09A1 provides a quality assurance (QA) layer (indicating the abnormal pixels contaminated by snow, ice or clouds on the data) and a day of the year (DOY) layer (marking the dates of data acquisition for respective pixels). Utilizing the QA and DOY flags, we identified the abnormal VI pixel values that were affected by snow, ice or clouds. Then, we replaced the abnormal VI values during the start and end stages of the year using the nearest normal VI values and other abnormal VI values with linearly interpolated VI values of the two nearest normal VI values at both sides. Moreover, we identified the date (D f ) when the minimum VI (maximum PSRI) occurred in the first half of the reconstructed VI sequence and the date (D s ) when the minimum VI (maximum PSRI) occurred in the second half of the reconstructed VI sequence. All VI values from 1 January to D f and from D s to 31 December were further replaced by the average of VI values corresponding to D f and D s . Here, the minimum VI (maximum PSRI) value represents the initial background VI value [42]. Finally, we smoothed the VI time series with a five point moving average filter ( Figure 2). As the double logistic method is demonstrated to be particularly useful for describing annual VI time series for phenological monitoring [31,43], we fitted the smoothed VI time series with the double logistic function [28,32]. The corresponding formula is as follows:

Equations 2 References
where VI(t) is the VI value of a given pixel at time t; VI max and VI min are the maximum and minimum VI values in each annual time series; I and D denote the maximum falling and rising slopes at the inflection points of each VI curve and S and E denote the corresponding dates of I and D. n equals to 1 for PSRI, while n equals to 0 for other VIs. The EOS dates (day of year) were then identified using Equation (2) [32] at each pixel and for each year, as they are the dates of the last local minima curvature change rate of the fitted curves. All data processing was conducted with MATLAB software (Matlab R2014a, MathWorks Inc. Natick, MA, USA).
To compare the ability of different VIs to capture the spatial pattern of ground autumn phenology, we calculated the root mean square error [RMSE; Equation (3)] and Pearson correlation coefficient (r) between the spatial series of the ground-observed LF date at 27 stations and the VI-derived EOS date at the corresponding pixels. The comparison was made both for individual years as well as for a multi-year average.
where EOS i and LF i are the remotely-sensed end date of the growing season and the ground-observed leaf fall end date at station i, respectively. n is the number of stations. The best performing VI was selected by the lowest RMSE value and highest correlation (r) for multi-year average, and by the largest number of years with the lowest RMSE value and highest correlation for individual years. Then, we used the best performing VI to identify the EOS for each pixel over the entire study area.
To ensure the accuracy of identified EOS dates, we eliminated pixels with the following attributes: (1) The EOS date did not appear in the current year; (2) The annual maximum VI (minimum PSRI) was not detected from April to October. The satellite-derived EOS dates were resampled to a spatial resolution of 0.1 degrees using a moving window averaging method to match the spatial resolution of climatic data. Because of the removal of abnormal pixels, the total number of pixels in TDBF were different in different years. To ensure the same geographic range among the different years, we extracted the common pixels from 2000 to 2012, and obtained 4091 pixels for further study.

Spatial Phenology Models Based on Geo-Location Indicators and Climatic Factors
To examine the attribution of yearly and multi-year average EOS spatial patterns, multiple regression models were created among the EOS date and geo-location indicators (latitude, longitude and elevation) across the 4091 common pixels [16]. The elevations at each pixel were determined by Digital Elevation Model (DEM) data extracted from the Global 30 Arc-Second Elevation (GTOPO30) dataset, in which DEM data were resampled from the spatial resolution of 30 to 0.1 • using a cubic convolution method in ArcGIS 10.3 software. The multiple regression model is shown as follows: where EOS (x, y, z) is the EOS date (in day of year) at latitude x, longitude y and altitude z; a 0 is a constant; a x , a y and a z are regression coefficients denoting latitudinal sensitivity, longitudinal sensitivity and altitudinal sensitivity, respectively. Partial correlation analysis was also used to detect dependency of the EOS spatial pattern on each of the geo-location indicators (x, y and z), respectively. We used mean air temperature and accumulated precipitation spatial series during the optimum length period to correlate with the EOS spatial series. The basic hypothesis was that the station-to-station variation of a phenological event occurrence date over an area is mainly influenced by station-to-station variation of climatic factors within a particular length period (LP) of days during and before its occurrence over the area [5]. The LP in this study was determined as the time periods during which the pixel-to-pixel variation in daily mean temperature or accumulated precipitation affects the pixel-to-pixel variation of EOS most remarkably in a year. The computation process was implemented through the following steps: First, we defined the number of days between the earliest and latest dates of the EOS spatial series across all the pixels in a given year as the basic length period (bLP). Then, we calculated backwards (to an earlier time direction) the mean air temperature and accumulated precipitation spatial series across all the pixels in the year during bLP + mLP. Here, mLP was a moving length period prior to the earliest date of EOS spatial series, which moved with a step length of 1 day, namely, during the bLP+1 day, bLP+2 day, bLP+3 day, etc. The maximum mLP was set to 60 days, as the mean temperature and accumulated precipitation prior to this period were usually assumed to have limited effects on spatial variation of occurrence dates of a phenological event. Thus, the LP is defined as follows: Furthermore, we calculated the correlation coefficients between the EOS spatial series and mean temperature/accumulated precipitation spatial series during different testing LPs across all the pixels in a specific year. Finally, we obtained the optimum LP based on the highest correlation between the EOS spatial series and the corresponding mean temperature/accumulated precipitation spatial series. The temperature or precipitation within the optimum LP across all the pixels in the year was thought of as the major influence factor of EOS spatial variation. Hence, the optimal temperature/precipitation-based spatial phenology model for each year can be created as: where X is the mean temperature or accumulated precipitation during the optimum LP; a and b are the intercept and slope of the equation, respectively. To evaluate the relative importance of temperature and precipitation on controlling EOS spatial patterns, we further executed a partial correlation analysis between EOS spatial series and temperature/precipitation spatial series during the optimum LP. The significance levels of partial correlation were evaluated based on two-tailed significance tests.

Ground Validation of EOS Retrieved from Different Vegetation Indices
We used the spatial RMSE to measure the station-to-station differences between satellite-derived EOS dates and ground-observed LF dates in each year. Results show that the medians and ranges of RMSE for EOSs retrieved from all six VIs were between 14 days and 23.2 days, and between 5.6 days and 11 days, respectively, in which either the median or range of RMSE for PSRI-derived EOS was the smallest (Figure 3). To evaluate the spatial variation consistency between satellite-derived EOS dates and ground-observed LF dates, we conducted a correlation analysis between EOS spatial series and LF spatial series across all stations/pixels on a multi-year mean and for each year. On average, EOS dates retrieved from all six VIs were significantly (p < 0.001) and positively correlated with LF dates, although the satellite estimates were systematically earlier ( Figure 4). In addition, relatively high correlations (r ≥ 0.7) were found between all EOS estimates and LF. The correlation coefficient between the PSRI-derived EOS and LF was also the highest (0.76; Figure 4f), followed by that of LSWI (0.74; Figure 4b). For each year, the correlation analysis showed a consistently positive relationship. The number of years with a significantly positive correlation (p < 0.01) ranged from 7 years (GRVI) to 12 years (PSRI) out of the 13 years in total ( Table 2). The PSRI-based EOS performed better than other VIs in all cases, in terms of yielding relatively higher correlation coefficients.
In summary, an outstanding performance of the PSRI-derived EOS was demonstrated by its consistence with ground-observed LF. In contrast, the performance of the GRVI-derived EOS appeared to be the poorest. Performances of EOS retrieved from the other four VIs fall in between these two. Among the four VIs, the NDVI-derived EOS was the most accurate (with the smallest RMSE; Figure 4a) and the LSWI-derived EOS could more consistently reflect the spatial variation of ground-observed LF (with largest correlation coefficient; Figure 4b). Finally, we selected the PSRI-derived EOS for the further analysis.

Spatial Patterns of EOS and Their Relation to Geo-Location Factors
Based on PSRI data from 2000 to 2012 at each pixel, we extracted the EOS over the entire TDBF area. The multi-year mean EOS ranged from 272 DOY to 346 DOY, and occurred earlier at higher latitudes than at lower latitudes. Specifically, the earliest EOS was detected at the northern parts of the study area, namely in the Da Hinggan Ling Mountains and Xiao Hinggan Ling Mountains, whereas the latest EOS was found in the south, such as in the Qin Ling Mountains and Taihang Shan Mountains ( Figure 5). This multiple regression model explained 80% of the spatial variance of multi-year mean EOS, with a simulation error of 7.1 days. Partial correlation analyses showed significantly (p < 0.01) negative relationships of multi-year mean EOS with latitude, longitude and altitude. Overall, the latitudinal effect on EOS was stronger (r x = −0.81) than the longitudinal effect (r y = −0.46) and altitudinal effect (r z = −0.5). Similar relationships were also detected in all years when yearly EOS was fitted (Table 3). Table 3. Multiple regression equation parameters, evaluation indicators and partial correlation coefficients among growing season end date, latitude, longitude and altitude. a 0 is a constant; a x , a y and a z are regression coefficients denoting latitudinal sensitivity, longitudinal sensitivity and altitudinal sensitivity, respectively; EOS is the mean growing season end date across the study area; R 2 is the goodness of fit of the multiple regression model; r x , r y and r z are partial correlation coefficients between the growing season end date and latitude/longitude/altitude after controlling for other two geo-location indicators, respectively. Moreover, both the absolute values of yearly regression equation slopes along the altitude (a z ) and the yearly partial correlation coefficients between EOS and the altitude (r z ) indicated significant decrease trends from 2000 to 2012 ( Figure 6). That is, either the altitudinal sensitivity of EOS or the effect of altitude on EOS became weaker and weaker from 2000 to 2012.

Spatial Relationship between EOS and Climatic Factors
To investigate climatic controls of the EOS spatial pattern, we used both temperature-based and precipitation-based spatial phenology models to simulate and predict the EOS spatial series for each year. The temperature-based spatial phenology models showed a significant (p < 0.001) positive relationship between EOS and mean temperature during the optimum LP in each year-namely, the higher the autumn temperature at a pixel, the later the EOS date at the pixel, and vice versa. These models explained 53% to 73% of the EOS spatial variance, and yielded RMSEs of 8.7 (in 2012) to 11.3 days (in 2001) between the observed and predicted dates. Slopes of the linear models indicate that the fitted EOS varied with LP mean temperatures spatially at a rate of 2.4 days/ • C (in 2010) to 3.6 days/ • C (in 2003 and 2011; Figure 7). Similarly, the precipitation-based spatial phenology models also displayed a significant (p < 0.001) positive relationship between EOS and accumulated precipitation during the optimum LP in each year. This indicated that the more autumn accumulated precipitation at a pixel, the later EOS at the pixel, and vice versa. The precipitation-based spatial phenology models explained 10% to 46% (p < 0.001) of the EOS spatial variance, which were much smaller than the temperature-based spatial phenology models. In addition, the precipitation-based models yielded larger RMSEs than temperature-based spatial phenology models, which ranged from 12.7 (in 2010) to 16.7 days (in 2003; Figure S1). Yearly partial correlation coefficients between EOS and temperature during the optimum LP were all significant (p < 0.001), and ranged from 0.62 (in 2001) to 0.80 (in 2006 and 2012), while partial correlation coefficients between EOS and precipitation during the optimum LP were significant (p < 0.001) in 9 of 13 years, and ranged from −0.12 (in 2002) to 0.48 (in 2011). It should be noted that absolute values of partial correlation coefficients between EOS and accumulated precipitation were much smaller than those of partial correlation coefficients between EOS and the mean temperature (Table 4). Table 4. Yearly partial correlation coefficients between EOS and the mean temperature/accumulated precipitation during the optimum length periods. r t denotes the partial correlation coefficient between the growing season end date and the mean temperature after eliminating the influence of accumulated precipitation; r p denotes the partial correlation coefficient between the growing season end date and accumulated precipitation after eliminating the influence of the mean temperature.

Discussion
This study assessed the ability of different remotely sensed vegetation indices in monitoring autumn phenology. Previous studies showed that different VIs are more sensitive to different aspects of forest canopies. For instance, NDVI, EVI and WDRVI are more sensitive to canopy structure and color [21,39,40], but GRVI to greenness [41], LSWI to canopy moisture content [44] and PSRI to carotenoid and mesophyll cell structure [45]. As each type of vegetation indices shows obvious seasonal characteristics that reflect local vegetation dynamics, all indices can be used for detecting autumn phenology in deciduous forests. In this study, we found that the more widely used NDVI and EVI did not perform better than other VIs in monitoring deciduous forest autumn phenology [21,23,28], and that PSRI best represented the spatial distribution of the ground-observed autumn phenology. This may be due to the fact that PSRI can capture not only the change of mesophyll cell structure but also the change of the relative proportion of carotenoid and chlorophyll in leaves, giving it greater sensitivity to leaf senescence than other VIs [45].
The multiple regression model of multi-year mean remotely sensed EOS dates showed similar latitudinal and altitudinal gradients with models of multi-year mean estimated leaf coloring dates for all three representative deciduous broadleaf species in China [46]. That is, autumn phenology occurs earlier as latitude and altitude increase. Moreover, the longitudinal gradient of EOS was consistent with that of estimated leaf coloring dates for two of three species (except Fraxinus chinensis), namely, autumn phenology occurs earlier as longitude increases. Although the two modeling regions did not overlap completely, consistent responses of both EOS date and leaf coloration date on geo-location indicators evidence that remotely sensed autumn phenology can reflect overall spatial variations of autumn phenology on the ground. The significantly decreasing trend of altitudinal sensitivities (in advancing days per 100 m upward) for EOS from 2000 to 2012 might be explained by larger delay rates of EOS dates at high elevations than at low elevations, considering global change. These kinds of uneven delay rates of autumn phenology at different elevations have been also detected in Slovenia for in situ European beech leaf coloration dates [47].
In this study, we employed a spatial phenology model that was originally developed using ground observations [5] to simulate multi-year mean and yearly spatial patterns of remotely sensed autumn phenology. Unlike the previously revealed spatial relationship between discrete ground autumn phenological date and mean temperature during the optimum LP [5], the remotely sensed EOS dates allowed for the characterizing of similar spatial relationships with mean temperature (delayed autumn phenology with an increase of the optimum LP mean temperature) in a continuous manner. The consistency of temperature control on the spatial pattern of remotely sensed and ground-observed autumn phenology suggests that the spatial phenology model is applicable to both ground-observed and remotely-sensed phenological simulations. Sensitivities of remotely-sensed autumn phenology to mean temperature during optimum LPs ranged from 2.4 days/ • C to 3.6 days/ • C in this study (Figure 7), which is also similar to previous estimates (2.2 days/ • C to 3.2 days/ • C) using leaf fall dates of ground trees in northern China [5].
As anticipated, our study showed that temperature-based models can explain more spatial variations of EOS than precipitation-based models in the temperate deciduous forests. Our results further confirm that the spatial variation of remotely sensed EOS in temperate eastern China are controlled mainly by spatial variations of temperature [48], rather than precipitation. Moreover, the partial correlation coefficients between spatial EOS and spatial accumulated precipitation during the optimum LP indicate inconsistent relations (positive or negative) among years (Table 4). Thus, the underlying reasons for precipitation controls on the spatial pattern of EOS need to be further studied.

Conclusions
This study assessed and compared the reliability of six representative vegetation indices in monitoring autumn phenology, and analyzed spatial characteristics of growing season end date and its climatic controls in the temperate deciduous broadleaf forest area of northeastern China. The main conclusions are as follows: 1.
PSRI-derived EOS can more accurately capture spatial variations of ground-observed leaf fall dates of dominant trees on the multi-year average and in each year than the other five vegetation indices. The spatial variation of PSRI-derived EOS dates correlate significantly and positively with the spatial variation of ground based leaf fall dates in each year.

2.
Multi-year mean PSRI-derived EOS represent a latitudinal sensitivity of −2.98 days per degree northward, a longitudinal sensitivity of −1.03 days per degree eastward and an altitudinal sensitivity of −1.5 days per 100 m upward. The altitudinal sensitivity of EOS became weaker and weaker from 2000 to 2012, which might be attributable to larger delaying rates of EOS dates at high elevations than at low elevations under recent climate warming.

3.
Temperature-based phenological models can explain more spatial variations of EOS with less simulation errors than the precipitation-based models in the study area. Thus, the spatial variation of remotely sensed EOS are controlled mainly by spatial variations of temperature, rather than precipitation.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/11/13/1546/s1. Figure S1: Spatial correlation and regression analyses between PSRI-derived growing season end date and accumulated precipitation during the optimum length period (LP) in each year across the study area; Table S1: Geographical coordinates and elevation of phenological stations and selected plant species at each station.