Spring and autumn phenology across the Tibetan Plateau inferred from normalized difference vegetation index and solar-induced chlorophyll fluorescence

ABSTRACT Plant phenology is a key parameter for accurately modeling ecosystem dynamics. Limited by scarce ground observations and benefiting from the rapid growth of satellite-based Earth observations, satellite data have been widely used for broad-scale phenology studies. Commonly used reflectance vegetation indices represent the emergence and senescence of photosynthetic structures (leaves), but not necessarily that of photosynthetic activities. Leveraging data of the recently emerging solar-induced chlorophyll fluorescence (SIF) that is directly related to photosynthesis, and the traditional MODIS Normalized Difference Vegetation Index (NDVI), we investigated the similarities and differences on the start and end of the growing season (SOS and EOS, respectively) of the Tibetan Plateau. We found similar spatiotemporal patterns in SIF-based SOS (SOSSIF) and NDVI-based SOS (SOSNDVI). These spatial patterns were mainly driven by temperature in the east and by precipitation in the west. Yet the two satellite products produced different spatial patterns in EOS, likely due to their different climate dependencies. Our work demonstrates the value of big Earth data for discovering broad-scale spatiotemporal patterns, especially on regions with scarce field data. This study provides insights into extending the definition of phenology and fosters a deeper understanding of ecosystem dynamics from big data.


Introduction
Spatial heterogeneity is among the most outstanding features of physical geography, and the understanding of this heterogeneity requires extensive data covering a broad range of the Earth's surface (Guo, 2017;Nativi et al., 2015). However, the acquisition of such extensive data had been difficult until recent decades, when space-based Earth observations emerged and has since grown quickly (Guo et al., 2020). Indeed, the recent rapid growth of Earth observation data has pushed physical geography, as well as many other geoscience subjects, into a "Big Data Era" (Guo, 2017;Guo et al., 2017Guo et al., , 2020Hampton et al., 2013;LaDeau, Han, Rosi-Marshall, & Weathers, 2017). Nonetheless, challenges remain on how to better inform geographical patterns, processes, and their underlying mechanisms with big data.
Space-based Earth observations are particularly useful for regional studies of the Tibetan Plateau (TP), where the high elevation and harsh environment make access to ground-based monitoring extremely difficult (Piao et al., 2019b;Yao et al., 2019). Rising over 4000 m above the sea level on average, the TP has experienced one of the most rapid rates of land warming in recent decades (Chen et al., 2015;Hansen, Ruedy, Sato, & Lo, 2010;Yao et al., 2019). This fast warming, together with the high climate sensitivity of the TP's alpine ecosystems, drives rapid changes in its ecosystem structures and functions (Klein et al., 2014;Piao et al., 2019a). Plant phenology has long been used as a valuable ecosystem indicator sensitive to climate change (Piao et al., 2019a;Richardson et al., 2013). Investigations based on remote sensing data suggest a significant advancement trend for the start of the growing season (SOS) and a delay for the end of the growing season (EOS) on the TP since the 1980s Cong, Shen, & Piao, 2017;Piao et al., 2011). This result is consistent with observations over many other northern ecosystems (Jeong, Ho, Gim, & Brown, 2011). Whether or not this advancing SOS trend persists into the 21 st century remains an active topic of debate in the literature (Chen, Zhu, Wu, Wang, & Peng, 2011;Shen et al., 2013;Wang et al., 2013;Yi & Zhou, 2011;Yu & Xu, 2010;Zhang, Zhang, Dong, & Xiao, 2013). Inconsistent and even contrary findings have been reported, with such inconsistence potentially arising from differences among different studies in the selection of satellite products, data processing, and phenology retrieving algorithms (Shen et al., 2014;Wang et al., 2015;Zhang et al., 2013). Therefore, these debates bring the urgent need to verify the validity of phenology extracted from big data, by using new datasets and integrated comparisons of different algorithms.
The recent emergence of satellite-based solar-induced chlorophyll fluorescence (SIF) brings a new set of tools for monitoring vegetation photosynthetic activities (Chen et al., 2021;Frankenberg, Butz, & Toon, 2011;Joiner, Yoshida, Vasilkov, Yoshida, & Corp, 2011). As SIF is emitted from the chlorophyll of photosynthesis-active plants, it is arguably less impacted by cloud, snow cover, or land degradation that may contaminate greennessbased vegetation indices (Zhang, Joiner, Alemohammad, Zhou, & Gentine, 2018a). Since plant phenophase dates such as SOS and EOS are proxies for the start and end of growing season over a year. SIF products as the measurements of annual photosynthetic activities may also have the potential to track plant phenology changes. However, SIF-based phenology research is still in its infancy, and prudent cautions are needed in its development and application, particularly with regard to the appropriate algorithms to retrieve phenology dates from SIF data time series (Zhang, Commane, Zhou, Williams, & Gentine, 2020). Furthermore, it is unknown whether the SIF phenology shows similar responses to environmental variations as the NDVI phenology.
In this experiment, we used the contiguous solar-induced chlorophyll fluorescence (CSIF), a gridded SIF product fusing the Orbiting Carbon Observatory-2 (OCO-2), SIF data, and Moderate Resolution Imaging Spectroradiometer (MODIS) MCD43C4 C6 reflectance data (Zhang et al., 2018a), to investigate spatiotemporal changes in vegetation SOS and EOS on the TP from 2001 to 2018. The primary goal of our research is twofold: (1) to examine the suitability of using SIF products for large-scale vegetation phenology studies on the TP; and (2) to resolve the debate as to whether the phenology changes detected for the 1980s and 1990s has continued into the recent two decades. Therefore, we focused on the 2001-2018 period, because phenology changes in the TP since the 2000s are most controversial in current literature (Klein et al., 2014). We also chose to focus on this time period because SIF products are only available for that period. We used five common algorithms, which have been widely used in previous studies, to derive SOS and EOS from the time series data of CSIF (see Methods). Furthermore, as a comparison, we also repeated the analyses with MODIS NDVI data. For convenience, we hereafter denote SOS and EOS based on the SIF product as SOS SIF and EOS SIF , respectively; and SOS and EOS based on the NDVI product as SOS NDVI and EOS NDVI , respectively.

Data
The SIF time series dataset used in this study was downloaded from the CSIF project (Zhang et al., 2018a), which was calculated by a machine learning method using data of OCO-2 SIF and MCD43C4 C6 reflectance. This CSIF dataset covers the global land surface with a 4-day temporal resolution and 0.05° spatial resolution, and is available from 2001 to 2018 (https://doi.org/10.17605/OSF.IO/8XQY6) at the time when we accessed it. The data of NDVI was obtained from the product of MODIS MOD13C1 (Collection 6), which has a 16-day temporal resolution and 0.05° spatial resolution, and is available from 2000 to 2020 (https://lpdaac.usgs.gov/products/mod13c1v006/). The vegetation classification map of the TP was adopted from a digitalized 1:1000,000 vegetation map of China (Editorial Board of Vegetation Map of China, 2001, (http://westdc.westgis.ac.cn), which was used to remove pixels that lacked seasonal dynamics.
Climate data, including temperature and precipitation, were downloaded from the China Meteorological Forcing Dataset (CMFD, He et al., 2020;Yang & Huang, 2019), which has a 3-h temporal resolution and 0.1° spatial resolution. This data set is available from January 1979 to December 2018 (http://data.tpdc.ac.cn/en/data/8028b944-daaa-4511-8769-965612652c49/). The time period we adopted was from 2001 to 2018, according to the time range overlap of all datasets.

Preprocessing of datasets
We focused on natural vegetation with clear seasonal rhythms, thus excluding subtropical and cultivated areas following the vegetation classification map. We also removed areas with sparse vegetation cover or bare land, defined as pixels that have a multi-year average NDVI value under 0.1 (Zhou et al., 2001). To simplify the computation and ensure the comparability, we converted the spatial resolution of all datasets to 0.05° using the nearest neighbor interpolation (Gao et al., 2020). Furthermore, to avoid spurious start and end dates of the growing season, we imposed a thermal growing season, which was defined as periods with temperatures greater than 0°C in early spring and late autumn, to constrain these phenological dates . In addition, the climatic factors with a temporal resolution of 3 h were resampled to daily mean temperature, and daily sum precipitation.

Algorithms used for phenological data retrieval
We employed five commonly used methods to derive the spring and autumn phenophases from both the SIF and NDVI time series data: Hants-Mx, Polyf-Mx, Spln-MP, Dlog and Plog (Cong et al., 2012;Liu et al., 2016a;Zhang et al., 2020). The growing season length (GS) was the difference between EOS and SOS.
(1) Hants-Mx: The harmonic analysis of time series (Hants) function (Cong et al., 2012;Liu et al., 2016b) was used to fit SIF and NDVI data with different frequencies (2 or 4), depending on their thermal growing season length (Eqn 1). Any negative values were reassigned as zero. The function was used for fitting raw data of every pixel with an interval of 16 or 4 days, depending on the different datasets each year. The function removed instances of cloud contamination and outliers, and then fitted the continuous smoothed time series data. The maximum and minimum rate of change of fitted growing season data corresponded to SOS and EOS, respectively (Eqn 2).
where x is the day of year, a is the fitted coefficient. n and ω i are set to 1 and 2/4 π based on the thermal-growing season, if the thermal growing season length is more than 90, the ω i is 2 π, otherwise, it is set to 4 π. φ i is the maximum SIF value. The maximum and minimum of SIF ratio are the SOS and EOS, respectively.
(2) Polyf-Mx: To remove the interference of outliers or spikes, we first used the Savitzky-Golay filtering (SG filter) to smooth the raw data (Eqn 3, Chen et al., 2004). The ensemble mean value outside the thermal growing season was set as the background value. We then used sixth-and eighth-order polynomial functions to fit SIF and NDVI data to acquire the smoothed daily data (Eqn 2, Piao, Fang, Zhou, Ciais, & Zhu, 2006;Zhang et al., 2020). This process is similar to Hants-Mx for extracting phenological dates.
where x is the day of year, a is the fitted coefficient. n is set to 6 for NDVI and 8 is set for SIF in our study . The maximum and minimum of SIF ratio are the SOS and EOS, respectively (Eqn 2).
(3) Spln-MP: Similar to Polyf-Mx, we first removed outliers by SG filtering (Chen et al., 2004). Then, we used a cubic smoothing spline function to smooth raw data and acquired a continuous dataset for each year (Eqn 4, Cong et al., 2012;Liu et al., 2016b). We then calculated the date that reached 50% of the difference between the fitted maximum and minimum values on both sides of the curve; these dates were set as the SOS and EOS dates, respectively (Eqn 5).
where x is the day of year, and a is the fitted coefficient. The SOS and EOS occur when the ratios reach 50% of the maximum of SIF.
(4) Dlog: We first used a double logistic function to fit the raw data (Eqn 6). This method produced a good fit for the growing season curve. In contrast to the above-mentioned three methods, function parameters in Dlog were the SOS and EOS dates (Julien & Sobrino, 2009). The parameters were optimized using the Levenberg-Marquardt algorithm (Moré, 1978).
where x is the day of year, a, b, c and d are the fitted coefficients. SIF min and SIF max are the annual minimum and maximum values of SIF, respectively.
(5) Plog: We used a piecewise logistic function to fit the raw data. Similar with Dlog, this algorithm used two logistic functions to fit the first half curve and the second half curve, respectively (Eqn 7, Moré, 1978;Zhang et al., 2003).
Eqn7 where x is the day of year, a 1 , a 2 , b 1 and b 2 are the fitted coefficients. c 1 and c 2 are the background SIF values of non-growing season. SIF max is the annual maximum value of SIF. b 1 and b 2 are the SOS and EOS, respectively.

Analysis
Firstly, we examined the uncertainties of different methods by calculating standard deviations of different algorithms and comparing them with ground observation. Due to the lack of long-term observations, only one ground observation of SOS is acquired from reference Li et al. (2016). With derived SOS and EOS dates from the SIF and the NDVI time series data, we first examined each data set's spatial and temporal patterns on the TP. We then investigated their uncertainties, which we derived from different algorithms by calculating the interalgorithm standard deviations (SD). We also compared the spatiotemporal differences in phenophases extracted from SIF and NDVI. Notably, we used the multi-year average values to show the spatial patterns of SOS and EOS, and defined the temporal change rate (β) as the slope of these variables over time (Eqn 8). We reserved the term "trend" for temporal changes that were significant at the p < 0.05 level.
where time is the years from 2001 to 2018 for each grid i. To further understand potential climatic drivers of the spatiotemporal patterns of derived phenophases, we performed a partial correlation analysis (Eqn 9) between phenological dates and climatic factors over 18 years for each pixel. Because both SOS and EOS are strongly controlled by their pre-season climate (Fu et al., 2015;Gusewell, Furrer, Gehrig, & Pietragalla, 2017), we performed a partial correlation analysis between these phenological dates and "pre-season" climate, with a 10-day interval decrement for each step before the SOS or EOS dates. For example, if the Julian date of SOS was 100, the loop time intervals were 10 (from 90 to 100), 20 (from 80 to 100), and 100 (from 0 to 100) days. The period with the maximum correlation coefficient (either positive or negative) for each climatic variable was defined as the most related preseason of each variable for each phenophase. All statistical analyses were reported at a significance level of p < 0.05 unless otherwise stated.
x 2 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where R x 1 y is the partial correlation coefficient of the factor (x 1 ), R yx 1 , R yx 2 and R x 1 x 2 are the simple correlation coefficients between x 1 and the dependent variable (y), x 2 and y, x 1 and x 2 , respectively. Figure 1 shows the ensemble mean values of SOS SIF and EOS SIF from the five algorithms averaged for the period from 2001 to 2018 across the TP. SOS SIF on the TP ranged from May to July, with a clear spatial pattern that was earlier in the east than in the west ( Figure  1a). EOS SIF mostly ranged from early September to early October ( Figure 1b). However, unlike SOS SIF , EOS SIF did not show a highly distinct spatial pattern, although it seemed to be earlier in the center of the plateau than in the eastern or western parts ( Figure 1b). Overall, SOS SIF on the TP spanned about 2.5 months, while the range of EOS SIF was much narrower, spanning only about 1 month across the entire plateau ( Figure 1a and Figure  1b). Furthermore, the SD of SOS SIF and EOS SIF computed by the five algorithms also increased from east to west, with the highest SD found in the semiarid western area (Figure 1c and Figure 1d). This distinct spatial pattern found in SOS SIF , as well as the less distinct spatial pattern found in EOS SIF , was also confirmed by examining the results from each individual algorithm (Fig. S1). However, marked differences were also found among results from each of these five algorithms (Fig. S1). In particular, SOS SIF derived with the Polyf algorithm was consistently earlier than that derived with other algorithms by 8.0-14.4 days on average over the entire TP (Figure 1e). Similarly, EOS SIF derived with the Polyf algorithm was later than that derived with other algorithms by 4.3-12.5 days (Figure 1f). EOS SIF derived with the Dlog algorithm was also later than that from the Hants, Spln, and Plog by 7.9, 6.6, and 8.1 days on average, respectively (Figure 1f). These inter-algorithm differences were also evident by looking into the distribution histograms of SOS SIF and EOS SIF from the individual algorithms (Figure 1e and Figure 1f).

Temporal changes in SIF-derived SOS and EOS
Next, we analyzed temporal changes in SOS SIF and EOS SIF from 2001 to 2018. As shown in Figure 2, the ensemble mean SOS SIF in most areas (67%) was advanced over the past two decades, with 12% of the pixels showing a significant trend (Figure 2a). Most of the significant advancing trend in SOS SIF was found in the northeastern and central TP ( Figure  2a). However, it was also notable that SOS SIF showed delayed changes in the western semiarid ecosystems, although the change was not statistically significant probably due to the limited number of years investigated (Figure 2a). For EOS SIF , its temporal changes  ranged between −0.8 and 0.8 days yr −1 , being negative (advanced) in the east but positive (postponed) in the west and the south (Figure 2b). However, most of the changes were insignificant (Figure 2b).
These spatial patterns of ensemble mean temporal changes in both SOS SIF and EOS SIF were also reproduced by each individual algorithm (Fig. S2). Importantly, results of the temporal changes from different algorithms were nearly identical for both SOS SIF and EOS SIF , evident by the small inter-algorithm SD in annual SOS SIF (or EOS SIF ) changes (<-0.2 d yr −1 , Figure 2c and Figure 2d) and their almost overlapped frequency distributions of annual change values (Figure 2e and Figure 2f). This high inter-algorithm similarity in temporal SOS SIF (or EOS SIF ) changes was in contrast to the clear inter-algorithm differences for the spatial pattern of multi-year averaged SOS SIF (or EOS SIF ).

Climatic determinants of SIF-derived SOS and EOS
Partial correlation analysis between phenophases (i.e. SOS SIF and EOS SIF ) and preseason temperature or precipitation, suggested that temperature generally had a negative correlation with SOS SIF (nearly 76% of the total grids (significant for 35% of the total grids)). Pixels with significant negative correlations between temperature and SOS SIF were mostly located in the central and eastern TP (Figure 3a). Positive correlations between preseason temperature and SOS SIF accounted for the remaining 34% (12% significant) of the grids, primarily in the southwestern TP (Figure 3a). The effect of temperature on EOS SIF was roughly equal between negative (51%, among which 19% being significant) and positive (49%, among which 16% being significant) (Figure 3c). It is interesting to note that in the eastern TP, temperature exerted a negative impact on both SOS SIF and EOS SIF ; while in the southwest, it turned positive for both SOS SIF and EOS SIF (Figure 3a and Figure 3c). By contrast, precipitation had a negative relationship with SOS SIF (70% (29%)) but a positive relationship with EOS SIF (73% (32%)) in most areas (Figure 3b and Figure 3d).

Spatiotemporal patterns of NDVI-based SOS and EOS and their climatic determinants
As a comparison, we also extracted the two phenophases with NDVI, that is, SOS NDVI and EOS NDVI , using the same five algorithms. Both the spatial distribution and the temporal change of SOS NDVI and EOS NDVI seemed to follow the same patterns as that of SOS SIF and EOS SIF , respectively, with the exception of the spatial pattern of EOS (Figure 4, S1 and S2). Compared to EOS SIF , which did not show a clear spatial pattern, EOS NDVI progressively advanced from east to west (Figure 4d). Similar to SIF-based phenophases, SOS NDVI showed an advancing trend in the east of the plateau, but no significant temporal changes in other regions ( Figure 4g); and EOS NDVI rarely showed any significant changes across the entire region ( Figure 4j). Furthermore, the spatial patterns of the partial correlations between climate and phenophases were similar for SIF-based versus NDVIbased phenophases (comparing Figure 5 to Figure 3, Fig. S5).
Looking into the difference between SOS SIF and SOS NDVI (and between EOS SIF and EOS NDVI ), we found that spring started earlier by SOS NDVI than by SOS SIF , and autumn ended later by SOS NDVI than by SOS SIF, for nearly all the pixels (97% and 99% of the total area, respectively; Figure 6). A difference of up to 30 days for both SOS and EOS was primarily found in the eastern TP (Figure 6a and Figure 6b). As a result, observations based on NDVI led to a much longer growing season (GS), on average 39.9 days, than that based on SIF (Figure 6c). In addition to the differences in SOS, EOS, and GS, by the two satellite products, we also found that temporal changes in their phenophases showed divergent patterns between SIF-and NDVI-based results (Figure 6d-Figure 6f). The divergence of temporal trends in SOS was mainly found in the central and eastern TP (Figure 6d). For EOS, although asynchronous temporal trends were also found in some regions, the discrepancy may not make any difference due to their non-significant temporal changes (Figs. 4k and Figure 6e).

Discussion
Using a satellite-based SIF product, we provided a comprehensive assessment on the spatial patterns of spring and autumn phenophases, and their temporal variations on the Figure 3. Spatial patterns of the partial correlation coefficients between climatic factors and CSIFderived SOS and EOS, which are calculated by the combined mean from the five methods. (a-d) partial correlation coefficients between preseason temperature/precipitation and SOS SIF and EOS SIF . Black dots mark significant temporal changes at p < 0.05. The inserted pie charts indicate the percentage of Pn (positively non-significant value), Ps (positively significant value), Nn (negatively non-significant value), and Ns (negatively significant value).
TP over the past two decades. Consistent with earlier studies (Shen et al., 2014), our results revealed a clear spatial pattern of progressively delayed spring phenology along the eastto-west gradient. This pattern was also confirmed by MODIS NDVI-based analysis, although SOS SIF was consistently higher (later) than SOS NDVI . However, for the autumn Figure 4. Spatiotemporal patterns of SOS and EOS inferred from MODIS NDVI, which are calculated by the combined mean from the five methods. (a, d) spatial patterns of multi-year averaged SOS NDVI and EOS NDVI ; (g, j) spatial patterns of the temporal change rate in SOS NDVI and EOS NDVI ; (b, e, h, k) the multialgorithm standard deviation (SD) corresponding to figs (a, d, g, j), respectively; (c, f, i, l) frequency distributions of SOS NDVI (or EOS NDVI ) with each of the five phenophase detecting algorithms, corresponding to figs (a, d, g, j), respectively. Vertical lines indicate average values of each algorithm. Black dots mark significant temporal changes at p < 0.05. The inserted pie charts indicate the percentage of Pn (positively non-significant value), Ps (positively significant value), Nn (negatively non-significant value), and Ns (negatively significant value). phenology, EOS SIF showed a rather heterogeneous spatial distribution without a clear gradient pattern. Surprisingly, EOS SIF in some western areas even ended later than the wetter eastern and central TP. By contrast, EOS NDVI still displayed a distinct advancement from the east to west, that is, autumn ended earlier in the west than in the east. Considering that SOS SIF agreed well with SOS NDVI on the spatial pattern, it is unknown why EOS SIF did not repeat the spatial pattern found in EOS NDVI , despite the same phenophase extraction algorithms applied to these two datasets. One possible explanation is that the seasonal cycle curve of SIF may not be as symmetric between spring and autumn as that of NDVI. SIF curves are more similar to that of NDVI in the spring than in the autumn (Jeong et al., 2017). Our further analysis also found SIF and NDVI in spring had a higher relationship than that in the autumn, for most areas (77%, Fig. S6). This asymmetry of SIF seasonal variations indicates that algorithms used to detect SOS may not perform well in detecting EOS. Moreover, as SIF is a better proxy of GPP and is more responsive to instant environmental changes (Zhang et al., 2018b), our results may also imply that the photosynthesis phenology in the autumn could be more complex than leaf Figure 5. Spatial patterns of the partial correlation coefficients between climatic factors and NDVIderived SOS and EOS, which are calculated by the combined mean from the five methods. (a-d) partial correlation coefficients between preseason temperature/precipitation and SOS NDVI and EOS NDVI . Black dots mark significant temporal changes at p < 0.05. The inserted pie charts indicate the percentage of Pn (positively non-significant value), Ps (positively significant value), Nn (negatively non-significant value), and Ns (negatively significant value).
senescence. Plants may have stopped photosynthetic activities well ahead of observed leaf senescence (Bauerle et al., 2012;Jeong et al., 2017;Zhang et al., 2020). For example, the autumn phenology of photosynthesis (or EOS SIF ) may be ended by surplus carbohydrates (Fu et al., 2014;Zani, Crowther, Mo, Renner, & Zohner, 2020). Furthermore, EOS SIF is more controlled by photoperiod, even when the leaf still remains green (Bauerle et al., 2012). Figure 6. Spatial distribution of the difference in spring and autumn phenophases and in growing season (GS) between results derived with MODIS NDVI data and with SIF data (NDVI-SIF). (a-c) the difference (DIF) between multi-year average SOS NDVI and SOS SIF (a), EOS NDVI and EOS SIF (b), GS NDVI and GS SIF (c); (d-f) the temporal change rate of SOS (d), EOS (e), and GS (f) from 2001-2018. Note in (d-f), the matrix legend indicates both the difference in temporal change rates (DCR) between NDVI-and SIF-based results, and their average temporal change rates (ACR). For example, the blue color in (d) indicates the average of SOS NDVI and SOS SIF has a negative (advancing) change rate, and SOS NDVI has a higher (more negative) change rate than that of SOS SIF . The inserted histograms indicate frequency distribution of the DIF or DCR.
Our results also indicated that the spring advancement had continued in the eastern TP, but not in other parts of the TP. Indeed, SOS SIF in the western TP even showed delayed changes, although such changes were not statistically significant. This finding was also supported by analyzing MODIS NDVI data, and is consistent with the conclusion of Shen et al. (2014). The distinct spatial patterns are mainly caused by their different driving forces: the eastern and western ecosystems are subject to cold and drought stresses, respectively; and in turn are more sensitive to temperature and moisture variations (Fig.  S5, Shen et al., 2014). Whether the observed spring phenology change will persist into a warmer future will also depend upon the precipitation change, especially in the western TP. As SIF is less impacted by snow compared to vegetation index-based satellite products (Wang et al., 2019a;Zhang et al., 2018a), our validation based on ground observation further demonstrated the high performance of SIF compared to NDVI among the five algorithms (Fig. S7). This finding may help resolve the debate as to whether the spring phenology has continued to advance into the 21 st century (Klein et al., 2014). We suggest that the answer is regional dependent, as the warmer eastern TP continued to witness a more advanced spring ( Fig. S3 and S4). For the autumn phenology, however, we did not see significant changes for both SIF and NDVI. Meanwhile, significant climate changes in autumn were still observed in many different locations within the TP, for instance, the preseason warming in autumn ( Fig. S3 and S4). Therefore, the sensitivity of autumn phenology to climate change may have decreased (Meng et al., 2018). On the other hand, leaf longevity and annual total photosynthesis capacity may also be limited (Green, Berry, Ciais, Zhang, & Gentine, 2020;Kikuzawa & Lechowicz, 2011;Rollinson, 2020;Zani et al., 2020), which may also contribute to the observed insignificant changes in both EOS SIF and EOS NDVI . In addition, although we found similar temporal trends and smaller standard deviations among different algorithms, we lack long-term autumn observations to validate these retrievals. Uncertainties of different algorithms may also cause the nonsignificant trend of autumn phenology changes. Future studies need more ground-based observations to verify the result validity derived from big Earth data.
In this study, we considered only temperature and precipitation impacts on phenology. In the TP, where solar radiation is extremely high (Norsang, Kocbach, Stamnes, Tsoja, & Pincuo, 2011;Wang & Qiu, 2009), this high radiation may inhibit plant activities (Cong et al., 2017). Such radiation impacts on photosynthetic phenology (SIF) and structural phenology (NDVI) may also be different. However, currently, it is difficult to know how different these impacts may be, particularly considering it is hard to distinguish the radiation effects from that of other confounding climatic and environmental factors on the TP.
Our work was also one of the few that compared key phenological dates with different algorithms (see also Cong et al. (2012); Shen et al. (2014)). The general consistency among SOS and EOS derived with almost all the algorithms proved the robustness of our findings. On the other hand, significant inter-algorithm differences did exist. In particular, the uncertainty measured by the standard deviation among these algorithms increased from the east to the west, similar to the finding of Shen et al. (2014). This may suggest that phenology retrieval from arid and semiarid grassland is more uncertain relative to other ecosystems. A similar phenomenon has been previously documented for the arid ecosystems of Northern Australia (Wang et al., 2019a). Arid and semiarid grasslands generally have a shorter growing season and sparser vegetation coverage, which may be responsible for their higher uncertainties in satellite derived phenological dates. SIF derived growing season length was consistently shorter than that derived from NDVI across most of the areas. This is because photosynthesis starts after the unfolding of plant leaves and stops before leaf senescence. However, the large difference in their derived growing season lengths between the two satellite datasets (on average 39.9 days) may also imply that the phenophase detection algorithms previously developed for vegetation index products may not be directly applied to SIF products, at least for accurately dating key phenophases. These algorithms need to be modified or to be re-parameterized for the purpose of SIF-based accurate phenology dating. To do so, we need to have them evaluated and calibrated with ground observation data. Therefore, at least in the case of the TP, while satellite big data provide us unprecedented opportunities to explore many interesting questions, more ground observation sites broadly distributed on the plateau are still needed for calibrating satellite products.
The similarity and difference between phenophases derived from SIF and NDVI could also shed light on mechanisms driving ecosystem phenology and their carbon cycle consequences. Importantly, in spring, the SIF-and NDVI-based phenology had similar climatic drivers, the dominant driver was transformed from temperature in the eastern TP to moisture in the western TP, in line with a previous study (Shen, Piao, Cong, Zhang, & Jassens, 2015). We also found that SOS NDVI started earlier than SOS SIF . Several reasons may explain this difference between SOS NDVI and SOS SIF . Firstly, the photosynthesis activity relies on leaf structure; therefore, SOS SIF lags behind SOS NDVI . Secondly, divergent sensitivities to environmental variations, such as light, could also lead to a mismatch between SOS derived from different satellite products . In autumn, more moisture could delay autumn phenology in most regions, but the effect of temperature became diversified for EOS SIF and EOS NDVI . Specifically, EOS SIF was advanced by warming in the east but delayed in the west; however, EOS NDVI was delayed in most regions. This divergency between EOS SIF and EOS NDVI indicates the complexity of autumn phenology. Also, an earlier study (Piao et al., 2008) has shown that autumn warming actually leads to carbon loss because of enhanced respiration. Adding to the finding of Piao et al. (2008), we showed that SIF-based autumn phenology ended earlier than NDVI-based. Therefore, vegetation in late autumn could have stopped photosynthetic carbon uptake, yet still maintain its leaves, which are costly for maintenance respiration. This further explains why autumn warming may result in net ecosystem carbon loss to the atmosphere.

Conclusions
The precise extraction of plant phenology from big data is the first critical step to estimate the effect of ongoing warming on natural ecosystems. Leveraging the SIF and the NDVI dataset along with the multi-algorithm of phenological extraction, our work suggests that SIF has a high potential as a reliable data source for phenological studies, especially the SOS SIF . In particular, the comparison between phenophases derived with SIF and with MODIS NDVI complements each other in providing a more complete understanding of vegetation phenology from the structural (leaves) and functional (photosynthesis) prospective. Our finding of regionally dependent spring phenology highlights the divergent responses of the TP alpine ecosystem phenology to climate change. With the use of extensive satellite observation data, this study also underscores the value of big Earth data in research addressing important ecological and geographical questions on regions with scarce field data.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This study was supported by the National Natural Science Foundation of China (41861134036, 41988101) and by the Second Tibetan Plateau Scientific Expedition and Research Program (Grant No.2019QZKK0405)

Notes on contributors
Dr. Shilong Piao is a full Professor of Peking University. His current research focuses on data-model integration to improve our ability for predicting terrestrial carbon cycle responses to global change. He has contributed to the 5 th and 6 th Assessment Reports of the IPCC. He is on the editorial board of Global Change Biology, National Science Review, Agricultural and Forest Meteorology and so on.

Dr. Fandong Meng is a Post-doctor in the Sino-French Institute for
Earth System Science, College of Urban and Environmental Sciences, Peking University. His main research interests focus on the responses of plant phenological sequences to global change based on ground observations and remote sensing, and their feedbacks.
Ling Huang is a PhD candidate in the Sino-French Institute for Earth System Science, College of Urban and Environmental Sciences, Peking University. Her current research interest mainly focuses on the effects of climate changes on water cycle based on numerical modelling and remote sensing.