## 4 Discussion
First, satellite-derived VIs cannot fully capture the influences of the interannual variability in environmental drivers and biotic factors on GPP [50]. For example, moisture in the root zone, an important ecological factor, can affect annual GPP in crops and seasonally dormant forests, which was not directly observed in the VIs [50]. Second, seasonal changes in the NIR reflectance do not necessarily indicate changes in leaf reflectance, leaf area or vegetation productivity [47]. Third, cloud cover curtails the retrievals of vegetation dynamics and the presence of ephemeral seasonal snow cover can limit the identification of the seasonal change of VIs [51], resulting in the deviation of the VIs at the annual scale. Finally, the interannual variability of GPP was relatively small for some ecosystems such as tropical forests and irrigated croplands [52], and the small ranges in annual VIs and GPP could lead to weaker VI-GPP relationships. At the site and biome levels, all-around VIs with significant correlation (P < 0.05) to tower GPP only existed in ENF, GRA, WSA, and WET (Table 1) probably due to less seasonal changes of the vegetation canopies in these four ecosystems. NDVI, EVI2, EVI\({}_{\text{BRDF}}\), EVI\({}_{\text{2 BRDF}}\), and NIR\({}_{\text{V, BRDF}}\) had a relatively better performance in approximating GPP except for SAV and CRO. These results indicate that there is no vegetation index with a strong relationship (R\({}^{2}\) >0.7) with tower GPP at the annual scale across all the biomes (Figure 3). At the biome level, the performances of the VIs in predicting GPP were comparable to each other (R\({}^{2}\) = 0.84-0.86, p < 0.0001). EVI explained slightly lower variance in GPP (R\({}^{2}\) = 0.73, p < 0.0001) probably due to the low values of EVI in open and closed shrulands. Interestingly, VIs generally had similar performance to the MODIS GPP product in estimating GPP at both monthly and annual timescales. Besides fPAR, three other input variables: PAR, temperature stress, and water stress are explicitly incorporated into the MODIS GPP algorithm. This demonstrates that VIs are generally good proxies of vegetation productivity and can be used to estimate GPP across a wide variety of biomes. Meanwhile, similar to the MODIS GPP product, VIs only explained \(\sim\)20-30% of the variance in GPP at the annual scale. This indicates that caution should be taken for interpreting the results on the interannual variability and long-term trends in VIs. Capturing the interannual variability and long-term trends in GPP remains a challenge to not only satellite-derived VIs but also light use efficiency and process-based ecosystem models [53; 54]. Both observational and modeling studies have reported that the increase of the diffuse radiation improved plant LUE [55; 56; 57]. The surface reflectance and the resulting VIs do not account for the fraction of diffuse radiation, and therefore the effects of diffuse radiation were not considered. Our results also showed that these VIs responded to temperature stress fairly well. In particular, three VIs: EVI2, EVI\({}_{\text{2 BRDF}}\), and NIR\({}_{\text{V, BRDF}}\) could capture temperature stress almost as well as the MODIS GPP algorithm that incorporates a temperature scalar. A late spring frost reduced vegetation productivity as approximated by NDVI in a beech dominated Mediterranean mountain forest [58]. Satellite-derived photosynthetic activity approximated by NDVI in the Swiss Alps revealed high-elevation growth enhancement and low-elevation growth suppression in response to the 2003 heat wave [59]. Although VIs were less sensitive to water stress than to temperature stress, VIs generally captured water stress as well as MODIS GPP for most biomes. The VIs respond to dry and wet spells because of their sensitivity to water stress and have been used to examine the responses of vegetation productivity to drought [60; 61].