Dynamic Responses of Atmospheric Carbon Dioxide Concentration to Global Temperature Changes between 1850 and 2010

Anthropogenic carbon dioxide (CO₂) emissions from fossil-fuel usage and land-use changes have been almost exponentially increasing since the Industrial Revolution (Fig. 1). Their accumulation in the atmosphere appears to be changing Earth's climate (IPCC, 2013), yet the full strength of anthropogenic CO₂ emissions for changing the climate has not been reached. Only 41%-45% of the CO₂ emitted between 1850 and 2010 remained in the atmosphere, while the rest was sequestered by lands and oceans (Jones and Cox, 2005; Canadell et al., 2007; Raupach et al., 2008; Knorr, 2009; also see Fig. 1 of this study). This largely constant ratio, generally referred to as the "airborne fraction" (denoted as "α" in this paper), was conventionally used to evaluate the efficiency of global carbon sinks in assimilating the extra CO₂ from the atmosphere (Jones and Cox, 2005; Canadell et al., 2007). However, recent studies indicate that the airborne fraction can be influenced by other factors and thus may not be an ideal indicator for monitoring changes in the carbon sink efficiency (Knorr, 2009; Gloor et al., 2010; Frölicher et al., 2013). The calculation of the airborne fraction also neglects important responses of the global carbon cycle to climate changes. Global surface temperature has increased by ∼1^°C since the beginning of the 20th century (Hansen et al., 1999; Brohan et al., 2006). Given the tight coupling between temperature and the carbon cycle, the warming alone can release a large amount of CO₂ from the land and the oceans into the atmosphere, redistributing carbon among these reservoirs (Keeling et al., 1995; Joos et al., 1999, 2001; Lenton, 2000; Friedlingstein et al., 2006; Boer and Arora, 2009, 2013; Rafelski et al., 2009; Frank et al., 2010; Willeit et al., 2014). Such climatic impacts need to be accounted for in order to diagnose and monitor changes in the global carbon cycle.

Figure 1.  Time series of global anthropogenic CO$_2$ emissions (red line), atmospheric CO$_2$ concentrations (green line), and the anomalous CO$_2$ fluxes induced by warming surface temperatures (gray shading) between 1850 and 2010. The top panel indicates the accumulated CO$_2$ fluxes or the total concentration changes, while the bottom panel shows them at annual steps. The thick and thin lines indicate long-term and interannual variations of the time series, respectively. The mathematical symbols are the same as in Eq. (1) and explained in the text. In both annual and accumulative cases, CO$_2$ emissions largely increase as an exponential function of time, while changes in the atmospheric CO$_2$ concentration are proportional to the corresponding emissions by a factor of about 0.41-0.45.

Previous studies have noticed that the effects of temperature on atmospheric CO₂ vary at different time scales, ranging from 1-2 ppm ^°C^-1 at the scale of years to 10-20 ppm ^°C^-1 over centuries or millennia (Woodwell et al., 1998). The long-term temperature sensitivity of atmospheric CO₂ is usually estimated by examining the correlations between ice-core CO₂ measurements and reconstructed paleoclimatic records before the industrial era, which comprise negligible anthropogenic influence (Woodwell et al., 1998; Scheffer et al., 2006; Lüthi et al., 2008; Lemoine, 2010; Frank et al., 2010). The short-term sensitivity is estimated from contemporary observations of temperature and atmospheric CO₂ (Keeling et al., 1995; Adams and Piovesan, 2005; Wang et al., 2013), which provide the necessary high temporal resolutions for the task. Because the contemporary atmospheric CO₂ records are strongly influenced by growing anthropogenic carbon emissions, the trends in the CO₂ and temperature time series often need to be removed before the data are used in correlation analysis (Wang et al., 2013). This "de-trending" process, however, inhibits the estimation of the long-term relationship between the two fields by conventional correlation-based methodology.

The timescale-dependent temperature sensitivity of atmospheric CO₂ concentration has also been suggested by modeling studies. By driving an earth system model of intermediate complexity with idealized periodic forcing, (Willeit et al., 2014) found that the simulated temperature sensitivities of atmospheric CO₂ range from 1-4 ppm ^°C^-1 at interannual and decadal scales to 5-12 ppm ^°C^-1 at centennial scales, which are generally consistent with observation-based estimates (e.g., Woodwell et al., 1998). However, model results are prone to uncertainties induced by different representations of carbon cycle processes and different initial state conditions in the simulations. For instance, the centennial temperature sensitivities of atmospheric CO₂ estimated from an ensemble of CMIP5 (Coupled Model Intercomparison Project Phase 5) models show a wide range from 7 to 30 ppm ^°C^-1 (Arora et al., 2013; Willeit et al., 2014). Observational constraints are thus required by models to derive more realistic results (Cox et al., 2013).

In this study, we explore a simple linear scheme to quantify the dynamic responses of global atmospheric CO₂ concentration to anthropogenic carbon emissions and global temperature changes based on observational records. Although the coupled climate-carbon system is nonlinear in nature, it can be linearized around a (relatively) steady point within a neighborhood in its state space (Khalil, 2001), which is evident in the fact that the atmospheric CO₂ concentration and the corresponding global climatology had been relatively stable for thousands of years before the industrial era (IPCC, 2013). The literature is rich with respect to the application of linear models to study the dynamics of the global carbon cycle or to diagnose its characteristics (e.g., Oeschger and Heimann, 1983; Maier-Reimer and Hasselmann, 1987; Enting and Mansbridge, 1987; Wigley, 1991; Jarvis, 2008; Boer and Arora, 2009, 2013; Gloor et al., 2010; Joos et al., 1996, 2013). In particular, these previous studies suggest that the characteristic impulse response function (IRF)——or more generally, the Green's function——of a complex carbon-cycle model to external disturbances of carbon emissions, can be captured by a few linear modes (Maier-Reimer and Hasselmann, 1987; Young, 1999; Joos et al., 1996, 2013). The states of atmospheric CO₂ calculated through convolutions of the simplified IRF and records of CO₂ emissions agree well with the corresponding simulations of the "parent" carbon cycle model (Wigley, 1991; Li et al., 2009). However, few studies (if any) have studied the dynamic responses of atmospheric CO₂ to disturbances of temperature changes using linear models, which is a main interest of this paper.

This study also intends to address a couple of important issues that, to the knowledge of the authors, have not been rigorously discussed in the literature on the development of simple diagnostic models for the global carbon cycle. One such issue is the determination of key state variables to include in the models so that they will not neglect important dynamic features of the real-world system. For instance, a common approach in the literature has been to represent the net carbon flux from the atmosphere to the surface as a linear term that depends on the anomalous atmospheric CO₂ concentration and the residence time of such anomalies (e.g., Gloor et al., 2010; Raupach et al., 2014). Sometimes, an additional term is also included to account for the residual fluxes from processes such as reforestation (Rayner et al., 2015). However, the surface carbon storage as a state variable (or variables) has tended to be neglected in these studies. The diagnostic framework developed in Boer and Arora (2009, 2013) explicitly considers the feedbacks of a warming climate on the global carbon cycle; but it also neglects the influence of the surface reservoirs on the atmosphere. Because variations of surface carbon reservoirs play a fundamental role in regulating carbon fluxes from the surface to the atmosphere, their absence in these previous models leads to significant restrictions on the simulated carbon cycle dynamics. Indeed, because of these restrictions, some previously reported results should be interpreted with caution (see below).

Another practical factor to decide upon in developing a diagnostic model is the complexity of the linear tool itself. Initial evaluation of the observational datasets used in this study indicates that they only allow the retrieval of a few independent system parameters (see below). Therefore, in the main text we only demonstrate our analytical framework by a simple two-box model that represents carbon exchanges between the atmosphere and the surface (i.e., land and ocean) reservoirs. Although such a "toy" model may sit at the lowest rank in the hierarchy of global carbon cycle models (Enting, 1987), its simplicity does not prevent it from shedding light on important and less well-known characteristics of the atmospheric CO₂ dynamics. The use of a simple model by no means implies a compromise in the scientific rigor of our findings, which we analytically verified with a generalized linear model comprised of an arbitrary large number of carbon reservoirs. For the sake of simplicity, these analytical proofs are not supplied with this paper, but will be published separately in the future.

Finally, throughout the analysis we also compare the results obtained from the two-box model to those from the more advanced Bern model (Siegenthaler and Joos, 1992; Enting et al., 1994; IPCC, 1996, 2001). The Bern model couples the atmosphere with a process-based ocean biogeochemical scheme (Siegenthaler and Joos, 1992; Shaffer and Sarmiento, 1995; Joos et al., 1999) and a multi-component terrestrial biosphere module (Siegenthaler and Oeschger, 1987). The original Bern model does not consider the effects of changing global temperatures on terrestrial ecosystem respiration, which plays an important role in regulating the variability of the global carbon cycle at interannual to multidecadal time scales (Woodwell et al., 1998; Rafelski et al., 2009; Willeit et al., 2014). Therefore, we revised the Bern model to account for the effects of temperature on terrestrial ecosystem respiration and subsequently recalibrated the model (see Appendix for details). The global carbon cycle processes described in the Bern model help us diagnose the biogeochemical mechanisms underlying the characteristics of the atmospheric CO₂ dynamics identified with our simple linear model.

The annual atmospheric CO₂ concentration data from 1850 to 1960 are based on the ice core CO₂ records from Law Dome, Antarctica (Etheridge et al., 1996), and those between 1960 and 2010 are compiled from the NOAA (National Oceanic and Atmospheric Administration) Earth System Research Laboratory (ESRL) (Conway et al., 1994; Keeling et al., 1995). We merged the data following the approach described in (Le Quéré et al., 2009) and calculated the annual CO₂ growth rate as the first-order difference of the yearly CO₂ concentrations. Long-term records of anthropogenic CO₂ emissions from fossil fuel burning and cement production compiled by (Boden et al., 2011), and those of land-use changes from (Houghton, 2003) were both downloaded from the Carbon Dioxide Information Analysis Center at Oak Ridge National Laboratory, TN, USA (http://cdiac.ornl.gov). Two sets of monthly surface temperature data are used, including GISTEMP (GISS Surface Temperature Analysis) from the NASA (National Aeronautics and Space Administration) Goddard Institute for Space Studies (Hansen et al., 1999) and the CRU-NCEP (Climatic Research Unit-National Centers for Environmental Prediction) climate dataset (Sitch et al., 2008; Le Quéré et al., 2009), available from 1901 to the present with spatial resolutions of 0.5^°× 0.5^° (CRU-NCEP) or 1^°× 1^° (GISTEMP). Monthly time series of temperature are aggregated globally and over the tropics (24^°N-24^°S), and smoothed with a 12-month running window to convert the monthly data to annual values. We calculated temperature anomalies relative to their 1901 to 1920 annual mean and assumed the 20-year mean temperature to be representative of temperature climatologies between 1850 and 1900. This assumption is reasonable, as suggested by analysis of other long-term coarse-resolution temperature datasets (Jones and Moberg, 2003; Brohan et al., 2006).

This study considers only the "fast" carbon flows between the atmosphere and the surface at time scales within hundreds of years (IPCC, 2001). Our two-box approach treats the atmosphere as one carbon reservoir ("box") and the world's land and oceans as the other one. Following the mathematical notation developed in Boer and Arora (2009, 2013), we describe the linearized two-box carbon system as follows: \begin{align} \dot{H}'_{ A}&=-B_{ A}H'_{ A}+B_{ S}H'_{ S}+\Gamma T'+\dot{E}'(1a) ,\\ \dot{H}'_{ S}&=+B_{ A}H'_{ A}-B_{ S}H'_{ S}-\Gamma T' , (1b)\end{align} where H _A and H _S denote carbon storages in the atmosphere and the surface reservoirs, with the subscripts "A" and "S" indicating the corresponding reservoirs, respectively. E is the accumulated anthropogenic CO₂ emissions since the industrial era. The three variables can be measured by the same unit of parts per million by volume (1\; ppm=∼2.13× 10¹⁵ grams of carbon, or 2.13 PgC). The prime symbol (e.g., "E'") indicates anomalies of a variable relative to its preindustrial steady-state level. The preindustrial emissions are assumed negligible so that E'≈ E (this assumption does not affect the results reported in this paper). The dot accent (e.g., "\(\dot E'\)") indicates the first-order derivative with regard to time, such that \(\dot E'\) represents the annual rate of CO₂ emissions (ppm yr^-1). The positive parameters B _A and B _S (yr^-1; note that "B" is the uppercase Greek letter "beta") describe the decaying rates of corresponding carbon anomalies. Their reciprocals (i.e., τ _A=1/B _A,τ _S=1/B _S) are often referred to as the response time of the carbon reservoirs (IPCC, 2001). T (^°C) denotes indices of global (or large-scale) surface temperatures and the parameter Γ (ppm yr^{-1 °}C^-1) indicates the rate at which carbon fluxes are released from the surface reservoirs to the atmosphere per unit change in temperature. The Γ parameter does not necessarily represent the responses of atmospheric CO₂ concentration to temperature changes, for the estimation of the latter must also take the dynamics of the system into account. This study assumes all three parameters (B _A,B _S, and Γ) to be constant, which is found to be a good approximation for the global carbon cycle over the time period under investigation (see below).

The inclusion of the term B _SH' _S in Eq. (1) reflects an important structural difference between this paper and the aforementioned previous studies (e.g., Boer and Arora, 2009, 2013; Gloor et al., 2010; Raupach et al., 2014; Rayner et al., 2015). This extra term helps complete the modeled carbon cycle dynamics to the first-order approximation. Carbon outflows from the atmosphere, B _AH' _A, are inflows to the surface (e.g., through photosynthesis in green vegetation and the dissolution of CO₂ in surface water), while carbon outflows from the surface, B _SH' _S, (e.g., through respiration and the outgassing of dissolved CO₂) are the inflows to the atmosphere. The effects of temperature changes, Γ T', revise the relative carbon balance between the atmosphere and the surface reservoirs. Human emissions of CO₂, on the other hand, represent an "external" source of CO₂ to the system by rapidly releasing carbon (e.g., fossil fuel burning) from reservoirs that were formed over millions of years and by permanently altering the structure of land surface carbon pools (e.g., land-cover/land-use changes). It is easy to check that neglecting the B _SH' _S term is equivalent to setting B _S to zero or setting τ _S to infinity, which implies that the surface reservoirs will never saturate. As will be discussed later, such negligence can introduce large biases to the simulated carbon cycle over long-term time scales.

Because mass (carbon) is conserved in the two-box model, Eqs. (1a) and (1b) are not independent. Adding the two equations together, we can easily see that $$ H'_{ A}+H'_{ S}=E' ,(1c) $$ which simply states that the anthropogenically emitted CO₂ either resides in the atmosphere or in the surface reservoirs (i.e., the land and the oceans). Substituting this relationship into Eq. (1a) to replace H' _S, we obtain $$ \dot{H}'_{ A}+(B_{ A}+B_{ S})H'_{ A}=\Gamma T'+B_{ S}E'+\dot{E}' . (2a) $$ Therefore, the dynamics of atmospheric CO₂ represented by the two-box model is determined by an ordinary differential equation of H' _A under the disturbances of anthropogenic emissions (E' and \(\dot E'\)) and the changing climate (T').

To determine the parameters of Eq. (2a) with observational records of H' _A, T' and E', we rearrange the equation, as follows, to construct a regression model: $$ \dot{E}'-\dot{H}'_{ A}=(B_{ A}+B_{ S})H'_{ A}-B_{ S}E'-\Gamma T' , (2b) $$ where \(\dot E'-\dot H'_ A\) represents the strength of annual carbon sinks. However, the direct regression analysis of Eq. (2b) is obscured by the "collinearity" among the regressors (Chatterjee and Hadi, 2006). Indeed, because H' _A≈α _LTE' (Fig. 1), substituting this long-term (indicated by the subscript "LT") airborne-fraction relationship into Eq. (2b) leads to $$ \dot{E}'-\dot{H}'_{ A}=[B_{ A}-\left(\dfrac{1}{\alpha_{ LT}}-1\right)B_{ S}]H'_{ A}-\Gamma T' , (2c) $$ which implies that only a combination of B _A and B _S can be estimated from the regression analysis. Our estimate of the overall coefficient B _A-(1/α _LT-1)B _S is 0.04 0.001 (yr^-1), while α _LT is separately estimated to be 0.41 (Fig. 1).

Additional information is required to resolve B _A and B _S. One source of such additional information comes from previous studies based on the measurements of carbon isotope ratios in wood and marine material (Revelle and Suess, 1957), from which the response time (τ _A) of atmospheric CO₂ is inferred to be on the order of 10 years. We can also extract information from process-based model studies. Because τ _A determines the initial decaying rate of the IRF of a global carbon cycle model (see the proof outlined in the next section), applying this result to analyze the ensemble IRFs reported in (Joos et al., 2013) suggests τ _A to be ∼14 years. In this study, we choose τ _A to be 12 years (B _A≈0.083 yr^-1) so that the IRF of our linear model closely matches with the Bern model during the initial decaying stage (see the next section). Subsequently, we estimate τ _S to be ∼34 years (B _S≈0.029 yr^-1).

The estimation of the Γ parameter in Eq. (2c) requires the choice of a large-scale temperature index that is representative of climate change and closely related with the global carbon cycle. Our previous study shows that the land surface air temperature in the tropics (24^°S-24^°N) is most strongly coupled with the interannual variations of the growth rate of atmospheric CO₂ by a sensitivity (Γ) of ∼1.64 0.28 ppm yr^{-1 °}C^-1 (Wang et al., 2013). Here, we find the same temperature sensitivity parameter to be valid for Eq. (2c). Indeed, because the system is linear, variations of all the variables in Eq. (2c) over different time scales must satisfy the equation separately. Because the interannual variations ("IAV") of the carbon emissions (\(\dot E'\) and E') and the atmospheric CO₂ concentration (H' _A) are relatively small (Fig. 1), neglecting them in Eq. (2c) leads to $$ \dot{H}'_{ A,IAV}\approx\Gamma T'_{ IAV} ,(2d) $$ which is the same linear relationship used in (Wang et al., 2013) and other previous studies (Adams and Piovesan, 2005). In addition, because the long-term trends in global temperature (T') and atmospheric CO₂ concentrations (r≈ 0.9) are significantly correlated (H'_A), estimating Γ directly from Eq. (2c) is also subject to the influence of the collinearity between the two fields. This practical concern also makes it more reasonable to estimate Γ based on Eq. (2d) than Eq. (2c).

With the above parameters determined, we use the two-box model to simulate changes in atmospheric CO₂ concentration between 1850 and 2010 from historical records of temperature and carbon emissions (Fig. 2). The results reproduce the evolution of the observed CO₂ time series to a high degree of accuracy, capturing more than 96% of the variability (i.e., r²>0.96) of the latter (Fig. 2). The standard deviations (σ) of the differences between the simulations and the measurements since 1960 are ∼0.9 ppm for the atmospheric CO₂ concentration and ∼0.4 ppm yr^-1 for its growth rate, respectively (Fig. 2). These results are highly comparable to those simulated with the revised Bern model (Fig. 2) or other sophisticated climate-carbon models reported in the literature (e.g., Joos et al., 1999; Lenton, 2000; Friedlingstein et al., 2006), demonstrating that the atmospheric CO₂ dynamics in the past one and half centuries can be properly approximated with suitable lower-order linear models.

Figure 2.  Simulations of the observed atmospheric CO$_2$ concentrations (top panel) and growth rates (bottom panel) from anthropogenic CO$_2$ emissions and land-surface air temperature data using the two-box model ("2box") and the revised Bern model ("Bern"). Close-up views of the simulations between 2000 and 2010 (lightly shaded) are shown in the inset figures. The atmospheric CO$_2$ concentration in 1850 (i.e., 284.7 ppm) is used as the initial condition for the model integration. Long-term mean temperature before 1901 is assumed to be stable and represented by the 1901-20 mean. Other model parameters used in these simulations are explained in the main text (the two-box model) or Appendix (revised Bern model).

Figure 3.  Disturbance-response functions of the atmospheric CO$_2$ concentration simulated by the two-box model ("2box") and the revised Bern model ("Bern"). The top panel shows the responses of atmospheric CO$_2$ concentration to an impulse increase (of 100 ppm) in anthropogenic CO$_2$ emissions and the bottom panel shows the corresponding responses to a step increase (of 1$^◦$C) in surface temperatures.

We first check the model’s responses to an impulse disturbance of anthropogenic CO₂ emissions. Shown in Fig. 3, the initial atmospheric CO₂ anomaly decays relatively fast, as 60%-70% of the emitted CO₂ is absorbed by the surface reservoirs within 20 years of the disturbance. However, the rate of carbon assimilation by the land and the oceans slows down significantly in the following decades and eventually becomes neutral as the system approaches a steady state, suggesting that 15%-25% of the simulated CO₂ anomaly will likely stay in the atmosphere for thousands of years (Fig. 3). These results are consistent with the findings from fully coupled climate-carbon models (Archer et al., 2009; Cao et al., 2009; IPCC, 2013; Joos et al., 2013).

The IRF of the linear models can be analytically characterized. For the model of Eq. (2a), when the system approaches a (new) steady state after the disturbance, all the time derivatives (\(\dot E'\) and \(\dot H'_ A\)) will be zero. Assuming that temperature does not change during the process, we easily obtain the steady state of H' _A as $$ H'_{ A}=\dfrac{B_{ S}}{B_{ A}+B_{ S}}E'=\dfrac{\tau_{ A}}{\tau_{ A}+\tau_{ S}}E' , (3a) $$ or more generally, $$ \dfrac{H'_{ A}}{\tau_{ A}}=\dfrac{H'_{ S}}{\tau_{ S}} , (3b) $$ where the mass-conservation relationship represented by Eq. (1d) is used in the derivation. Therefore, the extra CO₂ added to the "fast" carbon cycle by anthropogenic emissions will be partitioned between the atmosphere and the surface corresponding to the response times (τ) of the reservoirs, respectively (Revelle and Suess, 1957). Because τ _S>τ _A, a majority of the emitted CO₂ will eventually be absorbed by the surface carbon reservoirs (Fig. 3). If τ _S approaches infinity, the proportion of the emitted CO₂ that stays in the atmosphere would decrease to zero within a few centuries, which is, however, unrealistic (IPCC, 2013). This discrepancy highlights a key limitation in neglecting the influence of surface carbon storage in previous diagnostic models discussed in section 3.

The rates at which the atmospheric CO₂ anomaly decays are determined by the solutions (i.e., eigenvalues) to the characteristic equation of the system. For a two-box system like Eq. (2a), the problem is particularly simple because the only non-zero eigenvalue (Λ) is $$ \lambda=B_{ A}+B_{ S} ,(4) $$ and the solution of Eq. (2a) is therefore $$ H'_{ A}=\dfrac{B_{ A}}{B_{ A}+B_{ S}}\exp[-(B_{ A}+B_{ S})t]+\dfrac{B_{ S}}{B_{ A}+B_{ S}} . (5a) $$ A helpful observation of Eq. (5a) is that, when t« 1/(B _A+B _S), the solution can be approximated by $$ H'_{ A} \approx \dfrac{B_{ A}}{B_{ A}+B_{ S}} [1-(B_{ A}+B_{ S})t]+\dfrac{B_{ S}}{B_{ A}+B_{ S}}=1-B_{ A}t \approx \exp(-B_{ A}t) . (5b) $$ That is, H' _A initially decays at a maximum rate of B _A.

Next, we consider the system's responses to disturbances induced by changes in surface temperatures. Unlike anthropogenic CO₂ emissions, changes in temperature do not add additional CO₂ to the "fast" carbon cycle, but only redistribute carbon between the atmosphere and the surface [Eqs. (1a) and (1b)], and so the system will recover to its initial steady state once the temperature anomaly is removed. However, increases in temperature are persistent under climate change scenarios. Therefore, we examine the long-term responses of atmospheric CO₂ to a step change in temperature, which is easily determined from Eq. (2a) as: $$ H'_{ A}=\dfrac{\Gamma}{B_{ A}+B_{ S}}T' .(6a) $$ Because B _A>B _S, for quick estimates we can also use $$ H'_{ A}\approx\dfrac{\Gamma}{B_{ A}}T'=\tau_{ A}\Gamma T' . (6b) $$

Equations (6a) or (6b) indicates that, on the first order, the long-term temperature sensitivity of atmospheric CO₂ concentration is largely determined by the product of the Γ parameter and the response time (τ _A) of the atmosphere carbon reservoir. Based on the model parameters retrieved in the preceding section, we estimate that atmospheric CO₂ will rise by ∼15 ppm for an increase of 1^°C in temperature over multidecadal to centennial time scales (Fig. 3). This result generally agrees with the estimates inferred from (reconstructed) temperature and atmospheric CO₂ records over the Little Ice Age [∼20 ppm ^°C^-1 (Woodwell et al., 1998)], the past millennium [1.7-21.4 ppm ^°C^-1 (Frank et al., 2010)], or estimates from coupled model experiments [5-12 ppm ^°C^-1 (Willeit et al., 2014)].

The relationships represented by Eqs. (3b), (5b) and (6b) can be generalized to arbitrary high-order linear systems. The uncertainties associated with these results——especially the long-term responses of atmospheric CO₂——need to be emphasized. One key source of the uncertainties is that the model's parameters are not fully determined by the observations of the global climate-carbon system. As discussed in section 4, the estimation of the model parameter B _S depends on the choice of B _A, which is only loosely constrained by the prior knowledge. General analysis indicates that this situation only worsens in higher-order (N-box) systems as the number of system parameters increase by the order of N² (also see Joos et al., 1996). It is possible for us to choose another pairing of B _A and B _S, or a higher-order linear model, so that the derived disturbance response functions better approximate those of the Bern model (Fig. 3). However, tuning the model in this fashion has only cosmetic effects on the results and does not reduce the associated uncertainties. Furthermore, in the real world, the climate system and global carbon cycle are not independent, but tightly coupled. A comprehensive assessment of the long-term fate of anthropogenic CO₂ emissions in the atmosphere must account for the effects of the associated changes in global temperature, which is beyond the scope of this study.

The above analysis strongly suggests that the appropriate representation of the effects of temperature on the carbon cycle in our linear model helps improve the model's accuracy in approximating the observed dynamics of the atmospheric CO₂ across multiple time scales. To illustrate, we further rearrange Eq. (2c) to obtain $$ \dot{E}'-\dot{H}'_{ A}+\Gamma T'=\left[B_{ A}-\left(\dfrac{1}{\alpha_{ LT}}-1\right)B_{ S}\right]H'_{ A} .(7) $$ On the left-hand side of the equation, the term "\(\dot E'-\dot H'_ A\)" is usually used to measure the net strength of annual global carbon sinks. However, because the warming temperature also releases carbon from the surface into the atmosphere (Γ T'), this extra source of CO₂ has to be absorbed by the global carbon sinks. By accounting for the effects of temperature changes, the term "\(\dot E'-\dot H'_ A+\Gamma T'\)" thus defines the gross global carbon sinks.

Figure 4.  Global annual carbon sinks (ppm yr$^-1$) as a function of atmospheric CO$_2$ concentration from 1850 to 2010. The green dots indicate the observed "net" carbon sinks and the red dots indicate the "gross" carbon sinks that accounted for the effects of temperature changes (Eq. 7). The differences between the gross and the net carbon sinks (shaded area) indicate the extra carbon fluxes released into the atmosphere as a result of warming temperatures (Fig. 1). The gray arrow ("$H_ A,0$") indicates the estimated atmospheric CO$_2$ level (284.7 ppm) that was stable at preindustrial CO$_2$ emission rates and climate conditions. The slopes between the global annual carbon sinks and corresponding changes in atmospheric CO$_2$ concentration (relative to $H_ A,0$) were interpreted as the carbon sequestration efficiency of global land and ocean reservoirs in the literature (Gloor et al., 2010; Raupach et al., 2014). However, such an interpretation is valid only when the long-term airborne fraction ("$\alpha_ LT$") of anthropogenic CO$_2$ emissions is approximately constant (see discussion in the main text).

Examining Eq. (7) with the observational data shows that both the net and gross carbon sinks have been steadily increasing in response to the rising atmospheric CO₂ concentration in the past 160 years, reaching ∼2.5 ppm yr^-1 and ∼4.0 ppm yr^-1, respectively, in 2010 (Fig. 4). The gross carbon sinks have a near direct linear relationship (with a constant slope of ∼0.04 yr^-1; r=0.98) with the atmospheric CO₂ concentrations throughout the entire data period. In comparison, the relationship between the apparent carbon sinks and the CO₂ concentrations is slightly nonlinear, with its slope decreasing from ∼0.03 yr^-1 in 1960 to ∼0.02 yr^-1 in 2010. Therefore, our linear approximation approach would not be able to achieve the same high accuracy if the effects of temperature on the carbon cycle were not correctly represented. Note that the slopes of these linear relationships are sometimes interpreted as the efficiency of surface carbon reservoirs in sequestering annual CO₂ emissions (Gloor et al., 2010; Raupach et al., 2014). Figure 4 shows that although the gross carbon sequestration rates B _A-(1/α _LT-1)B _S of the surface reservoirs changed little, the net "efficiency" (the ratio between \(\dot E'-\dot H'_ A\) and H' _A) of the system has slowed by 28%-33% in the past five decades. This finding is essentially the same as reported in (Raupach et al., 2014), but our analysis emphasizes that this declining carbon sequestration rate mainly reflects the impacts of climate changes on the global carbon cycle. Also, it must be cautioned that since the coefficient represented by B_A-[(1/α _LT-1)B _S] is influenced by the long-term AF factor (α _LT), it is not an intrinsic characteristic of the carbon cycle system. The interpretation of the coefficient as "carbon sink efficiency" is only meaningful when α _LT is constant, which is largely valid in the historical records but subject to changes in the future (IPCC, 2013).

The constancy of the model parameters B _A, B _S and Γ needs further discussion. Previously, Boer and Arora (2009, 2013) developed a linear diagnostic framework to quantify the coefficients of CO₂ concentration-and temperature-related carbon cycle feedbacks (i.e., B _A and Γ) with the Canadian Earth System Model. Their results show long-term trends in both parameters as CO₂ concentration and temperature increase through the 21st century (2009, 2013). To explain the discrepancies between this study and Boer and Arora (2009, 2013), we notice that the term B _SH' _S was neglected in the previous studies. Therefore, the "B _A" parameter estimated in their studies is actually the bulk coefficient, B _A-[(1/α _LT-1)B _S], associated with the CO₂ concentration anomaly (H' _A) in Eq. 2(c). As discussed earlier in this paper, this coefficient is sensitive to future changes in the long-term AF factor (α _LT). Also, changes in atmospheric CO₂ concentration and global temperature are much higher under the future emission scenarios studied in Boer and Arora (2009, 2013) than the historical period investigated in this study. The changes of the parameters (B _A and Γ) reported in Boer and Arora (2009, 2013) may also partially reflect the fact that the global carbon cycle is gradually pushed away from its linear resilience zone in their model experiments.

To explain the biogeochemical meaning of the Γ parameter, our previous analysis (Wang et al., 2013) suggests that it mainly reflects the temperature sensitivity of respiration of land-surface carbon pools (biomass and soil carbon). This explanation is supported by the simulations of the Bern model in this study, in which terrestrial carbon sinks have much stronger responses to temperature changes than the oceanic counterpart (not shown). Furthermore, both our simulations and those from the literature (e.g., Canadell et al., 2007; Le Quéré et al., 2009) indicate that the total carbon storage in the land-surface reservoirs remains largely stable between 1850 and 2010, which is a necessary condition for Γ to be constant. For instance, because terrestrial carbon uptake accounts for 50%-60% of the global net sinks in our simulations, the accumulated terrestrial net carbon sinks are about 71-85 ppm in 2010, representing a 7%-8% increase in the total terrestrial carbon storage (∼1040 ppm as of 1850). At the same time, the accumulated terrestrial carbon losses through land-use changes are about 74 ppm in 2010 based on the dataset of (Houghton, 2003). These results suggest that the net changes in the total terrestrial biomass and soil carbon are (relatively) small during the past 160 years, providing further justification for our linear modeling approach.

Finally, because of the buffering effect of the ocean carbonate chemistry (Revelle and Suess, 1957), the responses of oceanic carbon uptake to changes in the atmospheric CO₂ concentration are generally estimated to be small compared with its terrestrial counterpart (e.g., Le Quéré et al., 2009). Our analysis thus suggests that the increasing atmospheric CO₂ concentration must have promoted carbon assimilation by the terrestrial biosphere (Ballantyne et al., 2012), most likely through the CO₂ fertilization effect (Long, 1991; Körner and Arnone, 1992; Oechel et al., 1994; Long et al., 2004) and the associated ecological changes (Keenan et al., 2013; Graven et al., 2013). Indeed, because the surface warming rapidly releases a proportion of the assimilated carbon back to the atmosphere (Fig. 4) (Piao et al., 2008; Wang et al., 2013), the increased turnover rate may have obscured the evaluation of the magnitude of the CO₂ fertilization effects, which we found in calibrating the Bern model (see Appendix). In other words, the gross CO₂ fertilization effect of terrestrial vegetation is likely higher than previously thought (Schimel et al., 2014).

This paper develops a simple linear model to describe carbon exchanges between the atmosphere and the surface carbon reservoirs under the disturbances of anthropogenic CO₂ emissions and global temperature changes. We show that, with a few appropriately retrieved parameters, the model can successfully simulate the observed changes and variations of the atmospheric CO₂ concentration and its first-order derivative (i.e., CO₂ growth rate) across interannual to multi-decadal time scales. The results are highly comparable to those obtained with more sophisticated models in the literature, confirming that the simple linear model is capable of capturing the main features of atmospheric CO₂ dynamics in the past one and half centuries.

A distinct advantage of our linear modeling framework is that it allows us to analytically, and thus most directly, examine the dynamic characteristics of the (modeled) carbon cycle system. Our analyses indicate that many such characteristics are closely associated with the response times of the atmosphere and surface carbon reservoirs. For instance, the response time of the atmosphere determines the initial decaying rate of an impulse of CO₂ emitted into the atmosphere, and plays a major role in connecting the short-term and long-term temperature sensitivities of atmospheric CO₂ concentration. The ratio between the response times of the atmosphere and the surface reservoirs also determines the proportion of the CO₂ emissions that stays in the atmosphere at long-term time scales. Unfortunately, the collinearity exhibited by the observed time series of CO₂ emissions and atmospheric CO₂ concentrations has obscured the determination of the response times for individual surface reservoirs, inducing uncertainties to the estimated long-term responses of the global carbon system.

Our model results have important biogeochemical implications. They highlight that the responses of the global carbon cycle to recent anthropogenic and climatic disturbances are still within the resilience zone of the system, such that annual (gross) terrestrial and ocean carbon sinks linearly increases with the atmospheric CO₂ levels. On the one hand, the elevated atmospheric CO₂ concentration must have enhanced land carbon uptakes through the "fertilization" effects and the associated ecological changes. On the other hand, the enhanced gross carbon uptakes are partially offset by the increases in global surface temperatures, which accelerate the release of carbon from the surface reservoirs into the atmosphere. As a result, the "net" efficiency of global land and oceans in sequestering atmospheric CO₂ may have slowed by ∼30% since the 1960s, although the airborne fraction of CO₂ emissions remains largely constant.

Finally, and importantly, we emphasize that the linear approximation of the global carbon cycle discussed in this paper is conditioned on the preindustrial (quasi) steady state of the system. The global climate-carbon system is clearly nonlinear beyond this scope (Archer et al., 2009), which can establish different steady states over glacial/interglacial time scales (Sigman and Boyle, 2000). A major concern stemming from climate change is that, because the post-industrial anthropogenic disturbances on the global carbon cycle are so strong and rapid, they may abruptly alter the pace at which the natural climate-carbon system evolves, and drive the system into a different state at a drastically accelerated rate (IPCC, 2001). Our results clearly indicate that the rising atmospheric CO₂ concentrations and the associated increases in global temperature have significantly intensified the global carbon cycle in the past one and half centuries. Although such intensification of the carbon system seems to be within the linear zone as of now, its resilience may be weakened, or lost, in the future. As the anthropogenic CO₂ emissions continue to increase and the global temperature continues to warm, scientists generally expect surface-in particular, terrestrial-carbon reservoirs to saturate and their CO₂ sequestration efficiency to decrease, such that the responses of the global carbon cycle to the anthropogenic disturbances will eventually deviate from their original path. With this concern in mind, the simple linear model developed in this study may serve as a useful tool to monitor the early signs of when the natural carbon system is pushed away (by anthropogenic disturbances) from its linear zone.

Calibrations of the Bern Carbon Cycle Model

The Bern model is a coupled global carbon cycle box model (Siegenthaler and Joos, 1992;Enting et al., 1994) that was used in previous IPCC (Intergovernmental Panel on Climate Change) Assessment Reports to study changes in atmospheric CO₂ concentration under different emissions scenarios (IPCC 1996,2001).It couples the High-Latitude Exchange/Interior Diffusion--Advection (HILDA) ocean biogeochemical model (Siegenthaler and Joos, 1992; Shaffer and Sarmiento, 1995; Joos et al., 1999) with an atmosphere layer and a multi-component terrestrial biosphere model (Siegenthaler and Oeschger, 1987).The HILDA model describes ocean biogeochemical cycling through two well-mixed surface layers in low and high latitudes,a well-mixed deep ocean in the high latitudes, and a dissipative interior ocean in the low latitudes. Ocean tracer transport is represented by four processes:~(1) eddy diffusion within the interior ocean (k, 3.2\times 10^-5 m² s^-1); (2) deep upwelling in the interior ocean (w, 2.0\times 10^-8 m s^-1), which is balanced by lateral transport between the two surface layers as well as the downwelling in the polar deep ocean; (3) lateral exchange between the interior ocean and the well-mixed polar deep ocean (q, 7.5\times 10^-11 s^-1); and (4) vertical exchange between the high-latitude surface layer and the deep polar ocean (u, 1.9\times 10^-6 m s^-1) (Shaffer and Sarmiento, 1995). The effective exchange velocity between surface ocean layers and the atmosphere in both low and high latitudes is assumed to be the same (2.32\times 10^-5} m s^-1) (Shaffer and Sarmiento, 1995). Ocean carbonate chemistry (e.g., the Revelle buffer factor) is based on the formulation given by Sarmiento et al. (1992). In addition, we implemented the influence of SST on the partial pressure of dissolved CO₂ in seawater with a sensitivity of \sim 4.3{\%} ^◦C^-1 (Gordon and Jones, 1973; Takahashi et al.,1993; Joos et al., 2001). The changes in global mean SST is approximately 0.8^◦C--1.0^◦C from the 1850s to 2000s (Rayner et al., 2003; Brohan et al., 2006), slightly lower than that of the tropical land-based air temperature (\sim 1.0^◦C) but with a trend resembling the latter (Rayner et al., 2003; Jones and Moberg, 2003; Hansen et al., 2006). For simplicity, therefore, we used the long-term trend of the tropical land air as a proxy for the corresponding trend in global SST.

The terrestrial biosphere in the Bern model is represented by four carbon compartments (ground vegetation, wood, detritus, and soil) with prescribed turnover rates and allocation ratios. The global net primary production (NPP), the influx to the biosphere, is assumed to be 60 PgC yr^-1 at the preindustrial level; and the effect of CO₂ fertilization on NPP (i.e., the \beta-effect) is described with a logarithmic function with a \beta parameter of 0.38 (Enting et al., 1994). The original Bern model does not consider the effects of changing global temperatures on terrestrial ecosystem respiration, which have been suggested to play an important role in regulating the variability of the global carbon cycle at interannual to multidecadal time scales (Rafelski et al., 2009; Wang et al., 2013).Therefore, we implemented the effects of temperature on terrestrial ecosystem respiration in the Bern model with an overall sensitivity (Q₁₀) of \sim 1.5 (Lenton, 2000; Davidson and Janssens,2006; Wang et al., 2013). We also changed the preindustrial CO₂ concentration to 285 ppm in the Bern model to reflect the findings obtained from the observations (Fig.~4 of the main text).

We calibrated the Bern model so that the model outputs fit the observed atmospheric CO₂ data most favorably. Because no major revisions were made to the ocean carbon cycle module (HILDA),we focused mainly on calibrating the biosphere module. With the original biosphere model parameters, the simulated atmospheric CO₂ concentrations were found to be distinctly higher than observations, reaching \sim 411 ppm in 2010. These results are induced because rising temperatures enhance respiration in the model, reducing the net land carbon sinks to an unrealistic \sim 0.5 ppm yr^-1 in 2010. To balance the temperature-enhanced respiration, we need to increase the \beta parameter from 0.38 to 0.64 to incorporate a higher rate of gross biosphere carbon uptake, as enhanced by CO₂ fertilization (Long et al., 2004) and the associated ecological changes Keenan 2013}. With the \beta parameter set at 0.64, the simulated global terrestrial NPP increased by 14% from its preindustrial level and reached \sim 69 PgC yr^-1 in 2010, which qualitatively agrees with recent estimates inferred from isotope measurements (Welp et al., 2011). As such, the recalibrated Bern model is able to simulate accurately the observed changes/variations in atmospheric CO₂ concentration and growth rate in the past 160 years (Fig.~2 of the main text). The simulated ocean and land components of global carbon sinks are also consistent with estimates found in previous studies (e.g., Canadell et al., 2007; Le Quéré et al., 2009).