## 4 Discussion
### Uncertainties in the Retrieved Phenology  
Global vegetation phenological events, including SOS, EOS, and LOS, were mapped from remote sensing imagery, with their spatial distribution patterns and trends highlighted in the last two decades. The flux-based observations were used to validate the remote-sensing-derived phenological metrics. We acknowledge that several factors can lead to uncertainties in the extraction of the phenological dates. First, the fitting models of vegetation photosynthesis level can introduce extra noise in some cases. In this study, the curves of vegetation photosynthesis level were fitted using harmonic functions, which, despite being applied widely, are sensitive to missing data, possibly even resulting in degraded data quality [18; 19]. In addition, mapping vegetation phenology at the continental or global scale requires locally adjusted extraction models to provide an accurate result [22]. A number of methods, including the amplitude threshold, first-order derivative (FOD), second-order derivative, relative changing rate, third-order derivative, and curvature change rate, were proposed to retrieve SOS and EOS [14]. Although our approach, using FOD, can largely overcome the limitation of the predefined threshold value and is believed to be able to dynamically fit a wide range of conditions, it is not robust for vegetation cover types that show too slow changes or a flat shape in photosynthesis intensity during a growth cycle. Another difficulty in extracting vegetation phenology is related to multiple growing seasons (MGSs), which needs to be addressed [8; 14]. MGSs are particularly prominent in croplands characterized by the rotation of crops. In this study, the dates of the year (DOY) corresponding to the maximum and minimum FOD were taken as the SOS and EOS events for vegetation cover presenting a single growth cycle. Additionally, a flexible range of 10% relative difference was used to determine the first and last FOD points if the values were within that range for MGS cases. Nevertheless, the hard-coded threshold of 10% is likely unable to fully address the issue. Both flat-shaped and MGS patterns in vegetation growth cycles are interesting topics that deserve more attention in future studies. Figure 9: Comparison of the temporal trends of annual GPP and LOS for the significant positive correlation area, which accounts for 20.2% of the total vegetated area (see Table 1). Second, a single indicator, i.e., GPP, might not be sensitive enough to indicate the changes in vegetation phenology. Multiple indicators sensitive to reflect the changes in vegetation at different growth stages may be combined to better reflect the changes in vegetation growth stages. One of the indicators is the optimized soil-adjusted vegetation index (OSAVI), which was found to be a strong predictor of EOS for evergreen needleleaf forests [52]. Solar-induced chlorophyll fluorescence and vegetation optical depth may provide alternative methods to extract vegetation phenology [53; 54; 55]. Our study provides high-level statistics of the correlations between different ecosystem types, such as cryptlands, forests, grasslands, savanna, and shrub, and climate variables (Table 1), which is useful to understand the response of vegetation phenology to climate changes for different ecosystems. However, such statistics may not reveal detailed differences within each type. Lastly, the relatively coarse datasets applied for global-scale vegetation phenology mapping may introduce noise into the final results. The spatial resolution of MODIS data (500 m) is likely to hide the spatial detail in the patterns of vegetation phenology. Spatial resolution is an important factor in phenology extraction, and spatially aggregated GPP over a large area, i.e., 500 m \(\times\) 500 m, represented by a single pixel can easily smooth out finer location-dependent GPP across an area [56]. The coarse spatial resolution of MODIS data also leads to difficulty in the validation from field-based observations.