Evaluation of the Latest MODIS GPP Products across Multiple Biomes Using Global Eddy Covariance Flux Data

The latest MODIS GPP (gross primary productivity) product, MOD17A2H, has great advantages over the previous version, MOD17A2, because the resolution increased from 1000 m to 500 m. In this study, MOD17A2H GPP was assessed using the latest eddy covariance (EC) flux data (FLUXNET2015 Dataset) at eighteen sites in six ecosystems across the globe. The sensitivity of MOD17A2H GPP to the meteorology dataset and the fractional photosyntheticallyactive radiation (FPAR) product was explored by introducing site meteorology observations and improved Global Land Surface Satellite (GLASS) Leaf Area Index (LAI) products. The results showed that MOD17A2H GPP underestimated flux-derived GPP at most sites. Its performance in estimating annual GPP was poor (R2 = 0.62) and even worse over eight days (R2 = 0.52). For the MOD17A2H algorithm, replacing the reanalysis meteorological datasets with the site meteorological measurements failed to improve the estimation accuracies. However, great improvements in estimating the site-based GPP were gained by replacing MODIS FPAR with GLASS FPAR. This indicated that in the existing MOD17A2H product, the errors were originated more from FPAR than the meteorological data. We further examined the potential error contributions from land cover classification and maximum light use efficiency (εmax). It was found that the current land cover classification scheme exhibited frequent misclassification errors. Moreover, the εmax value assigned in MOD17A2H was much smaller than the inferred εmax value. Therefore, the qualities of FPAR and land cover classification datasets should be upgraded, and the εmax value needs to be adjusted to provide more accurate GPP estimates using MOD17A2H for global ecosystems.


Introduction
GPP, the total photosynthetic uptake of carbon by vegetation, plays a key role in understanding the carbon balance between the atmosphere and biosphere [1,2].The magnitude of GPP also defines the ability to produce critical ecosystem services, notably fuel, food, fiber and materials for the construction industry [3].Therefore, quantifying of the spatial-temporal characteristics of GPP has received great attention in the field of global change studies [4][5][6].
GPP values can be estimated from satellite remote sensing products [7][8][9].The first consistent, near-real-time GPP dataset allowing the estimation of global vegetation in an eight-day interval at 1-km resolution is the MOD17A2 MODIS product [10,11].The MODIS GPP data correlated well with the GPP records from twelve flux towers across North America, indicating that MOD17A2 is a reliable data source [12].Even so, there are still many potential error sources in MODIS GPP products arising from input data, the parameters describing the biophysical properties of vegetation and the algorithm itself [13,14].Firstly, the spatial resolution of meteorological reanalysis data is quite different from that of MODIS products, which introduces inaccuracies in estimating atmospheric conditions at a scale that matches heterogeneities in the land surface [15].Secondly, estimation of the FPAR cannot be accurately achieved because the MODIS Leaf Area Index (LAI) is poorly correlated with ground measurements of the vegetation canopy [16].Thirdly, the land cover classification provided by MODIS is not always accurate, and the misclassifications are quite common for similar land cover classes [17].Finally, the use of a constant value attributed to the maximum value for light use efficiency (LUE) cannot be justified for the same biome type [18].All sources of error mentioned above (meteorological data, FPAR product, land cover classification result and LUE value) require separate examination on an individual basis [19].
The analysis of errors arising from the MODIS GPP product is difficult because of inherent limitations in directly measuring GPP values.For the present, the eddy covariance (EC) technique is the frequently-used approach [20][21][22].By performing comparisons between MODIS GPP products and EC flux-derived GPP measurements, much validation work has been conducted to evaluate the accuracies of previous MOD17A2 products [23,24].GPP values from MOD17A2 were comprehensively compared with GPP measurements from flux towers at fifteen research sites across North America by Heinsch et al. (2006) [25].Their results showed that MOD17A2 GPP on average overestimated the site-based GPP by 20-30%.Moreover, the effectiveness of MOD17A2 GPP in Africa was also evaluated by using site meteorological measurements and flux-derived GPP at twelve flux tower sites.The results suggested that GPP values from MOD17A2 agreed with those from the flux towers in wet sites, but were much smaller in dry sites [26].
For MODIS GPP products, many validation efforts have promoted updates and continuous improvements to MOD17A2 [27].Note that the latest MOD17A2H product has essential advantages over the previous version (MOD17A2) because the resolution increased for the first time from 1000 m to 500 m.The major feature of MOD17A2H is its innovation in two sources of input data: meteorological data and FPAR products [28].Firstly, GMAO (Global Modeling and Assimilation Office) reanalysis data substitutes DAO (Data Assimilation Office) reanalysis data to estimate the actual values of meteorological elements.Secondly, the latest MODIS FPAR product, MOD15A2H, substitutes the previous version, MOD15A2, to estimate the actual values of FPAR and to enable improvement in spatial resolution from 1000 m to 500 m.
It was uncertain whether the improvements of MOD17A2H over MO17A2 could estimate GPP more accurately, which will attract great attention from researchers around the world.In this context, the evaluation of MOD17A2H GPP becomes an essential work.Currently, FLUXNET has released the latest global flux dataset, FLUXNET2015 [29].It is clear that the latest EC flux dataset can provide a more accurate reference standard for the validation of MOD17A2H GPP.Therefore, it is of great significance to verify the latest MOD17A2H product using the latest FLUXNET2015 dataset.
In this study, we evaluated MOD17A2H GPP using the FLUXNET2015 dataset at eighteen sites in six biome types around the world (Table 1).Not only the eight-day and annual MODIS GPP values were systematically verified against the flux-derived GPP measurements, but also the role of each input variable was carefully analyzed.The objectives of this study were to assess the performance of MOD17A2H on a global scale by comparisons with EC flux-derived GPP measurements and to explore sources of error among all input parameters (meteorological data, FPAR product, land cover classification result and LUE value).

MOD17A2H GPP Algorithm
Based on the rationale originally described by Monteith (1972) [30], the MOD17A2H GPP algorithm can be expressed as: where FPAR is the fraction of PAR (MJ•m −2 •d −1 ), and PAR can be calculated from the total shortwave solar radiation (SWRad, MJ•m −2 •d −1 ): The magnitude of LUE ε (gC 1) is calculated as: where is the maximum LUE specified by the Biome Properties Look-Up Table (BPLUT) based on biome type, and T s and VPD s are the attenuation scalars influenced by daily minimum temperature (T min , • C) and daytime vapor pressure deficit (VPD avg , Pa), respectively [31][32][33].
The MODIS images from 2001-2014 were downloaded from the NASA Land Processes Distributed Active Archive Center [34].Based on the geographical coordinates of the flux tower sites, 5 km × 5 km cut-outs centered over every site location were extracted to represent: (1) land cover classification (MCD12Q1, Collection 5.1, annual product with 500-m spatial resolution); (2) FPAR (MOD15A2H, Collection 6, eight-day product with 500-m spatial resolution); and (3) GPP (MOD17A2H, Collection 6, eight-day product with 500-m spatial resolution).

GMAO Meteorological Reanalysis Data
The MODIS GPP Collection 6 (MOD17A2H, C6) data used in this paper implement Global Modeling and Assimilation Office (GMAO) Reanalysis data for direct meteorological inputs.The current version of GMAO is an hourly time step dataset with 0.5 latitude degree by 0.625 longitude degree resolution.For every flux tower site, hourly meteorological data including air temperature (T air , • C), specific humidity (SH, kg•kg −1 ), surface pressure (Pres, Pa) and SWRad (W•m −2 ) during 2001-2014 were downloaded from NASA Goddard Space Flight Center [35].Through data preprocessing, key parameters including the daily average air temperature (T avg , • C), daily minimum air temperature (T min , • C), daytime average VPD (VPD avg , Pa) and daily total SWRad (MJ•m −2 •d −1 ) were obtained at a 0.5 • × 0.625 • spatial resolution.Using a spatial nonlinear interpolation scheme developed by Zhao et al. (2005) [12], we further interpolated the GMAO data down to the 500-m MOD17A2H pixel spatial resolution.In order to make comparisons between reanalysis meteorological data and observed meteorological data, the same spatial nonlinear interpolation method was also employed to interpolate the GMAO data to EC flux tower locations [36].

GLASS LAI Data
FPAR is a critical component of the MOD17A2H GPP algorithm.In this study, we introduced the GLASS LAI product to assess the influences of FPAR on the effectiveness of MO17A2H GPP.In previous studies [37][38][39], it had been proven that the GLASS LAI product could provide better estimates of LAI than existing MODIS and CYCLOPES LAI datasets.For each EC flux tower site, the GLASS LAI subsets (eight-day interval at a spatial resolution of 0.05 • ×0.05 • ) during 2001-2014 were downloaded from the Generation and Applications of Global Products of Essential Land Variables website [40].We interpolated the GLASS LAI product down to the 500-m MOD15A2H pixel spatial resolution using a spatial nonlinear interpolation method (the same as the GMAO data).Based on the Beer's law approach, we converted GLASS LAI (m 2 •m −2 ) to GLASS FPAR by following equation: where K is the light extinction coefficient of canopy, which is set to 0.5 [41].
From the FLUXNET-Fluxdata website, we downloaded the newest FLUXNET2015 dataset to obtain site meteorological data and GPP measurements during 2001-2014 for each flux tower site.The daily meteorological data (T avg , T min , PAR and VPD avg ) were calculated from the half-hourly site meteorological measurements, and the daily GPP values were calculated from the half-hourly flux-derived GPP.We then obtained the eight-day site-based GPP by accumulating the daily calculations to ensure that our data remained consistent with the 8-day interval of the MOD17A2H GPP products.

Experimental Design
To distinguish the impacts of meteorological data and FPAR on MOD17A2H GPP, we conducted the following four groups of experiments to recalculate GPP under dissimilar circumstances: (1) GMAO meteorology data and MOD15A2H FPAR were chosen as meteorology inputs and FPAR inputs, respectively (MOD_GMAO GPP) (the same as the standard MOD17A2H algorithm); (2) using site meteorological data as meteorology inputs and MOD15A2H FPAR as FPAR inputs (MOD_Tower GPP); (3) GMAO meteorology data were used as meteorology inputs, and GLASS FPAR were introduced as FPAR inputs (GMAO_GLASS GPP); (4) using site meteorological data and GLASS FPAR to estimate GPP (Tower_GLASS GPP).At last, the four types of GPP estimates were compared with GPP values measured by EC flux tower sites (Flux_Tower GPP), respectively.

Analytical Methods
Three statistical indices, including the coefficient of determination (R 2 ), bias and root mean square error (RMSE, %), were calculated to assess accuracies and uncertainties.The value of R 2 denoted the amount of variation in the measurements that could be predicted by the estimates, and the value of bias expresses the differences between the estimates and the measurements.The value of RMSE was computed by: where M was the mean value of the measurements and E i and M i are the estimated and measured values, respectively.

Eight-Day GPP
The eight-day Flux_Tower GPP were compared with the results of MOD_GMAO GPP, MOD_Tower GPP, GMAO_GLASS GPP and Tower_GLASS GPP (Figure 1, Table 2).In general, the MOD_Tower GPP algorithm (R 2 = 0.52) did not result in obvious improvements in fitting the Flux_Tower GPP when comparing with the MOD_GMAO GPP algorithm (R 2 = 0.52).However, the GMAO_GLASS GPP algorithm (R 2 = 0.63) and the Tower_GLASS GPP algorithm (R 2 = 0.65) made considerable improvements, and the Tower_GLASS GPP algorithm exhibited the best estimating effectiveness.It should be noted that all of the above four algorithms underestimated the Flux_Tower GPP to varying degrees.meteorological data and GLASS FPAR to estimate GPP (Tower_GLASS GPP).At last, the four types of GPP estimates were compared with GPP values measured by EC flux tower sites (Flux_Tower GPP), respectively.

Analytical Methods
Three statistical indices, including the coefficient of determination (R 2 ), bias and root mean square error (RMSE, %), were calculated to assess accuracies and uncertainties.The value of R 2 denoted the amount of variation in the measurements that could be predicted by the estimates, and the value of bias expresses the differences between the estimates and the measurements.The value of RMSE was computed by: where M was the mean value of the measurements and Ei and Mi are the estimated and measured values, respectively.

Eight-Day GPP
The eight-day Flux_Tower GPP were compared with the results of MOD_GMAO GPP, MOD_Tower GPP, GMAO_GLASS GPP and Tower_GLASS GPP (Figure 1, Table 2).In general, the MOD_Tower GPP algorithm (R 2 = 0.52) did not result in obvious improvements in fitting the Flux_Tower GPP when comparing with the MOD_GMAO GPP algorithm (R 2 = 0.52).However, the GMAO_GLASS GPP algorithm (R 2 = 0.63) and the Tower_GLASS GPP algorithm (R 2 = 0.65) made considerable improvements, and the Tower_GLASS GPP algorithm exhibited the best estimating effectiveness.It should be noted that all of the above four algorithms underestimated the Flux_Tower GPP to varying degrees.Table 2. Biome-based comparison of eight-day GPP derived from the site estimates with those derived from the MODIS algorithms (MOD_GMAO GPP, MOD_Tower GPP, GMAO_GLASS GPP and Tower_GLASS GPP).ENF: evergreen needleleaf forests; EBF: evergreen broadleaf forests; DBF: deciduous broadleaf forests; MF: mixed forests; GRA: grasslands; CRO: croplands.From the perspective of individual sites, the four algorithms (MOD_GMAO GPP, MOD_Tower GPP, GMAO_GLASS GPP and Tower_GLASS GPP) underestimated the Flux_Tower GPP more obviously (Figure 2).The Flux_Tower GPPs were underestimated by above four algorithms at 15 of the 18 sites.For the remaining three sites (AU-ASM, AU-Whr and US-Wkg), GMAO_GLASS GPP and the Tower_GLASS GPP algorithms still underestimated the site-based GPP.In detail, AR-Vir, BR-Sa3 and PA-SPs suffered the highest degree of underestimation by the four algorithms (Figure 2a,g,p, respectively), and CA-Gro, CZ-BK2 and FR-Pue suffered the lowest degree of underestimation (Figure 2h,I,m, respectively).With respect to the estimating effectiveness of the four algorithms, the MOD_GMAO GPP, MOD_Tower GPP, GMAO_GLASS GPP and Tower_GLASS GPP algorithms fitted relatively effectively for one site (US-Wkg), two sites (AU-Rig and AU-Whr), ten sites (AU-ASM, BE-Lon, BE-Vie, CA-Gro, CZ-BK2, DK-Sor, FR-Gri, FR-Pue, IT-Ren and PA-SPs) and four sites (AR-Vir, DE-Hai, GH-Ank and US-CRT), respectively.None of the four algorithms successfully fitted the Flux_Tower GPP for BR-Sa3 (Figure 2g).The GMAO_GLASS GPP algorithm had a relative advantage in terms of bias (−14 gC•m −2 •8-d −1 ).
The eight-day GPP values of MOD_GMAO GPP, MOD_Tower GPP, GMAO_GLASS GPP and Tower_GLASS GPP were also compared with the Flux_Tower GPP at different biome types (Table 2).As the analyses of individual sites showed, all above algorithms underestimated the Flux_Tower GPP for all six biome types.In detail, deciduous broadleaf forests (DBF) suffered the highest degree of underestimation, mainly because the four algorithms seriously underestimated the Flux_Tower GPP of DE-Hai (Figure 2j) and DK-Sor (Figure 2k).In contrast, the degree of underestimation was low for grasslands (GRA).This could be largely attributed to the fact that the four algorithms effectively estimated the Flux_Tower GPP of CZ-BK2 (Figure 2i) and US-Wkg (Figure 2r).With respect to the estimating effectiveness of the four algorithms, the GMAO_GLASS GPP and Tower_GLASS GPP algorithms were superior to the MOD_GMAO GPP and MOD_Tower GPP algorithms across all six biome types.The GMAO_GLASS GPP algorithm provided the smallest values of bias for all biomes, and the Tower_GLASS GPP algorithm exhibited the highest values of R 2 for most biomes, except croplands (CRO); for croplands, no algorithm fitted the Flux_Tower GPP more effectively than the GMAO_GLASS GPP algorithm because it provided the highest value of R 2 (0.

Meteorology
In this study, the MOD_GMAO GPP algorithm followed the standard MOD17A2H GPP algorithm, which used the GMAO meteorology reanalysis data.The GMAO meteorology reanalysis data were replaced with the site meteorological data for the MOD_Tower GPP algorithm.The MOD_GMAO GPP and MOD_Tower GPP algorithms were thus compared to determine the influences of the meteorological data on the GPP estimates.The results showed that the MOD_Tower GPP algorithm did not estimate eight-day GPP more reliably than the MOD_GMAO GPP algorithm (Figure 1a,b; Table 2).Conversely, estimates of annual GPP from the MOD_Tower GPP algorithm were inferior to thoseof the MOD_GMAO algorithm (Figure 3a,b; Table 3).
We directly compared the daily GMAO meteorology reanalysis data and the site meteorological data (Figure 4).The two datasets were highly significantly correlated for T avg (R 2 = 0.97) and T min (R 2 = 0.94) (Figure 4a,b), which implied that the GMAO meteorology reanalysis data could be reliably used to estimate the on-site temperature conditions.Additionally, VPD avg was also statistically significantly correlated (R 2 = 0.88) (Figure 4c), which indicated that the on-site moisture variations could also be dependably estimated from the GMAO meteorology reanalysis datasets.With respect to PAR, the GMAO meteorology reanalysis data generally overestimated the site meteorological data (Figure 4d), and the correlation (R 2 = 0.85) was stronger than that of VPD avg , but weaker than those of T avg and T min .

FPAR
Unlike the MOD_GMAO GPP algorithm, the GMAO_GLASS GPP algorithm followed the MOD17A2H GPP algorithm by substituting GLASS FPAR (derived from GLASS LAI datasets) for MODIS FPAR (derived from MOD15A2H products).Therefore, comparing these two types of algorithms will be helpful for understanding the impacts of FPAR on the GPP estimates.In the case of the eight-day GPP estimates, the effectiveness of the GMAO_GLASS GPP algorithm for estimating the Flux_Tower GPP was much better than that of the MOD_GMAO GPP algorithm (Figure 1a,c; Table 2).Additionally, the GMAO_GLASS GPP algorithm achieved substantial improvements over the MOD_GMAO GPP algorithm for the annual GPP estimates (Figure 3a,c; Table 3).The GMAO_GLASS GPP algorithm estimated the Flux_Tower GPP far more effectively than the MOD_GMAO GPP algorithm for EBF site (Tables 2 and 3).
In our study, comparisons were made between GLASS FPAR and MODIS FPAR (Figure 5, Table 4).Overall, the correlation between the two types of FPAR was not strong (R 2 = 0.45).The correlations between GLASS FPAR and MODIS FPAR were also weak for every single biome type, except ENF; the correlation between them for EBF was lowest (R 2 = 0.02).Obviously, GLASS FPAR overestimated MODIS FPAR at most sites.This overestimation problem was most serious at GH-Ank (bias = 0.479), PA-SPs (bias = 0.375) and BR-Sa3 (bias = 0.292).Moreover, GLASS FPAR also overestimated MODIS FPAR for all six biome types.

FPAR
Unlike the MOD_GMAO GPP algorithm, the GMAO_GLASS GPP algorithm followed the MOD17A2H GPP algorithm by substituting GLASS FPAR (derived from GLASS LAI datasets) for MODIS FPAR (derived from MOD15A2H products).Therefore, comparing these two types of algorithms will be helpful for understanding the impacts of FPAR on the GPP estimates.In the case of the eight-day GPP estimates, the effectiveness of the GMAO_GLASS GPP algorithm for estimating the Flux_Tower GPP was much better than that of the MOD_GMAO GPP algorithm (Figure 1a,c; Table 2).Additionally, the GMAO_GLASS GPP algorithm achieved substantial improvements over the MOD_GMAO GPP algorithm for the annual GPP estimates (Figure 3a,c; Table 3).The GMAO_GLASS GPP algorithm estimated the Flux_Tower GPP far more effectively than the MOD_GMAO GPP algorithm for EBF site (Tables 2 and 3).
In our study, comparisons were made between GLASS FPAR and MODIS FPAR (Figure 5, Table 4).Overall, the correlation between the two types of FPAR was not strong (R 2 = 0.45).The correlations between GLASS FPAR and MODIS FPAR were also weak for every single biome type, except ENF; the correlation between them for EBF was lowest (R 2 = 0.02).Obviously, GLASS FPAR overestimated MODIS FPAR at most sites.This overestimation problem was most serious at GH-Ank (bias = 0.479), PA-SPs (bias = 0.375) and BR-Sa3 (bias = 0.292).Moreover, GLASS FPAR also overestimated MODIS FPAR for all six biome types.Significance levels: ** p < 0.01.

Land Cover Classification
We obtained the pixel values within a range of 5 km × 5 km centered over each site, as well as the numbers of corresponding pixel values from the MCD12Q1 image data during 2001-2013.MODIS land cover classification results at each site were then calculated from the statistics of the

Land Cover Classification
We obtained the pixel values within a range of 5 km × 5 km centered over each site, as well as the numbers of corresponding pixel values from the MCD12Q1 image data during 2001-2013.MODIS land cover classification results at each site were then calculated from the statistics of the pixel values.
Additionally, the actual land cover type of each site was obtained by consulting site descriptions provided by the FLUXNET2015 dataset.
Remote Sens. 2017, 9, 418 12 of 20 pixel values.Additionally, the actual land cover type of each site was obtained by consulting site descriptions provided by the FLUXNET2015 dataset.We also compared MODIS land cover classification results and actual land cover types in this study (Figure 6).It was found that MCD12Q1 misclassified eleven of the eighteen sites.As shown in Figure 6, the misclassified sites included AR-Vir (mistook ENF for savannas), AU-ASM (mistook ENF for open shrublands), AU-Rig (mistook GRA for CRO), AU-Whr (mistook EBF for woody savannas), CZ-BK2 (mistook GRA for ENF), DE-Hai (mistook DBF for MF), DK-Sor (mistook DBF for CRO), FR-Pue (mistook EBF for MF), IT-Ren (mistook ENF for MF), PA-SPs (mistook GRA for savannas) and US-Wkg (mistook GRA for open shrublands).

Light Use Efficiency
The constant values of ε max used in the MOD17A2H GPP algorithm varied for the different types of land cover; each biome type had an individual ε max value.By using Equations ( 1) and (3) in Section 2.1.1,we recalculated the inferred ε max for each land cover type and then directly compared them to the MOD17A2H ε max (Figure 7, Table 5).The results showed that the inferred ε max was higher than MOD17A2H ε max by 31% for ENF, 21% for EBF, 72% for DBF, 61% for MF, 31% for GRA and 25% for CRO.Moreover, from the perspective of a single site, the inferred ε max was higher than MOD17A2H ε max at almost all sites except AU-Whr (bias = 0.093).

Light Use Efficiency
The constant values of εmax used in the MOD17A2H GPP algorithm varied for the different types of land cover; each biome type had an individual εmax value.By using Equations ( 1) and (3) in Section 2.1.1,we recalculated the inferred εmax for each land cover type and then directly compared them to the MOD17A2H εmax (Figure 7, Table 5).The results showed that the inferred εmax was higher than MOD17A2H εmax by 31% for ENF, 21% for EBF, 72% for DBF, 61% for MF, 31% for GRA and 25% for CRO.Moreover, from the perspective of a single site, the inferred εmax was higher than MOD17A2H εmax at almost all sites except AU-Whr (bias = 0.093).

Discussion
In this study, Flux_Tower GPP from eighteen sites were used to verify four groups of GPP estimates recalculated following the MOD17A2H GPP algorithms (i.e., MOD_GMAO GPP, MOD_Tower GPP, GMAO_GLASS GPP and Tower_GLASS GPP).It was found that the four groups of GPP estimates deviated from the flux-derived GPP values to varying degrees.The extent of the deviations varied depending on the site location, the type of biome and the time interval.We can attribute these deviations to different factors.Undoubtedly, the effectiveness of GPP estimation critically depends on the precision of every input data involved in the MOD17A2H GPP algorithm.Therefore, the meteorological data, FPAR product, land cover classification result and ε max value were all potential sources of error [43].The accuracy of the flux-derived GPP value and the reliability of the MOD17A2H GPP algorithm itself should also be considered [44].

Impact of Meteorological Data on MODIS GPP
There are three meteorological factors (PAR, T min and VPD avg ) involved in the MOD17A2H GPP algorithm, which suggests that meteorology reanalysis data could be the main source of error in the GPP estimates [44].For the previous MOD17A2 GPP product, the meteorological input data were selected from the DAO meteorology dataset (1.00 • × 1.25 • spatial resolution).In contrast, the latest MOD17A2H product uses an updated version of the GMAO meteorology dataset (0.5 • × 0.625 • spatial resolution).However, the spatial resolution of GMAO meteorology data is still greatly inferior to that of the MODIS GPP product (500 m).Although a strict interpolation was completed in the data pre-processing stage, the errors caused by the original low resolution could not be eliminated.Moreover, the interpolation process may have resulted in additional uncertainties.Analytic comparisons between the GMAO meteorology reanalysis dataset and the site-based meteorological observations indicated that GMAO data were effective for estimating on-site temperature elements and moisture conditions (Figure 4a-c).
Compared with the previous research results about MOD17A2 meteorology data, the GMAO data adopted in the current MOD17A2H product achieved substantial progress in estimating actual values of meteorological elements [45,46].This demonstrated the considerable effectiveness of the efforts that MODIS GPP has taken in the version upgrading process to reduce the error derived from meteorological data.Note that as the source of original energy for GPP, PAR plays an extremely important role among the four meteorological elements involved in the MOD17A2H algorithm.It can be seen from this work that the effectiveness of GMAO data in estimating PAR was less than that of GMAO data in estimating the temperature and moisture elements (Figure 4d), which is mainly because the GMAO data largely overestimated the site-measured PAR.
In this study, we replaced the GMAO meteorology reanalysis dataset with the site-based meteorological observations to recalculate the MOD17A2H GPP algorithm and gained new GPP estimates (MOD_Tower GPP).Unexpectedly, the eight-day GPP estimates recalculated using the site-based meteorological observations did not result in obvious improvements (Figure 1a,b).Conversely, the annual GPP estimates recalculated using the site-based meteorological observations were not as accurate as those calculated using the GMAO meteorology reanalysis dataset (Figure 3a,b).The implication of these results is that an improvement in meteorological input data will not necessarily improve the reliability of GPP estimation by the MOD17A2H GPP algorithm, emphasizing the requirement of identifying additional sources of uncertainty.

Impact of FPAR Product on MODIS GPP
To assess the influences of FPAR on the accuracy of GPP estimation, we introduced an improved FPAR dataset (GLASS FPAR) and designed the GMAO_GLASS GPP algorithm.Our results showed that the greatest improvements were gained in estimating the eight-day GPP values when the MODIS FPAR was replaced by the GLASS FPAR (Figure 1a,c).Similarly, the precision in estimating annual GPP was also improved substantially after substituting GLASS FPAR for MODIS FPAR (Figure 3a,c).By comparing two groups of FPAR data, we found that the values of MODIS FPAR were generally smaller than that of GLASS FPAR (Table 4).Indeed, the MOD17A2H GPP (MOD_GMAO GPP algorithm) generally underestimated the flux derived GPP values (Tables 2 and 3).The improved GLASS FPAR data thus balanced the underestimation resulting from the use of MODIS FPAR to some extent.
Admittedly, GLASS FPAR adopted in the GMAO_GLASS GPP algorithm was improved from the previous MOD15A2 LAI/FPAR product, which had a spatial resolution of 1 km.The MOD_GMAO GPP algorithm adopted the latest MOD15A2H FPAR product, whose spatial resolution increased to 500 m.However, our results showed that the GMAO_GLASS GPP algorithm estimated the site-based GPP much more effectively than the MOD_GMAO GPP algorithm.This indicated that what had been done in MODIS GPP products to reduce the error derived from FPAR during the version upgrading process was still far from sufficient.
The qualities of the meteorological data and the FPAR product were simultaneously improved in the Tower_GLASS GPP algorithm.We simply substituted the site-based meteorological observations for the GMAO meteorology reanalysis dataset and replaced the MODIS FPAR with GLASS FPAR.The results showed that the effectiveness of eight-day GPP estimation had the greatest improvements (Figure 1a,d).However, this replacement strategy did not produce the most effective algorithm in terms of annual GPP estimates (Figure 3a,c,d).This may indicate that the accuracy of GPP estimates could be enhanced by simultaneously improving the qualities of multiple inputs of the MOD17A2H algorithm; but it did not imply that improvements could be gained simply by replacing the input data because the overall gain in algorithm performance was not simply the sum of the individual gains.More detailed analysis needs to be conducted to evaluate the interaction effects of the joint variation of the model inputs in future research [1].
Compared with MOD17A2, MOD17A2H improved in terms of meteorological data and FPAR.Therefore, four groups of experiments (i.e., MOD_GMAO GPP, MOD_Tower GPP, GMAO_GLASS GPP and Tower_GLASS GPP) were designed to verify the contributions of the upgraded meteorology dataset and FPAR product to the estimation accuracy of GPP.It was discovered that the GMAO_GLASS GPP and Tower_GLASS GPP algorithms achieved enormous performance gains over the MOD_GMAO GPP algorithm in terms of flux-derived GPP estimates and that the performance of the GMAO_GLASS GPP algorithm was almost the same as that of the Tower_GLASS GPP algorithm.This may indicate that the error from FPAR was larger than that from the meteorological data in the current version of the MODIS GPP product.

Impact of Land Cover Classification Result on MODIS GPP
The MOD17A2H GPP algorithm uses the MCD12Q1 product as its land cover classification scheme.The accuracy of the land cover classification directly determines the veracity of ε max and further influences the effectiveness of the GPP estimates.Therefore, the influences of land cover misclassification on GPP estimation are non-negligible.
In our study, we found that MCD12Q1 misclassified most of the sites (Figure 6).These misclassifications directly led to the misestimating FPAR by MOD15A2H and the misjudgments of ε max by BPLUT and thus affected the MOD17A2H GPP estimates for corresponding sites.

Impact of LUE Value on MODIS GPP
The value of ε max is biomespecific, representing the maximum light use efficiency (LUE) of the corresponding vegetation in the process of photosynthesis.For a given biome type, the value of ε max is constant and assigned by the BPLUT.Although the BPLUT has corrected and updated ε max values, the colossal diversity of terrestrial environments and climatic conditions does not allow a value of ε max at one location to be used reliably in another location, even for the same biome type [47].The optimized ε max that we deduced for the eighteen sites in this study greatly differed from the MOD17A2H ε max that was specified by the BPLUT (Figure 7, Table 5).Because the ε max value used in the MOD17A2H GPP algorithm sets the foundation for the actual conversion of light energy, the misestimation of its value will inherently reduce the accuracy of GPP estimations.

Uncertainties, Errors and Accuracies
In this study, our evaluation of MOD17A2H GPP assumed that the values of GPP derived from EC flux tower sites are reliable, even though it is well known that the flux-derived GPP values encompass many uncertainties [48].The ecosystem respiration (R eco ) is an essential variable when the flux-derived GPP value is calculated.However, accurate estimation of R eco is difficult, and misestimation can produce random and systematic errors in estimating GPP.Moreover, mismatches in scale between MODIS pixels and tower flux footprints lead to further uncertainties.Additionally, researchers have discovered structural errors in the MODIS GPP algorithm [49][50][51].For example, it was found that the MODIS GPP algorithm cannot correctly identify the contribution of leaves in both insolated and shaded locations to the calculation of total GPP when compared with GPP estimates from a two-leaf process-based model [52].

Conclusions
In our study, we systematically assessed both eight-day and annual MOD17A2H GPP using the latest FLUXNET2015 dataset at eighteen sites in six different ecosystems across the globe during 2001-2014.The sensitivities of the MOD17A2H algorithm to the meteorology reanalysis dataset and FPAR product was evaluated by introducing site-based meteorological observations and an upgraded GLASS FPAR product.We further explored the potential error contributions from land cover classification and ε max .Our principal conclusions are summarized as follows: (1) The effectiveness of the standard MOD17A2H product (i.e., MOD_GMAO GPP) for estimating annual GPP was poor (R 2 = 0.62) and even worse over eight days (R 2 = 0.52).Replacing the GMAO meteorology reanalysis dataset with the site-based meteorological observations (i.e., MOD_Tower GPP) failed to improve the correlations with the flux-derived eight-day GPP (R 2 = 0.52) and annual GPP (R 2 = 0.56).However, great improvements in estimating the flux-derived annual GPP (R 2 = 0.79) were gained by replacing MODIS FPAR with GLASS FPAR (i.e., GMAO_GLASS GPP).This may indicate that in the current version of the MODIS GPP product, the error from FPAR was larger than that from the meteorological data.When the meteorological dataset and FPAR product were replaced by upgraded data simultaneously (i.e., Tower_GLASS GPP), the accuracy in estimating the flux-derived eight-day GPP (R 2 = 0.65) was significantly improved.(2) There are four possible sources of error related to the input data of MOD17A2H algorithm: meteorological reanalysis dataset, FPAR product, land cover classification result and the ε max value.The GMAO meteorology reanalysis dataset commonly overestimated the site-measured PAR, and MODIS FPAR underestimated the upgraded FPAR in most flux tower sites.In addition, MCD12Q1 exhibited frequent misclassification errors, and the MOD17A2H ε max values were much smaller than the inferred ε max values for all six ecosystems discussed in our study.
From the above analyses, the latest MODIS GPP product, MOD17A2H, still has some defects in estimating the actual value of GPP.We suggest that FPAR and land cover classification products need to be improved, and ε max should be adjusted in the next versions of datasets for GPP estimation.
35) and the lowest values of both bias (−5 gC•m −2 •8-d −1 ) and RMSE (89%).Remote Sens. 2017, 9, 418 7 of 20 MOD_Tower GPP algorithms across all six biome types.The GMAO_GLASS GPP algorithm provided the smallest values of bias for all biomes, and the Tower_GLASS GPP algorithm exhibited the highest values of R 2 for most biomes, except croplands (CRO); for croplands, no algorithm fitted the Flux_Tower GPP more effectively than the GMAO_GLASS GPP algorithm because it provided the highest value of R 2 (0.35) and the lowest values of both bias (−5 gC•m −2 •8-d −1 ) and RMSE (89%).

Figure 2 .
Figure 2.(a-r) Time series of eight-day GPP derived from the tower estimates (i.e., Flux_Tower GPP) and four different experimental groups (i.e., MOD_GMAO GPP, MOD_Tower GPP, GMAO_GLASS GPP and Tower_GLASS GPP) at all eighteen sites.The full name for each site is listed in Table1.The units of the mean and bias are gC•m −2 •8-d −1 , and the RMSE is expressed as a percentage.Significance levels: ** p < 0.01.

Figure 2 .
Figure2.(a-r) Time series of eight-day GPP derived from the tower estimates (i.e., Flux_Tower GPP) and four different experimental groups (i.e., MOD_GMAO GPP, MOD_Tower GPP, GMAO_GLASS GPP and Tower_GLASS GPP) at all eighteen sites.The full name for each site is listed in Table1.The units of the mean and bias are gC•m −2 •8-d −1 , and the RMSE is expressed as a percentage.Significance levels: ** p < 0.01.

Figure 4 .
Figure 4. Scatter plots of the daily GMAO meteorology data against the daily site meteorology data for: (a) T avg ; (b) T min ; (c) VPD avg ; and (d) PAR.Significance levels: ** p < 0.01.

Figure 6 .
Figure 6.University of Maryland (UMD) land cover classification (MCD12Q1 Collection 5.1, Land Cover Classification Type 2) for the eighteen sites.The red columns indicate the actual land cover at each site.

Figure 6 .
Figure 6.University of Maryland (UMD) land cover classification (MCD12Q1 Collection 5.1, Land Cover Classification Type 2) for the eighteen sites.The red columns indicate the actual land cover at each site.

Table 4 .
Comparison of the FPAR derived from the GLASS LAI datasets (GLASS FPAR) with those from the MOD15A2H products (MODIS FPAR) for various biomes.

Table 4 .
Comparison of the FPAR derived from the GLASS LAI datasets (GLASS FPAR) with those from the MOD15A2H products (MODIS FPAR) for various biomes.