Hilbert problems for the geosciences in the 21st century

. The scientiﬁc problems posed by the Earth’s ﬂuid envelope, and its atmosphere, oceans, and the land surface that interacts with them are central to major socio-economic and political concerns as we move into the 21st century. It is natural, therefore, that a certain impatience should prevail in attempting to solve these problems. The point of this review paper is that one should proceed with all diligence, but not excessive haste: “festina lente,” as the Romans said two thousand years ago, i.e. “hurry in a measured way.” The paper traces the necessary progress through the solutions to the ten problems: A uniﬁed framework is proposed to deal with these problems in succession, from the shortest to the longest time scale, i.e. from weeks to centuries and millennia. The framework is that of dynamical systems theory, with an empha-sis on successive bifurcations and the ergodic theory of nonlinear systems. The main ideas and methods are outlined and the concept of a modeling hierarchy is introduced. The methodology is applied across the modeling hierarchy to Problem 5, which concerns the thermohaline circulation and its variability.


The problems of the Earth's fluid envelope and their societal relevance
Humanity has had for centuries a disruptive as well as a beneficial effect on its local environment: urban pollution and changes in flora and fauna over large fractions of the Earth's surface go back to the rise of major civilizations in Africa, the Americas, Asia and Europe several millennia ago. Increasing industrialization and the spread of industrial methods to agriculture, forestry and fisheries at the end of the 2nd millennium raise the possibility of the disruptive effects becoming global in the next centuries. In order to assess whether, to what extent, and in which ways we are modifying our global environment, it is essential to understand how this environment functions. We take, therefore, a planetary view of the Earth's climate system, of the pieces it contains, and of the way these pieces interact. This will allow us to understand how we might be acting on the individual pieces, and thus, on the whole. The global climate system is composed of a number of subsystems: atmosphere, biosphere, cryosphere, hydrosphere and lithosphere, each of which has distinct characteristic times, from days and weeks to centuries and millennia (Ghil and Childress, 1987;Trenberth, 1992). The atmosphere has a characteristic time of days-to-weeks in terms of the life cycle of extratropical weather systems. The global mixing of atmospheric trace gases, on the other hand, takes one or more years. The meanders and rings of the major wind-driven currents are the oceanic counterpart of weather systems; their characteristic time is as short as months to years. Temperature and salinity contrasts, on the other hand, drive the oceans' overturning circulation; its characteristic time is as long as centuries to millennia. Snow cover and sea ice have a huge seasonal cycle, as well as sub-and interannual variability, while continental ice sheets take many millennia to build up and at least centuries to collapse.
Each subsystem has therewith its own internal variability, all other things being constant, over a fairly broad range of time scales. These ranges overlap between one subsystem and another. The interactions between the subsystems thus give rise to climate variability on all time scales.
Can we hope to predict with confidence and eventually control in a rational way the effects of human intervention in this complex system? For humankind to survive and develop in a sustainable way in such a complex global environment as the climate system, we must at least understand the most basic workings of this system.
We outline here the rudiments of the way in which dynamical systems theory is starting to provide such an understanding. This understanding proceeds through the study of successively more complex flow patterns, for each time scale mentioned in the Abstract and in this section so far. The main features of dynamical systems theory (e.g. Smale, 1967;Arnol'd, 1983) that are important for the study of climate have been summarized by  and Ghil and Robertson (2000). They involve essentially bifurcation theory (e.g. Guckenheimer and Holmes, 1983) and the ergodic theory of dynamical systems (Eckmann and Ruelle, 1985).
Even a sketchy review of recent progress on all ten problems goes well beyond the scope of this paper. We thus give only an idea of the unified methodology that can help solve all ten. The methodology is illustrated by applying it to a subset of the problems. Problems 1-4 are merely touched upon, while Problems 5 and 9 are treated in somewhat greater detail. In Sect. 2, we describe the climate system's dominant balance between incoming solar radiation and outgoing terrestrial radiation. This balance is consistent with the existence of multiple equilibria of surface temperatures (Held and Suarez, 1974;Ghil, 1976;North et al., 1981). Such multiple equilibria are also present for other balances of climatic actions and reactions. Thus, the thermal driving of the mid-latitude westerly winds is countered by surface friction and mountain drag (Charney and DeVore, 1979;. These forces play an important role in solving Problem 1 and, hence, Problem 2. Multiple equilibria typically arise from saddle-node bifurcations of the governing equations (Ghil, 1994). Transitions from one equilibrium to another may result from small and random pushes, a typical case of minute causes having large effects in the long term.
In Sect. 3, we describe the ocean's overturning circulation between cold regions, where water is heavier and sinks, and warm regions where it is lighter and rises. Understanding this circulation, the forces that affect it and its resulting variability is the objective of Problem 5. The effect of temperature on the water's mass density and hence, motion, is in competition with the effect of salinity: density increases, through evaporation and brine formation, compete further with decreases in salinity and density through precipitation and river run-off. These competing effects can also give rise to two distinct equilibria (Stommel, 1961;Marotzke et al., 1988). In the present-day oceans, a thermohaline circulation prevails, in which the temperature effects dominate. About 50 million years ago, a halothermal circulation may have formed, with salinity effects dominating (Broecker et al., 1985;Kennett and Stott, 1991). In a simplified mathematical setting, these two equilibria arise by a pitchfork bifurcation that breaks the problem's mirror symmetry (Quon and Ghil, 1992;Thual and McWilliams, 1992).
On shorter time scales, of decades-to-millennia (Martinson et al., 1995), oscillations of intensity and spatial patterns in the thermohaline circulation seem to be the dominant mode of variability (Chen and Ghil, 1995). We show how interdecadal oscillations in the ocean's circulation arise by Hopf bifurcation (Quon and Ghil, 1995;Chen and Ghil, 1996).
In Sect. 4, we discuss the implications of multiple equilibria and interdecadal oscillations for our understanding of the effects that human activities might have on the climate system (Problem 9). The system's predictability in the absence of such effects (Problems 2 and 4) is presented (Lorenz, 1963a(Lorenz, , 1969Ghil et al., 1985Ghil and Jiang, 1998). Tentative conclusions are drawn about the identification and optimization of human effects on the climate system (Problem 10).

Connections to Hilbert's mathematical problems
Before launching into the ambitious enterprise described so far, we need to clarify the way in which these problems of atmospheric, oceanic and climate dynamics differ from, but also resemble, those posed by David Hilbert (1862-1943see Weyl, 1944, andReid, 1970) one hundred years ago for mathematics. First of all, in mathematics, one knows at least when a problem is solved, i.e. a conjecture is proved or disproved, once and for all. Errors may creep into a very lengthy and complicated mathematical proof, but the correctness of such a proof is, with notable exceptions, decidable at least in principle. The exceptions arise from the fundamental issues raised by Hilbert's (1900) first two problems that concern the completeness and compatibility of the arithmetic axioms (see also Gödel, 1940). This is not the case for problems in the physical sciences, in general, and those of the geosciences, in particular. Hilbert's sixth problem concerned, in fact, the mathematical treatment of the axioms of physics (see Corry, 1997), but progress in the sciences, unlike in mathematics, occurs to a large extent by inductive rather than deductive reasoning. In the sciences, a problem, such as Problems 3 and 4 listed in the Abstract, can only be answered based on observations available at a particular time. Given such observations, one solution among two or more that are being proposed may appear more satisfactory than another. As additional evidence in the form of observations or detailed computer simulations (Ghil and Robertson, 2000) becomes available, new solutions can be found and a better one can be selected (Kuhn, 1970).
The problems mentioned here are thus proposed as longterm challenges to geoscientists working on the Earth's fluid envelope. It is hoped that better solutions than those currently available will be found over the next decades. In particular, it is very important to have a coherent and self-consistent idea of a satisfactory set of solutions to Problems 1 through 9 in order to arrive at a well informed solution to Problem 10. The similarity with Hilbert's problems lies primarily in the two interrelated facts that (i) the problems cannot be solved in complete separation from each other; and (ii) ideas and methods developed in order to solve one problem can be immensely useful in solving others, as well.
2 Energy-balance models and the modeling hierarchy

Climate dynamics and the global environment
Climate dynamics is a modern scientific discipline that encompasses and extends beyond atmospheric and ocean dynamics. This view of the classical discipline of climatology first emerged about 40 years ago (Pfeffer, 1960). At this turn of the century, climate dynamics studies the variability of the atmosphere-ocean-cryosphere-biosphere-lithosphere system on time scales longer than the life span of individual weather systems and shorter than the age of our planet.
When defined in these broad terms, the variability of the climate system is characterized by a power spectrum that has three components. The first is a "warm-coloured" broadband component, with power increasing from high to low frequencies. The second is a line component associated with purely periodic forcing, both annual and diurnal. The third represents a number of broad peaks that might arise from external forcing that is not purely periodic (e.g. orbital changes or solar variability), from internal oscillations, or from a combination of the two (Mitchell, 1976;Ghil and Childress, 1987, Ch. 11;Ghil and Le Treut, 1999).
Understanding the climatic mechanism or mechanisms that give rise to a particular broad peak or set of peaks represents a fundamental goal of climate dynamics across all the time scales implicit in our ten problems. The regularities are of interest in and of themselves, for the order they create in our sparse and inaccurate observations. They also facilitate a prediction for time intervals comparable to the periods associated with a given regularity (Ghil and Childress, 1987, Sec. 12.6;Ghil and Jiang, 1998).
The climate system is highly complex, with its main subsystems having very different characteristic times, and the specific phenomena involved in each one of the climate problems listed in the Abstract are quite diverse. It is inconceivable, therefore, that a single model could successfully be used to incorporate all the subsystems, capture all the phenomena, and solve all the problems. Hence, the concept of a hierarchy of climate models, from the simple to the complex, has been developed about a quarter of a century ago (Schneider and Dickinson, 1974).
The methods of dynamical systems theory have been applied first to simple models, starting about forty years ago (Stommel, 1961;Lorenz, 1963a, b;Veronis, 1963Veronis, , 1966. More powerful computers now allow their application to fairly realistic and detailed models of the atmosphere (Strong et al., 1995;Keppenne et al., 2000), ocean (Jiang et al., 1995;Speich et al., 1995;Chen and Ghil, 1996;Dijkstra, 2000;Dijkstra and Ghil, 2001). We start, therefore, by presenting such a hierarchy of models.
This presentation is interwoven with that of the successive bifurcations that lead from simple to more complex solution behaviour for each climate model. Useful tools for comparing model behaviour across the hierarchy and with observations are provided by ergodic theory. Among these, advanced methods for the analysis and prediction of uni-and multivariate time series play an important role (Ghil and Vautard, 1991;Plaut et al., 1995;Ghil et al., 2001). Their applications to Problems 3 and 5 will also be mentioned here.

Radiation balance and energy-balance models (EBMs)
At present, the best developed hierarchy exists for atmospheric models; we summarize this hierarchy following Ghil and Robertson (2000). Atmospheric models were originally developed for weather simulation and prediction on the time scale of hours to days (Thompson, 1961). Currently, they serve in a stand-alone mode or coupled to oceanic and other models to solve all of the Problems 1-10 mentioned in the Abstract.
The first rung of the modeling hierarchy for the atmosphere is formed by zero-dimensional (0-D) models, where the number of dimensions, from zero to three, refers to the number of independent space variables used to describe the model domain, i.e. to physical-space dimensions. Such 0-D models essentially attempt to follow the evolution of global surface-air temperature T as a result of changes in global radiative balance (Crafoord and Källén, 1978;Ghil and Childress, 1987, Sect. 10.2): The other arrows show the hysteresis cycle that global temperatures would have to undergo for a transition from the upper stable branch to the lower one and back. The angle γ gives the measure of the present climate's sensitivity to changes in insolation (after Ghil and Childress, 1987).
Here R i and R o are incoming solar radiation and outgoing terrestrial radiation. The heat capacity c is that of the global atmosphere, plus that of the global ocean or some fraction thereof, depending on the time scale of interest: one might include only in c the ocean mixed layer when interested in subannual time scales, but the entire ocean could be included in c when studying paleoclimate. The rate of change of T with time t is given by dT /dt, while Q 0 is the solar radiation received at the top of the atmosphere, σ is the Stefan-Boltzmann constant, and µ is an insolation parameter, equal to unity for present-day conditions. To have a closed, selfconsistent model, the planetary reflectivity or albedo α and greyness factor m have to be expressed as functions of T ; m = 1 for a perfectly black body and 0 < M < 1 for a grey body, such as planet Earth.
There are two kinds of one-dimensional (1-D) atmospheric models, for which the single spatial variable is latitude or height, respectively. The former are so-called energy-balance models (EBMs: Budyko, 1969;Sellers, 1969) which consider the generalization of the model (2.1) for the evolution of surface-air temperature T = T (x, t), for example, Here the terms on the right-hand side can be functions of the meridional coordinate x (latitude, co-latitude, or sine of latitude), as well as of time t, and temperature T . The horizontal heat-flux term D expresses heat exchange between latitude belts; it typically contains first and second partial derivatives of T with respect to x. Hence, the rate of change of local temperature T with respect to time also becomes a partial derivative, ∂T /∂t.
The first striking results of theoretical climate dynamics were obtained in showing that slightly different forms of Eq. (2.2) could have two stable steady-state solutions, depending on the value of the insolation parameter µ (see Eq. (2.1b)) (Held and Suarez, 1974;Ghil, 1976;North et al., 1981). This multiplicity of stable steady states or physically possible "climates" of our planet can be explained in its simplest form in the 0-D model (2.1). This simple explanation resides in the fact that for a fairly broad range of µ-values around µ = 1.0 the curves for R i and R o as a function of T intersect in 3 points. One of these corresponds to the present climate (highest T -value), and another one corresponds to an ice-covered planet (lowest T -value); both of these are stable, while the third one (intermediate T -value) is unstable. To obtain this result, it suffices to assume that α = α(T ) is a piecewise-linear function of T , with high albedo at low temperature, due to the presence of snow and ice, and low albedo at high T , due to their absence, while m = m(T ) is a smooth, increasing function of T that attempts to capture in its simplest form the "greenhouse effect" of trace gases and water vapor (Ghil and Childress, 1987, Ch. 10).
The bifurcation diagram of such a 1-D EBM is shown in Fig. 1. It displays the model's mean temperature T as a function of the fractional change µ in the insolation Q at the top of the atmosphere. The 'S'-shaped curve in the figure arises from two back-to-back saddle-node bifurcations. The normal form of the first such bifurcation iṡ (2.3a) Here X stands for a suitably normalized form of T andẊ ≡ dX/dt is the rate of change of X, while µ is a parameter that measures the stress on the system, in particular, a normalized form of the insolation parameter. The upper-most branch corresponds to the steady-state solution X = + √ µ of Eq. (2.3a) and is stable. It matches rather well the Earth's present-day climate for µ = 1.0, more precisely, the steady-state solution T = T (x; µ) of the full 1-D EBM (not shown) matches closely the annual mean temperature profile from instrumental data over the last century (Ghil, 1976).
The intermediate branch starts out at the left as the second solution, X = − √ µ, of Eq. (2.3a) and is unstable. It blends smoothly into the upper branch of a coordinate-shifted and mirror-reflected version of (2.3a), for example, (2.3b) This branch, X = X o + √ µ 0 − µ, is also unstable. Finally, the lowermost branch in Fig. 1 is the second steady-state solution of Eq. (2.3b), X = X o − √ µ 0 − µ, and it is also stable. It corresponds to an ice-covered planet the same distance from the Sun as the Earth.
The fact that the upper-left bifurcation point in Fig. 1 is so close to present-day insolation values created great concern in the climate dynamics community in the mid-1970s, when these results were first obtained. Indeed, much more detailed computations (see Sect. 2.3 below) confirmed that a reduction of about 2-5% of insolation values would suffice to precipitate Earth into a "deep freeze" (Wetherald and Manabe, 1975). The great distance of the lower-right bifurcation point from present-day insolation values, on the other hand, suggests that one would have to nearly double atmospheric opacity for the the Earth's climate to jump back to more comfortable temperatures.

Other atmospheric processes and models
The 1-D atmospheric models in which the details of radiative equilibrium are investigated with respect to a height coordinate z (geometric height, pressure, etc.) are often called radiative-convective models (Manabe and Strickler, 1964;Ramanathan and Coakley, 1978;Charlock and Sellers, 1980). This name emphasizes the key role that convection plays in vertical heat transfer. While these models preceded historically EBMs as rungs on the modeling hierarchy, it was only recently shown that they too can exhibit multiple equilibria (Li et al., 1997;Rennó, 1997;Ide et al., 2001). The word (stable) "equilibrium," here and in the rest of this paper, refers simply to a (stable) steady-state of the model, rather than to true a thermodynamic equilibrium.
Two-dimensional (2-D) atmospheric models also have two types according to the third space coordinate which is not explicitly included. Models that resolve explicitly two horizontal coordinates on the sphere or on a plane tangent to it tend to emphasize the study of the dynamics of large-scale atmospheric motions (see Sec. II in Ghil and Robertson, 2000). They often have a single layer (Charney and DeVore, 1979;Legras and Ghil, 1985) or two (Lorenz, 1963b;Reinhold and Pierrehumbert, 1982). Those that resolve explicitly a meridional coordinate and height are essentially combinations of EBMs and radiative-convective models and emphasize therewith the thermodynamic state of the system, rather than its dynamics (Saltzman and Vernekar, 1972;MacCracken and Ghan, 1988;Gallée et al., 1991;Berger et al., 1998).
Yet another class of "horizontal" 2-D models is the extension of EBMs to resolve zonal, as well as meridional surface features, in particular, land-sea contrasts (Adem, 1970;North et al., 1983;Chen and Ghil, 1996). We shall see in Sect. 3.2 how such a 2-D EBM is used, when coupled to an oceanic model, to help solve Problem 5.
Additional types of 1-D and 2-D atmospheric models are discussed and references to these and to the types discussed above are given by Schneider and Dickinson (1974) and by Ghil and Robertson (2000), along with some of their main applications. Finally, to encompass and resolve the main atmospheric phenomena with respect to all three spatial coordinates, general circulation models (GCMs) occupy the pinnacle of the modeling hierarchy (Randall, 2000). Ghil and Robertson (2000) also describe the separate hierarchies that have grown over the last quarter of a century in modeling the ocean, and the coupled ocean-atmosphere system. More recently, an overarching hierarchy of Earthsystem models that encompass all the subsystems of interest (atmosphere, biosphere, cryosphere, hydrosphere and lithosphere) has been developing. Eventually, the partial results about each subsystem's variability, as outlined in this section and the next one, will have to be verified from one rung to the next of the Earth-system modeling hierarchy. No reliable solution to Problem 10 is conceivable, at least at the present time, without the full use of such a hierarchy.

Climate sensitivity
The results of climate simulations with GCMs, whether atmospheric or coupled, are often interpreted in terms of the understanding gained from 0-D or 1-D EBMs. Wetherald and Manabe (1975), using a simplified sectorial GCM, confirmed the dependence of mean zonal temperature on the insolation parameter µ (the normalized "solar constant") obtained for 1-D EBMs (see Ghil and Childress, 1987, Ch. 10). These authors could not confirm the presence of the intermediate, unstable branch or of the "deep-freeze" stable branch in Fig. 1 with their GCM, due to the model's computational limitations. But the parabolic shape of the upper, presentdaylike branch near the upper-left bifurcation point in our figure (Eq. 2.3a) was well supported by their GCM simulations.
In fact, the sensitivity tan γ = (dT /dµ)| µ=1.0 of global temperature T to changes in µ near the present-day climate (see Fig. 1) equals about 1K per 1% change in the insolation for both EBMs and GCMs. Many studies of climate-change response to increases in greenhouse trace-gas concentrations use, therefore, a linearized version of Eq. (2.1), to interpret the results of detailed GCM simulations. Such interpretations focus on the roles of the different feedbacks λ i , positive (λ i < 0) or negative (λ i > 0), and heat sources, Q j > 0, or sinks, Q j < 0 (e.g. Schlesinger and Mitchell, 1987;Cess, et al., 1989;Li and Le Treut, 1992). These interpretations are limited, therefore, to equilibrium or quasi-equilibrium responses. The complete study of the human effects on changes in climatic variability requires the use of models that reproduce natural variations in a satisfactory manner, rather than simulating only equilibria and simple transients (Ghil and Le Treut, 1999). We proceed, there- Table 1. Thermohaline circulation (THC) oscillations (adapted from Ghil, 1994)
fore, with a description of internal variability that arises in the context of solving Problem 5.

Theory and simple models
Historically, the thermohaline circulation (THC) was first among the climate system's major processes to be studied using a very simple mathematical model. Stommel (1961) formulated a two-box model and showed that it possessed multiple equilibria. A sketch of the Atlantic Ocean's THC and its interactions with the atmosphere and cryosphere on long time scales is shown in Fig. 2. These interactions can lead to climate oscillations with multi-millennial periods, such as the Heinrich events (Ghil, 1994, and references therein), and are summarized in the figure's caption, following . An equally schematic view of the global THC is provided by the widely known "conveyor belt" diagram (e.g. Broecker, 1991). The latter diagram does not commonly include the THC's interactions with water in its gaseous and solid phases, which our figure does include.
Basically, the THC is due to denser water sinking, lighter water rising, and water-mass continuity closing the circuit through near-horizontal flow between the areas of rising and sinking. This is roughly the oceanic equivalent of the atmosphere's Hadley circulation, with two notable differences: i) The ocean water's density ρ is a function of temperature T and salinity S, while that of the air depends on temperature and humidity.
ii) Water sinks in and near fairly concentrated regions of intense convection, currently located primarily in high latitudes, and rises diffusely over the rest of the ocean. Air, on the other hand, does rise most intensely in cumulus towers, but overall, the areas of net rising and sinking air in a Hadley cell are quite comparable in extent when viewed on the synoptic and planetary scales.
The effects of temperature and salinity on the ocean water's density, ρ = ρ(T , S), oppose each other: the density ρ decreases with increasing T and increases with increasing S. It is these two effects that give the thermohaline circulation its name, from the Greek words for T and S. In high latitudes, ρ increases as the water loses heat to the air above and, if sea ice is formed, as the water underneath is enriched in brine. In low latitudes, ρ increases due to evaporation but decreases due to heat flux in the ocean.
For the present climate, the temperature effect is commonly assumed to be stronger than the salinity effect. Ocean water is observed to sink in certain areas of the high-latitude North Atlantic and Southern Oceans, with very few, limited areas of deep-water formation elsewhere, and to rise everywhere else. Thus, the adjective "thermohaline" indicates that T is more important than S by placing "thermo" before "haline". During some remote geological times, deep water may have formed in the global ocean near the equator; such an overturning circulation of the opposite sign to that prevailing today has been dubbed halothermal, S before T (e.g. Kennett and Stott, 1991;see, however, Bice and Marotzke, 2001). The quantification of the relative effects of T and S on the oceanic water masses' buoyancy in high and low latitudes is far from complete, especially for paleocirculations; the association of the latter with salinity effects that exceed the thermal ones is thus rather tentative. Stommel (1961) considered a two-box model, with two pipes connecting the two boxes. He showed that the system of two nonlinear, coupled ordinary differential equations (ODEs), which govern the temperature and salinity differences between the two well-mixed boxes has two stable steady-state solutions, distinguished by the direction of flow in the upper and lower pipe. Stommel's paper was primarily concerned with distinct local convection regimes, and hence vertical stratifications, for example, in the North Atlantic and Mediterranean (or Red Sea). Today, we primarily think of one box as representing the low latitudes and the other one as representing the high latitudes in the global THC. The next step in the hierarchical modeling of the THC is that of 2-D meridional plane models (see Sect. I.B of Ghil and Robertson, 2000), in which the temperature and salinity fields are governed by coupled nonlinear partial differential equations with two independent space variables, for example, latitude and depth. Given boundary conditions for Changes in the radiation balance R in − R out are due, at least in part, to changes in the extent of the northern hemisphere (NH) snow and ice cover, V , and how they affect the global temperature, T ; the extent of the southern hemisphere ice is assumed constant, to a first approximation. The change in hydrologic cycle expressed in the terms P rain − P evap for the ocean and P snow − P abl for the snow and ice is due to changes in ocean temperature. Deep-water formation in the North Atlantic Subpolar Sea (North Atlantic Deep Water: NADW) is affected by changes in ice volume and extent, and regulates the intensity C of the THC; changes in Antarctic Bottom Water (AABW) formation are neglected in this approximation. This, in turn, affects the system's temperature, and is also affected by it (after ). such a model that are symmetric about the equator, as are the equations themselves, one expects a symmetric solution in which water either sinks near the poles and rises everywhere else (thermohaline), or sinks near the equator and rises everywhere else (halothermal). These two symmetric solutions would correspond to the two equilibria of Stommel's (1961) box model.
In fact, Fig. 3 shows that symmetry breaking can occur, leading gradually from a symmetric two-cell circulation (Fig. 3a) to an antisymmetric one-cell circulation (approximately achieved in Fig. 3c). In between, all degrees of dominance of one cell over the other are possible, with one such intermediate state shown in Fig. 3b. A situation lying somewhere between Figs. 3b and c seems to resemble most closely the meridional overturning diagram of the Atlantic Ocean in Fig. 2.
This symmetry breaking can be described by a pitchfork bifurcation (e.g. Guckenheimer and Holmes, 1983): (3.1) Here X stands for the amount of asymmetry in the solution, so that X = 0 is the symmetric branch, and µ is a parameter that measures the stress on the system, in particular, a normalized form of the buoyancy flux at the surface. For µ < 0 the symmetric branch is stable, while for µ > 0, the two branches X = ± √ µ inherit its stability.
Thus, Figs. 3b and c both lie on a solution branch of the 2-D THC problem for which the left cell dominates, i.e. North Atlantic Deep Water extends to the Southern Ocean's polar front, as it does in Fig. 2. According to Eq. (3.1), another branch exists whose flow patterns are mirror images in the rectangular box's vertical symmetry axis (the "equatorial plane") of those in Figs. 3b and c. The existence of this second branch was verified numerically by Quon and Ghil (1992;their Fig. 16). Thual and McWilliams (1992) considered more complex bifurcation diagrams for a similar 2-D model and showed the equivalence of such a diagram for their 2-D model and a box-and-pipe model of sufficient complexity (see also Dijkstra and Ghil, 2001). Bryan (1986) was the first to document the transition from a two-cell to a one-cell circulation in a simplified ocean GCM with idealized, symmetric forcing, in agreement with the three-box scenario of Rooth (1982). Manabe and Stouffer (1999), among others, questioned, however, the realism of more than one stable THC equilibrium by using coupled ocean-atmosphere GCMs. The situation with respect to the THC's pitchfork bifurcation (3.1) is thus subtler than it was with respect to Fig. 1 for radiative equilibrium. In Sect. 2, atmospheric GCMs confirmed essentially the EBM results. Climbing the rungs of the modeling hierarchy for the THC, however, yields contradictory results that are still in need of further clarification.

Bifurcation diagrams for GCMs
Internal variability of the THC, with smaller and more reg- (b, c) increasingly asymmetric solutions as the forcing is increased (from Quon and Ghil, 1992). ular excursions than the huge and completely irregular jumps associated with bistability, was studied intensively in the late 1980s and the 1990s. These studies placed themselves on various rungs of the modeling hierarchy, from Boolean delay equation models (so-called "formal conceptual models": Darby and Mysak, 1993) through box models (Welander, 1986) and 2-D models (Quon and Ghil, 1995) to ocean GCMs. A summary of the different kinds of oscillatory variability found in the latter appears in Table 1. Additional GCM references for these three types of oscillations are given by McWilliams (1996). Such oscillatory behaviour seems to match more closely the instrumentally recorded THC variability (see Sect. 3.3 below), as well as the paleoclimatic records for the recent geological past (see Ghil, 1994), than bistability. The interaction of the (multi-)millennial oscillations with variability in the surface features and processes shown in Fig. 2 is discussed by Ghil (1994). Chen and Ghil (1996), in particular, studied some of the interactions between atmospheric processes and the THC. They used a so-called hybrid coupled model to depict a (horizontally) 2-D EBM (see Sect. 2.3) coupled to a rectangular-box version of the North Atlantic rendered by a low-resolution ocean GCM. This hybrid model's regime diagram is shown in Fig. 4a. A steady state is stable for high values of the coupling parameter λ ao or of the EBM's diffusion parameter d. Interdecadal oscillations with a period of 40-50 years are self-sustained and stable for low values of these parameters.
The self-sustained THC oscillations in question are characterized by a pair of vortices of opposite sign that grow and decay in quadrature with each other in the ocean's upper layers. Their centers follow each other counterclockwise through the northwestern quadrant of the model's rectangular domain. Both the period and the spatio-temporal characteristics of the oscillation are thus rather similar to those seen in a fully coupled GCM with realistic geometry (Delworth et al., 1993).
The transition from a stable equilibrium to a stable limit cycle, via Hopf bifurcation, in Chen and Ghil's hybrid coupled model is shown in Fig. 4b. The physical characteristics of the oscillatory instability that leads to the Hopf bifurcations have been described in further detail by Colin de Verdière and Huck (1999), using both a four-box oceanatmosphere and a number of more detailed models.

Interannual and interdecadal climate variability
Human intervention in the workings of the climate on the global scale (Problem 9) is most likely to occur on time scales comparable to those of major socio-economic changes. The latter typically take place at this turn of the century on interannual to interdecadal time scales. The solution to Problem 9 is thus clearly predicated, at least, on satisfactory solutions to Problems 3 and 5.
The best known climatic regularities on the interannual time scale are the quasi-biennial and low-frequency (4-5-year) component of the El-Niño/Southern-Oscillation (ENSO) phenomenon (Rasmusson et al., 1990;Neelin et al., 1998). A particularly appealing explanation of these two spectral peaks that also allows for the broadband spectrum in the 1-10-year band is given by the Devil's staircase mechanism (Jin et al., 1994(Jin et al., , 1996Tziperman et al., 1994;Saunders and Ghil, 2001). This broadband spectral background is consistent with the occurrence of major warm events in the tropical Pacific every 2-7 years.
This mechanism is most easily explained in two steps. First, it is by now well accepted that the coupled oceanatmosphere system in the tropical Pacific exhibits a selfsustained oscillation with a period of 2-3 years for meanannual conditions. Second, this self-sustained oscillation in-teracts with the seasonal cycle in insolation and it frequencylocks, preferentially to an integer multiple of the forcing period: 2, 3, 4 or 5 years. Jumps between these broad integer steps in the Devil's staircase, as well as between narrower steps associated with rational periods, such as 4/3 years (≈ 16 months), occur as the result of atmospheric, higher-frequency "noise." This mechanism is described further and compared with other scenarios of interannual variability and with the observational evidence in Sec. III of Ghil and Robertson (2000); see also Ch. 5 of Dijkstra (2000).
A fairly satisfactory, albeit partial, solution to Problem 3 is thus at hand, insofar as we understand ENSO's spectral regularities. No complete consensus has been achieved, however, as yet, on the relative role of weather noise and coupledmode nonlinearities in ENSO's irregularity, i.e. the irregular occurrence of major warm events. Ghil and Vautard (1991) described statistically significant regularities in the interdecadal band of climatic variability. Broad peaks of roughly 13-15 and 23-27 years have been since confirmed in various records, both instrumental (Plaut et al., 1995;Moron et al., 1998) and paleoclimatic. The evidence for these interdecadal peaks is reviewed by Ghil et al. (2001). The interdecadal oscillations in the THC reviewed in Sect. 3.2 above are a plausible mechanism for these peaks, but not yet generally accepted as such. We are thus making good progress on Problem 5, but do not expect a fully selfconsistent solution in less than 2-3 decades.

Concluding remarks
We have cast a bird's-eye view on selected problems among the ten formulated in the Abstract. This view has clearly emphasized strong interconnections between them, both pairwise and all around. Unified methods of solution are provided by dynamical systems theory. Each problem requires the use of a model hierarchy for its solution.
To conclude, we briefly address Problem 9. More precisely, we ask whether the impact of human activities on the climate is observable and identifiable in the instrumental records of the last century-and-a-half and in recent paleoclimate records? The answer to this question depends on the null hypothesis against which such an impact is tested.
The current approach that is generally pursued assumes essentially that past climate variability is indistinguishable from a stochastic red-noise process (Hasselmann, 1976), whose only regularities are those of periodic external forcing (Mitchell, 1976). Given such a null hypothesis, the official consensus of IPCC (1995) tilts towards a global warming effect of recent trace-gas emissions, which exceeds the cooling effect of anthropogenic aerosol emissions.
Atmospheric and coupled GCM simulations of the tracegas warming and aerosol cooling buttress this IPCC consensus. The GCM simulations used so far do not, however, exhibit the observed interdecadal regularities described at the end of Sect. 3.3. They might, therewith, miss some important physical mechanisms of climate variability and are, there-a) Regime diagram b) Bifurcation diagram Fig. 4. Dependence of THC solutions on two parameters in a hybrid coupled model (HCM); the two parameters are the atmosphereocean coupling coefficient λ ao and the atmospheric thermal diffusion coefficient d. (a) Schematic regime diagram. The full circles stand for the model's stable steady-states, the open circles for stable limit cycles, and the solid curve is the estimated neutral stability curve between the former and the latter. (b) Hopf bifurcation curve at fixed d = 1.0 and varying λ ao ; this curve was obtained by fitting a parabola to the model's numerical-simulation results, shown as full and open circles (from Chen and Ghil, 1996). fore, not entirely conclusive.
As northern hemisphere temperatures were falling in the 1960s and early 1970s, the aerosol effect was the one that caused the greatest concern. As shown in Sect. 2.2, this concern was bolstered by the possibility of a huge, highly nonlinear temperature drop if the climate system reached the upper-left bifurcation point of Fig. 1.
The global temperature increase through the 1990s is certainly rather unusual in terms of the instrumental record of the last 150 years or so. It does not correspond, however, to a rapidly accelerating increase in greenhouse-gas emissions or a substantial drop in aerosol emissions. How statistically significant is, therefore, this temperature rise, if the null hypothesis is not a random coincidence of small, stochastic excursions of global temperatures with all, or nearly all, the same sign?
The presence of internally arising regularities in the climate system with periods of years and decades suggests the need for a different null hypothesis. Essentially, one needs to show that the behaviour of the climatic signal is distinct from that generated by natural climate variability in the past, when human effects were negligible, at least on the global scale. As discussed in Sects. 2.1 and 3.3, this natural variability includes interannual and interdecadal cycles, as well as the broadband component. These cycles are far from being purely periodic. Still, they include much more persistent excursions of one sign, whether positive or negative in global or hemispheric temperatures, say, than does red noise. Ghil and Jiang (1998) showed that on the seasonal-tointerannual time scale, climate predictability is greatly enhanced by the presence of these regularities. This is the case not only when the predictability is measured against the classical benchmark of a red-noise process of first-order autoregressive type (Hasselmann, 1976;Mitchell, 1976;IPCC, 1995). It is also the case when measured against a purely chaotic, albeit deterministic process, such as that governed by Lorenz's (1963a) model.
The answer to the global warming problem on the interdecadal time scale on which we are posing this problem is comprised of two parts. First, we need to describe and understand climate variability on this time scale as well as on the interannual time scale, in particular its regularities, i.e. solve Problems 3 and 5. Once that is done, a prediction with known error bars on the interdecadal time scale will be possible (see Plaut et al., 1995).
With these results in hand, we should be able to proceed to the second part of the answer. Can we identify with measurable certainty deviations of the current record from predictions based on past natural variability? If so, such deviations have to be attributed to new causes. The "suspects" clearly include human effects, and attribution to them will become thereby both easier and more reliable.
At the same time, one hopes that the applications of dynamical systems theory to the global socio-economic system and to population dynamics will have made considerable progress (Keilis-Borok and Sánchez Sorondo, 2000). With this in mind, a rational approach to predicting and controlling the coupled system formed by all living beings on this planet and their physico-chemical environment would become feasible by the end of this century.