Evaluation of MODIS Gross Primary Production across Multiple Biomes in China Using Eddy Covariance Flux Data

MOD17A2 provides near real-time estimates of gross primary production (GPP) globally. In this study, MOD17A2 GPP was evaluated using eddy covariance (EC) flux measurements at eight sites in five various biome types across China. The sensitivity of MOD17A2 to meteorological data and leaf area index/fractional photosynthetically active radiation (LAI/FPAR) products were examined by introducing site meteorological measurements and improved Global Land Surface Satellite (GLASS) LAI products. We also assessed the potential error contributions from land cover and maximum light use efficiency (εmax). The results showed that MOD17A2 agreed well with flux measurements of annual GPP (R2 = 0.76) when all biome types were considered as a whole. However, MOD17A2 was ineffective for estimating annual GPP at mixed forests, evergreen needleleaf forests and croplands, respectively. Moreover, MOD17A2 underestimated flux derived GPP during the summer (R2 = 0.46). It was found that the meteorological data used in MOD17A2 failed to properly estimate the site measured vapor pressure deficits (VPD) (R2 = 0.31). Replacing the existing LAI/FPAR data with GLASS LAI products reduced MOD17A2 GPP uncertainties. Though land cover presented the fewest errors, εmax prescribed in MOD17A2 were much lower than inferred εmax calculated from flux data. Thus, the qualities of meteorological data and LAI/FPAR products need to be improved, and εmax should be adjusted to provide better GPP estimates using MOD17A2 for Chinese ecosystems.


Introduction
Gross primary production (GPP), which is defined as the overall photosynthetic fixation of carbon by plants, is an important variable in studies of the carbon balance between the atmosphere and biosphere [1,2].GPP is also the basis for essential ecosystem services such as food, fiber, fuel and construction materials [3].Thus, the quantification of GPP has become a topic of wide concern in global change studies [4][5][6][7].
By using optical and near-infrared spectral wavelengths, GPP can be estimated from satellite remote sensing [8,9].In particular, the Moderate Resolution Imaging Spectroradiometer (MODIS) primary production products (MOD17A2) are the first regular, near-real-time GPP datasets for the repeated monitoring of global vegetation at a 1-km resolution every eight days [10].Zhao et al. (2005) [11] had demonstrated that MODIS GPP fitted well with GPP derived from 12 flux towers over North America, indicating MOD17A2 were reliable products.However, for MODIS GPP products, there are still exist many potential sources of errors that arise from the input datasets, the parameters used to describe the biophysical behavior of vegetation, and the algorithm itself [12].First, there is a large disparity between the spatial resolution of meteorological data and the resolution of MODIS products, which will provide inaccurate atmospheric conditions at scales consistent with land surface heterogeneities [13].Second, the MODIS leaf area index (LAI) has a poor correlation with the ground measurements, which will lead to an erroneous estimation for fraction of photosynthetically active radiation (FPAR) [11].Third, the accuracy of MODIS land cover classification is not very satisfactory, and most mistakes are between similar classes [14].Finally, it does not conform to reality to assign a constant value of maximum light use efficiency (LUE) to the same biome type [5].As a result, each of these error sources (i.e., meteorology, LAI/FPAR, land cover and LUE) requires a corresponding validation procedure and must be examined separately [15].
The error analysis of MODIS GPP products is a challenging task because of the difficulty in making direct measurements of GPP values.One widely used approach uses the eddy covariance (EC) technique, which measures the fluxes of carbon, water and energy between the atmosphere and land [16,17].A number of validation efforts have been established to evaluate the accuracy of MOD17A2 products that use time series comparisons between MODIS-based and EC flux tower-based GPP for one or more 1-km 2 cells centered on towers [18,19].Turner et al. (2003) [20] evaluated the MODIS GPP products at two sites: a temperate forest site and a boreal forest site.Their results showed that, relative to the flux tower measurements, MODIS overestimated the GPP by 35% at the boreal forest site, but the MODIS estimates were comparable to the tower results for the temperate forest site.Heinsch et al. (2006) [14] carried out a comprehensive evaluation of MOD17A2 by comparing MODIS-based GPP to flux tower-based GPP at 15 research sites in six different biome classes across North America.The authors reported that the MODIS GPP products overestimated the tower-based calculations by 20%-30% on average.They further found that the use of MODIS GPP products with DAO (NASA'S Data Assimilation Office) meteorology overestimated the annual GPP, whereas the use of tower-specific meteorology in the MODIS GPP calculations led to underestimates.The performance of MODIS GPP in Africa has also been evaluated using in situ measurements of meteorology and flux tower GPP for 12 sites.The study indicated that MOD17A2 agreed well with the tower-based GPP for wet sites, whereas the estimates were too small for dry sites [21].
Based on EC flux data provided by the Chinese FLUX Observation and Research Network (ChinaFLUX), several studies have assessed the performances of MODIS GPP in Chinese different biome types [22,23].For example, Zhang et al. (2008) [24] evaluated MODIS GPP by using estimated GPP from EC flux measurements over an alpine meadow on the Tibetan Plateau.Their results showed that the mean annual MODIS GPP accounted for 1/2-2/3 of the flux-based GPP in the study region.Liu et al. (2015) [25] reported that MODIS GPP performed poorly for evergreen forests but provided accurate estimates for grassland and mixed forests.However, these studies explored error sources mainly focused on one or two aspects, while comprehensive and detailed error analyses are still needed.Thus, we not only methodically evaluated the eight-day, seasonal and annual MODIS GPP against tower GPP, but also carefully analyzed each input.In this paper, we evaluated MOD17A2 GPP using EC flux tower data at eight sites in five various biome types across China (Figure 1 and Table 1).The objectives of our study were to evaluate the performance of MOD17A2 in China through comparisons with GPP measured at EC flux towers and to examine the potential error contributions of all input variables (meteorology, LAI/FPAR, land cover and LUE) used in the algorithm.

MODIS GPP Algorithm
The GPP calculation used in the MODIS GPP algorithm (MOD17A2) is based on the original logic of Monteith (1972) [32] that relates gross photosynthesis to the amount of photosynthetically active radiation (PAR) absorbed by plants over a growing season.As described in detail by Running et al. (2000) [33], the algorithm was developed as follows: where FPAR is the fraction of PAR absorbed by the vegetation canopy, and PAR is determined as a fraction of the downward solar shortwave radiation (SWRad): PAR " SWRad ˆ0.45 The magnitude of LUE ε in Equation ( 1) is calculated as where ε max is the maximum LUE obtained from a Biome Properties Look-Up Table (BPLUT) on the basis of biome type, and T s and VPD s are the attenuation scalars from cold temperature (low daily minimum temperature, T min ) and water stresses (high daily vapor pressure deficit, VPD), respectively [34,35].Following the MODIS GPP algorithm, essential MODIS product images from 2003 to 2005 were downloaded from the Level 1 and Atmosphere Archive and Distribution System (LAADS) website [36].Using the latitudes and longitudes of the flux tower sites, we obtained 5 km ˆ5 km cutouts centered over each tower location that represented the: (1) land cover classification (MCD12Q1, C5.1, 500-m resolution, annual product); (2) LAI and FPAR (MOD15A2, C5, 1-km resolution, eight-day product); and (3) GPP (MOD17A2, C5.5, 1-km resolution, eight-day product).
In particular, MCD12Q1 is used in the MOD15A2 and MOD17A2 calculations.Unlike the MOD15A2 algorithm, which uses MCD12Q1 Land Cover Classification Type 3 (LAI/FPAR scheme) to calculate FPAR, the MOD17A2 algorithm uses MCD12Q1 Land Cover Classification Type 2 (University of Maryland land cover classification scheme) to map biome-specific physiological parameters (ε max , maximum and minimum air temperature and VPD) according to the BPLUT [37].

NCEP-DOE Reanalysis II Data
The MODIS GPP Collection 5.5 (MOD17A2, C5.5) data used in this paper implement National Center for Environmental Prediction-Department of Energy (NCEP-DOE) Reanalysis II data for direct meteorological inputs.For each flux tower site, six-hour data composed of air temperature (T air , ˝C), downward shortwave radiation (SWRad, W¨m ´2), surface pressure (Pres, Pa), and specific humidity (SH, kg¨kg ´1) for 2003-2005 were downloaded from the Earth System Research Laboratory (ESRL) [38].From these data, we obtained the daily average air temperature (T avg , ˝C), daily minimum air temperature (T min , ˝C), daily average VPD (VPD avg , Pa), and daily total SWRad (MJ¨m ´2¨d ´1) at a 1.9 ˝ˆ1.9 ˝resolution.To interpolate the 1.9 ˝ˆ1.9 ˝resolution NCEP-DOE data down to the 1-km MODIS pixel resolution, the daily data for each pixel were interpolated using a spatial nonlinear interpolation scheme following Zhao et al. (2005) [11].

‚
Site Meteorological Data Site meteorological data (PAR, air temperature and VPD) during 2003-2005 were obtained directly from the flux towers or were derived from meteorological station data provided by the China Meteorological Administration (CMA) [39].Based on the half-hourly flux tower measurements, which are composed of air temperature (T air , ˝C) and either PAR (µmol¨m ´2¨s ´1) or downward shortwave radiation (SWRad, W¨m ´2), we obtained the daily T avg , T min and PAR.In addition, the daily average e a (i.e., the actual air vapor pressure, Pa) and daily average RH (i.e., relative humidity, %) were selected from the meteorological station datasets to calculate VPD avg [40].The site meteorological data (T avg , T min , PAR and VPD avg ) were then directly compared to the NCEP-DOE Reanalysis II data to evaluate how well the reanalysis datasets represented the local climatic conditions.

GLASS LAI Data
For the MODIS GPP algorithm, the FPAR is an essential input.To examine the potential error contributions from LAI/FPAR, the Global Land Surface Satellite Leaf Area Index (GLASS LAI) products, which are an improved LAI dataset, were introduced into our study [41].The GLASS LAI algorithm uses general regression neural networks (GRNNs) to retrieve LAI data.In addition, the GRNNs are trained by the fused LAI values from time-series multi-sensor remote sensing data, including MODIS and CYCLOPES LAI products and the corresponding reprocessed MODIS reflectance values [42].It has been demonstrated that the GLASS LAI algorithm is able to produce temporally continuous LAI datasets with significantly improved accuracies compared with current MODIS and CYCLOPES LAI products [43].In this paper, the GLASS LAI subsets (0.05 ˝ˆ0.05 ˝resolution, eight-day interval) during the 2003 to 2005 period centered over each of the flux tower locations were obtained from the Generation and Applications of Global Products of Essential Land Variables website [44].
Using a simple Beer's Law approach [45], the GLASS LAI data can be converted to FPAR data by: where K is the canopy light extinction coefficient, which is set to 0.5.To explore the influence of FPAR on GPP, we replaced the MOD15A2 FPAR with the above FPAR, which originated from GLASS LAI, to recalculate MODIS GPP in this study.

EC Flux Data
Eight EC flux tower sites across China were included in this study (Figure 1).The sites cover a diversity of biome types and climate regimes, including mixed forests, evergreen needleleaf forests, evergreen broadleaf forests, croplands and grasslands [29].A brief description of the sites is presented in Table 1, and detailed information is available at the ChinaFLUX website [46].
The EC flux data measured at the eight flux tower sites for 2003-2005 were downloaded from the ChinaFLUX website.The daily values, which were composed of the net ecosystem exchange (NEE) and ecosystem respiration (Re), were used to calculate GPP: where the GPP values are in units of gC¨m ´2¨d ´1 [47].For consistency with the time interval of the MODIS GPP products, we then aggregated eight-day flux derived GPP from the calculated daily values.

Experimental Design
To differentiate the effects of the meteorological elements and LAI/FPAR on the GPP estimates, we designed the following four experimental groups to recalculate GPP under different scenarios.
(1) Following the standard MOD17A2 products completely, NCEP-DOE Reanalysis II meteorology data and MOD15A2 LAI/FPAR products were used to estimate GPP (MOD_NCEP GPP); (2) Site meteorology data and MOD15A2 LAI/FPAR products were used to estimate GPP (MOD_Tower GPP); (3) NCEP-DOE meteorology data and GLASS LAI/FPAR products were used to estimate GPP (NCEP_GLASS GPP); (4) Site meteorology data and GLASS LAI/FPAR products were used to estimate GPP (Tower_GLASS GPP).Finally, all the GPP estimates were directly compared with the EC flux derived GPP (Flux_Tower GPP) to assess their accuracies and uncertainties.

Analytical Methods
In this study, the computation results for eight-day, seasonal and annual GPP had been analyzed, respectively.We firstly calculated the daily GPP, and then integrated the daily GPP into eight-day, seasonal and annual GPP by accumulating.The sensitivity of MOD17A2 algorithm to meteorological data and LAI/FPAR products were examined by introducing site meteorological measurements and improved GLASS LAI products.We also assessed the potential error contributions from land cover and ε max .Four statistic indices were used to evaluate the accuracies and uncertainties: (1) the coefficient of determination (R 2 ), which represents the amount of variation in the observations explained by the simulations; (2) bias, which quantifies the difference between the simulations and the observations; (3) the root mean square error (RMSE, %), which is computed as where Ó is the average value of the observed data, and S i and O i are the simulated and observed values, respectively; and Equation (4) the relative error (RE, %), which is
From the perspective of a single site, in CBS, QYZ and YC, all four algorithms (MOD_NCEP GPP, MOD_Tower GPP, NCEP_GLASS GPP, and Tower_GLASS GPP) underestimated the Flux_Tower GPP to varying degrees (Figure 3, Table 2).The MOD_Tower GPP were the most underestimated (´13 gC¨m ´2¨8-day ´1 for CBS, ´18 gC¨m ´2¨8-day ´1 for QYZ and ´22 gC¨m ´2¨8-day ´1 for YC), and the NCEP_GLASS GPP were the least underestimated (´1 gC¨m ´2¨8-day ´1 for CBS, ´2 gC¨m ´2¨8-day ´1 for QYZ and ´14 gC¨m ´2¨8-day ´1 for YC).However, for NMG, all four algorithms overestimated the Flux_Tower GPP.The MOD_NCEP GPP were the most overestimated (3 gC¨m ´2¨8-day ´1), and the NCEP_GLASS GPP and Tower_GLASS GPP were the least overestimated (1 gC¨m ´2¨8-day ´1).For the other sites, some algorithms overestimated the Flux_Tower GPP, whereas some underestimated them.For CBS, QYZ and DHS, the Tower_GLASS GPP algorithm most effectively estimated the Flux_Tower GPP (R 2 = 0.92, R 2 = 0.87 and R 2 = 0.75, respectively), and the greatest improvements were obtained over the MOD_NCEP GPP algorithm (R 2 = 0.76, R 2 = 0.52 and R 2 = 0.30) in particular.However, for YC and NMG, the Tower_GLASS GPP algorithm was not as effective as the other three algorithms for estimating the Flux_Tower GPP.In the case of YC, the performance of the Tower_GLASS GPP algorithm (R 2 = 0.77) was close to that of the most effective algorithm (MOD_NCEP GPP, R 2 = 0.81).However, for NMG, the gap between the Tower_GLASS GPP (R 2 = 0.64) and the best estimating algorithm (MOD_Tower GPP, R 2 = 0.75) was large.For HB and DX, all four algorithms effectively estimated the Flux_Tower GPP.In particular, the MOD_NCEP GPP (R 2 = 0.95, Bias = 0 gC¨m ´2¨8-day ´1) and MOD_Tower GPP (R 2 = 0.89, Bias = 0 gC¨m ´2¨8-day ´1) algorithms performed best for HB and DX, respectively.Regarding XSBN, however, no algorithm succeeded in estimating the Flux_Tower GPP.From the perspective of a single site, in CBS, QYZ and YC, all four algorithms (MOD_NCEP GPP, MOD_Tower GPP, NCEP_GLASS GPP, and Tower_GLASS GPP) underestimated the Flux_Tower GPP to varying degrees (Figure 3, Table 2).The MOD_Tower GPP were the most underestimated

Seasonal GPP
The tower-based GPP measurements (Flux_Tower GPP) were compared with the results achieved by the four algorithms (i.e., MOD_NCEP GPP, MOD_Tower GPP, NCEP_GLASS GPP, and Tower_GLASS GPP) over each season (Figures 4-7 and Table 3).Overall, the four algorithms effectively estimated the Flux_Tower GPP in the spring, autumn and winter.In particular, the algorithms had the best estimating effectiveness for the winter (the values of R 2 ranged from 0.87 to 0.92).However, the overall estimating effectiveness of the four algorithms was poor for the summer (values of R 2 ranged from 0.36 to 0.52).The details are as follows: (1) Spring (March-May): The estimation performances of MOD_NCEP GPP (R = 0.64) and Tower_GLASS GPP (R 2 = 0.65) were excellent (Figure 4, Table 3).The MOD_NCEP GPP, MOD_Tower GPP and Tower_GLASS GPP underestimated the Flux_Tower GPP to varying degrees (´39 gC¨m ´2¨3-month ´1, ´122 gC¨m ´2¨3-month ´1 and ´80 gC¨m ´2¨3-month ´1, respectively), primarily because of the considerable underestimation for QYZ and YC (Figure 4a,b,d).

Annual GPP
The tower-based GPP measurements (Flux_Tower GPP) were compared with the results of the four algorithms (i.e., MOD_NCEP GPP, MOD_Tower GPP, NCEP_GLASS GPP, and Tower_GLASS GPP) during 2003-2005.(Figure 8, Table 3).For all the sites, the annual average value of the Flux_Tower GPP was 1331 gC¨m ´2¨year ´1, and those of the four algorithms were 1170 gC¨m ´2¨year ´1, 827 gC¨m ´2¨year ´1, 1410 gC¨m ´2¨year ´1 and 972 gC¨m ´2¨year ´1, respectively.The NCEP_GLASS GPP algorithm accurately estimated the measurements for CBS, QYZ and HB (Figure 8c) and therefore had the smallest bias (79 gC¨m ´2¨year ´1).However, for the MOD_Tower GPP and Tower_GLASS GPP algorithms, the Flux_Tower GPP for CBS, HB, QYZ, XSBN and YC were severely underestimated (Figure 8b,d), which led to large biases (´504 gC¨m ´2¨year ´1 for MOD_Tower GPP and ´359 gC¨m ´2¨year ´1 for Tower_GLASS GPP).In terms of the estimating effectiveness, the MOD_NCEP GPP algorithm had the greatest performance because it yielded the largest R 2 (0.76) and smallest RMSE (31%) and RE (34%).Significance levels: ** p < 0.01.

Meteorology
In this paper, the MOD_NCEP GPP algorithm followed the standard MOD17A2 algorithm, which used the NCEP-DOE Reanalysis II meteorology data.However, the NCEP-DOE Reanalysis II meteorology data were replaced with the site meteorology data for use with the MOD_Tower GPP algorithm.Thus, the MOD_NCEP GPP and MOD_Tower GPP algorithms were compared to study the impacts of the meteorological data on the estimation results.An analysis of the seasonal and annual GPP estimates showed that the MOD_Tower GPP algorithm with the site meteorology data was inferior to the MOD_NCEP GPP algorithm (Table 3).However, in the case of the eight-day GPP estimates, the estimating effectiveness of the MOD_Tower GPP algorithm was better than that of the MOD_NCEP GPP algorithm, and the improvements were seen at most of the sites (Table 2).
The correlations between the daily NCEP-DOE Reanalysis II meteorology data and the daily site meteorology data were analyzed (Figure 9).The results demonstrated significant correlations between the two types of data in terms of T avg (R 2 = 0.92) and T min (R 2 = 0.90) (Figure 9a,b).This meant that the meteorological reanalysis datasets could be used to effectively estimate the on-site temperature variations.However, in the case of VPD avg , the correlation between the NCEP-DOE Reanalysis II meteorology data and the site meteorology data was poor (R 2 = 0.31) (Figure 9c) because the meteorological reanalysis datasets failed to estimate the on-site moisture conditions.Regarding PAR, the NCEP-DOE Reanalysis II meteorology data commonly overestimated the site meteorology data (Figure 9d).The correlation (R 2 = 0.62) was stronger than that of VPD avg but weaker than those of T avg and T min .
Remote Sens. 2016, 8, 395 15 of 24 temperature variations.However, in the case of VPDavg, the correlation between the NCEP-DOE Reanalysis II meteorology data and the site meteorology data was poor (R 2 = 0.31) (Figure 9c) because the meteorological reanalysis datasets failed to estimate the on-site moisture conditions.Regarding PAR, the NCEP-DOE Reanalysis II meteorology data commonly overestimated the site meteorology data (Figure 9d).The correlation (R 2 = 0.62) was stronger than that of VPDavg but weaker than those of Tavg and Tmin.

LAI and FPAR
Unlike the MOD_NCEP GPP algorithm, the NCEP_GLASS GPP algorithm followed the MOD17A2 algorithm by substituting GLASS_FPAR (i.e., the FPAR derived from GLASS LAI datasets) for MODIS_FPAR.Therefore, comparing the MOD_NCEP GPP and NCEP_GLASS GPP algorithms was helpful for understanding the influences of FPAR on the estimating effectiveness.Similar to the MOD_Tower GPP algorithm, the NCEP_GLASS GPP algorithm was outperformed by the MOD_NCEP GPP algorithm in terms of the annual GPP estimates (Table 3).However, regarding the seasonal estimates, the effectiveness of the NCEP_GLASS GPP algorithm for estimating the Flux_Tower GPP was better than that of the MOD_Tower GPP algorithm, which was vastly inferior to the MOD_NCEP GPP algorithm.In particular, no algorithm was more accurate than the NCEP_GLASS GPP algorithm for estimating the Flux_Tower GPP in the summer.In addition, the NCEP_GLASS GPP algorithm achieved substantial improvements over the MOD_NCEP GPP algorithm for the eight-day estimates (Figure 2a,c).It estimated the Flux_Tower GPP for CBS, QYZ and DHS far more effectively than the MOD_NCEP GPP algorithm (Table 2).

LAI and FPAR
Unlike the MOD_NCEP GPP algorithm, the NCEP_GLASS GPP algorithm followed the MOD17A2 algorithm by substituting GLASS_FPAR (i.e., the FPAR derived from GLASS LAI datasets) for MODIS_FPAR.Therefore, comparing the MOD_NCEP GPP and NCEP_GLASS GPP algorithms was helpful for understanding the influences of FPAR on the estimating effectiveness.Similar to the MOD_Tower GPP algorithm, the NCEP_GLASS GPP algorithm was outperformed by the MOD_NCEP GPP algorithm in terms of the annual GPP estimates (Table 3).However, regarding the seasonal estimates, the effectiveness of the NCEP_GLASS GPP algorithm for estimating the Flux_Tower GPP was better than that of the MOD_Tower GPP algorithm, which was vastly inferior to the MOD_NCEP GPP algorithm.In particular, no algorithm was more accurate than the NCEP_GLASS GPP algorithm for estimating the Flux_Tower GPP in the summer.In addition, the NCEP_GLASS GPP algorithm achieved substantial improvements over the MOD_NCEP GPP algorithm for the eight-day estimates (Figure 2a,c).It estimated the Flux_Tower GPP for CBS, QYZ and DHS far more effectively than the MOD_NCEP GPP algorithm (Table 2).
Comparisons were made between GLASS_LAI and MODIS_LAI as well as GLASS_FPAR and MODIS_FPAR (Figure 10 and Table 5).The results showed extremely significant correlations for the average values of LAI (R 2 = 0.93) and FPAR (R 2 = 0.91).MODIS_FPAR underestimated GLASS_FPAR in all of the eight sites except in NMG and DX.In particular, CBS, HB and DX exhibited the strongest correlations between MODIS_FPAR and GLASS_FPAR, and DHS, QYZ and XSBN exhibited relatively poor correlations (Figure 10b).From the perspective of the biomes, the correlations between MODIS_FPAR and GLASS_FPAR were strong for MF and poor for EBF.Comparisons were made between GLASS_LAI and MODIS_LAI as well as GLASS_FPAR and MODIS_FPAR (Figure 10 and Table 5).The results showed extremely significant correlations for the average values of LAI (R 2 = 0.93) and FPAR (R 2 = 0.91).MODIS_FPAR underestimated GLASS_FPAR in all of the eight sites except in NMG and DX.In particular, CBS, HB and DX exhibited the strongest correlations between MODIS_FPAR and GLASS_FPAR, and DHS, QYZ and XSBN exhibited relatively poor correlations (Figure 10b).From the perspective of the biomes, the correlations between MODIS_FPAR and GLASS_FPAR were strong for MF and poor for EBF.

Land Cover
The four algorithms used the MCD12Q1 Land Cover Classification Type 2 as their land cover classification scheme.The results of the land cover classification directly determine the value of the maximum light use efficiency (εmax) and further influence the estimation results.Based on the user's guide for flux data provided by ChinaFLUX, we obtained the actual land cover type of each site.By comparing the MCD12Q1-based classification results with the actual land cover types at all eight

Land Cover
The four algorithms used the MCD12Q1 Land Cover Classification Type 2 as their land cover classification scheme.The results of the land cover classification directly determine the value of the maximum light use efficiency (ε max ) and further influence the estimation results.Based on the user's guide for flux data provided by ChinaFLUX, we obtained the actual land cover type of each site.By comparing the MCD12Q1-based classification results with the actual land cover types at all eight sites, we found that MCD12Q1 correctly classified most of the sites (Figure 11).The only mistake was that MCD12Q1 misclassified the QYZ site to the MF category (the actual land cover was ENF).
Remote Sens. 2016, 8, 395 17 of 24 sites, we found that MCD12Q1 correctly classified most of the sites (Figure 11).The only mistake was that MCD12Q1 misclassified the QYZ site to the MF category (the actual land cover was ENF).

Light Use Efficiency
For the MODIS GPP algorithm, the values of the maximum light use efficiency (εmax) depend on the types of land cover.Each land cover type corresponds to a constant value of εmax.Equations ( 1) and (3) in Section 2.1 were utilized to obtain the inferred εmax for each land cover type, which were then compared with MOD17A2 εmax directly (Figure 12 and Table 6).The results show that the inferred εmax were higher than MOD17A2 εmax by 60% for MF, 74% for ENF, 11% for EBF, 143% for Crop and 37% for Grass, respectively.However, in the case of a single site, the inferred εmax were not invariably higher than MOD17A2 εmax.For example, the inferred εmax for DHS and NMG were less than the MOD17A2 εmax.

Light Use Efficiency
For the MODIS GPP algorithm, the values of the maximum light use efficiency (ε max ) depend on the types of land cover.Each land cover type corresponds to a constant value of ε max .Equations ( 1) and (3) in Section 2.1 were utilized to obtain the inferred ε max for each land cover type, which were then compared with MOD17A2 ε max directly (Figure 12 and Table 6).The results show that the inferred ε max were higher than MOD17A2 ε max by 60% for MF, 74% for ENF, 11% for EBF, 143% for Crop and 37% for Grass, respectively.However, in the case of a single site, the inferred ε max were not invariably higher than MOD17A2 ε max .For example, the inferred ε max for DHS and NMG were less than the MOD17A2 ε max .

Discussion
In this paper, GPP measurements from eight EC flux towers (Flux_Tower GPP) were used to verify four groups of GPP values calculated from the MODIS GPP algorithms (i.e., MOD_NCEP GPP, MOD_Tower GPP, NCEP_GLASS GPP and Tower_GLASS GPP).The results show that the four groups of GPP values deviated from the tower-based GPP values to varying degrees.The extent of the deviations varied depending on the time interval, the site and the biome type used in the calculation.The deviations can be attributed to many potential factors.Clearly, the accuracy of each parameter input in the MODIS GPP formula will influence the accuracy of the GPP estimates.That is, the meteorology, LAI/FPAR, land cover and εmax data are all error sources [48].Furthermore, the accuracies of the tower-based GPP measurements and the reasonableness of the MODIS GPP algorithm itself should also be taken into account [49].

Discussion
In this paper, GPP measurements from eight EC flux towers (Flux_Tower GPP) were used to verify four groups of GPP values calculated from the MODIS GPP algorithms (i.e., MOD_NCEP GPP, MOD_Tower GPP, NCEP_GLASS GPP and Tower_GLASS GPP).The results show that the four groups of GPP values deviated from the tower-based GPP values to varying degrees.The extent of the deviations varied depending on the time interval, the site and the biome type used in the calculation.The deviations can be attributed to many potential factors.Clearly, the accuracy of each parameter input in the MODIS GPP formula will influence the accuracy of the GPP estimates.That is, the meteorology, LAI/FPAR, land cover and ε max data are all error sources [48].Furthermore, the accuracies of the tower-based GPP measurements and the reasonableness of the MODIS GPP algorithm itself should also be taken into account [49].

Impact of Meteorology on MODIS GPP
The PAR in the MODIS GPP algorithm together with T min and VPD avg , which influence LUE, are all meteorological factors.Thus, the meteorological inputs may be the largest sources of error in the GPP estimates [50].The NCEP-DOE Reanalysis II meteorology dataset used in the C5.5 version of MOD17A2 product has a low resolution.Although a strict interpolation was used to make the resolution comparable to the MOD15A2 product, the errors caused by the original low resolution could not be eliminated completely.In addition, the interpolation process produces new uncertainties.From comparisons between the meteorological reanalysis data and the site meteorology data, it can be clearly observed that the NCEP-DOE Reanalysis II meteorology data are ineffective for estimating PAR, and the outcome is even worse for VPD (Figure 9c,d).In fact, the MODIS GPP algorithm does not directly take the soil moisture into account.However, soil moisture plays a key role in GPP estimation [51].To compensate for this defect, the VPD was added to the algorithm as a proxy for the soil moisture; however, errors will be unavoidably introduced.In particular, the estimation error will increase when the deviations between the VPD estimates and measurements are large.
In this paper, we replaced the NCEP-DOE Reanalysis II meteorology data with the site meteorology data to recalculate the MODIS GPP formula and obtained new GPP values (i.e., MOD_Tower GPP).Unexpectedly, the GPP estimates over eight days calculated using the site meteorology data did not result in obvious improvements in the accuracy of the GPP estimates (Figure 2a,b).On the contrary, the estimates of the seasonal and annual GPP calculated using the site meteorology data were not as accurate as those calculated using the meteorological reanalysis data (Table 3).This implies that improved meteorological inputs will not necessarily enhance the estimation effectiveness of the MODIS GPP algorithm, highlighting the need to analyze other sources of error.

Impact of LAI/FPAR on MODIS GPP
In addition to improving the quality of the meteorological inputs, we also introduced an improved LAI dataset (GLASS LAI) to evaluate the impacts of LAI/FPAR on the GPP estimates (i.e., NCEP_GLASS GPP).The results show that after replacing MODIS FPAR with GLASS FPAR, substantial improvements were achieved in estimating the eight-day GPP values (Figure 2a,c).The effectiveness in estimating GPP during the summer was also improved (Figure 5a,c).For Cropland and Grassland, the GPP estimates became more accurate after FPAR was improved (Table 4).Based on the comparisons, it was found that the MODIS LAI/FPAR were generally less than the GLASS LAI/FPAR (Table 5).Indeed, the MOD17A2 GPP algorithm generally underestimated the tower-based GPP (Table 3).Therefore, the improved LAI/FPAR compensated for the underestimation caused by MODIS LAI/FPAR to some extent, which resulted in a higher effectiveness of the GPP estimation.
In this paper, we further estimated GPP by improving the quality of the meteorological inputs and the LAI/FPAR data jointly (i.e., Tower_GLASS GPP).We simply replaced the NCEP-DOE Reanalysis II meteorology data with the site meteorology data and replaced MODIS FPAR with GLASS FPAR.
The results demonstrate that the replacements resulted in the greatest improvements in estimating the eight-day GPP values (Figure 2a,d).However, this type of replacement strategy was not very effective in estimating the seasonal and annual GPP (Table 3).This implies that the GPP estimation accuracy can be greatly improved by simultaneously enhancing the quality of the multiple inputs of the MODIS GPP algorithm.However, it does not mean that improvements can be achieved simply by substituting the input parameters because the overall performance gains are by no means the sum of the individual gains.In the future, we will thoroughly study how to use the synergies among the individual gains to increase GPP estimation accuracy.

Impact of Land Cover on MODIS GPP
The impact of land cover classification on GPP estimation is non-negligible.In this paper, the land cover at site QYZ was of type ENF, but it was misclassified to type MF by MCD12Q1 (Figure 11).This misclassification directly led to the misestimation of the FPAR by MOD15A2 and the misjudgment of ε max by BPLUT and thus affected the MOD17A2 GPP estimate for site QYZ.

Impact of LUE on MODIS GPP
Regarding the MODIS GPP algorithm, ε max represents the maximum light use efficiency of the vegetation in photosynthesis.The value of ε max will vary with the type of biome.For a given biome type, ε max is set to a constant by BPLUT.Encouragingly, BPLUT corrected and updated the value of ε max several times.However, because of the immense diversity of earth surface environments and climate conditions, assigning a constant value of ε max to the same biome type does not conform to the truth [52].The improved ε max that we inferred for the eight sites in this paper greatly differed from the MOD17A2 ε max that was specified by BPLUT (Figure 12 and Table 6).Because the value of ε max provides the foundation for the actual light use efficiency, the misestimation of ε max will inherently decrease the accuracies of the GPP estimates.

Uncertainties, Errors, and Accuracies
Note that our evaluation of MODIS GPP is based on the assumption that the GPP values measured by the EC flux tower are the ground truth.However, many uncertainties exist in the tower-based GPP measurements [53].The tower GPP data used in this paper were calculated as the difference between the net ecosystem exchange (NEE) and ecosystem respiration (R eco ).The precise estimation of R eco is difficult and can lead to systematic and random errors in estimating GPP.Furthermore, uncertainties also arise due to scale mismatches between the tower flux footprints and MODIS pixels.In addition, the MODIS GPP algorithm itself is also a potential error source in GPP estimation.Recently, many studies have examined the structural errors of the MODIS GPP algorithm [54][55][56].For example, Zhang et al. (2012) [57] compared the MODIS GPP product with estimates from a two-leaf process-based model.Their results showed that the MODIS GPP algorithm cannot properly treat the contribution of shaded/sunlit leaves to the calculation of the total GPP.
We should also notice that the number of observations in the regression analysis could affect our estimation results.The seasonal and annual GPP reduced a lot the number of observations in relation to the eight-day data.This could be one reason why the eight-day data gave, in general terms, more accurate estimations.

Conclusions
In this study, we methodically evaluated the eight-day, seasonal and annual MODIS GPP using EC flux measurements at eight sites in five various biome types across China from 2003 to 2005.The sensitivity of MOD17A2 algorithm to meteorological data and LAI/FPAR products were examined by introducing site meteorological measurements and improved GLASS LAI products.We also assessed the potential error contributions from land cover and ε max .Each of these validation steps would help to isolate and identify sources of error.The main conclusions can be summarized as follow: (1) The standard MOD17A2 product (i.e., MOD_NCEP GPP) performed better at estimating the annual GPP (R 2 = 0.76) and agreed well with the tower GPP during the autumn (R 2 = 0.80) and winter (R 2 = 0.92).However, the effectiveness of the MOD_NCEP GPP algorithm for estimating GPP over eight days was poor (R 2 = 0.55) and even worse during the summer (R 2 = 0.46).In addition, the MOD_NCEP GPP algorithm was ineffective when estimating the annual GPP for mixed forests, evergreen needleleaf forests and cropland.(2) Replacing the NCEP-DOE Reanalysis II meteorology data with the site meteorology data (i.e., MOD_Tower GPP) only slightly improved the correlation with the tower GPP over eight days (R 2 = 0.56).However, substantial improvements in estimating the tower GPP over eight days (R 2 = 0.65) and during the summer (R 2 = 0.52) were achieved by substituting GLASS_FPAR for MODIS_FPAR (i.e., NCEP_GLASS GPP).For cropland, the GPP estimates were more accurate after the FPAR data were improved.When the meteorology inputs and FPAR data were simultaneously replaced with improved data (i.e., Tower_GLASS GPP), the effectiveness in estimating the tower GPP was improved significantly for mixed forests and evergreen needleleaf forests.(3) There are four potential error sources related to the inputs of the MOD17A2 algorithm: meteorology, LAI/FPAR, land cover and ε max .The NCEP-DOE Reanalysis II meteorology data failed in estimating the tower measured VPD (R 2 = 0.31), and MODIS_FPAR underestimated the improved FPAR data at most sites.Although MCD12Q1 succeeded in classifying most of the sites correctly, the values of MOD17A2 ε max were much smaller than the optimized ε max values for all five biome types discussed in this paper.
From above analysis, we suggest that the qualities of the meteorological data and LAI/FPAR products need to be improved and the BPLUT parameters should be adjusted to provide better GPP estimates using MOD17A2 for Chinese ecosystems.In future research, additional high-resolution LAI/FPAR and GPP products should be considered, such as MOD15A2H and MOD17A2H, which have been recently released.

Figure 3 .
Figure 3.Time series of eight-day GPP derived from the tower estimates (i.e., Flux_Tower GPP) and four different experimental groups (i.e., MOD_NCEP GPP, MOD_Tower GPP, NCEP_GLASS GPP and Tower_GLASS GPP) at: (a) CBS forest site; (b) QYZ forest site; (c) DHS forest site; (d) XSBN forest site; (e) YC cropland site; (f) HB grassland site; (g) NMG grassland site; and (h) DX grassland site.The full name for each site is listed in Table1.

Figure 3 .
Figure 3.Time series of eight-day GPP derived from the tower estimates (i.e., Flux_Tower GPP) and four different experimental groups (i.e., MOD_NCEP GPP, MOD_Tower GPP, NCEP_GLASS GPP and Tower_GLASS GPP) at: (a) CBS forest site; (b) QYZ forest site; (c) DHS forest site; (d) XSBN forest site; (e) YC cropland site; (f) HB grassland site; (g) NMG grassland site; and (h) DX grassland site.The full name for each site is listed in Table1.

Figure 9 .
Figure 9. Scatterplots of the daily NCEP-DOE II meteorology against the daily tower meteorological data for: (a) T avg ; (b) T min ; (c) VPD avg ; and (d) PAR.Significance levels: ** p < 0.01.

Figure 10 .
Figure 10.Comparison of the average (a) LAI and (b) FPAR from the MODIS datasets and the GLASS LAI products.Significance levels: ** p < 0.01.

Figure 10 .
Figure 10.Comparison of the average (a) LAI and (b) FPAR from the MODIS datasets and the GLASS LAI products.Significance levels: ** p < 0.01.

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

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

Table 2 .
Comparison of eight-day GPP derived from the tower estimates with those derived from the MODIS algorithms (MOD_NCEP GPP, MOD_Tower GPP, NCEP_GLASS GPP, and Tower_GLASS GPP) at eight sites.The full name for each site is listed in Table1.

Table 3 .
Comparison of the seasonal and annual GPP derived from the tower estimates with those derived from the MODIS algorithms (MOD_NCEP GPP, MOD_Tower GPP, NCEP_GLASS GPP, and Tower_GLASS GPP).

Table 3 .
Comparison of the seasonal and annual GPP derived from the tower estimates with those derived from the MODIS algorithms (MOD_NCEP GPP, MOD_Tower GPP, NCEP_GLASS GPP, and Tower_GLASS GPP).

Table 4 .
Biome-based comparison of the annual GPP derived from the tower estimates with those derived from the MODIS algorithms (MOD_NCEP GPP, MOD_Tower GPP, NCEP_GLASS GPP, and Tower_GLASS GPP).

Table 5 .
Comparison of the LAI and FPAR derived from the GLASS LAI datasets with those from the MODIS products for various biomes.

Table 5 .
Comparison of the LAI and FPAR derived from the GLASS LAI datasets with those from the MODIS products for various biomes.

Table 6 .
Comparison of the inferred εmax and MOD17A2 εmax at the eight sites.

Table 6 .
Comparison of the inferred ε max and MOD17A2 ε max at the eight sites.