Vegetation structural change since 1981 significantly enhanced the terrestrial carbon sink
LAI data analysisA global LAI time series from 1981 to 2016 at 8 km resolution and 16-day (1981–2000) and 8-day (2001–2016) intervals was produced using Advanced Very High Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data with the GLOBMAP algorithm19. The global distribution of the temporal trend of LAI over this period is shown in Fig. 1. About 74.2% of the land surface shows an increasing trend, among which 52.7% is significant at p = 0.05 level (two-tailed). Globally, the increase in the maximum LAI in the peak growing season is about twice as large as that in the annual mean, suggesting that growing season lengthening is not the main factor explaining the greening trend. This general increase in LAI resulted from the combined effects of various drivers including CO2, climate, and nitrogen deposition over the same period, and therefore provides a new base for separating the effects of these drivers on vegetation structure and growth. The focus of this study is on the increase of LAI on plant growth and the land carbon sink. Considering the uncertainty in the temporal trend of this LAI time series, we used two other LAI products: GLASS LAI20 and GIMMS LAI3g21 (Supplementary Fig. 1 and Table 1).
Analysis of the residual land sinkDriven by one of LAI time series at a time as well as climate, CO2, soil, and nitrogen deposition data, BEPS is used to simulate GPP, AR, and HR at daily time intervals for each pixel. The sum of simulated net ecosystem productivity (NEP), taken as GPP-AR-HR, for all land areas is compared with the global residual land sink (RLS) reported by the GCP4, which is computed as the sum of emissions from fossil fuel consumption, cement production, and land-use change minus the sum of CO2 accumulated each year in the atmosphere and ocean. The modeled annual NEP as average from the BEPS model forced by each of the three LAI products closely follows the trend and interannual variability of the residual land carbon sink (Fig. 2), although it does not capture well extremely low and high values in some years, such as 1987, 1991, 2000, 2002, and 2009. Over the 1981–2016 period, the modeled accumulated NEP is 95 ± 5 Pg C, whereas the accumulated RLS is 94 ± 5 Pg C.
Attribution of the land sinkThe increase in atmospheric CO2 concentration is modeled to be the dominant diver for the land sink enhancement during the 1981–2016 period. The CO2 concentration increase after 1981 alone enhanced the global land sink by 44.6 Pg C accounting for 47.0% of the total accumulated sink enhancement after 1981 (Fig. 3). The simulated global total net primary productivity (NPP) increased by 11.6% with the increase of CO2 concentration from 340.13 p.p.m. in 1981 to 404.20 p.p.m. in 2016, whereas climate, LAI, and N deposition remain at baseline values. The sensitivity of simulated total NPP in the northern hemisphere to atmospheric CO2 concentration (β-factor3) was 18.6%/100 p.p.m. during 1981–2016. This β-factor is about 6.1% smaller than the average value of prognostic models participating in the Coupled Model Intercomparison Project 5 (CMIP5)3. Our β-factor excludes the effects of CO2 fertilization on LAI, whereas CMIP5 models lump vegetation LAI responses to CO2 increase into the β-factor. If the influence of LAI change is included in our β-factor calculation, it would be 31.8%/100 p.p.m. for the northern hemisphere, whereas the global average value is 27.0%, because LAI increased much less in the southern hemisphere than in the northern hemisphere. Both BEPS and CMIP5 models include the influences of CO2 increase on stomatal conductance affecting photosynthesis, transpiration, and soil moisture. Our β-factor is considerably larger than the values inferred from remote-sensing NPP models3, because these models are mostly based on empirical LUEs that do not explicitly consider enhanced LUE at higher CO2 levels24 (see Supplementary Fig. 4 for further comparisons). CO2 fertilization enhanced the carbon sink in all regions (Fig. 4), especially in regions with high NPP (15°S–30°N and 45°N–60°N) (Supplementary Fig. 5).
GPP modeling methodsThe BEPS model17,18 used in this study is a process-based diagnostic model driven by remotely sensed vegetation parameters, including LAI, clumping index, and land cover type, as well as meteorological and soil data. It simulates photosynthesis, energy balance, and hydrological and soil biogeochemical processes at daily time steps17,28. For GPP simulation, BEPS uses the leaf-level biochemical model29 with a two-leaf upscaling scheme from leaf to canopy:17
Full carbon cycle modeling methodsBEPS includes modules to calculate HR and NEP28. Based on a modified Century model38,39, it stratifies the biomass carbon stock into four pools (leaf, stem, coarse root, and fine root pools) and the soil carbon stock into nine pools (surface structural litter, surface metabolic litter, soil structural litter, soil metabolic litter, coarse woody litter, surface microbe, soil microbe, slow, and passive carbon pools). These carbon pools are initialized in the following way. The model first calculates NPP in 1901 using N deposition and CO2 concentration in 1901, a random year of climate data selected in the period from 1901 to 1910, seasonally variable LAI averaged over the 1982–1986 period, and default C:N ratios for all carbon pools. The nine soil and four biomass carbon pools are then estimated under the assumption that the carbon cycle of terrestrial ecosystems was in dynamic equilibrium in 1901. With this assumption, all carbon pools are determined by solving a set of equations describing the dynamics of carbon pools40. For the period from 1901 to 1980, the model is run using historical data of N deposition, CO2 concentration, and climate, and the average LAI during 1982 to 1986. Due to lack of data, we assume that LAI in 1982–1986 represents that in 1901–1981. If LAI increased in this period, NPP in 1901–1910 would be overestimated, leading to larger soil carbon pools in 1901 and smaller NEP in 1981–2016. To address this issue, we conducted a set of simulations by extending the LAI time series to 1901 according to atmospheric CO2 concentration with the rate of LAI change with CO2 determined using 1981–2016 data. We find that the enhancement of the land sink due to LAI change in 1981–2016 decreases from 12.4% to 10.5% relative to the accumulated sink in the same period when the possible LAI increase from 1901 to 1981 is considered. This decrease is due to higher accumulated NEP by 6.8 Pg C during 1981–2016, resulting from lower initial soil carbon pools at 1901 when LAI is smaller than our simulations without considering LAI change over 1901–1981. This set of simulations suggests that the impact of possible LAI changes prior to 1981 on the role of LAI after 1981 is within a few percent and does not affect our conclusion on the significance of LAI increase after 1981 in enhancing the land sink. Although extrapolating LAI according to CO2 concentration is overly simplistic, it may be considered as setting the upper bound of the possible error due to LAI changes prior to 1981, because climate change could have been negative on plant growth and LAI. The decomposition of soil carbon and mineralization of soil nitrogen are regulated by soil temperature, moisture, texture, and chemical property of soil carbon pools. Nitrogen available for vegetation growth consists of the total of mineralized and deposited nitrogen. The uptake of nitrogen by vegetation is simulated according to temperature, total amounts of soil carbon and nitrogen, and vegetation demand. The absorbed nitrogen is allocated daily to different vegetation carbon pools based on the C:N ratios and NPP allocation coefficients. The nitrogen content of leaves is used to adjust the parameter Vcmax at 25 °C, which consequently affects the photosynthesis rates of sunlit and shade leaves in the Farqhuar model29. AR consists of maintenance and growth respiration, and maintenance respiration depends on foliage, stem and root biomass, and temperature, whereas growth respiration is taken as a fraction (25%) of GPP. NEP of each modeling grid equals GPP-AR-HR.
LAI dataLAI is an input into the BEPS model for the simulation of the carbon flux. Three LAI time series, GLOBMAP-V2, GLASS, and LAI3g are used in this study and are shown in Supplementary Fig. 1 in comparison with other LAI time series. The GLOBMAP _V2 product is the basis for our simulations, whereas GLASS and LAI3g are used to assess uncertainties in carbon sink estimation due to the choice of LAI products. GLOBMAP_v2 over the period from 1981 to 2016 was generated through fusing LAI inverted from MODIS reflectance data with AVHRR GIMMS NDVI data. LAI from 2001 to 2016 was first derived from the MOD09A1 C6 land surface reflectance and the associated illumination and view angles based on the GLOBCARBON LAI algorithm19,41, which was developed on the basis of the 4-Scale geometric optical model42. This algorithm explicitly considers the effects of the bidirectional reflectance distribution function on reflectance over the canopy as measured by the sensors41. For the fusion of MODIS and AVHRR remote-sensing data, the relationships between GIMMS NDVI and MODIS LAI were established pixel by pixel over a period (2001–2006)19 that they overlap. Then the AVHRR LAI from 1981 to 2000 was generated using these relationships, to ensure the temporal consistency between these two sensors. The spatial resolution of the LAI series is 0.072727° × 0.072727° and temporal resolution varies from 16 days (1981 to 2000) to 8 days (2001 to 2016). In the simulation, these 16 days and 8 days LAI values were interpolated into daily values. GLASS and LAI3g have the similar temporal coverage. All three long-term LAI series used in this study have similar magnitudes, because they have all considered the three-dimensional canopy structure, as characterized by the clumping index, in their retrieval algorithms43. Both GLOBMAP_V2 and GLASS used the same global clumping index map30, whereas LAI3g considered clumping in a different way21. For accurate simulation of sunlit and shaded leaf area and GPP, both LAI and clumping index are needed. All these three products used the processed AVHRR data (GIMMS). The issues with possible artifacts and errors in the AVHRR data series (GIMMS), such as sensor degradation, sensor intercalibration, orbital drift causing changes in sun-target-view geometry, distortions by clouds, and abnormal aerosol absorption by two major volcanic eruptions, have been fully considered and rectified to a large extent, ensuring the useful signals in the trend being extracted44. The GIMMS time series has quality flags with values from 0 to 6. GLOBMAP used only the top quality 0 or 1. Depending on the strength of quality control, different LAI products could show different interannual variations, although they are consistent in their increasing trends (Supplementary Table 1). LAI data have been assimilated into diagnostic ecosystem models to optimize plant carbon allocation, stock, and residence time, as well as carbon-use efficiency16,45. In our study, we similarly optimized some of these parameters through a spin-up procedure and also used LAI data to calculate the long-term global carbon cycle.
Meteorological dataMeteorological data required to force the BEPS model include daily maximum and minimum temperatures, downward solar radiation, relative humidity, and precipitation. These data are interpolated from the 0.5° × 0.5° CRUNCEP V8.0 dataset, which is a combination of CRU monthly climatology and 6 hourly NCEP reanalysis meteorological data46. Daily maximum and minimum temperatures, downward solar radiation, and precipitation are directly retrieved from the CRUNCEP V8.0 dataset. Relative humidity is calculated from temperatures, specific humidity, and pressure from the CRUNCEP V8.0 dataset. The 0.5° × 0.5° meteorological data are interpolated into 0.072727° × 0.072727° resolution using a bilinear interpolation method. These data are the same as those used by the TRENDY models.
Soil dataFractions of clay, silt, and sand are retrieved from the harmonized global soil database (http://www.fao.org/nr/lman/abst/lman_080701_en.htm) and are used to determine soil physical parameters, including wilting point, field capacity, porosity, hydrological conductance, exponent of the moisture release equation, and so on.
Nitrogen deposition dataThe yearly global nitrogen deposition data at 0.5° × 0.5° resolution over the period from 1960 to 2009 are estimated from tropospheric NO2 column density retrieved from Global Ozone Monitoring Experiment and Scanning Imaging Absorption Spectrometer for Atmospheric Cartography sensors, meteorological data, and NOx emission inventory data47. For the years from 2010 to 2016, nitrogen data are extrapolated using the estimated nitrogen data over the period from 2000 to 2009. For the period from 1901 to 1959, nitrogen data are extrapolated based on the change rates of nitrogen deposition over the period from 1960 to 1969. The 0.5° × 0.5° nitrogen deposition data are interpolated into the 0.072727° × 0.072727° resolution using a bilinear interpolation method. The N deposition dataset used in this study is compared with that used by TRENDY models (Supplementary Fig. 2). These two datasets are similar before 1990, as both are based on measurements, but the increasing trends after 1990 are different, because the data used by the models are based on linear extrapolation from 1990 to 2050, at which the nitrogen deposition is estimated based on projected anthropogenic sources and other assumptions48, whereas satellite measurements from 2000 to 2009 are used in our dataset47 and could follow the realistic trend more closely than the linearly extrapolated trend used by the models. Over the 1981–2016 period, the total N deposition is 301 Tg N in our study, whereas it is 403 Tg N in TRENDY. The difference could be due to the overall decrease in N deposition in North America and other regions in this period.
Uncertainty assessmentThe uncertainty of model results shown in Figs. 2 and 3 is estimated based on differences among simulated results using the three LAI products. This uncertainty is of a similar magnitude to those estimated from the uncertainties of model parameters that influence the simulated temporal trends, because bias errors are mostly constrained by the atmospheric CO2 concentration and the main interest of this study is to attribute the land sink to the various drivers through their influences on the temporal trend. The uncertainty for the CO2 fertilization effect, e.g., is mostly caused by the uncertainties in the slope and intercept of the linear relationship between stomatal conductance and photosynthesis rate. The trend of NEP against climate is mostly controlled by the sensitivities of AR and HR to temperature.
Fiedlingstein, P. et al. Uncertainties in CMIP5 climate projections due to carbon cycle feedbacks. J. Clim. 27, 511–526 (2014).
Langley, J. A. & Hungate B. A. Plant community feedbacks and long-term ecosystem responses to multi-factored global change. AoB Plants 6, https://doi.org/10.1093/aobpla/plu035 (2014).
Smith, W. K. et al. Large divergence of satellite and Earth system model estimates of global terrestrial CO2 fertilization. Nat. Clim. Change 6, 306–310 (2016).
Le Quere, C. Global carbon budget 2017. Earth Syst. Sci. Data 10, 405–448 (2018).
Huntzinger D. N. et al. Uncertainty in the response of terrestrial carbon sink to environmental drivers undermines carbon-climate feedback predictions. Sci. Rep. 7, 4765, https://doi.org/10.1038/s41598-017-03818-2 (2017).
Piao, S. et al. Evaluation of terrestrial carbon cycle models for their response to climate variability and to CO2 trends. Glob. Change Biol. 14, 2015–2039 (2013).
Weider, W. R., Cleveland, C. C., Smith, W. K. & Todd-Brown, K. Future productivity and carbon storage limited by terrestrial nutrient availability. Nat. Geosci. 8, 441–444 (2015).
Clark, P. U. et al. Consequences of twenty-first century policy for multi-millennial climate and sea-level change. Nat. Clim. Change 6, 360–369 (2016).
Schimel, D. et al. Observing terrestrial ecosystems and the carbon cycle from space. Glob. Change Biol. 21, 1762–1776 (2015).
Bayer, A. D., Pugh, T. A. M., Krause, A. & Arneth, A. Historical and future quantification of terrestrial carbon sequestration from a greenhouse-gas-value perspective. Glob. Environ. Change 32, 153–164 (2015).
Arora, V. K. et al. Carbon concentration and carbon climate feedbacks in CMIP5 Earth System Models. J. Clim. 26, 5289–5314 (2013).
Thomas, R. Q., Canham, C. D., Weathers, K. C. & Goodale, C. L. Increased tree carbon storage in response to nitrogen deposition in the US. Nat. Geosci. 3, 13–17 (2009).
Chen, J. M. & Black, T. A. Defining leaf area index for non-flat leaves. Plant Cell Environ. 15, 421–429 (1992).
Zhu, Z. et al. Greening of the Earth and its drivers. Nat. Clim. Change 6, 791–795 (2016).
Jiang, C. et al. Inconsistencies of Interannual Variability and trends in long-term satellite leaf area index products. Glob. Change Biol. 23, 4133–4146 (2017).
Exbrayat, J. F. et al. Understanding the land carbon cycle with space data: current status and prospects. Surv. Geophys. 40, 735–755 (2019).
Chen, J. M., Liu, J., Cihlar, J. & Guolden, M. L. Daily canopy photosynthesis model through temporal and spatial scaling for remote sensing applications. Ecol. Model. 124, 99–119 (1999).
Chen, J. M. et al. Effects of foliage clumping on global terrestrial gross primary productivity. Glob. Biogeochem. Cycles 26, GB1019 (2012).
Liu, Y., Liu, R. & Chen, J. M. Retrospective retrieval of long-term consistent global leaf area index (1981-2010) maps from combined AVHRR and MODIS data. J. Geophys. Res. Biogeosci. 117, G04003 (2012).
Xiao, Z. et al. Long-time-series global land surface satellite leaf area index product derived from MODIS and AVHRR surface reflectance. IEEE Trans. Geosci. Remote Sens. 54, 5301–5318 (2016).
Zhu, Z. et al. Global data sets of vegetation leaf area index (LAI) 3 g and Fraction of Photosynthetically Active Radiation (FPAR) 3g derived from Global Inventory Modeling and Mapping Studies (GIMMS) Normalized Difference Vegetation Index (NDVI3 g) for the period 1981-2011. Remote Sens. 5, 927–948 (2013).
Ryu, Y. et al. Integration of MODIS land and atmosphere products with a coupled-process model to estimate gross primary productivity and evapotranspiration from 1 km to global scales. Glob. Biogeochem. Cycles 25, GB4017 (2011).
Potter, C. S. et al. Terrestrial ecosystem production: a process model based on global satellite and surface data. Glob. Biogeochem. Cycles 7, 811–841 (1993).
De Kauwe, M. G., Keenan, T. F., Medlyn, B. E., Prentice, I. C. & Terrer, C. Satellite based estimates underestimate the effect of CO2 fertilization on net primary productivity. Nat. Clim. Change 6, 892–893 (2016).
Zaehle, S. Terrestrial nitrogen-carbon cycles interactions at the global scale. Philos. T. Roy. Soc. B 368, https://doi.org/10.1098/rstb.2013.0125 (2013).
Churkina, G. et al. Synergy of rising nitrogen deposition and atmospheric CO2 on land carbon uptake moderately offset global warming. Glob. Biogeochem. Cy. 23, GB4027 (2009).
Song, X. P. et al. Global land change from 1982 to 2016. Nature 560, 639–643 (2018).
Ju, W. M. et al. Modelling multi-year coupled carbon and water fluxes in a boreal aspen forest. Agr. For. Meteorol. 140, 136–151 (2006).
Farquhar, G. D., von Caemmerer, S. & Berry, J. A. A biochemical model of photosynthetic CO2 assimilation in leaves of C3 species. Planta 149, 78–90 (1980).
He, L., Chen, J. M., Pisek, J., Schaaf, C. B. & Strahler, A. H. Global clumping index map derived from the MODIS BRDF product. Remote Sens. Environ. 119, 118–130 (2012).
Ball, J. T. An analysis of stomatal conductance, Ph.D Thesis, Standford University (1988).
Feng, X. et al. Net primary productivity of China’s terrestrial ecosystems from a process model driven by remote sensing. J. Environ. Manag. 85, 563–573 (2007).
Liu, Y. et al. Evapotranspiration and water yield over China’s landmass from 2000 to 2010. Hydrol. Earth Syst. Sci. 17, 4957–4980 (2013).
Schwalm, C. R. et al. A model-data intercomparison of CO2 exchange across North America: results from the North American Carbon Program site synthesis. J. Geophys. Res. Biogeosci. 115, https://doi.org/10.1029/2009jg001229 (2010).
Sprintsin, M., Chen, J. M., Desai, A. & Gough, C. M. Evaluation of leaf-to-canopy upscaling methodologies against carbon flux data in North America. J. Geophys. Res. Biogeosci. 117, https://doi.org/10.1029/2010jg001407 (2012).
Wang, Q., Tenhunen, J., Falge, E., Bernhofer, C. & Granier, A. Simulation and scaling of temporal variation in gross primary production for coniferous and deciduous temperate forests. Glob. Change Biol. 10, 37–51 (2004).
Matsushita, B. & Tamura, M. Integrating remotely sensed data with an ecosystem model to estimate net primary productivity in East Asia. Remote Sens. Environ. 81, 58–66 (2002).
Ju, W. M., Chen, J. M., Harvey, D. & Wang, S. Future carbon balance of China’s forests under climate change and increasing CO2. J. Environ. Manag. 85, 538–562 (2007).
Parton, W. J. et al. Observations and modeling of biomass and soil organic-matter dynamics for the grassland biome worldwide. Glob. Biogeochem. Cy. 7, 785–809 (1993).
Chen, J. M. et al. Spatial distribution of carbon sources and sinks in Canada’s forests. Tellus B 55, 622–642 (2003).
Deng, F., Chen, J. M., Plummer, S., Chen, M. Z. & Pisek, J. Algorithm for global leaf area index retrieval using satellite imagery. IEEE Trons. Geosci. Remote Sens. 44, 2219–2229 (2006).
Chen, J. M. & Leblanc, S. G. A four-scale bidirectional reflectance model based on canopy architecture. IEEE Trans. Geosci. Remote Sens. 35, 1316–1337 (1997).
Chen, J. M. in Comprehensive Remote Sensing, Vol. 3 (ed. Liang, S.) Pages 53–77, ISBN 9780128032213, https://doi.org/10.1016/B978-0-12-409548-9.10540-82018 (Elsevier, Oxford, 2017).
Tucker, C. J. et al. An extended AVHRR 8-km NDVI dataset compatible with MODIS and SPOT vegetation NDVI data. Int. J. Remote Sens. 26, 4485–4498 (2005).
Bloom, A. A., Exbrayat, J. F., van der Velde, I. R., Feng, L. & Williams, M. The decadal state of the terrestrial carbon cycle: global retrievals of terrestrial carbon allocation, pools and residence times. Proc. Natl Acad. Sci. USA 113, 1285–1290 (2015).
Viovy, N. CRUNCEP data set, available at: ftp://nacp.ornl.gov/synthesis/2009/frescati/temp/land_use_change/original/readme.htm, last access (2016).
Lu, X. H. et al. Estimated global nitrogen deposition using NO2 column density. Int. J. Remote Sens. 34, 8893–8906 (2014).
Dentener, F. J. Global maps of atmospheric nitrogen deposition, 1860, 1993, and 2050. Data set. Available on-line [http://daac.ornl.gov/] from Oak Ridge National Laboratory Distributed Active Archive Center, Oak Ridge, Tennessee, USA. https://doi.org/10.3334/ORNLDAAC/830 (2006).