Exploring Simple Algorithms for Estimating Gross Primary Production in Forested Areas from Satellite Data

Abstract

Algorithms that use remotely-sensed vegetation indices to estimate gross primary production (GPP), a key component of the global carbon cycle, have gained a lot of popularity in the past decade. Yet despite the amount of research on the topic, the most appropriate approach is still under debate. As an attempt to address this question, we compared the performance of different vegetation indices from the Moderate Resolution Imaging Spectroradiometer (MODIS) in capturing the seasonal and the annual variability of GPP estimates from an optimal network of 21 FLUXNET forest towers sites. The tested indices include the Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), Leaf Area Index (LAI), and Fraction of Photosynthetically Active Radiation absorbed by plant canopies (FPAR). Our results indicated that single vegetation indices captured 50–80% of the variability of tower-estimated GPP, but no one index performed universally well in all situations. In particular, EVI outperformed the other MODIS products in tracking seasonal variations in tower-estimated GPP, but annual mean MODIS LAI was the best estimator of the spatial distribution of annual flux-tower GPP (GPP = 615 × LAI − 376, where GPP is in g C/m²/year). This simple algorithm rehabilitated earlier approaches linking ground measurements of LAI to flux-tower estimates of GPP and produced annual GPP estimates comparable to the MODIS 17 GPP product. As such, remote sensing-based estimates of GPP continue to offer a useful alternative to estimates from biophysical models, and the choice of the most appropriate approach depends on whether the estimates are required at annual or sub-annual temporal resolution.

1. Introduction

Understanding the mechanisms of global climate change and improving our predictions of possible future climates requires accurate quantifications of carbon uptake by terrestrial vegetation. Gross Primary Production (GPP), the amount of carbon assimilated by vegetation through photosynthesis over time, fundamentally influences the carbon cycle and human habitability of our planet [1]. While the mechanisms regulating GPP are well understood and modeled at the leaf and canopy level, regional and global GPP models still need to address large uncertainties in the estimates [2].

Historically, scientists have tried to explain the spatial and temporal patterns of GPP by exploiting its relation with vegetation biophysical metrics such as Leaf Area Index (LAI) [3], foliar standing crop [4], the Fraction of Photosynthetically Active Radiation absorbed by plant canopies (FPAR) [5], or with vegetation indices like the integrated Normalized Differential Vegetation Index (NDVI) [6]. Following the recognition that drought conditions suppress summer photosynthesis [7], process-oriented models were developed to estimate GPP with a combination of meteorological data and satellite-observed biophysical parameters, such as FPAR [8] and LAI [9]. These process-oriented models, while allowing simulations of ecosystems under future climate scenarios, present large uncertainty when used for current estimations of GPP because the required input parameters are not available with adequate accuracy at the global scale. The possibility of estimating GPP based on meteorological data only (e.g., Kato et al. [10]), combinations of satellite-observed biophysical parameters (e.g., Wu et al. [11]; Xiao et al. [12]), or combinations of meteorological data and simple biophysical processes (e.g., temperature × (precipitation − evapotranspiration) [13]), has also been explored.

For global-scale studies of GPP, satellite data have the advantage of providing continuous, repetitive, and consistent observations of vegetation dynamics. The Advanced Very High Resolution Radiometer (AVHRR, available since 1981 and originally designed for weather monitoring) has been widely used for long-term analysis of global vegetation dynamics [14,15] but suffers from orbital drift, water vapor contamination and lack of onboard calibration [16]. In 1999 NASA launched the Moderate Resolution Imaging Spectroradiometer (MODIS) onboard the TERRA platform with the specific objective of monitoring Earth’s surface and atmosphere environment. MODIS data have now been used to generate a decade’s worth of standard land data products including vegetation indices, biophysical properties of vegetation, and carbon fluxes [17].

Concurrent with MODIS, participants in the global FLUXNET network of micrometeorological towers with eddy covariance instrumentation have recorded land-atmosphere exchanges of carbon dioxide, water vapor, and energy at hundreds of diverse sites [18,19]. Although uncertainties exist in the FLUXNET observations, the network has provided unprecedented estimates of GPP at high temporal resolution for almost a decade.

The development of these two data streams—MODIS and FLUXNET—has led to extensive tests of the ability of remote sensing to approximate GPP, especially at local scales. Among the MODIS products, the Enhanced Vegetation Index (EVI)—an optimized vegetation index developed from NDVI [20]—has shown a strong seasonal correlation with GPP [21–24]. Sims et al. [23] found a high correlation between GPP and EVI in deciduous forests but not in evergreen forests. Huete et al. [24] extended the relationship between GPP and EVI to the monsoon forests of Southeast Asia. Xiao et al. [21] indicated that EVI could track FPAR by measuring photosynthetically active vegetation while excluding the senescent parts of vegetation. Some biophysical models have already tried to use EVI in place of NDVI to estimate vegetation carbon production [21,25–27].

Although past research seems to have established the seasonal correlation of GPP with EVI at individual flux tower sites, questions still remain whether or not the relationship is consistent across sites, whether it is superior to the other existing satellite metrics of vegetation abundance and structure (NDVI, LAI, and FPAR), and whether it lends itself to estimating annual GPP. Although there are now a number of process-based models for GPP estimation, the use of a single satellite-derived metric still has great appeal for ease of application, assessing the effects of uncertainties caused by poor quality input data, and simplifying the exploration of the underlying mechanisms of vegetation responses to climate change and variability. In this study we analyzed and compared how well four MODIS vegetation products relate to short-term and annual GPP estimated at a network of ideal-quality FLUXNET sites.

2. Data and Methods

Our overall approach was to compare a suite of remotely-sensed vegetation metrics to GPP estimated at a quality-screened subset of the FLUXNET sites. In the following sections, we describe: (1) the selection and processing of MODIS data stored as subsets for the FLUXNET sites, (2) the conceptual basis of the MODIS products, (3) the processing and selection of FLUXNET data, and (4) our analytical approach.

2.1. MODIS Land Products Subsets

Collection 5 MODIS Land Products Subsets for FLUXNET sites are readily available from the Distributed Active Archive Center (DAAC) at Oak Ridge National Laboratory [28]. These subsets represent MODIS product cutouts of 7 km by 7 km around each FLUXNET tower, greatly simplifying the task of MODIS data retrieval for the FLUXNET sites.

We retrieved Collection 5 MOD13Q1 (NDVI and EVI) and MOD15A2 (LAI and FPAR), and Collection 5.1 MOD17A2 (GPP) (hereafter referred to as MOD13, MOD15, and MOD17, respectively) for the period 2001 to 2008 (

Table 1. Description of Moderate Resolution Imaging Spectroradiometer (MODIS) product subsets.

	Products	Spatial Resolution	Compositing Period	Citation
MOD13Q1	NDVI and EVI	250 m	16 day	Huete et al. [20]
MOD15A2	LAI and FPAR	1 km	8 day	Myneni et al. [30]
MOD17A2	GPP/NPP	1 km	8 day	Zhao et al. [31]

). We used MOD17—a process-oriented model driven by FPAR and meteorological data—as a comparison with our simpler, single-metric approaches of estimating GPP. Unfortunately, MOD13 is only archived at 250 m, while the other products have a spatial resolution of 1 km. To make the subsets comparable, we averaged the sixteen MOD13 pixels over the 1 km area surrounding each tower. For each site, we then selected the central pixel from the 7 km by 7 km subset to represent the FLUXNET tower footprint and assembled the 2001–2008 MODIS products time series. The choice of working with the central 1km pixel is justified by the fact that the geolocation accuracy of MODIS is about 50 m (1 σ) at nadir [29] and because for most FLUXNET sites the representative vegetation cover extends for less than 1 km.

We only used data from years in which MODIS data were at least 80% complete (i.e., years for which more than 37 and 19 valid data per year were available for 8-day and 16-day composite periods, respectively). To compare across MODIS products, we did not use the quality flags (except for the snow flag), because the categories and the procedures defining the MODIS quality assurance fields vary for each MODIS product. Instead of using the quality flags, we applied the maximum value composite (MVC) method [32], and set all the non-processed data values to the corresponding maximum. If all the data in the MVC window were not processed, we labeled the pixel value as missing data. To assess the stability of the results with respect to the composite period, we processed the MODIS products to 8-day, 16-day, and 32-day composites.

2.1.1. MOD13 NDVI and EVI

MOD13 is the 250 m, 16-day composite NDVI and EVI product, derived from MOD09 reflectances. Even though vegetation indices are a unitless metric of vegetation greenness, they are widely used as monitoring tools because of the simplicity of their formulas. NDVI is defined as:

$$\mathit{NDVI}=\frac{{\rho }_{\mathit{NIR}}-{\rho }_{\mathit{red}}}{{\rho }_{\mathit{NIR}}+{\rho }_{\mathit{red}}}$$

(1)

where ρ_NIR is near infrared reflectance (841–876 nm) and ρ_red is red reflectance (620–670 nm) .

Among all vegetation indices, NDVI is the most widely used for monitoring and modeling vegetation dynamics [33,34]. However, NDVI-based monitoring of vegetation dynamics suffers from several drawbacks: it saturates over dense vegetation, is frequently contaminated by clouds and aerosols, and is sensitive to soil background reflectance. EVI [20] was introduced as a partial solution to these limitations of NDVI, and is calculated as:

$$\mathit{EVI}=G\frac{{\rho }_{\mathit{NIR}}-{\rho }_{\mathit{red}}}{{\rho }_{\mathit{NIR}}+C_1{\rho }_{\mathit{red}}-C_2{\rho }_{\mathit{blue}}+L}$$

(2)

EVI empirically corrects for aerosols with the reflectance in the blue band (ρ_blue, 459 to 479 nm) and adds L (1) as a background brightness correction factor. G (2.5) is the gain factor and C₁ (6) and C₂ (7.5) are band-specific atmospheric resistance correction coefficients.

We used the MOD13 quality assurance fields to identify the presence of snow, in which case we set the values of both NDVI and EVI to zero.

2.1.2. MOD15 LAI and FPAR

MOD15 is the 1 km, 8-day composite LAI and FPAR product [35]. The main algorithm for retrieving LAI and FPAR is based on a lookup table (LUT) inversion approach. The biome-specific LUTs, generated using a three-dimensional radiative transfer model [30,35], consists of spectral surface reflectances as a function of LAI, soil reflectances and view geometry. The inversion results in minimizing the difference between the observed surface reflectances and the simulated LUT entries to generate an average LAI value from a solution set [30]. If the main algorithm fails, a back-up algorithm based on the relationship between NDVI and LAI/FPAR is used in its place (in this study, the main algorithm had been applied to 72.8% of the LAI/FPAR data). Summarizing LAI validation studies, Yang et al. [35] showed that the MOD15 LAI product had an R² of 0.87 when compared to measured LAI from 29 field campaigns, with an average overestimation of 12% (RMSE = 0.66). In a separate validation effort over the tropics, Aragão et al. [36] reported an overestimation by a factor of 1.18 at Tapajos, Brazil. The new Collection 5 (C5) MOD15 data, used in this study, was an improvement over the earlier collections, especially for the forest biomes [35,37,38]. A recent study by De Kauwe et al. [37] shows that the C5 MOD15 data has higher spatial consistency (R² = 0.42) than Collection 4 (R² = 0.28), when compared to upscaled in situ measurements of mixed coniferous forests. We screened the C5 MOD15 product by assigning the LAI and FPAR of each 8-day composite to zero when the quality assurance fields indicated that snow was present over a pixel. We calculated annual mean LAI and FPAR as for NDVI and EVI.

2.1.3. MOD17 GPP

MOD17 is the 1 km 8-day average GPP and is based on a light use efficiency (LUE) model [31]:

$$\mathit{GPP}={\epsilon }_{\text{max}}*\mathit{PAR}*\mathit{FPAR}*f(\mathit{VPD})*f(\mathit{Tmin})$$

(3)

where ε_max is the biome-specific maximum light use efficiency (g C/MJ), PAR is photosynthetically active radiation (MJ), and potential GPP is linearly scaled between 1 (no limitation) to 0 (no GPP) using biome-specific functions for Vapor Pressure Deficit (f(VPD)) and minimum temperature (f(Tmin)) [39]. FPAR is from MOD15 and includes filling of cloud-contaminated pixels, while meteorological inputs are from the Goddard Earth Observing System Model (GEOS) of the NASA Data Assimilation Office (DAO) [40]. In this study, MOD17 was used as a reference dataset since it had been well validated with flux tower data [31,39].

2.2. FLUXNET Dataset

We used standardized Level 4 FLUXNET data provided by AmeriFlux and CarboEurope networks. The Level 4 FLUXNET data include meteorological and eddy covariance flux data and estimated GPP (note that we use “GPP” interchangeably with “Gross Ecosystem Production”, though a small difference in photorespiration exists between the two). Level 4 data contain GPP derived from estimated Net Ecosystem Exchange (NEE) by different methods (two different storage term calculations and two different gap-filling methods). The storage term required for the estimation of NEE was either provided by each tower principal investigator (original) or estimated by a standard method assuming that CO₂ concentration measured at eddy flux measurement height is constant within the canopy (standard). The gaps were filled either by the Marginal Distribution Sampling method [41] or by the Artificial Neural Network method [42]. We used the original Marginal Distribution Sampling GPP because the original storage term is more reliable than the standard method at tall tower sites [43] and the Artificial Neural Network approach may fill gaps that are too long. The u* filtering was applied to NEE following Reichstein et al. [41] and the time-series were finally gap-filled at the daily, weekly, and monthly time scale.

The partitioning of the measured NEE fluxes into GPP and Total Ecosystem Respiration (TER) was based on a temperature-dependent exponential model of TER estimated based on nighttime NEE data. GPP was finally obtained as the sum of NEE and TER according to Reichstein et al. [41]. Results of this partitioning method have been recently validated at a large number of FLUXNET sites against alternative estimates based on light response curves [44].

Scaling from point measurements at the flux tower to the 1 km MODIS resolution is not trivial because the heterogeneity of the flux tower footprint can lead to a time-variant measurement bias. Ideally, we should use flux towers whose eddy covariance footprints are homogeneous over the MODIS pixel, but this condition is not often met. In some studies (e.g., Turner et al. [45]), finer resolution satellite data have been used to bridge between flux tower data and coarse resolution data to account for the heterogeneity of the tower sites. Here, to avoid scaling issues, we selected flux tower sites where the land cover extended homogenously roughly over a 500 m radius circle around the tower (assessed by visual interpretation of Google Earth imagery). Thus, we only selected AmeriFlux and CarboEurope forest flux tower sites (

Table 2. Site description of flux towers used in this study. Forest types are: deciduous broad leaf forest (DB), evergreen needleleaf forest (EN), and evergreen broad leaf forest (EB). Note that the years given in the table refer to the available data before applying the criteria described in Sections 2.1 and 2.2, and thus the actual number of years used in the analysis may be smaller.

	Site Name	Abbreviation	Country/State	Latitude (°N)	Longitude (°E)	Forest Type	Year	Citation
Non-tropics	Hainich	DE-Hai	Germany	51.079	10.452	DB	2001–2007	Kutsch et al. [47]
	Tharandt	DE-Tha	Germany	50.964	13.567	EN	2001–2003	Grünwald et al. [48]
	Wetzstein	DE-Wet	Germany	50.454	11.458	EN	2002–2008	Rebmann et al. [49]
	San Rossore	IT-SRo	Italy	43.730	10.287	EN	2001–2004	Chiesi et al. [50]
	Loobos	NL-Loo	Netherlands	52.168	5.744	EN	2001–2008	Dolman et al. [51]
	UCI-1964 burn site	CA-NS3	Canada	55.912	−98.382	EN	2001–2005	McMillan et al. [52]
	UCI-1989 burn site	CA-NS6	Canada	55.917	−98.964	EN	2001–2005	McMillan et al. [52]
	Howland Forest (west tower)	US-Ho2	Maine, US	45.209	−68.747	EN	2001–2004	Davidson et al. [53]
	Morgan Monroe State Forest	US-MMS	Indiana, US	39.323	−86.413	DB	2001–2006	Schmid et al. [54]
	Missouri Ozark	US-MOz	Missouri, US	38.744	−92.200	DB	2004–2006	Gu et al. [55]
	Metolius Intermediate Pine	US-Me2	Oregon, US	44.452	−121.557	EN	2002–2007	Irvine et al. [56]
	Metolius First Young Pine	US-Me5	Oregon, US	44.437	−121.567	EN	2001–2002	Irvine et al. [56]
	North Carolina Loblolly Pine	US-NC2	North Carolina, US	35.803	−76.668	EN	2005–2006	Domec et al. [57]
	Niwot Ridge Forest	US-NR1	Colorado, US	40.033	−105.546	EN	2001–2007	Monson et al. [58]
	Mize	US-SP2	Florida, US	29.765	−82.245	EN	2001–2004	Clark et al. [59]
	Univ. of Mich. Biological Station	US-UMB	Michigan, US	45.560	−84.714	DB	2001–2006	Schmid et al. [60]
Tropics	Santarem Km83 Logged Forest	BR-Sa3	Brazil	−3.018	−54.971	EB	2001–2003	Saleska et al. [61]
	Bukit Soeharto	ID-Buk	Indonesia	-0.833	117.050	EB	2002	Huete et al. [24]
	Palangkaraya	ID-Pal	Indonesia	−2.345	114.036	EB	2002–2003	Hirano et al. [62]
	Mae Klong	TH-Mae	Thailand	14.575	98.858	DB	2003–2004	Huete et al. [24]
	Sakaerat	TH-Sak	Thailand	14.493	101.922	EB	2002–2003	Huete et al. [24]

and Figure 1) that satisfied the following requirements:

• The tower was surrounded by homogeneous land cover for a radius of at least 500 m;
• Temporally overlapping coverage with the MODIS subset data existed;
• The gap-filled ratio was less than 20% for all 12 months in a year.

The only tropical forest AmeriFlux tower meeting the above requirements was in Tapajos, Brazil. To increase the sample size for tropical forests and improve the robustness of our analysis, we added four AsiaFlux tropical sites from Southeast Asia: Bukit Soehart, Palangkaraya, Mae Klong, and Sakaerat. Only annual GPP data were published for these sites in Hirata et al. [46], and were derived from nighttime respiration estimated via exponential respiration model as in Reichstein et al. [41]. Following our full screening, we assembled a total of 21 FLUXNET sites with corresponding MODIS subset data for the years 2001 to 2008. For brevity, we refer to the sites with the codes defined in Table 2.

To assess the degree to which the selected FLUXNET sites were representative of the climatic range of the global forested regions, we first generated a 1982 to 2005 average global climatology of temperature and precipitation. We used the Climate Research Unit (CRU) TS3.1 dataset [63] and the MCD12C1 V005 [64] IGBP classification in 2001 (aggregated to a 0.5 degree resolution map by specifying the most frequent land cover in each pixel). We calculated the FLUXNET sites’ annual mean temperature and precipitation during the flux measurement period and compared the two datasets in the temperature-precipitation domain following Wang et al. [65]. The selected sites ranged from −5 °C to 25 °C in mean annual temperature and from 200 mm/year to 2,300 mm/year in mean annual total precipitation, largely covering the climatic variability of global forested areas (Figure 2).

2.3. Reduced Major Axis Regression Analyses

We examined the relationship between flux tower GPP and MODIS products with reduced major axis regression analysis, which minimizes the sum of the product of deviations from the model along both the x and y axes [66], both at the short-term and at the annual time scale. The correlation was evaluated with the coefficient of determination R². The p-value testing the null hypothesis that the slope of the regression line is equal to zero was calculated using two-tailed t-tests.

2.4. Short-Term Analysis

The first question we asked is which MODIS product(s) best explains the short-term (i.e., intra-annual) variations in flux tower GPP. We averaged daily flux tower GPP estimates over each composite periods (8-day, 16-day, and 32-day), and applied MVC to each MODIS product to the corresponding composite periods. We then calculated the R² between the flux tower average GPP and corresponding satellite-derived composite for each available year from 2001 to 2008, and then averaged the R² of all years. R² between flux tower GPP and satellite composites were also calculated after stratifying by forest type class (i.e., deciduous, evergreen, non-tropical, and tropical). In the same way as done with average R², we also calculated average root mean square error (RMSE) using the reduced major axis regression.

2.5. Annual Analysis

The second question we asked is which MODIS product(s) best explains the annual variations in flux tower GPP. We calculated the annual mean of each MODIS product as the sum of the available values divided by the total number of available values. We then analyzed the regressions between flux tower annual total GPP and the annual mean of each MODIS product for the data organized into the following groups: all sites, non-tropical sites only, and tropical sites only. Similarly to the short-term analysis, we varied the composite period from 8, 16 to 32-days. We did not composite MOD17 because this would have caused to neglect its biophysically-induced reductions (e.g., reductions induced by drought periods rather than by cloud or aerosol contamination) and therefore lead to an overestimation of annual GPP.

Although annually integrated vegetation indices are sometimes used for this kind of analyses (e.g., Goward et al. [6]), we preferred to use annual mean vegetation indices. Satellite data are plagued by frequent missing data due to cloud contamination, high aerosol, etc., causing the annually integrated indices to be affected by the number of valid data. Also, the use of mean values facilitates the comparisons across deciduous and evergreen forest types.

3. Results

3.1. Short-Term Correlation between MODIS Products and Flux Tower GPP

Short-term GPP was more highly correlated with EVI than with NDVI, LAI, or FPAR (

Table 3. Average R² between each MODIS product and flux tower GPP for each composite period. The values in parenthesis are root mean square error (RMSE; gC/m²/day) calculated by reduced major axis regression analysis. Average R² and RMSE are the average of all the R² and RMSE calculated from each year of each site, respectively. n is number of sites used for the calculations.

		Composite Period	All (n = 21)	Deciduous (n = 5)	Evergreen (n = 16)	Non-Tropical (n = 16)	Tropical (n = 5)
MOD13	NDVI	16-day	0.44 (2.08)	0.64 (2.36)	0.38 (1.98)	0.51 (2.14)	0.13 (1.78)
		32-day	0.47 (1.94)	0.65 (2.31)	0.41 (1.80)	0.54 (2.00)	0.16 (1.68)
	EVI	16-day	0.55 (1.70)	0.78 (1.73)	0.47 (1.69)	0.64 (1.70)	0.11 (1.67)
		32-day	0.54 (1.67)	0.77 (1.80)	0.46 (1.62)	0.62 (1.72)	0.20 (1.42)
MOD15	LAI	8-day	0.38 (2.21)	0.60 (2.56)	0.30 (2.08)	0.44 (2.31)	0.13 (1.82)
		16-day	0.43 (1.86)	0.69 (2.05)	0.34 (2.02)	0.50 (2.09)	0.13 (1.75)
		32-day	0.51 (1.70)	0.69 (2.09)	0.38 (1.77)	0.56 (1.89)	0.10 (1.72)
	FPAR	8-day	0.31 (2.42)	0.54 (2.81)	0.23 (2.28)	0.35 (2.57)	0.14 (1.81)
		16-day	0.36 (2.26)	0.62 (2.45)	0.27 (2.20)	0.41 (2.37)	0.12 (1.77)
		32-day	0.40 (1.90)	0.62 (2.38)	0.33 (1.91)	0.47 (2.19)	0.12 (1.63)
MOD17	GPP	8-day	0.68 (1.31)	0.72 (1.95)	0.67 (1.11)	0.82 (1.17)	0.12 (1.91)

), in agreement with other studies that have shown high correlation between EVI and GPP [21–24]. The correlations with NDVI and EVI were minimally affected by changing the compositing period from 16-days to 32-days, demonstrating that the auto-correction feature of vegetation indices to non-vegetation signals through normalization of the reflectances is effective. On the other hand, LAI and FPAR showed improved fitting when they were composited to 32-days. Because of their inability of capturing short-term changes in GPP (i.e., temporary reductions in GPP due to adverse weather conditions), both MOD13 and MOD15 products showed weaker correlation for evergreen forests than for deciduous forests. For all the products, no correlation was found with GPP of tropical forests (Table 3). This can be attributed to the small seasonal changes observed in tropical forest leaves, which are smaller than the interannual variation in leaf phenology (e.g., Huete et al. [24]). EVI generally showed smaller RMSE than NDVI, LAI, and FPAR.

For all analyzed sites, the reference data of MOD17 GPP was most highly correlated with short-term flux tower GPP, both for deciduous and evergreen forests (Table 3). This result can be explained by the fact that MOD17 GPP parameters have been well calibrated and tested at many FLUXNET sites [31,39,67]. This result also supports the need for meteorological input data in addition to the satellite metrics for better capturing the seasonal variations in photosynthesis (Running et al. [9]).

We plotted 32-day composites of MODIS products against flux tower GPP estimates for all the sites (Figure 3). NDVI and FPAR clearly showed saturation at values of 0.75 and above. EVI showed the highest correlation with flux tower GPP, but the modest R² of 0.62 implied low accuracy in the estimation of short-term GPP through one equation using satellite metrics only. Theoretically, the four metrics should be correlated with each other. Indeed, Figure 3(a,b) shows a similar scatter even though NDVI and FPAR were calculated through very different approaches. This similarity can be explained by the well-known linear relationship between NDVI and FPAR [68]. FPAR and LAI are correlated through an exponential relationship usually described by Beer’s law, which explains the differences between the scatters of Figure 3(c) and Figure 3(d). However, even correlating the logarithmic transformation of GPP with MODIS NDVI and FPAR did not significantly increase the R².

3.2. Annual Analysis between MODIS Products and Flux Tower GPP

3.2.1. NDVI and EVI vs. Annual Mean GPP

When all sites were analyzed at a 32-day resolution, annual mean NDVI and EVI showed comparable correlations with annual GPP (Figure 4). When the GPP of tropical and non-tropical sites were compared separately to EVI, however, the correlation was much lower and non-significant (

Table 4. Summary of R² between each MODIS product and flux tower annual GPP. The stars indicate p-value < 0.05 (*) and p-value < 0.001 (**).

		Composite Period	All Sites	Non-Tropical Sites	Tropical Sites
MOD13	NDVI	16-day	0.50**	0.48**	0.30
		32-day	0.57**	0.42**	0.32
	EVI	16-day	0.56**	0.01	0.01
		32-day	0.54**	0.01	0.11
MOD15	LAI	8-day	0.78**	0.43**	0.70*
		16-day	0.83**	0.57**	0.68*
		32-day	0.88**	0.68**	0.66*
	FPAR	8-day	0.56**	0.46**	0.64*
		16-day	0.58**	0.59**	0.55*
		32-day	0.62**	0.63**	0.57*
MOD17	GPP	8-day	0.81**	0.59**	0.01

and Figure 4(b)). In contrast, NDVI showed moderate correlation both for all sites and non-tropical sites, but the linear relationship between NDVI and annual GPP collapsed at NDVI greater than 0.75. Similarly to the short-term analysis, compositing to 32-day did not significantly change the results (Table 4).

3.2.2. LAI and FPAR vs. Annual Mean GPP

LAI had the highest correlation with annual GPP, regardless of the compositing period length (Table 4 and Figure 5(a)). LAI showed high correlation with tower GPP both for non-tropical sites (R² 0.43 to 0.68) and tropical sites (R² 0.66 to 0.70). Increased compositing length tended to increase R² for non-tropical sites and decrease R² in tropical sites. The R² values for the comparisons between FPAR and annual GPP were similar, but slightly lower than those for LAI (Table 4). At FPAR greater than 0.7, flux tower GPP was generally underpredicted by the regression at non-tropical sites (Figure 5(b)). As for LAI, when compositing to longer periods, FPAR correlations increased at non-tropical sites and tended to decrease at tropical sites (Table 4).

While several metrics showed promise in estimating annual flux tower GPP, annual mean LAI derived from 32-day composites had the highest correlation, explaining 88% of measured variance. The simple model is:

$$\mathit{GPP}=615*\mathit{annual}\mathrm{\hspace{0.17em} }\mathit{mean}\mathrm{\hspace{0.17em} }\mathit{LAI}-376$$

(4)

where annual GPP is in g C/m²/year. It is noteworthy that we tried the same analysis with max LAI, but max LAI was poorly correlated with annual GPP, as previously shown by other studies [10,13].

3.3. Assessment of the Simple GPP Model with MOD17 GPP and Climate

As a broader evaluation of the use of a simple approach to estimate annual GPP, we then applied equation 4 to MOD15 LAI product globally and compared it to MOD17 GPP (both subset for forested areas, aggregated to 0.5 degrees, and calculated as a 2001 to 2006 climatology). We used the monthly 0.25 degree Collection 5 MOD15 LAI product from Boston University [37] and 1 km global MOD17 Collection 5.1 annual GPP from the University of Montana [69]. We compared the single-metric GPP with MOD17 GPP for Evergreen Needleleaf Forest (ENF), Evergreen Broadleaf Forest (EBF), Deciduous Needleleaf Forest (DNF), Deciduous Broadleaf Forest (DBF), and mixed forest. Figure 6 shows that the best-performing single-product GPP model (Equation (4)) and the more complex MOD17 GPP produced comparable climatological GPP estimates with similar standard deviations, with the average LAI-based GPP (Equation (4)) only slightly lower than MOD17 GPP, except for deciduous broadleaf forest.

Finally, to assess whether our LAI-based GPP estimates were consistent with expectations relatively to climate, we plotted global annual GPP against mean temperature, annual precipitation, and mean diurnal temperature range (DTR) derived from the CRU TS 3.1 datasets. Our results for global forested areas show that the LAI-based GPP behaves as expected, with the highest values found in conditions of high temperature and precipitation and moderate DTR, while the lowest values are found in cold and dry climates (Figure 7).

4. Discussion

4.1. Why is Annual Mean LAI the Best Satellite-Derived Estimator of Annual GPP?

LAI showed the highest correlation with annual GPP among the satellite metrics presented in this study (Table 4). This strong connection between annual GPP and annual mean LAI has previously been documented through ground measurements by Kira et al. [70] and Waring et al. [71], who showed high correlations between GPP and Leaf Area Duration (LAD = growing season LAI × the number of growing season months). The relationship can be explained biophysically with the observation that larger canopies and longer growing seasons allow more carbon to be assimilated. Conversely, as more carbon is fixed as GPP, more carbon is allocated to leaf production. Although annual mean LAI correlates well with annual GPP, seasonal LAI correlates poorly with seasonal GPP because it does not respond to short term stresses such as anomalous heat waves.

The high correlation across different forest types can be explained by the functional convergence hypothesis [72]. The functional convergence hypothesis states that plants tailor their photosynthetic capacity to the amount of the most limiting resource so that optimal use is made of all resources [72]. The hypothesis has been well examined at the leaf level and is used to explain the convergence of leaf traits, such as the high correlation between photosynthetic capacity and leaf nitrogen content [73,74]. At the landscape or global level, the functional convergence hypothesis can be extended to explain the high correlation between absorbed PAR (APAR) and net primary production, and leads to assuming a constant LUE (Field, 1991). Goetz et al. [75] concluded that, because of the wide variability of carbon use efficiency (NPP/GPP), it is the LUE of GPP that should converge rather than the LUE of NPP. Since APAR is a function of LAI, the functional convergence theory should apply also for the relation between annually integrated GPP and LAI. Hence, we expect the annual GPP and LAD to be tightly correlated, even though over shorter time periods there can be a discrepancy between GPP and LAI, as shown in Table 3. Similar results were found when calculating GPP and LAI using prognostic ecosystem models [65].

Vegetation indices plots (Figure 3) showed non-linear relationship with GPP and a very different scatter from LAI plots (Figure 4). Vegetation indices do not monotonically increase with LAI mainly because of (1) canopy structural parameterization in the radiative transfer model of the MOD15 algorithm and (2) higher sensitivity of LAI to the reflectances in the saturation range of vegetation indices [35].

4.2. Why was the Performance of EVI Different between the Short-Term and Long-Term Spatial Analysis?

Our results suggest that simple models are viable approaches for estimating GPP, but that fundamentally different approaches should be used depending on the spatio-temporal scale of the analysis. Among the four satellite-derived metrics, EVI showed the highest correlation with short-term GPP (Table 3) but not with annual GPP (Table 4). The reason for the higher correlation of EVI than NDVI with short-term measurements of GPP can be attributed to (1) its broader dynamic range compared to NDVI, which is due to non-saturation at high LAI, (2) empirical auto-correction for atmospheric contamination of the reflectance, and (3) consistent retrievability without need of backup algorithms. Other studies also showed that EVI can track the seasonal change in GPP and other phenological events better than NDVI for deciduous forest [21].

However, when assessed at the level of either non-tropical or tropical sites to annual mean GPP, EVI had almost no correlation with annual GPP, while NDVI had significant correlation (Table 4). This result indicates that the conversion factor from EVI to GPP varies across forest types. That is not surprising considering that vegetation indices are a metric that is only theoretically proportional to leaf amount at any single point [76]. There is no clear biophysical justification for which the GPP/VI ratio should be constant among different biomes. Huete et al. [77] also reported that the slopes of the linear regression between monthly EVI and monthly GPP depended on biome type. In this study, the correlations between annual maximum VIs and annual maximum GPP were 0.47 (p < 0.001) for NDVI and 0.17 for EVI, indicating that the GPP/EVI ratio was site-specific rather than constant. This explains why annual mean EVI showed no correlation with annual GPP (Table 4).

To check if these results depended on the assumptions made for winter LAI in the presence of snow, and hence the calculation of the mean values, we set the value of the snow-covered pixels to the annual minimum rather than zero. While this increased the correlations of NDVI and EVI (R² = 0.63 for NDVI and R² = 0.62 for EVI), it did not significantly affect the correlation of LAI and FPAR. However, these results did not change how well those metrics were correlated because increases in the winter values proportionally increased the annual mean values. Thus, our results were not dependent on setting the snow values to zero.

4.3. Assessment of the LAI Empirical Relationship with Other Datasets

Our analysis shows that the simple model for estimating annual GPP from LAI has comparable accuracy—at global scales—to the MOD17 GPP product (Figure 6). LAI-based estimates are also consistent with the expected influence of climate on GPP (Figure 7), with estimates increasing with precipitation and temperature and average DTR. High DTR is indicative of dry conditions, while lower DTR is associated with low radiation. Both conditions are unfavorable for vegetation photosynthetic activity, as are extreme cold and low precipitation [78]. The fact that our GPP estimates broadly follow these patterns lends additional confidence that the simple approach does not generate spurious predictions. Our estimates of GPP for both tropical (2,332.9 g C/m²/year) and non-tropical (755.6 g C/m²/year) forest are also comparable to the ensemble median of 6 models by Beer et al. [2], further supporting that our results are reasonable.

At finer temporal resolutions, or for interannual variability, important differences may exist between the LAI-based GPP and MOD17, and we do not argue the two approaches are necessarily interchangeable. Indeed, in spite of errors and/or biases in input meteorology, upstream MODIS products, and structural issues with the algorithm itself [31,79], MOD17 is more highly correlated to the short-term measured GPP than any of the single MODIS products. But for research focused on broad-scale patterns of annual GPP variability, the simple model does appear to provide a useful alternative, and one that is consistent with theoretical principles developed decades ago.

4.4. Effects of Uncertainties in MODIS Products on Their Relationship with Annual GPP

In Section 2 we discussed the uncertainty introduced in the MODIS products by possible cloud contamination of the pixels and how to minimize it by temporally compositing the datasets.

In the case of MOD15, the propagation of algorithm uncertainties presents additional issues in the estimation of GPP. As described in Section 2.1, the algorithm uncertainty of LAI and FPAR was estimated through ground validation studies. Here, we assess the influence of MOD15 uncertainty on the GPP estimate through comparison with CYCLOPES [80], the other existing global LAI dataset. CYCLOPES LAI has a 1/112 degree spatial resolution and 10-day temporal sampling. The algorithm utilizes a neural network approach to retrieve LAI from SPOT VEGETATION reflectance data. We processed the CYCLOPES dataset to 30-day composites similarly to the MOD15 LAI 30-composites and calculated the correlation with annual GPP. In contrast to the high correlation between MODIS annual mean LAI and GPP, CYCLOPES annual mean LAI had an R² of only 0.02 when compared to the annual flux tower GPP. Thus, annual GPP estimates appear to be sensitive to the uncertainty of the LAI products. The differences between MOD15 and CYCLOPES were mostly pertaining to: (a) Definition of LAI, with CYCLOPES measuring effective LAI while MODIS algorithm retrieves true LAI accounting for leaf clumping; (b) Nadir view TOA reflectances are input to the CYCLOPES LAI algorithm, while spectral surface reflectances with varying view angles are input into the MODIS algorithm; (c) Early saturation of the CYCLOPES LAI product, with maximum values limited to LAI = 4 while MODIS retrieves LAI of values up to 7 (includes both main algorithm and backup algorithm retrievals); (d) Classification errors due to differences in land cover mapping. These algorithm differences suggest caution when applying the proposed GPP estimation method to algorithm dataset different than MOD15.

Another potential source of uncertainty could be caused by the spatial resolution of the datasets. We used 250 m resolution data of MOD13, 1 km resolution data of MOD15 and MOD17, and flux tower point GPP data. We avoided the uncertainty caused by within-pixel spatial heterogeneity by applying a strict selective process to use only data for towers located over homogenous land covers, as described in Section 2. Also, Sims et al. [23] showed that comparisons between flux tower GPP and EVI were largely insensitive to the selection of a 3 km by 3 km grid versus a single 1 km pixel. Thus, we believe that our analysis was hardly affected by spatial resolution effects.

Unfortunately, the number of tropical forest flux sites was very small compared to boreal or temperate forest sites. In this study, we used 4 sites in South East Asia and 1 site in South America. More data for tropical forest are needed to improve the GPP estimation over these regions, likely altering the regression for the tropics.

5. Conclusions

We analyzed the correlation of MODIS NDVI, EVI, LAI and FPAR products with short-term and annual GPP measured at eddy covariance flux towers. Our results show that EVI is useful for the analysis of short-term variations in tower-estimated GPP (i.e., at the time scale of weeks to months) in agreement with other studies. Some diagnostic models started to utilize the EVI-GPP seasonal correlation, but caution should be exercised when mapping geographic patterns of GPP using EVI, because the EVI/GPP relationship is not constant across forest types. In our study, annual mean MODIS LAI was the best metric for estimation of the spatial patterns of annual GPP (GPP in g C/m²/year = 617 × LAI − 376). Annual mean NDVI and EVI were less correlated to tower-estimated GPP, but annual mean NDVI was a better predictor than annual mean EVI of annual GPP in non-tropical areas, suggesting that the GPP/EVI ratio is site specific. Extrapolating the relationship between annual mean LAI and annual GPP at the tower sites to forests globally produced reasonable GPP estimates that were also consistent with known climatic influences on GPP. The results allowed the rediscovery of the ground measured LAI and tower-estimated GPP relationship, as explained by strongly interrelated biophysical processes and the functional convergence hypothesis.

There are several concerns related to the uncertainty of MODIS products. These analyses were based on MODIS products, and different satellite products may be not able to get such high correlation between annual mean LAI and annual GPP. For the uncertainty caused by different resolution analysis, we minimized the influence on GPP estimate by carefully choosing the dataset not to include contaminated or heterogeneous data.