## 1 Introduction  
Terrestrial gross primary productivity (GPP), the integral of organic carbon fixed through photosynthesis over all of the leaves, is the beginning of the biogeochemical cycle and a primary input for vegetation carbon [1; 2]. About half of GPP will be incorporated into new plant tissues (e.g., leaves, roots), and the remainder will be released back into the atmosphere [3; 4]. Therefore, quantitative estimates of the magnitude and variability of GPP can enhance further comprehension about the biogeochemical dynamics and the carbon cycle of terrestrial ecosystems. Numerous ecosystem models, such as the Lund-Potsdam-Jena Dynamic Global Vegetation Model (LPJ-DGVM) [5], have been widely used as a practical method for quantifying temporal and spatial variation of GPP at regional, continental, and global levels. However, considering the spatial-temporal heterogeneity of ecosystem processes, it might be difficult to acquire accurate input parameters of models, thus possibly resulting in low-confidence output results [2; 6]. Moreover, the simulation errors caused by the uncertainties of models, both in process mechanisms and related parameters settings, can accumulate gradually and cause results deviating from the true values [7]. Satellite remote sensing data could offer broad spatial-temporal distribution information on vegetation coverage, as well as the state variables of models (e.g., biomass, leaf area index) [8]. Thus, it can be considered as a good option for integration with process-based models to reduce uncertainties in model estimates. Data assimilation (DA) is one technique widely used for combining observation data with process-based models and incorporates real-time data into models constantly to generate a better estimation [7], which can provide a favorable method for accurately analyzing dynamics of ecosystem, expending the use of observations in ecosystem studies. Through assimilation, observed variables are applied to adjust the model state variables; hence, several related variables will be optimized. The ensemble Kalman Filter (EnKF) algorithm has gained wide attention owing to its simple conceptual formation, optimum performance, and easy implementation [9; 10]. Based on the Monte Carlo method, EnKF designs an ensemble of state variables, and the average of the ensemble can be used as the optimal estimation of state variables. Many studies have adopted EnKF for integrating remote sensing data with models in meteorological and hydrological fields [11; 12; 13; 14]. [PERSON] et al. [15] assimilated observed soil moisture into the LPJ-DGVM by using the EnKF method to construct a data assimilation system for optimizing LPJ performance. In the field of land surface vegetation, [PERSON] et al. [16] employed EnKF in optimizing model parameters that significantly improve the simulation of GPP, latent heat (LE), and sensible heat (SH), especially during dry periods. [PERSON] et al. [17] demonstrated that the assimilation of Moderate Resolution Imaging Spectroradiometer (MODIS) reflectance data into an ecosystem model performed better in improving the estimates of GPP and net ecosystem productivity (NEP). As a sequential assimilation method, EnKF is naturally designed to deal only with sequential information; therefore, it lacks the temporal smoothness constraint included in the four-dimensional (4D) variational data assimilation (4 DVar) method. Research on coupling 4 DVar with EnKF has been conducted to fully utilize their respective advantages [18; 19; 20; 21]. Recently, a hybrid method known as the proper orthogonal decomposition (POD)-based ensemble 4D variational data assimilation method (PODEn4 DVar), was proposed by [PERSON] et al. [21].