The declining uptake rate of atmospheric CO2 by land and ocean sinks

. Through 1959–2012, an airborne fraction (AF) of 0.44 of total anthropogenic CO 2 emissions remained in the atmosphere, with the rest being taken up by land and ocean CO 2 sinks. Understanding of this uptake is critical because it greatly alleviates the emissions reductions required for climate mitigation, and also reduces the risks and dam-ages that adaptation has to embrace. An observable quantity that reﬂects sink properties more directly than the AF is the CO 2 sink rate ( k S ), the combined land–ocean CO 2 sink ﬂux per unit excess atmospheric CO 2 above preindustrial levels. Here we show from observations that k S declined over 1959–2012 by a factor of about 1 / 3, implying that CO 2 sinks increased more slowly than excess CO 2 . Using a carbon–climate model, we attribute the decline in k S to four mechanisms: slower-than-exponential CO 2 emissions growth ( ∼ 35 % of the trend), volcanic eruptions ( ∼ 25 %), sink responses to climate change ( ∼ 20 %), and nonlinear responses to increasing CO 2 , mainly oceanic ( ∼ 20 %). The ﬁrst of these mechanisms is associated purely with the trajectory of extrinsic forcing, and the last two with intrinsic, feedback responses of sink processes to changes in climate and atmospheric CO 2 . Our results suggest that the effects of these intrinsic, nonlinear responses are already detectable in the global carbon cycle. Although continuing future de-creases in k S will occur under all plausible CO 2 emission scenarios, the rate of decline varies between scenarios in non-intuitive ways because extrinsic and intrinsic mechanisms respond in opposite ways to changes in emissions: extrinsic mechanisms cause k S to decline more strongly with increasing mitigation, while intrinsic mechanisms cause k S to decline more strongly under high-emission, low-mitigation scenarios as the carbon–climate system is perturbed further from a near-linear regime.


Introduction
The properties of natural land and ocean CO 2 sinks have major implications both for climate mitigation goals and for adaptive responses. The CO 2 airborne fraction (AF, the fraction of total anthropogenic CO 2 emissions from fossil fuels and net land use change that accumulates in the atmosphere) determines the fraction of emissions that contribute to rising atmospheric CO 2 concentrations, with the remainder (the sink fraction, SF = 1−AF) being absorbed by land and ocean sinks.
Since the commencement of high-quality atmospheric CO 2 measurements in 1958, the AF has averaged about 0.44 (Canadell et al., 2007;Knorr, 2009;Le Quéré et al., 2009;Ballantyne et al., 2012), with significant interannual variability (Keeling and Revelle, 1985). This fact is one of the most important attributes of the contemporary carbon cycle, with major policy implications both for the climate mitigation challenge and also for adaptation to climate change. The natural CO 2 sinks that absorb more than half of all anthropogenically emitted CO 2 represent a massive ecosystem service to humankind (Millennium Ecosystem Assessment, 2005), with implications that are directly quantified by the AF (Raupach et al., 2008). The mean AF and its possible past and future trends are therefore important, for both biophysical and policy reasons.
The basic reason for the approximate past constancy of the AF is well known: a constant or zero-trend AF would be expected under a "LinExp" idealisation of the carbon cycle, in which land and ocean CO 2 sinks increase linearly with excess CO 2 above preindustrial concentrations (assumption "Lin") and total anthropogenic CO 2 emissions increase exponentially ("Exp") (Bacastow and Keeling, 1979;Hofmann et al., 2009;Gloor et al., 2010;. This idealisation is a reasonable first approximation to the past behaviour of the carbon cycle, where total CO 2 emissions (both annual and cumulative) have increased roughly exponentially for more than a century, and sinks have increased roughly linearly with excess CO 2 .
However, the LinExp idealisation -despite its utility in explaining the observed approximate constancy of the AF -is imperfect even for the past, and is likely to become more so in the future. Several analyses (Canadell et al., 2007;Raupach et al., 2008;Le Quéré et al., 2009) have detected a small increasing trend in the AF since 1958 at a mean relative growth rate around 0.2 to 0.3 % year −1 , with significance (probability of positive trend) in the range 0.8 to 0.9. Methodological issues have been raised to question this result, concerning trend detection methods (Knorr, 2009), data (Ballantyne et al., 2012;Francey et al., 2010) and uncertainty analyses (Ballantyne et al., 2012). Nevertheless, multiple studies find results for the magnitude and significance of the AF trend that are in approximate agreement when consistent definitions are used (Canadell et al., 2007;Knorr, 2009;Le Quéré et al., 2009;Ballantyne et al., 2012).
Although the AF and its trend are important metrics of the behaviour of the carbon cycle with direct policy implications, they do not provide unambiguous information about the behaviour or "efficiency" of CO 2 sinks. The main reason is that AF mixes information about sinks and anthropogenic emissions, since AF = 1 − sinks/emissions. For example, under the LinExp idealisation, an AF trend appears when emissions increase non-exponentially, even when sinks are linear in excess CO 2 and thus have a constant efficiency .
Recognising these issues with the interpretation of the AF, we here provide an attribution of observed CO 2 sink behaviour over the last 50 years by using a novel observable diagnostic, the CO 2 sink rate. This is a more direct measure of sink efficiency than the AF and offers complementary insights into carbon cycle behaviour, although the two quantities are related and are both obtained from the same observations. We compare the observed past behaviours, likely future trajectories and diagnostic properties of the sink rate and the AF.

CO 2 mass balance and airborne fraction
The atmospheric CO 2 mass balance is or in more compact form Here c A = 2.127 [CO 2 ] − [CO 2 ] q is the excess CO 2 in Pg C (where [CO 2 ] is the CO 2 mixing ratio in ppm, and [CO 2 ] q = 278 ppm is [CO 2 ] at preindustrial equilibrium); a prime denotes differentiation with respect to time (t), so that c A = dc A /dt is the atmospheric CO 2 accumulation (in Pg C year −1 ); f E = f Foss + f LUC is the total CO 2 emission flux, the sum of emissions from fossil fuels and other industry (f Foss ) and from net land use change (f LUC ); and f ↓S = −f L − f M is the total (land plus ocean) CO 2 sink flux, the negative sum of the land-air (f L ) and ocean-air (f M , marine) exchange fluxes. All fluxes are in units of Pg C year −1 and are positive upward (surface to atmosphere), except for f ↓S , which is positive downward to denote a CO 2 sink. The airborne and sink fractions are the dimensionless quantities The AF is often alternatively defined as an "apparent airborne fraction" c A /f Foss (Oeschger et al., 1980); the definition used here allows the anthropogenic contribution to CO 2 growth from net land use change to be distinguished from the terrestrial carbon sink (Raupach et al., 2008;Le Quéré et al., 2009;Gasser and Ciais, 2013).

CO 2 sink rate
The CO 2 uptake rate by land and ocean sinks (k S , henceforth called the CO 2 sink rate) is the combined land-ocean CO 2 sink flux (f ↓S ) per unit mass of excess atmospheric CO 2 above preindustrial concentrations (c A ). It is defined (Raupach 2013) by and has dimension 1/time. This definition has two simple (and related) physical interpretations: k S is the sink strength per unit excess CO 2 , and equivalently the instantaneous fractional rate of decrease in excess CO 2 caused by sinks alone. Both lead to the interpretation of k S as a measure of "sink efficiency" . Several important properties follow from the definition of k S . First, Eqs. (2) and (5) imply that Therefore, k S (like the AF and SF) can be readily observed using basic data on global CO 2 emissions and concentrations. The simple relationship between k S , the SF and the AF is Second, the definition of k S can also be written as with k L = −f L /c A and k M = −f M /c A . Thus, k S can be split into additive components k L and k M for land and ocean sinks, respectively, as for the LF and OF (Eq. 4). Similar decompositions are possible for the regional sub-components of the land and ocean sinks. Third, k S depends directly on the sink flux and only indirectly (and weakly) on emissions trajectories through their effect on excess CO 2 . By contrast, the AF is directly affected as much by a change in emissions as a change in sinks . Therefore k S is a more appropriate diagnostic for sink properties than the AF.
Fourth, a trend in k S indicates a difference in the relative growth rates of sinks and excess CO 2 . The relative growth rate (RGR) of a quantity is its absolute growth rate normalised by its mean, with dimension 1/time; thus, the RGR of a time series X(t) is RGR(X) = d(ln X)/dt ≈ X −1 dX/dt , where angle brackets denote expected values. The evaluation and properties of RGRs are summarised in Appendix A; one important property is that the RGR of a product (or quotient) is the sum (or difference) of the RGRs of its factors (Eq. A3). Combining this with Eq. (5), the relative growth rate of k S can be written as Thus, k S increases (has a positive RGR) when the RGR for sinks exceeds that for excess CO 2 , and vice versa. Fifth, k S constitutes an observable weighted mean of the multiple timescales governing the global carbon cycle, describing their composite effect at any one time on excess atmospheric CO 2 . It is well known that there is no single lifetime for atmospheric CO 2 , because the carbon cycle includes multiple processes with timescales from days to millennia (Archer et al., 2009). Describing these processes is a fundamental challenge for carbon cycle modelling. A linearised, multi-pool carbon cycle model is equivalent to a pulse response function for atmospheric CO 2 (the airborne fraction after time t of a pulse of CO 2 into the atmosphere) of the form of a sum of exponentials: G(t) = a m exp(−λ m t), where the sum is over a set of modes m with turnover rates λ m and weights a m Joos et al., 2013;. The modes are a set of independent carbon pools z m (t), superpositions of physical carbon pools, that sum to the atmospheric excess carbon c A (t). It can be shown (see Appendix B, Eq. B8) that where b m is the fraction of c A in mode m, summing over m to 1. Thus, k S is a weighted sum of the turnover rates λ m , where the weights b m are time-dependent in linearised, pulse-response-function models of the carbon cycle, and the rates λ m are also time-dependent in nonlinear models. Equation (10) shows how k S aggregates the effects of multiple processes with different rates to determine the net drawdown rate of atmospheric CO 2 by sinks at any particular time.
Sixth, under the LinExp idealisation defined in Sect. 1, both k S and the AF are constant in time also Bacastow and Keeling, 1979, for the AF only). Conversely, neither k S nor the AF are constant if CO 2 sinks are nonlinear in excess CO 2 (departure from "Lin") or emissions are non-exponential (departure from "Exp").

Estimation of trends
Monthly trajectories for AF and k S from January 1959 to December 2012 (henceforth 1959.0-2013.0) are shown in Fig. 1. These were obtained from collated data (Le  on CO 2 emissions from fossil fuels (f Foss ) and net land use change (f LUC ), together with global atmospheric CO 2 concentrations; see Appendix C1 for details and references to primary sources. To quantify trends in AF and k S , we used several different data treatments (Appendix C2) and trend estimation methods (Appendix C3). Our measure of trend is the relative growth rate (Appendix A).
For AF, our best trend estimate ( Fig. 2 and Table 1) is RGR (AF) = 0.24 ± 0.20 % year −1 (±1σ , P = 0.89) about a mean AF of 0.44 over 1959.0-2013.0, where ±1σ denotes a 1-standard-deviation confidence interval, and the significance (P ) is the probability of positive trend. Both the trend and its significance are comparable with earlier studies cited in the Introduction, when consistent definitions are used; in particular, the statistical significance of the AF trend is found by all studies (including this one) to be less than 95 %, between "likely" and "very likely" in the standard terminology of the Intergovernmental Panel on Climate Change (IPCC, 2007).
For k S , the best trend estimate ( Fig. 2 and Table 2) is RGR (k S ) = −0.91 ± 0.17 % year −1 (±1σ , P > 0.999 for negative trend), about a mean k S of 0.028 (= 1/36) year −1 . The observed decreasing trend in k S is statistically robust and "virtually certain" in IPCC terminology, in contrast with the AF trend. The above uncertainty estimates for trends in AF and k S reflect variability associated with CO 2 growth rate, but not uncertainties in CO 2 emissions from fossil fuels (f Foss ) and net land use change (f LUC ). As described in Appendix D, these uncertainties were assessed by repeating the estimation of RGR (AF) and RGR (k S ) with 3 alternative f Foss trajectories (Francey et al., 2010;Gregg et al., 2008;Guan et al., 2012) (Fig. D1) and 11 alternative f LUC trajectories (Le Quéré et al., 2009) (Fig. D3). The resulting trend estimates are statistically indistinguishable from our best estimates.
The RGR of a time series X(t) > 0 can be estimated in two ways, as RGR(X) = d(ln X)/dt and RGR(X) ≈ X −1 dX/dt (Appendix A, Eqs. A1 and A2, respectively). The former definition is the more fundamental because it precisely preserves identities for RGRs of products and quotients (Eq. A3). However it is not usable in practice if the series X(t) includes any negative values (for which ln X is undefined), so the latter approximate definition must be   Table 1 caption.
used instead, even though the resulting RGR estimates do not exactly satisfy Eq. (A3). The RGR estimates in Tables 1  and 2 are obtained with Eq. (A2) because the input monthly time series change sign. Later we also use estimates from Eq. (A1), in circumstances where it is important that RGR product and quotient identities be satisfied. In practice the difference between the two RGR estimates is less than the statistical uncertainty in either.

Approach and model
We attribute trends in AF and k S by using a nonlinear carbon-climate model that approximately reproduces observed trends in AF and k S in its full form. By progressively simplifying the model to eventually reach the LinExp idealisation in which all trends are zero, the contributions of different factors to observed trends can be identified.
In general, it must be noted that attribution of an observed effect from multiple processes to individual process contributions is necessarily a modelling exercise (UNFCCC, 2002), with results that are model-dependent and not directly verifiable by observations unless the system can be manipulated experimentally. However, the present approach of progressively removing processes to reach a known analytic approximation has the advantage that two points in a "process space" are well characterised: the real world (which needs to be approximately reproduced by the model for any attribution exercise to be effective) and the simple analytic idealisation. Our choice of model for this exercise is determined by the requirements that (1) it can approximately reproduce observed trends AF and k S , and (2) it can be reduced to the LinExp idealisation by formal linearisation.
The model is the Simple Carbon-Climate Model (SCCM), a globally aggregated model of the carbon-climate system (Harman et al., 2011;Raupach et al., 2011). Model state variables comprise one atmospheric CO 2 store, two land carbon stores, four perturbation carbon stores in the ocean, the atmospheric concentrations of four major non-CO 2 greenhouse gases (CH 4 , N 2 O and two representative halocarbons), and three perturbation global temperatures representing heat stores with different turnover rates . Radiative forcing by aerosols is incorporated via a simple parameterisation based on a proportionality with f Foss (t), with a time-dependent coefficient. Carbon in the ocean mixed layer is modelled using a pulse response function that emulates the mixing dynamics of www.biogeosciences.net/11/3453/2014/ Biogeosciences, 11, 3453-3475, 2014 several complex ocean circulation models (Joos et al., 1996). The model ocean-atmosphere CO 2 flux incorporates full, nonlinear ocean carbonate chemistry (Lewis and Wallace, 1998). The model terrestrial biosphere includes a nonlinear dependence of terrestrial net primary production (NPP) on CO 2 concentration to account for CO 2 fertilisation of plant growth, and a nonlinear dependence of heterotrophic respiration on temperature. The effect of volcanic activity on the terrestrial carbon cycle (Jones and Cox, 2001) is included through an enhancement factor for terrestrial NPP that is proportional to a global volcanic aerosol index (Ammann et al., 2003), tested using recent major eruptions. SCCM does not resolve interannual variability associated with short-term climate fluctuations, regionally specific processes, and climate effects on the carbon cycle beyond those captured by a response to global temperature. In exchange for these simplifications, an important benefit for this work is that SCCM can be linearised analytically , allowing linearisation to be included explicitly as a simplifying step.

Model-data comparisons
SCCM satisfactorily reproduces observed trends in AF and k S over 1959.0-2013.0, as shown in Fig. 3 with a comparison between model predictions (red bars) and RGR estimates using Eqs. (A1) and (A2) (black and grey bars, respectively). Of these, the like-with-like comparison is between red and black bars, for which RGRs are calculated for both the model and observations (respectively) using annual data with Eq. (A1). The grey bars show the best RGR estimates (with uncertainties) from Tables 1 and 2, calculated with Eq. (A2); the difference between black and grey bars does not exceed ±1σ confidence intervals.
Comparisons of model predictions against observed time series of CO 2 , temperature, AF and k S over the period of high-quality CO 2 observations from 1959 onward indicate satisfactory performance for the purpose of attributing trends over this period (Figs. 4 and 5, right panels). In particular, the model reproduces the observed perturbations in AF and k S due to major volcanic eruptions (indicated by dots in Fig. 5). However, the model does not reproduce interannual climate variability related to El Niño-Southern Oscillation (ENSO) and other interannual climate modes. Also, model performance against data prior to 1959 is weaker than after 1959, for reasons including data quality, lack of account for interdecadal variability in air-ocean CO 2 and heat exchanges, and lack of incorporation of the full range of forcing factors into the model. Figure 3 shows the effects on the modelled trends in AF and k S of successive simplification by removing processes from the model, while leaving all model parameters unchanged.  Note that RGR (k S ) is negative and is plotted with reversed sign. Black bars show observed trends estimated from annual data using Eq. (A1) in Appendix A; these observed trend estimates are directly comparable with model estimates, which were obtained in the same way. Grey bars with ±1σ confidence intervals are the best trend estimates from Tables 1 and 2, obtained with Eq. (A2). There are small differences between the estimates (within ±1σ confidence intervals). Reasons for the use of these two different estimation methods are given in the text and in Appendix A.

Process attributions
The first simplification (V1 to V2, where V1 is the full model) is linearisation of the model carbon cycle, using the tangent-linear form of SCCM. This removes all nonlinear dependences of CO 2 fluxes and radiative forcing on carbon stores and temperatures, but retains linearised interactions among these quantities. The result is a reduction in the magnitude of the k S trend by ∼ 20 % (noting that RGR (k S ) is negative), most of the reduction being due to removal of nonlinearities associated with the dependence of ocean-air CO 2 exchange on atmospheric CO 2 .
The next simplification (V2 to V3) is carbon-climate decoupling, performed by removing all dependences of CO 2 fluxes on temperature through terrestrial NPP, heterotrophic respiration and ocean chemistry (recalling that linearised versions of these interactions were retained in the step from V1 to V2). This simplification also removes all effects of non-CO 2 gases on the carbon cycle, since these are mediated entirely by temperature in this model. This step reduces the magnitude of RGR (k S ) by another ∼ 20 % of its full-model value.
The third simplification (V3 to V4) is removal of the effects of volcanism on terrestrial NPP. This causes another ∼ 25 % reduction in the magnitude of RGR (k S ).
The last simplification (V4 to V5) is replacement of real total CO 2 emissions (f Foss + f LUC ), which depart from  Figure 4. Total CO 2 emissions (f E , top row) and SCCM predictions for CO 2 concentration (middle row) and temperature (bottom row), with analytic scenarios for future emissions of CO 2 and non-CO 2 gases (CH 4 , N 2 O, CFCs) such that the all-time cumulative total CO 2 emission Q takes values from 1000 to 3000 Pg C. Scenarios and model details (including treatment of aerosols) are given in . Left panels show plots against time from 1800 to 2200; right panels zoom in to the period 1950-2020 to compare model with data. This figure is a variation with added detail of Fig. 6 in . exponential growth , by an exponential trajectory with the same mean growth rate over 1850-2011. This removes the remaining ∼ 35 % of the magnitude of RGR (k S ). After all four simplification steps, the k S trend is reduced to zero in the model, consistent with the theoretical requirement of the LinExp idealisation.
We have repeated this progressive model simplification experiment with different orderings for process removal, finding that the above result is independent of ordering to a very good approximation. Not surprisingly, model simplification causes the agreement between model and observations to weaken progressively as processes are removed (Figs. 6 and 7).
The sequence of effects of progressive model simplification is not as simple for the AF trend as for k S (Fig. 3). For RGR (AF), the largest single change is brought about by removal of volcanic effects: this step alone eliminates the observed positive trend in AF, in accord with other recent findings . In Sect. 5.3 we investigate the reasons for the different responses of RGR (k S ) and RGR (AF) to progressive model simplification, and therefore the different attributions of the observed trends to processes.

Extrinsic and intrinsic mechanisms
We have attributed the decline in k S to four mechanisms. One of these -departure of emissions from exponential growthis "extrinsic", arising from the trajectory of external anthropogenic forcing of the carbon-climate system. Two othersnonlinear carbon cycle responses to CO 2 and carbon-climate coupling -are "intrinsic", arising from process feedbacks in the system. Volcanic effects are both extrinsic and intrinsic, involving feedbacks triggered by non-anthropogenic forcing from volcanic aerosols.
The primary extrinsic mechanism operates thus: when CO 2 emissions increase more slowly than exponentially, the fast-response, low-capacity modes of the carbon cycle saturate more rapidly than slow modes, so the weights b m in Eq. (10) decrease with time for faster modes and reciprocally increase for slower modes, causing k S to decrease. This mechanism is associated with sink capacities through the effect of CO 2 emissions trajectory on the distribution of carbon among the ocean, land and atmospheric stores. It can be described purely by linear theory .
In contrast, the primary intrinsic mechanisms arise from internal feedbacks. Many (though not all) of these are fundamentally nonlinear: prime examples are the dependences of ocean-atmosphere CO 2 fluxes and terrestrial NPP on CO 2 ,   (11) and (12), respectively. All growth rates are evaluated with Eq. (A1) and annual data to ensure that the product and quotient rules in Eq. (A3) are satisfied exactly.

Sign
Term RGR 23 and the dependence of heterotrophic respiration on temperature. These feedbacks have the net effect of decreasing the turnover rates λ m in Eq. (10) with increasing CO 2 and temperature, hence decreasing k S . Figure 4 (left panels) shows the trajectories of total CO 2 emissions, excess CO 2 and temperature under a set of analytically specified future emission scenarios for all modelled forcing agents; see Appendix B of  for details. The scenarios are characterised by the parameter Q, the cumulative total CO 2 emission (f E = f Foss + f LUC ) integrated from preindustrial times to the far future when emissions decline to zero. This ranges from high emissions (Q = 3000 Pg C) to strong mitigation (Q = 1000 Pg C). Also shown in Fig. 4 is a "StopEmission" projection for the case where all anthropogenic emissions of CO 2 and other forcing agents are stopped instantly at the present time, taken as 2013.0 for these calculations . The temperature projection for this scenario shows a rapid rise in temperature of about 0.5 K as the net cooling radiative forcing from aerosols is switched off suddenly, followed by a slow decline.

Future behaviours of AF and k S
The corresponding future behaviour of AF and k S is shown in Fig. 5 (left panels). Under a high-emissions scenario the AF remains close to its present value for a century or more into the future, while under a strong-mitigation scenario the AF declines fairly rapidly, becoming negative after the time of peak CO 2 . For the StopEmission scenario the AF is undefined from the present time onward.
There is continuing decline in k S under all scenarios. For the StopEmission scenario this decline is initially very rapid, with k S decreasing to less than half of its present value over a few years as fast modes saturate. Among the scenarios with Q ranging from 1000 to 3000 Pg C, the future rate of decline in k S varies surprisingly little. This occurs because extrinsic and intrinsic mechanisms respond in opposite ways to changes in emissions: extrinsic mechanisms cause k S to decline more strongly with increasing mitigation, as emissions trajectories fall progressively further below exponential growth. In contrast, intrinsic mechanisms cause k S to decline more strongly under high-emission, low-mitigation scenarios as the carbon-climate system is perturbed further from a near-linear regime and rates for individual sink processes decrease. The net result of these opposing influences is that our projected future values (in 2100) of the composite drawdown timescale 1/k S range from ∼ 120 to ∼ 180 year, for scenarios from emissions-intensive to strong-mitigation (Fig. 5).

Implications of differing trends in AF and k S
The proportional effects of the four process-removal steps (Fig. 3, Sect. 4.3) are not the same for RGR (AF) as for RGR (k S ), because these two growth rates depend in different ways upon changes in carbon-cycle fluxes and stores. From Eq. (7), RGR (k S ) can be written as showing that the growth rate of k S is a linear combination of the growth rates for SF (= 1 − AF), total emissions (f E ) and excess CO 2 (c A ). The equivalent expression for RGR (AF) is Thus, the growth rate of the AF is a linear combination of the growth rates of atmospheric CO 2 accumulation (c A , the time derivative of the excess CO 2 ) and emissions. The AF increases when CO 2 accumulation grows faster than emissions, and vice versa. Table 3 shows the contributions to RGR (k S ) and RGR (AF) of the terms in Eqs. (11) and (12), respectively. All growth rates are evaluated using Eq. (A1) to ensure that the product and quotient rules in Eq. (A3) are satisfied exactly. Over 1959.0-2013.0, RGR (k S ) was negative mainly because the growth rate of excess CO 2 c A /c A was significantly larger than the growth rate of total emissions, with a smaller contribution from the growth rate of SF. A positive growth in AF occurred because the growth rate of the CO 2 accumulation c A /c A exceeded the growth rate of emissions. These different dependencies indicate that there is no direct relationship between the growth rates of AF and k S , so it is not surprising that the process contributions to the two growth rates are very different (Table 3).
This discussion highlights the different insights obtained from absolute and relative growth rates. CO 2 sinks have unquestionably increased in absolute magnitude since 1959 (Ballantyne et al., 2012); recent work  has focussed on quantifying this absolute trend in units of Pg C year −2 . However, the increase in sinks (f ↓S ) has been  accompanied by increases in total emissions (f E ) and atmospheric accumulation (c A ) (the other terms in the atmospheric CO 2 mass balance, Eq. 2), and also by continuing growth in the excess CO 2 concentration itself (c A ) (Le Quéré et al., 2009; Fig. C1). Therefore, it is important to investigate not only the absolute growth rate in sinks but also the growth rates of these quantities relative to each other, or which quantities are "winning the race". Trends in AF, SF and k S answer this question, with the help of Eqs. (9), (11) and (12). Over 1959.0-2013.0, the positive sign of RGR (AF) indicates that c A grew faster than f E , the negative sign of RGR (SF) indicates that f ↓S grew slightly slower than f E , and the negative sign of RGR (k S ) indicates that f ↓S grew slower than excess CO 2 (c A ). This provides a simple rationale for the significance of the relative growth rates.

Model-independent and model-dependent findings
Many of the findings of this work are wholly or nearly independent of the particular simple model used here (SCCM); rather, they are based on observations or simple analytic inferences. The past decline in k S follows from observations of CO 2 emissions and accumulation. Attribution of a significant fraction of this decline to extrinsic mechanisms (associated with the effect of emissions trajectory on the distribution of carbon among stores, and thence on sink capacities) is based on robust linear theory, effectively a pulse-response-function description of the global carbon cycle . A continued future decline in k S from extrinsic mechanisms alone is expected from the same linear theory. Other aspects of our findings are model-dependent, including the precise fractional attributions of the decline in k S among mechanisms (Fig. 3), which may change as more sophisticated carbon-cycle models are brought to bear on the attribution problem. However, it is important to recognise that any model used in this way, whether simple or complex, must first pass two basic tests: that it includes the processes to be attributed, and that it can reproduce observations well enough for attribution to be possible. These tests are more fundamental than the level of complexity of the model.
As an illustration of the challenge, Figs. 8 and 9 (respectively for AF and k S ) compare the mean and relative growth rate over 1959.0-2013.0 between data, SCCM and the 11 models in the C4MIP intercomparison , in both uncoupled and coupled modes. For the growth rate of AF, 7 out of 11 C4MIP models in coupled mode predict the wrong sign (negative rather than positive as observed), and for k S , the negative growth rate is underestimated by all C4MIP models in coupled mode. For the AF, a similar comparison has been presented previously (Le Quéré et al., 2009); the comparison for k S is presented here for the first time. The reasons for these discrepancies may partly lie with the fact that the C4MIP protocol did not incorporate volcanism, which has only recently been found to have a large influence on carbon-cycle trends over decades .

Conclusions
The implications of this work can be summarised as follows. First, the trajectories of AF and k S provide different insights into the behaviour of the carbon cycle: trends in the AF indicate differences in the relative growth rates of excess CO 2 accumulation and anthropogenic CO 2 emissions (Eq. 12), while trends in k S indicate differences in the relative growth rates of sinks and excess CO 2 concentration (Eq. 9). Immediate implications of the observed decline in k S over 1959.0-2013.0 are that CO 2 sinks increased more slowly than excess CO 2 , and that the sink efficiency (the sink strength per unit excess CO 2 ) decreased. Second, k S constitutes an observable weighted mean of the multiple rates λ m of processes controlling the global carbon cycle, describing their combined effect on excess atmospheric CO 2 through land and ocean sinks (Eq. 10). Over 1959.0-2013.0, the composite drawdown timescale 1/k S increased from ∼ 30 to ∼ 45 years, and is projected to increase further in future. Therefore the mix of carbon-cycle timescales contributing to drawdown of CO 2 by sinks has shifted observably towards longer scales.   in both uncoupled and coupled modes (orange and green bars, respectively). Lower panel: relative growth rate RGR(k S ), estimated using Eq. (A2). See Fig. 8 caption for note on consistency of SCCM trend estimates between Figs. 9 and 3. Third, we attribute the observed decline in k S to four mechanisms: slower-than-exponential CO 2 emissions growth (∼ 35 % of the trend); volcanic eruptions (∼ 25 %); responses of CO 2 sinks to climate change (∼ 20 %); and nonlinear responses to increasing CO 2 , mainly associated with ocean CO 2 sinks (∼ 20 %). The first of these mechanisms is "extrinsic" (associated with the trajectory of external forcing of the carbon-climate system and its effect on sinks through the distribution of carbon among ocean, land and atmospheric stores). The last two are "intrinsic" (associated with feedback responses of sink processes to changes in climate and atmospheric CO 2 ). Volcanic effects include both non-anthropogenic extrinsic forcing and intrinsic feedback responses.
Fourth, the observed decline in k S is projected to continue under all realistic emissions scenarios (Fig. 5). This implies that the shift of the mix of carbon cycle timescales towards longer scales, already evident in past observations, is projected to continue in future. By contrast, future trends in AF are much more strongly dependent on emissions scenarios, with the AF becoming negative under strong-mitigation scenarios. , 11, 3453-3475, 2014 www.biogeosciences.net/11/3453/2014/ Fifth, our model-based attribution suggests that the effects of intrinsic mechanisms (carbon-cycle responses to CO 2 and carbon-climate coupling) are already evident in the carbon cycle, together accounting for ∼ 40 % of the observed decline in k S over 1959.0-2013.0. These intrinsic mechanisms encapsulate the vulnerability of the carbon cycle to reinforcing system feedbacks. By comparison, the extrinsic, sinkcapacity mechanisms are much easier to describe and are captured by pulse-response-function models of the global carbon cycle. An important open question is how rapidly the intrinsic mechanisms and associated feedbacks will contribute to further decline in k S under various emission scenarios.

Biogeosciences
Finally, the approach of progressive model simplification used here can be applied to attribute trends in k S with other suitable models. While our attribution is necessarily restricted to processes resolved in the simple model used here, a more complex model could attribute trends to more finely resolved processes such as regional contributions to land and ocean sinks Sitch et al., 2013). The approach ensures that all contributions sum to the full model trend in k S . Such an attribution would show not only how different regions contribute to the global ecosystem service provided by land and ocean carbon sinks, as quantified by their additive contributions to the global sink flux or global sink rate k S , but also how these contributions are changing in different ways in response to both extrinsic (forcing) and intrinsic (feedback) influences. For a process X(t) with all X(t) > 0, the relative growth rate is where angle brackets denote expected values. However, time series for noisy processes X(t) often have values of both signs -for example, monthly time series of AF and k S . In this case, Eq. (A1) must be approximated as Equation (A1) is the more fundamental of the two definitions because it has the important property of automatically satisfying the following identities for relative growth rates of products, quotients and powers of any processes X(t) and Y (t): Relative growth rates calculated with Eq. (A2) do not automatically satisfy these identities. However, in practice for the the AF and k S , the difference between Eqs. (A1) and (A2) is less that the statistical uncertainty in RGR from either method. Equation (A2) is used with monthly data for the "bestestimate" RGR values in this paper, because most monthly series have values of both signs. For the attribution of contributions to RGR (k S ) and RGR (AF) (Fig. 3, Table 3), Eq. (A1) is used with annual data to ensure that the attributed contributions sum to the total trend.

as a weighted mean of turnover rates
Here it is shown that the sink rate k S is a weighted mean of the turnover rates contributing to a pulse response function for atmospheric CO 2 , following previous work (Raupach, 2013) with simplified notation. At a high level of generality, a linearised, multi-pool model of the carbon cycle is where c i (t) is the excess carbon (the perturbation above a preindustrial equilibrium state) in pool i, f i (t) is the anthropogenic carbon input into pool i, and K = K ij is a square transfer matrix describing inter-pool transfers. This is a coupled dynamical system that can be solved readily by method of normal modes. The approach is to transform the system to a new frame where the state variables c i (t) become "normal modes" satisfying independent, uncoupled equations. In the new frame, the atmospheric excess carbon pool c A can be written as a sum of rescaled normal modes z m (t): The modes z m (t) are linear superpositions of the excess carbon pools c i (t), governed by where λ m is the turnover rate for mode m, and a m is a weight (summing over m to 1) specifying the fraction of total emissions to the atmosphere (f E = f Foss + f LUC ) entering mode m. The rates λ m are the eigenvalues of the transfer matrix K.
The solution for c A (t) is then given by where is a pulse response function (PRF) for atmospheric CO 2 (the fraction of an instantaneous pulse of CO 2 into the atmosphere that remains airborne after time t) taking the form of a sum of decaying exponential terms with decay rates λ m . One of the decay rates is often taken as zero, so that G(t) = a 0 + a 1 exp(λ 1 t) + . . .
Summing Eq. (B3) over modes m, the excess atmospheric CO 2 is governed by From Eq.
(2), the atmospheric CO 2 budget (with the total CO 2 sink expressed using the definition of the sink rate k S , Eq. 5) is Equating the last terms in Eqs. (B6) and (B7), it follows that Hence k S is a weighted mean of the turnover rates λ m for different modes. The weights b m are the fractions of c A appearing in the modes m, and from Eq. (B2), these weights sum to 1. The weights b m depend on time in general, because the modes z m grow at different rates λ m . If emissions f E (t) were steady, then z m for faster modes with larger λ m would saturate to the equilibrium value f E(steady) /λ m more rapidly than z m for slower modes. This would cause b m to give progressively higher relative weight to slower modes as time advances, so that k S would decrease. In the case where emissions increase exponentially, k S is constant in time, like the AF. An exponentially increasing trajectory f E (t) is the only case leading to constant k S and AF   2. AF(m) is a simple, untreated monthly AF time series: AF(m) = ( c A / t)/(f Foss + f LUC ) with t = 1 month and discretisation to yield 12 month-centred estimates per year (at 2009 + 1/24, 2009 + 3/24, . . . , 2009 + 23/24). Linear interpolation between annual data points was used to estimate emissions at intervening times, and linear interpolation between monthly data points was used similarly for concentrations.
3. AF(m, s) is a version of the monthly series AF(m) with 15-month running-mean smoothing applied to time series of dc A /dt before calculation of the AF. This removes most high-frequency (faster than annual) variability (Raupach et al., 2008).
4. AF(m, n) is a monthly AF series without smoothing but with noise reduction by removal of the fluctuating component linearly correlated with El Niño-Southern Oscillation (ENSO) and volcanic aerosol indices. These together account for about half the variance in dc A /dt from fluctuations shorter than a decade (Raupach et al., 2008). Contributions to the other half of the variance include climate modes other than ENSO, nonlinear effects, and regionally specific effects.
5. AF(m, s, n) is a monthly AF series with both 15-month smoothing as in AF(m, s) and noise reduction as in AF(m, n). Results are insensitive to the order in which smoothing and noise reduction are applied.
For the CO 2 sink rate k S , similar data treatments were used. This yielded five series: k S (a), k S (m), k S (m, s), k S (m, n) and k S (m, s, n).

C3 Trend estimation methods
The four trend estimation methods were as follows.
1. Linear regression: simple least-squares linear regression overestimates the confidence in the estimated trend, yielding a spuriously low CI (confidence interval) and spuriously high P value (probability of positive trend) when a series is temporally autocorrelated, as for all our series.
2. Stochastic method: this method estimates the confidence interval for the trend with account for the temporal autocorrelation of the time series (Canadell et al., 2007;Le Quéré et al., 2007). For a time series X(t), steps are as follows: (a) the trend X T is found by conventional least-squares linear regression, yielding a trend line The lagged autocorrelation function of the residual (X − X T ) is fitted with an autoregressive (AR) model (Box et al., 1994). To represent the autocorrelation function for monthly data adequately, an AR model of order 20 is used, noting that a good fit to the autocorrelation function is desirable and does not translate into overfitting in the final result because of the stochastic nature of the method. (c) An ensemble of 1000 stochastic realisations of the data is generated with mean trend X T and residuals correlated as in the AR model. Members of this ensemble have a similar mean trend and autocorrelation function to the original series, but vary stochastically among realisations. (d) The trend for each member of this ensemble is determined by least-squares linear regression, yielding 1000 estimates of the trend (x 1 ). (e) The probability density function (PDF) of trend estimates x 1 is calculated, yielding trend statistics. The P value is the fraction of the 1000 estimates of x 1 that is greater than 0 for a positive trend, or less than 0 for a negative trend.
3. Bootstrap method: both regression and stochastic methods suffer from sensitivity to the choice of start and end times, yielding different results if start or end times are shifted by a few months. The bootstrap method overcomes this problem. An ensemble of time series is constructed by selecting continuous subseries from the original series X(t) with randomised start and end times, subject to the condition that the minimum record length of each ensemble member is at least a fraction f Bts of the complete record. The ensemble is constructed with replacement, so this is a "bootstrap" method. The members of the ensemble are not independent, but represent different possible realisations of the observational constraints determining the length of the series. The choice of f Bts is a compromise between the requirements that (a) ensemble members include most of the original series and (b) the variations in start and end times be enough to randomise their effects. The latter requirement demands that the omitted portions of the record typically encompass several integral timescales for the series X(t). We used f Bts = 0.8. The bootstrap method reduces the influence of the choice of start and end times and therefore yields an improved estimate of mean trend, but provides no estimate of uncertainty information (CI or P value) because the ensemble members are not independent.
4. Combined method: here both the stochastic and bootstrap methods are applied together. Steps are as follows: (a) an ensemble of continuous subseries with randomised start and end times is selected as in the bootstrap method. (b) The trend (x 1 ) for each subseries is determined by linear regression. (c) The lagged autocorrelation function for the entire series is found, as in the stochastic method. (d) Using this autocorrelation function and the trend for each subseries, a stochastic ensemble of series is generated. (e) The trend for each ensemble member is found by linear regression. (f) Statistics of the ensemble of trend estimates are found as in the stochastic method. The combined method provides our best estimates, because it combines the benefits of the bootstrap method for estimation of trend and the stochastic method for confidence interval.  (2010)

Appendix D: Implications of uncertainties in emissions
The uncertainty estimates for RGR (AF) and RGR (k S ) in Fig. 2 and Tables 1 and 2 reflect the stochastic variability associated with CO 2 growth rate, but not the uncertainty in data on CO 2 emissions from fossil fuels and other industrial processes (f Foss ) and from net land use change (f LUC ). These are assessed as follows. Uncertainty in CO 2 emissions from fossil fuels and other industrial processes (f Foss ): the uncertainty in f Foss is estimated as ±6 % (Andres et al., 2012;Marland, 2008). If this is random, the uncertainty propagated into RGR (AF) and RGR (k S ) is very small. However, some studies have suggested systematic or strongly temporally autocorrelated biases for some countries, notably an underestimate of up to 20 % in the late 1990s to early 2000s for China (Gregg et al., 2008). This would also be consistent with a suggested underestimate of global f Foss + f LUC for this period (Francey et al., 2010). Also, it has been suggested recently that there are significant uncertainties in Chinese emissions, particularly since 2005, from discrepancies between national and summed provincial accounts (Guan et al., 2012).
To assess the consequences of these possible revisions to f Foss , we computed RGR (AF) and RGR (k S ) using three alternative f Foss series (Fig. D1), in which f Foss was increased (1) between 1998 and 2003, (2) between 1993 and 2003, and (3) from 2000 onward using revised Chinese emissions based on provincial rather than national data (Guan et al., 2012). The resulting trends RGR (AF) and RGR (k S ) (Fig. D2), computed using data treatment (m, s, n) and the combined trend detection method, are slightly smaller in magnitude than the best estimates with primary f Foss data, but the differences are not statistically significant. Therefore, our conclusions are unaffected by any of the three possible revisions to f Foss (all of which are still speculative). Uncertainty in CO 2 emissions from net land use change (f LUC ): it is well known that uncertainty in f LUC is significant (Houghton, 2010(Houghton, , 2003  and three alternative f Foss trajectories shown in Fig. D1 (pale grey bars). Trends are estimated using Eq. (A2). All estimates are computed using data treatment (m, s, n) and the combined trend detection method, as for best estimates ( Fig. 2 and Tables 1 and 2). Error bars show ±1σ confidence intervals. For RGR (AF), all P values (probability of positive trend) exceed 0.69; for RGR (k S ), all P values (probability of negative trend) exceed 0.996. uncertainty in AF trend estimates (Le Quéré et al., 2009;Raupach et al., 2008). The error on all f LUC estimates is large, typically ±50 %. For estimation of both RGR (AF) and RGR (k S ), uncertainties arising from systematic biases in f LUC data (both in level and trend) are more important than uncorrelated random errors in annual estimates. The primary data used here (Fig. C1) Figure D3. Alternative trajectories for CO 2 emissions from net land use change (f LUC ): (black) primary data; (coloured dots) alternative trajectories described in Table D1. Upper and lower panels show "Recent Fall" and "Recent Const" extrapolations, respectively. on high-resolution remote sensing imagery, while latest estimates for Indonesia are based on data for 2003 and 2006, in contrast with 2005 estimates based on forecasts for this period. Recent studies in parts of the world that have dominated global deforestation inventories in past decades, including Brazil (Nepstad et al., 2009;Regalado, 2011) and Indonesia (Hansen et al., 2008), support the hypothesis that f LUC has declined significantly through 2000-2009.
To assess the implications of uncertainty in f LUC for trends AF and k S , we replaced the primary f LUC data with 11 alternative annual time series from other assessments ( Fig. D3 and Table D1). These alternative series are not independent, being based on just three sources of land cover data (FAO (FAO, 2010), SAGE (Ramankutty and Foley, 1999) and HYDE (Goldewijk, 2001)) and several carbon cycle models. Nevertheless, these series represent presently available estimates of global f LUC from numerous investigators. All alternative f LUC series end between 1990 and 2000, so each series was extrapolated in time by assuming either a linear decrease in f LUC from 2000 onward at 0.03 Pg C year −1 per year consistent with our primary f LUC data or a constant f LUC from the end of each alternative f LUC series (respectively denoted "Recent Fall" and "Recent Const" in Fig. D3).  Figure D4 (source: aaTrendEstimation.Results.V34.xls) Figure D4. Estimates of RGR (AF) using primary data for f LUC (blue bars) and 11 alternative f LUC trajectories (Table D1 and Fig. D3, pale grey bars). Trends are estimated using Eq. (A2). P values are shown in Table D1. All estimates are computed using data treatment (m, s, n) and the combined trend detection method, as for best estimates in Fig. 2 and Tables 1 and 2. Error bars show ±1σ confidence intervals. Upper and lower panels show RGR (AF) using "recent fall" and "recent const" extrapolations, respectively.
Because of the likely decline in f LUC since 2000, "Recent Fall" is the more likely scenario.
The resulting trends RGR (AF) (Fig. D4) and RGR (k S ) (Fig. D5), computed using data treatment (m, s, n) and the combined trend detection method, are slightly smaller in magnitude than the best estimates with primary f LUC data, but the differences are not statistically significant. This indicates that uncertainty in f LUC does not significantly affect our results.  (Table D1 and