Tighten the Bolts and Nuts on GPP Estimations from Sites to the Globe: An Assessment of Remote Sensing Based LUE Models and Supporting Data Fields

: Gross primary production (GPP) determines the amounts of carbon and energy that enter terrestrial ecosystems. However, the tremendous uncertainty of the GPP still hinders the reliability of GPP estimates and therefore understanding of the global carbon cycle. In this study, using observations from global eddy covariance (EC) ﬂux towers, we appraised the performance of 24 widely used GPP models and the quality of major spatial data layers that drive the models. Results show that global GPP products generated by the 24 models varied greatly in means (from 92.7 to 178.9 Pg C yr − 1 ) and trends (from − 0.25 to 0.84 Pg C yr − 1 ). Model structure differences (i.e., light use efﬁciency models, machine learning models, and process-based biophysical models) are an important aspect contributing to the large uncertainty. In addition, various biases in currently available spatial datasets have found (e.g., only 57% of the observed variation in photosynthetically active radiation at the ﬂux tower locations was explained by the spatial dataset), which not only affect GPP simulation but more importantly hinder the simulation and understanding of the earth system. Moving forward, research into the efﬁcacy of model structures and precision of input data may be more important for global GPP estimation.


Introduction
Terrestrial gross primary production (GPP) or the total photosynthetic uptake of carbon by plants plays a critical role in maintaining the global carbon balance between the biosphere and atmosphere. However, the estimation of terrestrial GPP by existing models remains highly uncertain, with global estimates ranging widely between 92.7 to 168.7 Pg C yr −1 [1][2][3]. This large uncertainty poses a serious obstacle to quantifying and understanding the global carbon cycle [4]. It is broadly agreed that, to reduce the

Materials and Methods
In this study, using observations from global EC flux towers, we appraised 8 common satellite-data-based GPP models and the quality of major spatial data layers that drive the models. In order to utilize remote sensing data as directly as possible and minimize error propagation, 2 new GPP models were developed and compared with the 8 common satellite-data-based GPP models. Finally, 24 existing global GPP products, including the global GPP product generated from the new model that showed better performance, were compared and analyzed for similarity and differences.

Description of Existing GPP Models Used in this Study and the New Models
A total of 8 common satellite-data-based GPP models and 22 global GPP products were included in this study. The 8 existing models were the eddy covariance-light use efficiency (EC-LUE) model [9], the vegetation-indices (VI) model [25], the temperaturegreenness (TG) model [26], the vegetation photosynthesis model (VPM) [11], the carbon fixation model (CFIX) model [27], the greenness-radiation (GR) model [28], the alpine vegetation (AVM) model [29], and the moderate resolution imaging spectroradiometer (MODIS) model [30]. At the site scale, we reproduced the 8 common satellite-data-based GPP models and compared them with the two new models in this study.
The two new models were adaptations of the EC-LUE model, which included corrections for cloudiness and CO 2 concentration, with two derivations for determining water stress: One from the evaporative fraction (i.e., the LUE-EF model) and another from the normalized difference water index (NDWI) (i.e., the LUE-NDWI model), as described below.

Description of the LUE-EF Model
This model was developed mainly based on the principles of the EC-LUE model [9]. Specifically, the regulation of water on GPP is represented by the evaporative fraction (EF), taking advantage of the newly available EF products [31]. In addition, two new modifiers of GPP were added to the original EC-LUE model. The first modifier considers the impact of cloudiness on GPP. The other modifier addresses the fertilization effect of increased CO 2 concentration in the atmosphere.
The LUE-EF model can be expressed as follows: where PAR is incident photosynthetic active radiation (MJ/m 2 ) over a period of time; FPAR is the fraction of PAR absorbed by the vegetation; F CI is the regulation of cloudiness on GPP; F CO2 is the regulation scalar of atmospheric CO 2 concentration; LUE (MAX) is the maximum light use efficiency; and T S and W S are regulation scalars respectively for temperature and water stress on GPP, from which the minimum value is taken, following the Leibig law [9], and the expression of Ts is the same as that of EC-LUE model. The determination of model parameters was done as: 1.
FPAR is in practice approximated by EVI [11], since photosynthetically active vegetation is estimated as a ratio α of EVI, set to be α = 1: 2.
Most previous models underestimates of GPP on cloudy days mainly because photosynthesis can be increased by diffuse radiation under cloudy conditions [6,32]. The regulating effect of cloud cover on GPP was expressed by a cloudiness index (CI) as follow: Remote Sens. 2021, 13, 168 4 of 20 where CI is the ratio of PAR to potential PAR (PPAR) [6], and PPAR is PAR that reaches the upper atmosphere. Using the FLUXNET2015 dataset, the coefficients were determined to be a (= 2.9) and b (= 1.2); 3.
For calculating the influence of atmospheric CO 2 on GPP, we employed the algorithm in the Frankfurt Biosphere Model (FBM): where CCL is the internal CO 2 concentration of leaves, and it assumed to be 70% of atmospheric CO 2 concentration [33]. ∆(T) is the CO 2 compensation point for gross photosynthesis and photorespiration at temperature T ( • C) [34]:

4.
The regulation scalar of water on GPP, W S , was expressed as the evaporative fraction (EF) of the total sensible and latent heat [9]: where LE is latent heat flux (W m −2 ), and H is sensible heat flux (W m −2 ).

Description of the LUE-NDWI Model
The LUE-NDWI model can be expressed as follows: It can be seen that the only difference between LUE-EF and LUE-NDWI is the expression for water stress. The NDWI, strongly related to vegetation water content [35,36], can be a very good proxy for vegetation water stress. In addition, NDWI data can be obtained through satellite observation directly and the EF obtained by calculating sensible heat and latent heat, cannot be directly observed. After examining measurements from many flux towers, we found that the following nonlinear function can be used to represent W S , the regulation scalar of water stress on GPP, using NDWI: Using the 'nls' method of parameter optimization adjustment, the coefficients were determined to be a (= 0.35), b (= 2.14), and c (= 0.086). The W S_NDWI values vary between 0 and 1, with values beyond the bounds set to 0 or 1, respectively. The use of W S_NDWI can be very convenient for applications from local to global scales as the NDWI fields can be directly derived from satellite observations.

Data for Evaluation of Models and Spatial Data Products
For an evaluation of models and remotely sensed data products, we used eddy covariance (EC) flux tower data from the FLUXNET2015 dataset (https://fluxnet.fluxdata. org) [37]. Our study included data from 151 EC tower sites that belonged to the following 12 terrestrial biomes: Croplands (CRO), closed Shrublands (CSH), deciduous broadleaf forest (DBF), deciduous needleleaf forest (DNF), evergreen broadleaf forest (EBF), evergreen needleleaf forest (ENF), grasslands (GRA), mixed forests (MF), open shrublands (OSH), savannas (SAV), permanent wetlands (WET), and woody savannas (WSA). We used two criteria to filter the data, following [6]: (1) If more than 20% of the data in a given year was missing, the whole year was discarded, and (2) after this first step of processing, EC-tower sites with records for less than two years were completely discarded. EC-towers details are in Table A1. To simulate global GPP, the following spatial data products were used: (1) Meteorological data fields-radiation, air temperature, latent heat flux, and sensible heat fluxderived from the second Modern-Era Retrospective analysis for Research and Applications (MERRA-2), (2) MODIS satellite products including enhanced vegetation index (EVI) and normalized difference wetness index (NDWI), and (3) global atmospheric CO 2 concentration product from the Earth System Research Laboratory Global Monitoring Division (https://www.esrl.noaa.gov/gmd/dv/site/).

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 (α = 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.

Comparison of Model Performance at the Site Scale
Globally among all the 10 models evaluated, the LUE-EF model had the highest correlation coefficient of r = 0.86, followed by the LUE-NDWI model with r = 0.82 (Figure 1a). The LUE-EF model had the smallest normalized RMSE (0.50), indicating that the difference between the LUE-EF model and EC-tower GPP was the smallest among all models. The LUE-NDWI model had the second smallest normalized RMSE (0.52). The normalized RMSE for the other models was larger than the two new models. Model LUE-EF simulated the amplitude of the variations close to the data amplitude of EC-towers (SD ratio = 0.93). When grouped by latitudinal zone (Figure 1b-e) the improvements of the new models were more apparent for tropical and northern temperate zones. In the temperate zones, the distribution of models in the Taylor diagram was relatively concentrated, whereas in the tropical and boreal areas there were larger differences among models. When grouped by biome ( Figure A1), the new LUE-EF model showed advantages of fit in simulating daily GPP for most biomes, both in terms of correlation coefficients and RMSE. For example, correlation coefficients were highest for LUE-EF in deciduous broadleaved and evergreen needleleaf forests, wetlands, and grasslands (DBF, ENF, WET, and GRA), with LUE-NDWI being the second highest in the latter three of those. It is important to note that these biomes also have the largest number of EC-towers. It is relevant to note that many biomes are underrepresented in the current EC-tower network, such as closed shrublands or deciduous needle-leaved forests (CSH and DNF), each represented by only one and two EC-towers. The Taylor diagram shows that the models in CSH, DNF, and SAV are more dispersed, which means that model performances in these biomes vary greatly. The distribution of double mass curves ( Figure 2) show that those of the LUE-EF were the most concentrated around the 1:1 correspondence line among the 10 models compared, which indicates its greater ability to simulate the patterns in temporal variability of GPP. Figure 2 also shows the distributions of relative bias in ratio (PB-percentage biases), which ranged ± 0.4 for all models, indicating that the models had a large heterogeneity (i.e., the PB of each EC-tower is quite different) in simulating the temporal change of GPP across sites. The biases were however narrower for the LUE-EF and LUE-NDWI models, respectively, containing 120 and 110 of the EC-towers within ± 0.2 in their PB, the largest number of EC-towers amongst all models. In contrast, the number of sites with PB ± 0.2 did not exceed 100 sites for each of the other models. This indicates that the two new models had a stronger ability to capture the spatial (smaller dispersion of the double mass curves) and temporal (smaller deviation of the double mass curve from the 1:1 reference line) changes of GPP.

Biases in Remote Sensing Data Products and Consequences on Global GPP Estimation
First, various biases were found when the spatial datasets that feed the models for global GPP simulations were evaluated at the site scale (Figures 3 and A2). For example, the spatial PAR dataset only explained 57% of the observed PAR variation at the EC-towers, and the slope and intercept were 1.2 and 0.57, respectively, indicating that the PAR data fields overestimated PAR as a whole and slightly underestimated PAR at the low value. The determination coefficients of the global datasets of CO 2 , LE, and H at the EC-towers were less than 20% (R 2 < 0.2), and only temperature data was efficient in representing site-scale variation (R 2 = 0.89).
Second, the biases in the spatial datasets had a significant impact on GPP simulations was found through orthogonal regression. Before correcting these biases, the simulated GPP by the LUE-EF and LUE-NDWI models explained only 49% and 61% of the EC-tower GPP variation, and the slopes of the linear regression between simulated and towerestimated GPP were 1.54 and 1.31, respectively, and the corresponding intercepts were −2.09 and -1.23. These results indicate that both models overestimated GPP as a whole, but underestimated low GPP values (Figure 3b,f). After correcting the biases in the spatial datasets, the R 2 of LUE-EF and LUE-NDWI models improved to 0.80 and 0.79, with the slopes closer to 1 (1.20 and 1.18 values, respectively) and the intercepts closer to 0 (−0.59 and −0.91, respectively) (Figure 3c,g). The results also indicated that the LUE-NDWI model was less sensitive to the biases in the spatial data fields than the LUE-EF model, as shown by the smaller differences in R 2 before and after data correction. This can probably be attributed to the fact that the LUE-NDWI relies on NDWI, a factor that can be derived directly from remote sensing data and thereby less prone to error propagation than the LUE-EF model. Comparison of (a) spatial PAR (photosynthetic active radiation) and site PAR, (b) tower GPP and LUE-EF GPP (uncorrected spatial data), (c) tower GPP and LUE-EF GPP (corrected spatial data), (d) LUE-EF GPP with uncorrected and LUE-EF GPP corrected spatial data, (e) spatial EVI (500 m resolution) and spatial EVI (10 km resolution), (f) tower GPP and LUE-NDWI GPP (uncorrected spatial data), (g) tower GPP and LUE-NDWI GPP (corrected spatial data), and (h) LUE-NDWI GPP with uncorrected and LUE-NDWI GPP corrected spatial data. All comparisons are based on site scale. Corrected spatial data were generated from the uncorrected ones after correcting their biases using orthogonality regression (see Section 2.5).
At a global scale, the biases in spatial data inputs had a great impact on the simulated GPP even for the less sensitive model LUE-NDWI. Figure 4a shows the global average annual GPP distribution from 2000 to 2018, simulated by the LUE-NDWI model using corrected input data layers. The spatial pattern of GPP agrees well with previous studies. However, the impact of data biases on the spatial pattern of simulated GPP was obvious and not uniform across space (Figure 4b). The area overestimated is much larger than the underestimated area when the data biases were not attended, and the area fractions with GPP biases at (−50%)-(−30%), (−30%)-(−10%), 10-30%, and 30-50%, were 8%, 19%, 27%, 31%, and 15%, respectively. After data correction, the area of GPP serious reduction occurs in the mountain systems of the Tibetan plateau in Asia, northern Africa, and South America region. Growth trends of GPP were observed in Australia, northwest North America, and Siberia. The global annual average GPP estimated by the LUE-NDWI, after input data correction, was about 125.6Pg C yr −1 . Without data correction, the LUE-NDWI model would overestimate global GPP by 18% (Figure 5a,b). The corresponding global growth rate of GPP decreased from 0.34 to 0.17 Pg C yr −1 after input data correction (Figure 5c).

Comparison of 24 Global GPP Products
A comparison of 24 global GPP products is shown in Figure 5a. Large differences can be seen from these models with long-term GPP averages varying from 92.7 to 178.9 Pg C yr −1 , with more GPP estimates concentrated in the range of 120−130 Pg C yr −1 (Figure 5a,b). The number of times that a model's GPP stays within the interquartile range (IQR) of all 24 GPP products throughout the years can be used as an indicator of model performance. Result showed that the annual GPP values simulated by the following models stayed within their annual IQRs every year: Revised_ECLUE, VPM, SVR, LPJ_GUESS, and LUE-NDWI (after correcting input data biases), suggesting a better performance in simulating GPP. Although not in the category above, FLUXCOM RF, BESS, CABLE POP, and CLASS CTEM models have more than 10 years in IQR. It is worth noting that the GPP from the LUE-NDWI model, after data correction, was the closest to the median GPP value of the 24 global models.
The trends of GPP simulated by the 24 models also varied greatly from −0.25 to 0.84 Pg C yr −1 (Figure 5c). The trends of the machine learning models were smaller than those of other models. Most models showed positive trends, only the revised-ECLUE and CLASS-CTEM models showed downward trends, and some models demonstrated no significant trends (i.e., MODIS, FLUXCOM_ANN, FLUXCOM_MARS, and FLUXCOM_RF).

Adequacy of Model Structure in Representing Processes
The understanding of the processes involved in GPP is fundamental to building a reliable GPP model. For example, we found that a failed incorporation of the effect of clouds on GPP in some existing models significantly underestimated GPP in areas with frequent cloudy cover ( Figure A1). Under clear sky conditions, the upper canopy leaves are close to light saturation, while the lower canopy leaves are shaded and have limited light [51]. In contrast, under cloudy conditions, a higher proportion of the light in the form of diffuse radiation can reach the lower parts of the canopy, thus increasing the total photosynthetic use of PAR by vegetation [31]. Some studies indicated that a 1% increase in diffuse radiation induces a 0.94% increase in GPP [27,28]. In addition, many studies have shown that CO 2 fertilization has a significant effect on vegetation production, a dominant factor contributing to the 31% increase in global GPP since 1990 [52,53]. Nevertheless, many LUE models have not explicitly accounted for the effect of increasing atmospheric CO 2 concentration [6,31]. In our study, after incorporating the impacts of both cloud cover and CO 2 , the performance of the LUE-EF and LUE-NDWI models improved compared with the original EC-LUE model: R 2 improved from 0.61 to 0.68 for LUE-NDWI and from 0.61 to 0.74 for LUE-EF, respectively.
Water availability is an important factor that affects GPP [13]. In this study, we adopted two alternative structural expressions to represent the impact of water stress on GPP. First, the evaporative fraction (EF) of total energy, closely related to the Bowen ratio, was used in the LUE-EF model. The relevance of LUE-EF is grounded in the fact that less energy used for ecosystem evaporation (i.e., smaller evaporative fraction) implies a stronger water limitation [33], which has largely been verified using flux tower measurements [54]. Nevertheless, the application of the LUE-EF model is hindered by the derivation of EF that required multiple steps and input data layers and is therefore prone to error propagation [55]. In order to get a direct measure of water stress and therefore minimize error propagation, from satellite observations we replaced the EF using NDWI in the LUE-NDWI model, based on the evidence that NDWI is closely related to the plant water content and thus a good proxy for plant water stress [56]. The direct use of NDWI in the LUE-NDWI model makes it ideal for mapping GPP on regional to global scales.

Input Data Biases and Possible Impacts on GPP Simulations
Evaluating the quality of input data and understanding the impact of data biases on GPP simulation are prerequisites for improving GPP simulation accuracy. In this study, we found that the data from the spatial data fields at the EC-tower sites had various systematic deviations that seriously affected GPP estimates. For example, the spatial data fields of PAR explained only 57% of the PAR variation observed at the EC-tower sites. Reasons for the data biases are mainly rooted in the influence of sensor errors and atmospheric factors (e.g., cloud and snow) [50,57,58]. In addition, attempts unifying data from different spatial and temporal resolutions also bring biases as we found that data at different resolutions sometimes had poor correlations. The reason for the existence of resolution mismatch (i.e., spatial product data pixels cannot be unified in temporal or spatial scale) is mainly caused by mixed pixels and/or different time scales (e.g., daily, eight-day, or monthly data), compared to the spatio-temporal resolution that applies to ground conditions [51,59]. In general, spatial data biases are an important cause of the uncertainty of GPP simulation, which is an important reason for the large differences in GPP estimated among existing GPP models [60,61].
A comparison of existing global GPP products shows a huge variation from 92.7 to 178.9 Pg C yr −1 , which is not only affected by input data biases, but also by the model structure [13,62]. One should realize that input data biases affect the model outputs differently, depending on the model structure. In our study, it was found that the LUE-EF model is more susceptible to data errors, while the LUE-NDWI model is less affected by data biases (Figure 3d). The sensitivities of the GPP models to data deviations should be systematically investigated in future research as the sensitivities have not been effectively evaluated and compared.

Improving GPP Simulation Capability: The Ways Forward
A reduction in model complexity and error accumulation should be a major consideration in improving GPP estimation. For example, the evaporative fraction parameter of the LUE-EF model involves more steps than the direct use of NDWI from satellite in the LUE-NDWI model, and consequently the risk of the error accumulation increases. The direct observation and continuous recording of NDWI as a remote sensing product is one of the main advantages of the LUE-NDWI model, which minimizes the risk of error propagation of the LUE-NDWI model at a regional to global scale. Therefore, the LUE-NDWI model is more practical and attractive in spatially-explicit simulations of GPP. In summary, the simple forms of GPP models and widely and readily available inputs, as compared with more complex global models, make them more practical for applications over large areas and better suited for attribution and uncertainty analysis [53,63].
Vastly different GPP products, as shown by the means, trends, and interannual variabilities of GPP, generated by the 24 models suggest our current ability in simulating global GPP is not encouraging (Figure 5a). For nearly 40 years of commitment to the global simulation of GPP, there does not seem to have a clear direction of improvement in GPP estimation as shown by the vast differences among GPP models in simulating global GPP. GPP model development is not explicitly directed, despite the constant emergence of new models. Our research indicates that both the structures of the models and input data are error prone. Therefore, it is necessary to optimize model structures as well as sufficient validation and calibration of the input data to improve GPP simulation.

Conclusions
In this study, we first developed two new GPP models, driven by convenient remote sensing data and climate variables, with an emphasis on the representations of GPP responses to cloud, CO 2 , and water stress. Then, the performances of the new models along with other commonly used GPP models were appraised using GPP estimates from global eddy covariance (EC) flux towers. The quality of major spatial data layers that drive the models at the global scale was also evaluated at the EC-towers.
Two new models (LUE-EF and LUE-NDWI) showed advantages in simulating daily GPP for most biomes. Meanwhile, they had a stronger ability in capturing the spatial and temporal changes of GPP and the improvements of the new models were more apparent for tropical and northern temperate zones. Then, a new global GPP product was generated using the LUE-NDWI model. The LUE-EF model was not applied on 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. Various biases were found in the spatial datasets that feed the models, affecting the simulation of GPP from site to global scales. Addressing these biases is a high priority for earth system science. In addition, large differences were found from 24 global models with long-term GPP averages varying from 92.7 to 178.9 Pg C yr −1 . The newly developed LUE-NDWI model (after correcting input data biases) produced a global mean annual GPP that is very close to the mean level of all the models within the inter-quantile range of GPP from all 24 models. Moving forward, a reduction in model complexity and error accumulation, validation, and calibration of the input data fields are key processes of improving GPP simulation.
The newly developed new remote sensing-based GPP model had excellent performance and could be used to estimate GPP across a range of scales. Improving the quality of input data fields should be a major research component in reducing the uncertainty in GPP simulations on regional to global scales. The existent biases of spatial data not only affect GPP simulation but more importantly hinder the simulation and understanding of the earth system. Effective correction of spatial data is critical for reducing the uncertainty in GPP simulations and research on improving data quality should be encouraged. In this regard, our research cautions using currently available spatial data for relevant research.

Conflicts of Interest:
We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service, and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled, "Tighten the bolts and nuts on GPP estimations from sites to the globe: An assessment of LUE models and supporting data fields".

Appendix A
Existing models have estimation bias for GPP under different cloud cover conditions, and cloud cover index can adjust this phenomenon ( Figure A3). Under mostly cloudy conditions, it can be found that the overall distribution of PB in all models is biased to the negative side, which indicates that the model underestimates GPP under mostly cloudy conditions. Under partly cloudy conditions, the existing models PB is more evenly distributed around the unbiased line, yet systematically giving lower values of GPP than on clear days. Some models are more affected than others by this cloudiness-induced bias, with VPM, VI, and GR models being those most affected by GPP underestimation on cloudy days, whereas other models had less overall deviation. In order to improve model performances, the cloudiness correction factor (Equation (3)) dependent on cloud cover index (CI) was added and compared with the original model. We found that the cloudiness factor clearly improved GPP prediction for most models, except for GR, MODIS, and AVM models, which showed no obvious improvements. For example, the original MODIS model underestimates GPP under partly cloudy and mostly cloudy conditions, with no obvious deviation under sunny conditions. However, after adding CI, GPP underestimates appear under all three cloud cover conditions.
Comparing the fitting performance of the models at different scales can reflect the ability of the model to reflect spatio-temporal changes ( Table 2). As expected, the R 2 of the model was larger on the coarser spatio-temporal scales. On one hand, for the same time scale, the R 2 of the biome scale was higher than that of the site scale. Similarly, for a same spatial scale, the R 2 were increasingly larger from the daily, yearly, and to the scale of many years. This indicates that the model captures the spatial and temporal variation of GPP, but the capture capability is different. In addition to the biome_years scale, the R 2 of the LUE-NDWI model was 0.94 less than VPM model, and the R 2 of the LUE-EF and LUE-NDWI models were higher than that of other models. It can be seen that two new modes had advantages in capturing spatial and temporal changes. Second, on the scale of site_year, the slope of LUE-EF and LUE-NDWI models was 1 and extremely significant (p_value < 0.01), R 2 was also higher, 0.85 and 0.75 (p_value < 0.05), indicating that the LUE-EF and LUE-NDWI models explained 85% and 75%, respectively, of the EC-tower GPP. In addition, the slope of LUE-NDWI model in site_years was also 1 (p_value < 0.01), indicating that the model had the optimal simulation effect on these scales, that is, the applicability of the model was the most reliable.    The color dots represent the models in the corresponding legend. The Taylor diagram is a polar graph in which the cosine of the angle between the X-axis is the correlation coefficient between the GPP of the model and EC−tower. The radial direction is the ratio of model to EC−tower GPP standard deviation. The grey arcs represent the RMSE normalized by standard deviation for each model. The n is the number of EC-towers. Where (a) − (i) contains 12 terrestrial biomes: CRO (croplands), CSH (closed Shrublands), DBF (deciduous broadleaf forest), DNF (deciduous needleleaf forest), EBF (evergreen broadleaf forest), ENF (evergreen needleleaf forest), GRA (grasslands), MF (mixed forests), OSH (open shrublands), SAV (savannas), WET (permanent wetlands), WSA (woody savannas).