Underestimates of Grassland Gross Primary Production in MODIS Standard Products

Abstract

As the biggest carbon flux of terrestrial ecosystems from photosynthesis, gross primary productivity (GPP) is an important indicator in understanding the carbon cycle and biogeochemical process of terrestrial ecosystems. Despite advances in remote sensing-based GPP modeling, spatial and temporal variations of GPP are still uncertain especially under extreme climate conditions such as droughts. As the only official products of global spatially explicit GPP, MOD17A2H (GPP_MOD) has been widely used to assess the variations of carbon uptake of terrestrial ecosystems. However, systematic assessment of its performance has rarely been conducted especially for the grassland ecosystems where inter-annual variability is high. Based on a collection of GPP datasets (GPP_EC) from a global network of eddy covariance towers (FluxNet), we compared GPP_MOD and GPP_EC at all FluxNet grassland sites with more than five years of observations. We evaluated the performance and robustness of GPP_MOD in different grassland biomes (tropical, temperate, and alpine) by using a bootstrapping method for calculating 95% confident intervals (CI) for the linear regression slope, coefficients of determination (R²), and root mean square errors (RMSE). We found that GPP_MOD generally underestimated GPP by about 34% across all biomes despite a significant relationship (R² = 0.66 (CI, 0.63–0.69), RMSE = 2.46 (2.33–2.58) g Cm⁻² day⁻¹) for the three grassland biomes. GPP_MOD had varied performances with R² values of 0.72 (0.68–0.75) (temperate), 0.64 (0.59–0.68) (alpine), and 0.40 (0.27–0.52) (tropical). Thus, GPP_MOD performed better in low GPP situations (e.g., temperate grassland type), which further indicated that GPP_MOD underestimated GPP. The underestimation of GPP could be partly attributed to the biased maximum light use efficiency (ε_max) values of different grassland biomes. The uncertainty of the fraction of absorbed photosynthetically active radiation (FPAR) and the water scalar based on the vapor pressure deficit (VPD) could have other reasons for the underestimation. Therefore, more accurate estimates of GPP for different grassland biomes should consider improvements in ε_max, FPAR, and the VPD scalar. Our results suggest that the community should be cautious when using MODIS GPP products to examine spatial and temporal variations of carbon fluxes.

1. Introduction

Gross primary productivity (GPP), which is also known as the rate of photosynthesis, is the biggest carbon flux of terrestrial ecosystems [1]. This carbon flux plays an important role in the terrestrial carbon cycle. GPP is also the basis for ecosystem services such as food, fuel, and wood products [2]. The ability to accurately track the spatial and temporal variability of GPP is fundamental for understanding the biogeochemical dynamics of terrestrial ecosystems [3,4]. Therefore, it is critical for us to accurately estimate GPP and further understand the trends and variations of global and regional carbon uptake. However, there still exist considerable uncertainties in GPP estimation, which has attracted plenty of attention [5,6,7].

Satellite remote sensing provides an unprecedented and practical opportunity to estimate ecosystem GPP at large scales by using a diagnostic approach. Numerous remote sensing models have been proposed including the Global Production Efficiency Model (GLOPEM) [8], the Carnegie-Ames-Stanford Approach model (CASA) [9], the Vegetation Photosynthesis Model (VPM) [10,11], the light use efficiency model (EC-LUE) [12], and the net photosynthesis model (PSN) [13]. The Moderate Resolution Imaging Spectroradiometer (MODIS) primary production products (MOD17A2) based on the PSN model are open access datasets with high temporal and spatial resolutions that allow for the monitoring of global GPP at the 1-km resolution every eight days [14]. Previous studies using different light use efficiency (LUE) models estimated global GPP ranging from 105 to 177 P g C a⁻¹ [15,16] with a considerable uncertainty. Yuan et al. [17] compared seven LUE models and they found that only two models (EC-LUE and CFlux) showed higher correlations between GPP derived from eddy covariance (GPP_EC) and modeled GPP and performed better in simulating inter-annual variability of GPP than others. Furthermore, the performances of different models are not the same among all vegetation types or biomes [18,19]. As the only official global GPP products (MOD17A2/MOD17A2H), the MODIS GPP (GPP_MOD) products have been widely used in addressing scientific questions on terrestrial carbon uptake and carbon cycle [20,21]. For example, MOD17A2 has been used to understand the drought impact on global carbon uptake of terrestrial ecosystems [22].

However, there is still a great discrepancy about the reliability of GPP_MOD when compared with in-situ data collected in different locations [23]. Due to the unavailability of directly measured GPP data, the validation of MODIS GPP products is still challenging. The eddy covariance technique is an effective approach to evaluate GPP on the landscape scale. A previous study found that the GPP_MOD was reliable and matched well with GPP_EC in the forest biome at Vancouver Island, Canada [24]. However, there are many studies that have highlighted the uncertainties of GPP_MOD data in different ecosystem types. For example, Turner et al. [25] pointed out that the GPP_MOD underestimated GPP at hardwood forest sites and, yet, agreed well with GPP_EC data at boreal forest sites. In terms of savanna ecosystems, Leuning et al. [26] found that the GPP_MOD performed differently for tropical wet/dry savannas. It overestimated GPP during the dry season and arid summer but matched very well during the two wet seasons. Consistent with that in Leuning et al. [26], GPP_MOD underestimated GPP at dry sites located in the Sahel region, according to a study at 12 African sites of savanna ecosystems [27] while the product agreed well with GPP_EC at wet sites. It showed that generally GPP_MOD tend to be underestimated in drier forest ecosystems or savanna ecosystems. However, these studies were mostly based on the previous versions of the MODIS GPP product (Collection 4 or Collection 5), but the latest version (Collection 6) significantly improved the value of the maximum light use efficiency (for example, from 0.68 g C m⁻² d⁻¹ MJ⁻¹ to 0.86 g C m⁻² d⁻¹ MJ⁻¹ for grass) and the product footprint sizes (from 1 km to 500 m) [28]. Whether performance has improved for collection 6 or not in these ecosystems has not been studied on a global scale.

Compared to the forest and savanna biomes, validation of the MODIS GPP in grassland has also been studied extensively. Zhang et al. [29] found that MODIS data underestimated the mean annual GPP by 30% to 50% in the alpine meadow sites in China. Zhu et al. [30] also found that grassland GPP (MOD17A2) was underestimated at the Haibei grassland site during the summer. In the US, Doughty et al. [31] also found that GPP_MOD performed poorly at the improved pasture sites, the native pasture sites, and the winter wheat sites when compared with GPP_EC. However, these analyses were based primarily on limited spatial and temporal coverage of in situ observations. Regional-scale and global-scale evaluations of GPP_MOD for grassland ecosystems has not yet been conducted. With the increasing global distribution of eddy flux towers and the formation of networks (e.g., FLUXNET [32], AmeriFLux [33], EuropeanFLux [34], AsiaFLux [12]), a large scale validation of GPP products has become feasible [24,35]. Recently, the latest global flux dataset of FLUXNET2015 [36] was released and it provided a great opportunity for the validation of GPP_MOD for grassland ecosystems.

In this study, we aimed to conduct a systematic validation of GPP_MOD for grassland ecosystems by using the FLUXNET2015 dataset. The overarching goals of this study are to: (1) examine the performance of GPP_MOD across a network of grassland flux sites, (2) assess the uncertainties of GPP_MOD under different grassland biomes (temperate, tropical, and alpine), and (3) discuss the potential reasons that cause the underestimation of GPP in the MOD17A2H products.

2. Materials and Methods

2.1. GPP_MOD Algorithm

The MOD17A2H GPP product (Collection 6) is a standard global product with a 500-m spatial resolution and eight-day temporal resolution [22]. It is based on the light use efficiency approach, which calculates GPP by using the amount of photosynthetically active radiation (PAR) absorbed by vegetation over a growing season [37]. The algorithm was developed below.

$$GP{P}_{MOD}=\epsilon \times FPAR\times PAR$$

(1)

where ε is the actual light use efficiency, PAR is the photosynthetically active radiation (MJ m⁻²) per unit time, and FPAR is the fraction of PAR absorbed by vegetation canopy.

$$\epsilon ={\epsilon }_{max}\times TMI{N}_{scalar}\times VP{D}_{scalar}$$

(2)

where ε_max (g C m⁻² d⁻¹ MJ⁻¹) is the maximum light use efficiency, which is given in a Biome Parameter Look-up Table (BPLUT) for each land cover type in the PSN model. The TMIN_scalar and VPD_scalar are environmental stress factors of temperature (daily minimum temperature, Tmin, °C) and water (maximum daily vapor pressure deficit, VPD, Pa) and are parameterized according to Equations (3) and (4).

$$TMI{N}_{scalar}=\{\begin{array}{cc}1& TMIN>TMI{N}_{max}\\ \left(TMIN-TMI{N}_{min}\right)/\left(TMI{N}_{max}-TMI{N}_{min}\right)& TMI{N}_{min}\le TMIN\le TMI{N}_{max}\\ 0& TMIN<TMI{N}_{min}\end{array}$$

(3)

$$VP{D}_{scalar}=\{\begin{array}{cc}0& VPD>VP{D}_{max}\\ \left(VP{D}_{max}-VPD\right)/\left(VP{D}_{max}-VP{D}_{min}\right)& VP{D}_{min}\le VPD\le VP{D}_{max}\\ 1& VPD<VP{D}_{min}\end{array}$$

(4)

where TMIN and VPD are the daily minimum temperature (°C) and the average vapor pressure deficit (Pa), TMIN_max and VPD_max are the daily minimum temperature and the average vapor pressure deficit at which ε = ε_max, and TMIN_min and VPD_min are the daily minimum temperature and average vapor pressure deficit at which ε = 0 [28,38]. These parameters were determined, according to the BPLUT. The NASA Global Modeling and Assimilation Office (GMAO) provided input data for the GPP_MOD algorithm including global meteorological simulations of incoming PAR, daily maximum/minimum temperature, and daily maximum/minimum VPD. We downloaded MODIS Collection 6 data for 2000 to 2014 and we extracted GPP_MOD from pixels with the 500-m resolution in which the flux towers were located. GPP_MOD for each site was transformed from monthly totals to monthly means and were converted from kg C m⁻² to g C m⁻².

2.2. CO₂ Eddy Flux and Meteorological Data

We used the FLUXNET2015 datasets (http://www.fluxdata.org) to evaluate the performance of GPP_MOD. The FLUXNET2015 datasets include carbon dioxide (CO₂), water vapor, energy fluxes, shortwave and longwave radiations, photosynthetically active radiation, temperature, and precipitation for each flux tower site. We selected all grassland flux tower sites across the world with more than five years of observations between 2000 and 2014 using the 500-m resolution. The selection included 15 grassland sites with eddy covariance towers. All of these were included in this study (Figure 1 and

Table 1. Basic information of FLUXNET eddy covariance tower sites selected in this study.

Site ID	Site Name	LAT	LON	IGBP Class	Data Range
CH-Cha	Chamau	47.2102	8.4104	Temperate grassland	2005–2014
CH-Fru	Fruebuel grassland	47.1158	8.5378	Temperate grassland	2005–2014
CH-Oe1	Oensingen grassland	47.2858	7.7319	Temperate grassland	2002–2008
DE-Gri	Grillenburg	50.9495	13.5125	Temperate grassland	2004–2014
NL-Hor	Horstermeer	52.2404	5.0713	Temperate grassland	2004–2011
US-SRG	Santa Rita Grasslan	31.7894	−110.8277	Temperate grassland	2008–2014
US-Wkg	Walnut Gulch Kendall grasslands	31.7365	−109.9419	Temperate grassland	2004–2014
AU-Dap	Daly River avanna	−14.0633	131.3181	Tropical grassland	2007–2013
AU-Stp	Sturt Plains	−17.1507	133.3502	Tropical grassland	2008–2014
US-Var	Vaira Ranch-Ione	38.4133	−120.9507	Tropical grassland	2001–2014
AT-Neu	Neustift	47.1167	11.3175	Alpine grassland	2002–2012
CZ-Bk2	Bily Kriz grassland	49.4944	18.5429	Alpine grassland	2006–2012
IT-Mbo	Monte Bondone	46.0147	11.0458	Alpine grassland	2003–2013
RU-Sam	Samoylov	72.3733	126.4978	Alpine grassland	2002–2014
IT-Tor	Torgnon	45.8444	7.5781	Alpine grassland	2008–2014

). The meteorological data including precipitation, air temperature, solar radiation, and VPD were used in this study. We used the gap-filled GPP data derived from the Marginal Distribution Sampling (MDS) method for gap-filling and the night-time partitioning method for the net ecosystem exchange (NEE) partitioning [39]. The daily GPP_EC values were transformed to monthly GPP_EC values (g C m⁻² day⁻¹) for a comparison with the GPP_MOD.

2.3. Evaluation of Model Performance

In this study, GPP_MOD was evaluated against GPP_EC. First, according to the climate characteristics of the grassland sites, we merged the 15 grassland sites into three grassland biomes, i.e., temperate, tropical, and alpine grasslands. Second, we evaluated the model performance by comparing the seasonal dynamics of GPP_MOD and GPP_EC by using the linear regression model (GPP_EC = a × GPP_MOD + b). During the evaluation, modeled data (GPP_MOD) are plotted on the x-axis and observation data (GPP_EC) on the y-axis. This avoids estimating spurious biases in the representation [40]. The coefficient of determination (R²) was used to evaluate the models’ explanatory abilities for variances in GPP. Third, the root mean squared errors (RMSE) were calculated to quantify the agreement between GPP_EC and GPP_MOD during the plant-growing season (GPP > 1 g C m⁻² day⁻¹).

$$RMSE=\sqrt{\frac{1}{n}\times {{\displaystyle \sum }}_{n=1}^n{\left(X_i-Y_i\right)}^2}$$

(5)

where n is the total number of sample points and X_i and Y_i represent the observed and simulated values, respectively.

We used the bootstrapping method to provide confidence intervals necessary to determine whether these indicators (slopes, RMSE, and R²) are different or not and repeated sampling for 4000 times in the bootstrapping analysis.

2.4. Estimation of ${\epsilon }_{max}$

The ${\epsilon }_{max}$ is an important parameter in the GPP_MOD algorithm. The ${\epsilon }_{max}$ for vegetation types can be obtained from an analysis of gross ecosystem exchange of CO₂ and photosynthetic photon flux density (PPFD) at an eddy flux tower site [41,42]. In this study, we used a nonlinear model between GPP and PPFD (at half-hour time step) data to estimate the ${\epsilon }_{max}$ for each grassland site [43,44,45]. The abnormal observations (including zero values or extremely high values) have been excluded before fitting the function. The model is described below [42,43].

$$NEE=\frac{\alpha \times PPFD\times GE{E}_{max}}{\alpha \times PPFD-GE{E}_{max}}-R$$

(6)

where $\alpha$ is the apparent quantum yield and is assumed to be the ${\epsilon }_{max}$. R is the ecosystem respiration, GEE_max is the maximum gross ecosystem exchange (GEE), and NEE is the net ecosystem exchange.

3. Results

3.1. Comparison of GPP_MOD and GPP_EC

Figure 2 illustrates the seasonal variation of GPP_MOD and GPP_EC at all sites. GPP_MOD and GPP_EC exhibit consistency in the magnitudes and seasonal variations at most of the grassland sites. GPP values started near zero in the winter, began to increase in the spring, reached its peak during the summer, and decreased quickly after its peak (Figure 2). The GPP_EC in most sites were underestimated by the GPP_MOD. For example, the GPP_EC were underestimated by 58.48%, 56.24%, and 53.25% in AT-Neu, RU-Sam, and IT-Tor, respectively (Figure 2a), ((l), and (j)). The coefficients of determination (R²) between GPP_MOD and GPP_EC varied from 0.17 (CI, 0.08–0.29) at the RU-Sam site to 0.83 (0.77–0.89) at the DE-Gri site with all being statistically significant at p < 0.05 (Figure 3) and all these analyses were tested with a 95% confidence interval (

Table 2. Statistics of the comparison between monthly GPP_MOD and GPP_EC for all sites by using the linear regression model. The 95% confidence interval using the bootstrapping method is in parenthesis.

Site ID	Slope	R2	RMSE
AT-Neu	1.89 (1.72, 2.07)	0.78 (0.72, 0.83)	4.05 (3.52, 4.55)
AU-Dap	0.92 (0.60, 1.24)	0.28 (0.14, 0.45)	3.31 (2.88, 3.70)
AU-Stp	0.89 (0.70, 1.07)	0.53 (0.34, 0.70)	1.04 (0.83, 1.24)
CH-Cha	1.14 (1.03, 1.25)	0.79 (0.70, 0.86)	3.64 (3.29, 4.00)
CH-Fru	1.27 (1.14, 1.39)	0.78 (0.71, 0.84)	3.19 (2.78, 3.57)
CH-Oe1	1.14 (0.94, 1.34)	0.61 (0.48, 0.74)	3.19 (2.72, 3.66)
CZ-Bk2	1.04 (0.93, 1.15)	0.82 (0.74, 0.88)	1.31 (1.05, 1.54)
DE-Gri	1.15 (1.06, 1.24)	0.83 (0.77, 0.89)	1.90 (1.61, 2.18)
IT-Mbo	1.05 (0.92, 1.18)	0.67 (0.57, 0.76)	2.36 (1.87, 2.82)
IT-Tor	1.30 (1.01, 1.58)	0.50 (0.33, 0.67)	2.81 (2.23, 3.34)
NL-Hor	0.86 (0.71, 1.00)	0.59 (0.48, 0.68)	2.05 (1.72, 2.36)
RU-Sam	0.58 (0.37, 0.79)	0.17 (0.08, 0.29)	0.89 (0.74, 1.03)
US-SRG	0.50 (0.37, 0.63)	0.40 (0.18, 0.59)	0.68 (0.58, 0.77)
US-Var	0.91 (0.76, 1.07)	0.63 (0.46, 0.76)	1.72 (1.47, 1.97)
US-Wkg	1.59 (1.41, 1.77)	0.70 (0.61, 0.78)	0.63 (0.49, 0.75)
ALL	1.22 (1.18, 1.26)	0.66 (0.63, 0.69)	2.46 (2.33, 2.58)

). The linear regression between GPP_MOD and GPP_EC showed that GPP_MOD performed well at DE-Gri (R² = 0.83 (0.77–0.89), RMSE = 1.90 (1.61–2.18) g Cm⁻² day⁻¹), CZ-Bk2 (R² = 0.82 (0.74–0.88), RMSE = 1.31 (1.05–1.54) g Cm⁻² day⁻¹), CH-Cha (R² = 0.79 (0.70–0.86), RMSE = 3.64 (3.29–4.00) g Cm⁻² day⁻¹), CH-Fru (R² = 0.78 (0.71–0.84), RMSE = 3.19 (2.78–3.57) g Cm⁻² day⁻¹), and AT-Neu (R² = 0.78 (0.72–0.83), RMSE = 4.05 (3.52–4.55) g Cm⁻² day⁻¹), but it performed poor at RU-Sam (R² = 0.17 (0.08–0.29), RMSE = 0.89 (0.74–1.03) g Cm⁻² day⁻¹), AU-Dap (R² = 0.28 (0.14–0.45), RMSE = 3.31 (2.88–3.70) g Cm⁻² day⁻¹), US-SRG (R² = 0.40 (0.18–0.59), RMSE = 0.68 (0.58–0.77) g Cm⁻² day⁻¹), and IT-Tor (R² = 0.50 (0.33–0.67), RMSE = 2.81 (2.23–3.34) g Cm⁻² day⁻¹) (Figure 3). Moreover, GPP_MOD had a R² of 0.66 (0.63–0.69) and RMSE of 2.46 (2.33–2.58) g Cm⁻² day⁻¹ for all grassland tower sites (Figure 4).

3.2. Model Performances in Different Grassland Biomes

We also considered the seasonal dynamic of GPP_MOD by analyzing the relationship between GPP_MOD and GPP_EC within each grassland biome. The 15 flux sites were divided into three grassland biomes including the temperate, tropical, and alpine biomes. All analyses were tested with a 95% confidence interval (

Table 3. Statistics regarding the comparison between GPP_EC and GPP_MOD for temperate, tropical, and alpine grassland biomes using a linear regression model. The 95% confidence interval derived by the bootstrapping method is in parenthesis.

Grass Type	Slope	R2	RMSE
Tropical grassland	1.02 (0.83, 1.21)	0.40 (0.27, 0.52)	2.45 (2.13, 2.76)
Temperate grassland	1.25 (1.19, 1.30)	0.72 (0.68, 0.75)	2.40 (2.25, 2.54)
Alpine grassland	1.24 (1.17, 1.32)	0.64 (0.59, 0.68)	2.55 (2.30, 2.79)

). Figure 5 shows the different relationships between GPP_EC and GPP_MOD for each biome type. From Figure 5, we find GPP_MOD had better performance in a temperate grassland (R² = 0.72 (0.68–0.75), RMSE = 2.40 (2.25–2.54) g Cm⁻² day⁻¹) than an alpine grassland (R² = 0.64 (0.59–0.68), RMSE = 2.55 (2.30–2.79) g Cm⁻² day⁻¹) or a tropical grassland (R² = 0.40 (0.27–0.52), RMSE = 2.45 (2.13–2.76) g Cm⁻² day⁻¹) (Figure 5). Furthermore, MODIS GPP products underestimated GPP about 4% for the tropical grassland, 29% for the temperate grassland, and 41% for the alpine grassland, respectively.

4. Discussion

4.1. Underestimation of MODIS GPP in Grasslands and Comparison with Previous Studies

This study showed that GPP_MOD had a promising ability to explain GPP_EC variance (R² = 0.66 (0.63–0.69) and RMSE = 2.46 (2.33–2.58) g Cm⁻² day⁻¹) at grassland sites, but GPP_MOD consistently underestimated GPP_EC (Figure 4). We also found that different grassland biomes had different performances. Specifically, GPP_MOD performed better in the temperate grassland sites than in the alpine and tropical grassland sites (Figure 5).

This study indicated that GPP_MOD substantially underestimated GPP in grasslands, which agreed with previous studies [29,30,31]. For instance, Zhang [29] used GPP_EC measurements to verify GPP_MOD over an alpine meadow on the Tibetan Plateau, which further pointed out that MOD17 GPP products underestimated the mean annual GPP by 30% to 50%. Turner et al. [46] also reported that GPP_MOD underestimated GPP at desert grassland sites. Xiao [47] evaluated GPP_MOD by using GPP data from eddy covariance flux towers and their results showed that the GPP_MOD products underestimated GPP by 45% at the Walnut Gulch Kendall grassland site. Despite the wide use of GPP_MOD to analyze inter-annual variation at regional scales in previous studies [48,49], the results from this study and the previously mentioned studies implied additional caution should be given when using GPP_MOD for inter-annual and intra-annual variation analyses of GPP in grasslands. We found that GPP_MOD underestimates observed GPP and that GPP_MOD explained 66% (63%–69%) of the variance in GPP_EC (R² = 0.66 (0.63–0.69), RMSE = 2.46 (2.33–2.58) g Cm⁻² day⁻¹).

In addition, we found that the GPP_MOD could explain more GPP variance for the temperate and alpine grassland sites than for the tropical grassland sites (Figure 5). Thus, GPP_MOD performed better in lower GPP situations (e.g., drought condition and temperate grassland biome), which further illustrated that the GPP_MOD tended to underestimate GPP.

4.2. Attributing Underestimation in Grassland GPP and Its Implications

In this study, GPP_MOD explained GPP_EC by 72% (68%–75%) at temperate grassland sites, 64% (59%–68%) at alpine grassland sites, and 40% (27%–52%) at tropical grassland sites (Figure 5). Potential reasons for the uncertainty in MODIS GPP products could be due to different factors such as maximum LUE, meteorological data, FPAR, and land cover/vegetation types [50].

These varied performances may be related to the maximum light use efficiency (ε_max) [51]. For the MODIS GPP algorithm, ε_max was determined by using the Biome Properties Look-Up (BPLUT) for given biome types, which was for 0.86 g C MJ⁻¹ grass [38]. The underestimated ε_max values could be the major reason for underestimates of GPP in the grasslands [51]. Based on the nonlinear model (Equation (6)), we estimated the ε_max for each grassland site in this study (

Table 4. The comparison of the ${\epsilon }_{max}$ between the estimated values based on the nonlinear model (${\epsilon }_{max}$-EST) and the values in the biome properties look-up table (${\epsilon }_{max}$-BPLUT) used in the MOD17 GPP algorithm. “RMSE after” and “RMSE before” referred to the RMSE between the GPP_EC and GPP_MOD with ${\epsilon }_{max}$-EST and ${\epsilon }_{max}$-BPLUT, respectively.

Site ID	${\mathit{\epsilon }}_{\mathit{m}\mathit{a}\mathit{x}}$-EST	${\mathit{\epsilon }}_{\mathit{m}\mathit{a}\mathit{x}}$-BPLUT	RMSE After	RMSE Before	${\mathit{\epsilon }}_{\mathit{m}\mathit{a}\mathit{x}}$ bias
	g C/MJ	g C/MJ	g Cm−2 day−1	g Cm−2 day−1	g C/MJ
AT-Neu	1.71	0.86	2.30	4.05	0.85
AU-Dap	0.87	0.86	3.30	3.31	0.01
AU-Stp	0.73	0.86	0.98	1.04	−0.13
CH_Cha	0.97	0.86	3.19	3.64	0.11
CH-Fru	1.28	0.86	2.24	3.19	0.42
CH-Oe1	1.13	0.86	2.64	3.19	0.27
CZ-Bk2	2.57	0.86	7.42	1.31	1.71
DE-Gri	1.09	0.86	1.53	1.90	0.23
IT-Mbo	0.98	0.86	2.24	2.36	0.12
IT-Tor	1.13	0.86	2.55	2.81	0.27
NL-Hor	0.93	0.86	2.04	2.05	0.07
RU-Sam	0.93	0.86	0.89	0.89	0.07
US-SRG	0.77	0.86	0.57	0.68	−0.09
US-Var	0.76	0.86	1.52	1.72	−0.10
US-Wkg	1.07	0.86	0.60	0.63	0.21

). In most sites (except for AU-Stp, US-SRG, and US-Var), the estimated values of ε_max were larger than that used in the MOD17 GPP algorithm, which can partly explain the underestimation of the MODIS GPP in these grassland sites. Moreover, when the ε_max used in the MOD17 GPP algorithm (${\epsilon }_{max}$-BPLUT) were simply replaced with the estimated ε_max (${\epsilon }_{max}$-EST), the RMSEs were smaller than before (except for the site of CZ-Bk2) (Figure 6). This demonstrated that the ε_max have an important influence on the result of the GPP estimation based on the MOD17 GPP algorithm. Previous studies have also suggested different ε_max for different grassland biomes. For example, 1.31 g C MJ⁻¹ for alpine grassland and 1.21 g C MJ⁻¹ for tropical grassland [52,53]. Sjöström et al. [27] increased ε_max for grasslands from 0.86 g C MJ⁻¹ to 2.01 g C MJ⁻¹, which improved the ability for GPP_MOD to explain the variance of GPP_EC from 25% to 74%.

There are also potential uncertainties from FPAR, which could affect the accuracy of GPP_MOD. The GPP_MOD algorithm directly used the MODIS FPAR products (MOD15A2H), which was derived from a radiative transfer model. The MODIS FPAR usually overestimated the “green” FPAR with a positive offset for barren land and was saturated at high GPP values. This behavior further affected ε and GPP estimation [54]. In this study, GPP_MOD underestimated GPP_EC for grassland ecosystems, which may be caused by the uncertainty of FPAR and the methods used to simulate FPAR. For example, Liu et al. [55] used a three-dimensional formulation of the radiative transfer process in the canopy and assessed FPAR_canopy performances of scaled EVI (FPAR_chl1), NDVI, scaled NDVI (FPAR_chl2), and EVI. The results showed that the FPAR_canopy of scaled EVI (FPAR_chl) improved the accuracy of GPP_MOD for grasslands.

Another previous study [55] compared MODIS-FPAR and GLASS-FPAR (derived from the product of the Global Land Surface Satellite (GLASS)). MODIS-FPAR performed well with a higher R² and a lower bias for mixed forests and cropland. However, GLASS-FPAR performed better in grasslands with a higher R² and a lower bias [55]. Furthermore, Kucharik et al. [33] and Wu et al. [51] pointed out that the adaptive Savitzky-Golay filtered FPAR data could have better performance compared with the FPAR data based on the linear interpolation approach used in the MODIS algorithm. The adaptive SG method could maintain some intrinsic seasonal variability and keep reliable values of FPAR. Further uncertainties also arise due to the canopy photosynthesis saturation effects on the FPAR [56]. In this study, GPP_MOD products underestimated GPP in tropical grassland sites because of the saturation of canopy photosynthesis. A previous study showed that the MODIS light use efficiency approach does not account for saturation of canopy photosynthesis under clear sky conditions in a tropical forest site [50]. Saturation can result in the underestimation of GPP with low values of FPAR but overestimation with high FPAR values [50]. The tropical grassland received solar radiation for extended periods and it was easy for canopy photosynthesis saturation to occur [57] especially in periods of drought. Previous studies also reported that saturation of canopy photosynthesis might occur on clear-sky days at the hourly and daily timescales [25,58,59].

The species composition of grassland ecosystems could also cause uncertainty in GPP estimates. Grasslands are dynamic ecosystems that often have diverse mixtures of C3 and C4 species, which have different light use efficiencies. Furthermore, the ratio of C3 and C4 species within a grassland can have high intra-annual and inter-annual variability [31]. Such community composition data is not available for all the FLUXNET sites. The effects of the composition of C3 and C4 grasses on GPP simulation is an interesting question, but it is beyond the scope of this specific study. However, a recent GPP product considering the C3 and C4 photosynthesis pathways showed improved performances in grassland ecosystems [60]. The water stress scalar could also have uncertainty when using VPD as a proxy of water stress in the PSN model. To decrease uncertainties, we suggest using long-term meteorological observation data or remote sensing-based water stress indicators [61] to quantify the water stress scalar especially in drought conditions. Furthermore, fundamental improvements could be made to GPP_MOD products. For instance, the Soil Moisture Active Passive (SMAP) mission has two to three day temporal fidelity and an enhanced (≤9 km) spatial resolution that provide accurate global mapping of the freeze–thaw state and the surface soil moisture [62].

4.3. Model Performances under Drought and Non-Drought Conditions

GPP_MOD could have different performances under drought conditions. We used the Standardized Precipitation Index (SPI) to identify the periods of drought. The SPI is designed to be a spatially invariant indicator of drought, which recognizes the importance of time scale in the analysis of water use [63]. In this study, the precipitation from the FluxNet sites was used to calculate SPI for each month to determine the status of drought [64]. If SPI was less than −0.5, it was defined as drought [65]. We analyzed the relationships between monthly GPP_MOD and GPP_EC in both non-drought and drought conditions. The results showed that GPP_MOD had higher R² (0.72 (0.66–0.79)) and lower RMSE (2.33 (2.05–2.60) g Cm⁻² day⁻¹) in drought conditions than in non-drought conditions with lower R² (0.64 (0.61–0.68)) and higher RMSE (2.48 (2.34–2.62) g Cm⁻² day⁻¹), which indicates that GPP_MOD had better explanatory capabilities for GPP_EC variances under drought conditions than under non-drought conditions (

Table A1. The relationship between GPP_EC and GPP_MOD for the temperate, tropical, and alpine grassland biomes in both non-drought and drought conditions. The 95% confidence interval using the bootstrapping method is in parenthesis.

	Type	Slope	R2	RMSE
Drought conditions	Tropical grassland	1.05 (0.26, 1.85)	0.37 (0.01, 0.86)	3.55 (2.39, 4.49)
	Temperate grassland	1.22 (1.12, 1.32)	0.79 (0.73, 0.85)	2.14(1.86, 2.42)
	Alpine grassland	1.20 (1.04, 1.36)	0.68 (0.57, 0.78)	2.36 (1.83, 2.84)
	ALL	1.19(1.10, 1.28)	0.72 (0.66, 0.79)	2.33 (2.05, 2.60)
Non-drought conditions	Tropical grassland	1.01 (0.81, 1.20)	0.40 (0.28, 0.54)	2.31 (1.99, 2.62)
	Temperate grassland	1.26 (1.20, 1.32)	0.70 (0.66, 0.74)	2.45 (2.28, 2.61)
	Alpine grassland	1.26 (1.17, 1.35)	0.63 (0.58, 0.68)	2.59 (2.32, 2.87)
	ALL	1.23 (1.18, 1.28)	0.65 (0.61, 0.68)	2.48 (2.34, 2.62)

).

Although several studies have reported that the standard MODIS GPP products did not accurately estimate carbon uptake during drought conditions [24,66,67]. The role of VPD in determining GPP has been examined [68]. For example, Turner et al. [46] pointed out that a high VPD scalar could reduce GPP values during a dry period. Another study also found that VPD had a weak positive correlation with GPP. Furthermore, VPD sometimes failed to capture drought events [69]. In this study, however, it is notable that GPP_MOD had a better performance in drought conditions than in non-drought conditions for grasslands.

We analyzed the relationships between the VPD-based water scalar and GPP_MOD and found a higher R² between GPP and VPD-based water scaler in drought conditions (Figure 7) even though they had weak correlations in both drought and non-drought conditions (26% vs. 8%, Figure 7b,c). The MOD17A2H GPP products of grassland have better performance in drought conditions, which could be attributed to the W_scalar. However, there could be some uncertainty in the drought identification since SPI calculation was based on a short-term precipitation record in some sites. More studies are still needed in the future for a thorough analyses. The weak correlation between W_scalar and GPP suggest that the water scalar could be further improved in the grasslands.

5. Conclusions

The GPP of grassland ecosystems plays a vital role in carbon sequestration, food production, and biodiversity [70]. In this study, we investigated and evaluated the performance and robustness of GPP_MOD at grassland sites across the globe. We found that: (1) GPP_MOD validation at 15 global eddy covariance sites suggested a high explanation capability of GPP variances (66%) but GPP_MOD generally underestimated GPP about 34% for grassland ecosystems. (2) GPP_MOD performed better for the temperate grassland (R² = 0.72 (0.68–0.75), RMSE = 2.40 (2.25–2.54) g Cm⁻² day⁻¹) than the alpine (R² = 0.64 (0.59–0.68), RMSE = 2.55 (2.30–2.79) g Cm⁻² day⁻¹) and the tropical grasslands (R² = 0.40 (0.27–0.52), RMSE = 2.45 (2.13–2.76) g Cm⁻² day⁻¹). The inconsistent underestimates of GPP_MOD for the three grassland biomes may be rooted in the MODIS GPP algorithm, which underestimated ε_max for the different grassland biomes. The varied performances of the GPP standard products in different grassland types implied that the parameterization of the MOD17A2H GPP products (e.g., ε_max, FPAR, and water scalar) could be improved to better capture changes in GPP.