## 1 Introduction  
Terrestrial vegetation gross primary productivity (GPP) is a measure of the amount of carbon absorbed by plants through photosynthesis. It plays an indispensable part in quantifying the global carbon cycle. Accurately quantifying GPP is essential for improving our understanding of ecosystem carbon dynamics [1], climate change [2], crop production [3], and human welfare [4]. Satellite-based light use efficiency (LUE) models have become an important method for quantifying regional and global GPP, since they can provide spatially and temporally continuous information derived from remote sensing data [5; 6; 7; 8]. Numerous satellite-based models have been developed to estimate GPP at regional and global scales [7; 9; 10], and several global GPP products have been generated based on these models [11; 12; 13]. Current GPP products are mainly based on medium or coarse resolution satellite datasets and cannot reveal the spatial heterogeneity of GPP [14; 15]. For example, Global Land Surface Satellite (GLASS) GPP [13] has a spatial resolution of 500 m \(\times\) 500 m, which limits its use in applications that require GPP with finer spatial resolution, such as monitoring deforestation and restoration [16], quantifying GPP in urban areas [17], and estimating crop yield [18]. These GPP products rely on the biome-specific parameters applied from a biome parameter look-up table (BPLUT). The BPLUT parameters are provided at coarse spatial resolutions and thus do not represent the full characteristics within a given biome, especially in heterogeneous areas, due to the mismatch between the coarse spatial resolution and the eddy covariance footprints. Although BPLUT-based methods allow global estimation of GPP, the coarse inputs and BPLUT approach can introduce considerable uncertainties into GPP simulations [19]. One potential method to solve this problem involves incorporating footprints into model parameterization for better constraining parameters, making the parameters more representative [20]. Another involves using remote sensing data as inputs with high spatial resolution in the modeling [15]. Eddy covariance (EC) measurements are commonly used to validate, develop, and optimize terrestrial models because they can provide standardized and high-quality flux data within their footprint area. Although the footprints concept has been applied in some data assimilation studies [21; 22], it has often been neglected in model parameterization due to the spatial mismatch between footprints and remote sensing data [20]. In addition, flux tower sites are generally located in areas with a heterogeneous surface, and it is difficult to assess whether flux-tower observations are equal to model assumptions because models are not generally optimized to the resolution of the tower footprint [20]. Surface spatial heterogeneity has a strong influence on flux measurement, and therefore an accurate representation of the EC flux footprint is important for understanding the carbon cycle models [21]. In the cases of complex heterogeneity on land surface and EC tower sites, it is useful to evaluate their footprints and integrate the footprints into the models, thereby ensuring that the model parameters are more representative and improving model estimates of GPP. Most studies have parameterized and evaluated LUE models for GPP estimation using spatially-averaged MODIS data in windows of 6 \(\times\) 6 grid cells centered round tower sites, fully covering the tower footprints [1; 23]. MODIS data are provided in grid cells of about 10\({}^{5}\) m\({}^{2}\)\(\sim\)10\({}^{6}\) m\({}^{2}\), which is close to the spatial scale of most tower footprints. However, data in windows of 3 \(\times\) 3 grid cells cannot capture all the details of a footprint, such as the shape and trajectory of the fluxes, because flux footprints cover a relatively small area (10\({}^{3}\) m\({}^{2}\)-10\({}^{6}\) m\({}^{2}\)).