Regional contribution to variability and trends of global gross primary productivity

Terrestrial gross primary productivity (GPP) is the largest component of the global carbon cycle and a key process for understanding land ecosystems dynamics. In this study, we used GPP estimates from a combination of eight global biome models participating in the Inter-Sectoral Impact-Model Intercomparison Project phase 2a (ISIMIP2a), the Moderate Resolution Spectroradiometer (MODIS) GPP product, and a data-driven product (Model Tree Ensemble, MTE) to study the spatiotemporal variability of GPP at the regional and global levels. We found the 2000–2010 total global GPP estimated from the model ensemble to be 117 ± 13 Pg C yr−1 (mean ± 1 standard deviation), which was higher than MODIS (112 Pg C yr−1), and close to the MTE (120 Pg C yr−1). The spatial patterns of MODIS, MTE and ISIMIP2a GPP generally agree well, but their temporal trends are different, and the seasonality and inter-annual variability of GPP at the regional and global levels are not completely consistent. For the model ensemble, Tropical Latin America contributes the most to global GPP, Asian regions contribute the most to the global GPP trend, the Northern Hemisphere regions dominate the global GPP seasonal variations, and Oceania is likely the largest contributor to inter-annual variability of global GPP. However, we observed large uncertainties across the eight ISIMIP2a models, which are probably due to the differences in the formulation of underlying photosynthetic processes. The results of this study are useful in understanding the contributions of different regions to global GPP and its spatiotemporal variability, how the model- and observational-based GPP estimates differ from each other in time and space, and the relative strength of the eight models. Our results also highlight the models’ ability to capture the seasonality of GPP that are essential for understanding the inter-annual and seasonal variability of GPP as a major component of the carbon cycle.


Introduction
Terrestrial ecosystems play a critical role in the global carbon cycle (Le Quéré et al 2016). Gross primary productivity (GPP), the carbon uptake by terrestrial ecosystems through plant photosynthesis, is the largest global CO 2 flux (Le Quéré et al 2016) and the major driver of many ecosystem processes. Therefore, it is important to understand the spatiotemporal variability of GPP for obtaining reliable estimates of terrestrial ecosystems capacity to serve as a major reservoir for carbon, especially in light of continued buildup of atmospheric carbon dioxide in the atmosphere for the rest of this century .
The spatial and temporal variations of global GPP are controlled by climate conditions, vegetation types and their spatial distribution, the nutrient availability, and other factors such as land-use practice that affect the distribution and composition of ecosystems (Ahlström et al 2015). Modeling studies suggest large differences in different approaches used to obtain estimates of global GPP (Anav et al 2013), largely controlled by seasonal, decadal and longer time variability and change in different regions, globally (Ahlström et al 2015). Many studies have estimated GPP in key regions of the world: for example, Lee et al (2013) estimated GPP in Amazonia using Greenhouse gases Observing SATellite (GOSAT) measurements; Liu et al (2014) estimated GPP in China using five GPP models and Moderate Resolution Imaging Spectroradiometer (MODIS) observations; Jung et al (2008) estimated GPP in Europe (EU) using terrestrial ecosystem models; Nightingale et al (2008) and Chen et al (2011) estimated GPP in the USA using MODIS observations; and Friedlingstein et al (2010) estimated GPP in the Sahel region of Africa using ORCHIDEE model. In spite of this large body of research, there is limited information on the relative contributions of these regions to the seasonal, inter-annual and longer time variability of total global GPP.
One significant challenge is that GPP cannot be directly measured (Ma et al 2015). The most common methods for estimating regional and global GPP include using remote sensing data, statistical interpolation, and process-based carbon cycle model simulations. For example, the MODIS GPP estimate is derived from satellite-based observations and a lightuse efficiency algorithm (Running et al 2004); the Model Tree Ensemble GPP (MTE GPP) product is derived from eddy covariance flux measurements by the global FLUXNET network and the parameterization of the relationship between GPP and the explanatory variables (Jung et al 2011, Beer et al 2010. Generally, good spatiotemporal correlation between site-level inferred GPP and MODIS-based estimates as well as MTE products has been reported in literature (Sjöström et al 2013, Zhu et al 2016, Gebremichael and Barros 2006, Turner et al 2006, and these products have been widely used for analyzing spatiotemporal variation of global GPP, and for benchmarking the process-based model estimates of GPP. With common simulation protocols, and consistent input data, model intercomparison projects (MIPs) have been launched to analyze terrestrial ecosystem model outputs (e.g. GPP) in a systematic fashion, improve models' estimates, and facilitate model improvement. The recent Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) provides a framework for developing and using such a common modeling protocol by a number of models, and offer a unique opportunity to use their simulated results to better understand the spatial, inter-annual and seasonal variation of GPP at regional and global scales (Warszawski et al 2014). These models serve as an important and alternative way to estimate GPP variation for a range of spatial and temporal scales, and in response to multiple environmental factors (Ahlström et al 2015, Xia et al 2015. However, as indicated by previous MIPs, such as the Multi-scale synthesis and Terrestrial Model Intercomparison Project (MsTMIP) (Huntzinger et al 2013) and TRENDY (Sitch et al 2008), the magnitude and patterns of model-estimated GPP varies largely across the models, indicating limited understanding of the fundamental underlying ecological process.
In this study, we used a combination of GPP products from MODIS and MTE GPP, and those simulated by eight ISIMIP Phase 2a (ISIMIP2a) terrestrial biome models to understand the contribution of terrestrial ecosystems to carbon cycle under the historical conditions, for the 1971-2010 period. Our main objectives were to: (1) report the ISIMIP2a model-estimated GPP and evaluate the ISIMIP2a models' performance against MODIS and MTE-GPP; (2) examine the spatial, inter-annual and seasonal variability of GPP at global and regional scales; and (3) assess the relative contribution of major regions to the global GPP.

Data
The MODIS GPP (MOD17) (Zhao and Running 2010) is the first satellite-based modeled dataset for monitoring vegetation productivity at the global scale. A light use efficiency model is the core of MODIS GPP algorithm (Running et al 2004). The MODIS GPP data is available for 8 day, monthly and annual timescales at 1 km and 0.05 • × 0.05 • spatial resolution since 2000, and we analyzed the data for 2000-2010 period in this study. The MODIS GPP products were resampled to 0.5 • × 0.5 • resolution by area-weighted averaging of the values in each 10 × 10 pixel window from the original 0.05 • × 0.05 • product.
The MTE-GPP product (Jung et al 2009(Jung et al , 2011 was upscaled based on the globally distributed FLUXNET eddy-covariance tower measurements with a 'Model Tree Ensembles' machine-learning algorithm. The MTE-GPP is currently available as a monthly We used eight ISIMIP2a biome models: CARAIB, DLEM, JULES, LPJmL, LPJ-GUESS, ORICHIDEE, VEGAS, and VISIT from the ISIMIP Phase 2a project. A summary of underlying processes for these models is presented in table 1. We used the simulated GPP from these eight models for 1971-2010 period, all at 0.5 • × 0.5 • spatial resolution. We calculated the mean GPP from the eight models (hereafter 'ENSEMBLE'). The models' simulations were performed using the common ISIMIP2a protocol, including time variant CO 2 concentrations, climate forcing, and land use change data. All models used the same climate forcing data from Global Soil Wetness Project 3 (GSWP3). More details of ISIMIP2a models are in the supplementary material available at stacks.iop.org/ERL/12/105005/mmedia, and more information about the models' simulation protocol and input data are available at www.isimip.org/protocol/#isimip2a.

Analysis
We divided the global land area into 12 regions We calculated the temporal trend of GPP for the annual global and regional GPP (area-weighted) using linear least square regression method. F-statistic was used to test the significance of the linear trend. GPP anomaly was calculated as the departure of the longterm mean of the detrended annual-GPP. The GPP inter-annual variability (IAV) for each region and the entire globe was calculated as the standard deviation of the detrended annual GPP at each of these spatial scales. We compared the IAV of ISIMIP2a model estimated GPP against MODIS and MTE by calculating the Pearson's correlation coefficient (R), with the detrended annual GPP. For each region and the entire globe, we calculated the mean monthly GPP to analyze the seasonal variations, which are also compared across ISIMIP2a models, MODIS and MTE with the Pearson's correlation coefficient. We also calculated each region's relative contribution to global GPP trend, seasonality and IAV based on the method from Ahlström et al (2015) (see supplementary material for details). The Matlab R2016b was used for all statistical analysis.

Spatial pattern of global GPP
The spatial pattern of GPP from MODIS, MTE and ENSEMBLE generally agree well (figure 2). The highest GPP values are in tropical regions (e.g. Amazonia, Central Africa and Southeast Asia) because of the wet and humid climate; and, the lowest GPP values are mainly in the arid regions of North Africa, West Asia and South  America, as well as the cold and snow-covered regions (e.g. Greenland), due to the unfavorable conditions for plant photosynthesis. However, the representation of this variability by eight models was different, with JULES and DLEM producing exceptionally high GPP in tropical regions, especially in Amazonia and the Sahel, while CARAIB producing high GPP in the temperate regions.
Over the 2001-2010 period, the eight models' estimates diverge in total global GPP, ranging from 106 (ORCHIDEE) to 134 (VISIT) Pg C yr −1 (figure 3(a), figure S1a), with a standard deviation of 13 Pg C yr −1 . The ENSEMBLE estimate of total global GPP is 118 Pg C yr −1 , which is higher than MODIS estimate of 112 Pg C yr −1 but lower than MTE estimate of 120 Pg C yr −1 (figure 3(a)). Based on the ENSEM-BLE estimates (figure 3(b)), the regional sum of GPP in Tropical Latin America is the highest among all 12 regions (23 Pg C yr −1 ), followed by Southern Africa (17 Pg C yr −1 ) and Temperate Asia (13 Pg C yr −1 ). The top three regions together contribute about 45% to global GPP (figure 3(c)). However, GPP estimates per unit area are the highest in tropical regions (Tropical Africa, Tropical Asia and Tropical South America) with > 2000 g C m −2 yr −1 (figure 3(d)). The regional contributions to global GPP from MODIS and MTE show similar patterns to the ISIMIP2a models based estimates.

Trend and regional contributions
The   5(a)). Inter-annual changes of global and regional GPP. The anomalies were calculated as the difference between annual GPP and the long-term mean over the period for each data source. The numbers located at the top-left corners indicate the linear trends of the ENSEMBLE (black), MTE (blue) and MODIS (red). * indicate the trend is significant (p < 0.01). Figure 6 illustrates the mean seasonal cycle of global and regional GPP. At the global scale, all datasets show similar seasonal cycles. For MODIS and ENSEMBLE, global GPP starts from as low as about 50 g C m −2 month −1 in the Northern Hemisphere (NH) winter months (i.e. December, January and February) to high values of about 100 g C m −2 month −1 in the NH summer months (i.e. June, July and August). The amplitude of MTE GPP is about 10 g C month −1 higher. The NH regions generally have similar seasonal patterns but higher amplitudes (as large as ∼150 g C m −2 month −1 ) than the global results, but the amplitudes of GPP seasonality in Semi-arid and Arid Asia-Europe-Africa are much lower than the other NH regions. There are no strong seasonal variation in tropical regions such as Tropical Latin America, Southern Africa and Tropical Asia. The Southern-Hemisphere (SH) regions of Temperate South America and Oceania show opposite patterns, which offset the higher seasonal amplitude of the NH regions resulting in the lower amplitude in seasonality of global GPP. We found a strong correlation (R > 0.9) between the seasonality of ISIMIP2a models, MODIS and MTE GPP at the global scale (figures S2(a) and (b)). This varies at the regional level, however, the correlation coefficients are high in boreal and temperate regions, such as Boreal North America, Temperate North America, Temperate South America, Europe, Boreal Euro-Asia and Temperate Asia, but low in the other semi-arid, arid and tropical regions. The correlation is especially low in the Tropical Latin America, Tropical Africa and South Africa regions. In addition, the models do not agree completely on the seasonal phase of GPP in reference to MODIS and MTE products. Based on the ENSEMBLE estimates, each of the regions in the NH, including Boreal North America, Temperate North America, Europe, Boreal Euro-Asia and Temperate Asia, contributes to the global GPP seasonal variation (figure 5(b)) ranging from 17% to 29%. Southern Hemisphere regions (Temperate South America, South Africa and Oceania) together contribute about −14% to global GPP variability. The Semi-arid and Arid Asia-Europe-Africa and tropical regions have small contributions (−1% to 1%) to seasonal variations of global GPP. The relative regional contributions calculated from MODIS, MTE and each of the eight models support similar conclusions.

Inter-annual variability and regional contributions
During the 2000-2010 period, MODIS data suggest that the interannual variability (IAV) of global GPP was 0.8 Pg C yr −1 , while MTE and ENSEMBLE estimates are nearly twice this amount (1.42 and 1.48 Pg C yr −1 , respectively; figure S1(d)). Among the regions, the ENSEMBLE generally shows the highest IAV in Boreal North America, Tropical Latin America, Semi-arid and Arid Asia-Europe-Africa, Tropical Africa, South Africa, Boreal Euro-Asia and Oceania; MODIS shows relatively higher IAV than MTE in all the regions except in Europe and Semi-arid and Arid Asia-Europe-Africa. The IAV estimates vary among ISIMIP2a models. For example, the IAV of GPP at the global scale estimated from JULES, LPJmL and VISIT are larger than 2 Pg C yr −1 , significantly higher than the 0.8 Pg C yr −1 from DLEM. At the regional scale, for example, according to LPJ_GUESS, LPJmL and VISIT the GPP in Oceania varies from 0.94 to 1.01 Pg C yr −1 , while it is 0.37 and 0.26 Pg C yr −1 , respectively, according to DLEM and JULES (figure S1(d)).
The correlation between the interannual variation of ISIMIP2a models, MODIS and MTE GPP at the global scale are not as good as that of the seasonality (figures S2(c) and (d)). The correlation coefficients were 0.41 and 0.57 for ENSEMBLE vs. MODIS and ENSEM-BLE vs. MTE, respectively. However, the correlation coefficients are relatively high in boreal and temperate regions such as Boreal North America, Temperate North America, Semi-arid and Arid Asia-Europe-Africa, South Africa, Temperate Asia, and exceptionally high in Oceania. In contrast, the correlations are weak and even negative in some tropical regions, for example, ISIMIP2a models estimated interannual variation is  1982andISIMIP models: 1971-2010).
negatively correlated with MODIS and MTE estimates in Tropical Africa and Tropical Africa, respectively. Figure 5(c) shows the ENSEMBLE estimates of GPP for Oceania region that explains the largest fraction (39%) of IAV for global GPP, followed by South Africa (32%), Tropical Asia (10%) and Tropical Latin America (9%) regions. In contrast, the contributions of Temperate North America, Temperate South America and Semi-arid and Arid Asia-Europe-Africa regions to IAV of global GPP are negative (−2% to −4%). The remaining regions contribute relatively smaller and positive fractions to IAV of the global GPP, with Boreal Euro-Asia contributing the least (<1%). The relative regional contributions to the IAV of global GPP estimated from MODIS and MTE do not completely agree with ENSEMBLE mean values, nor with each other (figure 5(c)). For example, South America contributes the most to the IAV of global GPP according to MODIS and MTE data (46% and 22%, respectively); and MODIS data indicates the contributions from Tropical Latin America, Europe and Boreal Euro-Asia are negative, while all regions except Boreal Euro-Asia contribute positively based on the MTE data. There are also significant uncertainties among the estimated relative regional contributions from the eight ISIMIP2a models, especially in Tropical Latin America, South Africa and Oceania regions.

Discussion
In this study, we examined the spatial pattern, trends, and inter-annual and seasonal changes of global and regional GPP as simulated by the eight ISIMIP2a terrestrial biome models, and benchmarked them against observation-based MODIS and MTE GPP products.
The spatial distribution and patterns of GPP based on the ISIMIP2a models ensemble mean, MODIS, and MTE all generally agree, and are consistent with other independent global GPP estimates reported in literature (Yuan et al 2010, Chen and Zhuang 2014, Anav et al 2015. The ISIMIP2a model ensemble mean is close to the mean annual GPP from MODIS and MTE, and consistent with GPP reported by Jiang and Ryu (2016) based on a mechanistic model and eddy flux data (122 ± 2.5 Pg yr −1 for 2001-2011). However, it is lower than the 130-169 Pg C yr −1 reported by the Earth system models (ESMs) simulations from the Fifth Climate Model Intercomparison Project (CMIP5) (Anav et al 2015); the 150-175 Pg C yr −1 suggested by atmospheric isotope measurements (Welp et al 2011); and the 146 Pg C yr −1 by using atmospheric CO 2 observations in a carbon cycle data assimilation system (Koffi et al 2012). These differences are probably due to a variety of factors such as the errors of the climate forcing in the coupled ESMs, lack of consideration of plant photorespiration in the isotope-based models, and the uncertainty of assimilation method used in the data assimilation systems (Anav et al 2015).
The temporal variation of global GPP based on ISIMIP2a models' estimates did not completely agree with MODIS and MTE estimates. The ISIMIP2a models show significant trends in the global and regional GPP in almost every region (except Oceania), while MODIS and MTE have few trends over their reported periods, although it is interesting to note that the MTE trend of 0.09 Pg C yr −2 is very close to the 1989-2008 soil respiration trend of 0.1 Pg C yr −2 estimated through a global synthesis (Bond-Lamberty and Thomson 2010). The overall lack of a strong GPP trend, however, is most likely due to the fact that neither MODIS nor MTE explicitly take account of CO 2 fertilization effect in their algorithms (De Kauwe et al 2016), despite the observed increase of atmospheric CO 2 concentration by about 1.7 ppm yr −1 (Conway et al 1994). Therefore, MODIS and MTE data do not appear to be suitable for trend analysis, as also indicated earlier by Jung et al (2009) andAnav et al (2015). In contrast, the increasing trends in model-estimated GPP appear to be consistent with other independent GPP indicators such as vegetation indices (de Jong et al 2012), and atmospheric carbonyl sulfide records (Campbell et al 2017).
Benchmarked against MODIS and MTE data, ISIMIP2a models well captured the global GPP seasonality, but they are relatively weaker in simulating the GPP seasonal cycle in the tropical regions primarily covered by tropical evergreen forests. In fact, most existing terrestrial biome models perform poorly in simulating tropical evergreen forest phenology, a key biophysical control of GPP seasonality (Kim et al 2012, Restrepo-Coupe et al 2017. Other recent studies also suggest that incorporating both the seasonal change of leaf quantity (i.e. leaf area) and quality (i.e. leaf photosynthetic capacity) can successfully explain seasonal GPP variation , Wu et al 2017. Including this new mechanism may help the models to better simulate tropical GPP seasonality.
The IAV of GPP from MODIS, MTE and the ISIMIP2a models are generally more consistent in the temperate than in the tropical regions (figures S2(c) and (d)). One possible explanation is that IAV from MODIS and MTE in tropical regions are not as reliable as their seasonality, because some of the year-to-year change of GPP sensitivity to climate are not captured in their algorithms (Piao et al 2013), and the satellitebased input data are in a relatively low quality due to dense and persistent clouds in these regions (Trenberth et al 2001). The uncertainty of climate forcing used in ISIMIP2a models are large due to less meteorological station observations in these regions than the rest of the world, thus the models could produce inaccurate GPP estimates in these regions (Chang et al 2017). Another possible reason could be the poor simulation of tropical phenology in the models (Pau et al 2011) as indicated above. Phenology controls the length of growing season, thus being a key determinant of annual GPP (Keenan et al 2012, Richardson et al 2013, Chen et al 2016. Many previous studies have suggested that accurate simulation of vegetation phenology is the key to capturing IAV of the carbon cycle (Keenan et al 2012. We urge the modeling community to improve the phenology algorithms in their models, and to evaluate how the improved ecosystems phenology, especially in tropical regions, can help improve the IAV of regional and global GPP. We found the regional contributions of terrestrial ecosystems to the global GPP trend, seasonality, and IAV to be generally consistent among independent methods and estimates (i.e. MODIS, MTE and ISIMIP2a models), and those reported in literature. Overall, all regions have increasing trends, and the NH regions contribute the most to the global increasing trend, mainly due to the lengthening of the growing seasons as a result of a warmer and wetter climate conditions that enhance plant growth (Piao et al 2007, Zhao andRunning 2010). Regions in the Southern Hemisphere have smaller increasing trends, and some of the models suggest negative trend in some regions (e.g. South Africa, Oceania). This is possibly due to the increasing evaporative demand, which leads to a drying trend in these regions (Zhao and Running 2010), and reducing soil moisture (Jung et al 2010) and vegetation greenness (Gobron et al 2010). The Northern Hemisphere dominates the world's land area, and has the largest forest ecosystems, therefore, it is not surprising that it contributes the most to the seasonal variability of global GPP; and, the SH regions generally contribute negatively to it. Some tropical regions, including Tropical Latin America and Tropical Africa, cover areas that extend across the equator, therefore their net contribution are close to zero.
Our results indicate that the Oceania and South Africa regions contribute the most to the IAV of global GPP, supported by all the data sources used in this study. The dominant land cover in these regions are semi-arid savannas, grasslands and shrublands that are very sensitive to climate conditions and its variability. These findings are consistent with those previously reported by (Ahlström et al 2015). The IAV of GPP for these vegetation types are especially sensitive to the change of precipitation. For example, due to the periodic El Niño and La Niña events, Australia has received record high precipitation in 1975, 2000 and 2010, and experienced exceptional droughts in 1972, 1994 and 2002 (according to the data from Australian Bureau of Meteorology, www.bom.gov.au/climate/change/). The influence of these extreme events/conditions are clearly reflected in the inter-annual variation of GPP in Oceania, and our finding are consistent with previous studies (Ahlström et al 2015, Zhang et al 2016, Zscheischler et al 2014. In addition, the ISIMIP2a models' estimates suggest that the tropical forests in Amazonia and Southeastern Asia (regions of Tropical Latin America and Tropical Asia in this study) are the second largest contributor to the IAV of global GPP, although MODIS and MTE data do not agree with this result. There remain some debates on satellite-observed tropical forests' response to extreme climate conditions (Morton et al 2014, Saleska et al 2007, thus we cannot conclude whether the models or MTE/MODIS data are more reliable, since largescale measurements of GPP do not exist, and terrestrial biome models show large uncertainties in these regions. The eight ISIMIP2a model-based GPP estimates used in this study, having common climate forcing, land use and CO 2 input data were expected to provide consistent GPP estimates, especially at inter-annual and seasonal levels. However, differences relating to the parametrization, model formulation and underlying assumptions that exist in these models appear to contribute to GPP variability, and the model ensemble approach is a good way to capture these uncertainties. It is not always clear how these assumptions propagate in the simulated parameters, especially if multiple processes are at play (Rafique et al 2015). Increasing complexity of models may lead to an increase in models' uncertainty associated with the introduction of new model parameters (Prentice et al 2015). In this study, models' simulations did not reveal a consistent set of processes behind the agreement/disagreement with MODIS and MTE GPP. Apparently, formulation of photosynthesis based on enzyme kinetics and light use efficiency (LUE) in the eight models performed equally well in simulating the inter-annual and seasonal variability of GPP. The addition of nitrogen and leaf to canopy details did not show significant influence on the GPP simulations. It is interesting that none of the models performed well in Tropical Latin America (Restrepo-Coupe et al 2017), suggesting that the phenological sub-models based on temperature and moisture may need further improvement. A lack of particular pattern in the models' performance does not mean that structural differences have no effect, however. For example, better parameter values of LUE and leaf to plant canopy scaling can largely improve the GPP estimates (Schaefer et al 2012). The number of parameters also differ widely among the eight models. Despite the differences among the eight models, the results highlight their ability to capture the importance of relative contributions of diverse ecosystems of the regions in determining the inter-annual and seasonal variability of global GPP.

Conclusion
The ISIMIP2a was established to foster model evaluation, and in the global biome models, how well they represent the role of terrestrial ecosystems in the carbon cycle under historical climate conditions. We used simulations from eight global terrestrial biome models participating in the ISIMIP2a to examine the spatial and temporal variability, the changes in twelve geographical regions, globally, and their contribution to the temporal variability and change in global and regional GPP. The simulated GPP from ISIMIP2a models, driven by common climate forcing, land use and CO 2 data, were 117 ± 13 Pg C yr −1 (model ensemble mean ± 1 standard deviation), which is in close agreement with independent global MODIS and MTE based GPP estimates. The model ensemble generally showed closer interannual variability for the higher latitudes (i.e. boreal and temperate regions) than the lower latitudes (i.e. tropical regions), while at seasonal scale, the models performed very well except in South Africa and tropical forests in Tropical Latin America and Tropical Africa regions, benchmarked against GPP estimates from MODIS and MTE. According to the model ensemble mean estimates, Tropical Latin America contributes the most to the global mean annual GPP; the regions in Asia contribute the most to the global GPP trend; the Northern Hemisphere regions dominate the global GPP seasonal variation; and Oceania is likely the largest contributor to the global GPP inter-annual variability. The results of this study provide some useful insights on the ability of the terrestrial vegetation models in capturing the spatiotemporal variability and change, and the relative contribution of twelve regions to global GPP, and its variability. These results can be used for model development/improvement, and predicting future changes in regional and global terrestrial ecosystems conditions and the models' performance as relate to the global carbon cycle, in a changing climate.