Integration of Sentinel-3 OLCI Land Products and MERRA2 Meteorology Data into Light Use Efﬁciency and Vegetation Index-Driven Models for Modeling Gross Primary Production

: Accurately and reliably estimating total terrestrial gross primary production (GPP) on a large scale is of great signiﬁcance for monitoring the carbon cycle process. The Sentinel-3 satellite provides the OLCI FAPAR and OTCI products, which possess a higher spatial and temporal resolution than MODIS products. However, few studies have focused on using LUE models and VI-driven models based on the Sentinel-3 satellites to estimate GPP on a large scale. The purpose of this study is to evaluate the performance of Sentinel-3 OLCI FAPAR and OTCI products combined with meteorology reanalysis data in estimating GPP at site and regional scale. Firstly, we integrated OLCI FAPAR and meteorology reanalysis data into the MODIS GPP algorithm and eddy covariance light use efﬁciency (EC-LUE) model (GPP MODIS-GPP and GPP EC-LUE , respectively). Then, we combined OTCI and meteorology reanalysis data with the greenness and radiation (GR) model and vegetation index (VI) model (GPP GR and GPP VI , respectively). Lastly, GPP MODIS-GPP , GPP EC-LUE , GPP GR , and GPP VI were evaluated against the eddy covariance ﬂux data (GPP EC ) at the site scale and MODIS GPP products (GPP MOD17 ) at the regional scale. The results showed that, at the site scale, GPP MODIS-GPP and GPP EC-LUE agreed well with GPP EC for the US-Ton site, with R 2 = 0.73 and 0.74, respectively. The performance of GPP GR and GPP VI varied across different biome types. Strong correlations were obtained across deciduous broadleaf forests, mixed forests, grasslands, and croplands. At the same time, there are overestimations and underestimations in croplands, evergreen needleleaf forests and deciduous broadleaf forests. At the regional scale, the annual mean and maximum daily GPP MODIS-GPP and GPP EC-LUE agreed well with GPP MOD17 in 2017 and 2018, with R 2 > 0.75. Overall, the above ﬁndings demonstrate the feasibility of using Sentinel-3 OLCI FAPAR and OTCI products combined with meteorology reanalysis data through LUE and VI-driven models to estimate GPP, and ﬁll in the gaps for the large-scale evaluation of via Sentinel-3 satellites.


Introduction
Gross primary production (GPP) usually refers to the total amount of CO 2 assimilated in photosynthesis [1]. GPP has been an important indicator for quantitatively describing the global carbon cycle and evaluating the sustainable development of terrestrial ecosystems [2,3]. Therefore, the real-time monitoring of changes in GPP is essential for the field of regional and global change studies [4]. The eddy covariance (EC) technique is the most widely used approach to measure CO 2 net uptake and can obtain GPP through different modeling methods [5]. However, the limitations of this technique: the limited number and solar-induced chlorophyll fluorescence (SIF). They proved the potential correlation between FAPAR products and GPP. Harris et al. [49] demonstrated that OTCI might be an effective substitute for MODIS to estimate GPP. Zhang et al. [50] combined OLCI FAPAR with in situ meteorological data to drive the MODIS GPP algorithm to estimate GPP across several biomes. The results indicated that the Sentinel-3 OLCI FAPAR products performed better than MODIS FAPAR products at the site scale in estimating GPP and demonstrated a significant correlation between OTCI and GPP. Nevertheless, the studies mentioned above were all based on the site scale, and the estimation of GPP on a large scale was not involved. Moreover, only a single LUE model was adopted, or only a simple relationship between in situ GPP and FAPAR was established. None of these studies utilized multiple LUE or VI models to estimate GPP.
This paper aims to demonstrate the potential of combining Sentinel-3 OLCI FAPAR and OTCI with meteorology reanalysis data in estimating GPP with multiple models at multiple scales. Specifically, this study has four objectives (1) integrating Sentinel-3 OLCI FAPAR and MERRA2 meteorology reanalysis data into the MODIS GPP algorithm and EC-LUE model to estimate GPP, (2) combining Sentinel-3 OTCI and MERRA2 meteorology reanalysis data with the GR and VI models to estimate GPP, (3) comparing the GPP obtained from four RS-driven models with on-site measured GPP and MODIS GPP products, and (4) applying Sentinel-3 OLCI FAPAR and meteorology reanalysis to estimate GPP at the regional scale and compare it with MODIS GPP products.

Research Sites and Region
The interest area covers regions from the upper-left corner at 45 • N and −124 • W to the bottom-right corner at 25 • N and −89 • W in North America ( Figure 1). The AmeriFlux eddy covariance flux tower data were used in the study. Since the full year of available measurements by the satellite products is from 2017, we chose sites with data after 2017. The eight EC flux sites cover eight biome types, including evergreen needleleaf forests (ENF), deciduous broadleaf forests (DBF), mixed forests (MF), closed shrublands (CSH), open shrublands (OSH), woody savannas (WAS), grasslands (GRA), and croplands (CRO). Since the EC flux sites with the biome types of evergreen broadleaf forests (EBF), deciduous needleleaf forests (DNF), and savannas (SAV) are not matched in time with Sentinel-3 OLCI, this research did not estimate the GPP under these biome types at the site scale. Additionally, due to rich biome types, California was selected as the study area for regional GPP estimation.  Sentinel-3 OLCI is an inheritor of the Medium Resolution Imaging Spectrometer (MERIS) and has many enhancements [51,57]. There are three main kinds of OLCI products available to the public (i.e., level-1B products, level-2 land products, and level-2 water products). The level-2 land product provides land and atmospheric geophysical products at full and reduced resolution, where full resolution is approximately 300 m on the ground and reduced resolution is approximately 1200 m on the ground. An OLCI Level-2 Land Full Resolution (OL_2_LFR) product includes an OGVI band and OTCI band. OGVI describes FAPAR in the plant canopy and OTCI represents the terrestrial chlorophyll index. Sentinel-3 OLCI products applied in this study covered January 2017 to December 2018 at the US-WCr site, US-Ton site, US-PFa site and US-Var site, respectively. However, due to the lack of site data, Sentinel-3 OLCI products from January 2017 to September 2018 at the US-Rls site and US-Rws site, January 2017 to March 2018 at the US-GLE site and April 2017 to December 2018 at the US-Bi2 site were used.
The Sentinelsat module offers a powerful Python API for searching and downloading the metadata of Sentinel satellite images from the Copernicus Open Access Hub (https: //scihub.copernicus.eu/dhus (accessed on 13 January 2021)) through executing condition and spatial queries (a region of interest can be created with a polygon on the GeoJSON website) and was used to batch download the OLCI online and offline products. The filtered OLCI products were automatically downloaded and stored locally as compressed files, and a total of about 5200 scene data were downloaded in this work. Then, we developed a complete set of OLCI product preprocessing processes, including mosaic, reprojection, resampling, and subset processes. The snappy module developed by the European Space Agency (ESA) and a custom script written in Python was applied to complete these data preprocessing operations. The traditional data processing software (i.e., SNAP and ENVI) either cannot batch process or takes a long time. Our python script can dramatically improve the efficiency of data processing through fully automatic data retrieval, batch processing and multithreading technology.
In order to run GPP estimation models at the same spatial resolution as OLCI products, we adopted the nonlinear spatial interpolation method developed by Zhao et al. [60] for downsizing the MERRA2 meteorology reanalysis data from a 0.5 • × 0.625 • resolution to a 300 m OLCI pixel resolution [43,59]. This method utilizes the fourth power of the cosine function and the weighted distance in the closest four grid cell for computing the value of every output pixel with an OLCI resolution [61]. The spatial interpolation scheme may improve the meteorological input data accuracy of every 300 m pixel, eliminating sudden changes between the two sides of the MERRA2 boundary [62].

EC Flux Data
The EC flux data of eight AmeriFlux sites were downloaded from the AmeriFlux data portal (https://ameriflux.lbl.gov/ (accessed on 13 January 2021)). These eight AmeriFlux sites are GLEES (US-GLE), Willow Creek (US-WCr), SE5 Aspen-5 CHEESEHEAD 2019 (US-PFs), RCEW Low Sagebrush (US-Rls), Reynolds Creek Wyoming big sagebrush (US-Rws), Tonzi Ranch (US-Ton), Vaira Ranch-Ione (US-Var), and Bouldin Island corn (US-Bi2), represented by ENF, DBF, MF, CSH, OSH, WSA, GRA, and CRO, respectively ( Table 1). The EC flux data have a half-hourly temporal resolution, and the daily GPP values were calculated as a sum of 30-min GPP fluxes. The REddyProc online tool (https://www.bgcjena.mpg.de/REddyProc/brew/REddyProc.rhtml (accessed on 13 January 2021)) [63,64] implements a standardized approach for processing EC data and was used to calculate GPP for the EC flux tower when GPP was not provided. Furthermore, the online tool filters out data collection during stable atmospheric conditions using u* (friction velocity) filtering, and gap-fills missing data using a marginal distribution sampling approach. Flux measurements (net ecosystem exchange, latent heat flux and sensible heat flux) and meteorology measurements (incoming shortwave radiation, vapor pressure deficit, relative humidity, air temperature, soil temperature and friction velocity) were obtained to calculate GPP.

MODIS Products
The most widely used GPP product is MOD17A2H, which uses updated Biome Property look-up tables and an updated version of the daily GMAO meteorological data as input data. Since the time resolution is 8 days, there are a total of 92 scene MODIS GPP products used within two years. The purpose of using MOD17A2H is to compare of the model prediction accuracy.
The MOD12Q1 land cover product with a 500 m spatial resolution, which is one of the MODIS products, was used in the MODIS GPP algorithm. Boston University's UMD classification scheme within the dataset was also employed for offering biome-specific information for GPP estimation models [6]. Moreover, a 1.5 km × 1.5 km area centered on every flux tower site was extracted to represent the flux footprint, the mean values of the 5 × 5 pixels and the mean values of the 3 × 3 pixels were the final GPP values for Sentinel-3 OLCI products and MODIS GPP products, respectively [5,50].

MODIS GPP Algorithm
Firstly, we used the MODIS GPP algorithm to estimate GPP. The MODIS GPP algorithm is applied based on the LUE model, which mainly uses the relationship between LUE and environmental factors (temperature, water vapor, and light) [28,72]. The algorithm can be expressed as follows: where ε max , TMIN max , TMIN min , VPD max , and VPD min are shown in Table 2. The biomespecific physiological parameters for each biome type are determined according to the MOD12Q1 UMD classification scheme and BPLUT (Table 3) [59]. TMIN and VPD scalars are simple linear ramp saturated at the maximum and minimum, respectively, and the range is from 0 to 1 [73]. FAPAR was obtained from the OLCI FAPAR product (OGVI). IPAR is equal to SWRad times 0.45 (Equation (1e)): where SWRad, T min and VPD are obtained from MERRA2 meteorology reanalysis data. VPD is calculated with T2MMEAN and RH (Equation (1f)) [73]: Daily minimum temperature at which ε = 0.0 VPD max Daylight average vapor pressure deficit at which ε = ε max VPD min Daylight average vapor pressure deficit at which ε = 0.0

EC-LUE
Yuan et al. [32] developed the EC-LUE model which is mainly driven by FAPAR, PAR, the air temperature, and the Bowen ratio of sensible to latent heat flux. Because of the weak simulation of sensible and latent heat fluxes on a large spatial scale, adopting the Bowen ratio for presenting water stress factors (W s ) will hinder EC-LUE model application on a large scale [75]. In this study, VPD was used to replace sensible and latent heat fluxes to calculate W s [21]. The GPP in EC-LUE is estimated as follows: where FAPAR is obtained from the OLCI FAPAR product (OGVI). IPAR is the same as in the MODIS GPP algorithm. ε max is the potential LUE without environmental stress and is set Remote Sens. 2021, 13, 1015 7 of 23 to 2.14 (gC/m 2 /d) [33]. T s and W s are the downward-regulation scalars for the individual influences of temperature and humidity on LUE. Min denotes the minimum value of T s and W s . T s indicates a limiting influence of temperature on vegetation photosynthesis, based on the algorithm of the temperature limiting factor in the Terrestrial Ecosystem Model (TME) [76,77]: where T is the average air temperature ( • C), and T min , T max , and T opt represent the minimum, maximum, and optimum air temperatures ( • C), respectively, for photosynthetic activities. If the air temperature is below T min or exceeds T max , Ts is set to zero. According to [32], this study set T min and T max to 0 and 40 • C, respectively, and T opt was determined to be 20.33 • C through nonlinear optimization. The W s can be expressed as follows: where VPD 0 is the half-saturation coefficient of the Equation (2c), obtained from previous research (Table 4) [21].

The Greenness and Radiation Model (GR)
Gitelson et al. first introduced the GR model [78]. The GR model estimates GPP by adopting the chlorophyll content, as well as incoming solar radiation. This model was successfully applied within irrigated and rainfed maize, wheat cropland, and forest [48]. The GR model can be expressed as follows: where Chl was replaced with OTCI in this study. The parameter m is a scalar decided by the model calibration. The calculation method of model calibration is described in Section 2.3.5.

The Vegetation Index Model (VI)
The VI model employed for GPP estimation was proposed by Wu et al. and applied in crop and deciduous forest ecosystems [42]. It assumes that the VIs are a reliable proxy of both LUE and FAPAR, so the product VIs * VIs * IPAR is used to estimate GPP. In this paper, VIs is same, so the VI model can be expressed as: where OTCI is regarded as the VI parameter of the model input, m is the same as in the GR model.

Calibration of the GR and VI Model
Model calibration is an important step for the operational application of the GR and VI models and greatly impacts model accuracy [48,79]. GR and VI model calibration mainly involves the calculation of scalar m. At present, calibration methods can be divided into adopting the half-data method (the rest of the data are used for model validation) [23] and the all-data method (all data are applied to model calibration and validation) [46,48]. In this work, to conduct a rigorous test of the model, we chose half-data to calibrate the model. The rest of the data were then applied when testing the model. The values of scalar m were obtained through establishing a directly linear relationship between in situ measured GPP and model-estimated GPP values, and scalar m is the slope of linear function.

Analytical Methods
Three statistic analytical indices can be adopted for evaluating the accuracies and uncertainties of all GPP models: (1) the coefficient of determination (R 2 ), which expresses the fit between model prediction results and true values; (2) the root mean square error (RMSE), which refers to sample standard deviation between predictions and measurements; and (3) bias, which reflects the differences between the estimates and the observation. RMSE and bias are calculated as follows: where M i and E i are the predicted GPP and measured GPP values, respectively. Notably, in subsequent chapters, GPP EC represents the on-site measured GPP values, and GPP MOD17 represents the GPP values from MODIS-GPP products. GPP MODIS-GPP , GPP EC-LUE , GPP GR , and GPP VI represent GPP estimation values obtained from the MODIS-GPP algorithm, EC-LUE model, GR model, and VI model, respectively.

Meteorology Variables
In this study, the MERRA2 meteorology reanalysis data were taken as site meteorological data for four GPP estimation models. For determining the meteorological data effects on GPP estimation, we conducted a direct comparison of daily MERRA2 meteorology reanalysis data and site meteorological data ( Figure 2). The two datasets have a strong correlation with T2M_MIN and T2M_MEAN (Figure 2a,b), suggesting that the MERRA2 meteorology reanalysis data can be reliably applied to estimate the on-site temperature conditions. The IPAR (Figure 2c) at the eight sites had a high correlation with the measured data. The R 2 was above 0.8, and the scatter distribution was very close to the 1:1 line, which indicated that the on-site IPAR could be reliably calculated from the MERRA2 meteorology reanalysis data. For VPD (Figure 2d), R 2 varied at different sites, and the range of R 2 was 0.51-0.95. Among them, VPD at the US-WCr (DBF) and US-Bi2 (CRO) obtained a relatively low accuracy, with R 2 = 0.54 and 0.62, respectively. The VPD across shrublands obtained the highest accuracy (R 2 = 0.95).
conditions. The IPAR (Figure 2c) at the eight sites had a high correlation with the measured data. The R 2 was above 0.8, and the scatter distribution was very close to the 1:1 line, which indicated that the on-site IPAR could be reliably calculated from the MERRA2 meteorology reanalysis data. For VPD (Figure 2d), R 2 varied at different sites, and the range of R 2 was 0.51-0.95. Among them, VPD at the US-WCr (DBF) and US-Bi2 (CRO) obtained a relatively low accuracy, with R 2 = 0.54 and 0.62, respectively. The VPD across shrublands obtained the highest accuracy (R 2 = 0.95).

Agreement between GPPMODIS-GPP, GPPEC-LUE, GPPMOD17 and GPPEC
We first compared the performance of GPPMODIS-GPP and GPPEC-LUE against that of the GPPEC at all sites for 2017 and 2018 ( Figure 3). For the MODIS-GPP algorithm, the results showed that the US-Ton site obtained the best performance (R 2 = 0.73), followed by US-GLE with R 2 = 0.72, US-Bi2 with R 2 = 0.68, US-WCr with R 2 = 0.67, US-PFa with R 2 = 0.66, and US-Rls with R 2 = 0.63 ( Table 5). The R 2 values between the GPPMODIS-GPP and GPPEC at the US-Var site and US-Rws site were lower than 0.5. In terms of the RMSE, the US-Bi2 site produced the maximum error (RMSE = 7.70 gC/m 2 /d) and the US-Ton site displayed the lowest error with RMSE = 1.15 gC/m 2 /d.

Agreement between GPP MODIS-GPP , GPP EC-LUE, GPP MOD17 and GPP EC
We first compared the performance of GPP MODIS-GPP and GPP EC-LUE against that of the GPP EC at all sites for 2017 and 2018 ( Figure 3). For the MODIS-GPP algorithm, the results showed that the US-Ton site obtained the best performance (R 2 = 0.73), followed by US-GLE with R 2 = 0.72, US-Bi2 with R 2 = 0.68, US-WCr with R 2 = 0.67, US-PFa with R 2 = 0.66, and US-Rls with R 2 = 0.63 ( Table 5). The R 2 values between the GPP MODIS-GPP and GPP EC at the US-Var site and US-Rws site were lower than 0.5. In terms of the RMSE, the US-Bi2 site produced the maximum error (RMSE = 7.70 gC/m 2 /d) and the US-Ton site displayed the lowest error with RMSE = 1.15 gC/m 2 /d.
For the EC-LUE model, the results showed that the US-Ton site obtained the best performance (R 2 = 0.74), followed by the US-Var with R 2 = 0.67, US-WCr with R 2 = 0.58, US-Rws with R 2 = 0.57, US-PFa with R 2 = 0.57, US-PFa with R 2 = 0.56, and US-Bi2 with R 2 = 0.53. However, insignificant correlations were obtained at the US-GLE site, and this situation was mainly influenced by the ENF when using EC-LUE model [32,75]. For MODIS-GPP products, the US-GLE site exhibited the best performance, with R 2 = 0.92. At the US-Rws site and US-Var site, GPP MOD17 performed relatively poorly, with R 2 = 0.48 and 0.46, respectively. These results were consistent with GPP MODIS-GPP . It is worth noting that the performance of GPP MOD17 was worse than GPP MODIS-GPP for the US-Ton, US-Var, and US-Bi2 site and worse than GPP EC-LUE for the US-Rws, US-Ton, and US-Var site. Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 24

Agreement between GPP GR , GPP VI and GPP EC
The scatter plots between the GPP GR , GPP VI , and GPP EC for every site are displayed in Figures 5 and 6. The performance of the two models is different at various sites. Therefore, the two models' applicability depends on the biome type. For example, insignificant correlations between the GR and VI models were obtained at the US-Rws site, US-Ton site and US-Rls site. Even for the same biome type, their performance was different. For instance, the GR model had a better performance than the VI model at the US-GLE site. Additionally, strong relationships between GPP GR , GPP VI and GPP EC were established at the US-WCr site, US-Bi2 site, US-PFa site, and US-Var site. In particular, the GR and VI models both obtained the best performance at the US-WCr site, with R 2 = 0.84 and 0.81, respectively. = 2.07 gC/m 2 /d) and the US-WCr site (RMSE = 4.05 gC/m 2 /d, bias = 1.71 gC/m 2 /d) were overestimated. Notably, GPPMODIS-GPP at the US-WCr site and GPPEC-LUE at the US-Ton site obtained a low bias but high RMSE. It was because the underestimation of GPP offset part of the deviation during the months of 2017, as well as the overestimation of compensation in 2018 or GPP underestimation during the growing season, alongside overestimation during the nongrowing season.

Agreement between GPPGR, GPPVI and GPPEC
The scatter plots between the GPPGR, GPPVI, and GPPEC for every site are displayed in Figures 5 and 6. The performance of the two models is different at various sites. Therefore, the two models' applicability depends on the biome type. For example, insignificant correlations between the GR and VI models were obtained at the US-Rws site, US-Ton site and US-Rls site. Even for the same biome type, their performance was different. For instance, the GR model had a better performance than the VI model at the US-GLE site. Additionally, strong relationships between GPPGR, GPPVI and GPPEC were established at the US-WCr site, US-Bi2 site, US-PFa site, and US-Var site. In particular, the GR and VI      Figure 7. Overall, the temporal dynamics of GPPGR and GPPVI are roughly the same. The GPPGR and GPPVI tracked the GPPEC well at the US-Bi2 site, US-PFa site, US-Var site, and US-WCr site. As shown in Table 6, GPPGR and GPPVI at the US-GLE site were overestimated, with bias = 3.18 gC/m 2 /d and 1.16 gC/m 2 /d, respectively. Additionally, GPPGR and GPPVI were underestimated at the US-Ton site, with bias = -0.75 gC/m 2 /d and -1.48 gC/m 2 /d, respectively. Notably, both GPPGR and GPPVI showed a poor relationship with GPPEC at the US-Rls site and US-Rws site. This situation was mainly due to the underestimated GPP values during the growing season and overestimating GPP values during the nongrowing season.  Figure 7. Overall, the temporal dynamics of GPP GR and GPP VI are roughly the same. The GPP GR and GPP VI tracked the GPP EC well at the US-Bi2 site, US-PFa site, US-Var site, and US-WCr site. As shown in Table 6, GPP GR and GPP VI at the US-GLE site were overestimated, with bias = 3.18 gC/m 2 /d and 1.16 gC/m 2 /d, respectively. Additionally, GPP GR and GPP VI were underestimated at the US-Ton site, with bias = −0.75 gC/m 2 /d and −1.48 gC/m 2 /d, respectively. Notably, both GPP GR and GPP VI showed a poor relationship with GPP EC at the US-Rls site and US-Rws site. This situation was mainly due to the underestimated GPP values during the growing season and overestimating GPP values during the nongrowing season.

Spatial−Temporal Consistency between GPP MODIS-GPP , GPP EC-LUE and GPP MOD17
We compared the spatial distribution of the annual mean GPP (gC/m 2 /d) and maximum daily GPP (gC/m 2 /d) from GPP MODIS-GPP and GPP EC-LUE across California in North America during 2017 and 2018 at a 300 m spatial resolution. The annual mean GPP was the average of each pixel throughout the year. The maximum daily GPP was the maximum value of each pixel throughout the year. Notably, the spatial distributions of GPP GR and GPP VI were not shown here due to the insufficient calibration data for upscaling to the regional scale. . The biggest discrepancy between the annual mean GPP value and maximum daily GPP value could be found in the central croplands, where the maximum daily GPP was high, and the annual mean GPP was moderate. Overall, the GPP MODIS-GPP and GPP EC-LUE showed similar spatial−temporal patterns to GPP MOD17 . Notably, compared with GPP MOD17 , GPP MODIS-GPP , and GPP EC-LUE relatively underestimated the GPP values for cropland.
37°N, −121°W-199°W). The biggest discrepancy between the annual mean GPP value and maximum daily GPP value could be found in the central croplands, where the maximum daily GPP was high, and the annual mean GPP was moderate. Overall, the GPPMODIS-GPP and GPPEC-LUE showed similar spatial−temporal patterns to GPPMOD17. Notably, compared with GPPMOD17, GPPMODIS-GPP, and GPPEC-LUE relatively underestimated the GPP values for cropland.

Discussion
The R 2 between GPPMODIS-GPP, GPPEC-LUE, and GPPEC showed the strong correlation between modelled GPP and measured GPP. The performance of RMSE varied at different

Discussion
The R 2 between GPP MODIS-GPP , GPP EC-LUE , and GPP EC showed the strong correlation between modelled GPP and measured GPP. The performance of RMSE varied at different sites. Compared with GPP MOD17 , GPP MODIS-GPP offered better GPP prediction at the US-Ton site and US-Bi2 site, and GPP EC-LUE offered better GPP prediction at the US-Ton site, US-Var site, and US-Rws site. At the US-Ton site, GPP MODIS-GPP and GPP EC-LUE obtained the best performance. Previous research showed that the savanna GPP estimation accuracy was smaller than other biome types because of the misclassification and sparse vegetation cover [50,74,80]. However, the woody savanna has a denser vegetation cover than the savanna [74,80,81]. This structure can reduce the influence of understory vegetation and improve GPP estimation accuracy [82]. Additionally, GPP MODIS-GPP and GPP EC-LUE also obtained a good performance at the US-WCr site and US-PFa site. This finding was consistent with previous studies [5,83]. For the US-WCr site, both GPP MODIS-GPP and GPP EC-LUE performed best in estimating GPP. The reasons for this can be summarized as follows. Firstly, DBF has obvious seasonal and phenological characteristics which can be monitored by satellites in real time [32]. Secondly, the vegetation structure of DBF is simple and the vegetation coverage varies greatly in different seasons, so the estimation of FAPAR is relatively simple and accurate [32]. Thirdly, DBF usually grows at mid and high latitudes with less cloudiness so that high-quality time-series images can be obtained [44]. For the US-PFa site, GPP MODIS-GPP and GPP EC-LUE represented the growth status of MF. However, compared with GPP MODIS-GPP , GPP EC-LUE was overestimated (RMSE = 3.23 gC/m 2 /d, bias = 2.07 gC/m 2 /d). The factor that caused this overestimation was the complex vegetation structure of the mixed forest, which had various ε max due to abundant plant species. Moreover, the leaf shape, tree height, and light energy absorption conditions are different among tree species. Therefore, a single ε max used in the input model causes certain errors [84,85]. The performance of GPP MODIS-GPP and GPP EC-LUE at two shrubland sites varied. For the US-Rls site, GPP MODIS-GPP and GPP EC-LUE both had a moderate correlation against GPP EC , with R 2 = 0.63 and 0.57, respectively. However, for the US-Rws site, the accuracy of GPP MODIS-GPP was the lowest, and the performance of GPP EC-LUE was better than that of GPP MODIS-GPP . The factor that caused this phenomenon was the estimation accuracy of GPP at shrubland sites, which varied in different areas [22,86]. The growth of shrubland is mainly affected by water, so more consideration should be given to water influence factors when estimating shrubland GPP [81,84]. Therefore, for shrubland in temperate and frigid zones, GPP estimation accuracy is higher, but in the tropics, affected by VPD and FAPAR, the accuracy of shrubland GPP estimation is lower. At the US-Var site, the performance of GPP EC-LUE was also better than GPP MODIS-GPP , and the accuracy of the GPP EC-LUE was the second highest (R 2 = 0.67). The accuracy of GPP MODIS-GPP (R 2 = 0.49) was consistent with that of GPP MOD17 (R 2 = 0.46). However, the performance of GPP MODIS-GPP at the US-GLE site was far better than that of GPP EC-LUE . For the US-GLE site, the accuracy of GPP MODIS-GPP was worse than the US-Ton site (R 2 = 0.72). For the US-Bi2 site, the relationship between GPP MODIS-GPP , GPP EC-LUE and GPP EC was moderate, and GPP MODIS-GPP and GPP EC-LUE obviously underestimated the cropland GPP (RMSE = 7.70 gC/m 2 /d, bias = −4.53 gC/m 2 /d and RMSE = 7.80 gC/m 2 /d, Bias = −4.15 gC/m 2 /d, respectively). This finding is aligned with previous research [50,87]. In crop ecosystems, different crops have different ε max and growth cycles, so a single ε max cannot represent well all types of crops [88,89]. Additionally, the rotation period varies for different crop types. Therefore, if the same input variables are used to estimate GPP for different crops, this can lead to a large deviation in GPP estimation [90].
The performance of GPP GR and GPP VI was diverse for different biome types. The correlation between GPP GR , GPP VI , and GPP EC was strongest at the US-WCr site, with R 2 = 0.84 and 0.81, respectively. Compared to the US-WCr site, the correlation between GPP GR , GPP VI and GPP EC was relatively weak at the US-GLE site. This result was consistent with previous research, where OTCI was highly correlated with tower GPP across deciduous forests and weakly correlated with tower GPP across evergreen forests [49,50]. The reason for this phenomenon could be that stress causes a decrease of the photosynthetic efficiency [91], the conical canopy structure, and the density of evergreen needleleaf trees. The consequent shadowing effect leads to the failure of OTCI for detecting subtle changes within the seasonal chlorophyll content of evergreen needleleaf forests [92,93]. Additionally, GPP GR and GPP VI showed moderate correlations with GPP EC at the US-PFa site. The factors that caused this can be considered from two aspects. Firstly, the vegetation structure of mixed forests is complex, and different vegetation types will affect the estimation of GPP. Secondly, mixed forest productivity will drop sharply when the temperature drops in a short period. However, the vegetation index cannot change significantly in this case, leading to its failure in capturing the impact of temperature change on GPP [47,94,95]. A strong correlation could also be observed between GPP GR , GPP VI , and GPP EC at the US-Bi2 site. Previous studies have found a similar relationship between OTCI or other chlorophyll indices and GPP [49,78,96]. Notably, GPP GR and GPP VI had the highest RMSE (4.21 gC/m 2 /d and 4.29 gC/m 2 /d, respectively) at the US-Bi2 site. The reason for this is that compared with GPP EC , GPP GR and GPP VI overestimated GPP in the nongrowing season and underestimated it in the growing season (Figure 7h). A previous study demonstrated that OTCI started to increase earlier than the time of corn sowing [49]. The high GPP GR and GPP VI obtained in the nongrowing season may represent fallow land, which is gradually colonized by weed species or affected by humidity changes, and not represent actual crop growth [49]. At the US-Var site, the correlation between GPP VI and GPP EC (R 2 = 0.66) was better than that between GPP GR and GPP EC (R 2 = 0.50). This difference may be related to the vertical and horizontal heterogeneity of grasslands [49]. C3 annual grasses dominated the US-Var site. Previous studies [49,97] showed that, compared with grasslands dominated by C4 species, C3 grasslands exhibited lower OTCI values during the peak growth period, explaining why GPP VI had a higher accuracy than GPP GR . However, GPP GR and GPP VI performed poorly in woody savannas and two shrubland sites. The probable reason for this was that the influence of the soil background on the low vegetation coverage at these sites led to the failure of OTCI in tracking GPP [50]. At the regional scale, GPP MODIS-GPP and GPP EC-LUE agreed quite well with GPP MOD17 . For the annual mean GPP, GPP MODIS-GPP and GPP EC-LUE were relatively high in the western coastal areas, where the biome types were mainly evergreen forests. The low GPP MODIS-GPP and GPP EC-LUE were obtained in grasslands. The finding aligned with a previous study [61], where evergreen forests exhibited strong photosynthesis and open shrublands were the least productive. Forest ecosystems have a relatively higher maximum daily GPP, and open shrublands have the lowest maximum daily GPP compared with other vegetation types. Notably, the maximum daily GPP across grasslands in the northeast was larger than in central areas, and the high sensitivity to soil moisture may be the cause of this phenomenon [61]. The inconsistency between the annual mean GPP and maximum daily GPP may be mostly due to the influence of temperature and rainfall on the length of the different growing seasons [61].
In addition, we further analyzed the uncertainty in simulating GPP using Sentinel-3 OLCI products and MERRA2 meteorology reanalysis data in this study. Several factors can influence the GPP estimation using the MODIS-GPP algorithm and EC-LUE model. First is the inconsistency of the spatial resolution between Sentinel-3 and meteorology reanalysis data. The errors may have been derived from the interpolation of MERRA2 meteorology reanalysis data from 0.5 • × 0.625 • to 300 m. Second is the coarse classification of vegetation types on land cover maps. Boston University's UMD classification scheme in the MOD12Q1 land cover dataset classifies croplands as one category and does not distinguish between C3 and C4 crops. C3 and C4 crops possess various photosynthetic pathways and light use efficiencies [88,89]. Because the C3/C4 mixing ratio for each cropland pixel is unknown, we simply used the single ε max and meteorological parameters to estimate GPP on a large scale. Therefore, GPP MODIS-GPP and GPP EC-LUE are underestimated or overestimated in croplands. Thirdly, the image data quality is also a significant factor influencing GPP estimation. Due to atmospheric contamination (i.e., clouds, aerosols), Sentinel-3 OLCI products have more missing data and low reliability during winter. Future work should be focused on spatial−temporal fusion techniques, which can further enhance the temporal and spatial continuity of Sentinel-3 OLCI products by fusing Sentinel-2 images with Sentinel-3 OLCI to achieve more accurate regional and large-scale GPP estimation [98,99].
Although our results showed that GPP MODIS-GPP , GPP EC-LUE , GPP GR and GPP VI did not always perform better than GPP EC or GPP MOD17 , it needs to be emphasized that the objective of our study was not to distinguish which model was superior for GPP estimation. Our study aimed to demonstrate the potential of integrating Sentinel-3 OLCI products and MERRA2 meteorology reanalysis data to estimate GPP at site and regional scales. Our results also demonstrated that the LUE model and VI-driven model were suitable for ESA's Sentinel-3 data. The Sentinel-3 Sea and Land Surface Temperature Radiometer (SLSTR) provides land surface temperature (LST) products with a higher temporal resolution. In future research, the integration of SLSTR LST products and OLCI FAPAR or meteorology reanalysis data in other models for estimating GPP should be conducted. For a wider application of Sentinel-3 OTCI products, future work requires more in situ data to explore the OTCI and GPP's relationship. Moreover, further investigations should focus on combining other remote sensing data or satellite-driven models with OTCIbased models. Previous studies have found that meteorological factors could improve VI-driven GPP estimation models [23]. Some existing satellite-driven models could make up for the insufficiency of GPP models based on physiologically driven spectral indices, such as MTCI [100,101]. Additionally, the Fluorescence Explorer satellite, the ESA's eighth Earth Explorer, is expected to be launched by 2022, which can offer important auxiliary data for using Sentinel-3 products to estimate GPP [33,102].

Conclusions
Within this research, we assessed the performance of two Sentinel-3 OLCI products (FAPAR and OTCI) combined with MERRA2 meteorology reanalysis data to estimate GPP through four GPP models (the MODIS GPP algorithm, EC-LUE model, GR model, and VI model) at the site and regional scales during 2017 and 2018. The following major conclusions can be drawn from the study: (1) The relationship between GPP MODIS-GPP , GPP EC-LUE and GPP EC is obvious for all sites. GPP MODIS-GPP exhibited the best performance with GPP EC at the US-Ton site (R 2 = 0.73), and the weakest performance for the US-Rws site (R 2 = 0.41). In addition, GPP MODIS-GPP and GPP EC-LUE were underestimated or overestimated at several sites, such as the US-Bi2 site, US-GLE site, US-WCr site and US-PFa site. Compared with GPP MOD17 , GPP MODIS-GPP was superior to GPP EC at the US-Ton site, US-Var site, and US-Bi2 site, and GPP EC-LUE was superior at the US-Ton site, and the US-Var site; (2) Compared with GPP EC , the performance of GPP GR and GPP VI varied across different biome types. Good performances were obtained across deciduous broadleaf forests, mixed forests, grasslands, and croplands, and low R 2 values were obtained across evergreen needleleaf forests, shrublands, and woody savannas; (3) At the regional scale, GPP MODIS-GPP and GPP EC-LUE both showed a very high spatial−temporal consistency with GPP MOD17 . The GPP value was high in the western coastal area and low in southwestern shrubland across California. The annual mean GPP value and maximum daily GPP exhibited a strong correlation with GPP MOD17 in broadleaf forests, mixed forests, grasslands and croplands, but a relatively weak correlation in evergreen needleleaf forests and shrublands.