Regional contribution to variability and trends of global gross primary productivity

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₂ 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 (Ciais et al 2013).

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 light-use 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.

2.1. Data

2.2. Analysis

We divided the global land area into 12 regions (figure 1) by considering their climate, land cover and geopolitical characteristics: Boreal North America (BNA), Temperate North America (TNA), Tropical Latin America (TLA), Temperate South America (TSA), Europe (EUR), Semi-arid and Arid Asia-Europe-Africa (SAAEA), Tropical Africa (TAF), South Africa (SA), Boreal Euro-Asia (BEA), Temperate Asia (TEA), Tropical Asia (TA) and Oceania (OCE). 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 long-term 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1. Spatial pattern of global GPP

3.2. Trend and regional contributions

The ENSEMBLE global GPP increased from 105 Pg C yr⁻¹ in 1971 to 118 Pg C yr⁻¹ in 2010, with a significant trend of 0.28 Pg C yr⁻² (p < 0.01, figure 4). During 1971–2010, all regions except Oceania show significant increasing trends in GPP, with the fastest rates of increase in Boreal Euro-Asia and Temperate Asia (both 0.04 Pg C yr^-2). The trend of the ensemble mean is similar to the trends of individual GPP from the eight ISIMIP2a models, although their magnitudes are different (figure S1(b)). For example, some models (e.g. CARAIB, LPJ_GUESS) have larger positive trend of global GPP than the others (e.g. VEGAS), (figure S1(b)). In contrast to the ISIMIP2a models, MTE and MODIS estimates of global GPP show smaller long-term change over their record periods. During 1982–2010, MTE global GPP increased at a small rate of 0.08 Pg C yr⁻², and only Europe, Semi-arid and Arid Asia-Europe-Africa and Temperate Asia regions show a weak but significant increasing trend (figure 4). In general, the MODIS product does not show a significant trend in any region during 2000–2010, except in Tropical Asia and Tropical Africa, which have exhibit a significantly decreasing trend of −0.09 Pg C yr⁻² and a significantly increasing trend of 0.04 Pg C yr⁻², respectively. In contrast, ENSEMBLE GPP during 2000–2010 showed a significant large increasing trend of 0.41 Pg C yr⁻², and only Temperate South America and Oceania had decreasing trends. MTE GPP during 2000–2010 generally showed similar but weaker trends as compared with ENSEMBLE GPP (figure S1(c)). Based on the ENSEMBLE GPP, Boreal Euro-Asia and Temperate Eastern Asia contribute most to the global GPP trend (each of them contributes about 15%), while Semi-arid and Arid Asia-Europe-Africa and Oceania contribute the least, only about 2% by each region (figure 5(a)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.3. Mean seasonality and regional contributions

3.4. Inter-annual variability and regional contributions

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⁻¹ for 2001–2011). However, it is lower than the 130–169 Pg C yr⁻¹ 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⁻¹ suggested by atmospheric isotope measurements (Welp et al 2011); and the 146 Pg C yr⁻¹ by using atmospheric CO₂ 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⁻² is very close to the 1989–2008 soil respiration trend of 0.1 Pg C yr⁻² 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₂ fertilization effect in their algorithms (De Kauwe et al 2016), despite the observed increase of atmospheric CO₂ concentration by about 1.7 ppm yr⁻¹ (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) and Anav 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 2016, 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 satellite-based 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, Richardson 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 and Running 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, Zhou et al 2014, Saleska et al 2007, 2016), thus we cannot conclude whether the models or MTE/MODIS data are more reliable, since large-scale 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₂ 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.

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₂ data, were 117 ± 13 Pg C yr⁻¹ (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.

This study was funded primarily by a Laboratory Directed Research and Development project sponsored by the Pacific Northwest National Laboratory managed by Battelle Memorial Institute for the US Department of Energy. The study was also partly supported by the National Aeronautics and Space Administration Carbon Monitoring System (CMS) and ACCESS programs under projects NNH12AU35I and NNH13AW58I. The ISIMIP project is grateful for funding from the German Federal Ministry of Education and Research (BMBF, grant no. 01LS1201A1). Catherine Morfopoulos, Richard Betts and Jinfeng Chang received support from the European Commission's 7th Framework Programme (EU/FP7) under Grant Agreement 603864 (HELIX). Tian, Pan, and Yang acknowledge funding support from US National Science Foundation (1210360, 1243232), National Key Research and Development Program of China (no. 2017YFA0604700).

Ben Bond-Lamberty https://orcid.org/0000-0001-9525-4633