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

Is the fraction of anthropogenically released CO2 that remains in the atmosphere (the airborne fraction) increasing? Is the rate at which the ocean and land sinks take up CO2 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 CO2 from the atmosphere is decreasing at approximately 0.54% yr−1.


Introduction
A part of the anthropogenically released CO 2 emitted to the atmosphere flows to the oceans (the ocean sink) and the terrestrial biosphere (the land sink).Approximately 45 % of released CO 2 stays in the atmosphere (the airborne fraction), while the two sinks take up approximately 24 % and 31 % of the CO 2 , 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 2 that ultimately remains in the atmosphere increases, and projections of future atmospheric CO 2 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 2 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 2 from the atmosphere is decreasing.Thus, analysing the behaviour of the sink rate can help predict the future uptake of CO 2 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 2 (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 2 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 Published by Copernicus Publications on behalf of the European Geosciences Union.
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 2 , and the CO 2 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 2 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.

The global carbon budget
The so-called global carbon budget is defined as where E ANT t is anthropogenically released CO 2 into the atmosphere, G t is growth of atmospheric CO 2 concentration, S O t is the flux of CO 2 from the atmosphere to the oceans (the ocean sink), and S L t is the flux of CO 2 from the atmosphere to the terrestrial biosphere (the land sink).In other words, Eq. (1) states that emissions of CO 2 should equal the fluxes of CO 2 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). 1 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 IM t for budget imbalance.Therefore, when considering actual data, the carbon budget is defined as (2) 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 IM t as part of an error term.The growth rate in atmospheric CO 2 data, G t , is from Dlugokencky and Tans (2018), the ocean sink data, S O t , are obtained from an ensemble of global biochemistry models, and the land sink data, S L t , 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 2 can be decomposed into two parts: where E FF t are emissions from fossil fuel burning, cement production, and gas flaring, while E LUC t are emissions from land-use change.Fossil fuel emissions, E FF t , are from Boden et al. (2018), while land-use change emissions, E LUC t , 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 2 in the atmosphere is constructed as where [CO 2 ] 1750 = 279 ppmv (parts per million by volume) and [CO 2 ] 1959 = 315.39ppmv are the concentrations of CO 2 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 preindustrial levels is given by the initial value in 1959 plus the cumulative sum of the growth in atmospheric CO 2 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 Biogeosciences, 16, 3651-3663, 2019 and the CO 2 sink rate as which is the flux of CO 2 from the atmosphere to the sinks (ocean plus land), normalized by the amount of CO 2 (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 respectively, with k 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 ANT t is the ratio of the growth of atmospheric CO 2 in period t to the amount of CO 2 emitted in period t.It is thus a measure of the fraction of emitted CO 2 that stays in the atmosphere.In contrast, the sink rate k S,t = (S O t + S L t )/C t is the ratio of the CO 2 flux in the sinks in period t to the total amount of CO 2 in the atmosphere (above pre-industrial levels).By writing S O t + S L t = k S,t C t , we can interpret the sink rate k S,t as the "efficiency" with which CO 2 flows from the atmosphere to the sinks, i.e. as the amount of CO 2 going into the sinks for an extra unit of CO 2 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.

Trend model specification
In this section, we consider several models for the datagenerating 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 where y t is a vector of observations at time t = 1, . .., 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 y t = (k O,t , k L,t ) , 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 E(x t |y 1 , . .., y n ; A, B, ξ , κ ), for t = 1, . .., 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 where β ∈ R is a fixed and unknown coefficient and η t is an iid Gaussian random variable with a mean of zero and variance σ 2 η .The solution to the difference equation in Eq. ( 6) is given as where T 1 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 1 , a deterministic linear term with slope β, and a random-walk component t−1 i=0 η t−i .Thus, T t can be interpreted as a longterm 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 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 σ 2 .
The local linear trend model can be cast in the state-space system Eq.( 5) where vectors and matrices are defined as www.biogeosciences.net/16/3651/2019/Biogeosciences, 16, 3651-3663, 2019 for t = 1, . .., 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.

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: t = where 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 zeromean and covariance stationary error process.
We consider our baseline local linear trend model of Sect. 3 for each of the objects, that is, for i = 1, 2, where the trend T (i) t is specified in Eq. ( 6) and with error (i) t .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 (i) t and the trend error term η (i) t , as well as the estimate of the slope parameter β, including the estimated standard error (SE) of β and the resulting t statistic, t-stat = β / SE( β).Based on these estimation results, we can formally test hypotheses of the type or, more relevantly, 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 0 in this case (p values 0.2711 and 0.4042 for the case of AF (1) t and AF (2) t , respectively).In other words, although the estimate β 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 y t − E(y t |y 1 , . .., y t−1 ; A, B, ξ , κ ), for t = 1, . .., 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 2 d , 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 2 , we have The reported diagnostic statistics and goodness of fit in Table 1 are satisfactory for the time series AF (1) t and AF (2) t .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 (2) t is well-described within our state-space framework, the extra error term ξ t = B IM t /E ANT t in AF (2) t , 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.
The state-space system allows both measures for the airborne fraction, AF (1) t and AF (2) t , 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 (1) t and T (2) t , both specified as Eq. ( 6), for the two time series.We have where the error terms t , for i = 1, 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 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 2 D , 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 χ 2 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 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.
model to be a good representation of the data.Furthermore, the slope is estimated to be positive in both cases (that is β > 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 AF (1) t and AF (2) t , respectively).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 A t , but with possibly different error terms (1) t and (2) t .In other words, we consider 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 β, 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).
5 Trend analysis of the CO 2 sink rate In this section, we analyse the CO 2 sink rate in the same way as the airborne fraction above.Analogously we can define two alternative versions of the sink rate: S,t = where now ξ t = B IM t /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 for i = 1, 2, where T (i) t 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 ξ t = B IM t /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 (2) S,t , respectively).We still consider a one-sided test as in Eq. ( 10), but now the relevant alternative hypothesis is H 1 : β < 0.
Analogously to the airborne fraction above, these data can be put in a joint uninformed system with two different trend components, and we have which can be compared with the model in Eq. ( 11).We report parameter estimates for the standard deviations σ (i) and σ (i) η , for i = 1, 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 2 D , and the Durbin-Watson (DW) test statistic for serial correlation.For details, see Table 1.In panel b, the trend coefficients (σ η and β) for AF (2) are the same as for AF (1) given the construction of the model (Eq.12).Finally, we consider the state-space system that imposes a common trend for both time series, T t , that is, 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 β = −0.00014,and this estimate is statistically significant: we reject the hypothesis H 0 : β = 0 in favour of H 1 : β < 0 at a 95 % confidence level (p value 0.0014).The mean of the sink rate (calculated using either data set k (1) S or k (2) S ) is 0.0258.It follows that we estimate the sink rate to be decreasing with approximately 0.00014/0.0258= 0.54 % yr −1 .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 reason- 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 2 D , 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.We report parameter estimates for the standard deviations σ (i) and σ (i) η , for i = 1, 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 2 D , and the Durbin-Watson (DW) test statistic for serial correlation; for details, see Table 1.In panel b, the trend coefficients (σ η and β) for k able is an interesting question beyond the scope of the present study.

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.
where the time series k O,t and k L,t are defined in Eq. ( 4), while the trend variables T O t and T L t 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 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 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 1 : β < 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 0 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, β = −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.

Discussion
Previous studies of the airborne fraction and the CO 2 sink rate have focused on detecting a single linear and deterministic trend in the data of the form a 0 + a 1 t, where a 0 and a 1 , are constants (Canadell et al., 2007a;Le Quéré et al., 2009;Knorr, 2009;Raupach et al., 2008Raupach et al., , 2014)).However, possible statistical difficulties in such analyses have been pointed out in Knorr (2009).For instance, a linear regression anal-ysis 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 σSlp > 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 ANT t might be autocorrelated.As shown in Tables 1-5, the diagnostic statis-Table 5. Analysis of ocean and land sink rates.
(a) Two trends and two observation series as in Eq. ( 16).

Parameter estimates
Correlation matrix ( ) Diagnostics We report parameter estimates for standard deviations σ (i) and σ (i) η , for i = 1, 2, 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 2 D 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).
tics 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 2 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 2 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 ANT t , 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 2 in the atmosphere (i.e.C t = G t ), which is not captured by the concurrent ratio G t /E ANT t .When studying the airborne fraction, it would perhaps be more reasonable to study an object taking this cumulative nature into account, e.g.G t / E ANT t = C t / E ANT t (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 O t and S L t are linear in concentrations C t ("Lin") and that emissions (E ANT t ) 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 2 concentrations: it is plausible that due to a rising level of CO 2 concentration, non-linear effects in the terrestrial CO 2 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(Cook, , 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 = 1959, 1960, . .., 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 2 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 ANT t .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.

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 2 and in the CO 2 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 2 from the atmosphere seems to be decreasing.The best estimate resulting from our state-space system is that the CO 2 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 −1 .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 2 .That is, although the land sink is itself increasing and thus continuing to take up a large part of anthropogenically emitted CO 2 , as also noted recently by, for example, Rayner et al. (2015), Keenan et al. (2016), andFerná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.
Data availability.The data used in this paper are available at the website of the Global Carbon Project (https://www.globalcarbonproject.org, Le Quéré et al., 2018).

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 2 may have risen to a point where a linear expansion of the logarithmic Bacastow and Keeling (1973) formula, describing the flux of CO 2 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 which implies that the flux of CO 2 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 2 is decreasing, i.e. that the flux of CO 2 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 L t and C t still eludes us, but (Bacastow and Keeling, 1973, p. 94) suggest that (in our notation) where α is a constant and C 0 = 591.30GtC is the amount of CO 2 in the atmosphere in pre-industrial times.Using this function, we can write a 2nd-order Taylor expansion: Thus, if C 0 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 L t and C t is reasonable.However, once C t becomes large compared to C 0 , this shows that the estimated sink rate will be found to be decreasing.To see this, use the Taylor expansion to write where . 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

Figure 1 .
Figure 1.Estimated trend T A t of the airborne fraction from the model 12.

S
given the construction of the model (Eq.15).

Figure 2 .
Figure 2.Estimated trend T S t of the CO 2 sink rate from Model (Eq.15).

Figure 3 .
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.

Table 2 .
The main difference from Table1is the inclusion of the estimated correlation matrix for ( The diagnostic test statistics are reasonable.In comparison with the univariate analysis, the goodness-of-fit values for R 2 d are slightly higher for the multivariate model.Hence we trust the t ).

Table 1 .
Univariate analysis of the airborne fraction.

Table 2 .
Multivariate analysis of the airborne fraction.Two individual trends as in Eq. (11).

Table 3 .
Univariate analysis of the CO 2 sink rate.

Table 4 .
Multivariate analysis of the CO 2 sink rate.
is decreasing in C t .In our data, we have C 1959 ≈ 80 GtC and C 2016 ≈ 267 GtC, resulting in C 1959 /C 0 ≈ 14 % and C 2016 /C 0 ≈ 45 %.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 1959 /C 0 ≈ 14 %, but is likely misspecified in the present, since C 2016 /C 0 ≈ 45 %.That is, if this model is accurate, then a decreasing (land) sink rate indicates that we have entered a regime of atmospheric CO 2 concentrations, where the linear approximation breaks down and higher-order terms should be taken into account.