## 2 Materials and Methods
### Global VPP Estimates from MODIS EVI and Land Surface Temperature  
A combined function of EVI and LST can well explain variation in GPP because of close correlations between EVI and light-use efficiency and between LST and environmental variables, such as radiation and vapor pressure deficit [35; 36]. This provides a sound theoretical basis for estimating VPP using a regression model driven by EVI and LST. Multiple linear regression models for different PFTs were developed based on the relationships between site-level VPP observations and MODIS EVI and LST data [24]. The EVI data from the MODIS/Terra Vegetation Indices 16-Day L3 Global 250-m SIN Grid (MOD13Q1) [37] and the LST data from the MODIS/Terra LST/Emissivity 8-day L3 Global 1-km SIN Grid V006 (MOD11A2) [38] were used during the process of model construction. For each site, the EVI and LST data were extracted from a single MODIS pixel covering the FLUXNET site. Low quality values for both EVI and LST data were removed based on quality assurance (QA) layers, such as pixel reliability values from the MOD13Q1 dataset for the EVI data and mandatory QA flags from the MOD11A2 dataset for the LST data. The EVI data were masked when their corresponding pixel reliability was equal to 2 or 3, and the LST data were removed when their corresponding mandatory QA flags were equal to \"10\". If the percentage of poor values was > 30%, or if a long period of poor values was observed over two sequential time steps (32 days for EVI data and 16 days for LST data) during DOY 65 to 334 (i.e., the period of active plant growth), the site-year of the data was not used. Following this process, both 16-day EVI and 8-day LST values were linearly interpolated into daily values. Then, the same was done with the daily GPP estimates, and both site-level daily EVI and LST values were smoothed using the Savitzky-Golay filter. A previous study has reported on the strong explanatory ability of the mean and standard deviation of seasonal EVI and LST for SOS and EOS modeling [24]. In addition, we hypothesize that inter-annual variations in SOG and EOG are strongly related to EVI and LST values at five key time points during the year that shape the curves of daily EVI and LST. Therefore, twenty-eight original independent variables were obtained from smoothed daily EVI and LST data for modeling site-level VPP estimates (Table S2). These include the mean, standard deviation, and coefficient of variation of seasonal EVI and LST [24; 31]. The EVI and LST values at five key time points, including the start of each season (spring, DOY 60; summer, DOY 152; autumn, DOY 243; winter, DOY 334 for the Northern hemisphere and spring, DOY 243; summer, DOY 334; autumn, DOY 60; winter, DOY 152 for the Southern hemisphere) and the maximum value during summer, were used in this study, as they shape the curves of daily EVI and LST [31]. A detailed description of these independent variables is shown in Table S2. We used a stepwise selection method to select the independent variables from the 28 original independent variables using significance levels of 0.05 and 0.1 (default values in SPSS software) as thresholds for adding and removing independent variables, respectively. Thus, we developed PFT-specific multiple linear regression models for estimating the SOG (Appendix A) and EOG (Appendix B) of each PFT. The PFT-specific multiple linear regression models were then applied to estimate global VPP. The 16-day EVI from the MOD13C1 v006 product and the 8-day LST from the MOD11C2 v006 product at a 0.05\({}^{\circ}\) spatial resolution were used for global-scale VPP estimation. The independent variables were extracted from the smoothed EVI and LST time series. The global PFT map is from the MCD12C1 v006 product based on the IGBP classification schemes. Finally, the independent variables were added to the multiple linear regression models to produce a distribution map of global SOG and EOG estimates for 2001-2017.