Global carbon cycle data

We used global carbon budget data from the Global Carbon Project1, in combination with diverse observational data sets from satellite remote sensing and distributed Earth observation networks, multiple prognostic DGVMs and a diagnostic modelling approach. The Global Carbon Project1 data set is a compilation of estimates of all major components of the global carbon budget, based on the combination of a range of data, algorithms, statistics and model estimates. This data set provided the long-term estimates of global emissions, along with estimates of the residual terrestrial sink from 10 DGVMs. Long-term atmospheric CO 2 concentrations were provided by the National Oceanic and Atmospheric Administration’s (NOAA) Earth System Research Laboratory ( http://www.esrl.noaa.gov). Annual emissions of CO 2 from the burning of fossil fuels, land use change and cement production from the Global Carbon Project1 were used in conjunction with the NOAA atmospheric CO 2 concentration data to calculate the annual airborne fraction.

Diagnosing changes in the growth rate of atmospheric CO 2

We used two methods to examine changes in the growth rate of atmospheric CO 2 over time: a statistical linear model of the growth rate as a function of emissions, atmospheric CO 2 concentrations and global CO 2 sinks46. Rayner et al.46 showed that the growth rate can be modelled as a linear function of atmospheric CO 2 concentration, and that deviations of the atmospheric growth rate from this linear model are an indication of changes in the sink strength of the biosphere. This is an informative approach as it incorporates the physical link between atmospheric CO 2 concentrations and global CO 2 sinks, thus allowing changes in that coupling to be studied. By examining the residuals of this linear model we identify a significant residual bias from 2002 to 2014, indicating a structural change in the relationship between atmospheric CO 2 concentrations, emissions and global CO 2 sinks (Supplementary Fig. 1). The second method, singular spectrum analysis (SSA), is used to extract the underlying temporal dynamics of both the growth rate and the airborne fraction.

Examining changes in atmospheric CO 2 with linear modelling

The growth rate can be modelled as a linear function of atmospheric CO 2 concentrations46. Deviations of the atmospheric growth rate from this linear model are an indication of changes in the sink strength of the biosphere. Such an approach is more informative than a simple linear model as a function of time, as it incorporates the physical link between atmospheric CO 2 concentrations and global CO 2 sinks, thus allowing changes in that coupling to be studied.

Formally, the global CO 2 sink can be described as a linear function of CO 2 concentration or, equivalently, CO 2 mass. Thus, we write

where M is the mass of CO 2 in the atmosphere and M 0 is the background or equilibrium mass of CO 2 in the atmosphere. Equation (1) can be simplified to

where B has units of per year and plays the role of an inverse residence time for excess carbon against the processes of land and ocean uptake.

Given that

where is the growth rate of atmospheric CO 2 . We can substitute equation (2) into equation (3) to get

Using observations of emissions from fossil fuel burning and land use change (F anthro ), along with the atmospheric CO 2 concentrations and growth rate, it is possible to estimate B and F 0 using standard statistical techniques.

The model thus constructed preserves the relationship between atmospheric CO 2 concentrations, emissions and the growth rate of atmospheric CO 2 , on the basis of a proportional response of the global CO 2 sinks. Any change in the strength of the global sinks can therefore be analysed by examining the residuals between the observed and predicted growth rate. We use the model, informed by the first 30 years of observations to test the hypothesis that the CO 2 growth rate maintains the same linear relationship with atmospheric CO 2 concentrations throughout. The residuals show a statistically significant deviation from being normally distributed around zero from 2002 (P<0.05, t-test, Supplementary Fig. 1). This shows that there is a pause in the growth rate of atmospheric CO 2 . Importantly, it also identifies changes in global CO 2 sinks as the cause. We confirm this change in the growth rate by using SSA to extract the low frequency mode of variability corresponding to 5 years (see below).

Examining changes in atmospheric CO 2 with SSA

SSA is a non-parametric spectral estimation method for extracting different modes of variability form a time series. The SSA method decomposes time series into a sum of components, each having a meaningful interpretation. Different modes of underlying variability can then be extracted by reconstructing the time series using only the eigenvalues relevant to the mode of variability in question. The name singular spectrum relates to the spectrum of eigenvalues in a singular value decomposition of a covariance matrix.

SSA first decomposes the time series into a set of empirical orthogonal functions and associated principal components in the spectral domain, following Takens’ embedding theorem47. Each component is dominated by a single oscillatory mode, and therefore has a simple representation in the frequency domain. This means that each empirical orthogonal function can thus be assigned a characteristic frequency. The functional separation from the decomposition step can be used to reconstruct the time series, either fully or only retaining specific modes of variability by using the relevant principal components. Thus, a time series can be finally described by a set of subsignals, X f , each of which belongs to a well-defined frequency bin. All SSA analyses in this paper were performed using a 5-year frequency.

Despite of orthogonal base functions, extracted subsignals are subject to uncertainty due to a degree of inseparability of closely related modes of variability48. To quantify this uncertainty, we use a surrogate technique (the Iterative Amplitude Adjusted Fourier Transform49), following Mahecha et al.50 This approach generates a set of surrogates for each residual, corresponding to the extracted 5-year subsignal. The subsignal X 5 is then re-extracted 100 times, giving an array of subsignal surrogates. The s.d. of the extracted subsignals quantifies the extraction uncertainty, and is shown in Fig. 1.

DGVMs of the terrestrial carbon cycle

DGVM output from two different model intercomparison projects was used in this analysis (Supplementary Table 1). As part of the Global Carbon Project1, transient runs from 10 DGVMs are available from 1959 to 2014 (www.globalcarbonproject.org/carbonbudget). These model simulations were assessed using observed atmospheric CO 2 growth rates (Supplementary Fig. 2) and used to explore the decadal change in the land sink from 1959 to 2014 (Fig. 2). See Le Quere et al.1 for a detailed description of the models. We also use global simulations from 8 DGVMs (Supplementary Table 1) run as part of the Trends in Net Land-Atmosphere Exchange (TRENDY-v1) project (http://dgvm.ceh.ac.uk/node/9). In this project, common input forcing data were prescribed for a series of model experiments from 1901 to 2010. Here we use two model experiments from TRENDY-v1, varying either CO 2 only (S1 (ref. 7): time-invariant climate; present-day land use mask) or climate only (S4 (ref. 51): time-invariant CO 2 ; present-day land use mask). For more details on the TRENDY model simulations see Sitch et al.7

Satellite-based estimates of the terrestrial carbon cycle

In addition to the DGVMs, we used satellite-based estimates of ecosystem GPP combined with ecosystem respiration estimates from a diagnostic coupled PR model to quantify the likely effect of changes in global vegetation (through the fAPAR), and changes in water availability (though Alpha) on enhanced terrestrial uptake, from 1982 to 2013.

The diagnostic coupled PR model is based on a new light use efficiency (LUE) model of photosynthesis developed from first principles, and a semi-empirical model of ecosystem respiration developed based on eddy-covariance flux data.

The mechanistic photosynthesis model proposed by Farquhar et al.52 successfully captures the biochemical controls of leaf photosynthesis and responses to variations in temperature, light and CO 2 concentration. According to the model, the gross photosynthesis rate, A, is limited by either the capacity of the Rubisco enzyme for the carboxylation of RuBP (ribulose-1,5-bisphosphate), or by the electron transport capacity for RuBP regeneration.

In the case of Rubisco limitation, the photosynthetic rate (A c ) is given by

where V cmax is the maximum rate of Rubisco activity, c i is the intercellular concentration of CO 2 , Γ* is the CO 2 compensation point in the absence of dark respiration and K is the Michaelis–Menten coefficient of Rubisco.

In the case of limitation by the electron transport capacity for RuBP regeneration, and assuming electron transport capacity is large (relative to V cmax ) such that the response of photosynthesis to light is linear under Rubisco limitation, the photosynthetic rate (A j ) is given by

where ϕ 0 is the intrinsic quantum efficiency and I is the absorbed light.

The co-limitation or coordination hypothesis, which is strongly supported by empirical evidence53,54, predicts that photosynthesis under typical daytime field conditions is close to the point where Rubisco- and electron transport-limited photosynthesis rates are equal (that is, equation (5)=equation (6))55. In other words, the photosynthetic capacity of leaves adjusts to acclimate to the typical daytime light levels to be neither in sufficient excess to induce additional, non-productive maintenance respiration nor less than required for full exploitation of the available light. Recent empirical support comes from Maire et al.53, who tested the coordination hypothesis with 293 observations for 31 species grown under a range of environmental conditions, and found that average daily photosynthesis under field conditions is close to the point, where the Rubisco and electron transport photosynthesis rates are equal.

The coordination hypothesis allows for the prediction of photosynthesis through equation (6) using a LUE approach. Indeed the success of LUE models generally in predicting photosynthesis can be explained by the co-ordination hypothesis. Importantly, it also allows for the effect of CO 2 on photosynthesis to be incorporated in such LUE models based on the first-principles understanding of the Farquhar et al.52 model.

By rewriting equation (6), substituting c i by the product of atmospheric CO 2 (c a ) and the ratio of leaf-internal to ambient CO 2 (χ=c i /c a ), GPP can be described as

where ϕ 0 is the quantum yield and I is the absorbed light (here derived from satellite fAPAR).

Γ* depends on temperature, as estimated through a biochemical rate parameter (x) as described by Bernacchi et al.56:

here R is the molar gas constant (8.314 J mol−1 K−1) and x 25 =4.22 Pa is the photorespiratory point at 25 °C. ΔH is the activation energy for Γ* (37,830 J mol−1) and T is the temperature in K.

χ depends on air temperature and the vapour pressure deficit (VPD; D), and can be estimated following the least-cost hypothesis54. This states that an optimal long-term effective value of χ can be predicted as a result of plants minimizing their total carbon costs associated with photosynthetic carbon gain, and explicitly expressed with the following model

where D is VPD, K is the Michaelis–Menten coefficient of Rubisco and η* is the viscosity of water relative to its value at 25 °C. D is estimated as the difference between saturated and actual vapour pressure. Saturated vapour pressure (e s ) is estimated as the averaged saturated vapour pressure at maximum and minimum temperature with the Clausius–Clapeyron relationship, which is well approximated by

The Michaelis–Menten coefficient of Rubisco (K) in equation (9) is given by

where K c and K o are the Michaelis–Menten coefficient of Rubisco for carboxylation and oxygenation, respectively, expressed in partial pressure units, and P o is the partial pressure of O 2 . K responds to temperature via K c and K o , which is also described by equation (8) with specific parameters: ΔH is 79.43 kJ mol−1 for K c and 36.38 kJ mol−1 for K o , x 25 is 39.97 Pa for K c and 27,480 Pa for K o . On the basis of equations (9)–(11), and assuming a typical value of χ 25 as 0.8 (at T=298.15 K and D=1 kPa), parameter β of equation (9) is estimated as 356.51.

The GPP model implicitly assumes that nutrient limitations on GPP are manifest in allocation to foliage and are therefore contained in the observed fAPAR, as has been reported in recent empirical observations57,58. This is consistent with the theoretical expectation and empirical evidence that CO 2 -induced enhancement of biomass growth is possible even under nutrient-limited conditions59, and findings of increased below-ground allocation, including root exudation, on less fertile soils60.

Ecosystem respiration (R eco ) in the diagnostic model is estimated via a photosynthesis-dependent respiration model61, which combines the joint influences of temperature, soil moisture and substrate availability on ecosystem respiration, and is designed for the diagnostic upscaling of R eco from observations at eddy-covariance towers to global scales. The model estimates monthly R eco as

where R 0 is the reference respiration rate at the reference temperature T ref (15 °C), E 0 is the activation energy and T 0 =−46.02 °C (ref. 61). k is the proportional contribution of GPP to ecosystem respiration through substrate availability at the reference temperature T ref , and is a measure of water availability, calculated as the ratio of actual to equilibrium evapotranspiration (see below). Conceptually, this model can be considered as the sum of a GPP-dependent term comprising autotrophic respiration and the fast-responding labile component of heterotrophic respiration, and a GPP-independent term standing for heterotrophic respiration of slower carbon pools. The two free model parameters (E 0 , k) were taken from the original study61, where 104 globally distributed sites from the FLUXNET network were used to derive plant functional type (PFT) specific parameters. Global PFT classifications were taken from the MODIS land cover product (MOD12, http://modis-land.gsfc.nasa.gov/landcover.html), curated at a resolution of 0.5° by the Global Land Cover Facility (http://glcf.umd.edu/data/lc/). For each 0.5° grid cell, we used the PFT that was most prevalent during the period 2000–2013. R 0 was estimated for each 0.5° grid cell analytically, by solving the combined equations (7) and (12) for equilibrium net carbon uptake under preindustrial conditions.

Diagnostic model simulations performed

To examine the role of changes in each of the model drivers (air temperature, atmospheric CO 2 , radiation, moisture availability and vegetation cover) used in our analysis, we ran multiple global 0.5° simulations from 1900 to 2013 with the PR model. For each simulation, we removed the long-term trend in all drivers but one. This allowed us to quantify the direct first-order effect of long-term changes in each. An fAPAR climatology was used pre-1981.

Diagnostic model forcing data

Global monthly gridded weather data at 0.5° were provided by the Climate Research Unit at East Anglia University (CRU TS3.21)62. The total available photosynthetically active radiation, VPD and the ratio of actual to equilibrium evapotranspiration (α) were calculated from insolation and CRU climate data using a simple process-based bioclimatic model (STASH63). The GIMMS 3G remotely sensed Normalized Difference Vegetation Index product64 provided monthly estimates of the fAPAR, an indicator of green vegetation cover, at 0.5°. Monthly fAPAR estimates are available from 1981 to the present.

Diagnostic model evaluation

We evaluated the efficacy of the PR model at multiple temporal and spatial scales. Evaluations performed include the following: a global comparison to interannual variability in the residual terrestrial sink from multiple DGVMs and estimates from the Global Carbon Project65 (Supplementary Fig. 2), site-level comparisons to time series of GPP and R eco from globally distributed individual sites in the La Thuile Fair Use FLUXNET data set (Supplementary Figs 5–8 and Supplementary Table 2); comparisons with the seasonal anomalies of GPP and R eco from 149 sites from the same data set (Supplementary Fig. 9); regional comparisons with seasonal changes in atmospheric CO 2 concentrations from the NOAA global sampling stations (Supplementary Fig. 10); and a latitudinal comparisons to an empirical upscaling estimate of global GPP66 (Supplementary Fig. 11).

To compare PR model estimates of NEP to the observations at the NOAA stations, we used the TM2 atmospheric transport model67 to integrate and transport detrended monthly values of NEP for each 0.5° grid cell to the station locations. We then calculated both the modelled and observed CO 2 seasonal cycle at the observation sites68.

Derivation of the sensitivity of GPP to CO 2

GPP can be described as a function of atmospheric CO 2 as

The sensitivity of GPP to c a can therefore be derived by taking the derivative of GPP with respect to c a , as

This requires the derivation of , which can be formulated as:

It therefore follows that

Taking χ=0.8 and Γ*=43 p.p.m. (at 25 °C), the sensitivity of GPP is thus calculated as 37% at current levels of atmospheric CO 2 (400 p.p.m.).

Examining the first derivative of the LUE model of GPP54,58 suggests a CO 2 sensitivity (β CO2 ) of GPP of 37% at current atmospheric CO 2 levels (400 p.p.m.), which is consistent with the observed response in FACE studies69. Our estimate of the change in GPP is also consistent with other process-based estimates7,8, but is larger than estimates from commonly used empirical upscaling techniques (Supplementary Fig. 12) as these do not account for the effect of increasing atmospheric CO 2 on photosynthesis.

Long-term change in the sensitivity of GPP to changing c a

The change in the sensitivity of GPP to c a can also be derived analytically.

If we denote the sensitivity (from equation (16)) as

Then we can calculate the partial derivative of β CO2 as

And the partial derivative of gamma as

Combining equations (18) and (19) gives

The negative coefficient implies that the sensitivity of GPP to c a , β CO2 , will decline with increasing c a . It is also clear from equation (20) that the response of β CO2 to c a is not linear, but decreases with c a . In other words, the magnitude of the sensitivity declination decreases with c a enhancement. Evaluating β CO2 at different atmospheric CO 2 concentrations shows a decrease in β CO2 from 37% under current CO 2 levels to 19% at double the current CO 2 levels (Supplementary Table 3).

Data availability

Data used in this study are available from the Global Carbon Project data archive at the Carbon Dioxide Information Analysis Center (http://cdiac.ornl.gov/GCP/). This includes global carbon budget data and long-term simulations from DGVMs. Additional simulations used are available from the Trends in Net Land-Atmosphere Exchange (TRENDY) project (http://dgvm.ceh.ac.uk/node/9). All other data that support the findings of this study are available from the corresponding author on request.

Code availability

The code used in this study is available from the corresponding author on request.