Analysis and Modeling of the Natural Variability of Climate (original) (raw)
1. Introduction
There is widespread agreement that the obliquity and precession of the earth’s orbit result in a significant portion of global temperature variance at frequencies near 40 000 and 20 000 yr. Aside from this “external” forcing of the climate system, there is no consensus on the physical processes governing the remainder of the variance. It is generally agreed that the variance resulting from orbital forcing is superimposed on a continuous “background” spectrum with a variance that increases with decreasing frequency (Kerr 1978). Most attempts to model this background spectrum assume a nonlinear response of the climate system to the orbital forcing (Imbrie and Imbrie 1980; Bhattacharya et al. 1982; Le Treut and Ghil 1983; Nicolis 1993; Hagelberg et al. 1994). Each of these studies has succeeded in reproducing only a limited portion of the spectrum of climatic variations.
Several studies have recognized that the low-frequency portion of the power spectrum, S(f), of temperature variations inferred from ice and ocean cores has the form of a Lorentzian distribution (Hasselmann 1971; Komintz and Pisias 1979; Lovejoy and Shertzer 1986):
This function is constant for very low frequencies and proportional to _f_−2 for f ≫ _f_0. Hasselmann (1971) and Komintz and Pisias (1979) have suggested that a stochastic model may be the most appropriate model for this spectrum since a Lorentzian spectrum can result from the sum of uncorrelated random pulses with a negative feedback mechanism limiting the variance at low frequencies.
At shorter timescales it is recognized that there is “persistence” or correlations in time in meteorological time series over a range of timescales. Persistence means that warm years (or weeks or months) are, more often than not, followed by warm years and cold years by cold years. Hurst et al. (1965) presented studies of these correlations using the rescaled-range technique (see Feder 1991 for an introduction). Hurst found that time series of annual mean temperature produced a power-law rescaled-range plot with an average exponent of 0.73. No persistence would yield an exponent of 0.5. Numerical studies have shown that series with a Hurst exponent of 0.73 have a power spectrum proportional to _f_−1/2 (Fox 1989; Higuchi 1990; Gomes da Silva and Turcotte 1994).
Understanding the natural variability of climate is one of the most important tasks facing climatologists. The 1995 Intergovernmental Panel on Climate Change (IPCC) report (1995) concluded that the “balance of evidence suggests a discernible human impact on the climate system.” This conclusion is based, however, on comparison with the natural variability of the climate system expected from general circulation model (GCM) control runs. GCM control runs exhibit significantly lower variability and a different frequency dependence than paleoclimatic data (IPCC 1995). Understanding more about the natural variability of climate is essential for an accurate assessment of the human influence on climate.
In this paper we present a model that provides a specific physical mechanism for the continuous spectrum from timescales of 1 day to 200 kyr. The model is based on internal, “natural” variability, produced by the stochastic nature of atmospheric turbulence. We present spectral analyses of paleoclimatic proxy data and instrumental data that support the observation of a low-frequency Lorentzian spectrum and a higher-frequency spectrum proportional to _f_−1/2. At very high frequencies we found a distinct difference between the power spectra of continental and maritime stations. Continental stations exhibit power spectra proportional to _f_−3/2 at timescales less than 1 month while the spectra of maritime stations remain proportional to _f_−1/2 down to timescales of 1 day. Our model is based upon an analytic approach to modeling the stochastic diffusion of heat in the atmosphere and ocean. The difference between continental and maritime stations results from the air mass above maritime stations exchanging heat with both the atmosphere above and the ocean below while the air mass above continental stations exchanges heat with only the atmosphere above it.
The model we present was first solved by van Vliet et al. (1980) to determine the power spectrum of variations due to the stochastic diffusion of heat in a metallic film in thermal equilibrium with a substrate. Temperature variations in the film and substrate occur as a result of fluctuations in the heat transport by electrons undergoing Brownian motion. The top of the film absorbs and emits blackbody radiation. In this paper we use van Vliet et al.’s model exactly as they presented it with the atmosphere as the metallic film and the ocean as the substrate. Turbulent eddies in the atmosphere are analogous to the electrons undergoing Brownian motion in a metallic film in contact with a substrate. The model studied by van Vliet et al. (1980), with thermal and diffusion properties appropriate to the atmosphere and ocean, yield a power spectrum with the same form as that of climatic variations recorded in instrumental records and inferred from ocean and ice cores from timescales of 1 day to 200 kyr. However, there is one major discrepancy between the model and observation: the atmospheric eddy diffusivity required to match the Vostok data is much smaller than the most reliable estimates.
2. Observations of climatic variations
In Fig. 1 we present of the power spectrum of the variations in temperature as a function of the frequency predicted by the model and supported by observational data. At low frequencies the power spectum is constant. Above f ≈ (40 kyr)−1 the power spectrum is proportional to _f_−2. Above f ≈ (2 kyr)−1 the power spectrum is proportional to _f_−1/2. At very high frequencies (above f ≈ (1 month−1) the spectrum varies as _f_−3/2 for continental stations and remains proportional to _f_−1/2 for maritime stations. The crossover frequencies quoted above are those observed in the climatological record.
Figure 2 shows the Lomb periodogram of the Vostok δ D record extending back 220 kyr converted into degrees Celsius by the conversion factor 5.6% (°C)−1 (Jouzel et al. 1987). We obtained the spectrum from M. Ghil and P. Yiou (1994, personnal communication). Jouzel and Merlinvat (1984) have concluded that the Vostok deuterium record is a proxy for local atmospheric temperature. It is not possible to directly estimate the power spectrum of the Vostok record with the FFT since the data are unevenly sampled. Press et al. (1992) suggest the use of the Lomb periodogram for such data. The periodogram shows a constant low-frequency region that changes to a region proportional to _f_−2 above f 1 (40 kyr)−1. Above f ≈ (2 kyr)−1 the power spectrum of temperature variations is ∝ _f_−1/2. To obtain a more consistent estimate of the power spectrum at high frequencies, we averaged the periodogram in logarithmically spaced bins above f = (10 kyr)−1.
In Fig. 3 we present the power spectra of temperature time series as a function of the frequency taken instrumentally at higher frequencies. We plot the average power spectrum of time series of monthly mean temperature from 94 stations worldwide with the yearly trend removed. The power spectra were computed by taking the modulus squared of the complex Fourier coefficients obtained by computing the fast Fourier transform using the Numerical recipes routine realft (Press et al. 1992). We computed the power spectra of all complete temperature series of length greater than or equal to 1024 months from the climatological database compiled by Vose et al. (1992). The yearly trend was removed by subtracting from each monthly data point the average temperature for that month in the 86-yr record for each station. All of the power spectra were then averaged at equal frequency values. The data yield a straight-line least squares fit with slope close to −0.5 indicating that S(f) ∝ _f_−1/2 in this frequency region. This is consistent with the results of Hurst et al. (1965), who analyzed correlations at instrumental timescales with the rescaled-range method. Numerical studies have shown that a time series that yields a power-law rescaled-range plot with Hurst exponent H = 0.73 has a power spectrum S(f) ∝ _f_−1/2 (Fox 1989; Higuchi 1990; Gomes da Silva and Turcotte 1994). Here, _f_−1/2 power spectra have also been observed in variations of atmospheric humidity on timescales of days to years (Vattay and Harnos 1994) and in tree-ring widths, river discharge, and precipitation up to several thousand years (Pelletier and Turcotte 1996, manuscript submitted to J. Hydrol.).
In Figs. 4a and 4b we present the average power spectrum of time series of daily mean temperature (estimated by taking the average of the maximum and minimum temperature of each day) from 90 maritime and 1000 continental stations, respectively, over 4096 days estimated by computing the fast Fourier transform as before. Maritime stations were sites on small islands far from any large landmasses. We chose all stations that fit this criterion and had no more than 2% missing record from the Global Daily Summary database compiled by the National Climatic Data Center (1994). Continental stations were any stations not on small islands. We filled in the missing days with the same value as the previous day. The yearly periodicity in these data were removed by subtracting from each station’s temperature time series a converged least squares fit [using Numerical recipes (Press et al. 1992) routine mrqmin] to a function of the form
where i is the number of the day in the year and T av, A, and ϕ were the fitting parameters. This procedure is a standard one for subtracting the yearly periodicity in a meteorological time series (Janosi and Vattay 1992). Continental stations exhibit a high-frequency region close to _f_−3/2. Maritime stations exhibit _f_−1/2 scaling up to the highest frequency.
A number of investigators have applied techniques of analysis developed in dynamical systems theory to meteorological time series. One of the most striking observations has been the calculation of low correlation dimensions from meterological and paleoclimatic time series [Nicolis and Nicolis 1984; see Islam et al. (1993) for a recent summary of these studies]. The correlation dimension is purportedly a measure of the number of degrees of freedom of the dynamical system (Nicolis and Nicolis 1984). That meteorological and climatic time series yield a low correlation dimension [dimensions between 3 and 5 are common values (Islam et al. 1993)] is generally considered to be a very surprising result since the dynamics of climate are assumed to be extremely complex, involving a large number of interacting components. The widespread belief that a low correlation dimension is indicative of a low-dimensional chaotic dynamical system has resulted in an emphasis on modeling the climate system as a forced nonlinear oscillator (Yiou et al. 1994). This belief is based on the assumption that stochastic models do not yield finite correlation dimensions as is the case for white noise. However, it has been shown that stochastic time series with power-law spectra yield low correlation dimensions (Osborne and Provenzale 1989). In particular, time series resulting from turbulent dispersion have been shown to yield a finite correlation dimension (Sanderson and Goulding 1990). Therefore, the stochastic model presented in this paper is consistent with the observation of low correlation dimensions for meteorological and paleoclimatic time series.
3. Parameterization of convective transport
The simplest approach to the turbulent transport of heat is to assume that the flux of heat is proportional to its gradient. The mixing of heat energy will then be governed by the diffusion equation. This parameterization is known as the flux-gradient parameterization with first-order closure (Garratt 1981). This approximation works well if the turbulence is dominated by small-scale eddies (much smaller than that of the mean gradient of potential temperature) where the convective action of turbulent eddies will be analagous to the molecular collisions responsible for molecular diffusion (Moffatt 1983). In the atmosphere and ocean these small-scale mixing processes are superimposed on a large-scale advection (Hadley and Walker circulation cells in the atmosphere and bottom water formation and large-scale upwelling in the ocean). This assumption is an accurate one for a neutral or stably stratified troposphere, but breaks down for a highly convective one (Garratt 1992). For studies of long-term climate change, the relevant transport mechanism or stratification is that which dominates in a long-term, global average. Diffusive transport has been assumed in analytical models of climate change in the atmosphere (Ghil 1983; North 1975). Both steady-state and transient tracer data in the atmosphere and ocean are well-modeled by a diffusion approximation with large-scale advection. For the ocean, Oeschger et al. (1975) and Munk (1966) have argued that diffusion combined with global-scale advection models the steady-state profiles of carbon and temperature, respectively. The vertical dispersion of tritium input into the oceans as a result of atmospheric nuclear testing is also well described by a diffusion model (Garrett 1984). In the atmosphere, the eruption of Mount Pinatubo and other recent volcanic events have provided us with a global tracer akin to bomb tritium in the ocean. Hofmann and Rosen (1987) have argued that vertical diffusion is the best model of the exponential decay of aerosol from El Chichon with a diffusivity of D ≈ 1 m2 s−1. Estimation of the vertical diffusion coefficient based on the steady-state distribution of aerosol in the middle atmosphere obtains a roughly constant diffusivity within the upper half of the troposphere (Parameswaran et al. 1993). Kida (1983) has traced the dispersion of air parcels in a hemispheric GCM and found the transport to be diffusive in the troposphere.
To test the hypothesis that advection-diffusion is a realistic first-order model of large-scale atmospheric transport, we have performed a correlation analysis of temperature and precipitable water vapor concentration vertically in the atmosphere measured with the TIROS operational vertical sounder (TOVS). We obtained 2 months (January–February 1985) of daily TOVS gridded temperature and water vapor concentration data with 1° horizontal spatial resolution at pressure levels 1000, 850, 700, 500, 400, 300, and 200 mbar. The time series were detrended with a least squares linear fit to remove the annual trend and leave only the high-frequency variations. For a fixed horizontal location and a fixed time lag, t o, the detrended temperature and water vapor concentration at one elevation, _h_1, was multiplied by that at another elevation, _h_2. These products were averaged for all pairs of data separated by the same distance r = |_h_1 − h_2| normalized by the product of the standard deviations of both time series to obtain the spatial cross-correlation function (for temperature fluctuations), Δ_T:
The cross-correlation function for a one-dimensional diffusion process is given by the solution to a one-dimensional diffusion process with a constant velocity advection, with velocity u, and diffusivity D (Voss and Clarke 1976):
The spatial cross-correlation as a function of vertical distance for time lags t = 1, 2, 3, 4, 5, and 6 days is presented in Figs. 5a and 5b for the temperature and precipitable water vapor concentration, respectively. The functions given by Eq. (1) for u = 0.01 m s−1 and D = 100 m2 s−1 are presented in Fig. 5b. The advection-diffusion correlation functions of Fig. 7b bear a close resemblance to the spatial cross-correlations of the TOVS data. The advection term is consistent with large-scale advection by Hadley and Walker circulations. Superimposed on the transient eddies that we model as diffusive there is a steady-state vertical advection of heat and water vapor by Hadley and Walker circulations. The globally averaged velocity of local vertical advection is 0.003 m s−1 (Peixoto and Oort 1992), close to the observed globally averaged velocity, 0.01 m s−1, from the TOVS analysis presented in Fig. 5. The diffusivity of D ≈ 100 m2 s−1 was with a least squares linear fit of ln[c(r, 1)] versus _r_2. We conclude from this analysis that an advection-diffusion model is a reasonable first-order model of large-scale atmospheric transport of heat and water vapor. In our modeling we will focus on the effect of internal noise and diffusive transport only on the local fluctuations of heat and water vapor. The advective transport would lead to only a periodicity in the time series of local temperature or water vapor with an average period of about 1 month (the return time for tropospheric circulation with u = 0.01 m s−1). No such periodicity was identified in the time series presented in section 2. Therefore it is most likely dominated by the stochastic variability of the system. These results are consistent with the work of Lindzen (1990), who argued that water vapor is delivered to the upper levels of the troposphere primarily by transport associated with convective systems of relatively small spatial scale.
Due to the stochastic nature of turbulent flow, a tracer will disperse a different way when placed in a turbulent flow at different times, although the mean (ensemble averaged) concentration may be described by the diffusion equation. Thus, a deterministic diffusion equation is inadequate to model turbulent transport. We will model the fluctuations in turbulent transport by adding a stochastic term to the heat flux of a deterministic diffusion equation. This term represents variations in convective activity from place to place and time to time in the turbulent atmosphere. Fluctuations in heat transport may play a dominant role in the natural variability of climate. The spatial distribution of heat energy in the atmosphere and oceans is a time integration of the flux of heat energy. Fluctuations caused by the random action of turbulent transport could lead to anomalously large heat storage in the upper atmosphere and ocean resulting in low temperatures in the lower atmosphere. In this paper, we investigate the fluctuations in local temperature expected from a stochastic diffusion model of turbulent heat transport and compare our predictions with instrumental records.
4. Fluctuation model of natural climatic variability
A stochastic diffusion process can be studied analytically by adding a noise term to the flux of a deterministic diffusion equation (van Kampen 1981):
where Δ_T_ is the fluctuation in temperature from equilibrium, σ is the vertical thermal conductivity, ρ is the density, c is the specific heat per unit mass, and the mean and variance of the noise are given by
η x, t(7)
and
η x, t η x t σ x T x_2_δ x x δ t t
Novikov (1963) has advocated this method for studying turbulent fluctuations.
Perhaps the simplest calculation possible with this equation is the power spectrum of temperature fluctuations in a layer of width 2_l_ of an infinite, one-dimensional, homogeneous space. It is also directly applicable to the difference between the power spectral behaviors of maritime and continental stations at high frequencies. The presentation we give is similar to that of Voss and Clarke (1976). The variations in total heat energy in the layer of width 2_l_ are determined by the heat flow across the boundaries. Figure 6a illustrates the geometry of the layer exchanging thermal energy with diffusing regions above and below it. A diffusion process has a frequency-dependent correlation length λ = (2_D_/f)1/2 (Voss and Clarke 1976). Two different situations arise as a consequence of the length scale, 2_l_, of the geometry. For very high frequencies, λ ≪ 2_l._ In that case, the fluctuations in heat flow across the two boundaries are independent. For low frequencies, λ ≫ 2_l_ and the boundaries fluctuate coherently. First we consider high frequencies. Since the boundaries fluctuate independently, we can consider the flow across one boundary only. The flux of heat energy is given by Eq. (4). Its Fourier transform is given by
where α = σ/ρ c is the vertical thermal diffusivity. The flux of heat energy out of the layer at the boundary at x = l (the other boundary is located at x = −l) is the rate of change of the total energy in the layer E(t): [dE(t)/_dt_] = J(l, t). The Fourier transform of this equation is
Therefore, the power spectrum of variations in E(t), S E(ω) = [|E(ω)|2], is
where the brackets denote an ensemble average. Since Δ_T_ ∝ Δ_E,_ S T(ω) ∝ _ω_−3/2 also.
If we include the heat flux out of both boundaries, the rate of change of energy in the layer will be given by difference in heat flux: [dE(t)/dt_] = J(l, t) − J(−_l, t). The Fourier transform of E(t) is now
Then
where θ = [ω/(ω o)]1/2, and ω o = D/2_l_2 is the frequency, where the correlation length is equal to the width of the layer. When λ ≪ 2_l,_ the above expression reduces to S T(ω) ∝ ω_−3/2. When λ ≫ 2_l, S T(ω) ∝ _ω_−1/2 (Voss and Clarke 1976).
In section 2 we presented evidence that continental stations exhibit a _f_−3/2 high-frequency region and maritime stations exhibit _f_−1/2 scaling up to the highest frequency. This observation can be interpreted in terms of the diffusion model presented above. The power spectrum of temperature variations in an air mass exchanging heat by one-dimensional stochastic diffusion is ∝_f_−1/2 if the air mass is bounded by two diffusing regions and is ∝_f_−3/2 if it interacts with only one. The boundary conditions appropriate to maritime and continental stations are presented in Figs. 6a and 6b, respectively. The maritime stations have a _f_−1/2 power spectrum up to the highest frequency because the air mass above a maritime station exchanges heat with both the atmosphere above and the ocean below. In maritime stations, fluctuations in heat energy are readily absorbed by both the ocean below and the atmosphere above (which can radiate heat into space). The fluctuation calculation appropriate for maritime stations is one in which the coherent fluctuations from two boundaries are considered as in the calculation of the _f_−1/2 spectrum. The air mass above continental stations exchanges heat energy only with the atmosphere above it, the lower boundary condition acting as an insulator. The calculation appropriate for continental stations is the one-boundary model that predicts the observed _f_−3/2 spectrum. At low frequencies, horizontal heat exchange between continental and maritime air masses limits the variance of the continental stations. This occurs at a timescale of 1 month.
At timescales up to 2000 yr, fluctuations in heat energy input into air masses from the ocean or variations in emitted radiation can be absorbed by the other reservoir. At lower frequencies, the atmosphere and ocean achieve thermal equilibrium. The ocean can no longer absorb thermal fluctuations input into the atmosphere from fluctuations in radiative emission at this timescale and vice versa. The variance in temperature of the atmosphere and ocean is then determined solely by the radiation boundary condition. The equation governing the average behavior of the system at these low-frequencies (averaged over the noise) is
where _c_′ is the heat capacity per unit mass of the ocean, _ρ_′ is the density of the ocean, w_2 is the depth of the ocean, and g is the radiative damping constant, the coefficient of linear expansion of the Stefan–Boltzmann radiation law around the equilibrium temperature of the upper atmopshere: g = 4_σ
T_3_o
≈ 1.7 W (_m_2 K)−1. Two different phenomena can occur depending on the vertical thermal conductivity of the atmosphere. If the conductivity is low compared to the timescale for radiative damping, τ = _c_′_ρ_′_w_2/g, the rate-limiting step of radiative equilibrium is the conduction of heat from the ocean through the atmosphere. The time constant of this thermal discharge (analagous to the electrical discharge of an RC circuit) is then τ = (RC)−1 = σ/_w_1_w_2_c_′_ρ_′. However, if the conductivity is high, the timescale of radiative equilibrium is determined by the timescale for radiative damping, τ = _c_′_ρ_′_w_2/g. Using the most reliable estimates for the thermal and diffusion properties of the atmosphere and the ocean (which we will introduce next), the conductivity is large and the low-frequency portion of the spectrum is governed by the radiative damping constant. This is problematic, however, because the timescale of radiative equilibrium is then τ = _c_′_ρ_′_w_2/g = 600 yr. Such a short timescale implies that no _f_−2 spectral region should exist at all. This is a major discrepancy between the model and the observed spectrum for which there is no clear explanation.
To show these quantitative features of the model more explicitly, we consider a coupled atmosphere–ocean model with an atmosphere of uniform density (equal to the density at sea level) in thermal contact with an ocean of uniform density. The height of our model atmosphere is the scale height of the atmosphere (height at which the pressure falls by a factor of e from its value at sea level). Figure 7 illustrates the geometry and constants chosen (where σ is the vertical heat conductivity, ρ is the density, c is the specific heat per unit mass, α is the vertical thermal diffusivity, and g is the thermal conductance of heat out of the earth by emission of radiation). Primed constants denote values for the ocean. The physical constants that enter the model are the thermal conductance by emission of radiation and the density, specific heat, vertical thermal diffusivity, and depth of the ocean and atmosphere. The density and specific heat of air and water are well-known constants. We chose an ocean depth of 4 km and an atmospheric height equal to the scale height of 8 km as used by Hoffert et al. (1980) in their climate modeling studies. The eddy diffusivity we employ for the ocean is 6 × 10−5 m2 s−1. This value has been obtained from Tritium dispersion studies (Garrett 1984). The vertical eddy diffusivity for the atmosphere we use is 100 m2 s−1, as obtained in our analysis of temperature and water vapor fluctuations measured by TOVS. It must be emphasized that there is large variability in estimates of the vertical eddy diffusivity in the atmosphere. Some studies have estimated it to be two orders of magnitude smaller for certain conditions. Pleune (1990) and Seinfeld (1986) have quoted 1 m2 s−1 as their estimate in stable air conditions. This eddy diffusivity implies an equilibration time of the tropospheric air column to be 2 yr. This value is roughly consistent with the 1-yr decay time of the Pinatubo and El Chichon aerosols (Hofmann and Rosen 1987; Rosen et al. 1994).
Since the timescale of horizontal diffusion in the atmosphere and ocean is so much smaller than the timescale of vertical diffusion, diffusion of heat into and out of a local air mass is one-dimensional. For this reason, we consider only the variations in local temperature resulting from heat exchange vertically in the atmosphere and ocean.
The equation for temperature fluctuations in space and time in the model is Eqs. (5)–(8):
with
η x, t(16)
η x, t η x t σ x T x_2_δ x x δ t t(17)
North and Cahalan (1981) analyzed a similar model of climate change with respect to predictability. They studied the diffusion equation in two dimensions as a model for heat transport horizontally in the atmosphere. They included a white noise term on the right side of the diffusion equation {they used η(x, t) where we use [∂η(x, t)/∂_x_]} to represent variations in heat transport by turbulent eddies. However, a noise term in the flux of temperature (rather than in the temperature itself) is more appropriate as a model for variations in convective activity.
The boundary conditions are that there be no heat flow out of the bottom of the ocean and continuity of temperature and heat flux at the atmosphere–ocean boundary:
At the top of the atmosphere we will impose a blackbody radiation boundary condition. Most (65%) of the energy incident on the earth is emitted as long-wavelength blackbody radiation from the H2O and CO2 in the atmosphere (Peixoto and Oort 1992). This heat emitted from the atmosphere is dependent on the temperature of the atmosphere at the point of emission according to the Stefan–Boltzmann law. It is common practice to assume that temperature variations from equilibrium are small (the global mean temperature has fluctuated by only about 8°C during the last glaciation). Within a linear approximation, the emitted temperature will be proportional to the temperature difference from equilibrium (Ghil 1983). The boundary condition at the scale height of the atmosphere (which we take to be representative of the average elevation where radiation is emitted from the atmosphere) is then
We will use the value g = 1.7 W (m2 K)−1 as used by Ghil (1983) and Harvey and Schneider (1985).
Van Vliet et al. (1980) used Green’s functions to solve this model. The Green’s function of the Laplace-transformed diffusion equation is defined by
where G is governed by the same boundary conditions as Δ_T._ This equation can be solved by separating G into two parts: G a and G b with x < _x_′ and _x_ > _x_′, respectively, where G a and G b satisfy the homogeneous (no forcing) diffusion equation with a jump condition relating G a and G b:
The power spectrum of the average temperature in the atmosphere in terms of G is given by van Vliet et al. (1980) as
where _G_1 stands for the solution to the differential equation for G where the source point is located in the atmosphere. Here, Re denotes the real part of the complex expression. Two forms of G_1_a and G_1_b are necessary for x located above and below _x_′, respectively, due to the discontinuity in the derivative of _G_1 created by the delta function (the jump condition). The solutions of _G_1 that satisfy the above differential equation and boundary conditions are
and
where
L = (α/i ω)1/2, and _L_′ = (_α_′/i ω)1/2. Performing the integration van Vliet et al. obtained
For very low frequencies, or timescales that are large compared to the diffusion time across both the atmosphere and the ocean,
Reducing Eq. (29),
which is the low-frequency Lorentzian spectrum observed in the Vostok data. The crossover frequency as a function of the constants chosen for the model is
As mentioned before, in order for this model to be consistent with the Vostok data, a vertical eddy diffusivity of 0.1 m2 s−1 would be necessary. Such a value appears to be much too small to be realistic.
At higher frequencies,
and then
as observed. The high- and low-frequency spectra meet at
which also agrees well with that observed in the Vostok data [_f_1 ≈ (2 kyr)−1].
The timescales of _f_1 and _f_0 correspond to thermal and radiative equilibration of the coupled atmosphere–ocean system, respectively. At timescales of 20 kyr, the entire atmosphere and ocean are in thermal equilibrium. Besides the Vostok data, observational support for this thermal equilibration time is provided by comparisons of Antarctic and Greenland ice cores. Bender et al. (1994) found that temperature variations in Antarctica and Greenland are correlated above timescales of 2 kyr and uncorrelated below it. This suggests that 2 kyr represents the timescale of global thermal equilibrium.
Besides the frequency dependence of the power spectrum, the model we have presented predicts that the distribution of temperature variations from equilibrium obeys a Gaussian distribution since the stochastic term obeys a Gaussian distribution function and the local temperature fluctuation is related to that stochastic term through a linear transformation. By definition, the probability density function is defined only for timescales in which the temperature fluctuation time series are stationary. Gaussian time series with power-law power spectra of the form S(f) ∝ f_−_β are stationary if β < 1 and nonstationary if β ≥ 1 (Gomes da Silva and Turcotte 1994). Thus, a unique probability density function exists only for very long timescales (greater than 100 kyr) where the power spectrum is constant (β = 0) and for the range of timescales in which the power spectrum obeys S(f) ∝ f_−_β with β = 1/2. Matteucci (1990) has computed the probability distribution function for climatic variations at very long timescales with the SPECMAP stack. He obtained a Gaussian distribution. Similarly, Janosi and Vattay (1992) have obtained a Gaussian distribution with monthly temperature datasets of several decades length with the annual variability removed.
Manabe and Stouffer (1996) have completed power spectral analyses of variations in local atmospheric temperature in control runs of a coupled atmosphere–ocean–land surface model. They computed the power spectrum of temperature time series in each surface grid point and then averaged the power spectra at equal frequency values, as in our observational power spectral analysis. They found different spectra for continental and maritime grid points. Maritime grid points exhibited power-law power spectra from timescales of 1 month to several hundred years with an exponent of close to −0.25. Continental grid points, however, showed flat spectra up to timescales of about 100 yr, in contrast to observations. Exploring the similarities and differences between the approach in this paper, GCM results such as those of Manabe and Stouffer (1996), and observations should enable us to learn more about this fundamental problem in earth science.
Time series analysis of paleoclimatic data often exhibit a dominant peak near 100 kyr. Although variations in the eccentricity of the earth’s orbit occur with this frequency, this variation is not expected to produce a linear influence on climate change since this orbital variation results in only a fraction of a percent change in the amount of radiation incident on the earth (Kerr 1978). Although there are nonlinear models that predict a 100 kyr periodicity, it is generally agreed that the underlying mechanism for this peak is not well understood (Kerr 1978). The model presented in this paper leaves the question open as it does not predict any periodicity. The only component of the system thought to have a characteristic timescale of 100 kyr is the cryosphere (Mitchell 1976). Perhaps the cryosphere can produce a 100 kyr peak in the power spectrum when forced by the background spectrum predicted by the model of this paper. We believe that studies incorporating the cryosphere into our model are an important extension of our work that may lead to new thoughts on the nature of the 100 kyr periodicity.
We have applied the same model presented in this paper to variations in the solar luminosity from timescales of minutes to months (Pelletier 1996). A stochastic diffusion model of the turbulent heat transfer between the granulation layer of the sun (modeled as a homogeneous thin layer with a radiative boundary condition) and the rest of the convection zone (modeled as a homogeneous thick layer with thermal and diffusion constants appropriate to the lower convection zone) predicts the same spectral form observed in solar irradiance data recorded by the active cavity radiometer irradiance monitor (ACRIM) project and observed in the climate spectra reported here. The timescales of thermal and radiative equilibrium of the solar convection zone based upon thermal and diffusion constants estimated from mixing-length theory match those observed in the ACRIM data.
5. Conclusions
We have presented evidence that the power spectrum of atmospheric temperatures exhibits four scaling regions. We presented a model originally due to van Vliet et al. (1980) proposed to study temperature fluctuations in a metallic film (atmosphere) supported by a substrate (ocean) that matches the observed frequency dependence of the power spectrum of temperature fluctuations well. The timescale of radiative equilibrium predicted by the model with the most reliable estimates of thermal and diffusion properties of the atmosphere and ocean is much smaller than that observed in the Vostok core, however. The difference between the high-frequency spectra of continental and maritime stations may be interpreted in terms of the fact that air masses above maritime stations interact with the atmosphere above and the ocean below while air masses above continental stations interact with only the atmosphere above. The model assumes that atmospheric turbulent transport can be modeled as a stochastic diffusion process. Evidence supporting this hypothesis was provided by a correlation analysis of TOVS temperature measurements vertically in the atmosphere.
Fig. 1.
Power spectrum of atmospheric temperature variance observed and predicted by the model presented in this paper as a function of the frequency in yr−1. The crossover frequencies observed in the climatological record are, from left to right, _f_0 = 1 (40 000 yr), _f_1 = 1 (2000 yr), and _f_2 = 1 (1 month).
Citation: Journal of Climate 10, 6; 10.1175/1520-0442(1997)010<1331:AAMOTN>2.0.CO;2
Fig. 4.
Averaged power spectrum of (a) 90 maritime daily temperature and (b) 1000 continental daily temperature time series with the annual variability removed as a function of frequency in yr−1. The crossover frequency for the continental spectra is _f_2 = 1 (1 month)−1.
Citation: Journal of Climate 10, 6; 10.1175/1520-0442(1997)010<1331:AAMOTN>2.0.CO;2
Fig. 5.
Spatial cross-correlation function of (a) TOVS temperature data, (b) TOVS precipitable water vapor concentration, and (c) advection–diffusion model with u = 0.01 m s−1 and D = 100 m2 s−1 at time lags t = 1, 2, 3, 4, 5, and 6 days.
Citation: Journal of Climate 10, 6; 10.1175/1520-0442(1997)010<1331:AAMOTN>2.0.CO;2
Fig. 6.
(a) Geometry of the 1D, infinite space diffusion calculation detailed in the text; (b) boundary conditions appropriate for the air masses above the ocean (maritime stations), where the ocean acts as a thermal conductor; and (c) boundary conditions appropriate for the air masses above the continents (continental stations), where the continents act as a thermal insulator.
Citation: Journal of Climate 10, 6; 10.1175/1520-0442(1997)010<1331:AAMOTN>2.0.CO;2
