## 3 Results and Discussion
### Uncertainties in Modeling Forest GPP  
Our results show that individual LUE models' performance varied drastically in GPP simulation across forest types at eddy covariance sites. Considering the averages of \(R^{2}\) and \(RMSE\), the VPM, EC-LUE, and MODIS products outperformed GLO-PEM, CASA, and C-Fix, while machine learning significantly improved the accuracy of estimation. The uncertainties in modeling forest GPP may come from ground-site and satellite observations, which were the model input parameters, and the structure of integrating temperature and water stress in different LUE models. First, the ground observations and subsequent processing procedures have uncertainties. GPP is not a variable that can be measured directly, and the uncertainties associated with deriving GPP products in FLUXNET2015 have been analyzed and evaluated previously in the literature [35; 63]. How to divide the NEE (Net Ecosystem Exchange) of CO\({}_{2}\) into GPP and plants' respiration and the data series' gap-filling methods are the major uncertainty source. There are two popular ways to calculate GPP from NEE measurements: daytime and nighttime fluxes method [35]. The daytime fluxes method uses a light response curve to simulate the daytime NEE, including the consideration of the temperature impact on respiration and the VPD constraint of vegetation photosynthesis [64]. The nighttime method extrapolates night-time plants respiration into the daytime by using a temperature response function that comes from long time observations [65]. Daytime and nighttime methods are popular in separating the NEE into respiration and GPP; nevertheless, both of them have advantages and disadvantages, so we used the average of FLUXNET2015's nighttime and daytime GPP products as the referenced value to validate the accuracy of our model estimation [34]. Gap-filling methods also introduce uncertainties in the calculation of GPP. The occasional instrumental failures, adverse weather conditions, and low turbulence data induced by advection issues cause missing data and require gap-filling [66; 67; 63]. The friction velocity (USTAR) indicates the turbulence strength, and a USTAR threshold was applied below which the eddy covariance NEE data were considered as bad data to be filled with predictions by a regression method [65]. These GPP derivation and gap-filling methods may introduce extra uncertainties in FLUXNET2015's GPP products. Another issue was the enclosure of the energy budget at eddy covariance sites because of wind variation, footprint size, and spatial heterogeneity. The imbalance in energy budget may induce an error of about 5% to 20%, although all FLUXNET sites were requested to obey a uniform data quality control process [35; 68; 69]. These uncertainties in eddy covariance measurements, meteorological observations, and subsequent data processing were introduced to individual LUE models and further brought into the fusion methods, which naturally impacted the accuracy of simulations. Second, the uncertainties in satellite observations contributed to the accuracy of LUE models that took remote sensing data as inputs. The errors in MODIS products used in this study have been reported [18; 36; 70] and contributed to the uncertainties in estimation of LUE models. The spatial heterogeneity was another error source. The spatial resolution of MODIS products is about 500 m, while the footprint depends on the height of the flux tower, which may be tens of meters or more, and the local terrain. Although MODIS pixels that covered the ground site were chosen carefully, mismatches still existed between MODIS pixels and footprints of ground EC sites. Third, the differences in LUE model structure can cause significantly different performance in GPP simulations [59; 71]. Temperature and water were the main constraints for most LUE models, but VPM included leaf phenology, and C-Fix considered CO\({}_{2}\) fertilization effects.