Climate Dynamics �2001) 18: 189±202 Ó Springer-Verlag 2001

A nonlinear impulse response model of the coupled carbon cycle-climate system (NICCS) Abstract Impulse-response-function (IRF) models are designed for applications requiring a large number of climate change simulations, such as multi-scenario climate impact studies or cost-bene®t integrated-assessment studies. The models apply linear response theory to reproduce the characteristics of the climate response to external forcing computed with sophisticated state-of-the-art climate models like general circulation models of the physical ocean-atmosphere system and three-dimensional oceanic-plus-terrestrial carbon cycle models. Although highly computer ecient, IRF models are nonetheless capable of reproducing the full set of climate-change information generated by the complex models against which they are calibrated. While limited in principle to the linear response regime (less than about 3 C global-mean temperature change), the applicability of the IRF model presented has been extended into the nonlinear domain through explicit treatment of the climate system's dominant nonlinearities: CO 2 chemistry in ocean water, CO 2 fertilization of land biota, and sublinear radiative forcing. The resultant nonlinear impulse-response model of the coupled carbon cycle-climate system (NICCS) computes the temporal evolution of spatial patterns of climate change for four climate variables of particular relevance for climate impact studies: near-surface temperature, cloud cover, precipitation, and sea level. The space-time response characteristics of the model are derived from an EOF analysis of a transient 850-year greenhouse warming simulation with the Hamburg atmosphere-ocean general circulation model ECHAM3-LSG and a similar response experiment with the Hamburg carbon cycle model HAMOCC. The model is applied to two long-term CO 2 emission scenarios , demonstrating that the use of all currently estimated fossil fuel resources would carry the Earth's climate far beyond the range of climate change for which reliable quantitative predictions are possible today , and that even a freezing of emissions to present-day levels would cause a major global warming in the long term.


Introduction
For comprehensive integrated assessment and other climate impact studies, computations of climate change are often required for a large number of greenhouse gas (GHG) emission scenarios. The most reliable instruments currently available for the estimation of anthropogenic climate change are coupled atmosphereocean general circulation models (AOGCMs) in combination with three-dimensional models of the carbon cycle and other, non-CO 2 greenhouse gases. However, for multi-scenario investigations, these models are prohibitively expensive in computation time. Ideally, a climate model designed for application in integrated assessment and climate impact studies should provide the desired climate-change information without excessive computational cost, while nevertheless approaching the reliability and detail of sophisticated, top-of-the-line climate models.
While AOGCMs process a huge amount of information on the three-dimensional ocean-atmosphere system, only a small subset of the data is normally required as output to characterize the resulting climate change. One is interested typically only in some vector xt representing, for example, the change in a set of twodimensional ®elds such as near-surface temperature, cloud cover, precipitation or sea level. As input characterizing the external anthropogenic forcing f t one is similarly concerned only with low-dimensional ®elds, or even a scalar, like the globally integrated input of fossilfuel carbon dioxide into the atmosphere.
Provided the change relative to a reference climate state is small, the response of xt to an arbitrary (but suciently small) forcing f t is given generally by a convolution with the climate system's linear impulse response function (IRF) R: The function Rt represents the response to a d-function forcing at time t 0. Once the IRF has been determined, for example by ®tting to a single climate change simulation with a sophisticated climate model (or, in practice, to separate response experiments for the individual physical climate and greenhouse-gas modules from which the full climate model is constructed), the simple convolution (IRF) model can be applied to any time-dependent forcing scenario without further reference to the sophisticated climate model against which it was calibrated. As long as one remains within the linear regime, the IRF model then serves as an exact substitute for the full model. In principle, IRF models can be designed to reproduce, without loss of information, any output from a sophisticated model, including annual cycles and derived quantities like extreme value statistics. They provide a highly ecient method of computing credible timedependent climate-change scenarios. For a single input variable, the cpu times are of the order of a second on a workstation. For a multidimensional input with n f independent degrees of freedom (as would be required, e.g., to describe spatially variable aerosol emissions), the cpu time increases linearly with n f , and one requires n f 1 GCM reference experiments (including a control run) to calibrate the model. In our applications, however, we shall consider only CO 2 emissions as input. Since CO 2 is well mixed in the atmosphere on the time scales relevant for climate change, the input can be characterized in this case by a single scalar variable representing the global integral of the CO 2 emissions.
According to the linear-response-®tting exercises for oceanic CO 2 uptake by Maier-Reimer and Hasselmann (1987), the linear response range for the carbon cycle is constrained to CO 2 concentrations less than about twice the preindustrial value of pCO 2;p = 280 ppm, corresponding to an equilibrium warming of less than about 3 C. This is consistent with the linear-response limits found by Hasselmann et al. (1993) in their analysis of the cold-start errors of global warming simulations with AOGCMs.
The goals of this study are twofold: to extend this range of applicability by including the main limiting nonlinear physical processes into an IRF-based model, and to generalize an earlier IRF model for global mean temperature, used by Hasselmann et al. (1997) in a coupled climate-socioeconomic model for the cost-bene®t analysis of optimal CO 2 emission paths, to spa-tially dependent ®elds and other climate variables (cloudiness, precipitation, and sea level). The following nonlinearities are considered: 1. The solubility of additional CO 2 in ocean surface water decreases with rising concentrations. This reduces the uptake of the mixed surface layer and thereby the downward transport from the mixed layer into the large deep-ocean reservoir. 2. The net primary production of land vegetation, which is believed to act as a sink for anthropogenic carbon, is assumed to respond logarithmically to increasing atmospheric CO 2 (Bacastow and Keeling 1973;Enting et al. 1994). This has been incorporated previously in a terrestrial biosphere model by Joos et al. 1996. 3. The radiative greenhouse forcing increases only logarithmically with increasing CO 2 concentrations, as the infrared absorption is already close to saturation in the principal CO 2 absorption bands. Caldeira and Kasting (1993) have pointed out that the higher CO 2 concentrations resulting from the decrease in solubility tend to be compensated by the weaker logarithmic radiative forcing in the nonlinear system. Thus IRF models of the combined carbon-cycle and physical ocean-atmosphere system give a better linear approximation of the net response of the system than the IRF models of each of the subsystems separately. However, we ®nd that the cancellation of nonlinearities is only partial. Furthermore, since the climate policy debate often focuses on CO 2 concentrations rather than global warming scenarios, it is desirable to model each of the subsystems as accurately as possible.
In summary, the nonlinear impulse-response model of the coupled carbon cycle-climate system (NICCS) presented in the following is an extended version of the impulse response function (IRF) climate model used in the structural integrated assessment model (SIAM) by Hasselmann et al. (1997), augmented by nonlinear ocean carbon chemistry, a simple IRF representation of the terrestrial biosphere adapted from Joos et al. (1996), a logarithmic formulation of the radiative greenhouse forcing, and spatial patterns of change in four impactrelevant climate variables.
Comprehensive climate models used to compute the climate response to anthropogenic CO 2 emissions normally consist of two modules: a carbon cycle module to compute the atmospheric concentration of CO 2 for given CO 2 emissions, and a coupled atmosphere-ocean general circulation model (AOGCM) to compute the climate change resulting from the change in atmospheric CO 2 concentration. Our IRF model similarly consists of two IRF modules: a carbon-cycle (CarC) IRF module calibrated against a three-dimensional ocean carbon cycle model and augmented by a terrestrial biosphere model, and a physical-climate-change (CliC) IRF module calibrated against an AOGCM. The net NICCS (nonlinear impulse-response model of the coupled carbon cycleclimate system), comprising the CarC and CliC IRF modules, can be run in a coupled or sequential mode. In the experiments discussed later we have run the two IRF modules sequentially, as we found the temperature feedback to be relatively small (see also Maier-Reimer et al. 1996). For greater transparency in illustrating other more important features of the model, we have therefore preferred to neglect this eect.
A more complete representation of the climate feedback on the carbon cycle would need to include also the impact of a change in the ocean circulation on the physical carbon pump (the downwelling of CO 2 enriched surface waters in the North Atlantic and Antarctica into the deeper ocean) and the biological pump (the downward transport of CO 2 through the rain of decaying plankton), see Maier-Reimer et al. (1996) and Sarmiento et al. (1998). These feedbacks were found to partially cancel each other (Maier-Reimer and Hasselmann 1987) and were not included in the carbon cycle model against which our CarC IRF was calibrated. Also not activated in the computed response of the parent ocean carbon cycle model to anthropogenic emissions (although included inz the model) were marine biological processes, as the biological pump is limited by nutrients rather than CO 2 and is thus insensitive to anthropogenic CO 2 emissions. Other feedbacks which we have neglected, largely because of lack of reliable information, concern the impact of changes in temperature, water availability and other climatic factors on the terrestrial biosphere.
The study is organized as follows: the next section describes the carbon-cycle IRF module, consisting of the ocean and terrestrial components. The atmosphereocean climate IRF module is presented in Sect. 3, while applications of the coupled IRF model NICCS are discussed in Sect. 4. Section 5 summarizes the principal conclusions. Details of the ocean carbon cycle IRF are given in the Appendix.

The carbon cycle module
The carbon-cycle IRF module consists of two components: the ocean carbon cycle and a land vegetation module.

The ocean carbon cycle
A number of 3D ocean carbon cycle models have been developed to compute the oceanic uptake of CO 2 , for example the Hamburg Model (HAMOCC1, Maier-Reimer and Hasselmann 1987) or the Princeton Model Sarmiento and Sundquist 1992;Siegenthaler and Sarmiento 1993). For changes in the atmospheric CO 2 concentration less than a factor of about 2, most 3D ocean carbon cycle models can be characterized by their linear IRFsR c t (Maier-Reimer and Hasselmann 1987).
Linear IRF models of the oceanic carbon uptake have been developed and applied e.g. by Siegenthaler and Oeschger (1978), Oeschger and Heimann (1983) and Siegenthaler (1983). We base our nonlinear IRF model on a linear impulse response representation of HAMOCC by Maier-Reimer and Hasselmann (1987). Their linear IRF model has previously been used for estimating future atmospheric CO 2 e.g. by Harvey (1989) and for incorporation in a coupled climate-socioeconomic model (SIAM: structural integrated assessment model) by Hasselmann et al. (1997).
Since the advective and diusive transport within the ocean is essentially linear (unless the circulation is signi®cantly changed through feedback from the climate change), the accuracy of the linear approximation is limited only by the nonlinear uptake of CO 2 through the ocean surface, which is governed by the chemical dissociation equilibrium relating the CO 2 partial pressure pCO 2 to the concentration of dissolved inorganic carbon (DIC) in the nearsurface water. At higher concentrations, additional carbon becomes less soluble, and thus a smaller amount of surface-water carbon for a given increase of pCO 2 is available for mixing down into the deep ocean by thermohaline overturning (Maier-Reimer and Hasselmann 1987;Joos et al. 1996).
A successful attempt to circumvent the limitation of ocean carbon cycle IRF models to small perturbations was made by Joos et al. (1996). An IRF representation was used to describe the linear mixing and transport processes within the ocean, while the nonlinear air-sea exchange was modelled by a dierential equation. The explicit formulation of the gas exchange not only extended the range of applicability to greater concentrations, but enabled also the model to be applied to all conservative tracers with sources and sinks in the atmosphere, like bomb radiocarbon or even (for small temperature changes) the oceanic heat uptake. Although requiring only modest CPU resources, the IRF substitute model reproduced the response of spatially resolving models to within a few percent, both for a wide range of carbon emission scenarios and for the uptake of bomb radiocarbon. However, the computational eciency of the model was compromised by two factors: the need for two nested time-step loops (due to the dierential treatment of the nonlinear air-sea exchange and the separate integral treatment of transport and diusion), and, for very high anthropogenic CO 2emission scenarios (pCO 2 > 2000 ppm), by the inelasticity of the air-sea exchange, which required very short time steps for the differential mixed-layer computations.

The composite atmosphere-plus-mixed-layer system
In principle, the limitation to small time steps can be overcome by using an implicit integration method. However, for a nonlinear system, this requires time-consuming iterations. Alternatively, the problem can be circumvented by regarding the mixed-layer-plusatmosphere subsystem as equilibrated with respect to CO 2 exchange. This is permissible if the relevant time scales of climate change are long compared with the equilibration time of the mixed layer-plus-atmosphere subsystem (termed in the following simply the composite layer). The transport of CO 2 through the surface needs then no longer to be modelled by a dynamical equation, so that the shortest time-scale is suppressed and the model can be integrated with signi®cantly longer time steps.
The composite-layer IRF can be obtained by supression of the shortest time scale component of the atmospheric IRF of the parent model (see Maier-Reimer and Hasselmann 1987), which describes the atmosphere-mixed layer equilibration process, with subsequent renormalization of the reduced model. However, in its standard convolution-integral form the composite-layer IRF model is not suitable for the incorporation of the non-linear chemistry governing the oceanic CO 2 uptake. For this purpose, the model needs to be translated into an equivalent dierential representation that is physically interpretable in terms of the carbon capacities of the two subsystems of the composite-layer, the atmosphere and mixed layer. This can be achieved by constructing a box-model analogue of the IRF model in the form of a cascade of layers which are coupled through carbon¯uxes proportional to the dierences in the layer concentrations. Anthropogenic CO 2 emissions are introduced into the uppermost or zeroth layer, which represents the composite atmosphere-plus-mixed-layer system. The CO 2 input into the composite layer is distributed quasi-instantaneously between the atmosphere and the mixed layer, and the composite layer is then coupled to the rest of the ocean via its mixed-layer subsystem, which is in contact with the next-deeper layer.
The cascade's parameters (layer thicknesses and Newtonian relaxation coecients) are chosen such that the uppermost (com-posite) layer's IRF matches the composite-layer IRF derived from the parent 3D model's atmospheric response (see Appendix for the model equations and tuning conditions). The decomposition of the uppermost (composite) layer into its atmospheric and mixed-layer subsystems is chosen such that the ratio of carbon uptake into the sublayers is in accord with the preindustrial mixed-layer buer factor (Revelle and Suess 1957) for small perturbations. Once the linear cascade has been tuned in this way, the atmospheric and mixed-layer fractions for a larger change in the composite-layer carbon content c 0 are computed as nonlinear functions of c 0 from the nonlinear chemical equilibrium governing the relation between partial pressure and total inorganic carbon concentration in sea water, following Maier-Reimer and Hasselmann (1987).

Calibration
Our ocean carbon-cycle IRF is a recent least-squares ®t to the HAMOCC3i (inorganic) response to a sudden increase of the atmospheric CO 2 concentration by 1% (2.78 ppm). The model was run without a biological pump and without CaCO 3 sediment interaction. The asymptotic airborne fraction (13%) is close to the value found for the 1987 HAMOCC1 IRF (14%; the small dierence can be attributed to the onset of a nonlinear eect due to the dierent impulse sizes used for the two calibrations).
The nonlinear IRF model was checked, using values for the chemical equilibrium constants corresponding to the present-day global-mean temperature, against the full HAMOCC3i's response for impulses in which the preindustrial atmospheric CO 2 was increased by 1%, 25%, 100%, and 300% ( Fig. 1 and Table 1).
For the largest impulse, the CO 2 uptake of the nonlinear IRF module is a few percent slower than in HAMOCC3i: the nonlinear retardation of the carbon uptake is slightly overestimated. However, small errors in this range are to be expected, as the nonlinear air-sea exchange in the 3D model is spatially dependent and cannot be accurately simulated by a one-dimensional model using only a single set of global-mean chemistry and mixed-layer parameters.
The direct reduction in oceanic carbon uptake by sea surface warming through temperature-related chemistry changes only is consistently estimated by dierent global warming-marine carbon cycle studies (e.g. Sarmiento et al. 1998;Matear and Hirst 1999;Joos et al. 1999;Plattner et al. 2001). However the total reduction of oceanic carbon uptake through all climate feedbacks combined, (including modi®cations of the thermohaline overturning circulation) is hard to predict even with respect to the sign. There is almost no agreement between dierent climate models. Predictions of the THC range from a complete switch-o, especially of the Atlantic THC, up to a slight enhancement. The majority of models predict a moderate reduction. Clearly, such a reduction, if not compensated by Antarctic deep water, would cause a weakening of the downward transport of CO 2 ; it would, however, also reduce the upwelling of nutrient-and DIC-rich water to the surface. In high latitudes, the reduction of deep mixing would make the conditions for biological production more favourable (as the production takes place in the euphotic zone) and thus enhance the regional downward transport through the biological pump. Globally, changes in the biological cycle could lead to a transient increase or also decrease of the oceanic carbon uptake. Acting in opposite directions, the combined climate-induced changes in circulation, chemistry and biology were found in various studies to induce rather marginal modi®cations of the oceanic CO 2 uptake in the range between À7% and +7% (Maier-Reimer et al. 1996;Joos et al. 1999;Matear and Hirst 1999;Plattner et al. 2001). A stronger decrease (À35% in the uptake rate at 4 Â CO 2 , and the integral amounting to À10% at 2 Â CO 2 , or À20% at 4 Â CO 2 ) was simulated by Friedlingstein et al. (2001). Thus, our uncertainties in the oceanic CO 2 uptake are comparable to uncertainties in the present understanding of the surface-to-deep transport rates.

The terrestrial biosphere module
A CO 2 sink of roughly 2 GtC/year in the global terrestrial biosphere is believed to approximately compensate carbon losses from deforestation and other land use changes, mainly in the tropics. The allocation is ascribed to accelerated plant growth due to the rising CO 2 concentrations (CO 2 fertilization) and nitrate fertilization. The eciency of the terrestrial carbon sink and the question whether it will counteract fossil-fuel emissions in the future is hotly debated. Despite numerous papers on this topic in recent years (see e.g. the review of Schlesinger 1993;IPCC 1990;Tans et al. 1990;Keeling and Shertz 1992;Friedlingstein et al. 1995;Keeling et al. 1996;Sellers et al. 1996;Knorr 1997;Gayler and Claussen 1997;Joos and Bruno 1998;Claussen et al. 1999;Ganopolski et al. 1998 and many others), the issue is still far from resolved. Terrestrial biosphere models of dierent complexity and spatial resolution have been mapped onto CO 2 uptake impulse response function models or equivalent box-type analogues (Meyer et al. 1999;Thompson and Randerson 1999).
We augmented our ocean carbon cycle module by a simple four-box representation of the terrestrial biosphere (Siegenthaler and Oeschger 1987), to account to ®rst order for changes in the terrestrial carbon storage under rising CO 2 . The global terrestrial net primary production (NPP) is assumed to be proportional to the logarithm of the atmospheric CO 2 concentration, and the respiratory CO 2¯u x back into the atmosphere is linear in the four reservoir contents (Joos et al. 1996;Kicklighter et al. 1999). Like other current terrestrial carbon cycle models, our terrestrial biosphere model neglects the eects of land-use changes and other interference and corresponding losses of biological diversity and productivity. It neglects furthermore complete NPP saturation even at high CO 2 levels, and it neglects an accelerated respirative return of carbon to the atmosphere as is expected in a warmer climate. It is possible that these additional feedbacks would reduce the terrestrial CO 2 uptake.
The aggregate model of the terrestrial carbon cycle is tuned to match estimates of the terrestrial carbon sink during the 1980s (Schimel et al. 1997). When driven by emission scenarios, its response in terrestrial CO 2 uptake was well within the range obtained with current, spatially resolved models (Kicklighter et al. 1999).

Fig. 1
Nonlinear impulse-response of the nonlinear ocean-carbon cycle IRF model (NO, solid lines) compared with its parent 3D model HAMMOC3i (dot-dot-dashed) to impulses increasing the preindustrial atmospheric CO 2 content by 1% (lowest), 100% (medium), and 300% (highest curve), respectively. Only the perturbations are shown, normalized to the impulse size Table 1 Amplitudes and time constants for the oceanic CO 2 uptake IRFR c : computed from a least-squares ®t to the HA-MOCC3i response to a sudden 1% increase of its preindustrial atmospheric CO 2 (see Fig. 1)

The climate change module
In addition to the carbon-cycle (CarC) IRF module, we require as second component of our IRF model an IRF representation of the physical coupled atmosphere-ocean climate system. This was calibrated against the Hamburg AOGCM, as described in Hasselmann et al. (1993, 1997) Cubasch et al. (1992. However, in contrast to these applications, we consider now not only the global mean temperature as climate-change index, but generalize the climate change (CliC) IRF representation to include four representative two-dimensional ®elds: near-surface temperature (T 2m ), cloud cover (clo), precipitation (pre), and sea level (sea). Schlesinger et al. (1998) proposed a procedure for combining a number of ®xed spatial patterns with time-dependent coecients in which the spatial patterns were derived from several equilibrium climate change simulations with a coupled atmospheric generalcirculation/mixed-layer-ocean model, while the trajectories of the corresponding time-dependent coecients were computed by an energy-balance-climate/upwelling-diusion-ocean model. Although this provides a number of spatial patterns and dynamic responses which can be combined to describe the net responses to dierent forcing mechanisms (greenhouse gases and direct sulfate aerosol forcing, for example), the spatial signals and temporal evolutions were derived from dierent models and were therefore not necessarily consistent. Furthermore, the parent models were strongly simpli®ed in at least one of their components (atmosphere or ocean) and were thus less reliable with respect to the net spatiotemporal response than a fully coupled AOGCM. However, the validity of a separation of variables into spatial patterns of change with associated time-dependent factors has been con®rmed by other authors, see Huntingford and Cox (2000) and citations therein. Huntingford and Cox (2000) reproduced decade-averaged global mean changes in a number of impact-relevant surface climate variables (for each month of the year) using a two-box model ®t to two 150-year/250-year greenhouse integrations with the Hadley Center AOGCM. The corresponding relative scaling patterns of regional changes were obtained by ®tting the AOCGCM's (decade-averaged monthly) global mean time series to the individual time series in each land grid cell by variation of one scaling factor per grid cell. Thus the temporal and spatial signals represent one common parent model and calibration experiment.
In the following, we pursue the same basic pattern-projection strategy as Huntingford and Cox (2000), deriving both the spatial and temporal information simultaneously from the same transient AOGCM simulation. We then reproduce the (dimensionally reduced) space-time dependent AOGCM signal with an IRF model. To enhance the signal-to-noise ratio we separated the space-time dependent AOGCM output ®elds into noise and signal components using an EOF analysis applied to a long transient experiment (850 years) with the parent AOGCM, and restricted the analysis furthermore to annual means.

Regional climate change signals
To extract the climate change signal from the AOGCM response, consisting of a superposition of the externally forced signal and the natural variability of the AOGCM, we represent the response (forced scenario minus control run) as the superposition of a set of EOFs (empirical orthogonal functions) f v i x with associated time-dependent scalar coecients, the principal components (PCs) p v i t. The EOF decomposition maximizes the fraction of total variance explained at any given expansion order.
The time-evolution of the coecients p v i t for which the climate change signal can be clearly distinguished from the background natural climate variability (See Cubasch et al. 1992;Santer et al. 1994) can then be represented by an IRF model in the same way as the mean-temperature in the case of a singleindex CliC IRF model.
The EOF patterns and corresponding PCs used for our IRF model were extracted from an 850-year transient AOGCM simulation with the periodically synchronously coupled models EC-HAM3 and LSG (Voss et al. 1998;Voss and Mikolajewicz 2000). The (equivalent) CO 2 concentration was prescribed as exponential growth up to the fourfold`preindustrial' level (330 ppm) at year 120, after which the concentration was kept constant.
The analysis was carried out for the annual means of nearsurface temperature, precipitation, cloudiness, and sea level. The time series of the ®rst six principal components for each of the four ®elds are shown in the four left panels of Fig. 6. The PCs of second and higher order show statistical¯uctuations around zero, without a clear signal. However, the ®rst PCs, p v 1 t, of all four variables start close to zero, but then clearly emerge from the noise during the course of the simulation. The signal growth closely follows the increase in the forcing during the ®rst 120 years, but continues to increase after the forcing is kept constant, although on a slower time scale. The signal is best discernible in the sea level rise, but can be clearly distinguished from the noise also in the three atmospheric variables. Extrapolation of the p S 1 t curve for sea level rise suggests that equilibrium would not be reached until well after a thousand years.

Impulse-response representation
Since the climate change signals of all four variables considered can be well captured by the ®rst EOFs, while the higher EOFs are indistinguishable from the noise, the regional climate change signals can be reproduced by an IRF representation of just the ®rst term in each of the expansions Eq. (2), The assumed logarithmic relation between atmospheric carbon load w and radiative forcing corresponds to the standard representation of the near-saturation of the principal CO 2 infra-red absorption bands (IPCC 1990;Myhre et al. 1998). For concentrations below the reference value of 2 ÂCO 2 , the logarithmic expression yields a slightly warmer equilibrium than the linear model, while at higher concentrations the warming is signi®cantly weaker.
To ®t the IRF function R v t to the AOGCM scenario simulation, the function was represented as a sum of exponentials: where P i a i 1, so that S v represents the model's asymptotic climate sensitivity to a CO 2 doubling.
The climate change signal patterns f v 1 x were normalized to have unit global means (Figs. 2±5). Thus, the p v 1 t time series represent global mean climate change signals, the patterns indicating where the change is larger (f v 1 x > 1) or smaller (f v 1 x < 1) than the global mean.
The time constants s v i and amplitudes a v i of R v were obtained by a least-squares ®t of the IRF model to the p v 1 time series (Table 2). Low signal-to-noise ratio did not permit determination of more than two time constants of the IRFs. Least-square ®t experiments on a coarse 2d s v 1 -s v 2 grid (with the a v i optimized at each s v 1 -s v 2 gridpoint) indicated that for the atmospheric variables (temperature, cloud cover, precipitation) the ®t quality (rms) is relatively little changed within a range of appropriate combinations of the two time constants (300±700 years for s v 1 and 12± 28 years for s v 2 ). However, for sea level the ®rst time constant was well determined at s v 1 800 years, with s v 2 lying in the range 20± 30 years. This can be explained by the dierent relative weightings of the short and long time scales for the atmospheric variables, for which the short-time relaxation terms dominate over the long-time terms by factors of 2±4, compared with sea level, for which the short-term contribution is less than 4% of the long-term contribution.
Although cloud cover and precipitation are important variables for impact studies, their climate change signals exhibit lower signal-to-noise ratios than the near-surface temperature and are thus less reliably determined (Fig. 6, right panels). However, all three variables represent simultaneous expressions of the total atmospheric response and may therefore be expected to exhibit similar time response characteristics. The responses would be identical, for example, if the atmosphere responds quasi-instantaneously to changes in the sea surface temperature, sea-ice cover and land moisture distribution, and the dynamic response characteristics of these variables with``memory'' can be represented by a single joint EOF pattern. We have accordingly tested the ®t of a single IRF to all three variables, using the t 2m IRF model (which is most closely constrained by the data), appropriately scaled by the individual sensitivities (À0:87% for clo and 0.15 mm/d for pre). A good ®t was achieved for both cloud coverage and precipitation, with rms errors only 3% greater than those of the independent best ®ts for each variable.
For greater numerical eciency of the NICCS model, in which the CarC & CliC IRF modules were coupled together, we also constructed for the CliC IRF an equivalent dierential box-model analogue that could be directly coupled to the dierential equivalent of the CarC IRF module. Both models could then be integrated within the same time-integration loop, avoiding also the second nested time-variable loop required for the standard integral formulation of the CliC IRF. We point out in conclusion that the limitation to the onedimensional representation of annual mean values of four selected variables is not dictated a priori by the model design. The approach may be readily generalized to higher-order EOFs, if these can be reliably distinguished from noise through suciently long integrations or Monte Carlo simulations, and it can be applied to any variable that is (directly or indirectly) provided by the parent AOGCM, such as seasonal variability or higher-moment statistics. The system's response to forcing mechanisms other than greenhouse gases, like sulfate or volcanic aerosols, or solar variability, can be similarly treated in terms of further linearly superimposable IRF models.

Nonlinear impulse response of the coupled carbon cycle-climate system
To illustrate the main dynamical features of our nonlinear IRF model, we computed the response of three variants of the model to three dierent d-function CO 2 -emission inputs, representing a sudden increase of the pre-industrial atmospheric CO 2 concentra-   (Fig. 7). For comparison, we also ran the parent ocean carbon cycle model for these cases (Fig. 1). The IRF model variants were: 1. The linear ocean carbon uptake module combined with linear radiative forcing (the linear convolution, or LC variant, as used in the impulse-response climate module of Hasselmann et al. (1997, dot-dashed lines in Fig. 7). 2. The nonlinear ocean IRF analogue together with the logarithmic radiation model (the NO variant, solid lines); 3. The same as the NO variant, but with the nonlinear ocean carbon cycle augmented by a simple CO 2 -fertilized terrestrial biosphere carbon pool, adapted from Joos et al. (1996) (the BJ variant, dashed lines).
The weakest CO 2 input, representing an increase of the initial atmospheric CO 2 concentration by 25%, or % 140 GtC, relative to the pre-industrial level, corresponds to the total accumulated anthropogenic emissions from early industrialization until the 1980s. The response of the ocean carbon cycle (NO) is still close to the linear case (LC). Inclusion of the land biosphere pool (BJ) leads to a faster decay initially, which slows down later, however, when the additional sequestered biospheric carbon starts returning to the atmosphere. The small asymptotic biospheric retention is determined by the equilibrium between the slightly increased NPP and the respirative decay of the additional carbon.
The temperature responses of all three model variants exhibit a relatively rapid adjustment to the sudden CO 2 increase initially, with time scales governed by the heat uptake of the ocean, mainly in the upper 1 km. This is followed by a slow temperature decrease mirroring the decay of the CO 2 concentration. The nonlinear model yields substantially larger temperature changes than the linear model, as the CO 2 concentrations remain well below the 2Âpreindustrial level, the break-even point at which the linear and logarithmic greenhouse forcing are the same. The enhanced logarithmic forcing relative to the linear forcing in this low-concentration range overcompensates the concentration drawdown by the land biosphere carbon pool.
The response in sea level is dominated by the extremely slow warming of the deep ocean. Thus the fast initial temperature response to the sudden CO 2 increase, which was governed mainly by the heat uptake in the main thermocline of the ocean, does not appear as a signi®cant signal in the sea level response. This is characterized rather by the long time scales describing the gradual relaxation of the CO 2 concentration to its equilibrium asymptotic value, both processes being determined by the rate of penetration of tracers (CO 2 and heat) into the deep ocean.
The intermediate impulse, representing an initial doubling of the CO 2 concentration relative to the preindustrial state, corresponds to the estimated accumulated emissions (560 GtC) for a typical business-as-usual emissions scenario some time near the . . . ; p v 6 . Right: ®t of IRF models to the p v 1 time series. Sign reversals are due to renormalization, for convenience, to unit global pattern mean and global-mean time series. Also shown for cloud coverage and precipitation changes are the appropriately rescaled temperature response curves (dashed) middle of this century. Although the oceanic uptake is already somewhat slower than in the ®rst experiment with weaker input, the overall climate response is not drastically changed. The dierence between the linear and logarithmic greenhouse modules has become smaller, thereby reducing the greenhouse forcing relative to the linear case and partially compensating the eect of the relatively higher CO 2 concentrations resulting from the slower nonlinear carbon uptake. In contrast, the largest CO 2 impulse, corresponding to a sudden CO 2 quadrupling (1650 GtC input), is suciently large to drive the oceanic carbon uptake well out of its linear regime, although even this input is still substantially smaller than estimates of the total fossil fuel resources, including anticipated but not yet discovered resources, of 4000 to more than 25000 GtC, See IPCC (1996. As the peak temperature response of 3 C is now above the breakeven point of the logarithmic radiative forcing, the nonlinear variant yields a weaker forcing than the linear variant. The eect is suciently strong to over-compensate the higher CO 2 concentrations of the nonlinear ocean uptake model, so that the peak warming is slightly lower than in the linear variant. Both the peak warming and the subsequent decay of the CO 2 concentration are retarded relative to the linear case. The impact of the nonlinearities is least pronounced in all three model variants in sea level, where the largest impulse produces only a weak retardation relative to the linear case.

Long-term emission scenarios
Typical scenarios of climate change are computed over time horizons of 100 years (IPCC 1992(IPCC , 2000NakicenovicÂ et al. 1998). It has been pointed out by several authors, in particular Cline (1992) (see also Hasselmann et al. 1997) that this time span is too short to cover the full range of the climatic consequences of today's policies, leading to dangerous underestimates of long-term climate change impacts. For many of the scenarios currently under discussion, the emissions have not ceased growing by the end of the 21st century, and even after the emissions begin to fall, the cumulative CO 2 input continues to rise. Because of the long residence time of CO 2 in the atmosphere, it is the cumulative emissions rather than the instantaneous emissions that govern climate change. The slow uptake of the CO 2 input by the oceans and the terrestrial biosphere and the large heat capacity of the ocean together produce an exceedingly long memory of the climate system extending over many centuries (See Fig. 7). This is further illustrated in Figs. 8 and 9, which show the CO 2 concentrations and climate change computed with the three IRF model variants described above, together with the CO 2 concentrations computed with the ocean carbon cycle parent model, for two representative 1000-year emission scenarios. The ®rst case corresponds to a long-range`Business as Usual' (BAU) scenario in which essentially all estimated fossil fuel resources are burnt in the course of a few centuries. The scenario corresponds to typical BAU scenarios (see IPCC 1996) for the twenty ®rst century, while over the entire time horizon, the total cumulated BAU emissions amount to 15 000 GtC, which lies in the middle range of estimates of total fossil resources (see above). In the second`Frozen Emissions' (FRE) scenario, the emissions are kept constant at the 1990 level of 5.5 GtC/year. The FRE scenario is representative of the cumulative emissions of typical``drastic-reduction scenarios'' (see IPCC 1996) for the period up to 2100.
In both scenarios, the largest changes in atmospheric CO 2 and climate occur well after the year 2100, with a millenium-time-scale decay of the climate signal even after the emissions have faded out. For the BAU scenario, the CO 2 concentrations reach extremely high values, between ten and twenty times higher than the preindustrial level, for which direct physiological damages to living organisms must be expected. The associated temperature changes are of the order of 10 C. However, even the FRE scenario yields temperature changes of the order of 5 C in the long term, of the same order as the warming since the last ice-age. Climate changes of this magnitude lie, of course, well outside the linear regime, in a range in which all climate models, including the parent models against which NICCS was calibrated, are no longer reliable. Thus, the computations should be interpreted only as an indication and warning of the major, basically unpredictable climate changes that can be anticipated if business-as-usual or insuciently restrictive climate policies are pursued over long periods.
The linear ocean carbon cycle IRF model (LC) severely underestimates the CO 2 concentrations predicted by the parent model (HAMOCC3i) for the BAU scenario. The concentrations are reduced by 25% already before the year 2100, while the peak concentration is reached two centuries too early and is too small by a factor of three. In contrast, the atmospheric CO 2 concentration computed with the nonlinear ocean carbon IRF module (NO) agrees with the parent model to within 10% during the entire Fig. 7 Nonlinear response of the coupled carbon cycleclimate model (NICCS) to sudden increases of the preindustrial atmospheric CO 2 concentration by 25%, (left), 100% (center), and 300% (right). From top: atmospheric CO 2 perturbation, global-mean near-surface air temperature change ( C), and global-mean sea level change (m). Each panel shows the response of three IRF model variants: the nonlinear ocean CO 2 model without (NO, solid lines) and with land biosphere (NB, dashed), in both cases coupled to the logarithmic greenhouse forcing climate module, and the coupled linear convolution models of oceanic CO 2 uptake and climate change (LC, dot-dashed) 1000-year BAU period, including even the extreme peak value of 5000 ppm. Note that in all runs the initial state was de®ned as the preindustrial state in the year 1800, so that the dierent model variants yield dierent CO 2 concentrations and climate states already today. This is most visible in the run with the land biosphere module.
Driving the ocean-chemistry module at a strongly reduced temperature (Northern Atlantic winter instead of global mean temperature) did not modify the results signi®cantly, while inclusion of the land biosphere module (BJ) shifted the concentrations down by about 15% in both scenarios. However, the terrestrial biosphere could have a stronger impact if the feedbacks through changes in temperature and water availability are included. In the BAU scenario, the impact of the terrestrial carbon sink on climate is particularly weak, since the logarithmic radiative forcing is insensitive to relative changes in CO 2 concentrations at large background concentrations. In the FRE climate response, the combined nonlinearities of ocean chemistry, land vegetation, and radiation happen to very nearly cancel.
In general, the net climate response of the IRF module was found to be rather robust with respect to details of the carbon cycle module, for example with regard to the direct temperature eect on ocean chemistry, the capacity of the terrestrial pool or modi®cation of the terrestrial biosphere through changes in climate. Although the CO 2 concentrations become more uncertain at higher levels, this is compensated in part by the decreased sensitivity of the climate response to changes in the CO 2 concentration as the CO 2 infrared absorption bands become more saturated.
A comparison of the changes in global mean temperature, as proxy for the atmospheric variables, with sea level (Fig. 9) shows that sea level responds much more slowly than the atmospheric variables, as found already in the impulse experiments. Since the dominant time scale of the sea level response (800 years) is large compared even with the multi-century growth time of the BAU concentrations, the sea level response for this scenario is similar to the response to a step-function increase in CO 2 concentration discussed earlier. The sea-level rise for the FRE scenario is approximately linear over the entire period.
In further experiments, the optimal emission-path computations of Hasselmann et al. (1997) were repeated using NICCS rather than the linear IRF model SIAM. As expected, the cost-bene®t analyses of Hasselmann et al. (1997) were found to be robust with respect to the relatively small nonlinear modi®cations of the climate response in the range of modest climate change (normally less than 3 C warming) occurring in the optimal emission solutions.

Summary
Integrated assessment of anthropogenic climate change requires cost-ecient models of the carbon cycle and the atmosphere-ocean climate system that approach nevertheless the reliability and credibility of complex, stateof-the-art 3D carbon cycle and general circulation and for the next 100 years (lower panel). Computations were made with three IRF model variants: the nonlinear ocean CO 2 model without (NO, solid lines) and with land biosphere (NB, dashed) and with the linear convolution model (LC, dot-dashed), and also with the parent 3D ocean carbon cycle model HAMOCC3i (dot-dot-dashed) Fig. 9 Annual-global-mean climate change for scenarios BAU und FRE computed using three IRF model variants: the nonlinear ocean CO 2 model without (NO, solid lines) and with land biosphere (NB, dashed), in both cases coupled to the logarithmic greenhouse forcing climate module, and the coupled linear convolution models of oceanic CO 2 uptake and climate change (LC, dot-dashed) models. As a convenient tool for this purpose, we have developed nonlinear impulse-response-function (IRF) representations of the response characteristics of the HAmburg Model of the Ocean Carbon Cycle (HAM-MOC3i) and the Hamburg coupled atmosphere-ocean general circulation model ECHAM3-LSG. Coupled together, the net IRF model NICCS (nonlinear impulse response model of the coupled carbon cycle-climate system) computes the atmospheric CO 2 concentration and the resulting changes in selected impact-relevant climate ®elds (near-surface temperature, cloud cover, precipitation and sea level) for a prescribed 1000 year CO 2 emission scenario within less than a second on a workstation. NICCS is thus a valuable instrument for providing for the integrated assessment community the detailed output information of state-of-the-art climate models without loss of reliability for modest climate change at greatly reduced computational cost.
The limitation of IRF models to modest perturbations (below CO 2 doubling and 3 C warming) for which the climate response can be approximately linearized was partially overcome in NICCS by explicit treatment of two dominant nonlinearities: the nonlinear inorganic carbon chemistry governing the CO 2 uptake in the ocean, and the logarithmic dependence of the radiative greenhouse forcing on the CO 2 concentration. This was augmented by a land vegetation carbon cycle module with a nonlinear formulation of net primary production. Although inclusion of these nonlinearities removes the more obvious shortcomings of linear response models, it must be stressed that many more nonlinearities arise in the real climate system (and state-of-the-art climate models) at higher climate-change amplitudes, and these can not be adequately reproduced in an IRF model. The lowest-order extension of a linear IRF model to a general quadratic response model, for example, requires already the calibration of a set of three-index response coecients rather than the set of standard two-index response matrices of a linear IRF model. Thus the introduction of just two dominant nonlinearites into the NICCS model should be regarded only as a stop-gap measure to obtain more realistic order-of-magnitude estimates of amplitudes, without claims to a realistic description of the nonlinear modi®cations of the climate-change response patterns.
Another shortcoming of NICCS is that it neglects feedbacks of greenhouse warming on the ocean carbon cycle resulting from changes in the ocean circulation, since the ocean carbon cycle IRF was calibrated against a 3D carbon cycle model with a prescribed ocean circulation ®eld. This shortcoming can be overcome by calibrating both modules of a coupled IRF model against a coupled model of the carbon cycle and the general atmosphere-ocean circulation system. However, previous studies with 3D ocean carbon cycle models (Maier-Reimer et al. 1996;Sarmiento et al. 1998) indicate that the feedbacks of global warming on downwelling transport, vertical mixing, solubility and the biological pump partly compensate each other for atmospheric CO 2 concentrations up to about 700 ppm, leaving only a small residual eect of global warming on the oceanic CO 2 uptake. Maier-Reimer et al. (1996) conclude that`the currently used modelling strategy of ®rst using a carbon cycle model for the transformation of anthropogenic emissions into pCO 2 and subsequently using the output as forcing for a physical climate model, appears justi®ed'. We have accordingly run NICCS in the sequential, decoupled mode, without consideration of global warming feedbacks on the carbon cycle.
At atmospheric CO 2 concentrations exceding 700 ppm, 3D ocean carbon cycle models indicate that the reduced ocean circulation and other climate feedbacks tend to slow down the CO 2 ocean uptake, so that NICCS probably underestimates the atmospheric CO 2 concentration in the high-emission scenarios. Fortunately, however, the uncertainties of high CO 2 concentrations map into smaller climate change uncertainties through the logarithmic dependence of the radiative forcing on the CO 2 concentration.
Another limitation of present IRF models is that they are unable to simulate an unstable transition of the climate system to a new state, such as a breakdown of the ocean thermohaline circulation, a destabilization of the West Antarctic ice shield, a run-away greenhouse eect triggered by the release of methane trapped in permafrost regions, or a large-scale disruption of terrestrial ecosysems. The various nonlinear, physical-biogeochemical processes involved in surprises of this kind are not yet well understood. They cannot be reliably simulated or predicted today, even with the most sophisticated climate models. Thus, there exist at present no suitable parent models against which an appropriately extended nonlinear IRF model could be calibrated.
Conceptually, many of these shortcomings can probably be overcome by a suitable generalization of the basic NICCS structure, once the governing processes are understood and the relevant sophisticated parent models needed for calibration have been developed. However, an important generalization of the present NICCS which is feasible already today is the inclusion of further climate change variables provided by the parent model, such as annual and diurnal cycles, the occurrence of extreme events and, generally, changes in the statistics of the internal spatiotemporal variability of the climate system. It is in these properties that the impact of future climate change will probably be felt most strongly.
Our examples of the application of NICCS to longterm CO 2 emission scenarios demonstrated that the estimated total fossil-fuel resources are more than sucient to carry the climate system into a range of extreme CO 2 concentrations and temperature increases far in excess of the bounds within which any climate model can presently provide reliable predictions. Even a freezing of CO 2 emissions at 1990 levels is unable to stabilize the CO 2 concentration and limit global warming to acceptable levels in the long term. However, the long memory of the climate system provides also an opportunity for the gradual transition to carbon-free energy technologies over several decades, without dislocations of the global economy (see Hasselmann et al. 1997). A repeat of the optimal emission path computations of Hasselmann et al. (1997), in which the linear climate module of their coupled climate-socioeconomic model SIAM was replaced by NICCS, con®rmed the robustness of the conclusions of these authors with respect to model details. NICCS has been used and is currently being applied in integrated assessment studies Petschel-Held et al. 1999;FuÈ ssel and van Minnen 2000) in investigations of climate change feedbacks onto the terrestrial carbon cycle Meyer et al. 1999), and as an educational tool developed for the EXPO2000 World Exhibition. It is available as a community model on the Internet.

Appendix: The ocean carbon cycle module
A.1 The linear impulse response function of the composite-layer A good ®t to the linear atmospheric responseR c t of the 3D ocean carbon cycle model of Maier-Reimer and Hasselmann (1987) to atmospheric CO 2 input can be obtained by a sum of four decaying exponentials plus a constant de®ning the asymptotic equilibrium state: with s 0 I 7 and X i A i 1 : 8 The shortest decay time s 4 can be interpreted as the compositelayer equilibration time, the associated amplitude A 4 representing the fraction of a given d-impulse carbon input into the atmosphere at time t 0 which becomes dissolved in the ocean surface layer within the time scale s 4 . The ratio of the impulse-added CO 2 content of the composite layer to the change of the atmospheric CO 2 content shortly after equilibration is accordingly 1=1 À A 4 . The linearized impulse response of the composite layer can thus be obtained from the IRF representation of the atmospheric response of the complete model of Maier-Reimer and Hasselmann (1987) by dropping the short-time scale term and subsequently renormalizing:  (9), the transfer matrix is given by The Newtonian transfer coecients g i i 1; . . . ; n À 1 and the layer thicknesses h i i 0; . . . ; n À 1 are tunable constants. Unfortunately, the relations between the parameters h i ; g i of the dierential analogue and the the parameters a k ; s k of the composite-layer IRF R c are nonlinear and cannot be derived in closed form. They must be determined by satisfying a set of tuning conditions derived from the analytical Green function solution of the linear model equations.
A.3 Tuning the dierential analogue To tune the parameters of the dierential analogue model to the parameters of the linear IRF model, we ®rst diagonalize D by expressing both the solution c i t and the forcing e i t in terms of the eigenvectors of D: where C ik is the ith-layer component of the eigenvector C k of D associated with the eigenvalue k k (for the Newtonian relaxation system described by Eqs. (12, 13), the k k are real and positive): X j D ij C jk k k C ik : 16 Comparison of the zeroth-layer solution for the case of a d-impulse forcing (et c d dt; r k t r d k dt , with constant c d ; r d k ), with the composite-layer impulse response function (9) yields the tuning conditions k k s À1 k 17 r d k C 0k c d a k k 0; . . . ; n À 1 : 18 The condition k 0 s À1 0 0 is satis®ed through the conservation of carbon by the analogue model, which requires P j D ij 0 and therefore a singular propagator, jDj 0.
The eigenvectors C k and eigenvalues k k , and thereby also the forcing representation in eigenvector coordinates r d k , depend on the layer thicknesses h i i 0; . . . ; n À 1 and diusion coecients g i i 1; . . . ; n À 1. These must be determined numerically such that the conditions (17) and (18) are ful®lled.
With given s 0 I and the renormalization condition P i a i 1, the composite-layer IRF Eq. (9) has six remaining independent parameters. For n 4, the analogue clearly contains the required four time constants (one of which is in®nite). However, it contains seven rather than six free tuning parameters (four layer thicknesses h i and three diusion constants g i ). The additional degree of freedom arises because the analogue model computes only the carbon content of the layers, not their concentrations. Thus the model is determined through R c only up to an arbitrary scaling factor: the transport matrix D is homogeneous in the ratios g=h, and the thicknesses of all layers can be changed by an arbitrary factor, provided the carbon exchange coecients are changed by the same factor.
The additional degree of freedom can be ®xed by the known relation between the carbon content and CO 2 concentration of the atmosphere. Applying the known linear-limit ratio of the CO 2 concentrations in the two subsystems atmosphere, mixed layer of the composite layer, one obtains for the layer thicknesses of the mixed-layer and composite-layer (expressed in equivalent water units), after some algebra: h s n p w p A oc m C RC p Á A 4 1 À A 4 ; h 0 n p w p A oc m C RC p Á 1 1 À A 4 ; 19 where n is the preindustrial Revelle buer factor (see later), w p the preindustrial atmospheric carbon content, m C the molar mass of carbon, and A oc the area of the world ocean.
A.4 Introduction of nonlinear chemistry The dierential cascade analogue to the IRF model has been introduced primarily as a mathematical tool to reproduce the response of the more sophisticated carbon-cycle model to anthropogenic forcing via the tuning to an intermediate IRF model, without reference to real physical processes or observations. However, we ascribe now a speci®c physical interpretation to the composite layer of the dierential analogue by regarding it as composed explicitly of the atmosphere and the mixed layer. Although we have used this terminology already in the analysis, the physical interpretation has had no mathematical implications so far apart from determining the free scaling parameter of the analogue model (which is irrelevant as long as we consider only the net carbon content of the composite layer, without attaching a physical signi®cance to the layer). By interpreting the composite layer physically, we may now extend the linear IRF model into the nonlinear domain by considering the nonlinear chemical equilibration between the atmosphere and the mixed layer. While the Newtonian¯ux q 1 from the zeroth (composite) layer into the ®rst oceanic layer is still linear at large partial pressures of CO 2 , the surface-layer carbon content anomaly c s becomes a nonlinear function of the composite-layer carbon content anomaly c 0 . The function c s c s c 0 can be computed from the nonlinear chemistry of the mixed layer, given the thicknesses of the mixed-layer and the atmosphere and the equilibrium ratios of the partial pressures of CO 2 in the atmospheric and ocean. Thus, the expression c 0 =h 0 in the analogue's dynamical Eqs. (12, 13) is no longer equal to the surface layer concentration. However, the thicknesses of the composite layer and its subsystems are chosen to ensure this is the case in the linear limit: lim c030 c s h s c 0 h 0 : 20 The evolution Eq. (13) therefore need to be reformulated for the two uppermost layers in terms of the explicit nonlinear relation c s c 0 : The linear Eq. (12) and (13) for _ c 2 and _ c 3 remain unchanged. For the numerical integration of the nonlinear analogue, the nonlinear relation c s c 0 is evaluated at each time step by computing the chemical equilibrium determined by the various chemical processes associated with the dissolution and dissociation of CO 2 in seawater.
Solving the system for several slightly perturbed values of pCO 2 relative to its preindustrial value pCO 2;p yields a ®nite-difference estimate of the preindustrial buer factor n p . With consistent values of pCO 2;p , RC p , and n p , the thickness of the composite layer, h 0 , (and thereby its volume, given the total surface area A oc of the world ocean) can be computed from Eq. (19). The carbon content anomaly in the mixed layer, c s , can then be found numerically for any given carbon content anomaly in the composite layer, c 0 .