## 2 Materials and Methods
### Evaluation of Spatial Data Products and Model Performance  
In this study, data from FLUXNET eddy covariance flux tower sites were retained for analysis after data quality control. Half of the EC-towers were used for model calibration and the other half for validation. Values of model parameters were estimated using the Gauss-Newton algorithm for nonlinear optimization as implemented in the 'nls' function in the R language (package'stats' v3.6.1) [38]. The model uncertainty at the site scale was assessed using correlation coefficient (r), normalized root-mean-square error (RMSE), and standard deviation (SD) of predicted and observed values, displayed by the Taylor diagrams. In addition, the ability to capture temporal changes of GPP is essential for GPP models. The double mass curve (cumulative GPP predicted by model versus cumulative GPP observed at EC-towers) computed per site can reflect the ability of models to simulate the temporal changes of GPP, and the percentage bias (PB) of each site was counted [39].  
It is critical that spatial data fields used to drive these models should represent the site-level conditions with high accuracy when these models are applied to estimate GPP at regional to global scales. In this study, the accuracy of spatial data products were evaluated against measurements collected at the flux tower sites. Specifically, orthogonal regressions were performed, comparing each of the spatial data products against their corresponding site-level measurements [40]. Their 1:1 correspondence and biases was assessed by a hypothesis tests (\(\alpha\) = 0.05) of the regression coefficients, with null hypothesis being slope = 1 and intercept = 0 [41]. Orthogonal regression can take into account the errors of both the independent and dependent variables at the same time [42]. These analyses were conducted using the'scipy.stats' and 'odr' packages in Python (v3) [43]. When generating global GPP products, the spatial data layers (i.e., uncorrected datasets) of driving variables were corrected according to the correction coefficients obtained at the site scale through orthogonal regression. The correction coefficients were the slope and intercept of linear regression between variable values measured at the flux towers and their equivalent from the spatial datasets.  
Finally, we evaluated a total of 24 global GPP products: Two from this study (LUE-NDWI with uncorrected and corrected spatial datasets, respectively), and 22 products from TRENDY ([[http://dgvm.ceh.ac.uk/](http://dgvm.ceh.ac.uk/)]([http://dgvm.ceh.ac.uk/](http://dgvm.ceh.ac.uk/))) [44] and other studies [3]. Specifically, the 22 global GPP products were EC-LUE, revised EC-LUE [6], MODIS [25], MOD17 C6 [45], PR [46], VPM [11], FLUXCOM-ANN, FLUXCOM-MARS, and FLUXCOM-RF from FLUXCOM [47], SVR [48], BESS [49], BEPS [50], and 10 TRENDY products. The performance of the LUE-EF model was not evaluated at a global scale because of the difficulty in estimating the water stress in the LUE-EF model (i.e., the evaporative fraction of net solar radiation), which cannot be readily derived from satellite observations. We analyzed the similarity and differences of the GPP products in magnitude, spatial patterns, and temporal trends.