A Comparison of OCO-2 SIF, MODIS GPP, and GOSIF Data from Gross Primary Production (GPP) Estimation and Seasonal Cycles in North America

Abstract

Remotely sensed products are of great significance to estimating global gross primary production (GPP), which helps to provide insight into climate change and the carbon cycle. Nowadays, there are three types of emerging remotely sensed products that can be used to estimate GPP, namely, MODIS GPP (Moderate Resolution Imaging Spectroradiometer GPP, MYD17A2H), OCO-2 SIF, and GOSIF. In this study, we evaluated the performances of three products for estimating GPP and compared with GPP of eddy covariance(EC) from the perspectives of a single tower (23 flux towers) and vegetation types (evergreen needleleaf forests, deciduous broadleaf forests, open shrublands, grasslands, closed shrublands, mixed forests, permeland wetlands, and croplands) in North America. The results revealed that sun-induced chlorophyll fluorescence (SIF) data and MODIS GPP data were highly correlated with the GPP of flux towers (GPP_EC). GOSIF and OCO-2 SIF products exhibit a higher accuracy in GPP estimation at the a single tower (GOSIF: R² = 0.13–0.88, p < 0.001; OCO-2 SIF: R² = 0.11–0.99, p < 0.001; MODIS GPP: R² = 0.15–0.79, p < 0.001). MODIS GPP demonstrates a high correlation with GPP_EC in terms of the vegetation type, but it underestimates the GPP by 1.157 to 3.884 gCm⁻²day⁻¹ for eight vegetation types. The seasonal cycles of GOSIF and MODIS GPP are consistent with that of GPP_EC for most vegetation types, in spite of an evident advanced seasonal cycle for grasslands and evergreen needleleaf forests. Moreover, the results show that the observation mode of OCO-2 has an evident impact on the accuracy of estimating GPP using OCO-2 SIF products. In general, compared with the other two datasets, the GOSIF dataset exhibits the best performance in estimating GPP, regardless of the extraction range. The long time period of MODIS GPP products can help in the monitoring of the growth trend of vegetation and the change trends of GPP.

1. Introduction

Gross primary production (GPP) is the total amount of carbon fixed by vegetation through photosynthesis, contributing the largest global carbon flux and driving ecosystem functions [1,2]. The carbon cycle of terrestrial ecosystems is complicated and the uncertainty of terrestrial ecosystem flux remains high, resulting in the so-called “mystery of missing carbon” [3,4]. This means that there is an imbalance between carbon sink and carbon sources [5]. GPP is a key index that can be used to estimate the carbon flux of terrestrial ecosystems [6,7]. Therefore, the accurate estimation of GPP is important for studies on the carbon cycle and climate change.

There are four methods commonly used to estimate GPP: (1) the eddy covariance (EC) technique [8], (2) the light use efficiency (LUE) model [9], (3) the process-based model [10], and (4) sun-induced chlorophyll fluorescence (SIF) [11,12,13,14]. The EC method is considered the most accurate of the four methods [15,16,17]. However, it requires the establishment of flux towers. At present, the distribution of flux towers is uneven and sparse, making it impossible to obtain GPP globally with a high accuracy. The LUE model is widely used to estimate GPP. The LUE approach can be expressed as the product of PAR (the photosynthetically active radiation received by the plant canopy), fPAR (the proportion of PAR absorbed by plant canopies), and LUE_P (light use efficiency). MODIS GPP products (e.g., MOD17A2, MYD17A2, MOD17A2H, and MYD17A2H) are estimated using the LUE model [18]. However, the LUE model depends on multiple types of ground parameters. Therefore, the accuracy of each parameter will affect the accuracy of the estimated GPP [19,20,21]. One of the process models commonly used is the soil-canopy observation of photosynthesis and energy (SCOPE) balance model. The SCOPE model is a one-dimensional vertical model. It is assumed that the canopy is homogeneous in the horizontal and vertical directions [22]. In addition, the photosynthesis, stomatal adjustment, and atmospheric disturbance module parameters need to be set with reference to the measured data. Setting a large number of parameters brings great uncertainty to the estimation of GPP. In recent years, SIF has brought new prospects for estimating GPP and monitoring vegetation growth states [11,12,23,24]. SIF is a by-product of photosynthesis [25]. Therefore, compared with the traditional vegetation index, SIF can be used to reflect the growth state of vegetation and helps estimate the GPP more accurately [26]. Many studies have investigated performances of estimating GPP using SIF products derived from greenhouse gases observing satellite (GOSAT) and global ozone monitoring experiment-2 (GOME-2) [27]. Because of the low spatial resolution of GOSAT and GOME-2 (GOME-2: 40 × 80 km; GOSAT: 10 km diameter), it is impossible to estimate GPP at an ecosystem scale [13,28]. With the launch of orbiting carbon observatory-2 (OCO-2), which provides SIF products with a fine spatial resolution (1.3 × 2.2 km), the spatial resolution of estimated GPP has been significantly improved [13]. However, the footprint of OCO-2 is very sparse and cannot match the flux tower well, which generates huge challenges in estimating the global GPP using OCO-2 SIF [19,23]. To solve this problem, Li and Xiao have developed a new SIF product based on OCO-2 SIF and meteorological data, named ‘GOSIF’, which has a time resolution of eight days and a spatial resolution of 0.05° [18].

The data of OCO-2 SIF, MODIS GPP, and GOSIF all have fine resolutions and have the potential to estimate GPP at an ecological scale. However, many results have demonstrated that MODIS GPP (MOD17A2) underestimates vegetation GPP, especially in Africa [29]. Other results have shown that MOD17A2 has a limited capacity for representing the temporal and spatial variation of annual GPP in croplands [8]. Furthermore, only a few studies have validated MODIS GPP (MYD17A2H) with in-situ GPP measurements. In addition, Li and Xiao only compared GOSIF data with OCO-2 coarse-resolution SIF data (1°) in seasonal cycles [18]. Consequently, there is an urgent need to evaluate those products comprehensively. Moreover, as introduced above, the three products are based on very different techniques. It is thus speculated that GPP estimates based on these products would demonstrate different characteristics. Therefore, it is necessary to conduct inter-comparisons to explore the performances of different products when estimating GPP at an ecological scale. In this study, we utilized in-situ measurements of EC flux towers as the truth value of GPP and established relationships between truth values of GPP and three products of SIF/GPP to evaluate their performances. Inter-comparison experiments were carried out from two perspectives, namely, the site and the vegetation type, in this study. The location of an EC flux site determines the geographic coordinates and meteorological parameters. Performances of satellite products are determined by the original signal quality, which is related to the observation geometry and the sun’s position. Further, meteorological parameters exert significant influences on generations of GPP [15]. This explains why it is important to evaluate the performance of satellite-observation-derived products in terms of the site. In addition, GPP and its seasonal cycle vary with the vegetation type. The response of SIF to GPP could also vary with the vegetation type. Hence, it is necessary to evaluate the performance of SIF-GPP models from the perspective of vegetation types.

Except for the above stated issues, there are another two topics to be discussed. Firstly, the extraction range of SIF data around the flux tower has always been a problem in correlating remotely sensed data and in-situ measurements. With a small extraction range, one often cannot obtain adequate SIF data to correlate with GPP. With a large extraction range, the spatial mismatch introduces additional errors, resulting in a reduced accuracy of GPP estimates [23]. Secondly, few studies have been dedicated to investigating the effects of different observation modes and bands of OCO-2 on the accuracy of estimating GPP.

The rest of this paper is organized as follows. Section 2 introduces the materials and methods. Section 3 introduces the primary results. Section 4 presents our discussion. Finally, the whole study is concluded in Section 5.

2. Materials and Methods

2.1. OCO-2 SIF Data and Flux Tower Data

OCO-2 was launched on July 2014, with the primary task of collecting space-based global measurements of atmospheric CO₂ column mixing ratios [30]. OCO-2 is equipped with three high-resolution spectrometers to measure reflected sunlight in the 0.76 μm O₂-A band and 1.61 μm and 2.06 μm CO₂ bands. Such a feature ensures the robust and accurate retrieval of SIF [11]. NASA developed the IMAP-DOAS algorithm to retrieve SIF. We have added the full width at half maximum (FWHM) in our paper, as follows: OCO-2 records high-resolution spectra in three bands: the O₂-A band (757–775 nm, FWHM = 0.042 nm), a weak CO₂ absorption band (1594–1627 nm, FWHM = 0.076 nm), and a strong CO₂ absorption band (2043–2087 nm, FWHM = 0.097 nm) [13]. OCO-2 has three observation modes: nadir, glint, and target. Typical OCO-2 measurements are collected alternately between nadir and glint mode with a repeat frequency of approximately 16 days. The instrument views the ground directly below the spacecraft in the nadir mode, tracks near the location with direct sunlight reflected in the glint mode, and collects a large number of measurements over calibration/validation sites in the target mode [31]. The viewing zenith angle (VZA) varies, depending on the observation mode applied, and has an impact on the SIF value [19]. Therefore, the observation mode could potentially impact the estimation of GPP. In this work, we used OCO-2 SIF data (OCO2_L2_Lite_SIF.8r). We selected all flux towers in North America, which have both SIF757 and SIF771 data (SIF at retrieval bands 757 and 771 nm) and GPP data spanning from September 2014 to 2018. To obtain more SIF data, we enlarged the area to within 25 km around each flux tower. In particular, the selection of EC sites was based on an assumption of threshold (60%) of International Geosphere-Biosphere Program (IGBP) classifications [32]. Hence, we calculated the maximum land cover of each flux tower. We chose an EC flux tower if the plant type observed at the flux site occupied at least 60% of the corresponding MODIS pixel area within 25 km [24]. However, OCO-2 SIF data of some vegetation types are too sparse, so we broadened this range (minimum is 46%), as shown in

Table 1. Site descriptions of the flux towers, including the site ID, vegetation type, longitude, latitude, year the data was used, maximum land cover (Max(LC)), and flux tower principal investigators (PIs) of each flux tower (ENF: evergreen needleleaf forests, DBF: deciduous broadleaf forests, OSH: open shrublands, GRA: grasslands, CSH: closed shrublands, MF: mixed forests, WET: permeland wetlands, CRO: croplands).

Site ID	Type	Lon (°)	Lat (°)	Year	Max (LC)	PIs
CA-SCC	ENF	−121.30	61.31	2014–2016	55	Manuel Helbig
CA-TP1	ENF	−80.56	42.66	2014–2017	75	M. Altaf Arain
CA-TP3	ENF	−80.35	42.71	2014–2017	73	M. Altaf Arain
CA-TP4	ENF	−80.36	42.71	2014–2017	73	M. Altaf Arain
CA-TPD	DBF	−80.56	42.64	2014–2015	73	M. Altaf Arain
US-EML	OSH	−149.25	63.88	2014–2018	46	Ted Schuur
US-KFS	GRA	−95.20	39.06	2014–2017	58	Nathaniel Brunsell
US-KLS	GRA	−97.59	38.77	2014–2017	59	Nathaniel Brunsell
US-Los	WET	−89.98	46.08	2014–2019	87	Ankur Desai
US-PFa	MFs	−90.27	45.95	2014–2019	83	Ankur Desai
US-PHM	WET	−70.83	42.74	2014–2018	47	Anne Giblin
US-Ro2	CRO	−93.09	44.73	2014–2016	63	John Baker
US-Ro4	GRA	−93.07	44.68	2014–2019	73	John Baker
US-Syv	MF	−89.35	46.24	2014–2018	89	Ankur Desai
US-Tw1	WET Wetlands	−121.65	38.11	2014–2018	52	Dennis Baldocchi
US-Tw3	CRO	−121.65	38.12	2014–2018	52	Dennis Baldocchi
US-Tw4	CRO Wetlands	−121.64	38.10	2014–2017	53	Dennis Baldocchi
US-Twt	CRO	−121.65	38.11	2014–2017	51	Dennis Baldocchi
US-UMd	DBF	−84.70	45.56	2014–2019	44	Christopher Gough
US-WCr	DBF	−90.08	45.81	2014–2018	78	Ankur Desai
US-ALQ	WET	−89.61	46.03	2015–2019	86	Brent Olson
US-GLE	ENF	−106.24	41.37	2014–2018	71	Bill Massman
US-RIs	CSH	−116.74	43.14	2014–2018	88	Gerald Flerchinger

. Finally, we chose 23 EC flux towers (Figure 1). The OCO-2 L2 data contains instantaneous SIF data and a daily correction factor. We used the daily correction factor to convert instantaneous SIF to daily SIF and removed daily SIF with values lower than 0.

The GPP values were provided by AmeriFlux, which is a network of PI (flux tower principal investigators)-managed sites measuring ecosystem CO₂, water, and energy fluxes in North, Central, and South America. Network data processing involves performing uniform data quality checks and transforming original measurements from individual sensors to ecologically or micro-meteorologically meaningful quantities. The EC data analysis procedure includes quality marks, filtering criteria, and gap-filling. These methods can reduce the bias and generate high-quality and standardized datasets. Among them, missing data are reported using –9999 as the replacing value. The L2 Flux data obtained from AmeriFlux include the net ecosystem exchange (NEE), total ecosystem respiration (Re), photosynthetically available radiation (PAR), wind speed, heat, water vapor, and temperature (T), with a time resolution of half an hour or an hour [33,34]. The NEE measurements are routinely partitioned into GPP and ecosystem respiration using a nigh-time partitioning method [35]. NEE is considered to be the primary source of observed interannual variability in atmospheric accumulation of carbon dioxide. Reis a process of oxidative decomposition of vegetation. For vegetation, NEE is governed by GPP and Re [36]. In this study, we used the flux tower GPP data spanning from 2014 to 2018 (

Table 2. Basic information of the OCO-2 sun-induced chlorophyll fluorescence (SIF), MODIS, and GOSIF datasets.

	OCO-2 SIF	MODIS	GOSIF
Temporal coverage	09/2014 to present	02/2000 to present	2000 to 2017
Spatial resolution	1.3 × 2.25 km2	500 m × 500 m	0.05°× 0.05°
Spatial sampling method	intermittent	continuous	continuous
Revisit period	16 days	8 days	8 days (and monthly, annual)

). There are eight biome types around the 23 flux towers: mixed forest (two sites), permanent wetlands (five sites), evergreen needleleaf forests (five sites), croplands (three sites), deciduous broadleaf forests (three sites), grasslands (three sites), open shrublands (one site), and closed shrublands (one site). For US-GLE and US-RIs, we used the night-time partitioning method to process the original GPP products [37]. The daytime ecosystem respiration Re_day was deduced from the relationship between Renight (equal to NEE exchange of the net ecosystem at night) and temperature. Then, GPP could be calculated by the daytime ecosystem exchange NEEday and Reday (Formula 1).

$$GPP={\mathrm{Re}}_{day}-NE{E}_{day}$$

(1)

2.2. GOSIF Data

The spatial coverage and temporal resolution of OCO-2 SIF products depend on the footprint of the satellite and the observation mode. Moreover, the sparse footprint of OCO-2 results in a poor match with the flux towers. The GOSIF dataset was developed to address this limitation. The GOSIF dataset is a global product of SIF, with a spatial resolution of 0.05° and a temporal resolution of 8 days spanning from 2000 to 2017 (http://globalecology.unh.edu/). The GOSIF dataset was developed by Li and Xiao [18]. For each grid cell, the enhanced vegetation index (EVI), land cover type (MCD12C1), air temperature, vapor pressure deficit (VPD), and photosynthetically available radiation (PAR) were used as the input for a cubist regression tree model. Among them, EVI was derived from the MODIS bidirectional reflectance distribution function (BRDF)-corrected reflectance product MCD43C4. Daily PAR, VPD, and air temperature were obtained from MERRA-2 [18]. The cubist regression tree model is an extension of Quinlan’s M5 model tree. Cubist creates rule-based predictive models from training samples, and each rule has an associated multivariate linear regression model [14]. Their modeled SIF showed a good correlation with the GPP of EC flux sites (R² = 0.73, p < 0.001). We extracted GOSIF data for 0.1 degree surrounding each flux tower during the period from 2014 to 2017. Further, we wanted to explore the effects of the radius on estimating GPP using GOSIF. We thus extracted GOSIF data within 25 km of each flux tower. Finally, we used linear interpolation to interpolate eight days to one day.

2.3. MODIS Data

MODIS is a sensor onboard Terra and Aqua satellites [38]. In this study, we used the MYD17A2H dataset from 2014 to 2017. MYD17A2H is a Level-4 standard product, which includes Gpp_500 m, PsnNet_500 m, and Psn_QC_500 m. We used the Gpp_500 m dataset, which has a spatial resolution of 500 m and temporal resolution of eight days. ‘Gpp_500m’ data is an eight-day summation GPP. To translate the numerical value of MODIS GPP into the sequestered carbon value, the appropriate scale factor (=0.0001) for the MYD17A2H data was used as recommended by MYD17 users guide [39]. We applied a scale factor of 0.0001 to each MYD17A2H pixel to obtain GPP products in units of kgCm⁻²day⁻¹. We extracted the MYD17A2H data within 5 and 25 km of each EC flux tower and averaged all valid data as the daily MODIS GPP. The daily GPP was interpolated from the eight-day GPP dataset using a linear-weighted method. Table 2 shows the spatial and temporal resolution and coverage of the MODIS, OCO-2 SIF, and GOSIF data.

2.4. Methodology

Firstly, we established linear relationships between GPP_EC and three remotely sensed products, namely OCO-2 SIF products, MYD17A2H GPP (MYD), and GOSIF datasets at a single tower. As stated in Section 2.1, the spectra of two bands can be used to retrieve SIF products. Previous studies have indicated that SIF757 products are superior to SIF771 products in estimating GPP [13]. Consequently, we utilized SIF757 products to correlate with GPP_EC at different flux towers at a daily timescale. Before correlating GPP_EC with remotely sensed products, it is necessary to ensure their spatial and temporal consistencies. The date of available OCO-2 SIF products was used to screen daily GPP_EC measurements. Eight-day MYD and GOSIF data were interpolated into daily data using linear interpolation. Then, temporal matching was conducted with respect to the date. Regarding the spatial registration, it is necessary to use an extraction range to determine which data can be included in experiments. Remotely sensed products located in cycles, of which the centers are flux towers and the radius is the extraction range, were screened to correlate with GPP_EC. Theoretically, a unified extraction range should be applied to different remotely sensed products to ensure fair inter-comparison experiments. However, spatial resolutions of OCO-2 SIF products, MYD17A2H GPP (MYD), and GOSIF datasets are different. The pixel heterogeneity of GOSIF and MYD data with finer spatial resolution will increase if the same extraction range is used. The amount of three kinds of data will directly affect the evaluation results of GPP. To ensure there are adequate data to establish robust relationships between GPP_EC and different remotely sensed products, we used different extraction ranges for three remotely sensed products. Besides, additional experiments using a unified extraction range are presented in the discussion to explore the effect of the extraction range. Owing to the sampling strategy and the retrieving method of OCO-2, there are many gaps and blanks in its SIF products. We thus utilized a large extraction range (25 km) to extract OCO-2 SIF products of nadir glint and target mode. Both MYD17A2H GPP (MYD) and GOSIF provide areal products with a full coverage over our study area. The spatial resolutions of MYD and GOSIF are 500 m and 0.05°(~5 km), respectively. Therefore, extraction ranges were set to 5 km and 0.1° for MYD and GOSIF, respectively. Finally, linear relationships between GPP_EC and three remotely sensed products were then established and were named SIF757- GPP_EC, MYD-GPP_EC, and GOSIF-GPP_EC. We used the determination coefficient (R²) to evaluate the relationships between SIF757-GPP_EC, MYD-GPP_EC, and GOSIF-GPP_EC. The closer R² is to 1, the better the ability of remotely sensed products to estimate GPP_EC is.

Moreover, we used SIF757-GPP_EC, GOSIF-GPP_EC, and MYD-GPP_EC models to generate estimates of GPP for each vegetation type, abbreviated as GPP_SIF, GPP_GOSIF, and GPP_MYD, respectively. Finally, we calculated the root mean square error (RMSE) of GPP_SIF, GPP_GOSIF, and GPP_MYD for different vegetation types. For each vegetation type, we analyzed the consistency of the seasonal cycle of the three GPP estimates and GPP_EC from 2015 to 2017. Afterwards, we utilized OCO-2 SIF products of different bands (SIF757 and SIF771) and observation modes (nadir, glint, and target) to establish SIF-GPP models, in order to explore the effects of bands and observation modes on performances of modeling GPP using OCO-2 SIF products. In addition, the explanation of abbreviations and acronyms in this paper are summarized in

Table 3. Explanation of Abbreviations and Acronyms.

Abbreviations and Acronyms	Explanation
GPP	Gross primary production
SIF	Sun-induced chlorophyll fluorescence
SIF757	OCO-2 SIF at 757nm retrieval band
SIF771	OCO-2 SIF at 771nm retrieval band
MYD	MYD17A2H/MODIS GPP
GPPEC	GPP of EC tower
GPPGOSIF	GPP estimated using GOSIF-GPPEC model
GPPSIF	GPP estimated using SIF757-GPPEC model
GPPMYD	GPP estimated using MYD-GPPEC model
Re	Total ecosystem respiration
NEE	Net ecosystem exchange
EVI	Enhanced vegetation index
PAR	The photosynthetically active radiation received by the plant canopy
fPAR	The proportion of PAR absorbed by plant canopies
LUEP	Light use efficiency
CRO	Croplands
CSH	Closed shrublands
DBF	Deciduous broadleaf forests
ENF	Evergreen needleleaf forests
GRA	Grasslands
MF	Mixed forests
OSH	Open shrublands
NDWI	Normalized difference water index
WET	Permeland wetlands

.

3. Results

3.1. Correlating MYD, GOSIF, and SIF757 with GPP_EC at a Single Tower

Table 4. Comparisons of correlation coefficients between the MODIS, GOSIF, SIF757, and EC gross primary production (GPP) at a single tower (parameter: R²).

Site ID	R2
	MYD-GPPEC	GOSIF-GPPEC	SIF757-GPPEC
Mixed Forest
US-PFa	0.57	0.79	0.99
US-Syv	0.62	0.86	0.67
Permanent Wetlands
US-Los	0.58	0.88	0.66
US-PHM	0.62	0.81	0.70
US-Tw1	0.37	0.70	0.11
US-Tw4	0.45	0.71	0.08
US-ALQ	–	0.13	0.41
Evergreen Needleleaf Forests
CA-SCC	0.74	0.87	0.74
CA-TP1	0.57	0.57	0.08
CA-TP3	0.64	0.64	0.56
CA-TP4	0.76	0.78	0.76
US-GLE	–	0.82	0.52
Croplands
US-Twt	0.31	0.52	0.26
US-Ro2	0.22	0.15	0.43
US-Tw3	–	0.24	0.18
Deciduous Broadleaf Forests
CA-TPD	0.79	0.88	0.49
US-UMd	0.42	0.76	0.52
US-WCr	0.62	0.82	0.53
Grasslands
US-KFS	0.30	0.37	0.51
US-KLS	0.15	0.16	0.56
US-Ro4	0.62	0.85	0.91
OSH
US-EML	0.69	0.70	0.69
CSH
US-RIs	0.70	0.73	0.36

illustrates the R² of MYD-GPP_EC, GOSIF-GPP_EC, and SIF757- GPP_EC models. On the whole, GOSIF (R² = 0.15–0.88, p < 0.001) and SIF757 (R² = 0.08–0.99, p < 0.001) were more strongly correlated with GPP_EC than MYD (R² = 0.15–0.79, p < 0.001). GOSIF demonstrates the strongest ability to estimate GPP. It can be concluded that SIF was better at estimating GPP than the LUE model at an ecological scale for a single flux tower. The low correlation of the MYD-GPP_EC model is attributed to the fact that MODIS uses the LUE model to calculate GPP. LUE_P is the critical input parameter of the LUE model, which is susceptible to the plant growth status and ecological environments. Therefore, theoretically, the light use efficiency varies with the environments and growth stages of vegetation. However, LUE_P is generally determined as a fixed parameter in the LUE model, according to specific vegetation types, resulting in large uncertainties of estimated GPP [20]. R² is missing in sites US-ALQ, US-GLE, and US-Tw3 because there is insufficient available MODIS GPP data to build linear relationships. Poor correlations of SIF757-GPP_EC at US-Tw1 and US-Tw4 are due to small amounts of SIF757 data. The footprints of OCO-2 could not match the flux towers well, even if the sampling range was extended to 25 km, resulting in a lack of high-quality data to establish a linear model at the US-Tw1 and US-Tw4 sites. Compared with OCO-2 SIF, the GOSIF provides a dataset covering the globe, thus avoiding the above problems.

3.2. Relationships of SIF757, GOSIF, and MYD with GPP_EC for Different Vegetation Types

We used SIF757-GPP_EC, GOSIF-GPP_EC, and MYD-GPP_EC models to estimate the GPP (GPP_SIF, GPP_GOSIF, and GPP_MYD, respectively) for each vegetation type and then established linear correlation between GPP_SIF, GPP_GOSIF, GPP_MYD, and GPP_EC, as shown in Figure 2. On the whole, GPP_GOSIF (R = 0.62–0.95, p < 0.001) was most strongly correlated with GPP_EC, compared with GPP_MYD (R = 0.67–0.92, p < 0.001) and GPP_SIF (R = 0.28–0.86, p < 0.001). The amount of SIF757 is small, which directly leads to a poor performance of the GPP_SIF-GPP_EC model, especially for (g) open shrublands (OSH), for which only seven observations are available. Moreover, the insufficient accuracy of a small portion of OCO-2 SIF757 products is also responsible for the performance of the GPP_SIF-GPP_EC model. Referring to Table 3, an interesting phenomenon is worth noting, in that GPP_MYD exhibits a poor correlation with GPP_EC in terms of a single flux tower. However, a good correlation between GPP_MYD and GPP_EC was witnessed in terms of the vegetation type, as shown in Figure 3.

Figure 3 shows the annual average GPP_EC, GPP_MYD, and GPP_GOSIF for different vegetation types. It can be concluded that the annual average GPP of evergreen needleleaf forests (ENF) is the highest, and that of OSH is the lowest. Compared with the annual average GPP_MYD, the annual average GPP_GOSIF is closer to the annual average GPP_EC, which suggests that the accuracy of GOSIF-estimated GPP is higher than that of MODIS GPP. It is worth noting that GOSIF overestimates GPP for deciduous broadleaf forests (DBF) and grasslands (GRA), implying that there is still necessary work to be done before estimating the GPP of those vegetation types using GOSIF datasets. On the contrary, MODIS GPP underestimated the GPP of eight vegetation types. In particular, for ENF and closed shrublands (CSH), more than half of the GPP values were underestimated by MODIS GPP. MODIS GPP has a long time series and can thus be used to observe changes in the vegetation coverage of long-term series. However, it is not suitable for quantitatively calculating carbon fluxes. The reason GPP_SIF is not shown in Figure 3 is that the data volume of SIF757 is too small. For most flux towers, the data volume of SIF757 is about 20 from 2014 to 2018. Because OCO-2 has sparse overpasses globally and short time series [19], we cannot calculate the annual average of GPP_SIF by using SIF757-GPP_EC model.

Table 5. Statistical comparisons of the GPP_GOSIF-GPP_EC, GPP_SIF-GPP_EC, and GPP_MYD-GPP_EC models from 2015 to 2017 for different vegetation types. Parameters: RMSE. The units are gCm⁻²day⁻¹.

Type	GPPGOSIF	GPPSIF	GPPMYD
CRO	3.121	4.578	3.326
GRA	3.154	3.402	3.253
MF	2.160	2.437	1.856
ENF	3.695	3.873	3.884
DBF	2.566	3.262	2.585
CSH	1.187	1.938	1.157
OSH	2.346	2.609	2.138
WET	3.981	2.399	1.245

illustrates the RMSE of the GPP_SIF-GPP_EC model, GPP_GOSIF-GPP_EC model, and GPP_MYD-GPP_EC model. For the eight vegetation types, the RMSE of the GPP_SIF-GPP_EC model is the largest. This can be attributed to the fact there are not sufficient OCO-2 SIF products to establish robust models. The RMSE of the GPP_MYD-GPP_EC model is lower than that of the GPP_GOSIF-GPP_EC model for mixed forests (MF), CSH, OSH, and permeland wetlands (WET), indicating that the accuracy of MODIS GPP estimates can be improved after averaging the MYD17A2H data of multiple study areas, according to the vegetation type.

3.3. Comparisons of Seasonal Cycles of Different GPP Estimates

We compared the seasonal trends of GPP_EC and that of SIF757 for eight vegetation types, as shown in Figure 4. Overall, SIF757 demonstrates a single peak state within one year, which is consistent with the growth pattern of vegetation. However, Figure 4 also shows that the amount of SIF757 data for each vegetation type from 2015 to 2017 is small. This is attributed to the fact that the OCO-2 footprint is sparse in space and does not match well with the flux tower. Therefore, it would be difficult to estimate the global distribution of GPP with a high accuracy using SIF757. Due to inadequate amounts of data, it is impossible to judge whether the seasonal cycle of the eight vegetation types is advanced or delayed. In Figure 4a, croplands (CRO), and in Figure 4d, ENF, SIF757 do not demonstrate a clear single peak.

The seasonal cycle of GPP_EC and MODIS GPP (MYD17A2H) for different vegetation types is shown in Figure 5. On the whole, the trends of MODIS GPP and GPP_EC are similar for eight vegetation types. However, MODIS GPP significantly underestimated GPP for (a) CRO, (b) CSH, (c) DBF, (d) ENF, (e) GRA, (f) MF, and (h) WET. This is because the LUE_P of these eight vegetation types is set too low in the LUE model, resulting in the low MODIS GPP value of the final estimate. Only for (g) OSH was the MODIS GPP value close to GPP_EC in 2016 and 2017. It can also be seen in Figure 5 that the MODIS GPP fluctuates greatly after reaching its peak in summer, especially for (c) CSH, (f) MF, and (g) OSH. MODIS GPP reaches its peak in February 2015 for (e) GRA and (h) WET, suggesting erroneous seasonal cycles of MODIS GPP products. However, MODIS GPP has a longer time series than OCO-2 SIF. MYD17A2H can thus be used to monitor changes in the vegetation cover.

We also evaluated seasonal cycles of GPP_EC and GOSIF data for each vegetation type from 2015 to 2017, as shown in Figure 6. The seasonal cycle of GOSIF was the most consistent with GPP_EC in all three datasets, verifying the conclusion that the linear correlation between GOSIF and GPP_EC is the highest. Especially for (c) DBF, (f) MF, and (g) OSH, the seasonal trend of GOSIF is exactly consistent with that of GPP_EC. However, the advanced seasonal cycle of (e) GRA and (h) WET and lagged seasonal cycle of (b) CSH suggest significant errors in GOSIF data for those vegetation types. However, there may be errors in the GPP_EC data of (a) CRO in 2017, which led to the low coincidence of seasonal cycles of GOSIF.

3.4. Effect of Observation Modes and Bands on Correlating OCO-2 SIF with GPP_EC

In order to explore the relationship between OCO-2 SIF and GPP_EC in different modes and bands, we analyzed the relationship between SIF and GPP_EC.

Table 6. Comparisons of correlation coefficients for SIF and EC GPP for the different observation modes and bands (parameter: R²).

Site ID	SIF757			SIF771
	Nadir	Glint	Target	Nadir	Glint	Target
Mixed Forest
US-PFa	0.99	0.32	0.67	0.95	0.48	0.59
US-Syv	0.67	0.57	–	0.71	0.49	–
Permanent Wetlands
US-Los	–	0.66	0.42	–	0.77	0.32
US-PHM	0.70		–	0.55	–	–
US-Tw1	0.11	0.01	–	0.38	0.01	–
US-Tw4	0.08	0.01	–	0.41	0.01	–
US-ALQ		0.41		–	0.32
Evergreen Needleleaf Forests
CA-SCC	0.74	–	–	0.48	–	–
CA-TP1	–	0.08	–	–	0.05	–
CA-TP3	0.56	0.20	–	0.44	0.20	–
CA-TP4	0.54	0.76	–	0.46	0.74	–
US-GLE	0.48	0.52		0.47	0.50
Croplands
US-Twt	0.26	–	–	0.12	–	–
US-Ro2	0.43	–	–	0.35	–	–
US-Tw3	0.18	0.06	–	0.35	0.07	–
Deciduous Broadleaf Forests
CA-TPD	–	0.49	–	–	0.48	–
US-UMd		0.52	–		0.16
US-WCr	–	0.49	0.53	–	0.22	0.50
Grasslands
US-KFS	–	0.51	–	–	0.47	–
US-KLS	0.38	–	–	0.56	–	–
US-Ro4	0.85	–	–	0.91	–	–
OSH
US-EML	0.69	0.15	–	0.41	0.29	–
CSH
US-RIs	–	0.36	–	–	0.40	–

shows the parameters of SIF757-GPP_EC and SIF771-GPP_EC linear models for the different OCO-2 observation modes. It can be concluded that the amount of SIF data in nadir mode and glint mode is greater than that of the target mode in North America [19]. Thus, because nadir and glint mode were observed every 16 days, more SIF data were obtained. Meanwhile, glint mode had the least number of observations. Moreover, the spatial resolution of data in nadir mode is the highest, so the correlation between SIF data in nadir mode and GPPEC is more efficient [32]. For different vegetation types, the R² of relationships between SIF and GPP varies with the observation mode. The amount of OCO-2 SIF data in nadir mode is the largest and there is a strong linear correlation with GPP_EC. Because nadir and glint mode were observed every 16 days, more SIF data were obtained. In general, SIF757 and SIF771 data showed a strong correlation with GPP in all observation modes [11,13,24]. The linear correlation between SIF757 and GPP is stronger than that between SIF771and GPP. It could be responsible for the fact that 757 nm is closer to the peak of the SIF spectrum than SIF771. For WET, the correlation between SIF and GPP of US-Tw1 and US-Tw4 of the glint mode is low, mainly due to the small amounts of SIF in glint mode.

3.5. The Impact of the Extraction Range on Relationships between the Remotely Sensed Product and GPP_EC

In this study, we extracted 0.1-degree GOSIF data around each flux tower, 5 km MODIS GPP data, and 25 km SIF757 data, to establish relationships with GPP_EC. Then, we compared the performance of GPP estimates from SIF757, MODIS GPP, and GOSIF. However, the different ranges of each data extraction result in a different proportion of heterogeneity between pixels, which may affect the results obtained when comparing the three data estimates’ GPP accuracy [19]. Therefore, on the one hand, in order to compare the accuracy of SIF757, MODIS GPP, and GOSIF, GPP was estimated in the same research range. On the other hand, in order to explore whether the extracted range would affect the relationship of SIF757, MODIS GPP, and GOSIF, we also extracted GOSIF and MODIS GPP data within 25 km around each flux tower. Then, the relationship of SIF757, MODIS GPP, and GOSIF in the same extraction range was analyzed.

Table 7. Correlation coefficients (R²) of GPP_GOSIF-GPP_EC, GPP_SIF-GPP_EC, and GPP_MYD-GPP_EC models with an extraction range of 25 km.

Type	GOSIF	SIF757	MYD
CRO	0.37	0.24	0.56
GRA	0.72	0.67	0.68
MF	0.87	0.71	0.82
ENF	0.74	0.37	0.71
DBF	0.91	0.60	0.87
CSH	0.69	0.36	0.58
OSH	0.68	0.69	0.39
WET	0.88	0.63	0.88

shows the linear relationship between GOSIF, SIF757, and MODIS and GPP_EC after extending the extraction range to 25 km. In general, the linear relationship between GOSIF and GPP_EC is the best, followed by MODIS and OCO-2 SIF. However, the correlation between GOSIF data and GPP_EC was low of GRO. This may be due to the CRO’s GPP_EC data in 2017 with errors, which we can see from Figure 6. On the whole, those results confirm that GOSIF yields the best performance in estimating GPP, regardless of the extraction range and vegetation types. However, the R² of three models slightly decreased to different extents after extending the extraction range, highlighting the significance of stringent temporal-spatial registration between remotely sensed products and ground-based in-situ measurements. In particular, for different vegetation types, the R² of SIF757 is slightly lower than MODIS GPP with GPP_EC. For a single flux tower, the R² of SIF757 is higher than that of MODIS GPP (Table 4). The main reason for this is that the amount of data for SIF757 is small. Moreover, the growth environment of different longitudes and latitudes is inconsistent, which leads to differences in the value of SIF757 of the same vegetation type. The average SIF757 data of different flux towers of the same vegetation type may contain errors.

4. Discussion

There is a significant linear relationship between SIF and GPP. We demonstrated that SIF757 is better than SIF771 in estimating GPP, which is consistent with many previous research results [12,23]. This is because 771 nm falls farther away from peak emission on the SIF spectrum and OCO-2 SIF757 is considered to have a higher level of retrieval precision than SIF771 [12]. Our results demonstrate that the linear relationship between SIF in nadir mode and GPP is slightly better than that for SIF in glint and target mode. We believe that SIF data in nadir mode has the highest spatial resolution. In this study, we can conclude that the ability of GOSIF and SIF757 to evaluate the GPP of a single tower is not significantly different. However, for different vegetation types, GOSIF was better at estimating GPP than SIF757 [19]. One of the main reasons for this is that the footprint of the OCO-2 is sparse. OCO-2 has a fine spatial resolution, which provides us with an unprecedented opportunity to explore the relationship between OCO-2 SIF and GPP of different vegetation types at an ecosystem scale. However, the amount of SIF data of OCO-2 around each flux tower is small. This causes the following problems. First, at present, the main research on SIF aims to estimate GPP by SIF. Many researches have demonstrated that there is a strong linear relationship between SIF and GPP compared with the traditional vegetation index [40,41]. There is a clear and simple mechanism for estimating GPP using SIF [42]. However, SIF estimates of the GPP have been limited by the satellite spatial resolution. The GOME-2 data covers the globe and has a long time series, but the noise of the data leads to a poor GOME-2 SIF data quality [43]. OCO-2 data has a high spatial resolution, but the data volume is small and the time series is short. Currently, there are no satellites dedicated to SIF observation, which will limit our ability to use SIF to estimate GPP. Therefore, we also cannot accurately estimate global carbon budgets. Secondly, many researchers have shown that the significance of the linear relationship between SIF and GPP for different vegetation types is different [13,44]. OCO-2 has a small amount of data, which affects our ability to explore SIF to estimate GPP for different vegetation types. Thirdly, in this study, the seasonal cycle of SIF757 presents a ‘single peak’. However, because the data volume of SIF is too small, we cannot compare whether the seasonal cycles of OCO-2 SIF and GPP_EC from 2015 to 2017 are advanced or delayed. Therefore, if we use OCO-2 SIF to explore the phenology of vegetation, it will be a great challenge [24]. Moreover, most flux tower data are not publicly available. Hence, we validated SIF data by using flux tower data on a global scale, which produces limitations. Above all, scientists using OCO-2 SIF data to estimate global GPP will also face huge challenges.

Although OCO-2 has a high temporal and spatial resolution, the dataset still cannot match flux towers well on an ecological scale [19,23]. On one hand, in order to get more OCO-2 SIF data in this study, we had to expand our study area to 25 km. However, the expansion of the research scope led to an increase in mixed pixels, and the increase in pixel heterogeneity may have had an impact on our research results. On the other hand, measuring GPP by EC is complicated with many factors influence the final measurement. For example, GPP values may be influenced by the temperature and tower height. GPP values obtained by EC technology vary in quality. Therefore, when establishing the linear relationship of GOSIF, SIF, and MODIS with GPP_EC, inaccurate GPP_EC values can also cause errors. For example, R² of GOSIF-GPP_EC of US-Ro2 and US-ALQ is low. The reason of our analysis is that the accuracy of GPP_EC value is low. In this study, we compared the SIF757-GPP_EC, GOSIF-GPP_EC, and MYD-GPP_EC within 25 km around each flux tower. However, SIF757, GOSIF, and MYD have different spatial resolutions, which means that even within the same extraction range, SIF757, GOSIF, and MYD obtain different amounts of data and the heterogeneity of pixels is different [23]. These factors also affect the accuracy of the results. Moreover, except for the uncertainty caused by the scale problem, the land cover and the quality of SIF757 are also sources of uncertainty of GPP_SIF.

The GOSIF data have a higher temporal and spatial resolution than the OCO-2 SIF dataset. Our results show that the GOSIF-GPP_EC model is an effective method for estimating GPP. However, the seasonal cycle of GOSIF of GRA and ENF occurred earlier than GPP_EC. This phenomenon was also observed in the MODIS GPP dataset. To calculate the GOSIF data, a large number of EVI and PAR values from the MODIS dataset were used. In addition, the cubist regression tree model’s input parameters from the MERRA-2 dataset have a coarse spatial resolution (0.5° × 0.625°), and the data were resampled to 0.05° [18]. The resampling process may cause discontinuities between the datasets, creating jagged or smoothing effects, and thus leading to errors in the predictions of SIF [45]. We can conclude that inconsistent seasonal cycles of GRA and ENF between GPP estimates and GPP_EC could be attributed to inaccurate data from the input model for MODIS or from the MERRA-2 data resampling. Therefore, we call for necessary improvements in the quality of input data for producing GOSIF datasets. When Li and Xiao developed GOSIF data, in order to improve the calculation efficiency of the model and to reduce the uncertainty of the data, only LAI, fPAR, NDVI, and the normalized difference water index (NDWI) were input into the cubist regression tree model [18]. However, change of land cover type has a great influence on the SIF value. The value of SIF varies with different vegetation types. Therefore, we suggest that land cover types should be added to the model input data to obtain more accurate GOSIF data. Based on our results, we believe that GOSIF data has great potential to estimate the GPP of different vegetation types globally. GOSIF data has a spatial resolution of 0.05 degrees and temporal resolution of eight days, which has important indicative significance for the study of the global carbon cycle using SIF data. At the same time, GOSIF data covering the whole world is of great significance to exploring the function of SIF in drought monitoring, vegetation phenology, and stress monitoring. MODIS GPP products are based on the LUE model. Many inaccurate ground parameters are input to the LUE model, resulting in errors of its GPP products [9,46]. In addition, few land cover classifications are incorporated into the MYD17A2H algorithm. At a 500 m resolution, the MODIS GPP dataset can cause the overlapping of land cover types rather than classifying individual pixels [30]. In this study, we averaged all MYD17A2H pixels within 25 km. To some extent, it can reduce errors and smooth the MYD17A2H data. However, the average MODIS GPP may contain many land types, which could introduce additional errors. Many studies have shown that MODIS data (MOD17A2) underestimate the GPP [47]. Our study also found that MYD17A2H data underestimate the GPP. The accuracy of MODIS GPP is largely determined by the accuracy of the LUE model and input parameters [48,49]. Therefore, we suggest improving the LUE model or controlling the quality of input parameters so that they have a higher precision.

In recent years, more and more researchers are exploring the wider application of SIF data of different vegetation types on the ecosystem scale [19,50,51]. It is a simple and effective method to use data of EC flux tower to verify the accuracy of their simulation results. However, EC flux towers are not evenly distributed on the globe, and data of most flux towers are not available to the public. Therefore, verification and use of SIF data is limited. What’s more, the problem that the range of the observation of the flux tower does not match the extraction range of the SIF data seriously affected the accuracy of using SIF to estimate GPP. Ideally, we believe that the type of vegetation around the flux towers is homogeneous and single. Yet that is not the case. In general, we believe that the smaller the extraction range of SIF data around the flux tower, the higher the accuracy of the result may be. Li et al. [23] extracted 2 × 2 km SIF data around tow flux towers and GPP data proved that OCO-2 SIF has a significant linear relationship with GPP. Lu et al. [50] extracted 10 km SIF data around flux towers found that a strong correlation between GPP and SIF. Wood et al. [26] and Verma et al. [15] made the evaluation of OCO-2 SIF using tower GPP in corn and grassland ecosystems within a radius of 10 to 25 km. Therefore, the extraction range of SIF varies for different research purposes. However, if the extraction range of SIF data around the flux towers is too large, the more vegetation types may be included and the higher pixel heterogeneity of SIF data will be. This will eventually affect the accuracy of results. Therefore, we calculated max(LC) to ensure the accuracy of the this study. With the launch of the OCO-3 (orbiting carbon observatory 3) satellite and the TROPOMI (TROPOspheric Monitoring Instrument) satellite with finer resolution [52]. An increasing amount of flux tower data are also public. The mismatch between extraction range of satellite’s SIF data and observation range of EC flux towers may be better solved in the future.

5. Conclusions

In this study, we evaluated the performance of the three remotely sensed products MYD17A2H, GOSIF, and OCO-2 SIF when estimating GPP. We modeled and compared remotely sensed product-derived GPP estimates with GPP_EC data from 23 flux towers and eight vegetation types (evergreen needleleaf forests, deciduous broadleaf forests, open shrublands, grasslands, closed shrublands, mixed forests, permeland wetlands, and croplands). The GOSIF dataset yields the best results in the three datasets, though it leans towards overestimating GPP for deciduous broadleaf forests and grasslands. The MYD17A2H dataset leans towards underestimating GPP by 1.157 to 3.884 gCm⁻²day⁻¹ for all vegetation types, especially for croplands, grasslands, and evergreen needleleaf forests. The insufficient amounts of OCO-2 SIF make it difficult to establish good relationships between OCO-2 SIF and EC GPP for different vegetation types. In terms of the seasonal cycle, GOSIF and MODIS GPP demonstrate a good consistency with GPP_EC. Seasonal cycles of MYD17A2H and GOSIF show advanced trends compared with GPP_EC for grassland and evergreen needleleaf forest. Though the amount of OCO-2 SIF data is small, its seasonal cycle basically follows the trend of GPP_EC. Observation modes of OCO-2 SIF exert an evident influence on estimating GPP. After applying a unified extraction range for the three remotely sensed products in the spatial registration process, the results still show that the ability of GOSIF data to estimate GPP is better than that of MYD17A2H and OCO-2 SIF. This also suggests better performances of SIF in estimating GPP than those of the LUE model. The MYD17A2H product has a long time series and is thus suitable for monitoring the seasonal cycle change of vegetation, as well as total amounts of GPP. The GOSIF dataset has a large amount of data with a high temporal and spatial resolution, which is useful for estimating global GPP and monitoring the vegetation growth status. The data volume of OCO-2 SIF is relatively small, but OCO-2 provides the first opportunity to estimate GPP at an ecological scale. With the development of spaceborne high-resolution spectrometer techniques, similar satellite-derived SIF products would demonstrate very promising prospects in monitoring GPP at an ecological scale.