Search tips
Search criteria 


Logo of scirepAboutEditorial BoardFor AuthorsScientific Reports
Sci Rep. 2017; 7: 44046.
Published online 2017 March 8. doi:  10.1038/srep44046
PMCID: PMC5341022

Robustness and Uncertainties of the “Temperature and Greenness” Model for Estimating Terrestrial Gross Primary Production


Terrestrial gross primary production (GPP) plays a vital role in offsetting anthropogenic CO2 emission and regulating global carbon cycle. Various remote sensing approaches for estimating GPP have attracted considerable scientific attentions, yet their robustness and uncertainties remain unclear. Here we evaluate the performance of the “temperature and greenness” (TG) model, a representative remote sensing model in estimating GPP, using the global FLUXNET GPP based on parameter sensitive analysis and optimization strategies. The results show that the minimum (xn) and optimum (xo) temperatures for photosynthesis are sensitive parameters but maximum temperature (xm) not. Optimized xn and xo differ largely from their defaults for more than half of 12 plant functional types (PFTs). Parameter optimization significantly improves the TG’s performance in forest ecosystems where temperature or solar radiation has significant contribution to GPP. For water-limited ecosystems where GPP are strongly dependent of EVI and EVI are sensitive to precipitation, parameter optimization has limited effects. These results imply that the TG model, and most likely for other kind of GPP models using same methodology, can’t be significantly improved for all PFTs through parameter optimization only, and other key climatic variables should be incorporated into the model for better predicting terrestrial ecosystem GPP.

Terrestrial gross primary production (GPP) is the major driver of global carbon cycle and it plays an important role in regulating the concentration of CO2 of the atmosphere by partly offsetting anthropogenic CO2 emissions1. However, direct measurements of GPP are not available, because no observational techniques are ready to quantify this process at the right spatial scale2.

Quantification of GPP at ecosystem level is mostly inferred from the measurements of net ecosystem production (NEP) between terrestrial ecosystems and the atmosphere using the eddy covariance (EC) equipment3. More than 950 site-years EC data have been archived in the international network of FLUXNET4 during the past three decades, which make the estimation of GPP possible at site level. However, estimate of GPP at a larger (for example regional, national or global) scale was difficult5. To retrieve GPP estimates at larger scales, site level of EC-inferred GPP have to be scaled up to spatial domains based on empirically statistical methods5,6, such as artificial neural networks7 or ensemble model trees5, or semi-empirical models including light use efficiency8 or water-use efficiency9 approaches. A common disadvantage of these approaches is the strong dependency on environmental (vegetation, soil or meteorological) variables6.

Another commonly used approach to quantify GPP based on land surface models (LSMs) or ecosystem models at different scales10. These models have explicit modules to simulate carbon cycle processes by describing plant physiological behaviour in relation to soil and atmospheric processes11. Models can be implemented alone (offline)12,13,14 or coupled with climate model (online, also termed Earth System Models)15,16,17. The significant advantage of these carbon models relies in the continuousness in both space and time, so that they have been often used to detect the trend and inter-annual variability of GPP in long term period and at larger scales18,19. However, LSMs require multiple driving data including meteorological variables, vegetation and soil maps, which are highly spatial heterogeneities and uncertainties.

Alternatively, GPP at larger scale can be estimated by the approach using satellite measurements of vegetation parameters that are directly related to plant photosynthesis (i.e., GPP)20. These parameters mainly include photosynthetically active radiation (PAR), normalized difference vegetation index (NDVI) or enhanced vegetation index (EVI), and leaf area index (LAI). Integrating one or more of these photosynthesis-related remote sensing measurements into an empirical or semi-empirical model was used to predict GPP. Amongst existed models, a classical type is termed light use efficiency (LUE) model, which describes GPP as a multiple product of PAR, fraction of PAR absorbed by the vegetation (fPAR), potential LUE (LUEmax) and the environmental constraints on LUEmax21,22,23,24,25,26. LUE-type models also require ancillary environmental variables such air temperature, vapor pressure deficit, soil moisture or canopy water content to constrain LUEmax26. In LUE-type models, fPAR was generally calculated as a linear function of NDVI27. However, NDVI was found to be easily saturated at moderate to high vegetation coverage28,29,30. Development of EVI, in certain extent, reduced the severity of early saturation31. EVI was found to perform better than NDVI in predicting GPP in high dense vegetation, such as evergreen and deciduous forests32,33 and cropland as well34. Based on strong correspondence of GPP on EVI over one or two week intervals35,36, a cluster of EVI-based models was developed to predict GPP32,33,35,36,37. Amongst the EVI-based models, temperature and greenness (TG) model received much attention due to its simplicity which uses purely remote sensed land surface temperature (LST) and EVI, free of any ancillary meteorological data, as drivers37. In the TG model, LST was considered to be representative of climatic variables and EVI include information of photosynthetic potential37. In addition, the TG model only consists of three temperature-related parameters (see Methods), which represent the minimum (xn), optimum (xo) and maximum (xm) LST values on photosynthesis, respectively. The TG model has been widely evaluated against with site level of derived GPP from EC38,39 and inter-compared with other LUE models39, it showed significant advantages in limited ecosystem types. However, the robustness of such a simply model in estimating GPP across a wide spectrum of ecosystems and associated uncertainties from model parameters have not yet been investigated.

Therefore, the first objective of this research is to conduct a comprehensive evaluation of the TG model for predicting GPP using the EC-derived GPP. The second objective is to deploy the potential of parameter optimization method in improving the performance of the TG model. Finally, this paper aims at investigating the robustness and uncertainties of the TG model among different biomes through introducing other climatic factors (e.g., precipitation and solar radiation) into a linearly multiple regression model. Such evaluation and investigation provided a common methodology to any other LUE-type models for predicting GPP, although the TG model was adopted as a case in this research.


Parameter sensitivities of the TG model

The TG model consists of three parameters, xn, xo and xm, which represent minimum, optimum and maximum temperature for plant photosynthesis37. Their values were set as 0, 30, and 50 °C, respectively, in the TG model37. The three parameters are referred to LST, rather than air temperature, which is more commonly used to drive productivity models37. The values of the three temperature-related parameters may differ largely across PFTs, so we set wide ranges of the three parameters (Table 1) to investigate their sensitivities using Morris method40. From the Morris screening method, two indices μ* and σ measure the influence of each parameter and its degree of involvement in non-linearities and interactions with other parameters, respectively. None of σ value was significantly larger than 2 μ* for all parameters and all PFTs, indicating that limited interactions and/or non-linearities existed (Fig. 1). Then we ranked the relative degree of sensitivities of three parameters based on their μ* values for each PFT and a larger μ* was considered as a more influential parameter (Fig. 1a). It is clearly that mostly influential parameters differ across PFTs. Parameter xn and xo were identified as mostly influential parameters of the TG model for all forests (including DBF, DNF, EBF, ENF and MF) and GRA, while xo was selected as the relatively influential parameters for other remaining PFTs (including CSH, CRO, WSA, WET, SAV and OSH). However, xm was not identified as the mostly influential parameter for any PFT (Fig. 1 and Fig. S1).

Figure 1
Mean elementary effects (μ*) (a) and standard deviation of elementary effects (σ) (b) derived from the Morris screening analysis for twelve PFTs. Each axis of the radar plot corresponds to the values of μ* or σ and different ...
Table 1
Default, ranges and optimized values for the three parameters of the TG model for 12 PFTs.

Optimization of parameters of the TG model

Theoretically, only those parameters selected as mostly influential by the Morris screening analysis were applied to MCMC for parameter optimization. Because there were only three parameters in the TG model and applying all three parameters to MCMC will not exert great burdens on computation, we included all three parameters (xn, xo and xm) during MCMC optimization.

To guarantee the validity of the optimized parameters from MCMC, convergence of the posterior parameters must be assured. A mostly direct way to diagnose the convergence is visually plotting the density of the posterior parameter. More like a Gaussian distribution of the posterior parameter, more possible the convergence can be assured. As shown in Fig. 2a, parameter xn was identified as convergent for ENF, EBF, DBF and WET. Parameter xo was considered as convergent for all PFTs except CRO, DNF and OSH (Fig. 2b). However, xm was identified as unconvergent during the range of 40–65 °C for all PFTs (Fig. 2c) due to their uniform or gamma distributions. The convergence of the posterior parameters optimized by MCMC agreed with the results diagnosed by running mean of the posterior parameters (Fig. S2). Overall, posterior distribution of parameter xo showed convergence in nine of twelve PFTs (excluding CRO, DNF and OSH), and that of parameter xn showed convergence in four of twelve PFTs (including ENF, EBF, DBF and WET) but posterior distribution of xm from MCMC did not show convergence for any PFT. It is worth mentioning that the parameters identified as convergent from MCMC (Fig. 2) were consistent with those parameters identified as mostly influential by the Morris screening analysis (Fig. 1a), which also indicated the advantage of applying all three parameters to MCMC.

Figure 2
Density of the posterior parameters for three parameters (a) for xn, (b) for xo and (c) for xm of the TG model. Lines with different colors represent different PFTs. Three vertically dashed lines represent the default values for parameter xn, xo and x ...

The optimized parameter value for those PFTs which were sensitive to xn were significantly different from the default value of xn (0 °C), and had wide ranges from −10 °C in ENF to 5 °C in WET (Fig. 2a). For parameter xo, the optimized value by MCMC differed largely across PFTs (Fig. 2b). In EBF, the optimized value for xo was 15 °C, half of the default value (30 °C) of the TG model. In WET, the optimized xo was 22.5 °C. In DBF, ENF, GRA and MF, the optimized xo by MCMC were about 25 °C. While, in the remaining three PFTs (CSH, SAV and WSA) which were sensitive to xo, the optimized values of xo ranged from 27.5 to 32.5 °C, and were close to the default xo (30 °C) of the TG model. Wide ranges of the optimized xo for LST by MCMC provided good support for previous report that optimum air temperature for canopy carbon flux varied from 7.5 to 30 °C and were highly correlated with mean summer temperature41. For parameter xm, no optimized value can be derived by the MCMC method during the range of LST between 40 and 65 °C (Fig. 3c), which was in good agreement with the analysis from the Morris screening method (Fig. 1).

Figure 3
Effects of MCMC optimization on performance of the TG model in twelve PFTs.

Effects of MCMC optimization on the TG model

Figure 3 shows the performance of the TG model before and after MCMC optimization. The default TG model (before optimization) performed large difference in different PFT. Correlation coefficient (R) between the observed GPP and the simulated GPP by the default TG model ranged from 0.55 in EBF to 0.96 in DNF, but were larger than 0.7 in ten of twelve PFTs (except EBF and WET), which indicated that the default TG model was able to explain more than 60% of variance of the observed GPP at monthly scale. After applying MCMC to the TG model, correlation coefficient (R) between data and model did not significantly increase in all twelve PFTs, but both normalized standard deviation (STD) and root mean squared difference (RMSD) decreased in EBF and ENF. STD decreased from 1.32 to 0.90, and RMSD decreased from 1.12 to 0.94 for EBF. For ENF, STD decreased from 1.01 to 0.86, and RMSD decreased from 0.67 to 0.62. Obviously, the effects of MCMC optimization on performance of the TG model were caused by the difference in the values of parameters xn and/or xo of the model (Fig. 2a and b). Optimized TG model by MCMC did not cause the increase in R between data and model, but caused the decrease in normalized STD and RMSD in EBF and ENF, which were also attributed by the nature that the maximum likelihood estimation was used as a cost function in MCMC.

Relative Contributions of EVI and Climatic Variables on GPP

It is clear that EVI had dominated effects on GPP in all PFTs except EBF in which RAD was found to have comparable effects on GPP (Fig. 4). Individual effect of EVI was able to contribute 13% and 28% of GPP variation in EBF and WET, respectively. In the remaining ten PFTs, EVI could contribute 39% (ENF) to 77% (CSH) of GPP variations. When introducing LST as a predictor to estimate GPP in MLR, LST showed its significant effects on GPP in forest PFTs (including DBF, DNF, EBF, ENF and MF) but insignificant in non-forest PFTs. Similarly, RAD was found to have significant effects to estimate GPP in forest PFTs (excluding DNF with limited measurements) which are mostly energy limited ecosystems and moderate effects on CRO and WSA as well. In contrast, PRC was found to have significant effects on GPP in WSA, SAV, OSH, CSH and CRO, in which GPP were mostly possible to be limited by water42.

Figure 4
Relative contributions of EVI (green), land surface temperature (LST, red), precipitation (PRC, blue) and radiation (RAD, purple) on GPP in twelve plant functional types (PFTs).

Given considering the four key variables (EVI, LST, PRC and RAD) together as predictors for estimating GPP using MLR, we found that the MLR using the four variables had comparable capability to estimate GPP as the TG model in ten out of twelve PFTs except EBF and WET. In EBF, MLR using four variables had worse ability to estimate GPP but better in WET.


Although the TG model was reported to have high accuracy of predicting GPP across a wide range of ecosystem types, previous evaluations were mostly focused on limited ecosystem types or few sites of data for an ecosystem type37,39. This study made utilize of a global FLUXNET GPP data to demonstrate that the TG model performed well (R2 > 0.7) for 7 of 12 PFTs (CSH, DFB, DNF, MF, OSH, SAV, and WSA), and moderately (0.5 < R2 < 0.7) for 3 PFTs (CRO, GRA and ENF), but poorly (R2 < 0.5) for EBF and WET. So it is necessary to understand why the TG model performed differently across PFTs.

Overall, the TG model predicted GPP well compared with the flux tower GPP in WSA and SAV where vegetation are sparse, which obviously conflicted with previous understanding37. It was deemed that solar elevation angel strongly affected the EVI values when vegetation is sparse37 which are cases for WSA and SAV. Our analysis showed that EVI was still well correlated with GPP in SAV and WSA, indicating that vegetation fraction coverage may be not an important source of errors. More importantly, EVI was proven to be very sensitive to PRC in SAV and WSA43, so LST was found to be have less contribution than PRC on GPP (Fig. 4). This is also evidenced in CSH, GRA and OSH. However, in DBF, EBF and ENF, RAD was found to have considerable contributions on GPP, indicating that GPP of forest ecosystems may be limited by RAD44. In WET, GPP was principally influenced by RAD45 and non-climatic variables such as leaf area index and water table depth45. As the number of GPP data for DNF was relatively limited, RAD did not show significant contribution on GPP.

The performance of a model can be determined by the model structure (behaviour), model parameters and input10. In this study, we mainly focused on the effects of model parameters on the performance of the TG model by optimizing the parameter based on MCMC. In terms of parameter values in the TG model, optimized values for both xo and xn differed largely from their default ones for more than half of 12 PFTs. However, parameter optimization using MCMC significantly improved the model’s performance for EBF and ENF only, where temperature and solar radiation are critically important for GPP. For other ecosystem types, such as DBF and WET, MCMC did not significantly improve the performance identified by a Taylor diagram, although the optimized parameters (xo and xn) values were largely different from their default. This suggests that the performance of the TG model can’t be improved by optimizing its parameters only, which also reflected that the TG model has large potential to improve in model structure and process. For example, solar radiation should be integrated into the model for forested ecosystems, and precipitation should be considered for water-limited ecosystems (WSA, SAV, and GRA). However, in some water-limited ecosystems, EVI was found to be very sensitive to precipitation43. Consequently, the TG model performed very well for these water-limited ecosystems (Fig. 2). However, we should note that in some PFTs, MCMC did not really improve the model performance but still induce a large discrepancy in the values of parameters. This suggests that the optimized parameter values may not correspond to true underlying values. Another reason is that the TG model was jointly determined by all parameters, rather than any individual one. Therefore, the optimized value for a single parameter by MCMC should not be over-interpreted.

Similar to previously studies using vegetation index to predict GPP46,47, LUE-type models showed robustness for some ecosystems and weakness as well for other ecosystems or when drought events presented. The essential causes were that the algorithms describing the environmental regulations, particularly water stress on GPP built in those models are largely different. By comparing the TG model (with or without parameter optimization) with a linearly multiple regression model, other climatic variables, such as radiation or precipitation, have significant effects on GPP in some ecosystems. In the TG model, LST was considered to be representative of climatic conditions (air temperature, vapor pressure deficit, or radiation). However, this prerequisite for the TG model does not hold for all ecosystems, evidenced by largely different “relative contribution” of each climatic variable to GPP (Fig. 4). This implies that the TG model, and most likely for other kind of GPP models using same methodology, should incorporate other key climatic variables for further improving the accuracy in simulating terrestrial ecosystem GPP.


Global GPP dataset

The analyses performed in this study are based on monthly GPP data from 155 flux towers, consisting of a total of 624 site-year datasets and representing a world wide spectrum of biomes and climate regions. These data covered 12 ecosystem types including cropland (CRO), closed shrubland (CSH), deciduous broadleaf forest (DBF), deciduous needleleaf forest (DNF), evergreen broadleaf forest (EBF), evergreen needleleaf forest (ENF), grassland (GRA), mixed forest (MF), open shrubland (OSH), savanna (SAV), wetland (WET), and wood savanna (WSA). Geophysical location of the flux tower sites are shown in the Fig. S3, and the number of each ecosystem type is listed in the Table S1. The majority of data were obtained from global FLUXNET (La Thule level 2 or level 4 products)48, and the rest of data were obtained directly from the site researchers. GPP data were inferred from net ecosystem production (NEP) observed by eddy covariance. Only site-years with small gaps (i.e., individual gaps in NEP of less than 5% of the entire annual record) were selected except in certain ecosystems of the boreal region where only growing season data were available. Missing values of NEP (and GPP) were gap-filled and partitioned by a publicly available marginal distribution sample (MDS) online tool3.

Satellite data

MODIS monthly enhanced vegetation index (EVI) (MOD13A3.005) and land surface temperature (LST) (MOD11C3.005) products for February 2000 to 2013 were obtained from the USGS repository ( Both EVI and LST dataset were produced globally over land at 1-km resolution. Due to the difficulty of precisely co-locating the pixels that directly correspond to the footprint of an EC flux site, a central 3 × 3 km window surrounding the flux tower was used to extract mean EVI and LST time series49. To reduce noise and uncertainties in the MODIS EVI and LST time series at each site, the singular spectrum analysis (SSA) was employed. In the monthly EVI and LST time series, a window length of 37 (i.e. 37 months) and 6 leading components were found to best captured the periodicity and simultaneously reduced random noises during reconstruction. To further minimize the contamination, the snow/ice flag in MOD13A3 quality assurance field was first used to remove the snow-covered EVI values; and then EVI values were further screened for effects of cold temperature (LST < −2.0 °C) using the MODIS LST50.

Temperature and greenness (TG) model

The temperature and greenness (TG) model simulates GPP as a proportionally linear function of the product of scaled EVI (EVIscaled) and LST (LSTscaled)37, i.e.,

An external file that holds a picture, illustration, etc.
Object name is srep44046-m1.jpg

where m is a scalar with unit of gC m−2 day−1.

The EVI scalar is simply defined as:

An external file that holds a picture, illustration, etc.
Object name is srep44046-m2.jpg

where 0.1 was the EVI value when GPP drops to zero.

Another scalar LSTscaled was defined as:

An external file that holds a picture, illustration, etc.
Object name is srep44046-m3.jpg

where min represent the minimum value of the two items in the square brackets, xn, xo and xm are three parameters, and represent minimum, optimum and maximum temperature for plant photosynthesis37. Their default values were set as 0, 30, and 50 °C, respectively37. In the TG model, LST represent climatic conditions in some extent, and EVI is used as a measure of vegetation greenness.

Ranking of parameter sensitivities

Although there are only three parameters in the TG model, it is still necessary to understand the sensitivity of the model to each parameter. Here we use the Morris method40,51 to screen out the most important parameters. The Morris screening is based on a one-at-a-time (OAT) approach, in which only one parameter is varied between two runs, allowing for calculating a local partial derivative of the output variable with respect to the input parameter. This impact of input parameter on the output is called elementary effect. The Morris method calculates the elementary effects of each parameter in several locations of the parameter space several times, so it is generally considered to be a “global” screening method. The mean (μ*) and standard deviation (σ) of all elementary effects for each parameter can be used to be representative of the behavior of this parameter in its entire range of variation. The two indices μ* and σ measure the influence of the parameter on output variable and its degree of involvement in non-linearities and/or interactions with other parameters, respectively.

In this study, the R package “sensitivity”52 was used to implement the Morris method for the TG model across 12 ecosystem types.

Bayesian approach for parameter optimization

Parameter optimization was performed within a Bayesian framework, which treats the unknown model parameters x as random variables with the posterior probability density function (pdf) p(x|y) given by

An external file that holds a picture, illustration, etc.
Object name is srep44046-m4.jpg

where p(x) denotes the prior distribution of x and p(y|x) signifies the likelihood function of x. Bayesian inference was performed using Markov Chain Monte Carlo (MCMC)53. For the application of Bayes theorem, an efficient algorithm called the Metropolis-Hastings algorithm54,55 was used to the posterior distribution p(x|y). The prior probability density function was first specified by providing a set of limiting intervals for the parameters. The likelihood function was then constructed assuming that errors in the observed data followed Gaussian distributions. In total, 105 iterations of probability distribution for each PFT were carried out to best satisfy convergence.

Statistical analysis

The Taylor diagram56 was used to describe the effects of parameter optimization on the performance of the TG model. The performance of the TG model was displayed in a single diagram, featured by correlation coefficient (R), standard error (STD), and root mean square difference (RMSD). The higher the R and the smaller the STD and RMSD, the better the agreement between model and data is. When comparing two simulations with different parameter values, i.e. the default TG model or the TG model after MCMC optimization, the longer the distance between the two simulation points, the greater effect of optimization on the model’s performance. To make the performance of the TG model across ecosystems, the Taylor diagram was shown in a normalized way, so “observations” are displayed in the point (R = 1, STD = 1, RMSD = 0).

To investigate why the TG model performed differently amongst PFTs, relative weight analysis (RWA) was used to quantify the relative contributions of EVI and three key climatic variables, LST, precipitation (PRC) and radiation (RAD) on GPP, supplementary to multiple linear regression (MLR) analysis. The function “rlw” in the R package “yhat” was used to implement RWA, and the R package “plotrix” was used to plot Taylor diagram.

Additional Information

How to cite this article: Dong, J. et al. Robustness and uncertainties of the “temperature and greenness” model for estimating terrestrial gross primary production. Sci. Rep. 7, 44046; doi: 10.1038/srep44046 (2017).

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Material

Supplementary Information:


This work was supported by the National Natural Science Foundation of China (Grant No. U1403382). We thank the China Scholarship Council for the financial support to Jiaqi Dong.


The authors declare no competing financial interests.

Author Contributions L.L. contributed to design this study. J.D. and H.S. performed the analysis. J.D., H.S. and L.L. wrote the manuscript. X.C., G.L. and Q.Y. contributed significantly to the discussion of results and manuscript refinement.


  • Cox P. & Jones C. Illuminating the Modern Dance of Climate and CO2. Science 321, 1642–1644 (2008). [PubMed]
  • Ma J., Yan X., Dong W. & Chou J. Gross primary production of global forest ecosystems has been overestimated. Scientific Reports 5, 10820 (2015). [PMC free article] [PubMed]
  • Reichstein M. et al. . On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algorithm. Global Change Biol. 11, 1424–1439 (2005).
  • Williams M. et al. . Improving land surface models with FLUXNET data. Biogeosci. 6, 1341–1359 (2009).
  • Jung M. et al. . Global patterns of land-atmosphere fluxes of carbon dioxide, latent heat, and sensible heat derived from eddy covariance, satellite, and meteorological observations. J. Geophys. Res.: Biogeosci. 116, G00J07 (2011).
  • Beer C. et al. . Terrestrial Gross Carbon Dioxide Uptake: Global Distribution and Covariation with Climate. Science 329, 834–838 (2010). [PubMed]
  • Beringer J. et al. . An introduction to the Australian and New Zealand flux tower network – OzFlux. Biogeosci. 13, 5895–5916 (2016).
  • Monteith J. L. & Moss C. J. Climate and the Efficiency of Crop Production in Britain [and Discussion]. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 281, 277–294 (1977).
  • Beer C. et al. . Temporal and among-site variability of inherent water use efficiency at the ecosystem level. Globle Biogeochem. Cycl. 23, GB2018 (2009).
  • Prentice I. C., Liang X., Medlyn B. E. & Wang Y. P. Reliable, robust and realistic: the three R’s of next-generation land-surface modelling. Atmos. Chem. Phys. 15, 5987–6005 (2015).
  • Bonan G. B. et al. . Improving canopy processes in the Community Land Model version 4 (CLM4) using global flux fields empirically inferred from FLUXNET data. J. Geophys. Res.: Biogeosci. 116, G02014 (2011).
  • Krinner G. et al. . A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system. Globle Biogeochem. Cycl. 19, GB1015 (2005).
  • Dai Y. et al. . The Common Land Model. Bull. Am. Meteorol. Soc. 84, 1013–1023 (2003).
  • Sitch S. et al. . Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model. Global Change Biol. 9, 161–185 (2003).
  • Zeng X., Shaikh M., Dai Y., Dickinson R. E. & Myneni R. Coupling of the Common Land Model to the NCAR Community Climate Model. J. Clim. 15, 1832–1854 (2002).
  • Wu T. et al. . Global carbon budgets simulated by the Beijing Climate Center Climate System Model for the last century. J. Geophys. Sci.: Atmosph. 118, 4326–4347 (2013).
  • Law R. M. et al. . The carbon cycle in the Australian Community Climate and Earth System Simulator (ACCESS-ESM1) – Part 1: Model description and pre-industrial simulation. Geosci. Model Dev. Discuss. 2015, 8063–8116 (2015).
  • Sitch S. et al. . Recent trends and drivers of regional sources and sinks of carbon dioxide. Biogeosci. 12, 653–679 (2015).
  • Anav A. et al. . Spatiotemporal patterns of terrestrial gross primary production: A review. Rev. Geophys. 53, 785–818 (2015).
  • Myneni R. B. et al. . Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data. Remote Sens. Environ. 83, 214–231 (2002).
  • Landsberg J. J. & Waring R. H. A generalised model of forest productivity using simplified concepts of radiation-use efficiency, carbon balance and partitioning. For. Ecol. Manage. 95, 209–228 (1997).
  • Potter C. S. et al. . Terrestrial ecosystem production: A process model based on global satellite and surface data. Globle Biogeochem. Cycl. 7, 811–841 (1993).
  • Xiao X. et al. . Satellite-based modeling of gross primary production in an evergreen needleleaf forest. Remote Sensing of Environment 89, 519–534 (2004).
  • Running S. W. et al. . A continuous satellite-derived measure of global terrestrial primary production. Bioscience 54, 547–560 (2004).
  • Prince S. D. & Goward S. N. Global Primary Production: A Remote Sensing Approach. J. Biogeogr. 22, 815–835 (1995).
  • Yuan W. et al. . Deriving a light use efficiency model from eddy covariance flux data for predicting daily gross primary production across biomes. Agric. For. Meteorol. 143, 189–207 (2007).
  • Myneni R. B. & Williams D. L. On the relationship between FAPAR and NDVI. Remote Sens. Environ. 49, 200–211 (1994).
  • Hatfield J. L. et al. . Leaf-area estimates from spectral measurements over various planting dates of wheat†. Int. J. Remote Sens. 6, 167–175 (1985).
  • Asner G. P., Scurlock J. M. O. & Hicke A., J. A. Global synthesis of leaf area index observations: implications for ecological and remote sensing studies. Global Ecol. Biogeogr. 12, 191–205 (2003).
  • Chen P.-Y., Fedosejevs G., Tiscareño-LóPez M. & Arnold J. G. Assessment of MODIS-EVI, MODIS-NDVI and VEGETATION-NDVI Composite Data Using Agricultural Measurements: An Example at Corn Fields in Western Mexico. Environ. Monit. Assess. 119, 69–82 (2006). [PubMed]
  • Huete A. et al. . Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sens. Environ. 83, 195–213 (2002).
  • Xiao X. M. et al. . Modeling gross primary production of temperate deciduous broadleaf forest using satellite images and climate data. Remote Sens. Environ. 91, 256–270 (2004).
  • Xiao X. M. et al. . Satellite-based modeling of gross primary production in an evergreen needleleaf forest. Remote Sens. Environ. 89, 519–534 (2004).
  • Zhang Q. Y. et al. . Estimation of crop gross primary production (GPP): I. impact of MODIS observation footprint and impact of vegetation BRDF characteristics. Agric. For. Meteorol. 191, 51–63 (2014).
  • Rahman A. F., Sims D. A., Cordova V. D. & El-Masri B. Z. Potential of MODIS EVI and surface temperature for directly estimating per-pixel ecosystem C fluxes. Geophys. Res. Lett. 32, 4 (2005).
  • Sims D. A. et al. . On the use of MODIS EVI to assess gross primary productivity of North American ecosystems. J. Geophys. Res. - Biogeosci. 111, G04015 (2006).
  • Sims D. A. et al. . A new model of gross primary productivity for North American ecosystems based solely on the enhanced vegetation index and land surface temperature from MODIS. Remote Sens. Environ. 112, 1633–1646 (2008).
  • Dong J. et al. . Comparison of four EVI-based models for estimating gross primary production of maize and soybean croplands and tallgrass prairie under severe drought. Remote Sens. Environ. 162, 154–168 (2015).
  • Wu C., Munger J. W., Niu Z. & Kuang D. Comparison of multiple models for estimating gross primary production using MODIS and eddy covariance data in Harvard Forest. Remote Sens. Environ. 114, 2925–2939 (2010).
  • Morris M. D. Factorial Sampling Plans for Preliminary Computational Experiments. Technometrics 33, 161–174 (1991).
  • Baldocchi D. et al. . FLUXNET: A New Tool to Study the Temporal and Spatial Variability of Ecosystem–Scale Carbon Dioxide, Water Vapor, and Energy Flux Densities. Bull. Am. Meteorol. Soc. 82, 2415–2434 (2001).
  • Haverd V., Smith B. & Trudinger C. Dryland vegetation response to wet episode, not inherent shift in sensitivity to rainfall, behind Australia’s role in 2011 global carbon sink anomaly. Global Change Biol. 22, 2315–2316 (2016). [PubMed]
  • Ma X. et al. . Spatial patterns and temporal dynamics in savanna vegetation phenology across the North Australian Tropical Transect. Remote Sens. Environ. 139, 97–115 (2013).
  • Nemani R. R. et al. . Climate-Driven Increases in Global Terrestrial Net Primary Production from 1982 to 1999. Science 300, 1560–1563 (2003). [PubMed]
  • Zhao L. et al. . Seasonal variations in carbon dioxide exchange in an alpine wetland meadow on the Qinghai-Tibetan Plateau. Biogeosci. 7, 1207–1221 (2010).
  • Zhang L.-X., Zhou D.-C., Fan J.-W. & Hu Z.-M. Comparison of four light use efficiency models for estimating terrestrial gross primary production. Ecol. Model. 300, 30–39 (2015).
  • Yuan W. et al. . Global comparison of light use efficiency models for simulating terrestrial vegetation gross primary production based on the LaThuile database. Agric. For. Meteorol. 192–193, 108–120 (2014).
  • Oak Ridge National Laboratory Distributed Active Archive Center (ORNL DAAC). FLUXNET Web Page. Available online [] from ORNL DAAC, 2015).
  • Sjöström M. et al. . Exploring the potential of MODIS EVI for modeling gross primary production across African ecosystems. Remote Sens. Environ. 115, 1081–1089 (2011).
  • Tan B. et al. . An Enhanced TIMESAT Algorithm for Estimating Vegetation Phenology Metrics From MODIS Data. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 4, 361–371 (2011).
  • Campolongo F., Cariboni J. & Saltelli A. An effective screening design for sensitivity analysis of large models. Environ. Model. Software 22, 1509–1518 (2007).
  • Pujol G., Iooss B. & A. J. Sensitivity: sensitivity analysis. R package 1.7 ed (2013).
  • Gelman A. et al. . Bayesian Data Analysis (CRC Press, 2013).
  • Metropolis N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H. & Teller E. Equation of state calculations by fast computing machines. J. Chem. Phys. 21 (1953).
  • Hastings W. K. Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57 (1970).
  • Taylor K. E. Summarizing multiple aspects of model performance in a single diagram. J. Geophys. Sci.: Atmosph. 106, 7183–7192 (2001).

Articles from Scientific Reports are provided here courtesy of Nature Publishing Group