Trend analysis of the airborne fraction and sink rate of anthropogenically released CO₂

Is the fraction of anthropogenically released CO₂ that remains in the atmosphere (the airborne fraction) increasing? Is the rate at which the ocean and land sinks take up CO₂ from the atmosphere decreasing? We analyse these questions by means of a statistical dynamic multivariate model from which we estimate the unobserved trend processes together with the parameters that govern them. We show how the concept of a global carbon budget can be used to obtain two separate data series measuring the same physical object of interest, such as the airborne fraction. Incorporating these additional data into the dynamic multivariate model increases the number of available observations, thus improving the reliability of trend and parameter estimates. We find no statistical evidence of an increasing airborne fraction, but we do find statistical evidence of a decreasing sink rate. We infer that the efficiency of the sinks in absorbing CO₂ from the atmosphere is decreasing at approximately 0.54 % yr⁻¹.

1 Introduction

A part of the anthropogenically released CO₂ emitted to the atmosphere flows to the oceans (the ocean sink) and the terrestrial biosphere (the land sink). Approximately 45 % of released CO₂ stays in the atmosphere (the airborne fraction), while the two sinks take up approximately 24 % and 31 % of the CO₂, respectively. (These percentages are calculated over the period 1959 to 2016 using the data described below; see, for example, Raupach et al., 2014, for similar estimates.) A key question is whether the airborne fraction is increasing or if it remains constant at around 45 %. An increasing airborne fraction implies that the share of anthropogenically released CO₂ that ultimately remains in the atmosphere increases, and projections of future atmospheric CO₂ levels need to take this into account (Gloor et al., 2010). Closely related is the question of whether the sinks will continue taking up CO₂ at the same rate (the sink rate) or if this rate is decreasing. A decreasing sink rate implies that the efficiency with which ocean and land sinks are absorbing CO₂ from the atmosphere is decreasing. Thus, analysing the behaviour of the sink rate can help predict the future uptake of CO₂ through the ocean and the land sink. The answers to the questions posed above are important for our understanding of the global carbon cycle and are relevant for policymakers and the public in general.

A series of papers argue that the airborne fraction of anthropogenically released CO₂ (mainly through fossil fuel emissions, cement production, and land-use change) is increasing (Canadell et al., 2007a; Le Quéré et al., 2009; Raupach et al., 2008; Rayner et al., 2015). Similarly, in Raupach et al. (2014), it is argued that, although the statistical evidence of an increasing airborne fraction is relatively weak, the evidence of a decreasing CO₂ sink rate is clearer. However, the methods in these studies have been criticized in, for example, Knorr (2009), Gloor et al. (2010), and Ballantyne et al. (2015). Indeed, by considering a longer data set and incorporating uncertainties into the data, Knorr (2009) found that the conclusion of an increasing airborne fraction was not warranted. Similarly, Ballantyne et al. (2015) argues that errors in the data can lead to erroneous conclusions regarding possible trends in the airborne fraction and in the sink rate.

In this paper, we address these statistical issues within the framework of a state-space system. It allows us to conduct statistical inference by taking explicit account of stochastic and deterministic trends in the data, transient shocks to the data (coming from, e.g. volcanic eruptions or strong El Niño events), and (potential) measurement errors. It also allows for the simultaneous incorporation of multiple data sets for the same object, which can improve the estimation of trends and increase reliability of parameter estimates. We find strong evidence for purely deterministic trends when we incorporate multiple measurements for the airborne fraction and the sink rate. These deterministic trends have a statistically significantly negative slope in the case of the sink rate and an insignificant slope in the case of the airborne fraction. These findings corroborate earlier findings in the literature, especially those of Raupach et al. (2014), but address the statistical concerns raised by Knorr (2009) and Ballantyne et al. (2015), among others. Finally, by analysing the ocean and land sink rates separately, we find no evidence of a decreasing ocean sink rate, but we do find evidence that the land sink rate is decreasing.

The paper is organized as follows. In Sect. 2 we state the fundamental equations of the global carbon budget, the definitions of the airborne fraction of anthropogenically released CO₂, and the CO₂ sink rate, which will motivate the specification of the state-space system. Section 3 introduces the state-space system used in the paper. In Sect. 4 we conduct a trend analysis of the airborne fraction within the proposed statistical framework. In Sect. 5 we carry out the corresponding analysis of the CO₂ sink rate, and in Sect. 6 we carry out the analysis of the land and ocean sink rates separately. Section 7 discusses the results, and Sect. 8 concludes the paper. The Supplement is available online.

2 The global carbon budget

The so-called global carbon budget is defined as

$$\begin{array}{}\text{(1)}& {\displaystyle }{E}_t^{\mathrm{ANT}}=G_t+S_t^O+S_t^L,\end{array}$$

where $E_t^{\mathrm{ANT}}$ is anthropogenically released CO₂ into the atmosphere, G_t is growth of atmospheric CO₂ concentration, $S_t^O$ is the flux of CO₂ from the atmosphere to the oceans (the ocean sink), and $S_t^L$ is the flux of CO₂ from the atmosphere to the terrestrial biosphere (the land sink). In other words, Eq. (1) states that emissions of CO₂ should equal the fluxes of CO₂ to the atmosphere, the ocean sink, and the land sink. We use the data set provided by the Global Carbon Project (Le Quéré et al., 2018).¹ All data are given in gigatonnes of carbon (GtC) and are recorded at a yearly frequency, beginning in 1959 and ending in 2016, resulting in 58 observations for each quantity in Eq. (1).

While the carbon budget is in principle always balanced for the physical quantities, in the sense that Eq. (1) always holds, this might not be the case when inserting actual data for emissions and sinks due to measurement errors in the data. For this reason, Le Quéré et al. (2018) introduce a residual term into the budget Eq. (1) to capture measurement error. It is denoted $B_t^{\mathrm{IM}}$ for budget imbalance. Therefore, when considering actual data, the carbon budget is defined as

$$\begin{array}{}\text{(2)}& {\displaystyle }{E}_t^{\mathrm{ANT}}=G_t+S_t^O+S_t^L+B_t^{\mathrm{IM}}.\end{array}$$

The sample mean of the budget imbalance over the observation period is not significantly different from zero and shows no sign of a trend (Le Quéré et al., 2018). These facts are important in the developments below, since they motivate treating $B_t^{\mathrm{IM}}$ as part of an error term.

The growth rate in atmospheric CO₂ data, G_t, is from Dlugokencky and Tans (2018), the ocean sink data, $S_t^O$, are obtained from an ensemble of global biochemistry models, and the land sink data, $S_t^L$, are estimated as the multi-model mean of several dynamic global vegetation models. See Le Quéré et al. (2018) for further information on the data. The anthropogenic emissions of CO₂ can be decomposed into two parts:

$${\displaystyle }{E}_t^{\mathrm{ANT}}=E_t^{\mathrm{FF}}+E_t^{\mathrm{LUC}},$$

where $E_t^{\mathrm{FF}}$ are emissions from fossil fuel burning, cement production, and gas flaring, while $E_t^{\mathrm{LUC}}$ are emissions from land-use change. Fossil fuel emissions, $E_t^{\mathrm{FF}}$, are from Boden et al. (2018), while land-use change emissions, $E_t^{\mathrm{LUC}}$, are averages of the results of the two bookkeeping models of Hansis et al. (2015) and Houghton and Nassikas (2017), updated as in Le Quéré et al. (2018). The time series of concentrations (above pre-industrial levels) of CO₂ in the atmosphere is constructed as

$$C_t=\mathrm{2.127}\cdot \left(\right[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{1959}}-\left[{\mathrm{CO}}_{\mathrm{2}}{]}_{\mathrm{1750}}\right)+\sum _{\mathit{\tau }=\mathrm{1}}^t{G}_{\mathit{\tau }},$$

where [CO₂]₁₇₅₀=279 ppmv (parts per million by volume) and [CO₂]₁₉₅₉=315.39 ppmv are the concentrations of CO₂ in the atmosphere in 1750 and 1959, respectively; see Raupach et al. (2014). The number 2.127 is the conversion factor from parts per million by volume to gigatonnes of carbon. In other words, the atmospheric concentration C_t above pre-industrial levels is given by the initial value in 1959 plus the cumulative sum of the growth in atmospheric CO₂ concentrations G_t, which result from the budget Eq. (1).

We follow Raupach (2013) and Raupach et al. (2014) and define the airborne fraction as

$${\mathrm{AF}}_t={\displaystyle \frac{G_t}{E_t^{\mathrm{ANT}}}},$$

and the CO₂ sink rate as

$$\begin{array}{}\text{(3)}& {\displaystyle }{k}_{\mathrm{S},t}={\displaystyle \frac{S_t^O+S_t^L}{C_t}},\end{array}$$

which is the flux of CO₂ from the atmosphere to the sinks (ocean plus land), normalized by the amount of CO₂ (above pre-industrial levels) currently in the atmosphere. We can also consider the individual components of the sink rate for ocean and land, which are given by

$$\begin{array}{}\text{(4)}& {\displaystyle }{k}_{O,t}={\displaystyle \frac{S_t^O}{C_t}},k_{L,t}={\displaystyle \frac{S_t^L}{C_t}},\end{array}$$

respectively, with $k_{\mathrm{S},t}=k_{O,t}+k_{L,t}$.

The airborne fraction and the sink rate are fundamentally different quantities. The airborne fraction AF${}_t=G_t/E_t^{\mathrm{ANT}}$ is the ratio of the growth of atmospheric CO₂ in period t to the amount of CO₂ emitted in period t. It is thus a measure of the fraction of emitted CO₂ that stays in the atmosphere. In contrast, the sink rate $k_{\mathrm{S},t}=(S_t^O+S_t^L)/C_t$ is the ratio of the CO₂ flux in the sinks in period t to the total amount of CO₂ in the atmosphere (above pre-industrial levels). By writing $S_t^O+S_t^L=k_{\mathrm{S},t}{C}_t$, we can interpret the sink rate k_S,t as the “efficiency” with which CO₂ flows from the atmosphere to the sinks, i.e. as the amount of CO₂ going into the sinks for an extra unit of CO₂ added to the atmosphere (Gloor et al., 2010; Raupach, 2013). We discuss the relationship between the airborne fraction and the sink rate further in Sect. 7.

3 Trend model specification

In this section, we consider several models for the data-generating process behind observations of the objects of interest defined in Sect. 2. Common to all models is that they can be cast in a state-space system of the form

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\mathit{y}}_t=\mathbf{A}{\mathit{x}}_t+{\mathit{\xi }}_t,\\ & {\mathit{x}}_{t+\mathrm{1}}=\mathbf{B}{\mathit{x}}_t+{\mathit{\kappa }}_t,\end{array}\end{array}$$

where y_t is a vector of observations at time $t=\mathrm{1},\mathrm{\dots },n$ with time series length n, and the system matrices A and B have appropriate dimensions. The vector x_t is usually referred to as the state vector, which can include deterministic and stochastic trends, and the error terms ξ_t and κ_t are both independent and identically distributed (iid) random vectors of appropriate dimensions and with a mean of zero. For example, when we need to model the airborne fraction alone, we have y_t=AF_t, and the state-space system represents a univariate dynamic model for the airborne fraction. When modelling the ocean and land sink rates jointly, we have ${\mathit{y}}_t=(k_{O,t},k_{L,t}{)}^{\prime }$, and the state-space system is a bivariate dynamic model. For given matrices A and B, and under the assumption of mutually and serially uncorrelated Gaussian errors ξ_t and κ_t (with their respective variance matrices Σ_ξ and Σ_κ), the state-space system is a linear Gaussian model. In such regular cases, an analytic formulation for the likelihood function is available and relies on the prediction error decomposition. Hence the parameters (variances and possibly covariances in Σ_ξ and Σ_κ) can be estimated by the maximum likelihood method. It requires the numerical optimization of the log-likelihood function that is evaluated via the Kalman filter. The resulting algorithm is initialized with specific starting values; we use a diffuse initialization as outlined in Chapter 5 of Durbin and Koopman (2012). The smooth estimate of the state process x_t can also be obtained by means of the Kalman filter together with a smoothing algorithm. The extracted state is effectively the conditional mean $\mathbb{E}\left({\mathit{x}}_t\right|{\mathit{y}}_{\mathrm{1}},\mathrm{\dots },{\mathit{y}}_n;\mathbf{A},\mathbf{B},{\mathbf{\Sigma }}_{\mathit{\xi }},{\mathbf{\Sigma }}_{\mathit{\kappa }})$, for $t=\mathrm{1},\mathrm{\dots },n$. Details of the state-space approach are given by Durbin and Koopman (2012), where both signal extraction and maximum likelihood estimation are discussed.

Our baseline model is the local linear trend (LLT) model. For a univariate time series y_t, we treat the underlying trend T_t as a stochastic process given by

$$\begin{array}{}\text{(6)}& {\displaystyle }{T}_{t+\mathrm{1}}=T_t+\mathit{\beta }+{\mathit{\eta }}_t,\end{array}$$

where β∈ℝ is a fixed and unknown coefficient and η_t is an iid Gaussian random variable with a mean of zero and variance ${\mathit{\sigma }}_{\mathit{\eta }}^{\mathrm{2}}$. The solution to the difference equation in Eq. (6) is given as

$$T_{t+\mathrm{1}}=T_{\mathrm{1}}+t\mathit{\beta }+\sum _{i=\mathrm{0}}^{t-\mathrm{1}}{\mathit{\eta }}_{t-i},t=\mathrm{1},\mathrm{2},\mathrm{\dots },n-\mathrm{1},$$

where T₁ can be treated as a fixed unknown coefficient (intercept or constant) or as a random variable. The solution shows that the trend component is made up of the starting value T₁, a deterministic linear term with slope β, and a random-walk component ${\sum }_{i=\mathrm{0}}^{t-\mathrm{1}}{\mathit{\eta }}_{t-i}$. Thus, T_t can be interpreted as a long-term trend in the time series and β as the slope of the deterministic part of the trend. We also considered a time-varying slope, β_t, but found no evidence supporting this generalization in either the airborne fraction or the sink rate. The observation equation for y_t is given by

$$\begin{array}{}\text{(7)}& {\displaystyle }{y}_t=T_t+{\mathit{ϵ}}_t,\end{array}$$

where T_t is given by Eq. (6) and ϵ_t captures deviations of the observed time series from the unobserved trend component. The deviations ϵ_t can be viewed as (i) actual (transient) disturbances of the physical systems arising from, for example, volcanic eruptions and El Niño events, and/or (ii) measurement errors arising from the way the data are collected. The random variable ϵ_t is assumed to be iid Gaussian with a mean of zero and variance ${\mathit{\sigma }}_{\mathit{ϵ}}^{\mathrm{2}}$.

The local linear trend model can be cast in the state-space system Eq. (5) where vectors and matrices are defined as

$$\begin{array}{rl}& {\mathit{x}}_t=\left(\begin{array}{c}{T}_t\\ \mathit{\beta }\end{array}\right),\mathbf{A}=\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\end{array}\right],\mathbf{B}=\left[\begin{array}{cc}\mathrm{1}& \mathrm{1}\\ \mathrm{0}& \mathrm{1}\end{array}\right],\\ & {\mathit{\xi }}_t={\mathit{ϵ}}_t,{\mathit{\kappa }}_t=\left(\begin{array}{c}{\mathit{\eta }}_t\\ \mathrm{0}\end{array}\right),\end{array}$$

for $t=\mathrm{1},\mathrm{\dots },n$. The state vector x_t consists of the two variables of interest: stochastic trend variable T_t and deterministic slope variable β. The state-space methods as discussed above can treat such mixed compositions of the state vector. We have illustrated how the state-space system can be used for a univariate time series. In the next sections, we also consider trend analyses based on multivariate time series models.

4 Trend analysis of the airborne fraction

It follows immediately from Eq. (2) that we can measure the airborne fraction AF_t in two alternative ways:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\mathrm{AF}}_t^{\left(\mathrm{1}\right)}={\displaystyle \frac{G_t^{\mathrm{ATM}}}{E_t^{\mathrm{ANT}}}},\\ & {\mathrm{AF}}_t^{\left(\mathrm{2}\right)}={\displaystyle \frac{E_t^{\mathrm{ANT}}-S_t^O-S_t^L}{E_t^{\mathrm{ANT}}}}={\mathrm{AF}}_t^{\left(\mathrm{1}\right)}+{\mathit{\xi }}_t,\end{array}\end{array}$$

where ${\mathit{\xi }}_t=B_t^{\mathrm{IM}}/E_t^{\mathrm{ANT}}$, since $E_t^{\mathrm{ANT}}-S_t^O-S_t^L=G_t+B_t^{\mathrm{IM}}$. Although the two quantities in Eq. (8) measure the same underlying object (the airborne fraction AF_t), they differ in practice because of a non-zero budget imbalance, i.e. ξ_t≠0. Our statistical analysis implies that ξ_t is a well-behaved zero-mean and covariance stationary error process.

We consider our baseline local linear trend model of Sect. 3 for each of the objects, that is,

$${\mathit{y}}_t={\mathrm{AF}}_t^{\left(i\right)}=T_t^{\left(i\right)}+{\mathit{ϵ}}_t^{\left(i\right)},$$

for $i=\mathrm{1},\mathrm{2}$, where the trend $T_t^{\left(i\right)}$ is specified in Eq. (6) and with error ${\mathit{ϵ}}_t^{\left(i\right)}$. Table 1 reports the output of the estimation, using the state-space system and the Kalman filter. The first part of Table 1 presents estimates of the standard deviations of the observation error term ${\mathit{ϵ}}_t^{\left(i\right)}$ and the trend error term ${\mathit{\eta }}_t^{\left(i\right)}$, as well as the estimate of the slope parameter β, including the estimated standard error (SE) of $\widehat{\mathit{\beta }}$ and the resulting t statistic, $t\text{-stat}=\widehat{\mathit{\beta }}/\text{SE}\left(\widehat{\mathit{\beta }}\right)$. Based on these estimation results, we can formally test hypotheses of the type

$$\begin{array}{}\text{(9)}& H_{\mathrm{0}}:\mathit{\beta }=\mathrm{0}\text{against}{H}_{\mathrm{1}}:\mathit{\beta }\ne \mathrm{0},\end{array}$$

or, more relevantly,

$$\begin{array}{}\text{(10)}& H_{\mathrm{0}}:\mathit{\beta }=\mathrm{0}\text{against}{H}_{\mathrm{1}}:\mathit{\beta }>\mathrm{0}.\end{array}$$

By using the normal approximation to the t distribution and for a 95 % confidence level, the critical value for the test Eq. (9) is 1.96, and for Eq. (10), it is 1.645. In the case of the airborne fraction, we are interested in testing Eq. (10). It is evident from Table 1 that we cannot reject H₀ in this case (p values 0.2711 and 0.4042 for the case of ${\mathrm{AF}}_t^{\left(\mathrm{1}\right)}$ and ${\mathrm{AF}}_t^{\left(\mathrm{2}\right)}$, respectively). In other words, although the estimate $\widehat{\mathit{\beta }}$ is positive, we cannot conclude, statistically at 95 % confidence, that the airborne fraction is increasing over time.

Table 1 also contains diagnostic statistics for the standardized prediction residual u_t based on

$${\mathit{y}}_t-\mathbb{E}\left({\mathit{y}}_t\right|{\mathit{y}}_{\mathrm{1}},\mathrm{\dots },{\mathit{y}}_{t-\mathrm{1}};\mathbf{A},\mathbf{B},{\mathbf{\Sigma }}_{\mathit{\xi }},{\mathbf{\Sigma }}_{\mathit{\kappa }}),$$

for $t=\mathrm{1},\mathrm{\dots },n$, and where Σ_ξ and Σ_κ are replaced by their respective maximum likelihood estimates. Under the assumption that the local linear trend model is correctly specified for the time series y_t, the residuals u_t are Gaussian iid (see Durbin and Koopman, 2012, p. 38). To verify these properties of u_t empirically, we consider two residual diagnostic statistics: the normality test statistic N of Jarque and Bera (1987) and the serial correlation test statistic of Durbin and Watson (1971). As a goodness-of-fit statistic, we consider the $R_d^{\mathrm{2}}$, which is a relative measure of model fit against a random-walk model. Since the statistic is defined in a similar way to the standard regression fit measure R², we have

$$\begin{array}{rl}& R_d^{\mathrm{2}}=\mathrm{1}-{\displaystyle \frac{{\sum }_{t=\mathrm{2}}^n{u}_t^{\mathrm{2}}}{{\sum }_{t=\mathrm{2}}^n\left[\right(y_t-y_{t-\mathrm{1}})-m{]}^{\mathrm{2}}}},\\ & m=(n-\mathrm{1}{)}^{-\mathrm{1}}\sum _{t=\mathrm{2}}^n(y_t-y_{t-\mathrm{1}}).\end{array}$$

The reported diagnostic statistics and goodness of fit in Table 1 are satisfactory for the time series AF${}_t^{\left(\mathrm{1}\right)}$ and AF${}_t^{\left(\mathrm{2}\right)}$. We may conclude from these results that the local linear trend model from Eqs. (6)–(7) provides an adequate description of the dynamic features in the time series. Since the AF${}_t^{\left(\mathrm{2}\right)}$ is well-described within our state-space framework, the extra error term ${\mathit{\xi }}_t=B_t^{\mathrm{IM}}/E_t^{\mathrm{ANT}}$ in AF${}_t^{\left(\mathrm{2}\right)}$, as introduced by the budget imbalance term in Eq. (8), is well-behaved. Hence the assumptions underlying the state-space system appear to be valid.

Table 1 Univariate analysis of the airborne fraction.

We report parameter estimates for the standard deviations σ_ϵ and σ_η and slope coefficient β together with its standard error (SE) and t statistic (t-stat). We further report the normality (N) test, the goodness-of-fit statistic $R_D^{\mathrm{2}}$, and the Durbin–Watson (DW) test statistic for serial correlation; all are computed for the standardized prediction errors u_t, which are obtained from the Kalman filter. The normality test N is the χ² distributed, with 2 degrees of freedom, statistic of Jarque and Bera (1987), with its 95 % critical value of 5.99; the statistic relies on the sample estimates of skewness and kurtosis of u_t. The goodness-of-fit statistic $R_d^{\mathrm{2}}$ is defined as $\mathrm{1}-ESS/DSS$, where ESS $={\sum }_{t=\mathrm{2}}^n{u}_t^{\mathrm{2}}$ and DSS $={\sum }_{t=\mathrm{2}}^n\left[\right({\mathit{y}}_t-y_{t-\mathrm{1}})-m{]}^{\mathrm{2}}$ with $m=(n-\mathrm{2}{)}^{-\mathrm{1}}{\sum }_{t=\mathrm{2}}^n({\mathit{y}}_t-y_{t-\mathrm{1}})$. The DW test statistic is developed by Durbin and Watson (1971), where also its critical values are tabulated. If DW =2 the sequence u_t is serially uncorrelated, if DW <2 there is evidence that the errors u_t are positively autocorrelated, and if DW >2 there is evidence that the errors u_t are negatively autocorrelated.

The state-space system allows both measures for the airborne fraction, AF${}_t^{\left(\mathrm{1}\right)}$ and AF${}_t^{\left(\mathrm{2}\right)}$, to be included in a single model with the purpose of improving the quality of the trend estimation and inference. We begin with an “uninformed” system using two different trend components, $T_t^{\left(\mathrm{1}\right)}$ and $T_t^{\left(\mathrm{2}\right)}$, both specified as Eq. (6), for the two time series. We have

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\mathit{y}}_t& =\left[\begin{array}{c}{\mathrm{AF}}_t^{\left(\mathrm{1}\right)}\\ {\mathrm{AF}}_t^{\left(\mathrm{2}\right)}\end{array}\right]=\left[\begin{array}{c}{G}_t^{\mathrm{ATM}}/E_t^{\mathrm{ANT}}\\ \mathrm{1}-(S_t^{\mathrm{OCEAN}}+S_t^{\mathrm{LAND}})/E_t^{\mathrm{ANT}}\end{array}\right]\\ & =\left[\begin{array}{c}{T}_t^{\left(\mathrm{1}\right)}\\ T_t^{\left(\mathrm{2}\right)}\end{array}\right]+\left[\begin{array}{c}{\mathit{ϵ}}_t^{\left(\mathrm{1}\right)}\\ {\mathit{ϵ}}_t^{\left(\mathrm{2}\right)}\end{array}\right],\end{array}\end{array}$$

where the error terms ${\mathit{ϵ}}_t^{\left(i\right)}$, for $i=\mathrm{1},\mathrm{2}$, are correlated and their correlation coefficient can be estimated by the method of maximum likelihood together with the other parameters. The estimation results for this model are presented in panel a of Table 2. The main difference from Table 1 is the inclusion of the estimated correlation matrix for $({\mathit{ϵ}}_t^{\left(\mathrm{1}\right)},{\mathit{ϵ}}_t^{\left(\mathrm{2}\right)})$. The diagnostic test statistics are reasonable. In comparison with the univariate analysis, the goodness-of-fit values for $R_d^{\mathrm{2}}$ are slightly higher for the multivariate model. Hence we trust the model to be a good representation of the data. Furthermore, the slope is estimated to be positive in both cases (that is $\widehat{\mathit{\beta }}>\mathrm{0}$). However, when testing the null hypothesis given in Eq. (10), we cannot reject the hypothesis that the slopes are zero (p values 0.3753 and 0.4895 for the case of ${\mathrm{AF}}_t^{\left(\mathrm{1}\right)}$ and ${\mathrm{AF}}_t^{\left(\mathrm{2}\right)}$, respectively).

Table 2 Multivariate analysis of the airborne fraction.

We report parameter estimates for the standard deviations ${\mathit{\sigma }}_{\mathit{ϵ}}^{\left(i\right)}$ and ${\mathit{\sigma }}_{\mathit{\eta }}^{\left(i\right)}$, for $i=\mathrm{1},\mathrm{2}$, correlation matrix for ϵ_t, and slope coefficient β together with its standard error (SE) and t statistic (t-stat). We further report the normality (N) test, the goodness-of-fit statistic $R_D^{\mathrm{2}}$, and the Durbin–Watson (DW) test statistic for serial correlation. For details, see Table 1. In panel b, the trend coefficients (σ_η and β) for AF⁽²⁾ are the same as for AF⁽¹⁾ given the construction of the model (Eq. 12).

Since the two quantities in Eq. (8) measure the same object, the airborne fraction, we now force the state-space system to recognize that these data are driven by the same underlying common trend, $T_t^A$, but with possibly different error terms ${\mathit{ϵ}}_t^{\left(\mathrm{1}\right)}$ and ${\mathit{ϵ}}_t^{\left(\mathrm{2}\right)}$. In other words, we consider

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\mathit{y}}_t& =\left[\begin{array}{c}{\mathrm{AF}}_t^{\left(\mathrm{1}\right)}\\ {\mathrm{AF}}_t^{\left(\mathrm{2}\right)}\end{array}\right]=\left[\begin{array}{c}{G}_t^{\mathrm{ATM}}/E_t^{\mathrm{ANT}}\\ \mathrm{1}-(S_t^{\mathrm{OCEAN}}+S_t^{\mathrm{LAND}})/E_t^{\mathrm{ANT}}\end{array}\right]\\ & =\left[\begin{array}{c}{T}_t^A\\ T_t^A\end{array}\right]+\left[\begin{array}{c}{\mathit{ϵ}}_t^{\left(\mathrm{1}\right)}\\ {\mathit{ϵ}}_t^{\left(\mathrm{2}\right)}\end{array}\right].\end{array}\end{array}$$

The output of the estimation of this system is shown in panel b of Table 2; the estimated common trend and the data are plotted in Fig. 1. A slight deterioration of the diagnostic statistics is to be expected when introducing a common trend into the system, but the diagnostic statistics are still such that we can accept Eq. (12) as a plausible model. For the estimate of the slope $\widehat{\mathit{\beta }}$, we find a larger t statistic in absolute value than in the uninformed model, indicating that the restriction to the common trend increases the precision of the estimates. An explanation of this finding is that the informed system used twice as many observations for estimating the trend compared to the uninformed system. The hypothesis test in Eq. (10) reveals that the estimate of the slope parameter, although again positive, is still not statistically different from zero (p value 0.2199).

Figure 1 Estimated trend $T_t^A$ of the airborne fraction from the model 12.

5 Trend analysis of the CO₂ sink rate

In this section, we analyse the CO₂ sink rate in the same way as the airborne fraction above. Analogously we can define two alternative versions of the sink rate:

$$\begin{array}{}\text{(13)}& k_{S,t}^{\left(\mathrm{1}\right)}={\displaystyle \frac{S_t^O+S_t^L}{C_t}},k_{S,t}^{\left(\mathrm{2}\right)}={\displaystyle \frac{E_t^{\mathrm{ANT}}-G_t}{C_t}}=k_{\mathrm{S},t}^{\left(\mathrm{1}\right)}+{\mathit{\xi }}_t,\end{array}$$

where now ${\mathit{\xi }}_t=B_t^{\mathrm{IM}}/C_t$ and where we used Eq. (2). As was the case for the airborne fraction, these two quantities measure the same underlying object (the sink rate, k_S,t) but differ in practice because of a non-zero budget imbalance, i.e. ξ_t≠0.

The basic (univariate) local linear trend model for each of these objects is then given by

$$y_t=k_{S,t}^{\left(i\right)}=T_t^{\left(i\right)}+{\mathit{ϵ}}_t^{\left(i\right)},$$

for $i=\mathrm{1},\mathrm{2}$, where $T_t^{\left(i\right)}$ is specified as in Eq. (6). When the model is cast in the state-space system, the parameters can be estimated for each of the data series individually. The estimation results are presented in Table 3. The diagnostic statistics are satisfactory, and we conclude again that the error term ${\mathit{\xi }}_t=B_t^{\mathrm{IM}}/C_t$ is well-behaved in the sense that the assumptions underlying the state-space system appear to be valid also for the alternative sink rate data, $k_{\mathrm{S},t}^{\left(\mathrm{2}\right)}$. Even though the estimates of the slopes are negative, we cannot reject the null hypothesis of β=0 (p values 0.2233 and 0.0761 for the case of $k_{S,t}^{\left(\mathrm{1}\right)}$ and $k_{S,t}^{\left(\mathrm{2}\right)}$, respectively). We still consider a one-sided test as in Eq. (10), but now the relevant alternative hypothesis is $H_{\mathrm{1}}:\mathit{\beta }<\mathrm{0}$.

Table 3 Univariate analysis of the CO₂ sink rate.

We report parameter estimates for the standard deviations σ_ϵ and σ_η and slope coefficient β together with its standard error (SE) and t statistic (t-stat). We further report the normality (N) test, the goodness-of-fit statistic $R_D^{\mathrm{2}}$, and the Durbin–Watson (DW) test statistic for serial correlation; all are computed for the standardized prediction errors u_t, which are obtained from the Kalman filter; for details see Table 1.

Analogously to the airborne fraction above, these data can be put in a joint uninformed system with two different trend components, and we have

$$\begin{array}{}\text{(14)}& {\displaystyle }{\mathit{y}}_t=\left[\begin{array}{c}{k}_{S,t}^{\left(\mathrm{1}\right)}\\ k_{S,t}^{\left(\mathrm{2}\right)}\end{array}\right]=\left[\begin{array}{c}(S_t^O+S_t^L)/C_t\\ (E_t^{\mathrm{ANT}}-G_t)/C_t\end{array}\right]=\left[\begin{array}{c}{T}_t^{\left(\mathrm{1}\right)}\\ T_t^{\left(\mathrm{2}\right)}\end{array}\right]+\left[\begin{array}{c}{\mathit{ϵ}}_t^{\left(\mathrm{1}\right)}\\ {\mathit{ϵ}}_t^{\left(\mathrm{2}\right)}\end{array}\right],\end{array}$$

which can be compared with the model in Eq. (11). The estimation results for this model are reported in panel a of Table 4. Although the slope estimates are negative, they are not significant (p values 0.3106 and 0.1947 for the case of $k_{S,t}^{\left(\mathrm{1}\right)}$ and $k_{S,t}^{\left(\mathrm{2}\right)}$, respectively).

Table 4 Multivariate analysis of the CO₂ sink rate.

We report parameter estimates for the standard deviations ${\mathit{\sigma }}_{\mathit{ϵ}}^{\left(i\right)}$ and ${\mathit{\sigma }}_{\mathit{\eta }}^{\left(i\right)}$, for $i=\mathrm{1},\mathrm{2}$, the correlation matrix for ϵ_t, and slope coefficient β together with its standard error (SE) and t statistic (t-stat). We further report the normality (N) test, the goodness-of-fit statistic $R_D^{\mathrm{2}}$, and the Durbin–Watson (DW) test statistic for serial correlation; for details, see Table 1. In panel b, the trend coefficients (σ_η and β) for $k_S^{\left(\mathrm{2}\right)}$ are the same as for $k_S^{\left(\mathrm{1}\right)}$ given the construction of the model (Eq. 15).

Finally, we consider the state-space system that imposes a common trend for both time series, $T_t^S$, that is,

$$\begin{array}{}\text{(15)}& {\displaystyle }{\mathit{y}}_t=\left[\begin{array}{c}{k}_{S,t}^{\left(\mathrm{1}\right)}\\ k_{S,t}^{\left(\mathrm{2}\right)}\end{array}\right]=\left[\begin{array}{c}(S_t^O+S_t^L)/C_t\\ (E_t^{\mathrm{ANT}}-G_t)/C_t\end{array}\right]=\left[\begin{array}{c}{T}_t^S\\ T_t^S\end{array}\right]+\left[\begin{array}{c}{\mathit{ϵ}}_t^{\left(\mathrm{1}\right)}\\ {\mathit{ϵ}}_t^{\left(\mathrm{2}\right)}\end{array}\right],\end{array}$$

which can be compared with model in Eq. (12). The estimation results are presented in panel b of Table 4. Similar to the analysis of the airborne fraction in the previous section, the diagnostic statistics are somewhat worse for the less flexible system with a common trend. However, the diagnostics are still satisfactory, while the goodness-of-fit statistics improved overall. The estimate of the slope is

$${\displaystyle }\widehat{\mathit{\beta }}=-\mathrm{0.00014},$$

and this estimate is statistically significant: we reject the hypothesis $H_{\mathrm{0}}:\mathit{\beta }=\mathrm{0}$ in favour of $H_{\mathrm{1}}:\mathit{\beta }<\mathrm{0}$ at a 95 % confidence level (p value 0.0014). The mean of the sink rate (calculated using either data set $k_S^{\left(\mathrm{1}\right)}$ or $k_S^{\left(\mathrm{2}\right)}$) is 0.0258. It follows that we estimate the sink rate to be decreasing with approximately $\mathrm{0.00014}/\mathrm{0.0258}=\mathrm{0.54}\mathit{\%}$ yr⁻¹. The estimated trend and the data are plotted in Fig. 2.

The state-space system is also well-suited for forecasting; see Durbin and Koopman (2012). Using the model in Eq. (15), we forecast the sink rate 25 years ahead in time. The output is presented in Fig. 3. For reference, the forecasts coming from an autoregressive model of order 1 (AR1) are also presented. The downward trend in the forecasts from the state-space model is the result of the negative estimate of β. Under current conditions, the forecast implies that in approximately 15 years, the sink rate will have declined to below 2 %. Conversely, the autoregressive model produces forecasts that converge to the mean of the original data series.

It is important to recognize that the validity of these forecasts are conditional on the assumption that the sink rate, k_S, is linear in concentrations. As seen from the analysis above, see also Fig. 2, this assumption has been approximately satisfied over the time period considered in this paper, but whether it will continue to be accurate is an open question (see Appendix A for a discussion of the possible future behaviour of the sink rate). The model-based forecasts of Fig. 3 should be seen in this light: these forecasts are obtained under the assumption that the sink rate will continue to be approximately linear in concentrations. Whether this assumption is reasonable is an interesting question beyond the scope of the present study.

Figure 2 Estimated trend $T_t^S$ of the CO₂ sink rate from Model (Eq. 15).

Figure 3 The blue solid line represents the data, while the red solid line represents the point forecasts from the Kalman filter with the unknown parameters estimated by maximum likelihood. The dashed red lines are 95 % confidence bands (±1.96 standard deviation) for the forecasts. The green line represents the forecasts from an autoregressive model of order 1.

6 Trend analysis of the ocean and land sink rates

We may conclude from the analysis in the previous section that the combined (land plus ocean) sink rate appears to be decreasing. To investigate this finding in more detail, we study two alternative models, which consider the two sink variables separately. The first model specifies local linear trends for ocean and land sink rates, i.e.

$$\begin{array}{}\text{(16)}& {\displaystyle }{\mathit{y}}_t=\left[\begin{array}{c}{k}_{O,t}\\ k_{L,t}\end{array}\right]=\left[\begin{array}{c}{S}_t^O/C_t\\ S_t^L/C_t\end{array}\right]=\left[\begin{array}{c}{T}_t^O\\ T_t^L\end{array}\right]+\left[\begin{array}{c}{\mathit{ϵ}}_t^{\left(\mathrm{1}\right)}\\ {\mathit{ϵ}}_t^{\left(\mathrm{2}\right)}\end{array}\right],\end{array}$$

where the time series k_O,t and k_L,t are defined in Eq. (4), while the trend variables $T_t^O$ and $T_t^L$ are specified as in Eq. (6). To inform the state-space system of the structure of the carbon budget, we also consider the model equations

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\mathit{y}}_t& =\left[\begin{array}{c}{k}_{O,t}\\ k_{L,t}\\ k_{\mathrm{S},t}\end{array}\right]=\left[\begin{array}{c}{S}_t^O/C_t\\ S_t^L/C_t\\ (E_t^{\mathrm{ANT}}-G_t)/C_t\end{array}\right]=\left[\begin{array}{c}{T}_t^O\\ T_t^L\\ T_t^O+T_t^L\end{array}\right]\\ & +\left[\begin{array}{c}{\mathit{ϵ}}_t^{\left(\mathrm{1}\right)}\\ {\mathit{ϵ}}_t^{\left(\mathrm{2}\right)}\\ {\mathit{ϵ}}_t^{\left(\mathrm{3}\right)}\end{array}\right].\end{array}\end{array}$$

This trivariate model equation can be cast in the state-space system (Eq. 5) as well. The model specification has two independent trend processes of the form Eq. (6) for land and ocean sinks. Since $k_{\mathrm{S},t}=k_{O,t}+k_{L,t}$, the time series k_S,t of combined ocean and land sinks must feature the sum of the two trend processes for the individual sinks as its trend process.

The estimation results for these two model specifications are presented in Table 5. The residual diagnostic statistics N and DW are satisfactory, but we are particularly interested in the estimates of the slope parameters. It seems that most of the decrease in the sink rate can be attributed to the land sink. The slope estimates of the trend driving the ocean sink rate are very close to zero and not statistically significant (p values 0.5227 and 0.5168, respectively). On the other hand, the slope estimates of the trend driving the land sink rate are negative for both specifications. In the first model (Eq. 16), we can reject the hypothesis that the slope of the trend driving the land sink rate is zero in favour of the one-sided alternative $H_{\mathrm{1}}:\mathit{\beta }<\mathrm{0}$ at a 95 % confidence level (p value of 0.0420). For the more informed model specification (Eq. 17), the estimation results are reported in panel b of Table 5. Here we can reject H₀ at a 90 % confidence level (p value of 0.0882). Further, the results show that the estimate of the slope parameter from the land sink rate is equal to the estimate of the slope parameter from the combined sink rate as in Sect. 5, that is, $\widehat{\mathit{\beta }}=-\mathrm{0.00014}$. In other words, it appears that the decrease in the combined sink rate studied in the previous section is entirely explained by the decrease in the land sink rate.

Table 5 Analysis of ocean and land sink rates.

We report parameter estimates for standard deviations ${\mathit{\sigma }}_{\mathit{ϵ}}^{\left(i\right)}$ and ${\mathit{\sigma }}_{\mathit{\eta }}^{\left(i\right)}$, for $i=\mathrm{1},\mathrm{2},\mathrm{3}$, correlation matrix for ϵ_t, and slope coefficient β together with its standard error (SE) and t statistic (t-stat). We further report the normality (N) test, the goodness-of-fit statistic $R_D^{\mathrm{2}}$ and the Durbin–Watson (DW) test statistic for serial correlation; for details see Table 1. In panel b, we have two trends and two sets of trend coefficients (σ_η and β) for k_O,t and k_L,t. The trend for k_S,t is a combination of the two given the construction of the model (Eq. 17).

7 Discussion

Previous studies of the airborne fraction and the CO₂ sink rate have focused on detecting a single linear and deterministic trend in the data of the form a₀+a₁t, where a₀ and a₁, are constants (Canadell et al., 2007a; Le Quéré et al., 2009; Knorr, 2009; Raupach et al., 2008, 2014). However, possible statistical difficulties in such analyses have been pointed out in Knorr (2009). For instance, a linear regression analysis of two or more non-stationary variables can yield invalid inference (Granger and Newbold, 1974). The approach of this paper is to consider the data in a state-space system. In this way, non-stationary components are explicitly modelled as unobserved trend components and inference is valid (e.g. Durbin and Koopman, 2012). Furthermore, the trend specification of the state-space system allows for both deterministic and stochastic trend components.

In some of the uninformed models (cf. Table 1, panel a of Table 2, and panel a of Table 4), we estimate ${\widehat{\mathit{\sigma }}}_{Slp}>\mathrm{0}$, and, thus, in these cases, we find evidence of the trend component varying in time. However, in our “informed” models with a single trend object for two alternative time series, the extracted trends are practically deterministic, that is, the estimates of σ_Slp in panel b of Tables 2 and 4 are near zero (cf. also Figs. 1 and 2). In conclusion, there is evidence that a simple deterministic trend fits both the airborne fraction and the sink rate data well, although this only becomes evident when incorporating two data sets for each of these objects.

Several studies have highlighted the need for accounting for noise in measurements of climate-related data (Knorr, 2009; Ballantyne et al., 2015). The state-space approach explicitly incorporates such noise in the framework as well. Ballantyne et al. (2015) argue that errors in $E_t^{\mathrm{ANT}}$ might be autocorrelated. As shown in Tables 1–5, the diagnostic statistics do not indicate that autocorrelated errors pose a serious problem in our specifications. Nevertheless, the state-space framework can incorporate autocorrelated errors in the measurement equation.

Why do we find statistical evidence of a decreasing CO₂ sink rate but no evidence of an increasing airborne fraction when these two quantities are closely linked and the data entering the analyses are the same? It was noted in Gloor et al. (2010) that the airborne fraction and the sink rate are actually not as closely linked as they may seem prima facie. In particular, an increasing airborne fraction does not necessarily imply a decreasing sink rate and vice versa (Gloor et al., 2010, Section 8). The findings of this paper support this claim by providing statistical evidence for a constant airborne fraction but at the same time for a decreasing sink rate. Secondly, the concept of an airborne fraction is that of a long-term quantity: the airborne fraction should represent the amount of anthropogenically released CO₂ that eventually stays in the atmosphere after other fluxes have been taken into account. However, the ratio of the concurrent fluxes, i.e. $G_t/E_t^{\mathrm{ANT}}$, is likely a very noisy measurement of this object. Also, as we saw above, it is reasonable to consider sink fluxes, and therefore indirectly G_t, as being dependent on the level of CO₂ in the atmosphere (i.e. $C_t=\sum G_t$), which is not captured by the concurrent ratio $G_t/E_t^{\mathrm{ANT}}$. When studying the airborne fraction, it would perhaps be more reasonable to study an object taking this cumulative nature into account, e.g. $\sum G_t/\sum E_t^{\mathrm{ANT}}=C_t/\sum E_t^{\mathrm{ANT}}$ (in fact, such specifications were often considered in earlier parts of the literature; e.g. Keeling, 1973; Bacastow and Keeling, 1979; Oeschger and Heimann, 1983; Enting and Pearman, 1986). However, cumulative statistics of this type would present other difficulties. The dominance of the long-term history may mask sudden changes, for example. In contrast, the sink rate S_t∕C_t, as a flow-to-stock ratio, is immediately compatible with the underlying theory, at least as long as the linear approximation of the relationship between S_t and C_t, such as was made in, for example, Gloor et al. (2010) and Rayner et al. (2015), is adequate.

What are possible physical explanations for the apparent decrease in the sink rate? Raupach (2013) argues that a necessary condition for a constant sink rate is that the so-called “LinExp” assumption holds, i.e. that the sink fluxes $S_t^O$ and $S_t^L$ are linear in concentrations C_t (“Lin”) and that emissions ($E_t^{\mathrm{ANT}}$) grow exponentially (“Exp”). Constancy of the airborne fraction rests on a similar LinExp argument. Since we find no statistical evidence that the airborne fraction, AF_t, and the ocean sink rate, k_O,t, are non-constant in time, it is unlikely that the Exp assumption is grossly violated over the observation period considered in this paper. In contrast, it was found above that the efficiency of the land sink, k_L,t, is decreasing. A plausible explanation of these findings is that the Lin assumption no longer holds for the land sink, for instance because the terrestrial sink could be slowly saturating (Canadell et al., 2007b). In Appendix A we give a formal argument for how this could lead to the findings documented above. In particular, we show that the findings of this paper can be explained by the land sink's response to high atmospheric CO₂ concentrations: it is plausible that due to a rising level of CO₂ concentration, non-linear effects in the terrestrial CO₂ carbon cycle have become noticeable. If this is indeed the case, it has obvious consequences for our understanding of the carbon cycle and should be a cause for substantial concern (Gloor et al., 2010, p. 7740). However, although this explanation of our findings is consistent with the data, we can not conclude that it is the only possible explanation. Further research into the underlying reasons for the decreasing sink rate would be very valuable and is left for future work.

It is possible that the analyses conducted above are influenced by external natural events such as the El Niño–Southern Oscillation (ENSO), volcanic eruptions, and the like (Frölicher et al., 2013). The state-space system used in this paper can explicitly account for such effects through the additive error terms ϵ_t (cf. Eq. 5). To verify that the approach is indeed robust to such external and transitory events, we have also conducted our analyses using 5-year average data. The findings from the estimated state-space system for these time series of averages confirm those reported above: in the joint estimation, we find no statistical evidence of a trend in the airborne fraction (p value of 0.3214), and we do find statistical evidence of a decreasing trend in the sink rate (p value of 0.00064). We conclude that the findings of this paper are not likely to be driven by external natural events such as ENSO and volcanic eruptions. We also considered 2-, 3-, and 4-year averages with similar results. We present details of this analysis in the Supplement (Sect. 3). To further check the robustness of the results, we examined whether there are any observations in the data set which are particularly influential. Statistically influential observations could be due to outliers, caused for instance by external natural events, such as the ones mentioned above. Using a statistic called Cook's distance (Cook, 1977, 1979; Atkinson et al., 1997), which is a measure of how influential a given observation is on the analysis, we did not find evidence of any one observation being particularly influential. Similarly, we tried estimating the slope parameter β after deleting the tth observation for each time point t in the sample, i.e. for $t=\mathrm{1959},\mathrm{1960},\mathrm{\dots },\mathrm{2016}$; the estimates of the slope parameter found in this way were very stable, which is further evidence of the robustness of the analyses to potential outliers and external events. Details can be found in the Supplement (Sect. 2.1).

This paper considers data recorded at a yearly frequency, while many of the previous studies of the airborne fraction and the sink rate use monthly data. The advantage of using monthly data is obvious: more observations. However, there are also some disadvantages. For instance, while the CO₂ concentration C_t (and therefore also the growth rate G_t) is recorded every month, these data contain a strong seasonal component induced by the photosynthesis–respiration cycle of terrestrial vegetation. This seasonality needs to be accounted for in some way; for instance, Raupach et al. (2014) smooth the data using a 15-month running mean. In contrast, some of the other data are recorded only yearly, for instance, the emission data available to us, $E_t^{\mathrm{ANT}}$. In this case, Raupach et al. (2014) use linear interpolation to get monthly estimates of emissions. Such transformations of the data, i.e. smoothing or interpolation, might introduce new and complicated errors, possibly invalidating the analyses. For these reasons, we prefer to work with yearly data.

8 Conclusions

We have argued that the state-space system can be a useful approach for analysing possible trends in the airborne fraction of anthropogenically released CO₂ and in the CO₂ sink rate. We have shown that deterministic and stochastic trend processes can be explicitly and jointly incorporated as unobserved components, allowing for valid inference, even when the observed time series are non-stationary. The state-space framework also allows for the incorporation of multiple data sets for the same object, which increases reliability of the resulting estimates.

We estimate a positive, yet statistically insignificant, slope in the data for the airborne fraction. Using two alternative time series for the sink rate and imposing a common trend, we obtain a significantly negative deterministic trend. Our analyses support the conclusions as set out by Raupach et al. (2014): the rate at which the combined (ocean plus land) sink takes up CO₂ from the atmosphere seems to be decreasing. The best estimate resulting from our state-space system is that the CO₂ sink rate, and therefore the efficiency with which the combined land and ocean sink is absorbing carbon from the atmosphere, is decreasing by 0.54 % yr⁻¹. We do not find evidence of this rate itself changing over time.

Finally, there is tentative evidence that the decrease in the sink rate is mainly driven by a weakening uptake in the land sink. This could be the result of non-linearities in the response of the terrestrial carbon sink to the level of atmospheric concentrations of CO₂. That is, although the land sink is itself increasing and thus continuing to take up a large part of anthropogenically emitted CO₂, as also noted recently by, for example, Rayner et al. (2015), Keenan et al. (2016), and Fernández-Martínez et al. (2019), the rate of this uptake appears to be decreasing. The statistical evidence for this is not strong, however, and we suggest that additional research must be conducted to further investigate this question.

Appendix A: Linear approximation of the relation of land sink and concentrations

In this Appendix, we argue that the levels of atmospheric concentrations of CO₂ may have risen to a point where a linear expansion of the logarithmic Bacastow and Keeling (1973) formula, describing the flux of CO₂ into the land sink, is no longer sufficient. Consequently, the Lin assumption of Raupach (2013) might be violated for the land sink, implying that 2nd-order effects may be driving the negative slope of the sink rate that we document in this paper.

From Eq. (4) we obtain the relation

$${\displaystyle }{S}_t^L=k_{L,t}\cdot C_t,$$

which implies that the flux of CO₂ to the land sink is linear in C_t, where k_L,t would then be treated as a constant. On the other hand, a decreasing k_L,t implies that the efficiency with which the land sink absorbs CO₂ is decreasing, i.e. that the flux of CO₂ to the land sink is non-linear in C_t and that this non-linearity is such that the efficiency is decreasing. These statements are consistent with simulation results from climate cycle models (Friedlingstein et al., 2006). Here we illustrate mathematically how such non-linearities can arise.

The precise relationship between $S_t^L$ and C_t still eludes us, but (Bacastow and Keeling, 1973, p. 94) suggest that (in our notation)

$${\displaystyle }{S}_t^L\approx \mathit{\alpha }\log (\mathrm{1}+C_t/{\mathcal{C}}^{\mathrm{0}}),$$

where α is a constant and 𝒞⁰=591.30 GtC is the amount of CO₂ in the atmosphere in pre-industrial times. Using this function, we can write a 2nd-order Taylor expansion:

$${\displaystyle }{S}_t^L\approx \mathit{\alpha }\log (\mathrm{1}+C_t/{\mathcal{C}}^{\mathrm{0}})\approx \mathit{\alpha }{\displaystyle \frac{C_t}{{\mathcal{C}}^{\mathrm{0}}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\alpha }{\left({\displaystyle \frac{C_t}{{\mathcal{C}}^{\mathrm{0}}}}\right)}^{\mathrm{2}}.$$

Thus, if 𝒞⁰ is large compared to C_t, this implies that the squared term in the above equation is small and thus that a linear specification between $S_t^L$ and C_t is reasonable. However, once C_t becomes large compared to 𝒞⁰, this shows that the estimated sink rate will be found to be decreasing. To see this, use the Taylor expansion to write

$${\displaystyle }{S}_t^L{\displaystyle }\approx k_{L,t}{C}_t,$$

where

$${\displaystyle }{k}_{L,t}={\displaystyle \frac{\mathit{\alpha }}{{\mathcal{C}}^{\mathrm{0}}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\alpha }}{{\mathcal{C}}^{\mathrm{0}}}}{\displaystyle \frac{C_t}{{\mathcal{C}}^{\mathrm{0}}}},$$

is decreasing in C_t. In our data, we have C₁₉₅₉≈80 GtC and C₂₀₁₆≈267 GtC, resulting in $C_{\mathrm{1959}}/{\mathcal{C}}^{\mathrm{0}}\approx \mathrm{14}\mathit{\%}$ and $C_{\mathrm{2016}}/{\mathcal{C}}^{\mathrm{0}}\approx \mathrm{45}\mathit{\%}$. In other words, the linear approximation to the Bacastow and Keeling model of the land sink flux might have been reasonable in the past, since $C_{\mathrm{1959}}/{\mathcal{C}}^{\mathrm{0}}\approx \mathrm{14}\mathit{\%}$, but is likely misspecified in the present, since $C_{\mathrm{2016}}/{\mathcal{C}}^{\mathrm{0}}\approx \mathrm{45}\mathit{\%}$. That is, if this model is accurate, then a decreasing (land) sink rate indicates that we have entered a regime of atmospheric CO₂ concentrations, where the linear approximation breaks down and higher-order terms should be taken into account.