Comparison of Phenology Estimated From Monthly Vegetation Indices and Solar-Induced Chlorophyll Fluorescence in China

Phenology is an important biological indicator for monitoring terrestrial ecosystems and global change. Solar-induced chlorophyll fluorescence (SIF) emitted by chlorophyll has been proven to characterize vegetation photosynthesis and phenology. In this study, we used monthly normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), and SIF products to qualitatively compare the effectiveness at detecting the phenological characteristics (SOS (start-of-season), EOS (end-of-season), and LOS (length-of-season)) over China during 2007–2013. The phenological characteristics determined by gross primary productivity (GPP) were applied as the reference to validate the phenological characteristics derived from NDVI, EVI, and SIF. The results demonstrated that the phenological characteristics derived from SIF were more consistent with that of GPP than VIs (NDVI and EVI) when considering all latitude grades, different elevation grades, and different land cover types in China. In the middle- and high-latitude regions, SOS derived from the vegetation indices (SOSVIs) did not deviate from those from GPP (SOSGPP) and SIF (SOSSIF), while in low latitudes, SOSVIs were about 20 d later than SOSSIF and SOSGPP. The VIs (EOSVIs) had a severe lag behind those of SIF (EOSSIF) in estimating the EOS at all latitudes. The EOSSIF had a deviation of fewer than 5 d compared with EOS estimated by GPP (EOSGPP), whereas the deviation of EOSVIs from EOSGPP was about 10–31 d across low to high latitude regions. The biases of SIF and VIs were due to the inconsistency between vegetation photosynthesis and leaf greenness. Also, VIs overestimated the LOS at all latitudes, the difference of LOS between estimated by NDVI and estimated by GPP was as long as 39 d in the high-latitude region. Our study suggests that SIF is suitable for estimating the phenological characteristics of vegetation regardless of different latitudes, elevation grades, and land cover types in China, providing a basis for SIF to study the vegetation phenological characteristics in a regional scope.


INTRODUCTION
Vegetation phenology is a significant factor in measuring the dynamic change of vegetation (Zhao et al., 2015) and influences such things as the terrestrial energy exchange and the global carbon cycle (Keenan et al., 2014;Tang et al., 2015). Research on vegetation phenology can further our understanding of the change of seasonal phenomena (Vrieling et al., 2013), agricultural production (Tao et al., 2006), and global change (Cleland et al., 2007). Vegetation phenology is especially strongly correlated with climatic fluctuation (Richardson et al., 2013). Global climate change seriously affects the growth of terrestrial vegetation and changes the distributions, growth cycles, and species composition of vegetation (Zhao et al., 2015). Therefore, accurate monitoring of phenology is of great significance for understanding the exchange of material and energy between vegetation and climate as well as for the accurate assessment of vegetation productivity and the global carbon budget (Piao et al., 2006;Yu et al., 2017).
Vegetation gross primary productivity (GPP) is the total photosynthetic CO 2 fixation of vegetation at the ecosystem level (Xia et al., 2015). Several studies reported that canopy phenology (using the camera) was consistent with the phenology from GPP observations in the field, such as Bartlett Experimental Forest and Howland Forest in the United States , the Lägeren Forest in northern Switzerland, and the Hainich Forest in central Germany (Ahrends et al., 2009). Therefore, GPP has been widely used as a reference for land surface phenology derived from satellite data analysis (Gonsamo et al., 2012;Jeong et al., 2017;Noumonvi et al., 2021).
In the past three decades, phenological vegetation information was monitored through time-series vegetation index data, which were obtained from satellite sensors, because these data have advantages in high temporal resolution and wide spatial coverage . Examples include the National Oceanic and Atmospheric Administration's (NOAA's) Advanced Very High-Resolution Radiometer (AVHRR) (Justice et al., 1985;Reed et al., 1994), Landsat (Melaas et al., 2013), MODerate Resolution Imaging Spectroradiometer (MODIS) (Zhang et al., 2003;Keenan et al., 2014;Walther et al., 2016), and Satellite Pour l'Observation de la Terre (SPOT) VEGETATION (Myneni et al., 1997;Zhou et al., 2001;Zhao et al., 2016). The normalized difference vegetation index (NDVI), which is associated with the greenness, seasonal, and interannual variabilities in most vegetation growth and activity, has been widely applied in the seasonal dynamic characteristics research on vegetation from the regional scale to the global scale (Tucker, 1979;Piao et al., 2006;Melaas et al., 2013). Compared with NDVI, the enhanced vegetation index (EVI) is less prone to saturation, which is caused by the luxuriant vegetation, and thus there are also many studies using EVI to observe vegetation growth (Sims et al., 2006;Testa et al., 2018). Although satellite vegetation indices (VIs)-based phenology is widely used to monitor seasonal growth and productivity of vegetation except for evergreen forests (Keenan et al., 2014;Testa et al., 2018), inconsistencies have been found between reflectance-based phenology and true vegetation photosynthesis based on ground observation (Churkina et al., 2005;Joiner et al., 2014;Cui T. et al., 2017;Jeong et al., 2017;Sun et al., 2018). In addition, NDVI-based phenological characteristics have also led to a longer length-of-season (LOS) than that derived by gross primary productivity due to an earlier estimated start-of-season (SOS) and later estimated end-of-season (EOS) in Harvard Forest (Gonsamo et al., 2012;Jeong et al., 2017;Lu et al., 2018). These discrepancies could be explained by the differences in temporal and spatial scales and dissimilarity between vegetation greenness and functions (Kross et al., 2011;Zhang et al., 2017). It has been suggested that the reflectance-based phenological detection method may not be an accurate proxy for plant photosynthesis (Frankenberg et al., 2011;Yang et al., 2014), and therefore, additional observations to help detect photosynthetic activity are necessary.
More recently, satellite remote sensing of solar-induced chlorophyll fluorescence (SIF) observations provides a new method for monitoring plant photosynthetic activities (Frankenberg et al., 2011;Guanter et al., 2014). SIF is a byproduct of the photosynthesis process of a plant, and approximately 1% of absorbed photosynthetically active radiation (APAR) is re-emitted at a longer wavelength (640-850 nm) by chlorophyll, which has two emission peaks: a red spectral region with a peak at approximately 685 nm and a near-infrared region (NIR) with a maximum at approximately 740 nm Lu et al., 2018). Using high spectral resolution spectrometers and the advanced retrieval methodology based on the exploitation of Fraunhofer lines, SIF can be retrieved from platforms such as the Greenhouse gases Observing SATellite (GOSAT) (Frankenberg et al., 2011;Guanter et al., 2012), Global Ozone Monitoring Instrument-2 (GOME-2) (Joiner et al., 2013;Zhang et al., 2016), Orbiting Carbon Observatory (OCO-2) (Frankenberg et al., 2014;Li and Xiao, 2019), and TROPOspheric Monitoring Instrument (TROPOMI) . Considering the longer time series of GOME-2 data, a number of studies applied it for long time-series research Cui Y. et al., 2017;Hu and Mo, 2020).
SIF data showed a significant positive linear relationship with GPP across the growing season, which suggests that SIF has good predictability for GPP from the canopy to satellite measurements (Damm et al., 2015;Yang et al., 2017;Li et al., 2018). Furthermore, some studies have shown that SIF could estimate the onset of spring and decline of photosynthesis in autumn more accurately in different biomes and regions, and thus can effectively estimate vegetation phenology Huang et al., 2021). Moreover, SIF had been demonstrated to be more sensitive for tracking short-term changes in GPP than reflectance-based vegetation indices (Damm et al., 2010;Yang et al., 2015;Luus et al., 2017;Sinha et al., 2017). However, most of the previous research using SIF to assess the vegetation phenology focused on the regions in relatively high latitude, in which the phenology is easily detected by SIF, few studies have compared vegetation indices and SIF in estimating phenological characteristics over a lengthy period at various latitudes. Also, few works estimated the phenology in different types of land cover or elevations in China (Jiao et al., 2020;Cheng et al., 2021).
In this study, we looked to quantify the improvement in detecting phenology by using SIF compared to NDVI and EVI products in different land cover types (evergreen forests, deciduous forests, shrublands, grasslands, and croplands), elevation grades (divided into 10 levels, according to the equal number of pixels [Supplementary Table S1]), as well as latitudes (high-(>40°N), middle-(30°N-40°N), and low-(<30°N) latitude regions). The reason for latitude grades division is in Synchronization of seasonal fluctuations between GPP and SIF, as well as VIs in each grid cell over China in China ( Figure 1). Specifically, we compared phenological characteristics derived from MODIS NDVI data with those derived from GOME-2 SIF data in deciduous forests, shrublands, grasslands, croplands, and evergreen forests using GPP data obtained from MODIS and the Max Planck Institute (MPI) as references. Considering that SIF data has been available since 2007, and GPP data accessed by MPI is only available until 2013, thus we limit the study year to 2007-2013. The main objectives of this research are twofold. One is to compare the performance of SIF-based and NDVIbased satellite remote sensing data in monitoring the satellite-GPP-based vegetation photosynthesis. The other is to study the phenological characteristics of SIF and NDVI in different regions (different latitudes, elevations, and land cover types), taking GPP as the reference.

Vegetation Indices Data
The monthly NDVI and EVI data from 2007 to 2013 were provided by the MODIS Terra Vegetation Indices dataset (MOD13C2) with a spatial resolution of 0.05°×0.05°. The MODIS NDVI dataset is based on specific spectral bands (red and infrared) designed to offer stable global vegetation leaves greenness monitoring. It has been widely used in land cover classification, carbon exchange modeling, and phenology studies (Shen et al., 2014;Aredehey et al., 2017;Biudes et al., 2021). In addition, EVI is more sensitive to high biomass areas and is minimally impacted by soil and atmosphere due to the increased blue band information (Huete et al., 2002;Jiang et al., 2008). Therefore, both NDVI and EVI were selected for phenological research in this study. In addition, NDVI and EVI were spatially resized to 0.5°×0.5°grid cells using the bilinear resampling method to make it consistent with the spatial resolution of SIF data.

Solar-Induced Chlorophyll Fluorescence Data
The SIF is acquired from the GOME-2 onboard the Eumetsat's MetOp-A platform, a sun-synchronous orbit satellite whose nominal nadir footprint size is 40 × 80 km and spectral resolution is 0.5 nm with an equator crossing time of 9:30 a.m. . The SIF signal is obtained primarily by the filling-in of solar Fraunhofer lines at wavelengths near 740 nm, which is the far-red fluorescence emission peak. A principal component analysis approach with a simplified radiative transfer model was used to separate the spectral signatures of atmospheric absorption and surface reflection from the SIF emission. The monthly SIF product had been cloud filtered and aggregated to 0.5°×0.5°spatial resolution by Joiner et al. (2013) and can be downloaded at http://avdc.gsfc.nasa.gov. In this study, we used version 2.7 Level 3 data from 2007 to 2013. The land cover was reclassified into evergreen forests, deciduous forests, shrublands, grasslands, croplands, and barren areas. The spatial resolution of the land cover map is processed to 0.5°× 0.5°. The study area was divided into the high-, middle-, and low-latitude regions by the bold and solid lines in the two panels.

Gross Primary Productivity Data
The GPP data were used as the phenological reference from two different sources. One source of GPP product was data-driven derived from the MPI (GPP MPI ). Specifically, it was estimated by a statistical model to obtain GPP estimates by upscaling global FLUXNET eddy covariance observations with the climate forcing data set CRUNCEPv6 using a machine learning method (Tramontana et al., 2016;Jung et al., 2017). GPP MPI has been widely used for GPP validation in many pieces of research (Frankenberg et al., 2011;Jung et al., 2011;Guanter et al., 2014;Jiang and Ryu, 2016). A more detailed description of GPP MPI is given in Tramontana et al. (2016). The monthly global GPP data set at the 0.5°resolution, covering 2007-2013, was downloaded from www.bgc-jena.mpg.de/geodb/projects/ Data.php. The second source of GPP data was the MODIS GPP product (GPP MOD ), which is calculated by using the MOD15A2 fraction of photosynthetically active radiation (fPAR) product and meteorological data, and applying the semi-empirical method (Running et al., 2004). In this study, we downloaded the MOD17A2 monthly product with a spatial resolution of 0.05°, which was accessed from the Numerical Terra dynamic Simulation Group (NTSG), covering the years 2007-2013. The data was downloaded from http://files.ntsg.umt.edu/data/NTSG_Products/ MOD17/GeoTIFF. Then, the bilinear resampling method was used to resize the spatial resolution of GPP to 0.5°×0.5°grid cell, which was consistent with the spatial resolution of SIF data.

Ancillary Data
The land cover data set was obtained from the MODIS dataset (MOD12Q1), the University of Maryland (UMD) land cover classification type was chosen and downloaded from http://files.ntsg.umt.edu/data/NTSG_Products/ MOD17/GeoTIFF/MOD12Q1/.
We chose the global long-term monthly means V4.01 data which covers 2007-2013, with the spatial resolution of 0.5°×0.5°.
Digital elevation data were generated from the CGIAR Consortium for Spatial Information (CGIAR-CSI) Shuttle Radar Topography Mission (SRTM) version 4.1 with a spatial resolution of 250 × 250 m. The elevation files have been mosaiced into a seamless near-global coverage (up to 60°north and south) and downloaded from https://cgiarcsi.community/.
Land cover data were spatially resampled to 0.5°×0.5°grid cells by the majority algorithm to be consistent with the spatial resolution of SIF and GPP MPI data, while elevation data were resampled by the mean resampling method.

Phenological Characteristics Extraction
To estimate vegetation phenological characteristics, the monthly SIF, VIs, GPP MPI , and GPP MOD data are smoothed and interpolated into a daily scale by the spline method (Chen et al., 2004), then normalized the daily scale data. In order to reduce noise and extract vegetation phenological characteristics from time-series data accurately, the Savitzky-Golay smoothing model (Chen et al., 2004;White et al., 2009;Brown et al., 2010;Zhao et al., 2016) was applied to the original SIF, VIs (NDVI and EVI), GPP MPI , and GPP MOD time series data. Smoothed and daily interpolated data have been applied to estimate fluctuation of vegetation growth or vegetation phenological characteristics with high accuracy (Jonsson and Eklundh, 2002;Piao et al., 2006;Zhang et al., 2013;Jeong et al., 2017;Chang et al., 2019;Wang et al., 2019). They have demonstrated that monthly data could provide a reliable result in phenological characterizing. The data used in this study were normalized by the annual maximum and minimum data to fill the magnitude gap among different data.
Then, the SOS, EOS, and LOS were estimated based on the daily SIF, VIs, GPP MPI , and GPP MOD . In this study, we defined SOS or EOS as a rapid continuous increase or decrease in remotely sensed time-series data after the long annual period of photosynthetic senescence (White et al., 2009). Thus, we adopted the threshold determined by White et al. (1997) FIGURE 2 | The method of estimating phenological characteristics. The "×" points are original data. The black solid line is a preprocessed seasonal curve that has been smoothed, interpolated, and normalized. The red horizontal dotted line is the 50% threshold of maximum and minimum of the preprocessed seasonal curve (the red vertical dotted line). The blue square dots are the SOS and EOS, and LOS is the difference between SOS and EOS.
Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 10 | Article 802763 which is the 50% point of annual normalized time series seasonal curve to calculate the phenology characteristics, because this threshold has been verified as being sensitive to the dynamics of underlying vegetation, canopy cover, and structure de Beurs and Henebry, 2010;Lu et al., 2018). We also calculated the LOS, which is the difference between the SOS and EOS. The schematic diagram of the phenological estimation method is shown in Figure 2, the intersection points of the 50% threshold line (the red horizontal dotted line) and the preprocessing seasonal curve (the black solid line) corresponded to SOS and EOS. Especially, the croplands in the middle-latitude region are planted twice a year, thus we estimate the first and last time to reach the 50% threshold as the SOS and EOS of these regions. Due to the data series being at the monthly temporal resolution, it may cause some deviation from the daily data in estimation phenology characteristics in the daily unit. The flux site GPP data, which has both monthly and daily data, were used to estimate the phenology characteristics, using Savitzky-Golay smoothing and spline interpolation, to determine the uncertainty of the data with different temporal resolutions, as Supplementary Figure S1 shows. It is found that after the preprocessing, the SOS and EOS derived from the daily data were 8 d different from the monthly data.
Resampling all the remotely sensed data to 0.5°spatial resolution may lead to some mixed pixels in the images, but the effect of the resampling was limited as Supplementary Figure  S2 shows. The EOS from the GPP MOD data was calculated with 0.05°(Supplementary Figure S2A) and 0.5°(Supplementary Figure S2B) spatial resolutions. We can find that the EOS distributions were not changed significantly before and after resampling. Supplementary Figure S2(C) shows the histograms of the calculated EOS. The histograms from different resolutions data were also similar, which indicates that the uncertainty from resampling preprocessing did not affect the phenology characteristics significantly.

Quantifying Synchronicity Methods
To quantify the synchronization of seasonal fluctuations of GPPs (GPP MPI and GPP MOD ), SIF, and VIs (NDVI and EVI), the mean absolute error (MAE, as Eq. 1) was conducted. In detail, we calculated MAE on monthly normalized SIF, NDVI, and EVI with GPP MPI and GPP MOD during spring and autumn (March, April, May, September, October, and November) in each grid cell from 2007 to 2013.

MAE
1 n n i 1 where y i represents the i-th GPP value in a grid cell, x i represents the i-th SIF, NDVI, and or EVI value in the grid cell at the same position, and n is the number of GPP-SIF(/NDVI/EVI) value pairs.

RESULTS
Synchronization of seasonal fluctuations between GPP and SIF, as well as VIs in each grid cell over China To explore the synchronization of the seasonal curves of GPPs and SIF, as well as GPPs and VIs in the period of vegetation leaf development and leaf senescence, the spatial distributions of MAE in spring and autumn are shown in Figure 3. The MAE values between GPPs and SIF ( Figure 3A for GPP MPI and Figure 3D for GPP MOD ), NDVI ( Figure 3B for GPP MPI and Figure 3E for GPP MOD ), as well as EVI ( Figure 3C for GPP MPI and Figure 3F for GPP MOD ) are shown. This figure shows lower MAE values obtained from SIF and GPPs than that from VIs and GPPs, which suggests that SIF curves can produce more similar phenological characteristics to GPP. Therefore, SIF could provide a basis for the application of SIF to indicate vegetation phenological characteristics. In addition, it was found that the synchronization between SIF and GPPs was differentiated in latitudes. SIF and GPPs showed high-level synchronicity in the high-latitude regions (>40°N),  medium-level synchronicity in the middle-latitude regions (30°N-40°N), and relatively low-level synchronicity in the lowlatitude regions (<30°N). The variation of the synchronizations provides a necessity to analyze the phenological characteristics by SIF at different latitudes. Therefore, the study area was divided into high-, middle-, and low-latitude regions, respectively. Meanwhile, the synchronicity between SIF and GPPs had some minor differences at the same latitude region, which might be caused by the influence of elevation and land cover types.

Seasonal Cycles of SIF and VIs at Different Latitudes
We analyzed the phenological seasonal cycles of different land cover types (deciduous forests, shrublands, grasslands, croplands, and evergreen forests) in different latitude regions (high-, middle-, and low-latitude regions) (Figures 4-6, and 9). The value of each point in these figures was the average monthly value of the corresponding land cover in the latitude region during the 2007-2013 period. We found that the seasonal dynamic curve of SIF was closer to that of GPPs, but the curve of VIs in autumn had a significant mismatch with the seasonal curves of two products of GPP in most cases. Figure 4 shows that the VIs curve ascended earlier than the other curves derived from GPPs in spring, while in autumn, it descended later than GPPs in the highlatitude region. Compared to GPPs, the SIF curve rose a little later than GPPs and fell together with GPPs.  Figure 5. While in contrast to the SIF, GPP MOD , and GPP MPI , the seasonal cycle of VIs was delayed in both spring and autumn in the low-latitude region ( Figure 6) on deciduous forests, shrublands, grasslands, and croplands. The difference between SOS SIF (SOS estimated by SIF) and SOS GPPs (SOS estimated by GPP MPI and GPP MOD ) (about 5 d) was smaller than that between SOS VIs (SOS estimated by NDVI and EVI) and SOS GPPs (about 18 d for EVI and about 24 d for NDVI). Additionally, the difference between EOS SIF and EOS GPPs (about 4 d) was also smaller than that between EOS VIs and EOS GPPs (about 10 d for EVI and about 28 d for NDVI). This phenomenon can also be seen in Figure 9. Furthermore, the mismatch between SOS and EOS derived from VIs, GPPs, and SIF, especially the advance of SOS VIs in the high-latitude region, and the serious delay of EOS VIs in all latitudes would lead to an overestimation of LOS VIs (LOS estimated by NDVI and EVI), especially the clear difference between LOS NDVI (LOS estimated by NDVI) and LOS GPPs (LOS estimated by GPP MPI and GPP MOD ) was about 39 d in high latitudes. In a word, considering all cases, the bias of SIF estimation of phenological characteristics was minimal. The annual changes of GPP MPI , GPP MOD , SIF, and NDVI as well as EVI of evergreen forests over 7 years (2007)(2008)(2009)(2010)(2011)(2012)(2013) are shown in Supplementary Figure S3. It was found that the changes in NDVI and EVI throughout the year are very apparent. The difference between the maximum and minimum values of NDVI did not exceed 0.15, as well as that of EVI, was no more than 0.2, and the upward and downward trends were not very distinct. Thus, it is difficult to distinguish SOS and EOS through VIs in the evergreen region. However, SIF showed a clear upward trend from March to May and a significant downward trend from September to November, which was the same as the changing trend of the two sources of GPP, and the phenological change trend could be distinguished. Therefore, the phenological characteristics estimated by SIF ( Figure 6E and LOS based on SIF and NDVI and the difference between their estimations in deciduous forests, shrublands, grasslands, and croplands in China over 7 years (2007)(2008)(2009)(2010)(2011)(2012)(2013) were mapped (Figure 7). In the high latitude areas, SOS SIF was later than SOS NDVI (SOS estimated by NDVI); in the low latitude areas, SOS SIF was earlier than SOS NDVI ; while in the middle latitude areas, SOS SIF was later than SOS NDVI in the west and earlier in the east, which may be due to elevation effect. Except for a few anomalies, the EOS SIF data was ahead of that of the EOS NDVI (EOS estimated by NDVI) data for all of China. The advance of SOS and the delay of EOS has led to an overestimation of LOS NDVI in most areas of China, except for croplands in the middle latitudes and a small area of Tibetan shrublands and grasslands. The spatial distribution of phenological characteristics estimated from EVI and the difference between SIF and EVI estimations were shown in Supplementary Figure S4. Compared with NDVI, the phenological characteristics estimated by EVI were closer to the phenological characteristics estimated by SIF. Overall, these results indicated that SIF and VIs are vastly different in estimating phenological characteristics.
The spatial distribution of the difference between the phenological characteristics estimated by GPP and those estimated by SIF, NDVI, and EVI are shown in Figure 8. The phenological characteristics of GPP were the average value of the two sources of GPP (GPP MOD and GPP MPI ) estimated phenological characteristics. It was found that the phenological characteristics estimated by SIF were more similar to those estimated by GPP, and the difference was within ±40 d. However, the difference between the phenological characteristics estimated by NDVI and those estimated by GPP was the largest. Most of the differences were delayed (that is, the difference was negative), with a difference of more than 40 d in many regions. The difference of the phenological characteristics estimated by EVI and GPP is a little smaller than that between NDVI and GPP. However, the estimated date was still delayed in many areas. Compared with VIs, SIF had the advantage to estimate EOS in that the difference between EOS SIF and EOS GPP was the smallest. EOS VIs had delayed estimation in almost the entire study area, distinguishing between −40 and −60 d in many areas. This also leads to an overestimation of LOS VIs , especially in high latitudes, and where the difference is around −50 d. Supplementary Figure S5 shows the frequency statistics of Figure 8, and it could also be found that SIF is more suitable than VIs in estimating vegetation phenological characteristics.

Phenological Characteristics Estimated
From SIF and VIs With GPP as a Reference in Different Land Covers Figure 9 represents the phenological characteristics derived from different datasets (SIF, NDVI, EVI, and GPP) in different land covers along latitudes. Different land cover types showed distinct SOS and EOS changes along with latitudes. For example, the SOS shifted to an earlier date and EOS delayed gradually for deciduous forests and grasslands from high to low latitudes and resulted in a longer LOS in lower latitudes. Nevertheless, for shrublands and croplands, the curves undulated intensively with the latitude change. For instance, the estimated SOS for croplands, especially in the middlelatitude region, showed a "concave". The reason may be intensive planting twice a year in central China, the phenology trend represents a two-peak pattern ( Figure 5D). In addition, since evergreen forest only appears in low latitudes, this study did not show the variation in middle and high latitudes.
However, the land cover types did not influence the difference in the estimation of phenological characteristics between SIF and VIs compared with GPP. Specifically, the difference between SIF and GPP is smaller than the difference between VIs and GPP in different land cover types. The differences of 7-years averaged SOS, EOS, and LOS between the NDVI and GPP, EVI, and GPP, as well as SIF and GPP are shown in Supplementary Figure S5. The phenological characteristics of GPP were the average value of the GPP MOD and GPP MPI estimated phenological characteristics. The frequency of difference values above the total and different land cover types by pixel were counted in this figure. The distribution curve conformed to an approximately normal distribution, and frequency distributions of different land cover types were similar to total frequency distribution. The differences between SIF and GPP phenological characteristics were more concentrated and closer to 0, while there was systematic bias between VIs and GPP phenological characteristics, and this caused the differences to be far from 0, especially in EOS. It manifested that SIF-estimated SOS, EOS, and LOS were all consistent with that estimated from GPP regardless of the land cover type. Conversely, the VIs-estimated SOS, EOS, and LOS were further away from those derived from GPP. The above results further proved that SIF is more available for estimating phenological characteristics in most land covers in China. Figure 10 shows the variation of SOS and EOS with elevation in different land cover types and different latitude regions. Overall, the phenological characteristics have the same tendency with the elevation rising, the SOS of the vegetation becomes later, the EOS is earlier, and the vegetation growth period is shortened. This phenomenon is particularly pronounced in the middle-and low-latitude regions, possibly because the temperature at these latitudes is strongly influenced by the elevation. In the middle-and low-latitude regions, the average elevation of the Qinghai-Tibet Plateau in the west of the study area can reach about 6000 m, while the elevation in the eastern coastal area is only about 0 m. The great variation of elevation would result in significant temperature differences, and the vegetation phenology will also change. In the high-latitude region, the variation of elevation is not clear, which may be due to the relatively small difference in elevation, and the lower temperature in high latitudes. However, the differences between SOS VIs and SOS SIF (or SOS GPPs ) and between EOS VIs and EOS SIF (or EOS GPPs ) have always existed. That is, the SOS and EOS derived from VIs deviated more from SIF and GPPs in most elevations.

DISCUSSION
The phenological changes of GPPs and SIF were consistent among the high-, middle-, and low-latitude regions, while the seasonal fluctuations of VIs were largely different from those of GPPs (Figures 4-6). In the high-latitude region, SOS VIs and SOS SIF were both close to SOS GPPs , but SOS VIs was always earlier than SOS SIF , whereas in the middle-latitude region, the seasonal curves of VIs rose synchronously with SIF and GPPs in spring phenology. However, in the low-latitude region, the VIs-based spring phenology has a lag to the SOS estimated by SIF and GPP. Meanwhile, VIs largely lagged behind SIF or GPPs in the estimation of autumn phenological characteristics at high to low latitudes, especially EOS NDVI . The estimated phenological characteristics by SIF in the highlatitude region were consistent with the results of previous research, in which NDVI showed an advance in spring and a distinct lag to GPP and SIF in autumn were found over northern forests (Jeong et al., 2017) or a specific forest such as Harvard Forest (Lu et al., 2018). Walther et al. (2016) also verified a similar seasonal change over deciduous broadleaf forests in the middle latitude along the east coast of North America, that is, NDVI lagged behind GPP, while SIF was synchronized with GPP. The distinct bias of VIs, compared with SIF, in estimating phenological characteristics might be due to the time inconsistency between photosynthesis and greenness. VIs generally indicated the greenness of vegetation, while SIF is closely related to vegetation photosynthesis. In spring, vegetation in the high-and middle-latitude regions might require time to photosynthesize and increase primary productivity by assimilating carbon after leaf emergence (Kikuzawa, 2003;Jeong et al., 2017). In addition, it caused a rapid increase in VIs values, but a slow increase in SIF and GPP values in the sprouting period. However, the result that SIF had an absolute advantage to estimate the spring phenological characteristics based on GPP in the low-latitude region is different from that in the high-and middle-latitude region and is still to be explained. On the other hand, EOS VIs have a lag at high to low latitudes, which may be due to the fact that VIs estimate vegetation phenology by tracking vegetation characteristics such as leaf area index and chlorophyll content and could not reflect slight phenological changes in time in autumn (Lu et al., 2018). In the case of autumn, carbon assimilation and SIF are mainly limited by sunlight or suitable temperature availability (Suni et al., 2003;Medvigy et al., 2013;Jeong and Medvigy, 2014;Wang et al., 2020), while vegetation indices are different from SIF, which are primarily controlled by the air temperature Yang et al., 2020). Thus, the leaves may drop after the photosynthesis has begun to decrease, which results in distinct lags in VIs, but a synchronized change with GPP in SIF (Daumard et al., 2010). Especially, EOS NDVI significantly lags behind other sources of data (GPPs, SIF, and EVI), possibly because NDVI was reported to have saturation effects in high biomass canopies (Walther et al., 2016) and it starts decreasing after a large number of leaves become senescent and drop. In addition, the phenological characteristics estimation of croplands in middle-latitude regions was only identified on the first and last date with the values of 50% thresholds as SOS and EOS. This may not be appropriate because the crops are always harvested two times per year. To detect the double peaks in midlatitude croplands is our task in future studies.
SIF has also been proven reliable to track GPP phenological fluctuations in evergreen forests, as shown in Supplementary Figure S3. This finding supports the viewpoint of Magney et al. (2019) that a good correlation between SIF and GPP in Colorado subalpine evergreen coniferous forests enables SIF and GPP to track each other consistently. The large seasonal variations in SIF yield highlight its unique ability to accurately track the seasonality of photosynthesis in different land covers, even in evergreen forests, which was not detected successfully by VIs.
Furthermore, since the temperature is one of the important climatic factors affecting vegetation phenology (Cong et al., 2013;Liang and Zhang, 2016;Lang et al., 2019;Han et al., 2020), the seasonal curve of temperature in the high-, middle-, and low-latitude regions in different land cover types is shown in Figures 4-6. We can find that the fluctuations of temperature were similar in the same latitude region in different land cover types, but the fluctuations of vegetation curves (GPPs, SIF, and VIs) were different, which indicated that the phenological characteristics of different vegetation have different sensitivities to temperature. For example, in the middlelatitude region, the seasonal fluctuations of deciduous forests were similar to the temporal change of temperature, whereas in shrublands and grasslands, the curves of GPPs, SIF, and VIs began to rise after about 2 months of temperature rise in spring. Moreover, by comparing the seasonal fluctuation of temperature and the curves of GPPs and SIF in various latitude regions, we found that the lower the latitude region, the higher the sensitivity of GPPs and SIF to temperature. These results are consistent with some previous studies (Doi and Takahashi, 2008;Gao et al., 2020). Therefore, the phenological characteristics calculated from the remotely sensed temperature have limitations in terms of the different responses of vegetation types to temperature, especially in the middle to high latitudes.
Overall, our research showed the fact that SIF tracked the variation in vegetation phenology more accurately than VIs in different latitudes, elevations, and land cover areas based on monthly data in China.

CONCLUSION
In this study, we compared and analyzed the discrepancies of phenological characteristics derived from GOME-2 SIF product and MODIS VIs products in China from 2007 to 2013. SIF showed more effectiveness in estimating the phenological characteristics than VIs in different latitudes, elevations, and land cover areas in China. The phenological curves of VIs and SIF in the high-and middle-latitude regions rose together in spring, while the VIs curve lagged far behind SIF in autumn. In the low-latitude region, in contrast to the SIF and GPP, the seasonal cycles of VIs were delayed in both spring and autumn.
Thus, SIF appears to be a better option to accurately estimate the characteristics of vegetation phenology and the phenological changes. This provides a basis for future research on the application scenarios of using SIF data to describe vegetation growth conditions and analyze vegetation phenological characteristics.
Moreover, the obviously different phenological changes between VIs and SIF in the highlatitude region, middle-latitude region, and low-latitude region provided evidence for future research on the relationship between vegetation photosynthesis and vegetation greenness across latitudes.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
XW, SL, and ZS designed the study; XW performed the analysis; XW and SL interpreted the data results; XW and ZZ drafted the original manuscript; SL, ZZ, and SZ acquired funding; and all authors contributed to the results and reviewed and edited the manuscript.