Response function analysis of carbon dioxide and climate using the Padé-Laplace technique

1.   Introduction

Impulse response functions are widely used as a characterisation of causal linear systems. Such impulse response functions have played various roles in the study of the carbon cycle. Notably the CO$_2$ impulse response function is central to the definition of Global Warming Potential (GWP), the factor that has been used in international agreements to convert emissions of other gases to a CO$_2$-equivalent ^[1]. Numerical fits to impulse response of models, commonly expressed as sums of exponentials, have been a common way of summarising model behaviour. Impulse response function descriptions for subsystems (e.g. ocean carbon mixing ^[2,3]) have often been used as components of more comprehensive models. Impulse response functions have also been used to characterise how global temperature responds to changes in radiative forcing, and, as illustrated below, can be useful in analysing the climate-carbon feedback loop.

In carbon cycle studies, Laplace transforms of impulse response functions have been used conceptually ^[4,5,6,7] to characterise generic aspects of the carbon cycle. However apart from some minor exceptions, noted in Section 3.4.1, the analysis using Laplace transforms has not been the basis of calculations in the transform domain. The present paper explores the possibility and utility of such transform domain computations. The technique that is used is the Padé-Laplace approximation developed by Claverie and colleagues ^[8,9], and used in areas such as circuit analysis ^[10]. The technique involves fitting low degree Padé approximants (ratios of polynomials in the transform variable) to Taylor series expansions of relations involving the transforms of impulse response functions. Such approximations can define "reduced model" expressions that give good approximations to more complicated relations.

The outline of the remainder of this paper is as follows: Section 2 references the mathematical basis that underlies this paper, including citing the computer code used for the numerical examples. Section 3 introduces the concept of impulse response functions using atmospheric CO$_2$ levels as an example. It reviews the properties of the Laplace transform and introduces the Padé-Laplace technique which makes it straightforward to perform quantitative calculations in the transform domain. Section 4 considers the carbon system of atmosphere, land and oceans and analyses how different response functions are related and how they relate to descriptions in terms of feedbacks from concentration changes. Section 5 considers response functions that characterise how global temperature responds to changes in radiative forcing, and the extent to which climatic influence on the carbon cycle creates a feedback loop that amplifies both the temperature response and the CO$_2$ response. Section 6 compares the parameterisations of the key responses to parameterisations used in various other studies. Section 7 reviews the various results, concentrating particularly on the utility of the Padé-Laplace formalism. Possible generalisations and directions for future studies are noted.

2.   Methods

The Padé-Laplace analysis, described in Section 3.4, is based on the widely-used Laplace transform, described in Section 3.2, along with key generic properties. Specific properties of the Laplace transform that are required for the Padé-Laplace analysis are given in Section 3.3.

The numerical calculations in the examples and figures were performed using R code ^[11] that was based on the routines ^[12] used by Enting and Clisby ^[13]. The coefficients that specify the various response functions used in the examples are given in Table 1.

3.   Impulse response functions and Laplace transforms

3.1.   Carbon dioxide

Impulse response functions were introduced into carbon cycle studies by Oeschger and Heimann ^[14] with the relation:

where $M_ {\rm{CO2}}(t)$ is amount of CO$_2$ (carbon dioxide in gigatonnes of carbon) in atmosphere at time $t$ so that $Q(t)$ is the perturbation from the initial amount. $R_ {\rm{X}}(t)$ is an impulse response function that specifies the proportion of CO$_2$ remaining in atmosphere at time $t$ after emission. $S(t)$ is CO$_2$ emissions, particularly anthropogenic emissions, from "outside" the system which may comprise atmosphere and/or ocean and/or land vegetation. The integration starts at $t_0$, a time of notional pre-industrial "equilibrium", with $S(t)$ taken as zero for earlier times.

The type of response function depends on what is defined as internal to the system (generally excluding direct human influence). Some examples of the generic $R_ {\rm{X}}(t)$ are:

● ${R_ {\rm{O}}}$ System is "atmosphere+ocean". This form was introduced by Oeschger and Heimann ^[14].

● ${R_ {\rm{L}}}$ System is "atmosphere + land biota" (i.e. including living and dead biomass and soil organic carbon). This response characterises the so-called "CO$_2$-fertilisation effect".

● ${R}$ System is "atmosphere + ocean + land biota".

● ${R_ {\rm{FB}}}$ System is "atmosphere + ocean + land biota", including climate-to-carbon feedback.

Relations between these various response functions are given by equations (4.4) and (5.8, 5.9) below. Enting ^[15] has noted that in many cases, when estimates of $R(t)$ are derived from models calibrated against 20th century data, the estimates were actually estimates of ${R_ {\rm{FB}}(t)}$ since the feedback from temperature change was occurring through the 20th century.

Impulse response functions can also be defined for subsystems within the carbon cycle, e.g. for the ocean mixed layer ^[3]. Such subsystem responses, converted to differential equations as described below (equation (3.15)), are sometimes used within more comprehensive models ^[2]. The various response functions are often expressed as sums of exponentials, most commonly obtained by numerical fits to model output.

The response function for CO$_2$ has been of considerable importance because it is used in the definition of GWP (Global Warming Potential). This has been used to define the relative importance of different greenhouse gases and convert non-CO$_2$ emission rates to a CO$_2$-equivalent. Originally GWP was notionally defined using $R(t)$. More recent usage is based on ${R_ {\rm{FB}}(t)}$, which includes feedbacks explicitly. The change in GWP values is relatively small, because, as noted above, estimates of $R(t)$ based on models calibrated against 20th century CO$_2$ may be effectively estimates of ${R_ {\rm{FB}}(t)}$ since the climate-to-carbon feedback was operating in the 20th century. Joos et al. ^[1] reported an extensive intercomparison of CO$_2$ response functions and a number of cases from that study are used in the examples below. Table 1 lists the coefficients in various response functions used in this paper, expressed as sums of exponentials.

3.2.   Laplace transforms

The Laplace transform of a function $f(t)$ is defined as

This paper uses the tilde notation to denote the Laplace transforms of functions of time. Earlier work, for example ^[4,5,7,16], followed Abramowitz and Stegun ^[17] in using upper case for functions of time and lower case for the transforms. Since many recent descriptions, especially on-line presentations, reverse this case distinction, the present paper follows ^[13] in order to avoid confusion and uses the tilde notation: ($\tilde{f}(p)$ as the Laplace transform of $f(t)$).

Some of the properties of the Laplace transform that are important in the present analysis are linearity

and the derivative relation

where $\tilde{\dot{f}}(p)$ denotes the Laplace transform of the time derivative of $f(t)$.

A very important property is

which means that impulse response expressions, such as (3.1), that are expressed as convolutions become simple products in the transform domain.

Derivatives of convolutions can be expressed in various ways. Taking $\tilde{k}(p) = \tilde{f}(p)\, \tilde{g}(p)$ one has

In the time domain, some of these equivalent forms can be regarded as being related by "integration by parts".

Two important connections between values in the time domain and values in the transform domain are:

and, subject to all poles in $p\, \tilde{R}(p)$ being in the left half plane,

Two important specific cases are:

and

where $\delta(t)$ is the "Dirac delta functtion".

For exponential growth from the indefinite past, with $S(t) = C \exp(c t)$

by putting $t" = t-t'$. The asymptotic airborne fraction, $\kappa(c)$, defined as either instantaneous, $\kappa = \dot{Q}(t)/S(t)$ or cumulative, $Q(t)/\int_{-\infty}^t S(t')\, dt'$ is given by

3.3.   Relations using Laplace transforms

Laplace transforms have been used for describing relations between aspects of the carbon cycle ^[4,5,16]. Their relation to the time domain can be regarded as comparable to the use of vector notation to express, in a compact way, relations between multiple components. This can be taken a step further by using Laplace transform expressions to derive techniques for time domain analysis ^[18,19], as is summarised in Section 5.2 below.

Specific calculations in the transform domain can be readily done for functions expressed as sums of exponentials:

which has the transform

This sums to a ratio of polynomials in $p$, with numerator and denominator of degrees $M-1$ and $M$ respectively. Thus $f(0)$ will be given by $\sum_j a_j$, consistent with (3.7).

Expressing a sum of exponentials as a sum of functions: $\tilde{f}(p) = \sum \tilde{f}_j(p) = \sum a_j/(p+\lambda_j)$ means that relations of the form $\tilde{f}(p)\, \tilde{g}(p)$ can be expressed as

This corresponds to a set of first order differential equations in the time domain, as is well known, without reference to Laplace transforms, e.g. ^[20].

A special case occurs if one of the $\lambda_j$ (which will taken to be $\lambda_1$) is equal to zero, as applies for most representations of the CO$_2$ response, $R_ {\rm{X}}(t)$. In this case analysis of $p\, \tilde{R}_ {\rm{X}}(p)$ is better suited to Padé-Laplace analysis. Property (3.8) means that $a_1$ will give the long timescale limit of $p\, \tilde{R}_ {\rm{X}}(p)$ which, as noted above, corresponds to the asymptotic airborne fraction. For an expansion of ${R}_ {\rm{X}}(t)$ in terms of $M$ exponentials plus a constant, $p\, \tilde{R}_ {\rm{X}}(p)$ will be the ratio of polynomials, both of degree $M$.

3.4.   The Padé-Laplace technique

The Padé-Laplace technique stems from the analysis by Yeramian and Claverie ^[8] with an extensive discussion by Claverie et al. ^[9].

The technique consists of deriving approximations as ratios of low-degree polynomials in the transform variable. Approximations of this kind are known as Padé approximants and characterised by $[M, N]$, the degrees of numerator and denominator respectively. The $[M, N]$ Padé approximant to a function $f(x)$ is the ratio of polynomials whose Maclaurin series agrees with the first $M+N+1$ coefficients in the expansion of $f(x)$. Conventionally, the polynomials are defined such that the constant term in the denominator is 1. Padé approximants can give good fits that converge more rapidly than the underlying series and can even extend series beyond their radius of convergence. The Padé-Laplace technique, as used here, consists of analysing response function relations in terms of Padé approximants to the Laplace transforms.

Some of the important aspects of the use of the Padé-Laplace technique for analysing response functions are:

● The various combinations of response function transforms can be expressed as ratios of polynomials if the inputs also have that form. However the degree of the resulting expressions may be higher than is justified by the knowledge of the system and lower-degree expressions from the Padé-Laplace technique may be more than adequate.

● The Padé-Laplace technique can be applied even for expressions that are not expressible as ratios of polynomials, e.g. those used by Enting ^[5].

● The calculations for combinations can be performed using the series expansions without having to express all intermediate expressions as ratios of polynomials.

● The technique is not restricted to the use of Maclaurin series, but can also use Taylor series expansions about some chosen $p_0$ as long as all input series are expansions about the same $p_0$.

● The choice of $p_0$ introduces an extra degree of freedom into the choice of approximation. The use of $p_0 = 0$ (i.e. Maclaurin series, as in the examples in ^[21]) is appropriate for the long timescales targeted for stabilising the human influence on climate. A choice of $p_0 \approx 0.02$ focuses the fits on the timescale of the growth in anthropogenic CO$_2$ emissions over the 20th century.

● More generally, it may be appropriate to impose additional constraints (and thus reduce the number of coefficients that are fitted) in order to reproduce known behaviour such as the limiting values. This possibility is not explored in this paper.

The aim of the Padé-Laplace formalism is to produce a "reduced model" description of a complex system. It needs to be noted that the choice of degrees of the approximant translates into a choice of the form of reduced model.

--- $[M-1, M]$: This case expands to a sum of partial fractions which transforms to a sum of exponentials. This represents the impulse response function of the reduced model which can, as noted above, be expressed as a set of differential equations.

--- $[M, M]$: This can be expressed as a constant, plus an approximant of degree $[M-1, M]$: and so the $[M, M]$ approximant transforms to a sum of exponentials plus a multiple of the Dirac $\delta$ function.

--- $[M, M-1]$: This can be expressed as a multiple of $p$, a constant and a sum of partial fractions. In the time domain the reduced model operation becomes differentiation, plus convolution with a sum of exponentials plus a multiple of the $\delta$ function.

3.4.1.   Precursors

An early case of approximating CO$_2$ responses in using low-degree Padé approximants in $p$ was given by Enting ^[5] in the analysis of the Box Diffusion Model of Oeschger et al. ^[22]. This work was not taken further at that time.

A more recent example of calculations in the transform domain was the analysis of greenhouse gas emission equivalence by Enting and Clisby ^[13]. This was based on the Forcing-Equivalent Index (FEI) defined by Wigley ^[23] which has emission histories equivalent if they have equal radiative forcing at all times. With ${a_ {\rm{CH4}}}$ and ${a_ {\rm{CO2}}}$ giving the radiative forcing per unit mass of CH$_4$ (methane) and CO$_2$ respectively, the FEI equivalence of CH$_4$ and CO$_2$ requires

which gives the FEI CO$_2$-equivalent of CH$_4$ emissions as

where $\left[ \tilde{R}_ {\rm{CO2}}(p) \right]^{-1}$ corresponds to the CO$_2$ deconvolution operator discussed below. Enting and Clisby discussed various proposed replacements of GWP in terms of how their behaviour matched $\left[ \tilde{R}_ {\rm{CO2}}(p) \right]^{-1} \times \tilde{R}_ {\rm{CH4}}(p)$ over the range of timescales, $p$.

4.   Carbon cycle relations

4.1.   The concentration feedback

A highly influential analysis by Friedlingstein et al. ^[24] expressed CO$_2$ changes as

with the factors $\beta_ {\rm{L}}$ and $\beta_ {\rm{O}}$ regarded a describing a feedbacks from concentration changes affecting land and ocean systems respectively. The term from temperature influence is given explicitly in equation (5.1) below.

Equation (4.1) represents emissions that are partitioned between atmosphere, land biota and oceans in the ratio: $1 :\beta_ {\rm{L}}:\beta_ {\rm{O}}$. An analysis by Oeschger et al. ^[25] shows that for exponentially growing emissions, these ratios depend on the growth rate, $c$, and should be written as the ratio: $1 :\beta_ {\rm{L}}(c):\beta_ {\rm{O}}(c)$. The relation (3.11) above indicates how responses specific to pure exponential growth can be interpreted as Laplace transform relations. These take the form

with $\tilde{R}_ {\rm{L}}(p)$ and $\tilde{R}_ {\rm{O}}(p)$ given by

The terminology of what is classified as feedbacks depends on what is regarded as the underlying system versus what is treated as external to the system. For a "null" case of the atmosphere only, accumulating net CO$_2$ fluxes, the "null" response has $\tilde{R}_ {\rm{null}}(p) = 1/p$. Wigley and Raper ^[26] effectively regarded the "base" carbon as atmosphere + ocean (characterised by $R_ {\rm{O}}(t)$) and referred to the land contribution as a feedback. Relative to a system of atmosphere + ocean + land, the influence of climate counts as a feedback. Recent usage has included the climate-to-carbon influence as part of the carbon system and defined GWPs in terms of $R_ {\rm{FB}}(t)$ ^[1].

The connection between an impulse response description and a description in terms of feedbacks, as shown by (4.2), seems to be little appreciated. For example, Gregory et al. ^[27] state that "The concentration to carbon feedback is negative; it has generally received less attention in the literature" [i.e. less than the temperature-to-carbon feedback]. Discussions of the carbon cycle in terms of concentration feedbacks generally ignore extensive studies of using response functions, e.g. ^[1], in connection with the definition of Global Warming Potentials.

4.2.   Combining feedbacks

Combining (4.3) and (4.2) gives

Enting et al. ^[6] gave a less direct derivation of this relation. This was not used numerically, although if $R_ {\rm{L}}$ and $R_ {\rm{O}}$ are expressed as sums of $M$ and $N$ exponentials respectively an explicit expression for $R$ as a sum of $M+N$ exponentials could be obtained. This level of detail is unlikely to be justified by the accuracy of estimates of $R_ {\rm{L}}$ and $R_ {\rm{O}}$. The Padé-Laplace formalism provides a systematic way of simplifying such expressions.

The Laplace transform also provides a convenient way of combining impulse response functions for subsystems into an overall system description. For example the study by Enting ^[5] shows how the carbon response for the ocean mixed layer can be incorporated into a linearised version of the Box Diffusion Model of Oeschger et al. ^[22]. These expressions can also be used with other estimates of the mixed layer response, e.g. ^[3]. Again, the Padé-Laplace technique provides a practical way of implementing such relations.

4.3.   Padé-Laplace analysis

To illustrate the Padé-Laplace technique we construct explicit representations of $\tilde{B}(p)$, using (4.2). The input is the response function "Iinit", expressed as a sum of 5 exponentials, taken from ^[6] and listed in Table 1. This case has been analysed with initial calculations using Maclaurin series ^[21]. If $R$ is expressed as a sum of $M$ exponentials, $p\tilde{R}(p)$ can be expressed exactly as an $[M-1, M-1]$ approximant, as can $\tilde{B}(p)$. $\tilde{B}(0) \neq 0$, as expected since it characterises long-term carbon partitioning into the oceans and land biota.

Figure 1 shows some examples of fitting $\tilde{B}(p)$, the Laplace transform of the concentration feedback. The solid line is $\tilde{B}(p)$ (a ratio of two 4th degree polynomials) derived from $R_ {\rm{Iinit}}$ from ^[6] using (4.2). The $[2, 3]$ Padé approximant fitted to the Maclaurin series of $\tilde{B}(p)$.

but is virtually indistinguishable from the full expression for $\tilde{B}(p)$ ^[21] and is not shown here. The dashed line is the $[1, 2]$ Padé approximant fitted to the Maclaurin series of $\tilde{B}(p)$ and gives a reasonable approximation. The dotted line and chain lines are the $[1, 2]$ and $[0, 1]$ Padé approximants fitted to the Taylor series of $\tilde{B}(p)$ expanded around $p = 0.02$. The $[1, 2]$ approximant (dotted line) gives a better approximation than the $[1, 2]$ approximant to the Maclaurin series. Even the $[0, 1]$ approximant, given by

gives a reasonable approximation to $\tilde{B}(p)$ except at very small $p$.

This analysis also forms the basis for analysing the deconvolution operator that inverts (3.1) to derive sources from concentrations. This can be expressed as

This gives a specific computational form of the generic relations from Enting and Mansbridge ^[4], confirming that differentiation, with its associated numerical noise amplification, is an inherent characteristic of this deconvolution. The various ways in which the second term can be expressed (cf. equation (3.6)) reflect alternatives given in ^[4].

5.   Climate

5.1.   Response and coupling to carbon

Friedlingstein et al. ^[24] expressed CO$_2$ changes as

(where equation (4.1) above considered only the first term) and

leading to a gain

These relations have provided a framework for characterisation of both observations and model studies such as intercomparisons ^[28,29].

The discussion in previous sections suggests that (5.1), (5.2) and (5.3) should be generalised to incorporate a dependence on timescale. This can be expressed as Laplace transform relations with $\alpha, \beta, \gamma$ replaced by $\tilde{A}(p), \tilde{B}(p), \tilde{\Gamma}(p)$, (the corresponding upper case Greek letters) as is done in relation (5.11). The Laplace transform analysis comes from Enting ^[7].

An impulse response for climate, giving the warming $W(t)$ as a response to radiative forcing, $F(t)$ can be written as

where $U(t)$ is an impulse response function characterising how global temperature responds over time to an instantaneous change in radiative forcing. Such impulse response expressions have been widely used (see ^[30,31] for example).

In describing climate-to-carbon feedback, Enting ^[7] split off a linearised form of the radiative forcing from CO$_2$, (as is done in (5.2)) expressing the Laplace transform of (5.4) as

where $F'(t)$, with transform $\tilde{F}'(p)$, denotes the radiative forcing from influences other than CO$_2$. Thus $\alpha$ generalises to

The coupling back to the carbon cycle was expressed as

where $H$ characterises the carbon flux resulting from warming (or other temperature change). This leads to

with gain given by

with $\tilde{H}(p) = - p \tilde{\Gamma}(p)$, as a generalisation of (5.3). The relation between the concentration and temperature feedbacks can be simplified by expressing the response as

illustrating the property that the gain from multiple linear feedbacks is additive.

5.2.   Application

Equation (5.7) can be written as

thus generalising (5.1) by incorporating a timescale, $p$, (with integration of the emissions, $S(t)$, transformed by multiplying the transform by $1/p$).

While (5.1) has lead to a number of studies (e.g. those cited in ^[19]) that sought to estimate $\gamma$ as a single parameter, independent of timescale, two more recent studies ^[18,19] incorporated a parameterised timescale dependence. The estimation by Rubino et al. ^[18] was based on the regression equation, from the time domain form of (5.7):

The analysis covered the pre-industrial period and so $S(t)$ was set to zero. The temperature series came from various paleo-temperature estimates.

$H$ was parameterized as

or

where $\gamma'$ is a scale factor defining the size of the initial pulse response and $\tau$ is a relaxation time for the effect to dissipate. A time constant of $\tau = 100$ years was used, guided by the timescales used for the climate-to-carbon influence in the MAGICC model ^[2]. The constant term $\theta_0$ is required in the regression both to counter the residual effects of initial conditions and to resolve the ambiguity in consistency of reference levels of $M(t_0)$ and the zero of $W(t)$. By fitting $M(t)$ rather than $Q(t)$, $\theta_0$ becomes an estimate of the "pre-industrial" CO$_2$ level implied by the zero of each particular paleo-temperature estimate.

Enting and Clisby ^[19] extended the analysis into the industrial period by including an additional term in the regression equation:

with $R(t)$ approximated by $R_ {\rm{Iinit}}(t)$ (see Table 1). The factor $\theta_s$ was included as an estimate of $1-G$ to account for the fact that $R_ {\rm{Iinit}}(t)$ was arguably closer to $R_ {\rm{FB}}$ than to $R$. The estimates of $\theta_s$ varied significantly depending on which temperature data set was fitted and on how much of the pre-industrial era was included in the fit. As noted below, parameterized representations of $1-G$ could provide a basis for refining the analysis that was based on a single parameter $\theta_s$. Another way in which the regression analysis from ^[19] should be refined is through refining the assumption of white noise error in the CO$_2$ data as is implicit in regression analysis of (5.12). Refined analysis needs to also include the noise structure in the paleotemperature data, especially given the finding of long-tailed (fractal) relations between CO$_2$ and temperature ^[32].

Both the study of pre-industrial times ^[18] and the study that included the industrial period ^[19] produced a wide range of estimates of $\gamma'$, depending on the time period that was analysed and the data set that was fitted. This may reflect the conclusion by Neukom et al. ^[33] that there were no globally coherent pre-industrial climate fluctuations in the last millennium.

5.3.   Decomposition

When response functions, or the models from which they are derived, are calibrated against 20th century data, they will be implicitly including the effects of climate-to-carbon feedback ^[15]. Thus many results will be estimates of $R_ {\rm{FB}}$ rather than $R$. In order to separate the effects, ^[21] took cases from Joos et al. ^[1,34] who gave an example where the feedback had been explicitly excluded. This made it possible to combine estimates of $\tilde{R}_ {\rm{no feedback}}$ and $\tilde{R}_ {\rm{feedback}}$ to give an estimate of the loop gain, $\tilde{G}(p)$.

This is shown as the solid line in Figure 2. The negative gain at very long timescales was an artifact of the original fits only using 1000 years of model output. Nevertheless, the [1, 2] approximant could follow the function from an unphysical starting point and give reasonable fits over the more relevant decadal to century timescales as could the $[2, 3]$ approximant ^[21].

In order to fit an approximation to the timescales actually fitted in ^[34], approximants were fitted to Taylor series expansions about $p = 0.02$. The $[0, 1]$ approximant is shown as the dotted line in Figure 2 and is given by

This produces an excellent fit over a wide range of timescales and extrapolates to $p = 0$ in a way that avoids the artifacts arising from the original fits of $\tilde{R}_ {\rm{no feedback}}$ and $\tilde{R}_ {\rm{feedback}}$. It could provide an improved representation of the factor $1-G$ that was fitted by the constant $\theta_s$ in (5.15) as used in ^[19]. Of particular note is the increase by a factor of 3 for the gain on century timescales compared to decadal timescales. Higher order approximants to the Taylor series show increasing amounts of influence from this anomalous small-$p$ behaviour of $\tilde{G}(p)$.

6.   Discussion: comparison with other approximations

6.1.   Framework

The expression of carbon-climate feedback in terms of $\alpha, \beta, \gamma$, as in equation (5.3) from ^[24] has proved highly influential. As described above, it can be generalised to include the dependence on timescale through $\tilde{A}(p), \tilde{B}(p), \tilde{\Gamma}(p)$. Some of the studies that go beyond treating all of $\alpha, \beta, \gamma$ as constant are noted in the following sections.

6.2.   Beta

Expressions such as (5.1) treat $\beta$ as a constant, rather than incorporating the more general timescale dependence of $\tilde{B}(p)$. Approximating $\tilde{B}(p)$ as a constant is equivalent to taking $p\, \tilde{R}(p)$ to be constant. This corresponds to assuming that the airborne fraction is constant, since in this case the asymptotic, cumulative and instantaneous airborne fractions are all equal. This corresponds to putting $R(t) \approx \kappa \, \Theta(t)$, where $\Theta(t)$ is the unit step function.

In 1983, with 25 years of direct measurements of CO$_2$ concentrations, Laurmann and Spreiter ^[35] suggested that for most purposes the approximation of a constant airborne fraction could be an adequate substitute for more complicated models. The exception was for deconvolution problems, as would be expected since such inverse problems are more sensitive to model error ^[4,36].

More recently, with 50 years of direct CO$_2$ observations and ice-core data extending the record to earlier times, the behaviour of the airborne fraction has been revisited. There were a series of disputed claims concerning whether or not the instantaneous airborne fraction was showing a significant medium term decrease and the possible reasons for such a change. In particular the question was whether a change indicated a non-linearity in carbon cycle response or whether it reflected departures from pure exponential growth over $(-\infty, t]$, and in particular the beginning of the industrial period. An overview by Raupach et al. ^[37] (which cites earlier stages of the discussion) concluded that the change reflected a combination of influences from slower than exponential CO$_2$ growth, influence of volcanic eruptions, sink responses to climate change and non-linearities in response to increasing CO$_2$.

The reason why the change in airborne fraction was hard to detect is that 20th century emissions grew with an $e$-folding time close to 50 years (32 year doubling). This has the consequence that observations primarily constrain knowledge of $\tilde{R}_ {\rm{X}}(p = 0.02)$ ^[15].

As discussed in more detail as a possible future direction (see Section 7 below) additional information can be obtained from analysis of data on the distribution of ${}^{14}$C (radiocarbon). As well as calculating the air:land:ocean carbon partitioning of $1 :\beta_ {\rm{L}}(c):\beta_ {\rm{O}}(c)$ for exponential growth at rate $c$, Oeschger et al. ^[25] calculated corresponding partitioning of perturbations in ${}^{14}$C. Similarly, Enting ^[5] explicitly considered ${}^{14}$C in terms of Laplace transforms. Such studies could give additional information regarding how $\tilde{B}(p)$ is partitioned into contributions from land, $\tilde{B}_ {\rm{L}}(p)$, and ocean, $\tilde{B}_ {\rm{O}}(p)$. The analysis of ${}^{14}$C responses is outside the scope of the present paper.

6.3.   Gamma

As noted in Section 5.2 above, a number of studies (e.g. those cited in ^[19]) sought to estimate $\gamma$ as a single parameter without any dependence on timescale. An early study by Woodwell et al. ^[38] had noted that the temperature influence on CO$_2$ was greater at decadal to century timescales than at interannual or millennial timescales. The relevant comparisons with ^[38] are in terms of

These are shown in Figure 3 for cases with with $- \, \tilde{\Gamma}(p) = \gamma' \frac{p}{p+1/\tau}$ for $\tau = 100$ and $\tau = 50$ and $1+ \tilde{B}(p)$ derived from Iinit as above. The case with $\tau = 50$ years shows less timescale dependence than found in the older study by Woodwell et al. ^[38]. Similarly, Enting and Clisby ^[19] used a Laplace transform analysis to revisit the estimates based on exponential smoothing by Bauska et al. ^[39] and concluded that the results were broadly consistent with $\tau \approx 125$ years.

Wang and Nemani ^[40] proposed a simplified model that linked the relaxation of responses to concentration and temperature change, thus linking the parameterisation of $\tilde{B}(p)$ and $\tilde{\Gamma}(p)$ as generalisations of $\beta$ and $\gamma$.

For the pre-industrial period, with $S(t) = 0$, equation (5.7) can be transformed to

relating the fluxes from temperature to the deconvolution of concentration data. This approach was used by Bauska ^[41].

7.   Conclusions

The examples presented here illustrate the power of the Padé-Laplace technique in analysing the carbon cycle and its coupling to climate change. The Laplace transform representation captures the dependence on timescales in a way that is often neglected in simple parameterisations. As noted in connection with $\beta$ and the airborne fraction, the neglect of timescale dependence will give most distortion in inverse problems such as those involved in fitting observational data. Thus, appreciation of the timescale dependence in the gain of the carbon-climate feedback loop will be important in relating model behaviour to current observations. Approximations using the Padé-Laplace technique provide a systematic computational basis for assessing the limitations of such simple parameterisations.

An important extension of this work would be to include an analysis of radiocarbon (${}^{14}$C). Before the availability of CO$_2$ concentration data from ice cores, much of the information about the carbon cycle came from radiocarbon data ^[42]. Broadly speaking, samples that showed the effects of atmospheric nuclear testing constrained the short timescale behaviour of the carbon cycle while material that had been isolated from the post-1950 atmosphere constrained long-term behaviour on timescales comparable to the ${}^{14}$C half-life. As noted above, Oeschger et al. ^[25] calculated growth-rate dependent expressions for ratios $1 :\beta_ {\rm{L}}(c):\beta_ {\rm{O}}(c)$, and these can be interpreted as Laplace transforms (see equation 3.11)). Oeschger et al. also calculated comparable expressions for the partitioning of radiocarbon. Enting ^[5] determined a number of radiocarbon responses of the box diffusion model, expressed as Laplace transforms. The use of ${}^{14}$C data may require extensions to the approach described here, and involve fitting Padé approximants with additional constraints beyond fitting a single Taylor series. In addition, simple expressions as sums of exponentials may be useful in characterising uncertainties.

Clearly, a major limitation of this analysis is that it is restricted to linear systems, while non-linearities in the earth system are becoming increasingly apparent. A recent analysis by Leach et al. ^[31] explores the possibility of analysing non-linear behaviour using impulse response functions with a scenario-dependent transformation of the timescale. Such an approach could extend the applicability of the analyses presented in the present paper.

Notation

Where appropriate, functions are listed in the form $f(t), \tilde{f}(p)$. Upper case Greek such as $\Gamma(t)$, and also $A(t)$ and $B(t)$, are used to represent impulse response functions, with $\tilde{\Gamma}(p)$, $\tilde{A}(p)$ and $\tilde{B}(p)$ as the respective Laplace transforms. The corresponding lower case Greek, $\gamma$, $\alpha$, $\beta$ respresent constants that approximate $\tilde{\Gamma}(p)$, $\tilde{A}(p)$ and $\tilde{B}(p)$. Units here are scaled to express CO$_2$ in terms of mass of carbon. Some other studies have used the same or similar symbols with CO$_2$ expressed in terms of mixing ratios.

$a_j$ Amplitude of $j$th term in expression of response as sum of exponentials.

$a_ \text{CO2}$ Radiative forcing per unit of additional CO$_2$ expressed as mass of carbon, and using a linearisation of logarithmic forcing relation.

$f(t), g(t), h(t), k(t), \tilde{f}(p), \tilde{g}(p), \tilde{h}(p), \tilde{k}(p),$ Generic functions used to express generic properties of the Laplace transform.

$F'(t), \tilde{F}'(p)$ Radiative forcing from non-CO$_2$ atmospheric constituents (and other changes).

$G(t), \tilde{G}(p)$ Response function specifying the gain around the carbon-climate feedback loop.

$H(t), \tilde{H}(p)$ Response function specifying the CO$_2$ flux, caused by temperature change, into Land and Ocean systems.

$M_ \text{CO2}$ Atmospheric CO$_2$ content, in Gigatonnes (Pg) of carbon.

$p$ Transform variable.

$Q(t), \tilde{Q}(p)$ Perturbation in atmospheric CO$_2$ content from notional pre-industrial level or other reference amount.

$R_ \text{X}(t), \tilde{R}_ \text{X}(p)$ Generic atmospheric CO$_2$ response. Specific cases are $R_ {\rm{L}}, R_ {\rm{O}}, R_ {\rm{FB}}$.

$t$ Time, notionally in years.

$U(t), \tilde{U}(p)$ Response defining warming from increase in radiative forcing.

$W(t), \tilde{W}(p)$ Temperature change ("warming"), relative to reference level.

$\alpha$ Parameterisation, as a single number, of the temperature increase from a change in CO$_2$.

$A(t), \tilde{A}(p)$ Response function giving the temperature change from a change in atmospheric CO$_2$ content. $A(t) = a_ {\rm{CO2}} U(t)$.

$\beta_ \text{L}, \beta_ \text{O}, \beta$ Parameterisation, as a single number, of the CO$_2$ exchange caused by atmospheric CO$_2$ change, for Land and Ocean systems and the combination $\beta = \beta_ {\rm{L}} + \beta_ {\rm{O}}$.

$B_ \text{L}(t), B_ \text{O}(t), B(t), \tilde{B}_ \text{L}(p), \tilde{B}_ \text{O}(p), \tilde{B}(p)$ Response function specifying the CO$_2$ exchange caused by atmospheric CO$_2$ change, for Land and Ocean systems and the combination $B(t) = B_ {\rm{L}}(t) + B_ {\rm{O}}(t)$.

$\gamma_ \text{L}, \gamma_ \text{O}, \gamma$ Parameterisation, as a single number, of the CO$_2$ exchange caused by temperature change, for Land and Ocean systems and the combination $\gamma = \gamma_ {\rm{L}} + \gamma_ {\rm{O}}$.

$\Gamma_ \text{L}(t), \Gamma_ \text{O}(t), \Gamma(t), \tilde{\Gamma}_ \text{L}(p), \tilde{\Gamma}_ \text{O}(p), \tilde{\Gamma}(p)$ Response function specifying the CO$_2$ exchange caused by temperature change, for Land and Ocean systems and the combination $\Gamma(t) = \Gamma_ {\rm{L}}(t) + \Gamma_ {\rm{O}}(t)$.

$\gamma'$ Scale factor used in the parameterisation of $H$, giving the size of the initial pulse response to a change in temperature, $W$.

$\theta$ Coefficients in regression fit. Specific cases $\theta_0$, $\theta_\gamma$, $\theta_s$ ^[19].

$\Theta(t)$ Unit step function.

$\kappa(p)$ Asymptotic airborne fraction as a function of $e$-folding rate, $p$.

$\lambda_j$ Growth rate of $j$th term in expression of response as sum of exponentials.

$\tau$ Recovery time parameter used in parameterisation of $H(t)$ ^[18,19].

Acknowledgments

The author wishes to thank Nathan Clisby and Cathy Trudinger for valuable comments on the manuscript. Helpful comments from the two referees are also acknowledged. The numerical calculations and plotting of the figures was performed using R, with the pade package from Avram Adler.

Conflict of interest

The author declares that he has no conflict of interest regarding the work in this paper.