## 1 Introduction  
Terrestrial gross primary productivity (GPP), the amount of carbon absorbed by terrestrial plants through photosynthesis, is the largest carbon flux between the terrestrial biosphere and the atmosphere. It also drives the terrestrial food chain and is the basis for agricultural and wood production. Quantifying GPP is therefore essential for assessing ecosystem carbon dynamics, climate feedbacks, agricultural productivity, and human welfare. Satellite-derived vegetation indices (VIs) are indicative of photosynthetic activity and have been widely used as proxies for GPP [1; 2; 3; 4] or to directly estimate GPP [5; 6; 7]. Ground-based GPP data from the eddy covariance (EC) flux towers are typically used as reference data to evaluate the performance of satellite-derived VIs in approximating or estimating GPP. Better understanding the relationships between the VIs and GPP for a wide variety of biomes at various timescales is essential for assessing ecosystem functioning, vegetation productivity, and carbon budgets at regional to global scales using satellite-derived VIs. The normalized difference vegetation index (NDVI) is perhaps the most widely used VI in satellite monitoring of vegetation productivity [8]. NDVI captures the contrast in reflectance between the red and NIR wavelengths and was first developed to estimate vegetation greenness [9]. It has been widely adopted as a proxy for GPP. For example, a close relationship between NDVI and GPP at the daily timescale was found for a maize cropland [10], a deciduous forest [11] and a pine forest [12]. However, NDVI is sensitive to soil background brightness and is influenced by soil texture and moisture conditions [13]; NDVI is also sensitive to scattering by the atmosphere from highly variable aerosols [14; 15] and also tends to saturate in multi-layered and dense canopies [16; 17]. The enhanced vegetation index (EVI) was proposed to reduce these limitations by using reflectance in the blue band to correct for effects of aerosols [17; 18]. Many studies reported that EVI had a better performance in estimating GPP for alpine meadow [19], maize [20] and wetlands [21] than NDVI. An EVI-based index without the blue band - EVI2 was later developed [22]. EVI2 is similar to the 3-band EVI especially when atmospheric effects are insignificant and the reflectance data are in good quality [22]. Compared with NDVI, this 2-band index was shown to be able to better capture subtle changes in vegetation condition and structure and to have higher correlation with GPP for maize and soybean ecosystems [10; 23] and a red pine forest [24]. EVI2 has also been shown to have a slightly higher correlation with GPP than did EVI in a forest ecosystem [25]. Although these satellite-derived VIs have been widely used to approximate or estimate GPP [5; 6; 7; 8; 26; 27], it remains unclear how strongly the VIs are related to GPP and how the VI-GPP relationship varies across VIs, biomes, timescales, and different climatic conditions. Prior to December 2015, the FLUXNET database available to the research community was the La Thuile dataset released in 2007. The duration of the GPP data from this database had only up to four years of overlap with the NDVI record from the moderate resolution imaging spectroradiometer (MODIS), which hindered the analysis for the VI-GPP relationship at the annual scale. The NDVI record derived from the Advanced Very High Resolution Radiometer (AVHRR) has a longer overlap with the La Thule database, but there is large scale-mismatch between the flux tower footprint and the AVHRR grid cell (\(\sim\)64 km\({}^{2}\)). The new FLUXNET2015 dataset, released in December 2015, consists of EC and meteorological data for 221 flux towers across the globe. Compared with its predecessor, the La Thuile dataset, the FLUXNET2015 dataset had much longer temporal coverage with some sites spanning the period from the early 1990s to 2014.