Earth system models underestimate carbon fixation by plants in the high latitudes
Large inter-model spread in greening sensitivityThe enhancement in NHL greenness throughout the observational period relates linearly to both increasing quantities, GDD0 and CO2 concentration, in general agreement between models and observations15,19,31. To avoid redundancy from co-linearity between the two driver variables, but retain their underlying time-trend and interannual variability (Supplementary Table 1), we use the dominant mode from a principal component analysis (PCA) of CO2 and GDD0 as the proxy driver (denoted ω) in subsequent analysis (Methods). Expressed in this compact form, greenness level (Fig. 2a) as well as greening sensitivity to ω (hereafter greening sensitivity, Fig. 2b) span a wide range across the multi-model ensemble. All models with low greenness levels (LAImax < 0.75 m2 m−2) tend to simulate low greening sensitivities (< 0.015 m2 m−2 LAImax per 1 unit ω), relative to observations. These models (NorESM1-ME, CESM1-BGC, and CanESM2) lack a representation of dynamic vegetation, i.e., do not allow plant functional type shifts in response to changing simulated climate, and/or show overly strong nitrogen limitations on plant growth and thus fail to capture GPP enhancement and its re-investment in green leaf area (Supplementary Table 2). The other four models behave randomly—some reproduce observed greenness levels (LAImax ~1.7 m2 m−2) but not the greening sensitivities (~0.045 m2 m−2 LAImax per 1 unit ω) and the others vice versa. Whether this is because these four models in common lack carbon–nitrogen interactions, or are missing some other key processes, is not known31, but the end result is a large range in model simulated greening sensitivity (hereafter LAImax sensitivity), during the historical period (0.022–0.075 m2 m−2 LAImax per 1 unit ω).
Emergent Constraint on projected increase of GPPWhat is known, however, is the strong linear relationship between modeled contemporaneous changes in LAImax and GPP arising from the combined radiative and physiological effects of CO2 enrichment in the range 1 × CO2 to 2 × CO2 (Supplementary Figs. 2 and 3). As a result, models with low LAImax sensitivity (Fig. 2b) project lower ΔGPP for a given increment of CO2 concentration, and vice versa. Thus, the large variation in modeled historical LAImax sensitivities (Fig. 2b) linearly maps to variation in ΔGPP at 2 × CO2 (Fig. 2c; r = 0.98, P = 0.0001), with the consequence that the uncertainty of the multi-model mean ΔGPP is large enough to undermine its value—e.g., 2.1 ± 1.91 Pg C yr−1 for 2 × CO2 in NHL. This linear relation in inter-model variation between ΔGPP at 2 × CO2 and historical LAImax sensitivities allows using the observed sensitivity as an EC on GPP estimation at 2 × CO2. Moreover, the probability contours about the best linear fit together with the uncertainty of observed sensitivity (blue and gray shadings in Fig. 2c) allow a robust characterization of the constrained estimate23, namely 3.4 ± 0.2 Pg C yr−1 for 2 × CO2 in NHL (Fig. 2d). This EC estimate is 60% larger than the multi-model mean value. Wenzel et al.23 reported a similar result for NHL (37 ± 9% versus 20–25% for relative GPP increase at 2 × CO2) and a somewhat smaller number for the extra-tropical vegetation in the northern hemisphere, both for the physiological effect only (Supplementary Fig. 4 shows that the radiative and physiological effects each contribute about half of the total GPP enhancement). Together, these results indicate that most models are largely underestimating photosynthetic carbon fixation, which is in contrast to previous studies5,7 that suggested an over-sensitivity of ESMs to atmospheric CO2. Below, we provide three independent lines of evidence, i.e., not using LAImax but atmospheric CO2 measurements, to buttress our EC estimate.
Independent lines of evidenceFirst, the seasonal cycle of CO2 concentration in the NHL, which shows a winter maximum due to respiratory processes and a late-summer minimum due to photosynthetic drawdown, may be considered as a proxy for NHL carbon exchange with the atmosphere15,16,17. Analyses of long-term measurements at NHL stations, Point Barrow (BRW, Alaska) and Alert Nunavut (ALT, Canada), reveal that this seasonal cycle has changed over time, dominated by a decreasing trend in the annual CO2 minimum (Fig. 3a, b). Nearly all of this change can be attributed to the land, as the trend in the abutting Arctic Ocean flux is ~15 times smaller (Fig. 3d; Methods). This strengthening of the seasonal swings of CO2 concentration relates to photosynthesis rather than respiration changes15,16,17 and thus features changes in GPP. So, if the EC estimate is closer to the true value of ΔGPP at 2 × CO2, then, models matching the EC estimate (e.g., MIROC-ESM) must also better simulate the changing CO2 seasonal cycle measured at the NHL stations, in comparison to models that over- (e.g., HadGEM2-ES) or underestimate (e.g., CESM1-BGC). Indeed, the MIROC-ESM best reproduces the average observed seasonal cycle, and critically, the change in summertime minimum over time at both stations, in comparison to the other models (Fig. 3a, b). None of the models reproduce the observed phase of the seasonal cycle, which suggests a recurring problem among models in their representation of vegetation phenology5. Nevertheless, the model that projects ΔGPP matching the LAImax-based EC estimate is also the one that best captures the changes in observed seasonal cycle suggesting that the EC estimate, rather than the corresponding multi-model mean, best represents the true value of ΔGPP at 2 × CO2. Thus, the multi-model mean is a large underestimate.
Observational LAI product (LAI3gV1)The new version (V1) of the LAI data set is an update of the widely used LAI3g data set26. It was generated using an artificial neural network (ANN) and the latest version (third generation) of the Global Inventory Modeling and Mapping Studies group (GIMMS) Advanced Very High Resolution Radiometer (AVHRR) normalized difference vegetation index (NDVI) data (NDVI3g). The latter has been corrected for sensor degradation, intersensor differences, cloud cover, solar zenith angle, viewing angle effects due to satellite drift, Rayleigh scattering, and stratospheric volcanic aerosols36. The ANN model was trained with overlapping data of NDVI3g and Collection 6 Terra MODIS LAI product37,38, and then applied to the full NDVI3g time series to generate the LAI3gV1 data set. This data set provides global and year-round LAI observations at 15-day (bi-monthly) temporal resolution and 1/12 degree spatial resolution from July 1981 to December 2016. Currently, it is the only data set that spans this long period. The quality of the previous version (V0) of the GIMMS LAI3g data set was evaluated through direct comparisons with ground-based measurements of LAI, indirectly with other estimates from similar satellite-data products, and also through statistical analysis with climatic variables, such as temperature and precipitation variability26. The LAI3gV0 data set (and related fraction vegetation-absorbed photosynthetically active radiation data set) has been widely used in various studies5,15,19,20,31,39,40,41. The new version LAI3gV1 used in our study is an update of that earlier version. For both, observational and CMIP5 data, LAI is defined as the one-sided green leaf area per unit ground area in broadleaf canopies and as one-half the green needle surface area in needleleaf canopies. It is expressed in units of m2 green leaf area per m2 ground area. In this study, we use the annual maximum value of LAI, LAImax, to quantify the greenness level of a surface. LAImax is less influenced by cloudiness and noise; accordingly, it is most useful in investigations of long-term greening and browning trends. The drawback of LAImax is the saturation effect at high LAI values42. However, this is less of a problem in high latitudinal ecosystems which are mostly sparsely vegetated, with LAImax values typically in the range of 2–3. The bi-monthly GIMMS LAI3gV1 data are merged to a monthly temporal resolution by averaging the two composites in the same month. Then, for model and observational data alike, the two-dimensional global fields are cropped to the northern high latitudinal band defined as 60°N to 90°N, averaged in space and temporally reduced to the annual maximum value. Although the AVHRR data underlying the LAI data in this study have corrections for various deleterious effects36, the data may still contain residual non-vegetation-related effects. Therefore, we sought confirmation of the greening trend19, on which the current study relies, from a more reliable but shorter record from the MODIS sensors37,38. These data are well calibrated, cloud-screened, and corrected for atmospheric effects, especially tropospheric aerosols. The sensor-platforms are regularly adjusted to maintain precise orbits. All algorithms, including the LAI algorithm, are physics-based, well-tested and currently producing the sixth generation data sets. The results, not shown here for brevity, illustrate global scale greening, across all latitudinal zones and broad vegetation classes. Zhang et al.43 also reported matching greening trends between the latest (Version 6) MODIS and AVHRR (Version 3) vegetation index data sets. We also found that the LAImax sensitivity derived with MODIS LAI data matched well with that obtained from the AVHRR LAI data (results not shown for brevity). Whether this indicates that the 17-year MODIS record from the period 2000 to 2016 captures information similar to the longer AVHRR record (1981–2016), or is simply a fortuitous occurrence, is not known, and deserves further study. In the present context, however, this adds confidence to the AVHRR LAI data used in our study.
Temperature data from ECMWF ERA-InterimEstimates of surface air temperature at 2 m height are from the widely used global atmospheric reanalysis product ERA-Interim by ECMWF44 (for details see https://www.ecmwf.int/en/research/climate-reanalysis/era-interim). The global temperature fields were retrieved at a resolution of 0.5° × 0.5° for monthly mean estimates derived from daily means. Other reanalysis products with similar specifications (NCEP/NCAR reanalysis, University of Delaware Air Temperature & Precipitation, and GHCN/CAMS reanalysis product) were also investigated. The differences among the various products were found to be minor.
CMIP5 models used in this studyIn this study, we analyze a set of the most recent climate-carbon simulations of seven ESMs participating in the fifth phase of the Coupled Model Intercomparison Project, CMIP528. The model data were obtained from the Earth System Grid Federation, ESGF (https://esgf-data.dkrz.de/projects/esgf-dkrz/). Seven ESMs provided output for the variables of interest for simulations esmHistorical, 1pctCO2, esmFixClim, and esmFdbk. The esmHistorical simulation spans the period 1850–2005 and was driven by observed conditions such as solar forcing, emissions or concentrations of short-lived species and natural and anthropogenic aerosols or their precursors, land use, anthropogenic as well as volcanic influences on atmospheric composition. The models are forced by prescribed anthropogenic CO2 emissions, rather than atmospheric CO2 concentrations. 1pctCO2 is an idealized fully coupled carbon/climate simulation initialized from steady state of the pre-industrial control run and atmospheric CO2 concentration prescribed to increase 1% yr−1 until quadrupling of the pre-industrial level. The simulations esmFixClim and esmFdbk and are set up as the 1pctCO2 with the difference, that in esmFixClim (esmFdbk) only the radiative effect from increasing CO2 concentration is included, while the carbon cycle sees the pre-industrial CO2 level (vice versa)28,45.
Historical simulation with MPI-ESM higher-resolution setupMPI-ESM-HR is the coupled high-resolution setup of the latest version of the Max-Planck-Institute Earth System Model MPI-ESM1.2, which is the baseline for the upcoming Coupled Model Intercomparison Project Phase 6 (CMIP6). Here, the atmospheric component ECHAM6.3 has 95 vertical levels and twice the horizontal resolution (~100 km) than the CMIP5 version. The ocean component MPIOM is set up on a tripolar grid at nominal 0.4° horizontal resolution (TP04) and 40 vertical levels. MPI-ESM1.2 includes the latest versions of the land and ocean carbon cycle modules, comprising the ocean biogeochemistry model HAMOCC and the land surface scheme JSBACH. The forcing components for the historical simulation are chosen from CMIP5 (Methods) as at the time the simulations were conducted CMIP6 forcing was not available46.
Atmospheric CO2 concentration dataMonthly means of atmospheric CO2 concentration at Point Barrow (71.3°N, 203.4°E) and Alert Nunavut (82.5°N, 297.7°E) are taken from the Global Monitoring Division measurement datasets (co2_brw_surface-insitu_1_ccgg_MonthlyData respectively co2_alt_surface-flask_1_ccgg_month) provided by the National Oceanic and Atmospheric Administration/Earth System Research Laboratory (NOAA/ESRL). Global monthly means of atmospheric CO2 concentration are taken from the GLOBALVIEW-CO2 product (obspack_co2_1_GLOBALVIEWplus_v2.1_2016_09_02; for details see https://doi.org/10.15138/G3259Z) also available at NOAA/ESRL.
Atmospheric CO2 inversion productsAtmospheric CO2 inversions estimate surface–atmosphere net carbon exchange fluxes by utilizing CO2 concentration measurements, a transport model and prior information on anthropogenic carbon emissions as well as carbon exchange between atmosphere and land respectively ocean47. We choose two products, which cover the longest time period (1980–2015) and are regularly updated, the Jena CarboScope32 (JENA, version s81_v3.8, for details see http://www.bgc-jena.mpg.de/CarboScope/s/s81_v3.8.html) and the Copernicus Atmosphere Monitoring Service33 (CAMS, version v15r2, for details see http://atmosphere.copernicus.eu/documentation-supplementary-products#greengas-fluxes) inversion systems. Both products provide monthly mean net flux estimates on a spatial resolution of 3.75° latitude and 5° longitude (JENA) and 1.875° latitude and 3.75° longitude (CAMS).
Calculation of growing degree days above 0 °C (GDD0)The global temperature fields from the reanalysis and model data are cropped to the northern high latitudinal band and averaged in space. The resulting one-dimensional time-series is converted to GDD above 0 °C by multiplying the days of each month with the respective monthly mean estimate if it is above 0 °C. Thus, we not only capture the warming signal, but also the prolongation of the growing season.
Dimension reduction using principal component analysisThe drivers GDD0 and atmospheric CO2 concentration vary co-linearly due to the radiative effect of increasing CO2 concentration in the NHL. Thus, it is problematic to conduct an accurate factor separation in terms of their respective contribution to increase in LAImax. However, the co-linearity suggests that a large amount of the signal is shared. Therefore, we conduct a PCA to apply dimension reduction48. The aim of the PCA is to find a linear combination of the driver variables that represents the one-dimensional projection with the largest possible variance. First, each driver time series xi is normalized by centering on its mean (subtracting x¯i
Estimation of historical LAImax sensitivityWe derive the historical LAImax sensitivity applying a standard linear regression model (fn)
Derivation of changes in NHL CO2 drawdown slopeGraven et al.17 showed that NHL CO2 drawdown mostly happens in June and July. ESMs, however, disagree on the phase, mainly due to a premature start of the growing season (Fig. 3a, b). As a consequence, the CO2 drawdown in models peaks earlier in the season. To obtain comparability for changes in CO2 drawdown strength, we calculate the first derivative of the CO2 concentration time series for the observational sites and each model individually. The annual minimum of the derivative in each time series reflects the months where the increase in photosynthetic CO2 fixation is strongest (CO2 drawdown slope). This procedure does not require a detrending of the atmospheric CO2 signal. For the BRW record, the 30 years of continuous overlap with the CMIP5 historical simulations were used to calculate the drawdown slopes (1974–2005). Due to the shorter overlap in the ALT record, 30 years of data from 1985 (start of measuring campaign) to 2015 were used for comparison with models. This is legitimate, because the CO2 concentration rate of increment for both periods are just about the same. Model time series are obtained from the near-surface CO2 concentration using the grid box in close proximity to each site. All yearly time series are slightly smoothed with a 2-year moving window. Changes are calculated from 5-year averages at the beginning and end of the record. Here, we only present a low-end, high-end and the closest-to-observation model from the greening EC analysis, because Wenzel et al.23 already reported the behavior of the entire CMIP5 ensemble in terms of simulating the NHL CO2 seasonal cycle.
Scaling of NPP estimatesWe scale and convert the EC estimate for changes in the GPP flux ΔFGPP,EC for a doubling of the pre-industrial CO2 level ([CO2]pi) to a NPP flux (ΔFNPP,EC) to obtain a comparable estimate to the atmospheric CO2 inversion datasets using
Comparison of C fluxes between Arctic Ocean and NHL landWe require the use of a fully coupled ESM to separate between land and ocean in terms of the sign, magnitude, and seasonal cycle of the respective net carbon exchange fluxes with the atmosphere. We have access to a spatially-high resolved historical run (10 realizations) of the MPI-ESM which has the ability to reproduce seasonality in the Arctic Ocean (Methods). The terrestrial carbon pools have not been brought into equilibrium due to computational limitations in these high-resolution simulations. Thus, we use simulations from the same model but at low spatial resolution (3 realizations), the CMIP5 esmHistorical simulation, to address land carbon exchange fluxes (Methods). The NHL land sink is approximately 2.5 times larger than the Arctic Ocean sink, on an annual basis. However, in terms of the change in carbon sink between the mid-1970s and early-2000s, the increase in CO2 uptake by the land is about 15 times larger than the ocean. Accordingly, the Arctic Ocean can be ignored when trying to explain changes during the recent past, i.e., BRW period of CO2 concentration measurements. During the months from May to September (may-to-sep) when photosynthetic CO2 drawdown is happening, the change in land sink is about 0.4 Pg C on an annual basis. Especially between May and July, the CO2 concentration is rapidly declining, i.e., photosynthesis prevails CO2 release processes. Thus, nearly the entire increase of 0.4 Pg C can be attributed to increasing NPP. The EC analysis shows that the MPI-ESM model is rather close to observations but generally underestimating greening sensitivity and thus also the GPP enhancement. These results are not provided as further proof of the EC estimate, although they are not contradictory, they are provided to compare the strength of NHL land and Arctic Ocean carbon sinks and why the ocean component can be neglected.
Bootstrapping for probability estimationWe apply bootstrapping to estimate the 68% confidence of the emergent linear relationship due to the small sample size of the CMIP5 ensemble. First, we randomly resample the data with replacement, where the size of the resample is equal to the size of the original sample N. Second, we compute the least-squares linear best fit for the resampled data. Third, we repeat this procedure m times (minimum m = 100) until the difference between the median best fit line of m − 1 and m computed regressions approaches zero (the actual threshold was set to a difference less than 1%). We derive the 68% confidence contours of equal probability based on the set of m random regression lines.
Calculation of probability density functionsWe derive a probability density function (PDF) for the observed sensitivity b′ (associated standard error σb, Methods)
Code availabilityThe code used in this study is available from the corresponding author upon request.
Quéré, C. L. et al. Global Carbon Budget 2017. Earth Syst. Sci. Data 10, 405–448 (2018).
Ciais, P. et al. in Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (eds. Stocker, T. F. et al.) 465–570 (Cambridge University Press, Cambridge and New York, 2013).
Zhao, F. & Zeng, N. Continued increase in atmospheric CO2 seasonal amplitude in the 21st century projected by the CMIP5 Earth system models. Earth Syst. Dyn. 5, 423–439 (2014).
Friedlingstein, P. et al. Uncertainties in CMIP5 climate projections due to carbon cycle feedbacks. J. Clim. 27, 511–526 (2013).
Anav, A. et al. Evaluating the land and ocean components of the global carbon cycle in the CMIP5 Earth system models. J. Clim. 26, 6801–6843 (2013).
Arora, V. K. et al. Carbon-concentration and carbon-climate feedbacks in CMIP5 Earth system models. J. Clim. 26, 5289–5314 (2013).
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).
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).
Campbell, J. E. et al. Large historical growth in global terrestrial gross primary production. Nature 544, 84–87 (2017).
Welp, L. R. et al. Interannual variability in the oxygen isotopes of atmospheric CO2 driven by El Nino. Nature 477, 579–582 (2011).
Koffi, E. N., Rayner, P. J., Scholze, M. & Beer, C. Atmospheric constraints on gross primary productivity and net ecosystem productivity: results from a carbon-cycle data assimilation system. Glob. Biogeochem. Cycles 26, GB1024 (2012).
Nemani, R. R. et al. Climate-driven increases in global terrestrial net primary production from 1982 to 1999. Science 300, 1560–1563 (2003).
Leakey, A. D. B. et al. Elevated CO2 effects on plant carbon, nitrogen, and water relations: six important lessons from FACE. J. Exp. Bot. 60, 2859–2876 (2009).
Thomas, R. T. et al. Increased light-use efficiency in northern terrestrial ecosystems indicated by CO2 and greening observations. Geophys. Res. Lett. 43, 11,339–11,349 (2016).
Forkel, M. et al. Enhanced seasonal CO2 exchange caused by amplified plant productivity in northern ecosystems. Science 351, 696–699 (2016).
Keeling, C. D., Chin, J. F. S. & Whorf, T. P. Increased activity of northern vegetation inferred from atmospheric CO2 measurements. Nature 382, 146–149 (1996).
Graven, H. D. et al. Enhanced seasonal exchange of CO2 by northern ecosystems since 1960. Science 341, 1085–1089 (2013).
Myneni, R. B., Keeling, C. D., Tucker, C. J., Asrar, G. & Nemani, R. R. Increased plant growth in the northern high latitudes from 1981 to 1991. Nature 386, 698–702 (1997).
Zhu, Z. et al. Greening of the Earth and its drivers. Nat. Clim. Change 6, 791–795 (2016).
Mao, J. et al. Human-induced greening of the northern extratropical land surface. Nat. Clim. Change 6, 959–963 (2016).
Cox, P. M. et al. Sensitivity of tropical carbon to climate change constrained by carbon dioxide variability. Nature 494, 341–344 (2013).
Wenzel, S., Cox, P. M., Eyring, V. & Friedlingstein, P. Emergent constraints on climate-carbon cycle feedbacks in the CMIP5 Earth system models. J. Geophys. Res. Biogeosci. 119, 794–807 (2014).
Wenzel, S., Cox, P. M., Eyring, V. & Friedlingstein, P. Projected land photosynthesis constrained by changes in the seasonal cycle of atmospheric CO2. Nature 538, 499–501 (2016).
Mystakidis, S., Davin, E. L., Gruber, N. & Seneviratne, S. I. Constraining future terrestrial carbon cycle projections using observation-based water and carbon flux estimates. Glob. Change Biol. 22, 2198–2215 (2016).
Screen, J. A. & Simmonds, I. The central role of diminishing sea ice in recent Arctic temperature amplification. Nature 464, 1334–1337 (2010).
Zhu, Z. et al. Global data sets of vegetation leaf area index (LAI)3g and fraction of photosynthetically active radiation (FPAR)3g derived from Global Inventory Modeling and Mapping Studies (GIMMS) normalized difference vegetation index (NDVI3g) for the period 1981 to 2011. Remote Sens. 5, 927–948 (2013).
Goetz, S. J., Bunn, A. G., Fiske, G. J. & Houghton, R. A. Satellite-observed photosynthetic trends across boreal North America associated with climate and fire disturbance. Proc. Natl. Acad. Sci. U.S.A. 102, 13521–13525 (2005).
Taylor, K. E., Stouffer, R. J. & Meehl, G. A. An overview of CMIP5 and the experiment design. Bull. Am. Meteorol. Soc. 93, 485–498 (2012).
Cook, B. I. & Pau, S. A global assessment of long-term greening and browning trends in pasture lands using the GIMMS LAI3g dataset. Remote Sens. 5, 2492–2512 (2013).
Park, T. et al. Changes in growing season duration and productivity of northern vegetation inferred from long-term remote sensing data. Environ. Res. Lett. 11, 084001 (2016).
Mahowald, N. et al. Projections of leaf area index in earth system models. Earth Syst. Dyn. 7, 211–229 (2016).
Rödenbeck, C., Houweling, S., Gloor, M. & Heimann, M. CO2 flux history 1982–2001 inferred from atmospheric data using a global inversion of atmospheric transport. Atmos. Chem. Phys. 3, 1919–1964 (2003).
Chevallier, F. et al. CO2 surface fluxes at grid point scale estimated from a global 21 year reanalysis of atmospheric measurements. J. Geophys. Res. 115, D21307 (2010).
Anav, A. et al. Spatiotemporal patterns of terrestrial gross primary production: a review. Rev. Geophys. 53, 785–818 (2015).
Ito, A. & Oikawa, T. A simulation model of the carbon cycle in land ecosystems (Sim-CYCLE): a description based on dry-matter production theory and plot-scale validation. Ecol. Model. 151, 143–176 (2002).
Pinzon, J. E. & Tucker, C. J. A non-stationary 1981–2012 AVHRR NDVI3g time series. Remote Sens. 6, 6929–6960 (2014).
Yan, K. et al. Evaluation of MODIS LAI/FPAR product collection 6. Part 1: consistency and improvements. Remote Sens. 8, 359 (2016).
Yan, K. et al. Evaluation of MODIS LAI/FPAR product collection 6. Part 2: validation and intercomparison. Remote Sens. 8, 460 (2016).
Piao, S. et al. Evidence for a weakening relationship between interannual temperature variability and northern vegetation activity. Nat. Commun. 5, 5018 (2014).
Poulter, B. et al. Contribution of semi-arid ecosystems to interannual variability of the global carbon cycle. Nature 509, 600 (2014).
Keenan, T. F. et al. Recent pause in the growth rate of atmospheric CO2 due to enhanced terrestrial carbon uptake. Nat. Commun. 7, 13428 (2016).
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 (2002).
Zhang, Y. et al. Reanalysis of global terrestrial vegetation trends from MODIS products: browning or greening? Remote Sens. Environ. 191, 145 (2017).
Dee, D. P. et al. The ERA-Interim reanalysis: configuration and performance of the data assimilation system. Q. J. R. Meteorol. Soc. 137, 553–597 (2011).
Taylor, K. E., Stouffer, R. J. & Meehl, G. A. A summary of the CMIP5 experiment design. PCDMI Rep. 33 (2009).
Müller, W. A. et al. A higher-resolution version of the Max Planck Institute Earth System Model (MPI-ESM1.2-HR). J. Adv. Model. Earth Syst. 10, 1383–1413 (2018).
Peylin, P. et al. Global atmospheric carbon budget: results from an ensemble of atmospheric CO2 inversions. Biogeosciences 10, 6699–6720 (2013).
Hannachi, A., Jolliffe, I. T. & Stephenson, D. B. Empirical orthogonal functions and related techniques in atmospheric science: a review. Int. J. Climatol. 27, 1119–1152 (2007).
Prentice, I. C. et al. in Climate Change 2001: the Scientific Basis. Contributions of Working Group I to the Third Assessment Report of the Intergovernmental Panel on Climate Change (eds. Houghton, J. T. et al.) 185–237 (Cambridge University Press, Cambridge and New York, 2001).
Zhang, Y., Xu, M., Chen, H. & Adams, J. Global pattern of NPP to GPP ratio derived from MODIS data: effects of ecosystem type, geographical location and climate. Glob. Ecol. Biogeogr. 18, 280–290 (2009).