Conceptual climate models are very simple mathematical representations of climate processes, which are especially useful because their workings can be readily understood. The usual procedure of representing effects of unresolved processes in such models using functions of the prognostic variables (parametrizations) that include no randomness generally results in these models exhibiting substantially less variability than do the phenomena they are intended to simulate. A viable yet still simple alternative is to replace the conventional deterministic parametrizations with stochastic parametrizations, which can be justified theoretically through the central limit theorem. The result is that the model equations are stochastic differential equations. In addition to greatly increasing the magnitude of variability exhibited by these models, and their qualitative fidelity to the corresponding real climate system, representation of unresolved influences by random processes can allow these models to exhibit surprisingly rich new behaviours of which their deterministic counterparts are incapable.
1. Introduction: model simplification and parametrization
Mathematical models of the climate system vary enormously in their levels of detail and complexity. The simplest are zero- or one-dimensional representations of the planet, and they calculate explicitly the dynamics of only a small number of (perhaps only one) prognostic variables (Schneider & Dickinson 1974; Saltzman 1978). At the other extreme are modern general circulation models (GCMs) of the coupled atmosphere–ocean system, which include comprehensive three-dimensional representation of the geophysical fluid dynamics in terms of millions or tens of millions of prognostic variables. Clearly the more complex models resemble the actual physical system more closely, whereas the more simplified models represent the climate system in only a generalized, or conceptual, way. However, because these conceptual models are sufficiently simple for their workings to be readily understood, they can be useful for gaining qualitative understanding of both the climate system, and the relationship between the models and the real system being modelled.
Any mathematical model of the climate requires simplifications relative to the real climate system. Of course, these simplifications are more numerous and severe for the simpler models. But regardless of the level of modelling detail, there are inevitably relevant processes occurring on time and space scales that are too small to be represented explicitly, which are not fully and quantitatively understood or both. The effects of such unresolved climate variables and climate processes are included in models in an approximate way, through functions—called parametrizations—of the resolved prognostic variables that summarize the effects of the unresolved processes. Fundamentally, these parametrizations are meant to represent the statistical characteristics of the effects of the unresolved processes on the model's prognostic variables.
By far the most common form for these parametrizations is representation of only the mean effects of the unresolved processes, with their other statistical properties being ignored. For example, figure 1 shows the relationship between the effects of unresolved processes on a resolved variable X in the idealized Lorenz (1996) model. The cloud of points represents individual ‘data’ values for the effects of unresolved variables on the time tendency of the resolved variable X, and the smooth curve is a regression function specifying (parametrizing) their mean effect, conditional on particular values of X. Such a parametrization could be implemented using an equation of the form(1.1)where R(X) represents the explicitly resolved dynamics and U(X) would be the smooth function of X indicated in figure 1, representing unresolved processes.
A conventional mean–response parametrization function such as U(X) in equation (1.1), which specifies a single conditional mean effect due to the unresolved scales, is deterministic; i.e. equation (1.1) contains no randomness, even though (according to the data scatter in figure 1) there is no unique or constant response from the unresolved scales for a given value of the resolved variable. By contrast, nature would provide a realization from the point cloud around U(X), rather than the value U(X) itself. Deterministic parametrizations have been justified heuristically by analogy to the macroscopic (‘resolved scale’) temperature response (of, say, the sensing surface of a thermometer) to a very large number of (‘unresolved’) interactions of the molecules of the medium. Because the number of microscopic interacting molecules is very large, the temperature is practically constant in time, even though the individual molecular kinetic energies, whose average effect yields the temperature, are highly variable. This conventional deterministic approach to parametrization has been called the ‘Reynolds/Richardson paradigm’ (Palmer et al. 2005), reflecting its historical genesis at the very beginning of modern fluid dynamical weather prediction.
This concept of a ‘thermodynamic limit’ and its applicability to the parametrization problem can be viewed in terms of the well-known central limit theorem from statistics (e.g. Lindgren 1968). One consequence of the central limit theorem is that the frequency distribution of the average of n independent random quantities (regardless of the distribution from which they have been drawn) tends towards the Gaussian (‘bell curve’) distribution as the number of averaged random effects increases. Furthermore, the variance of this distribution for the average decreases in proportion to 1/n. Craig & Cohen (2006) derived these same properties using a statistical mechanics approach, in which the unresolved small-scale influences are mass fluxes produced by independent (non-interacting) convective clouds in equilibrium with their large-scale forcing.
The solid curves in figure 2 illustrate these properties for the case in which the individual independent quantities to be averaged are drawn from an exponential distribution, shifted to the left so as to have zero mean, as shown in figure 2a. In this illustration, the probability density function in figure 2a represents the distribution for the effect of a single small-scale, unresolved event on a large-scale time tendency; figure 2b represents this distribution for the average effect of 10 such small-scale events. The dispersion of the probability distributions in figure 2 around their mean of zero corresponds to the scatter of points in figure 1 around the deterministic function U(X), for a particular value of X.
As progressively larger numbers of independent draws from the probability distribution in figure 2a are averaged (figure 2b–e), the shapes of the frequency distributions for their average effect approaches the Gaussian relatively quickly, even though the original distribution in figure 2a is very far from the Gaussian. On the other hand, extremely large numbers of these small-scale events must be averaged in order for the uncertainty regarding their net effect to vanish for practical purposes, which condition would be required for a deterministic (mean effect, only) parametrization of their influence to be justified. In practice, this zero-variance limit for parametrizations is not reached in climate models: if it could be, then the small-scale influences would be resolved, in effect, requiring no approximation of their contribution to the dynamics.
Note that independence of the small-scale variables is not necessary in order for the probability distribution for their average to approach the Gaussian, nor for the variance of that distribution to tend towards zero as the average is taken over progressively larger numbers of small-scale values. The probability densities in figure 2 (dashed curves) indicate corresponding results for moderately strongly correlated (ρ=0.8 for consecutive values in the n-member series) small-scale variables (Kotz & Adams 1964). These show that the distributions for the average of correlated values also tend towards the Gaussian, with decreasing variance as more values are averaged, but for both of these features the convergence is slower than in the case of independent small-scale effects. If the average effect of mutually correlated small-scale events on a large-scale tendency is to be represented by a parametrization, use of a deterministic (i.e. zero variance) parametrization is even less well justified.
2. Realistic variability in simple models: Hasselmann's paradigm
In a landmark paper, Hasselmann (1976) proposed that conventional mean–response parametrizations be extended to include also the variance associated with unresolved processes.
Integration of full GCMs of the climate system was only just beginning to become computationally feasible at that time, so most mathematical climate studies employed relatively simple conceptual models (Schneider & Dickinson 1974; Saltzman 1978). Among other attributes, these conceptual climate models highlight clearly the extreme simplification that may result from the use of deterministic parametrizations. Because most of these simple models do not exhibit chaotic dynamics, their variability is limited to adjustment towards their equilibrium solutions from non-equilibrium initial conditions, execution of fixed limit cycles (i.e. deterministic periodic solutions) or responses to variable external forcing. Attempts to simulate observed climate variability with deterministically parametrized conceptual models using realistic levels of external (e.g. orbitally varying solar) forcing had not been successful.
Hasselmann's insight was to realize that if the time scales of the resolved (‘climate’) variables and the to-be-parametrized unresolved (‘weather’) variables are well separated, then the central limit theorem implies that the effects of the unresolved scales on the resolved prognostic variables in a climate model can be represented as Gaussian white noise. In this paradigm, the governing equations for the climate model become stochastic differential equations, explicitly including randomness that represents effects of unresolved processes. Furthermore, this randomness may be specified using only its first two statistical moments because a Gaussian distribution is fully defined by its mean and variance. Hasselmann's paradigm extends the deterministically parametrized equation (1.1) to the stochastic differential equation(2.1)where σ is the square root of the Gaussian variance characterizing the uncertainty of the parametrized effects and z represents (unit variance) Gaussian white noise. When numerically integrated in discrete time, the effect of the term σz(t) is to randomly choose a value from a Gaussian distribution centred on U(X). This random variate corresponds to an individual value from the point cloud vertically above or below the relevant point on U(X) in figure 1. Although others had been thinking along similar lines at the time (Leith 1975; Lorenz 1975), Hasselmann (1976) was the first to provide a quantitative paradigm for stochastic parametrization in climate modelling.
Not surprisingly, a stochastically parametrized model of the sort shown in equation (2.1) can exhibit much more variability than its deterministic counterpart, equation (1.1). Hasselmann (1976) points out that, both conceptually and mathematically, a simple climate model with stochastic parametrization is analogous to Brownian motion, in which pollen grains exhibit (microscopic) motions that are the integral responses to their interactions with the molecules of the fluid in which they are suspended. In this analogy, the resulting random motions of the pollen grains might correspond to the variations in the average Earth temperature on the time scales of centuries or millennia, in response to the integrated effects of random synoptic (weather) disturbances in the general circulation of the atmosphere. Integration of a GCM of the atmosphere–ocean system would correspond to simulation of the dynamics of the individual fluid molecules, given fixed or only slightly perturbed positions of the pollen grains. Climate simulation using a simplified model having only non-stochastic parametrizations would correspond (in Hasselmann's words) ‘to determining the large-particle paths by considering only the interactions between the large particles themselves and the mean pressure and stress fields set up by the small-particle motions…’. Of course, the dynamics in this last analogy are drastically impoverished.
In one of the first studies to implement Hasselmann's paradigm, Lemke (1977) augmented a simple and well-studied zonally (i.e. around latitude circles) and annually averaged conceptual climate model (Budyko 1969; Sellers 1969) with a stochastic component. The conventional deterministic version of this model computes, for each resolved latitude, the energy balance as the sum of incoming and outgoing radiation (both as functions of the prognostic temperature), and an averaged meridional (i.e. north–south) heat flux convergence due to the net effect of advection by synoptic waves (parametrized as being proportional to the difference between the prognostic temperature at the latitude in question and the globally averaged temperature). These radiative and advective processes correspond to the terms R(T) and U(T), respectively, in equations (1.1) and (2.1). In addition, unresolved deviations from the average parametrized meridional heat flux divergence were represented by white noise forcing (σz(t) in equation (2.1)). Subsequently, Lorenz (1979) justified stochastic, non-diffusive meridional heat flux parametrizations on the basis of atmospheric observations.
A key quantity in equation (2.1) is the standard deviation (‘amplitude’) of the parametrized random synoptic forcing, σ. Setting this standard deviation equal to the average mid-latitude meridional heat flux (matching to the order of magnitude the available observations of heat transport by synoptic waves), Lemke (1977) obtained power spectra for the prognostic temperatures having the form shown in figure 3. Broadly consistent with observed climate spectra, apart from peaks attributed to external (prominently, orbital) variations in solar input, these spectra are flat at low frequencies and exhibit approximately an ω−2 dependence at higher frequencies. Depending on the choice of the magnitude of the modelled albedo dependence on temperature, Lemke (1977) generally found the transition between these two portions of the spectra to occur between 10−2 and 10−3 per year.
3. Stochastic parametrization and regime dynamics
Besides providing additional and realistic variability in simple models, stochastic parametrizations can also produce surprising new behaviours that are qualitatively different from those of their deterministically parametrized counterparts. This characteristic of stochastic parametrizations is illustrated here using a conceptual model for ice ages, which has been abstracted from Nicolis & Nicolis (1981), Sutera (1981), Benzi et al. (1982) and Nicolis (1982).
The model is a highly idealized zero-dimensional model for the global radiative energy balance, with the single prognostic variable being the globally and annually averaged temperature T(3.1)where C is an effective heat capacity. The Earth loses energy to space according to the well-known physics of blackbody radiation, in proportion to the fourth power of the prognostic temperature variable. Thus, the rate of energy loss −Eout corresponds to R(T) in equation (1.1). The factor ϵ is an average effective global emissivity (approx. 0.6). The rate of energy input is given by the proportion 1−α(T) of the annually and globally averaged solar flux Q (approx. 342 W m−2) that is absorbed, where α(T) is the average planetary albedo parametrized as a function of the prognostic temperature. Thus, the energy gain term Ein might be associated with the deterministic parametrization U(T) in equation (1.1).
One view of the workings of this simple deterministic model is shown in figure 4a, where Ein (dashed curve) and Eout (solid curve) are plotted as functions of temperature. The parametrized albedo function α(T) has been chosen in a way that Ein crosses the Eout curve in three places, yielding the three equilibrium solutions T−, T0 and T+. The solutions T− and T+ are stable equilibria because dT/dt>0 for T<T− and T0<T<T+, and dT/dt<0 for T−<T<T0 and T>T+. Thus, for any initial Tinit≠T0, the dynamics in equation (3.1) will ultimately arrive at one or the other of the stable solutions T− or T+. Although dT/dt=0 at T=T0, the slightest perturbation from this unstable solution will also tend towards one of the two stable solutions.
Another view of the dynamics and stability of this simple model is shown by the potential function Ψ in figure 4b, which is related to equation (3.1) through(3.2)The dynamics of equation (3.1) can be visualized by imagining a (two-dimensional) laboratory apparatus having the shape of the potential function in figure 4b, on which a steel ball can roll under the influence of gravity. Because the density of the steel is high, any effects of air turbulence in the laboratory have negligible effect on its motions. Unless the ball is initially balanced rather precisely at T0, it will promptly roll into one or the other of the potential wells (or ‘basins of attraction’) and stop at either T− or T+. From that point forward, dT/dt=0, and the model climate exhibits no further temperature variability, even though two temperature regimes, corresponding to the two potential wells, are physically possible.
Consider now this same simple climate model that has been extended with a stochastic parametrization to include variations in the global energy balance, which are not related to the global albedo according to the function 1−α(T) shown in figure 4a,(3.3)Because the Gaussian noise z(t) is additive, its effect does not depend on T, so that the potential function (equation (3.2)) for both equations (3.1) and (3.3) is the same. The dynamics expressed by equation (3.3) correspond to an air-filled, low-density ball rolling on the apparatus illustrated in figure 4b. Furthermore, the laboratory windows are open, so the motions of the ball are subject to turbulent accelerations that increase in intensity with increasing σ.
In this new situation, the ball still tends to roll towards the two stable solutions T− and T+, but nearly always remains in motion due to the erratic turbulent accelerations. If the wind is light (small σ), the random accelerations are small compared with the gravitational acceleration propelling the ball towards the bottoms of the wells, and the ball will usually reside very near T− or T+. In this small-σ situation, the ball will almost never switch from one well to the other. For moderate σ, the ball will move fairly frequently between the two wells, but will reside relatively near their centres more often than between or outside of them. If the wind is very strong (large σ), then the topography of the potential function apparatus makes little difference to the position of the ball, so the probabilities for its location are almost independent of position.
Figure 5 shows the probability density functions f(T) for the location of the ball in these three situations. In terms of the climate model, these correspond to the climatological distributions for globally averaged temperature, for the three levels of random forcing, σ. Quantitatively, these probability density functions could be computed according to (Nicolis & Nicolis 1981)(3.4)which becomes progressively flatter and more independent of the potential function as the variance of the stochastic forcing increases.
The double-peaked probability function for small σ in figure 5 may seem surprising since it implies occasional jumps between the two regimes, each of which requires a relatively large stochastic perturbation or sequence of perturbations. Such events become progressively rarer as σ decreases, but are not impossible unless σ=0 (i.e. for the deterministic climate model in equation (3.1)). It can be useful to characterize this tendency of the modelled system to jump between regimes using the average, or typical, time to achieve such a jump (Nicolis & Nicolis 1981; Sutera 1981)(3.5)where ΔΨ is the difference between the bottom of the potential well currently occupied and the local maximum Ψ(T0), specifying the height of the ‘barrier’. For sufficiently small σ2, this characteristic time could be large enough (in comparison with, say, the age of the Earth) that the transitions are theoretical possibilities only.
Simple climate models of the sort exemplified by equation (3.3) and figure 4a,b were originally constructed for conceptual investigation of the approximately 105-year ice age cycle. In this context, T+≈288 K might be associated with interglacials and T−≈278 K with ice age conditions. Jumps between the two regimes would then correspond to oscillations between these two climate states. However, working with a stochastically parametrized model qualitatively similar to equation (3.3), Benzi et al. (1982) found a purely red temperature spectrum, with no hint of ice age-like periodicities.
Consider now a modification to this simple climate model, in which the incoming solar radiation, Q, in equations (3.1) and (3.3), varies periodically in a way that mimics the Milankovich eccentricity cycle(3.6)where Q0 is the long-term average of the incoming solar flux and ω=10−5 per year. Figure 4c,d indicates the effect on equations (3.1) and (3.3) of a reduction of solar flux that might occur for t≈75 000 year. The effect is to lower the dashed curve for Ein in figure 4c relative to its level in figure 4a (which is reproduced in figure 4c with the grey dashed curve, for comparison). Intersections of the modified Ein function with Eout now define stable equilibria T− and T+ that are colder than their counterparts in figure 4a, while the unstable equilibrium solution T0 increases towards T+. Concurrently, the cold potential well Ψ(T−) becomes deeper relative to the well for the warm regime Ψ(T+). To see why this change in the shape of the potential function will occur, consider the situation where Ein is reduced to the point that it only just touches Eout at a warm equilibrium T+. At that point T0 merges with T+ and the previously existing ‘warm’ potential well shallows sufficiently for that local minimum in Ψ to disappear. The reverse effects occur for increases in solar input, e.g. at t≈25 000 year.
Figure 6a illustrates the cycle of potential functions that follows from periodic variation in the external solar forcing specified by equation (3.6). As solar input increases, the stable equilibrium temperatures T− and T+ increase, while the depth of the well representing the warm regime increases at the expense of the cold regime well (π/2). As the solar input subsequently decreases, the topography of the potential functions passes through its initial configuration (π), following which the stable equilibria cool further (3π/2), with the cold regime potential well deepening at the expense of the warm regime. The overall effect is that the nature of the model dynamics is being affected by the changes in the external solar input.
Imagine the deterministic dynamics of equation (3.1) as influenced by the periodic solar forcing, again in terms of the laboratory analogue. The track along which the steel ball may roll changes shape periodically, as shown by the cycle in figure 6a. If the ball begins at T=T+ at t=0, it will execute regular periodic excursions around this initial temperature, with period 1/ω=105 year, as the location of the maximum depth of the warm well periodically increases and decreases. However, the steel ball remains trapped in whichever well it begins, and its maximum temperature excursions are substantially smaller than the range of temperature variation between ice ages and interglacials, T+−T−. The most prominent climatic temperature excursions in recent geological history have occurred with a period of approximately 105 year, but the roughly 0.1% amplitude of the eccentricity variation in solar radiation is the smallest of the three Milankovich cycles. It appears to be inadequate to explain the ice ages, not only in this very simple model but also in other more realistic models.
Qualitatively different behaviour occurs when the stochastically parametrized model (equation (3.3)) is subjected to the periodic solar forcing of equation (3.6). In particular, the stochastic forcing σz(t) now allows solutions to make transitions between the interglacial and ice age regimes, resulting in probability distributions for T that range over a broad range of possibilities. Furthermore, according to equation (3.4) and as shown in figure 6b, these climatological temperature distributions f(T) execute a periodic cycle that follows the 105-year cycle of the potential function in figure 6a. At times of higher solar forcing (π/2), the warmer climate (around T+) becomes the more favoured state. Similarly, the colder climate (around T−) becomes more likely at times of lower solar forcing (3π/2).
The existence of these periodically shifting probability distributions is not by itself sufficient for the solutions T(t) of equations (3.3) and (3.6) to exhibit a realistic behaviour (e.g. a prominent spectral peak near ω=10−5 per year). In addition, the variance of the stochastic forcing in equation (3.3) must be of an appropriate magnitude, exemplifying a phenomenon that has come to be known as ‘stochastic resonance’ (Benzi et al. 1982). This is not a resonance in the usual sense of a coupling between two preferred frequencies, because the spectrum of the white noise forcing z(t) is completely flat: it has no preferred frequencies. Rather, stochastic resonance occurs when the variance σ2 of the stochastic forcing yields a characteristic transition time τ (equation (3.5)) that matches approximately the period ω−1 of the external forcing. Stochastic forcing that is too weak will result in transitions that occur too rarely, yielding solutions similar to those for the deterministic dynamics (equation (3.1)), although with a small amount of superimposed variability. Stochastic forcing that is too strong will tend to produce many more than two regime transitions over the course of a single period of ω−1 year.
When τ and ω−1 are well matched, the random transitions will tend to occur with the progression of the solar forcing cycle as the potential barrier ΔΨ decreases, while the simultaneous deepening of the preferred potential well will tend to retain the solution in the (new) preferred regime until the topography of the potential function changes again. The resulting time series for modelled temperature qualitatively resemble transitions between ice ages and interglacials (Benzi et al. 1982; Imkeller 2001), with the temperature response lagging the solar forcing by approximately π/4 (Nicolis 1982). This capacity for appropriately scaled random perturbations to amplify weak periodic forcings through stochastic resonance was first noted by Benzi et al. (1982) in the context of this conceptual climate model and has since seen numerous applications in various areas of physics (e.g. Gammaitoni et al. 1998).
The issue of regime switching among multiple climate equilibria is broader than the particular application to ice ages presented in this section, and conceptual models with stochastic parametrizations have been applied to regime behaviour in other settings as well. These include simple models of multiple equilibria in the North Atlantic thermohaline circulation (Stommel & Young 1993; Cessi 1994; Monahan 2002; Monahan et al. 2008) and in the large-scale atmospheric circulation (Egger 1981; Benzi et al. 1984; De Swart & Grassman 1987). Stochastic parametrization has also been found to improve the resemblance to the well-known three-dimensional Lorenz (1963) model of a two-dimensional principal components truncation of it (Selten 1995; Palmer 2001), in terms of both the individual trajectories and the overall bimodal ‘climatological distribution’.
This paper has reviewed a modest selection of the earlier and conceptually simpler work on stochastic parametrization of climate models, both for a historical perspective and as a pedagogical device. Not surprisingly, a prominent effect of including stochastic variability in parametrizations is to increase the variability exhibited by the models overall, to more realistic levels. This effect has been found also for stochastically parametrized versions of more complete and complex climate models (e.g. Lin & Neelin 2000; Sardeshmukh & Sura 2007). More recent work on stochastic parametrization in conceptual models has been summarized by Imkeller & von Storch (2001) and Imkeller & Monahan (2002).
The simplest approach to stochastic parametrization is through additive white (i.e. temporally independent) noise. A prominent time-scale separation between the resolved and parametrized processes implies that white noise will be a reasonable representation because successive averages of an autocorrelated process over sufficiently long periods (although still shorter than the resolved time scale) will be approximately independent of each other. When this condition is not met, representation of the unresolved parametrized effects using an autocorrelated stochastic process presents no problem in simulation (e.g. Egger 1981; Benzi et al. 1984; De Swart & Grassman 1987; Monahan 2002; Wilks 2005), although it does substantially complicate the calculation of analytical mathematical results (e.g. Egger 1981; Stommel & Young 1993). Similarly, in spatially explicit models, it has been found that spatially correlated random forcing may be appropriate (e.g. Lemke 1977; Penland & Sardeshmukh 1995; Buizza et al. 1999).
An important conclusion to be drawn from the very simple studies summarized here is that stochastic influences may fundamentally affect the simulation of regime behaviour in climate models, and in particular that rates and probabilities for regime transitions depend on the presence and magnitude of the stochastic input. Realistic modelling of climate changes may depend critically on realistic modelling of these regime changes. Palmer (1999) has argued that changes in the climate (including those due to anthropogenic influences) may manifest themselves in large part through changes of residence frequencies in its preferred regimes, noting that models with insufficient internal variability tend to exhibit overpopulation of the more stable regimes (Molteni & Tibaldi 1990). A more unsettling view is that climate change may manifest itself through an abrupt shift to a qualitatively different regime (e.g. Broecker 1997).
One contribution of 12 to a Theme Issue ‘Stochastic physics and climate modelling’.
- © 2008 The Royal Society