Impulse response functions of terrestrial carbon cycle models: method and application

To provide a common currency for model comparison, validation and manipulation, we suggest and describe the use of impulse response functions, a concept well‐developed in other fields, but only partially developed for use in terrestrial carbon cycle modelling. In this paper, we describe the derivation of impulse response functions, and then examine (i) the dynamics of a simple five‐box biosphere carbon model; (ii) the dynamics of the CASA biosphere model, a spatially explicit NPP and soil carbon biogeochemistry model; and (iii) various diagnostics of the two models, including the latitudinal distribution of mean age, mean residence time and turnover time. We also (i) deconvolve the past history of terrestrial NPP from an estimate of past carbon sequestration using a derived impulse response function to test the performance of impulse response functions during periods of historical climate change; (ii) convolve impulse response functions from both the simple five‐box model and the CASA model against a historical record of atmospheric δ13C to estimate the size of the terrestrial 13C isotopic disequilibrium; and (iii) convolve the same impulse response functions against a historical record of atmospheric 14C to estimate the 14C content and isotopic disequilibrium of the terrestrial biosphere at the 1° × 1° scale. Given their utility in model comparison, and the fact that they facilitate a number of numerical calculations that are difficult to perform with the complex carbon turnover models from which they are derived, we strongly urge the inclusion of impulse response functions as a diagnostic of the perturbation response of terrestrial carbon cycle models.


Introduction
The global carbon budget has changed significantly due to anthropogenic increases in atmospheric carbon dioxide (Keeling et al. 1989), and many suspect this increase may have wide-ranging effects on global climate as well as on basic terrestrial physiological and hydrological processes . But the magnitude of the potential change is uncertain (Moore & Braswell 1994;Sarmiento et al. 1995;Joos et al. 1996) as well as confounded by uncertainties in the scope of human activity (Houghton & Meira Filho 1995;Schimel 1995). A major question is whether the oceans or the terrestrial biosphere can mediate the recent increases in atmospheric carbon dioxide by increasing their carbon sequestration rates.
There is considerable evidence now that the terrestrial biosphere is, indeed, a significant sink of atmospheric carbon, as seen in analyses of atmospheric and ocean data (Tans et al. 1990;Ciais et al. 1995a;Ciais et al. 1995b), as well as in studies of the sensitivity of production to increased atmospheric CO 2 (Farquhar & Sharkey 1982;Farquhar & Wong 1987;Luo & Mooney 1995;Wullschleger et al. 1995;Amthor & Koch 1996;Luo et al. 1996) or nitrogen availability (Field & Mooney 1986;McGuire et al. 1992;McMurtrie & Wang 1993;Townsend et al. 1996;Holland et al. 1997;Vitousek et al. 1997a;Vitousek et al. 1997b). If production is actually increasing, it may lead to significant increases in carbon storage if the turnover of carbon is sufficiently slow (Taylor & Lloyd 1992;Friedlingstein et al. 1995;Lloyd & Farquhar 1996; Thompson et al. 1996;Townsend et al. 1996). The potential of a given ecosystem to maintain a significant carbon sink is given by the system's sink potential, or the product of NPP and the turnover time of the system, Pτ (Taylor & Lloyd 1992; Thompson et al. 1996), where P is production and τ is the turnover time.
Carbon turnover, however, unlike photosynthesis, is less well understood, due largely to its hidden, complex and heterogeneous nature (see review by Stevenson 1994). Soil chemical transformations are controlled by many things, including temperature, texture, moisture, biota, relief, and nutrient availability (Jenny 1980). Some soil carbon models, such as CENTURY (Parton & Scurlock JMO Ojima 1993;Schimel et al. 1994), Rothamsted (Jenkinson 1990;Jenkinson et al. 1991), Biome-BGC (Running & Coughlan 1988;Running 1990;Running & Hunt 1993) and FBM (Kindermann & Lü deke MKB Badeck 1993;Lü deke & Otto 1994;Lü deke & Otto 1995), have tackled these controls. For the sake of convenience, these models consign various classes of carbon compounds to a number of plant, detrital and soil carbon pools. There is little support for one model pool structure over another, and in many cases, this has inspired simplification of the models, reducing the controls to soil texture and climate (Schimel et al. 1994).

The impulse response approach
An impulse response function, defined as the system response to an impulse as a function of time since the impulse, can be derived from any model that is sufficiently linear and time-invariant (Nir & Lewis 1975;Lewis & Nir 1978). The rudiments of impulse response functions can be found in some of the very earliest literature analysing the carbon cycle, where the idea was to incorporate terrestrial functions into coherent, integrated and dynamic models of the entire global carbon cycle (Eriksson & Welander 1956;Eriksson 1971;Bolin & Rodhe 1973;Bolin 1975). They are also successful in describing the movement of tracers in atmospheric and ocean models (Siegenthaler & Oeschger 1978;Maier-Reimer & Hasselmann 1987;Sarmiento et al. 1992;Joos et al. 1996), and they have already been applied, at least in principle, to terrestrial carbon cycle models (Emanuel et al. 1981;Emanuel et al. 1984;Moore & Bolin 1986/1987Siegenthaler & Oeschger 1987;Moore & Braswell 1994;Joos et al. 1996).
Impulse response functions are powerful analytical tools. By convolving them against historical records of model tracer input (such as the concentration of atmospheric 13 CO 2 and 14 CO 2 , or against terrestrial NPP), one can efficiently calculate the isotopic disequilibria of 13 C or 14 C in the terrestrial biosphere (Enting et al. 1993;Enting et al. 1995;Fung et al. 1997), the 13 C and 14 C contents of various soil compartments (Balesdent 1987;Trumbore et al. 1995), or the controlling relationship between changes in NPP and the terrestrial carbon sink (Taylor & Lloyd 1992;Friedlingstein et al. 1995;Thompson et al. 1996;Townsend et al. 1996). Impulse response functions are not widely used, and they have not yet been applied to spatially resolved, complex models of terrestrial carbon turnover. Nevertheless, their potential utility in studies that explore the role of the terrestrial biosphere in the global carbon cycle is enormous (Joos et al. 1996).

This study
By demonstrating and refining this method we hope to clarify its use and encourage comparison and validation of complex carbon turnover models via their impulse response functions. We present the method and demonstrate its usefulness with respect to (i) the dynamics and diagnostics of continuous and discrete pool-based terrestrial carbon cycle models, (ii) the relative performance of an impulse response function compared to the model from which it was derived, and (iii) the utility of impulse response functions in calculating the reciprocal impacts of changes in atmospheric 13 C and 14 C on plant and soil isotopic content, as well as changes in biosphere function on atmospheric isotopic content.

Impulse response functions
Definitions. In practice, we find that different fractions of carbon leave the system at different times. The time it takes for a parcel of carbon to traverse a system is called its transit time or mean residence time, τ. Different parcels of carbon have different transit times. Each parcel also has an age, T, which is the time the carbon has spent in the system since sequestration via photosynthesis.
The steady-state distribution of carbon leaving the system as a function of transit time is called the impulse response of the system, Φ (τ). The integral of Φ (τ), with respect to τ, gives the total flux, Φ T , into and out of the system at steady state: (1) Note that Φ T does not vary with time, since Φ (τ) is a steady state function. Φ (τ) normalized by Φ T is the probability density function of transit times, τ, in the system: is also called Green's function, or kernel function, and is analogous to the 'unit hydrograph' used in groundwater precipitation response functions. φ (τ) can also be defined as (Lewis & Nir 1978): (i) the probability density function of transit times in the parcels of carbon entering the system, (ii) the outflow response of the whole system (in this case, the terrestrial biosphere) to an instantaneous input such as the Dirac delta function, or (iii) a weighting function describing the relative contribution of different inflows at times (t -τ) to the present outflow at time t. Inflow, in the case of the terrestrial biosphere, is gross (or net) primary production, while outflow is ecosystem (or heterotrophic) respiration. The first moment of φ (τ) with respect to τ is the mean transit time of carbon, τ: or, in other words, the average time required for a given parcel of carbon to traverse the biosphere from sequestration to eventual release through respiration. Each parcel of carbon in the biosphere has an associated age. The distribution of carbon in the terrestrial biosphere with respect to that age is the storage response of the system, Ψ(T), and can be alternatively defined as a weighting function giving the contribution of carbon sequestered at time (t -T) to the carbon presently stored in the system at time t. This function can be integrated with respect to T to give the total carbon storage of the terrestrial biosphere Ψ T at steady state: Ψ(T) normalized with respect to Ψ T is the probability density function of carbon storage with respect to T: and the first moment of ψ(T) with respect to T is the mean age, T, of carbon currently present in the biosphere: T and τ are equivalent when φ (τ) has an exponential form; that is, in a single-reservoir, first-order turnover model, the mean age and mean transit time are the same (Bolin & Rodhe 1973). An additional parameter, definable using terms described above, is the turnover time τ o ; it is usually expressed as the ratio of total carbon storage in the system to flux through the system: (7) Φ T at steady state. They subdivide steady-state systems into three different categories (Bolin & Rodhe 1973;Rodhe 1992): (i) T Ͻ τ ϭ τ o , in which carbon spends some finite time in the system before any appreciable fraction is respired (such as in the human population, where mean age may be around 30 but life expectancy is 70); (ii) T ϭ τ ϭ τ o , in which carbon is well -mixed and the probability of respiration remains constant irrespective of physical or temporal position within the system; and (iii) T Ͼ τ ϭ τ o , in which the probability of respiration declines the longer the carbon spends in the system (this is the situation most representative of that found in terrestrial ecosystems).
The age of each parcel, T, is the time since entering the system, while the transit time, τ, is the time required for a particle to traverse the system, irrespective of where it is found in the system. As soon as the particle departs the system, T ϭ τ. Making the distinction between the two is necessary for the development of the concepts presented here, but for practical purposes, it is necessary to switch to the sole use of τ for both τ and T. τ now represents the time since the particle entered the system: φ (τ) is the fraction of the initial impulse departing the system at time τ, and ψ (τ) is the fraction of the initial pulse remaining in the system at time τ after sequestration. T and τ will still be used for the mean age and mean transit time.
Additional diagnostics. At steady state, the quotient of the impulse response function and storage response function, k(τ), describes the probability of loss from the system as a function of age: Note the use of Φ(τ) and Ψ(τ), rather than φ (τ) and ψ (τ); in this case, we are interested in absolute carbon loss with respect to absolute carbon content. This function is called the decay function. k(τ) can be used only if the dynamics of the system do not change with time, i.e.
Indeed, this condition must hold for all functions derived using the impulse response method. The reciprocal of k(τ), τ o (τ), gives the instantaneous turnover time of carbon as a function of age τ Numerous Φ(τ) and Ψ(τ) functions from different cells within a spatially explicit model can be used as a whole by aggregating them with an area-weighted sum: where A x is the land area of cell x, and Φ (τ) and Ψ(τ) ( ( are the new aggregated impulse response and storage response functions, respectively. The functions can be converted to ( φ (τ) and ( ψ (τ) using (2) and (5).

Transformation of two models
The impulse response transformation: continuous vs. discrete. The transformation of carbon biogeochemistry models into their respective impulse response equivalents requires only that a pulse of carbon be applied to the model (e.g. through NPP, allocated to the different biomass pools) and followed through time as it is lost from the system through respiration. The fraction of the initial pulse remaining in the system after time τ is ψ (τ) while the fraction lost at time τ is φ (τ). In transformations of models with a discrete time step, impulse response functions can be defined discretely: In some models, the dynamics of carbon turnover change seasonally, a prohibited condition if we are using a single impulse response function. In such cases, it is necessary to define discrete, aggregated impulse response and storage functions that include loss from all subannum time steps into a single function. Thus, the time step of the model is forced to be annual. One problem with these aggregations is that many processes operate at time scales smaller than that chosen for the transformation time step such that carbon may be lost during the same time step as the pulse (i.e. carbon is respired immediately after sequestration, such as through plant respiration, at τ ϭ 0). This is especially important when linking changes in assimilation rate with changes in carbon storage. Operationally, Ψ τ is given by the carbon content of the system before the aggregated time step begins, and Φ τ is given by the carbon lost during the aggregated time step.
Single reservoirs within larger systems. The impulse response of individual reservoirs within a system may be calculated by monitoring them individually. But the carbon balance of a single pool, with respect to transit time and age, is slightly more complex than that of the entire system. First, the input to a single pool may be delayed since much of the carbon entering the pool may first pass through other pools, each with their own intrinsic turnover time. The input function of pool i can be defined in terms of the distribution of the carbon entering the pool with respect to its age, P i (τ), where age is defined with respect to the time since the carbon entered the system. We define two cases, one in which all carbon enters the pool with an age τ ϭ 0, called direct input, and the other where there is some disperse distribution of ages, called indirect input.
Second, carbon may depart the pool in one of two ways, by respiration or by transfer to another pool. The first case is denoted Φ r (τ), or the distribution of respiratory fluxes with respect to age (where age is the time since the parcel of carbon first entered the system), and the other is Φ t (τ), or the distribution of transfer fluxes to other pools with respect to age. It is useful to define a third function that gives the total loss from the system as a function of age: We can now write an equation for the mass balance of the reservoir as a function of τ: where dΨ(τ)/dτ is the change in carbon content of the pool as a function of age. For individual reservoirs, we must be careful how we define the fluxes. For instance, to calculate the mean transit time of a single reservoir, or the time required to traverse it, we use the first moment of the change in carbon content with age, divided by the total, time-integrated loss rate of the reservoir: For the mean age of carbon lost from the reservoir via respiration, we use: and for the mean age of carbon in the reservoir, we use: These distinctions are important when comparing model predictions with actual data. φ l (τ), φ r (τ) and φ t (τ) can be derived from Φ l (τ), Φ r (τ) and Φ t (τ), respectively, using (2).
NPP vs. GPP. Some models use NPP as their carbon input from the atmosphere (i.e. GPP minus autotrophic respiration), while others use just GPP. The difference in Fig. 1 The global terrestrial carbon turnover model of Emanuel et al. (1981). Values in boxes are the steady-state pool stocks in Pg C; values next to arrows are the annual flux rates from one pool to the next in Pg C y -1 . (a) The GPP-referenced version of the model where carbon is input via GPP. This model is described in eqns (A1-A3). Gross primary production (total plant assimilation or photosynthesis) is 113 Pg C y -1 . (b) The NPP-referenced model, derived in this paper and described in eqns (A 4-A6), with plant respiration removed and pool sizes of live biomass pools scaled to reflect overall change in input without a change in turnover time. Net primary production is 56 Pg C y -1 . In both models, the small flux (Ͻ 0.5) from the active soil carbon pool is ignored. model approach has implications for the shape of the impulse response function. If plant respiration, which is composed typically of very young carbon, is included in the impulse response function (as when using GPP as the carbon input), the distribution of fluxes will be skewed toward low residence times. Models that use NPP as their input include biomass pools with their own turnover times to introduce a delay between when carbon is assimilated and when it is delivered to the soil, but these pools do not respire carbon. In the absence of a clear consensus as to whether NPP or GPP is preferable, we make a distinction between the two whenever specific models are mentioned. Models that use GPP we call GPP-referenced and models that use NPP, NPP-referenced.
A simple box model of the terrestrial biosphere. To illustrate the impulse response method, we chose the simple 5-box model of global terrestrial carbon turnover developed by Emanuel et al. (1981). The model ( Fig. 1a) was chosen for four reasons: (i) it is simple, (ii) it is defined continuously, (iii) it is currently in use as a means of calculating the isotopic disequilibrium of 13 C (Enting et al. 1993;Francey et al. 1995), and (iv) it has been used as a benchmark for other models of carbon biogeochemistry (Wittenberg & Esser 1997). As such, it is an ideal tool for illustrating the impulse response concept as well as some of the pitfalls that may be encountered when using impulse response functions. The model of Emanuel et al. (1981) divides terrestrial carbon storage into five reservoirs, representing nontree ground vegetation, nonwoody tree parts, woody tree parts, detritus and active soil organic matter (Fig. 1a). It is a GPP-referenced model (indicated by the inclusion of live plant biomass pools with fluxes to the atmosphere) and is defined continuously on an annual basis by the steady state reservoir contents and annual flux rates between, into and out of pools (see Appendix). We calculate the impulse response function of this model by solving the model continuously and aggregating the solutions for each reservoir into a single storage response Fig. 2 The compartmentalized flow of carbon through live biomass, litter, and soil organic matter pools as described in the CASA biosphere model. For details, see Field et al. (1995), Randerson et al. (1996) and Thompson et al. (1996). function, Ψ(τ). Since Ψ(τ) is defined continuously, Φ(τ) is just the derivative of Ψ(τ). For each pool in the model, we also calculate ψ (τ), φ l (τ) and φ r (τ), as well as τ, T and τ o . Since each reservoir decays exponentially, φ l (τ) ϭ ψ (τ), and τ and T for each reservoir are equivalent. This is not the case for the model as a whole.
To test the difference between GPP and NPP-referenced models we convert this model into a NPP-referenced model by eliminating the plant respiration fluxes from the three biomass pools and adjusting the input fluxes accordingly, while maintaining their turnover times ( Fig. 1b; see Appendix). GPP (the gross annual flux of carbon from the atmosphere to the biosphere) in this model is 113 Pg C y -1 (Emanuel et al. 1981), while NPP (under our modifications) is 56 Pg C y -1 , a little less than half of GPP.
The CASA biosphere model. To illustrate the impulse response concept with a more complex, iterative model, we derived impulse response functions from CASA v1.2 (Fig. 2: Field et al. 1995;Randerson et al. 1996;Thompson et al. 1996;Malmströ m et al. 1997), a version in all ways identical to the model version used by Thompson et al. (1996). The CASA model is a global, spatially resolved model operating at the 1°by 1°scale on a monthly time step. All land pixels are modelled except those that are covered by ice (i.e. vegetation classes 1-12; Defries & Townshend 1994). The NPP model is driven by air temperature and precipitation data (Shea 1986), AVHRR NDVI data from 1990 with the FASIR correction (Los et al. 1994) 1990 monthly surface solar irradiance (Bishop & Rossow 1991), and a 1°by 1°soil texture classification (Zobler 1986). The soil carbon model is driven with carbon from the NPP submodel of CASA, as well as by air temperature, precipitation and soil texture. Global annual NPP under CASA using these data sets is 54.9 Pg C y -1 . This NPP is allocated to live roots, wood and leaves at a 1:1:1 ratio (or 1:1 ratio in systems with no wood: classes 7, 10 and 12), and the different biomass pools turn over on a first-order basis, constant over each of the vegetation classes defined by Defries & Townshend (1994). To account for the effects of herbivory on turnover, some carbon is shunted to the surface metabolic pool in proportion to annual NPP . The biomass turnover times were derived from an analysis by Kohlmaier et al. (1997), but autotrophic respiration is not calculated; this is an NPP-referenced model. Φ τ and Ψ τ were derived for each cell modelled by CASA on a discrete, annually aggregated basis, to a τ of 2000 years (it was not necessary, nor would it have been practical, to calculate Φ τ or Ψ τ to ϱ since Ψ τ becomes vanishingly small at high τ). Using (11) and (12), individual Φ τ 's and Ψ τ 's were globally aggregated on a land area weighted basis to provide global functions analogous to the functions derived from the model of Emanuel et al. (1981). The Φ τ and Ψ τ values of three points (centred at 2.5°S 59.5°W, 54.5°N 66.5°W, and 41.5°N 99.5°W) were saved for separate analyses (see below). Ψ τ is not defined for τ ϭ 0 since there is no carbon in the system before the pulse is applied, but this is not a concern since there is no respiration either.

Application of impulse response functions
Carbon sequestration in response to change in NPP. As Bolin et al. (1981) point out, with a time history of primary production and an impulse response function for the terrestrial biosphere, it is possible to calculate the time history of the sink. The gross flow of carbon from the biosphere to the atmosphere (the respiratory flow) is given by convolution integral of production against the normalized impulse response function: where P(t -τ) is the rate of primary production by the biosphere at some time (t -τ). The sink rate is then the difference between the incoming and outgoing fluxes: where N(t) is the net flux of carbon into the system at time t. φ (τ) can be used to calculate system respiration outside of steady state, since P(t -τ) does not need to be constant. Conversely, if N(t) is prescribed from some time t ϭ 0 to the present, we are given P at time t ϭ 0, and P does not vary as a function of time for all t Ͻ 0, then we can solve for the P(t) necessary to satisfy our prescribed N(t), either iteratively or with a Laplace transform. Thompson et al. (1996) reported that a long-term, 20% linear increase in global annual NPP over an initial value of 48.2 Pg C y -1 (a rate of increase of 0.1 Pg C y -1 ) was required for the terrestrial biosphere to sequester the µ 90 Pg estimated to have been sequestered from 1880 to 1990 (Houghton 1995;Sarmiento et al. 1995). Thompson et al. (1996) also concluded that compared to the magnitude of the NPP-dependent changes in carbon storage, changes in climate over the last 110 years had relatively little effect on long-term changes in respiration. As such, in the absence of significant deforestation or regrowth, change in production over that period would have to have been the primary control on carbon sequestration.
To test this assumption, as well as the CASA-derived, global impulse response function, we repeated their analysis, using the same initial conditions (48.2 Pg C y -1 NPP and steady-state carbon storage at and before 1880) and the same estimate for the terrestrial carbon sink from the study of Houghton (1995). Since it is impossible to incorporate changing climate into impulse response functions (due to the time invariant linear dynamics assumption), this exercise will also test the extent to which changing climate influences the impulse responsederived estimate of P(t).
The isodisequilibrium of 13 C. Changes in atmospheric 13 C content have been an important constraint in our estimates of the relative importance of oceans and the terrestrial biosphere as sinks for atmospheric CO 2 . Fossil fuel emissions, which are depleted in 13 C relative to the atmosphere, have diluted the atmosphere of 13 C (Fig. 3), creating a mismatch between carbon sequestered by the terrestrial biosphere some time in the past and carbon being sequestered now (Keeling 1979). This mismatch is called the isodisequilibrium, D 13 (Fung et al. 1997). It is formally defined as the 13 C flux from the atmosphere to the biosphere minus the 13 C flux from the biosphere to the atmosphere under conditions of steady state with respect to total carbon. It is commonly derived as the product of the gross flux of carbon between the biosphere and the atmosphere, P(t), and the difference between the δ 13 C of the atmosphere at time t and the δ 13 C of the atmosphere at some time (t -τ), where τ is the mean residence time of the biosphere (Tans et al. 1993;Ciais et al. 1995a;Ciais et al. 1995b). An important feature of the isodisequilibrium is that it increases with increasing τ.
Since the past-history of atmospheric δ 13 C is not a linear function of time and the biosphere does not follow simple exponential decay (as in the second category of Bolin & Rodhe 1973), it is more precise to calculate D 13 as the convolution of an impulse response function of the terrestrial biosphere against a past history of biosphereatmospheric 13 C flux: Fig. 3 Composite record of atmospheric δ A (t) from 1700 to the present using 1700-1953 ice core data (Friedli et al. 1986) 1978-81 atmospheric data (Keeling et al. 1989), and 1982-92 aggregate atmospheric data . The connecting line is a maximally rigid spline fit of the atmospheric observations and is used as a continuous estimate of past atmospheric 13 C in all analyses involving 13 C.
where δ A is the 13 C signature of the atmosphere, P is primary production (gross or net), and ε ab is the photosynthetic discrimination factor (usually -18‰ in most analyses; appropriate for C 3 plants). We assume ε ab to be constant through time, so it drops out of the equation. By eqn (21), D 13 is also a function of production. However, since D 13 is defined as the isotopic disequilibrium when total carbon is at steady state, P(t) can come out of the integral. Thus (21), for our purposes, can be rewritten: (22) Isodisequilibria are often reported divided by the gross production rate: For clarity, D 13 is the isodisequilibrium proper (a flux in units of Pg C‰ y -1 ) and D 13 (t)/P(t) is the isodisequilibrium forcing coefficient, ζ. The distinction between ζ 13 (equivalent to D b of Fung et al. 1997) and D 13 (equivalent to D b G b of Fung et al. 1997) is important: the term 'isodisequilibrium' has been used ambiguously for either one. D 13 is an index of the impact of the biosphere on the atmospheric 13 C budget, while the forcing coefficient is a normalized index of the isodisequilibrium, and has meaning for the atmospheric 13 C only in the context of its corresponding flux.
There are presently numerous estimates for the forcing coefficient, ranging from anywhere between 0.2‰ (Tans et al. 1993;Ciais et al. 1995b) and 0.23‰ (Enting et al. 1993) to 0.3‰ (Fung et al. 1997) and 0.41‰ (Wittenberg & Esser 1997). The wide variation in these estimates (Ͼ 50%) is the result of whether the model is GPP-or NPPreferenced. The kind of variation seen in the literature translates into large uncertainty with respect to estimates of the terrestrial carbon sink. We also calculate the isodisequilibria and forcing coefficients for the different models over the period from 1977 to 1992, as well the forcing coefficient of each cell modelled by CASA in 1992, comparable with plate 2 of Fung et al. (1997).
Calculating the 14 C content of stored and respired carbon. The 14 C content of stored and respired carbon has been shown in recent years to be a powerful tool for directly measuring rates of carbon turnover (Trumbore 1993;Trumbore et al. 1995). The impulse response method, and further revisions on the method, should allow experimenters and modellers to integrate these measurements with model performance. To use impulse response functions with radiocarbon data, it is necessary to construct a past history of atmospheric 14 C that takes into account the effect of age. The half-life of radiocarbon, 5730 years, is what dictates this correction.
The age and δ 13 C-corrected radiocarbon content of geophysical samples is given by the following relation (Stuiver & Polach 1977): where A SN is the activity of the sample corrected for δ 13 C, γ is 1/8267 y -1 (based on the 5730 years half-life), t is the year AD that the sample was created, and A ON is the δ 13 C corrected activity of the 1950 oxalic acid standard. The 14 C data that makes radiocarbon a useful tool for exploring plant and soil carbon dynamics is of recent origin, created in the atmosphere as a result of post WWII bomb testing. These data require no age correction. Older data, however, do require some correction assuming a constant, prebomb atmospheric ∆ of 0‰. For the northern hemisphere, we use the Georgian wine data of Burcholadze et al. (1989) from the beginning of the bomb era to 1985; data from 1986 to 1997 are derived from an assumption of a 9.7‰ per year decrease in ∆ estimated from trends in the data of Nydal & Lö vseth (1996) and Record and estimates of atmospheric ∆ A from 1920 to 1997 in the northern hemisphere (solid line) and the southern hemisphere (dotted line). All values are relative to a prebomb, age and δ 13 C corrected standard, which, by definition, has a ∆-value of 0‰. Prebomb era data are estimated assuming a steady state prebomb atmospheric 14 C content and a 14 C half-life of 5730 years. Levin et al. (1985) (S. Trumbore, pers. comm.). Southern hemisphere data up to 1993 are from the study of Manning & Melhuish (1994). 1994-97 southern hemisphere data are given the same values as the northern hemisphere (Fig. 4).
With the age correction, we calculate two parameters, the radiocarbon content of carbon respired from the terrestrial biosphere, ∆ R , and stored in the biosphere, ∆ C : The isotopic disequilibrium forcing coefficient of 14 C, ζ 14 , as in (23), is ∆ R minus ∆ A : ∆ C and ∆ R can be compared against bulk soil and respiratory measurements and against measurements of individual pools within a model (but only when the pool represents a measurable quantity). ζ 14 could be used in budgets of atmospheric 14 C, as for 13 C.

Site analysis
To demonstrate the usefulness of impulse response functions in summarizing important components of carbon turnover in different ecosystems, we perform the impulse response transformations of three cells in the CASA model in some depth. These are the same three cells analysed by Thompson et al. (1996) in their point analysis: a tropical forest cell centred at 2.5°S, 59.5°W, a boreal forest cell centred at 54.5°N, 66.5°W, and a grassland cell centred at 41.5°N, 99.5°W. A number of parameters, such as the mean age of stored carbon and the mean residence time, are calculated for each cell, but individual carbon pools are also analysed. Using impulse response functions for each pool in the system, we calculated for each the mean age of stored carbon, the transit or turnover time, and ∆ C (which on a pool by pool basis is the same as ∆ R ).

Results and discussion
Impulse response transformation of a simple box model Ψ(τ), Φ l (τ) and Φ r (τ) were calculated for the entire GPPand NPP-referenced Emanuel et al. (1981) models (Fig. 5) as well as for each pool within the GPP-reference model (Fig. 6). Since total system losses were solely respiratory, Φ l (τ), Φ r (τ) and dΨ(τ)/dτ for the whole system at τ Ͼ 0 were all equivalent (eqn 15). Since each pool in the model decayed exponentially, and because Φ r (τ) is a constant fraction of Φ l (τ), Φ r (τ), Φ l (τ) and Ψ(τ) were related to each other by a simple scalar. The major differences between the NPP-and GPPreferenced models were: (i) higher Ψ(τ) for low τ (up to τ ϭ 72 years) in the GPP-referenced model due to the fact that GPP is a larger pulse than NPP; and (ii) higher Ψ(τ) for high τ in the NPP-referenced model since plant respiration was eliminated, leaving older carbon cohorts to accumulate more carbon than their GPP counterparts. Furthermore, the change in slope of the log-y plots of Φ(τ) and Ψ(τ) (Fig. 5) shows that turnover rates varied with age: as the carbon shifted between different reservoirs, the slope was reduced. Ψ(τ) was calculated for each pool of the GPP-referenced model ( Fig. 6; see Appendix). Shifts in pool dominance were observed as the older pools accumulated carbon at the expense of younger ones (Fig. 6), consistent with the overall decline in system decomposition rate with age (as seen in Fig. 5). The majority of the carbon over the age of 10 years was found in woody tree parts or the active soil carbon reservoir, both of which, under the Emanuel model, decay slowly and contain the majority of total system carbon storage at steady state.  (Fig. 1a,b). Carbon storage Ψ(τ) and respiration Φ(τ) following a unit pulse of carbon. The solid line is Ψ GPP (τ), the dotted line is Φ GPP (τ), the long dashed line is Ψ NPP (τ), The nonconvergence between Ψ GPP (τ) and Ψ NPP (τ) and between Φ GPP (τ) and Φ NPP (τ) is due to small increases in the turnover time of the nontree ground vegetation and nonwoody tree parts pools with conversion to an NPPreferenced model. The insert is an expanded view of the first five years. As mentioned previously, each pool in the model of Emanuel et al. (1981) decays exponentially; thus, by definition, φ (τ) and ψ (τ) are equivalent. In fact, when any pool i is homogeneous, φ i (τ) ϭ ψ i (τ). This has a number of important implications. First, by (3) and (6), the mean transit time and mean age of each pool will be the same (Table 2); and second, assuming there is no significant discrimination for or against carbon isotopes in the decomposition process, the carbon stored in and © 1999 Blackwell Science Ltd., Global Change Biology, 5, 371-394 respired from the pool will be isotopically equivalent (eqn 23). The differences between the mean transit times of the different reservoirs (eqn 16) and their respective mean ages (eqn 17, Table 2) were the result of the additional time carbon spent in other pools before entering the pool in question (P(τ) is nonzero for τ Ͼ 0 (eqn 15). Only those pools in which carbon entered directly from the atmosphere (direct input) possessed equivalent turnover times and mean ages.

Φ(τ)
Steady state distribution of respired carbon with respect to transit time; impulse response function Φ T Steady state input and output flux rate between the system and its surroundings; (1) φ (τ) Probability density function of transit times; time-dependent outflow response to an instantaneous input; weighting function describing the relative contribution of different influxes at times (t -τ) to present outflow at time t; Green's function; kernel function; normalized impulse response function; (2) Ψ(T) Steady state distribution of carbon in the system with respect to age; storage response function Ψ T Steady state carbon storage in the system; (4) ψ(T), ψ (τ) Probability density function of carbon storage with age; normalized storage response function; (5) Φ x (τ) Impulse response function of cell x; (11) Ψ x (τ) Storage response function of cell x; (12) A x Surface area of cell x; (11) and (12) Φ (τ) Spatially aggregated impulse response function; (11) ( Ψ (τ) Spatially aggregated storage response function; (12) ( φ τ Discretized version of φ (τ); (13) ψ τ Discretized version of ψ (τ); (13) P i (τ) Distribution of carbon entering pool i as a function of time since sequestration (if all carbon enters directly from the atmosphere, P i (τ) ϭ δ(τ), the Dirac-delta function); (15) Φ r (τ) Steady state distribution of respiratory carbon loss from the reservoir as a function of time since sequestration; respiratory impulse response function; (14) Φ t (τ) Steady state distribution of carbon leaving the reservoir via transfer to another reservoir as a function of time since sequestration; transfer impulse response function; (14) Φ l (τ) Steady state distribution of all carbon leaving the reservoir (via either respiration or transfer) as a function of time since sequestration; Φ l (τ) ϭ Φ r (τ) ϩ Φ t (τ); (14) GPP Gross primary production; the rate of photosynthetic uptake of CO 2 from the atmosphere by plants NPP Net primary production; the difference between GPP and plant respiration GPP-referenced Input of carbon to system is via GPP and all carbon found in the system is GPP-derived NPP-referenced Input of carbon to system is via NPP and all carbon found in the system is NPP-derived k(τ) Decay function; instantaneous first-order decay constant of carbon in the system as a function of age; (8) τ o (τ) Instantaneous turnover time of carbon in the system as a function of age; the reciprocal of k(τ); (10) P(t) Net primary production (NPP) or gross primary production (GPP) as a function of time; (19) and (20)  A closer look at the age distributions of each pool (Fig. 6) suggests that some caution may in order when calculating radiocarbon ages from 14 C data. Although it is often assumed when calculating radiocarbon age in © 1999 Blackwell Science Ltd., Global Change Biology, 5, 371-394 soils (Balesdent 1987), several of the carbon pools do not have exponentially shaped storage response functions. For an exponential distribution to be appropriate the following criteria must be reasonably met: (i) carbon Table 2 Pool-by-pool analysis of the GPP-referenced carbon turnover model of Emanuel et al. (1981;Fig. 1a). Provided are the mean age T, mean transit time τ, turnover time τ o , and percentage of total carbon storage and respiration for each pool, using eqns 3, 6 and 7 and the solution to the GPP-referenced model of Emanuel et al. (1981)   entering the different pools comes directly from the atmosphere, (ii) transfer from one soil carbon pool to another is minuscule or nonexistent (i.e. reservoirs can only accept carbon directly from the atmosphere), and (iii) the pool must decay exponentially. In addition, the meaning of terms like turnover time, mean age and mean residence time should be properly distinguished; in many cases their values may differ significantly (ex: detritus/ decomposers; Table 2). τ o , T, and τ were calculated using (3), (6) and (7) with upper limits of integration that varied between 1 and 2000 (Fig. 7). For the GPP-referenced model, at 2000 years, τ o was 15, T was 70 and τ was 15 years. The differences between these numbers showed that overall the parent model did not decay exponentially, but more like the third case of Bolin & Rodhe (1973). At the end of the transformation the mean age of the NPP-referenced model was almost 10 years higher than that of the GPP-referenced model. In the NPP-referenced model, the turnover time was around 26 years but only 15 years in the GPPreferenced model. Slower turnover in the NPP-referenced model reflected the lack of plant respiration (and thus the lack of a fast initial return to the atmosphere).

Impulse response transformation of the CASA biosphere model
To illustrate this method on a spatially resolved, satellitebased NPP model, we determined the impulse response function of v1.2 of the CASA biosphere model (Figs 2 and  8,9,10). The advantage of using the globally aggregated impulse response functions Φ τ and Ψ τ determined in this ( ( way is that they are equivalent to the parent model and can be used like the simpler models of Joos et al. (1996), Emanuel et al. (1981) and Moore & Braswell (1994), but still capture spatial variability in decomposition and production. Similar functions can be produced for specific geographical regions, allowing integration into larger, more complex models of the entire global carbon cycle (Bolin et al. 1981;Enting et al. 1993;Enting et al. 1995). Examples of regions in which we might benefit having impulse response functions are the northern forests or the tropics. We suggest this method as a new standard for the 'simple' biosphere models used in global carbon cycle studies (Joos et al. 1996). The variation in τ o (τ) and k(τ) with cohort age for the entire globe (Fig. 9) demonstrates that carbon in the CASA model does not, taken as a whole, decompose exponentially (first-order decomposition would require constant τ o (τ) and k(τ) with τ). In the first two years, carbon turned over on a time-scale of 5 years; later on a time-scale of 35 years; and finally, as most of the carbon passed into the passive pool, the turnover time increased to hundreds of years. The plateauing observed here, where one pool or a group of similar pools at any given cohort age, dominates the turnover of the entire system, as has been previously observed (Schimel et al. 1994;McMurtrie & Comins 1996;Thompson et al. 1996). It has profound implications for transient studies, where the potential of a system to accrue carbon depends not only on the turnover time of the system at steady state, but also on the rate of change in production through time. Fast increases in production may lead to fast increases in respiration if turnover is rapid in the early stages of decomposition (Hungate et al. 1997).
The globally integrated turnover time, τ o , of the CASA model (eqns 7, 11 and 12), was 19.4 years ( Table 3). The global mean age of stored carbon T (eqn 6) and the global mean residence time of carbon τ (eqn 3) were 71.4 years and 19.2 years, respectively. The difference between the turnover time and the mean residence time was due the limited length of the transformation (2000 years), which was not sufficient for steady state. As seen with the simple model (Fig. 7) the longer the transformation is allowed to proceed, the more these numbers converge (Bolin & Rodhe 1973).
Differences in turnover time, mean age and mean transit time of the CASA model observed as a function of latitude (Fig. 10) indicates the fact that mean age and mean transit time have different controls. At high latitudes, the mean transit time is somewhat less than the turnover time by about 10% (since convergence of the two numbers requires more time where turnover times are high, such as in cool climates). Spatial variation in turnover and mean transit time appears to be driven by vegetation class (Defries & Townshend 1994;Randerson et al. 1996;Thompson et al. 1996; Table 3), which controls the presence or absence of wood, the turnover time of live biomass, and litter chemistry. Thus, mean residence and turnover times are dependent on the initial turnover characteristics of the system. The large swath of high turnover and mean residence times across Canada and northern Eurasia corresponds to the occurrence of vegetation class 4 (needleleaf evergreen forest, Defries & Townshend 1994). Mean age, on the other hand, is correlated with the long-term storage capacity of the soil, which is controlled primarily by climate and soil texture (Parton et al. 1993;Schimel et al. 1994).

Performance of an impulse response function against its parent model
To test impulse response functions as surrogates for their parent models, as well as their ability to operate under periods of changing turnover dynamics (a condition normally prohibited by the assumptions of the method), we repeated the analysis of Thompson et al. (1996) using the global impulse response function determined in the previous section (Fig. 8), assuming steady-state conditions at and before the beginning of the simulation, as well as an initial global annual NPP of 48.2 Pg C y -1 from Thompson et al. (1996), using (19) and (20).
We found close agreement (Ͻ 15% cumulative error) between the results from the impulse response function and the results from the Thompson et al. (1996) analysis (Fig. 11), but not enough agreement to justify the use of impulse response functions with these climate data sets for periods as long as 110 years. When performing analysis such as these, where it is known that climate or other factors may affect the turnover dynamics of the system, cumulative error should not be allowed to exceed 5-10%.
The change in NPP predicted by the impulse response method lacked the short-term variation seen in the work   Figs 8 and 9). Note the stronger latitudinal gradients in T relative to τ. Although it is hard to discern, at high latitudes τ undershoots τ o by as much as 10%.
of Thompson et al. (1996); high variation in the previous experiment was due to short-term variation in climate throughout the period, requiring high variation in NPP to match. It is uncertain that over timescales of decades this variation could be significant or even detectable given the large degree of uncertainty in the net terrestrial sink estimate. In addition, given the simplicity of this method, and the relative difficulty of the model protocol used by Thompson et al. (1996), limited analyses using impulse response functions over a period of a few decades are not only justifiable, but recommended. A good example of its utility would be in decadal deconvolution studies of atmospheric CO 2 and 13 CO 2 , where a link between terrestrial sequestration and changes in primary production could help constrain estimates of ocean and terrestrial biosphere sources and sinks. © 1999 Blackwell Science Ltd., Global Change Biology, 5, 371-394

Estimation of 13 C isodisequilibrium
Convolving an impulse response function, φ (τ), against past atmospheric 13 C content (Fig. 2) provides an efficient means of estimating the isodisequilibrium of 13 C, D 13 . We calculated the terrestrial D 13 and ζ 13 from 1977 to 1992 (Fig. 12) using three different impulse response functions: GPP-referenced from the Emanuel et al. (1981) model (Fig. 5); NPP-referenced from the Emanuel et al. (1981) model (Fig. 5); and NPP-referenced from the CASA model (Fig. 8). The GPP-referenced model from Emanuel et al. (1981) led to the highest average D 13 (21.9 Pg C ‰ y -1 ) for the period, and if incorporated into atmospheric 13 C budget analyses it would strongly shift the missing sink away from the terrestrial biosphere. The CASA model predicted the lowest average terrestrial Table 3 Carbon storage, net primary production (NPP), mean age T , mean residence time τ, and turnover time τ o of carbon after a 2000 years pulse response simulation, summarized for each of the vegetation classes defined in Defries & Townshend (1994). Values were calculated using the impulse response method by aggregating each vegetation class into their own impulse response functions (by area-weighted averaging, (11) Thompson et al. (1996) with varying climate, or using the protocol outlined in (19) and (20) and a globally aggregated impulse response function derived from the CASA model (Fig. 8). The cumulative increase in NPP predicted by the impulse response method from 1880 to 1990 undershoots the complete Thompson et al. (1996) model prediction by about 15%.
isodisequilibrium (16.6 Pg C ‰ y -1 ), and if incorporated would tend to shift the sink towards the terrestrial biosphere. Variation in D 13 for all three terms closely followed the negative rate of change in δ A (t) (Fig. 13), where D 13 was highest when δ A (t) declined most rapidly. Thus, the significant year-to-year variability in these estimates was mostly due to changes in δ A (t) at time t, rather than the long-term changes in δ A (t) of the years before. The highest average forcing coefficient, ζ 13 , belonged to the NPP-referenced Emanuel model (0.33‰), while the GPP-referenced Emanuel model predicted the lowest (0.19‰). Given the large variation observed in D 13 and ζ 13 , considerable care should be taken when determining what kind of model to use. Some confusion over this point has already arisen: Wittenberg & Esser (1997) assumed that the Emanuel model was NPP-referenced when they divided the isodisequilibrium (26.5 Pg C ‰ y -1 ; from Enting et al. 1993 a ) by 56 Pg C y -1 (the approximate NPP value for the Emanuel model) to derive a forcing coefficient of 0.47‰ for 1987. Consistency requires that they calculate a GPP-referenced forcing coefficient of about 0.23‰.
It should be emphasized that the high interannual a We were not able to reproduce the 13 C isodisequilibrium estimate of 26.5 Pg C ‰ y -1 for 1987 given by Enting et al. (1993). Using the piece-wise, linear δ 13 C record that they use in their study, we attempted to reproduce their results in three different ways: (i) using their method which gave an isodisequilibrium of 22.7 Pg C ‰ y -1 ; (ii) convolving a discrete version of Emanuel et al. (1981) GPP-referenced model against δ A (t) which gave a value of 23.0 Pg C ‰ y -1 ; and (iii) convolving an explicit integration of the Emanuel et al. (1981) impulse response function against the piece-wise δ A (t) of Enting et al. (1993) which gave a value of 23.1 Pg C ‰ y -1 . However, we did find that by using their method with the last section of the piece-wise, linear δ 13 C function changed so that it had the same slope as the previous term (the piece fom 1960 to 1980), we derived their published result (26.5 Pg C ‰ y -1 ).
variation in the isodisequilibrium, as calculated here, may have a significant impact on analyses using atmospheric 13 C and CO 2 data to calculate the terrestrial and ocean sinks. D 13 of the GPP-referenced Emanuel et al. (1981) model varied from as little as 17 Pg C‰ y -1 to as much as 28 Pg C‰ y -1 from 1977 to 1992. This variation is significant. Considering that the sink term of the budgetary equation (cf. Fung et al. 1997) is around 36 Pg C‰ y -1 (2 Pg C y -1 ϫ 18‰), ignoring this variation could lead to a significant underestimation of the variation in the distribution of land and ocean sinks. The efficiency of impulse response function in calculating D 13 and ζ 13 begs that instead of just assuming an interannually stable isodisequilibrium, it would be far better to include an impulse response function in the budgetary analysis. Indeed, if included, it would also be possible to solve for interannual variation in production using (19) and (20). Further effort should be made to use impulse response functions in these sorts of analyses. The global variation in ζ 13 (Fig. 13), from as little as 0.1‰ to more than 0.55‰ in the northern forests, highlights the strong spatial variability found in the terrestrial biosphere with respect to rates of τ o (Fig. 10). A globally uniform ζ 13 , as used by Enting et al. (1995), should be used with caution, especially since gross fluxes of carbon covary with turnover time . A high Fig. 13 The global distribution of ζ 13 in 1992 using (23) and each cell's φ (τ). The impulse response functions were convolved against δ A data from Fig. 3. ζ 13 superimposed on a large gross flux of carbon (high NPP) will lead to a large reduction in the time-decrease in atmospheric 13 C within that site's footprint and reduce the global need for a strong terrestrial photosynthesis to explain trends in atmospheric δ 13 C.

Estimation of stored and respired soil 14 C
If the carbon stored in the terrestrial biosphere is significantly older than the carbon respired (global T -71.4 years; global τ -19.2 years), then we would expect their corresponding 14 C ages to be quite different as well. However, a convolution of the northern hemisphere ∆ A data (Fig. 4) against Ψ τ and Φ τ from the CASA model (Fig. 8) indicates that these values are nearly the same: 170‰ and 175‰, respectively. The similarity is due, in this case, to the fact that ∆ A (t) does not change monotonically, but rather has a significant mode in 1964; the overlap of the two very different response functions, Ψ τ and Φ τ , across ∆ A (t) is such that (i) the peak is missed by the response function due to the higher probability of © 1999 Blackwell Science Ltd., Global Change Biology, 5, 371-394 younger carbon being represented (as in the case of Φ τ ), or (ii) it is smoothed out due to the high probability of much older carbon being represented (as in the case of Ψ τ ). The similarity between the two values is thus largely coincidental, but it does point out an important point, that despite large differences between τ and T, representative radiocarbon ages may not be divergent.
With northern and southern hemisphere ∆ A (t) data (Fig. 4), we calculated for 1995 the spatial distribution of bulk 14 C content in plants and soils, or the biospheric ∆ C , as well as the isodisequilibrium forcing coefficient of 14 C, ζ 14 (Fig. 14). The 'ghost' line along the equator of both images is due to the different ∆ A data used for the northern and southern hemisphere (Fig. 4). As far as we know, these are the first global estimates of spatially resolved ∆ C and ζ 14 .
With impulse response functions from three different cells in the CASA model, as well as separate functions for different pools within those cells, we made explicit comparisons between modelled and observed 14 C signatures. Three sites were chosen for this analysis: a tropical  (26) and historical ∆ A data (Fig. 4), as well as the 14 C forcing coefficient, ζ 14 , calculated by (27). wet forest site, a boreal evergreen needleleaf forest site, and a perennial grassland site. Several points (Table 4): (i) τ o varied by site, increasing in magnitude from the grassland site (9.83 years) to the boreal forest site (37.9 years), but T followed a different pattern, increasing from the tropical forest site (49.04 years) to the grassland site (141.5 years); (ii) the isodisequilibrium of 13 C increased with turnover time; and (iii) the turnover times of the different pools were typically less than their mean ages (the slow pools and microbial pools show this dramatically). Turnover time, in this case, is the time required to traverse the pool along, while mean age is the time required to traverse the entire system.
The last point demonstrates that an exponential decay assumption may not be valid when calculating the age of stored or respired carbon with radiocarbon. For example, soil inputs of above-and belowground woody biomass in the boreal site strongly skewed the age distribution of the soil pools away from exponential (extreme example: Table 4, boreal site slow pool, T vs. τ o ). Using an exponential decay assumption in this case would lead to an overestimation of the turnover time. There may be ways around this if properly modelled. Most promising is that with more detailed impulse response analyses of regional scale models, it may be possible to deconstruct currently Table 4 Carbon turnover statistics for three cells modelled under CASA, typical of a tropical forest (cell centred at 2.5°S, 59.5°W), a boreal forest (54.5°N, 66.5°W), and a perennial grassland (41.5°N, 99.5°W), as used in Thompson et al. (1996), from a 2000 years pulse response transformation. The model is NPP-referenced. The vegetation classification is by Defries & Townshend (1994) and soil texture classification is by Zobler (1986). 14 C and 13 C values are reported for 1996 using data from Figs 3 and 4. Turnover times for each pool are by (16 and  existing models such that they are more easily calibrated with 14 C data.

Uses and caveats
Impulse response functions derived from more complex carbon turnover models can be used (i) as surrogates for their parent models, (ii) as model diagnostics, or (iii) to determine the effect of the terrestrial biosphere on the global carbon cycle by convolving or deconvolving them against past changes in atmospheric or biospheric properties. It should be noted that these functions are only valid for linear systems with time-invariant dynamics; care must be taken when describing a potentially variable system with a single impulse response function, a problem pointed out by Joos et al. (1996) and Lewis & Nir (1978). A number of factors, including elevated nitrogen deposition and atmospheric CO 2 , may significantly influence terrestrial carbon dynamics and are worthy of further consideration in this context. Since impulse response functions are hysteric, any change in their shape may render their use over the historical period somewhat questionable. As such, one set of impulse response functions can go only so far to describe the terrestrial biosphere over a certain range of conditions.
The spatially aggregated impulse response functions derived from the CASA model are useful only when the initial geographical distribution of annual NPP (derived from NDVI) as well as terrestrial carbon turnover dynamics remain constant with time. On one level, this limits the usefulness of impulse response functions in comparisons of different models. But if the comparison is a steady state comparison, impulse response functions are ideal. These globally aggregated 1D models (Fig. 8) could become standards for comparison and validation in studies of present and future global carbon turnover. They may also become valuable tools in more complex analyses of the atmospheric budget of CO 2 , 13 CO 2 and 14 CO 2 .
Impulse response functions are useful comparative tools and may nicely complement other indices of model performance. The VEMAP (1995) exercise focused on changes in production and carbon storage under changing climate and doubled-ambient atmospheric CO 2 to determine the functional differences between three biogeochemistry models, TEM, Biome-BGC and CENTURY. Since impulse response functions are far more useful as prognostic variables than the mean transit time or mean age of a system, we urge that these functions be reported in any analysis of ecosystem response to perturbation. Impulse response functions are simple and require little CPU time. They have been used previously in descriptions of atmospheric, oceanic and terrestrial turnover (Bolin 1975), and as Joos et al. (1996) point out, although the interaction between different components of the global carbon cycle are sometimes nonlinear, their integrated function can be neatly described by a series of impulse response functions with explicit functions describing their interaction.

Conclusions
1 Impulse response functions facilitate diagnosis, comparison, and validation of terrestrial carbon biogeochemistry models. With them it is possible to calculate important flux-weighted parameters that are difficult to derive under complicated pool-orientated schemes. They can also be used as components of larger models describing atmosphere-ocean and atmosphere-biosphere exchange. Func-isodisequilibrium. The use of 14 C data with this method may prove a valuable tool for model validation.