Evaluating Parameter Adjustment in the MODIS Gross Primary Production Algorithm Based on Eddy Covariance Tower Measurements

How well parameterization will improve gross primary production (GPP) estimation using the MODerate-resolution Imaging Spectroradiometer (MODIS) algorithm has been rarely investigated. We adjusted the parameters in the algorithm for 21 selected eddy-covariance flux towers which represented nine typical plant functional types (PFTs). We then compared these estimates of the MOD17A2 product, by the MODIS algorithm with default parameters in the Biome Property Look-Up Table, and by a two-leaf Farquhar model. The results indicate that optimizing the maximum light use efficiency (εmax) in the algorithm would improve GPP estimation, especially for deciduous vegetation, though it could not compensate the underestimation during summer caused by the one-leaf upscaling strategy. Adding the soil water factor to the algorithm would not significantly affect performance, but it could make the adjusted εmax more robust for sites with the same PFT and among different PFTs. Even with adjusted parameters, both one-leaf and two-leaf models would not capture seasonally photosynthetic dynamics, thereby we suggest that further improvement in GPP estimaiton is required by taking into consideration seasonal variations of the key parameters and variables. OPEN ACCESS Remote Sens. 2014, 6 3322


Introduction
Gross primary production (GPP), which is the amount of light energy from the sun converted to chemical energy, determines the thermal, water and biogeochemical cycles in terrestrial ecosystems.However, even when eddy covariance (EC) flux data and remote sensing data are conjunctively introduced into various diagnostic models, uncertainties in modeled GPP are still vast.Current estimates of global GPP range between 102 and 165 petagrams of carbon per year, where uncertainties are likely from errors of input data, poor parameterization and model structure [1][2][3][4][5][6].Among the methods of reducing uncertainties in GPP simulation, parameter adjustment and structural modification are adopted, and the former method could compensate for the errors introduced by a model structure [2,7,8].
Many concepts were developed to simulate carbon assimilation [9][10][11].One of the widely accepted approaches is the light use efficiency model because of both its simple structure-which assumes that a fraction of the photosynthetically active radiation (PAR) absorbed by the vegetation canopy is used for plant primary production [9]-and its large amount of available input data, including EC measurements [12][13][14][15] and remotely sensed data [16][17][18].One application of this kind of model is the MODIS GPP algorithm [19].Its latest product MOD17A2 Collection 5 (https://lpdaac.usgs.gov/products) is forced by the National Center for Environmental Prediction-Department of Energy (NCEP-DOE) reanalysis II data and the default parameters are derived from the Biome Property Look-Up Table (BPLUT) [20].
Based on the algorithm and its product, many evaluations have been made.Early studies focused on uncertainties introduced by input data [5,16,21] and simulations using different input data [17,22,23].A direct comparison between the NCEP-DOC reanalysis II data and EC tower measurements at 12 African sites suggested that these two forcing data with different spatial resolutions were comparable in air temperature and atmospheric vapor pressure deficit (VPD), but the relationship was scattered in incoming PAR [14].Recently, many studies paid attention to structural errors of the MODIS GPP algorithm.Zhang et al. [24] compared the MODIS GPP product with the estimates using a process-based ecosystem model (Boreal Ecosystem Productivity Simulator).They noted that the MODIS GPP algorithm cannot properly treat the contribution of shaded leaves to canopy-level GPP.He et al. [13,15] developed a two-leaf light use efficiency model for improving the calculation of GPP and validated its application in six ecosystems.Previously, the two-leaf upscaling strategy had been well documented in the Farquhar model, as compared with the one-leaf strategy [25,26].Aside from model restructure, adjusting biome parameters is another way to improve simulation.Evaluations of MODIS GPP using EC data indicated that adjusting the maximum light use efficiency (ε max ), which was derived from the BPLUT [19,20], in the algorithm might be needed to better estimate GPP in both Africa [14] and northern China [27].As a key parameter in the algorithm, ε max was further affected by the scales of daily minimum temperature (T MIN ) and VPD [19].However, the T MIN function was observed not to constrain MOD17A2 GPP much [8,14] and both functions could be corrected using flux tower measurements [8].Moreover, some studies noted that optimizing the algorithm using the soil water content can improve agreement with the measurements [28], but only in dry regions [14,29,30] or throughout the growing season [16].
As mentioned above, both adjusting the key parameters and modifying the model structure can improve GPP estimation using the MODIS algorithm, and the former can compensate for errors introduced by the latter [2,7,8].Therefore, the best approach for improving the algorithm is open to debate, and benefits of parameter adjustment are needed to be validated for multiple biome types across different time scales and over relatively long time periods.Furthermore, less attention has been given to the effects of adopting different methods of parameter adjustment, such as using model-data fusion [31] and adding control factors [2].These issues must be clarified and resolved to reduce uncertainties in GPP estimates on regional and global scales.
We hypothesize that adjusting key parameter can improve estimates and can compensate for structural errors caused by adopting the one-leaf strategy in the MODIS GPP algorithm.The objective of this study is to evaluate capacity of parameter adjustment in the algorithm to improve estimates for multiple plant functional types.We used EC measurements to optimize key parameters in the algorithm.The selected EC towers represented nine plant functional types across six biomes in two main climate zones.Then we compared the adjusted models with the MOD17A2, the algorithm with default parameters, and a two-leaf Farquhar model across half hourly, daily, monthly and seasonal time scales.The Farquhar model is a default component of the Dynamic Land Model (DLM) [32,33].By using an existing dataset, this research proves a solid foundation for evaluating the MODIS algorithm's ability to estimate GPP across multiple plant functional types and a range of time scales.

Eddy-Covariance Data
We used the FLUXNET database (http://fluxnet.ornl.gov/) to calibrate the models and to validate GPP estimates.The dataset contains annual files of half-hourly meteorological and flux data from more than 400 EC sites across Europe (CarboEurope), America (AmeriFlux and Fluxnet-Canada), and Asia (AisaFlux and ChinaFLUX), etc.To reduce potential errors derived from the observations, we selected the EC towers according to the following criteria: (1) the site provides four or more years of continuous data as a part of the publicly accessible standardized Level 4 or 3 database; (2) a -site-year‖ is accepted for analysis if more than 90% of the half hours in a year contain non-missing values for the meteorological data (downwelling solar radiation, precipitation, wind speed, air temperature and relative humidity), the carbon flux data (net ecosystem exchange and ecosystem respiration (R eco )), and the energy fluxes (net radiation (R n ), ground heat flux (G), latent heat flux (LE) and sensible heat flux (H)); and (3) energy balance closure is evaluated for each site-year according to the ratio of the dependent flux variables (H + LE) against the independently derived available energy (R n − G) for each half hour [34].The values of the half-hourly energy balance closure ratio (H + LE)/(R n − G) deviate from the ideal closure (a value of 1) because random error exists, and the magnitude of the CO 2 uptake is less when the energy imbalance is greater [34,35].Thus we recorded the number of daytime half hours of which the ratio was in the range of 0.6-1.4,and then accepted a -site-year‖ when the accumulated number was greater than 60% of the total number of daytime half hours during the growing season.Uncertainties in the EC-measured GPP still exist because of the underestimation of R eco at night, gap filling algorithm uncertainty, partitioning uncertainty, random uncertainty, and threshold friction velocity uncertainty [36].We considered the supplied GPP as the -ground truth‖ [18] in this study.
Finally, 102 site-years, which represent six biome types across two main climatic environment zones (i.e., plant functional types, PFTs) [37] at 21 EC sites , were selected (Table 1).Seven of these sites are in boreal regions, and fourteen sites are in temperate regions.Five sites that have Mediterranean climates were characterized as being in the temperate zone because of the absence of Mediterranean forest in the PFTs we used [37].Though the three sites, CA-Ca1, CA-Ca2 and CA-Ca3, were located in adjacent areas, their planting years were 1949, 2000 and 1988, respectively, indicating different tree ages.This is a similar case with CA-Ojp and CA-Obs; their planting years were 1929 and 1879, respectively.We used two years of each site for model calibration, and another two consecutive years for validation.

MODIS 8-Day Average GPP Product
The MODIS GPP product MOD17A2 Collection 5 was designed to provide an 8-day average measure of the global terrestrial vegetation using MODIS land cover, vegetation product and surface meteorology at 1 km resolution [20,65].We used mean values of the 3-by-3 pixels, with the center pixel containing the tower location (Table 1).Because the footprint radius of most annual cumulative climatology data ranged from 0.70 to 1.5 km [66], the 3-by-3 pixel area is expected to represent the flux tower footprint well [18,67,68].

MODIS GPP Algorithm
MOD17A2 is calculated using a light use efficiency model with the one-leaf upscaling strategy based on the radiation conversion efficiency concept of Monteith [9] as follows [19,20,69]: where PAR is the incident photosynthetically active radiation-its value is assumed to be 45% of the incident shortwave radiation [19]; ε max is the maximum light use efficiency, which depends on the plant functional types; f(VPD) and f(T air ) are the scalars of vapor pressure deficit and air temperature, respectively, both of which are ranged from 0 to 1 to downscale the maximum light use efficiency of the canopy to the actual value.The scales were calculated as follows [8,13,15]: where VPD max and T MINmin are daily maximum VPD and daily minimum temperature at which the light use efficiency equals 0, and VPD min and T MINmax are daily minimum VPD and daily minimum temperature at which the light use efficiency is maximum.We replaced the daily minimum temperature in the temperature function [8,13,15] by EC-observed half-hourly temperature (T air ) to calculate GPP on half-hourly time scale.fPAR is the fraction of PAR being absorbed by the canopy, which was estimated using the Beer's law [21,70] as follows: ( where k is the light extinction coefficient, which is set to 0.5 [15,70,71], and LAI is the green area index. The soil moisture scalar (β t ) is another factor that impacts the light use efficiency of vegetation [16,17,29,30,72].We added this factor to Equation (1) in another simulation as follows [37]: (5) where w i is the plant-wilting factor for soil layer i, and r i is the fraction of roots in soil layer i.The soil profile is divided into fifteen layers, for which the depth of the layer increases exponentially with the soil layer number [37].The factor was calculated using DLM, because the soil moisture measurements were not provided by all site-years we selected and the values were simulated well by DLM [33].More details can be found in the Appendix and Oleson et al. [37]

Farquhar Model
The total canopy photosynthesis (A) in the two-leaf Farquhar model was calculated for sunlit and shaded parts by adopting sunlit and shaded leaf area indices (LAI sun and LAI sha ) separately, following [73]: (6) The net CO 2 assimilation rate at leaf level was expressed as follows [37,73]: min (7) where , , and are the Rubisco-limited rate, the light-limited rate, the export-limited rate and leaf dark respiration of sunlit or shaded leaf (i = 1 or 2), respectively.The Rubisco-limited rate was controlled by the soil water factor (Equation (A3)).Details can be found in Appendix.

Forcing Data
Off-line single point simulations with a 30 min time step were performed using observed meteorological data and land-surface data.Half-hourly meteorological data, including downwelling solar radiation (in W•m −2 ), precipitation (in mm), wind speed (in m•s −1 ), air temperature (in K), and relative humidity (in %), were measured at the EC towers.For these key model inputs, missing half-hourly values, which were due to periods of instrument failure, were gap-filled by linear interpolation for gaps of less than 2 hours.Larger gaps were filled by applying a simple interpolation technique of mean diurnal variation [74,75].Monthly LAI values for each site were extracted from a global LAI map based on 10-day synthesis VEGETATION images with 1 km spatial resolution taken in 2003.The values had been corrected based on a global clumping index map produced from the multi-angle observation of the POLDER 1, 2 and 3 sensors [4,76,77].We further corrected monthly LAI for each site using the ratio of the LAI max value (Table 1) against the extracted LAI value within the same month, in which LAI max was supplied by the biological information for each site.
Each site for the offline simulations using DLM were initialized by spinning-up for 200 years with repeat years using 1982-2001 atmospheric forcing dataset from the National Centers for Environmental Prediction reanalysis dataset [78] provided by National Center for Atmospheric Research.Although the years for which available supplementary land-surface data are available do not always correspond to the years being modeled, we assumed that the data are adequate for our photosynthesis modeling.We only utilized the biogeophysical module in DLM, thus the estimation was unaffected by biogeochemical (e.g., carbon-nitrogen coupling) uncertainties [2,79].

Parameter Selection and Optimization
We optimized some key biome-dependent parameters regarding carbon assimilation in the LUE, the LUE-SW and the Farquhar simulations.Adopting the parameter optimization algorithm by Chen et al. [38], we first identified which parameters are most sensitive to photosynthesis by randomly sampling parameters within their possible ranges and analyzing the response.The maximum light use efficiency and the leaf maximum carboxylation rate at 25 °C constrained by leaf nitrogen, which are significantly sensitive in the MODIS GPP algorithm and the Farquhar model, respectively, were selected to be optimized.Then, we applied the ensemble Kalman filter data-model synthesis approach, which encompasses both model parameter optimization and data assimilation, to optimize these parameters by minimizing the difference between observations and predications [80].These selected parameters were optimized at the site level to reduce errors introduced by plant functional type classification [81].To minimize VPD and temperature errors in the MODIS GPP algorithm (Equation ( 1)), we calculated VPD max and VPD min in Equation ( 2) and T MINmax and T MINmin in Equation (3) for each site following Kanniah et al. [8].

Experiment
We performed the following two simulations to document different methods of parameter adjustment applied to the MODIS GPP algorithm: LUE: A simulation with the MODIS GPP algorithm (Equaiton (1)) and optimized biome-parameters using EC measurements; LUE-SW: Addition of the soil water scale (Equation ( 5)) to the MODIS GPP algorithm.The parameters were optimized after the addition.
These two performances were evaluated using the MOD17A2 and the following two simulations forced by meteorological measurements: LUE def : A simulation with the MODIS GPP algorithm and the default biome-parameters supplied by BPLUT [20], which is aimed at testing default parameters forced by EC data; Farquhar: A simulation with the two-leaf Farquhar model (Equations ( 6) and ( 7)) to investigate structural error introduced by the one-leaf upscaling strategy and to validate compensation of parameter adjustment to the MODIS algorithm.

Model Performance
We quantified the model performance using statistical analysis based on the half-hourly GPP for each model-data pair.Model-data mismatch was evaluated using the bias, and the root-mean-square error (RMSE) [82][83][84], which are defined as follows: where P i and O i denote the predicated and observed values, respectively, and is the mean value of the observed data.
A final characterization of model performance uses the Taylor diagram [85], in which a single point indicates the linear correlation coefficient (R) and the ratio of the standard deviations of the prediction and the observation σ norm = σ p/ σ o ), along with the root-mean-square (RMS) difference of the two patterns on a two dimensional plot.An ideal model has a standard deviation ratio of 1.0 and a correlation coefficient of 1.0, i.e., the reference point on the x-axis.The Taylor skill (S) is a single-value summary of a Taylor diagram, where unity indicates perfect agreement with observations.More generally, each point for any arbitrary data group [85,86] can be scored as (10)

Model Parameters Variation
We optimized the key parameters of the LUE, LUE-SW and the Farquhar model (Table 2).For the selected 21 sites, the values of ε max ranged from 0.53 to 1.72 gC•MJ −1 in LUE.The highest ε max was exhibited in the boreal forest of deciduous broadleaf species.Considering the limitation of soil moisture, the optimized ε max in LUE-SW increased by 0.08 gC•MJ −1 on average.However, the statistical differences between the optimized parameters in each simulation (LUE or LUE-SW) and the default values (LUE def ) were not significant (p > 0.0 according to the Fisher's least significant difference test.The two-leaf Farquhar model calculated photosynthesis of sunlit and shaded leaves separately, but used the same values of the leaf maximum carboxylation rate at 25°C constrained by leaf nitrogen for both leaves (Equation (A3)).The average parameter was 37.16 μmol•m −2 •s −1 with a standard deviation of ±9.48 μmol•m −2 •s −1 for 21 sites.The absolute value of the standard deviation accounted for 25.51% of the average, which was less than the percentage of the absolute standard deviation in LUE (31.80% for ε max ) and was comparable with that in LUE-SW (25.56% for ε max ).These comparisons indicate that the key parameters in the LUE-SW and the Farquhar model were robust based on site-specific optimization.
We adopted the one-way analysis of variance (ANOVA) to determine whether there were differences among site-specific parameters in LUE, LUE-SW and the Farquhar model according to biome types or climate zones.As presented in Table 2, the differences were significant (p < 0.05) among nine plant functional types (i.e., biome types + climate zones) for all three simulations, but were not significant among biome types or climate zones, except for the biome-based category in the Farquhar model.These results suggest that it is necessary to specify parameters according to PFTs in both the MODIS algorithm and the Farquhar model to reduce errors introduced by parameter classification in regional simulation.

Model-Data Agreement on Half-Hourly and Daily Time Scales
Two methods of parameter adjustment, optimizing the key parameter ε max and adding the soil water factor, had the same effects on performances of the MODIS GPP algorithm.Examples are presented for one representative year of each site (Table 1) because the behavior is comparable from year to year in each simulation.On half-hourly time scale, σ norm values tended to increase linearly from 0.77 to 0.97 for LUE and from 0.80 to 0.97 for LUE-SW, in R of 0.74-0.92and 0.78-0.92,respectively (Figure 1a), indicating that adding the soil water factor to the algorithm could only improve GPP simulation slightly.The two-leaf simulation could well quantify canopy carbon assimilation, in which the average R and σ norm were 0.882(±0.042)and 1.007(±0.118),respectively, on a half-hourly time scale.LUE and LUE-SW still had same performances after accumulating half-hourly GPP into 8-daily average values (Figure 1b).Note that although the key parameters had been optimized for all three models at each site, the standard deviations of estimates by two one-leaf models were lower than those of observation systematically.The average σ norm values in LUE-SW decreased by 0.102 and 0.162 on half-hourly and 8-day average time scales, respectively, compared with the two-leaf model.These analyses suggest that parameter optimization is available for compensating the error caused by ignoring soil moisture, but is not for the model structural errors caused by the canopy upscaling strategy.Moreover, we quantified differences among MOD17A2, LUE def and three parameter-adjusted simulations (Figure 1b).Site-specific parameters made the simulations more robust than when default values for light use efficiency models were used overall.For most sites we selected, σ norm , which were ranged from 0.55 to 1.12, increased with R in MOD17A2.The excluded sites are the CA-NS1 (site number 10), FR-Pue (site number 12), IT-Cpz (site number 13), DK-Sor (site number 17), CA-Mer (site number 18), and CA-NS6 (site numbers 19), most of which are deciduous ecosystems.There was no consistency in sites with large or small σ norm between LUE def and MOD17A2.For instance, in the CA-N site, σ norm increased by 0.45 from LUE (0.76) to LUE def (1.21), and then by 0.70 from LUE def to MOD17A2 (1.91), which were similar for FR-Pue and CA-NS6.In the ES-ES1 site (site number 5), the value increased by 0.40 from LUE (1.19) to LUE def (1.59), but decreased by 0.25 from LUE to MOD17A2 (0.94).These results implied that poor parameterization is one reason for errors in MOD17A2, but not for all sites.Errors introduced by parameters could be compensated or intensified by uncertainty in meteorological and vegetation data.We further compared the daily simulations with the EC observations using the linear regression analyses for each PFT, as shown in Figure 2, whereas the daily observations were averaged into 8-day values to compare with MOD17A2.The slopes of the MODIS product varied largely from 0.49 to 1.50, and the R 2 ranged from 0.46 to 0.89.Adopting site-specific parameters and input data effectively improved the accuracy of simulations.Further improvements in both correlation and variability were achieved by using the two-leaf strategy.The biases of GPP simulated by the two-leaf Farquhar model were relatively small for nine PFTs, with the slopes of 0.81-1.07and R 2 of 0.67-0.92.Overall, model-data agreement across the selected 21 sites was better for the two-leaf model than for the one-leaf model.Systematic underestimation existed in the one-leaf models for all PFTs, and could not be compensated by adjusting key parameters.

Model-Data Agreement on Monthly and Seasonal Time Scales
We compared monthly and seasonal GPP variation among the simulations to explore the model's responses to varying weather conditions during different seasons.Table 3 shows that uncertainties of four estimates were large in the warm season.Although the MODIS product had relatively small biases from April to September compared with the other simulations, its standard deviation was 18.6(±15.9)times greater than the value of bias on average, and the RMSE ranged from 1.53 to 2.70 gC•m −2 •day −1 .Thus the low bias of MOD17A2 was caused by terms with opposite signs cancelling anthers.The standard deviations of the biases and the RMSE were relatively small in the other three simulations during the warm season, indicating that the biases could be considered as a measure of modeling uncertainties.LUE underestimated GPP for all seasons, especially in summer, when the bias was -1.28(±1.00)gC•m −2 •day −1 .The negative bias in summer was also apparent for LUE-SW, but it was slightly better.The bias in summer was canceled using the two-leaf Farquhar model with reduced RMSE.However, there were both overestimation in spring and underestimation in autumn for the Farquhar model overall, even though the leaf maximum carboxylation rate has been optimized for each site.
Photosynthesis in winter accounts for ~7% of yearly carbon assimilation in the selected 11 evergreen forests, and the contribution is as high as 16.6% at the ES-ES1 site.Thus, we further compared the seasonal GPP simulations of this biome type in Figure 3a.Ignoring the MODIS product, the site-level model-data agreement exhibited a low degree of variability from spring to fall, but the Taylor skill ranged from zero to unity in winter, indicating that photosynthesis dynamic could not be captured well during the cold season.The MODIS product for the evergreen forests had a relatively large Taylor skill, but agreed well with observations decreased for the deciduous forests and shrubs (Figure 3b) and the grasslands (Figure 3c) in all four seasons, especially summer, which contributed to the large RMSE indicated in Table 3.We compared the biases and RMSE of seasonally composite diurnal variations estimated by the LUE-SW and the Farquhar model to find further discrepancies between the simulations (Figure 4).Overall with site-specific parameters (Table 2), both models had similar half-hourly biases and RMSE in all four seasons, excluding obvious underestimations by LUE-SW during the morning and afternoon of summer.The negative biases were improved using the two-leaf Farquhar model with reduced RMSE, especially in the deciduous forests and shrubs.The average biases increased from −35.2(±24.0) to −0.8 ± .μgC•m −2 •s −1 and RMSE decreased from 43.3(±35.8)to 39.0(±32.4)μgC•m −2 •s −1 during daytime (6:00-18:00) for forests and shrubs in summer (Figure 4a,b).However, a significant improvement by the two-leaf model did not perform in the grasslands (Figure 4c).Note that the estimates by both models had great RMSE during the midday of spring and fall for all biome types, and during the midday of winter for evergreen forests.These results suggest that, whether model structures are complex or simple, the photosynthetic models could not well capture change of GPP consistently in spring and fall for forests and shrubs, and seasonal change for grasslands.

Uncertainties in Input Data and Parameters in MOD17A2
Many studies have demonstrated that the MODIS product normally underestimates GPP compared with the EC observations, including savanna and grassland in Africa [14], and forest, grassland and cropland in east Asia [12,13,27,87] and in North America [22,24,88].However, yearly overestimations were also found at some forest sites [8,13,89,90].In this study, overestimated and underestimated GPP were both observed at the selected sites (Figure 2), especially in the deciduous forests and shrubs (Figure 3b).It raises doubts as to the accuracy of the simulated temporal and spatial distributions of the MODIS GPP product at regional or global scale because of errors introduced by the reanalysis of meteorological data [65], the fraction of photosynthetically active radiation (fPAR in Equation ( 1)) [8,14,71], the land cover data [14,91] and the model structure [24,71].
Replacing the reanalysis meteorological and MODIS vegetation data [5,65] with the tower data (LUE def ) did not obviously improve GPP simulations (Figure 1b).Large and small σ norm were both observed in LUE def , and the sites with large errors were not consistent with those in MOD17A2 (Figure 1b).For some sites, such as ES-ES1, CA-Oas, CA-Mer, AT-Neu and IE-Dri, the performances were even worse in LUE def than in MOD17A2.Combining the tower observations with the default parameters in BPLUT [20] could compensate errors introduced by the biome-dependent parameters for some sites, such as ES-ES1, but would increase errors for others, as CA-NS1, FR-Pue and CA-NS6 (Figure 1b).Adjusting parameters using tower data could effectively improve GPP estimate by the MODIS algorithm, but systematic underestimations were observed (Figure 2), thereby indicating model structural errors.

Uncertainties in Parameter Sets in the MODIS GPP Algorithm
We adjusted the key parameter ε max in the models that based on the MODIS GPP algorithm (LUE and LUE-SW) at site level.The values in LUE-SW were greater than those in LUE on average, but both were not significantly different from the default values except for CA-Oas, DK-Sor, CA-Mer and AT-Neu, for which the value of ε max was greater than 1.3 gC•MJ −1 (Table 2).These sites include the biome types of shrub, grassland and deciduous forest.Overall, the optimized ε max values in both models were in the range of 0.55-2.8gC•MJ −1 , as reported by Garbulsky et al. [92], who calculated the global gross radiation use efficiency from the data provided by 35 EC flux sites and the MODIS fPAR data.Many studies have optimized ε max in the MODIS GPP algorithm based on site data [14,27].Those results have also demonstrated that the values for shrub and grassland were larger than those in the look-up table, and ε max reached ~1.56 gC•MJ −1 for deciduous forests [93].
Our values of ε max are acceptable, but an overestimated ε max could be induced by an underestimated fPAR value [14].In this study, fPAR was calculated based on a function of LAI and k (Equation ( 2)).Although we corrected the input LAI based on the site-measured LAI max (Table 1), uncertainties in the change of LAI with phenology still existed at the site level because the monthly data were extracted from a global LAI map [4].In Equation (2), we simplified the k as 0.5 across a range of biome types.However, this value is dependent on the leaf distribution [94] and the zenith angle, which ranged between 0.3 and 0.6 [95].The best estimates of k would be derived from stands with minimal clumping for each site [94,95].

Link between Parameter Sets and Model Structures
A revision of parameters in GPP simulation could cancel errors caused by model structure, including coupled stomatal conductance and carbon assimilation scheme [7], two-stream radiative transfer, and leaf photosynthesis [2].We evaluated the application of parameter adjustment in the MODIS land algorithm.The results demonstrated that the performances of parameter-optimized LUE and LUE-SW were identical in terms of GPP estimation (Figure 1), but the optimized ε max in the latter model had lower standard deviation than the former for sites within the same PFT (Table 2), even though the parameters categorized according to PFTs are not the best option.After optimizing parameters in simple photosynthesis and transpiration models using measurements from 101 EC flux towers, Groenendijk et al. [81] pointed out that a simple PFT classification could induce the uncertainties in the photosynthesis and water vapor flux estimates, and site-year parameters gave the best predictions.However, this parameter uncertainty could be reduced by adding control factors, such as considering the scale of soil moisture in the one-leaf model (Table 2).Furthermore, we quantified the relationships between annual average soil water factors with optimized ε max in LUE and in LUE-SW, excluding four sites for which ε max was greater than 1.3 gC•MJ −1 (Figure 5).The regression analysis shows that the ε max in LUE decreased linearly with the soil water factor.The reduction was as much as a half in LUE-SW, as evidenced by the slopes reducing from 0.77 to 0.38 and by the constants increasing from 0.32 to 0.73, especially for sites with low annual average soil water factors and large soil water changes.It should be noted that the annual rainfall values of our selected sites were all greater than 200 mm•yr −1 (Table 1).More validations were needed in those water-limited regions with annual precipitation less than 100 mm•yr −1 [8,96].
Moreover, parameter adjustment could not compensate all structural errors, such as the error caused by the canopy upscaling strategy.Although we optimized the parameters using data assimilation, GPP estimated by both LUE and LUE-SW still had large negative biases in summer (Table 3).The biases were reduced by 81.9% using the two-leaf Farquhar model.Schaefer et al. [36] compared 26 models with standard parameters at 39 EC flux tower sites across North America and found that the average GPP bias was −0.87 gC•m −2 •day −1 in summer (Figure 2a), which is between the biases of LUE-SW (−1.21 gC•m −2 •day −1 ) and the Farquhar model (0.17 gC•m −2 •day −1 ) in the same season that we estimated (Table 3).In their study, more than half of the models adopted the one-leaf upscaling strategy.An underestimation in the MODIS GPP product during summer was also reported by Zhang et al. [24], who evaluated different regions across the conterminous U.S. The reason for underestimation is that the one-leaf strategy ignores a large contribution of diffuse PAR to the shaded leaves, and the contribution is more efficiently absorbed by the canopy for photosynthesis than direct PAR.In the strategy of separating the canopy into two parts, the sunlit leaves achieve a high rate of light-saturated photosynthesis and have a low light use efficiency, whereas shaded leaves are only limited by the electron transport related to direct and diffused radiation, which lead to a high light use efficiency [25,73,97,98].Many efforts have been made to improve the structure of the MODIS land algorithm by considering photosynthesis rate saturation or light saturation.Propastin et al. [71] adopted a saturating function for light use efficiency adjustment that allowed for saturation of gross photosynthesis at a high irradiance.This modification improved the performance of the MODIS GPP algorithm for a tropical forest.Separating the canopy into sunlit and shaded leaves, He et al. [13,15] developed a two-leaf light use efficiency model based on the MODIS algorithm, and the model properly described differences in the light use efficiencies of sunlit and shaded leaves.Further improvement should focus on the photosynthesis models' capacities in capturing seasonal change.In this study, all models exhibited large uncertainties in winter and spring, even if the key parameters in the two-leaf Farquhar model have been well optimized (Table 3, Figures 3 and 4).This result is the same as the simulations from 26 models at 39 EC sites [36].Fortunately, many approaches could be referred to quantify seasonal variation of key variables in GPP estimates.A comparison between using the dynamic maximum velocity of carboxylation (Vc max in Equation (A3)) and the constant Vc max in the Farquhar models by Muraoka et al. [99] indicated that Vc max variation had remarkable effects on GPP, and an overestimate of 15% was caused by assuming Vc max to be constant in a cool-temperate deciduous broadleaf forest.Groenendijk et al. [100] upscaled the ecosystem parameters Vc max with LAI for 81 EC sites, but the seasonal variation of Vc max could not be sufficiently explained at the ecosystem scale.By seasonally changing the photosynthetic parameters in a Farquhar-type biochemical model, Zhu et al. [101] successfully reproduced the observed response in net assimilation rates at leaf scale.According to satellite data and the Biome-BGC terrestrial ecosystem model, Ichii et al. [102] suggested that proper setting of the root depths was important to simulate GPP seasonality in tropical forests.On the basis of these studies, a widely accepted concept is needed to improve seasonal change in GPP estimate.

Conclusions
By comparing GPP estimates with parameter adjustment in the MODIS algorithm across nine PFTs at half-hourly, daily, monthly and seasonal scales, our multisite study illustrates as follows: (1) Large bias was observed in the MODIS GPP product, especially in deciduous forests and shrubs and grasslands.Its uncertainties were affected by both input data and the look-up table values of ε max for individual PFTs.It is necessary to optimize the parameters in the look-up table used by the MODIS algorithm, but the optimized parameters should correspond to specific input data for applications, i.e., the optimized parameters cannot be applied to a simulation with changed driver data because errors from parameters and input data can accumulate.
(2) Optimizing the key parameter ε max in the MODIS GPP algorithm can compensate the errors caused by ignoring soil water factor at the site level, but the ε max values would have large uncertainties among sites within the same PFT and among the PFTs, especially for sites with low yearly average soil water factors.This result casts doubt on the accuracy of simulated spatial distribution of GPP yielded by the MODIS algorithm.Moreover, GPP was underestimated by the one-leaf models in summer, regardless of whether the soil water factor was considered, but could be improved by separating the canopy structure into sunlit and shaded parts.This result indicates that improving model structure is a better choice than only adjusting parameters.Photosynthetic dynamics in spring and fall for forests and shrubs and seasonal GPP change for grasslands could not be captured by both one-leaf and two-leaf models.Therefore, there is a need to improve seasonal and phenology variations of key parameters and variables in carbon assimilation calculation to reduce uncertainties in GPP simulation.

Figure 1 .
Figure 1.Performances of the GPP models for the 21 selected tower sites (Table 1).The statistics in the Taylor diagram were derived from the simulated and observed GPP of the representative year for each site: (a) half-hourly values and (b) 8-day average values.An ideal model has a standard deviation ratio (σ norm ) of 1.0 and a correlation coefficient of 1.0 (REF, the reference point).

Figure 2 .
Figure 2. Comparisons of the observed GPP and the MOD17A2: the GPP simulated by the LUE-SW and by the Farquhar model for different plant functional types: (a,b) needleleaf evergreen forests in temperate and boreal zones; (c) needleleaf deciduous forests in boreal zone; (d) broadleaf evergreen forests in temperate zone; (e,f) broadleaf deciduous forests in temperate and boreal zones; (g,h) broadleaf deciduous shrubs in temperate and boreal zones; and (i) grasslands.

Figure 3 .
Figure 3. Boxplots of Taylor skill (S) for daily GPP by models and seasons across (a) evergreen forests, (b) deciduous forests and shrubs, and (c) grasslands.Panels show the interquartile range (box), mean (square), median (solid line), range (whiskers), and outliers (cross).The models and seasons are sorted by the median Taylor skill.

Figure 4 .
Figure 4. Boxplots of bias and root-mean-square error (RMSE) for seasonally averaged diurnal composites simulated by the LUE-SW and the Farquhar model: (a) evergreen forests; (b) deciduous forests and shrubs; and (c) grasslands.

Figure 5 .
Figure 5.The annually averaged soil water factor (β t ) versus the optimized maximum light use efficiency (ε max ) used in the LUE and LUE-SW, excluding four sites with ε max greater than 1.3 gC•MJ −1 (gray symbols).The bar represents ±0.5 standard deviation of β t .

Table 1 .
Descriptions of the study sites.

Table 2 .
Model parameters a derived from the 21 selected tower sites for nine plant functional types.Parameter differences among plant functional types, biome types and climate zones were determined using the one-way analysis of variance.
a The terms ε max and are the maximum light use efficiency and the leaf maximum carboxylation rate at 25 °C constrained by leaf nitrogen, respectively.b Parameters obtained from Zhao and Running et al. [20].c Values in parentheses are standard deviations.

Table 3 .
Bias and root-mean-square error (RMSE) of daily GPP (in gC•m −2 •day −1 ) with respect to individual months and seasons.