## Abstract

The universal role of the nonlinear one-third subharmonic resonance mechanism in generation of strong fluctuations in complex natural dynamical systems related to global climate is discussed using wavelet regression detrended data. The role of the oceanic Rossby waves in the year-scale global temperature fluctuations and the nonlinear resonance contribution to the El Niño phenomenon have been discussed in detail. The large fluctuations in the reconstructed temperature on millennial time scales (Antarctic ice core data for the past 400 000 years) are also shown to be dominated by the one-third subharmonic resonance, presumably related to the Earth's precession effect on the energy that the intertropical regions receive from the Sun. The effects of galactic turbulence on the temperature fluctuations are also discussed.

## 1. Introduction

Because global climatology is in a state of data accumulation, obtaining systematic knowledge from the data is not an easy and straightforward process. There is only one way to further our understanding of the physics of the phenomenon: accurate and punctual analysis of data. It is normal at this stage that the models are rather crude. Moreover, the global climate, as we know it now, is a nonlinear phenomenon. The simplest energy balance model for the space-averaged surface temperature *T* is a nonlinear equation,
1.1
where *C* is the heat capacity. The first term on the right-hand side of the balanced equation represents the incoming solar energy, with *q* as the solar constant and *α* as the albedo, and the second term represents the outgoing infrared energy, with *ε* and *σ* as the emissivity factor and the Stefan constant, respectively.

A more complex nonlinear model describing the interaction between the surface energy balance and the mass balance of the cryosphere has been suggested in Saltzman *et al.* [1] (in dimensionless excess variables),
1.2
and
1.3
where *s* is the sea-ice extent, *ζ*=*T*−*s* and *T* is the mean ocean surface temperature, *c*_{1} and *c*_{2} are certain constants, and stands for solar forcing. Nicolis [2] showed that the system equations (1.2) and (1.3) can be considered as perturbations (for 1≫*s*) of a reference system
1.4
and
1.5
which is a Hamiltonian system for the conservative Duffing oscillator.

The nonlinear Duffing oscillator also appears in the conceptual model of the tropical Pacific oscillations (for positive upper ocean temperature anomalies *T* driving shallow thermocline depth anomalies *h*) [3,4]
1.6
and
1.7
The periodic term in this equation represents the solar forcing (*η* is ambient noise, and *γ*, *b* and *F* are certain geophysical constants; see below for more details). The appearance of the nonlinear Duffing oscillator in the models of the global climate processes for different time–space scales—from palaeoclimate equations (1.2)–(1.5) to the decadal and inter-year oscillations equations (1.6) and (1.7)—shows that this nonlinear oscillator can be considered as a simple prototype of nonlinear climate systems (the nonlinear Duffing oscillator also appears naturally in the solar dynamo models [5]). Basic symmetries of the systems can play a crucial role in this phenomenon (see below). For the nonlinear systems far from equilibrium, the oscillatory behaviour can appear near the bifurcation points (the dynamics generated by the Duffing oscillator can exhibit both deterministic and chaotic behaviour, depending on the parameter range). One of the specific nonlinear properties of the Duffing oscillator is the so-called one-third subharmonic resonance [6]. Let us imagine a forced excitable system with a large amount of loosely coupled degrees of freedom schematically represented by Duffing oscillators with a wide range of natural frequencies *ω*_{0} (it is well known [7] that oscillations with a wide range of frequencies are supported by oceans and the atmosphere; see Brindley *et al.* [8]),
1.8
where denotes the temporal derivative of *x*, *β* is the strength of nonlinearity, and *F* and *ω* are characteristic of a driving force. It is known [6] that, when *ω*≈3*ω*_{0} and *β*≪1, equation (1.8) has a resonant solution
1.9
where the amplitude *a* and the phase *φ* are certain constants. This is so-called one-third subharmonic resonance with the driving frequency *ω*. For the system of oscillators under consideration, an effect of synchronization can take place and, as a consequence of this synchronization, the characteristic peaks in the spectra of partial oscillations coincide [9]. It can be useful to note, for the global climate modelling, that the odd-term subharmonic resonance is a consequence of the reflection symmetry of the natural nonlinear oscillators (invariance to the transformation ; see Minobe & Jin [10]).

This nonlinear phenomenon can be observed directly in real data. Let us start from the well-known (in relation to the global warming) graph representing the monthly global temperature data shown in figure 1 (the instrumental monthly data are available at http://lwf.ncdc.noaa.gov/oa/climate/research/anomalies/index.html; see Smith *et al.* [11]). One can see that the temperature is strongly fluctuating. Actually, the fluctuations are of the same order as the trend itself (figures 1 and 2). While the nature of the trend has been widely discussed, the nature of these strong fluctuations is still quite obscure. It will be shown below that the nonlinear one-third subharmonic resonance also produces strong fluctuations in the global temperature for much larger *palaeoclimate* time scales, as shown in figures 3 and 4 (the reconstructed data are available at http://www.ncdc.noaa.gov/paleo/metadata/noaa-icecore-6076.html; see Kawamura *et al*. [12]). Again, while the nature of the trend has been widely discussed (in relation to the glaciation cycles), the nature of these strong fluctuations is still quite obscure. We also know [5] that, in the Sun itself, Nature uses the nonlinear one-third *subharmonic* resonance to amplify the hydromagnetic dynamo effects to the observed 11 year periodic solar activity. And again, global rotation and nonlinear waves are the main components of the nonlinear mechanism.

Such universality of this nonlinear mechanism makes it a very interesting subject for a thoughtful investigation. The problem has a technical aspect too: detrending is a difficult task for such strong fluctuations. Most of the regression methods have linear responses. At the nonlinear non-parametric *wavelet* regression, one chooses a relatively small number of wavelet coefficients to represent the underlying regression function. A threshold method is used to keep or kill the wavelet coefficients. At the wavelet regression, the demands for smoothness of the function being estimated are relaxed considerably in comparison with the traditional methods. These advantages make the nonlinear wavelet regression method an adequate tool for detrending the strongly fluctuating natural data.

## 2. One-third subharmonic resonance and Rossby waves

Figure 1 shows (as the dashed line) the instrumental monthly global temperature data (land and ocean combined) for the 1880–2009 period, as presented at the NOAA website (http://lwf.ncdc.noaa.gov/oa/climate/research/anomalies/index.html; see Smith *et al.* [11]). The solid curve (trend) in the figure corresponds to a wavelet (symmlet) regression of the data (see Scafetta & West [13]). Figure 2 shows the corresponding detrended fluctuations, which produce a statistically stationary set of data. We use a nonlinear non-parametric wavelet regression with a relatively small number of wavelet coefficients representing the underlying regression function. A threshold method is used to keep or kill the wavelet coefficients. In this case, in particular, the universal (VisuShrink) thresholding rule with a soft thresholding function was used. Figure 5 shows a spectrum of the wavelet regression detrended data calculated using the maximum entropy method (because it provides the optimal spectral resolution even for small datasets). One can see in this figure a small peak corresponding to a 1 year period and a huge, well-defined peak corresponding to a 3 year period. The fluctuations in oceanic temperature cause certain variations in the sea-surface height. These variations are intermixed with the sea-surface height variations caused by the oceanic planetary Rossby waves. The oceanic planetary Rossby waves play an important role in the response of the global ocean to the forcing [14], and they are of fundamental importance to ocean circulation on a wide range of time scales (it was also suggested that the Rossby waves play a crucial role in the initiation and termination of the El Niño phenomenon; see also below). Therefore, they present a favourable physical background for the global subharmonic resonance. It should be noted that the variable *h* in equation (1.7) represents western boundary Rossby wave reflection of its counterpart in equation (1.6). This provides a negative feedback to the *T* tendency in the eastern equatorial ocean via a nonlinear thermocline displacement: *h*+*bh*^{3} (its vertical gradient separating near-surface and deep-ocean mixed layers [3,4]). The solar periodic forcing plays an analogous role in this equation. For that reason, it is interesting to also look separately at global sea-surface temperature anomalies. These monthly data for the 1850–2008 period are available at http://jisao.washington.edu/data/globalsstanomts/ (see Smith *et al.* [15]). Figure 6 shows a spectrum of the wavelet regression detrended data. The spectrum seems to be very similar to the spectrum presented in figure 5. This is the one-third subharmonic resonance with the driving frequency *ω* corresponding to the *annual* North–South asymmetry of the solar forcing (the huge peak in figure 5 corresponds to the first term in the right-hand side of equation (1.9)).

Because the high-frequency part of the spectrum is corrupted by strong fluctuations (the Nyquist frequency equals 0.5 [month^{−1}]), it is interesting to look at the corresponding autocorrelation function *C*(*τ*) in order to understand what happens on the monthly scales. It should be noted that scaling of the defect of the autocorrelation function can be related to scaling of the corresponding spectrum,
2.1
Figure 7 shows the defect in the autocorrelation function on ln–ln scales in order to estimate the scaling exponent *α*≃0.6±0.04 (the straight line in this figure indicates the scaling equation (2.1)). The existence of oceanic Rossby waves was confirmed recently by NASA/CNES TOPEX/Poseidon satellite altimetry measurements. Corresponding to these measurements, a spectrum of the sea-surface height fluctuations, calculated in Zang & Wunsch [14] with *daily* resolution, is shown in figure 8. The dashed straight line in this figure is drawn in order to indicate correspondence to the scaling equation (2.1): 1+*α*≃1.6. The Rossby waves (together with the Kelvin waves) and a strong atmosphere–ocean feedback provide physical background for the El Niño phenomenon (see [3,4] and references therein). Figure 9 shows the spectrum for the wavelet detrended fluctuations of the so-called Global-SST ENSO index, which captures the low-frequency part of the El Niño phenomenon (the monthly data are available at http://jisao.washington.edu/data/globalsstenso/). The annual forcing can come from the oceanic Rossby waves (figure 8). To support this relationship, we show in figure 10 the defect in the autocorrelation function calculated using the wavelet detrended fluctuations of the Global-SST ENSO index. The ln–ln scales have been used in figure 10 in order to estimate the scaling exponent *α*≃0.6±0.03 (the straight line in this figure indicates the scaling equation (2.1)). Using these observations, one can suggest that the El Niño phenomenon has the one-third subharmonic resonance as a background.

## 3. Nonlinear palaeoclimate

Recent palaeoclimate reconstructions provide indications of nonlinear properties of the Earth's climate at the late Pleistocene [16] (the period from 0.8 Myr to present). A long-term decrease in atmospheric CO_{2}, which could result in a change in the internal response of the global carbon cycle to the obliquity forcing, has been mentioned as one of the principal reasons for this phenomenon [17]. At the present time, one can recognize at least two problems in the nonlinear palaeoclimate, which we will address in the present paper using recent data and speculations [18]. (i) Reconstructed air temperatures on millennial time scales are known to be strongly fluctuating (e.g. figures 3 and 4). While the nature of the trend has been widely discussed (in relation to the glaciation cycles), the nature of these strong fluctuations is still quite obscure. The spectral analysis of the wavelet regression detrended data reveals a rather surprising nature of the strong temperature fluctuations. Namely, the detrended fluctuations of the reconstructed temperature are completely dominated by the one-third subharmonic resonance, presumably related to the Earth's precession effect on equatorial insolation. (ii) The influence of galactic turbulent processes on the Earth's climate can be very significant for time scales less than 2.5 kyr.

Figure 3 shows reconstructed air temperature data (dashed line) for the 0–340 kyr period (as presented at http://www.ncdc.noaa.gov/paleo/metadata/noaaicecore-6076.html; Antarctic ice cores data; see also Kawamura *et al.* [12]). The solid curve (trend) in the figure corresponds to a wavelet (symmlet) regression of the data. Figure 4 shows corresponding detrended fluctuations, which produce a statistically stationary set of data. Figure 11 shows a spectrum of the wavelet regression detrended data. One can see in this figure a small peak corresponding to the period approximately 5 kyr and a huge well-defined peak corresponding to the period approximately 15 kyr. We also obtained analogous results (approx. 10% larger) from the ‘Vostok’ ice core data for the 0–420 kyr period (for a description of the data, see Petit *et al.* [19]).

The origin of the periodic energy input with the period approximately 5 kyr can be related to the dynamics of the energy that the intertropical regions receive from the Sun (equatorial insolation). Indeed, Berger *et al.* [20] found that a clear and significant 5 kyr period is present in this dynamic over the last 1 Myr. The amplitude of the 5 kyr cycle in the insolation decreases rapidly the further from the Equator. Using the fact that the double insolation maximum and minimum arise in the tropical regions in the course of 1 year, Berger *et al.* [20] speculated that this period in seasonal amplitude of equatorial insolation is determined by the fourth harmonic of the Earth's precession cycle. It should be noted that the idea of a significant role for the tropics in generating long-term climatic variations is rather a new one (see Berger *et al.* [20] for relevant references). In Hagelberg *et al.* [21], for instance, the authors speculated that the high-frequency climate variability (on millennial time scales) could be related to high sensitivity of the tropics to summer time insolation. Then, the oceanic advective transport could transmit an amplified response of tropical precipitation and temperature to high latitudes. The physical mechanism of this amplification is still not clear and the one-third subharmonic resonance discussed earlier could be a plausible possibility (in this respect, it is significant that we used the Antarctic data).

## 4. Galactic turbulence and the temperature fluctuations

Because the high-frequency part of the spectrum is corrupted by strong fluctuations (the Nyquist frequency equals 0.5 [(500 yr)^{−1}]), it is interesting to look at the corresponding autocorrelation function *C*(*τ*) and at the structure functions *S*_{p}(*τ*) (of different orders *p*) in order to understand what happens on the millennial time scales. The correlation function defect 1−*C*(*τ*) is proportional to the second-order structure function *S*_{2}(*τ*). Therefore, we can compare results obtained by these different tools. First, let us look at the autocorrelation function *C*(*τ*). Figure 12 shows a relatively small-time part of the correlation function defect. The ln–ln scales have been used in this figure in order to show a power law (the straight line) for the second-order structure function: *S*_{2}(*τ*)=〈|*x*(*t*+*τ*)−*x*(*t*)|^{2}〉,
4.1
This power law, ‘’, for the structure function (by virtue of the Taylor hypothesis transforming the time scaling into the space one [22]) is known for fully developed turbulence as Kolmogorov's power law.

Although the scaling interval is short, the value of the exponent is rather intriguing. This exponent is well known in the theory of fluid (plasma) turbulence and corresponds to the so-called Kolmogorov's cascade process. This process is very universal for turbulent fluids and plasmas [23,24]. In spite of the fact that a magnetic field is presumably important for interstellar turbulence, the Kolmogorov description can still be theoretically acceptable even in this area [25]. Moreover, Kolmogorov-type spectra were observed on scales up to kpc. In order to support Kolmogorov turbulence as the background of wavelet regression detrended temperature modulation, we also calculated structure functions *S*_{p}(*τ*)=〈|*x*(*t*+*τ*)−*x*(*t*)|^{p}〉 with different orders of *p*. In the classic Kolmogorov turbulence (at very large values of Reynolds number [26]),
4.2
(at least for *p*≤3). Figure 13 shows (as circles) the scaling exponent *ζ*_{p} against *p* for the same data as were used for the calculation of the autocorrelation function (figure 12). The bars show the statistical errors. The straight line in figure 13 corresponds to the strictly Kolmogorov turbulence with *ζ*_{p}=*p*/3 (equation (4.2)). One can see good agreement with the Kolmogorov turbulence modulation.

For turbulent processes on Earth and in the heliosphere, Kolmogorov's scaling with such large time scales certainly cannot exist. Therefore, one should think about a galactic origin of Kolmogorov turbulence (or turbulence-like processes [27]) with such large time scales. Let us recall that the diameter of the Galaxy is approximately 100 000 light years. This is not surprising if we recall the possible role of the galactic cosmic rays for the Earth's climate. Galactic cosmic ray intensity at the Earth's orbit is modulated by galactic turbulence [22]. On the other hand, the galactic cosmic rays can determine the amount of cloud cover (a very significant climate factor) on global scales through massive aerosol formation [28,29,30,31]. Thus, the galactic turbulence can modulate the global temperature fluctuations by the Kolmogorov scaling properties on millennial time scales. If one knows the characteristic velocity scale *v* for the Taylor hypothesis, one can estimate the outer space scale of the corresponding galactic turbulence as *L*≥2500 yr×*v*. However, it is not clear what estimate we should take for *v*. For instance, one could try the velocity of the solar system relative to the cosmic microwave background rest frame: *v*∼370 km s^{−1}. In this case, one obtains *L*∼1 pc. It should be noted that, in a recent paper [30], it is suggested that the typical outer scale for spiral arms can be as small as 1 pc, and in interarm regions the outer scale can be larger than 100 pc. Because the solar system and Earth are at the present time within the Orion Arm, this suggestion is in agreement with the above estimate. Although the suggestion made by Haverkorn *et al.* [30] is still under active discussion, the palaeoclimate consequences of this suggestion could be very interesting and we will discuss one of them here. Namely, when orbiting the galactic centre, the Solar System and Earth are in the interarm regions; the reverse of the Taylor hypothesis provides us with the outer time scale of approximately 2500×100 years. This time scale is larger than any known glaciation period (which are determined by the periods related to the Earth orbiting around Sun; see [32]). That strong *turbulent* fluctuations of the cosmic rays flux on such large time scales should prevent glaciation cycles from occurring when the Solar System is in the interarm regions. Deglaciation of the Earth related to the interarm regions was suggested in Shaviv [29] and explained by the difference in intensity of the cosmic ray flux in the spiral arms and in the interarm regions. It is difficult to estimate at the present time which of the two mechanisms is more efficient. Anyway, they are working to the same end and both are open to discussion.

## 5. Conclusion

In the complex natural systems under consideration, the fluctuations are often of the same order as those of the trend itself. This phenomenon is related to a nonlinear character of the fluctuations. Does this nonlinear phenomenon have a universal nature? The main three points of the above consideration are:

— the known (crude) models can be reduced to the Duffing oscillator;

— the one-third subharmonic resonance is a generic property of this oscillator; and

— the data indicate a crucial role of this resonance in the generation of large-scale fluctuation.

The next generation of models of the global climate should take into account the highly nonlinear nature of the Sun–climate interaction, which is not realized directly through the trends but through the large-scale fluctuations.

At this stage of the data accumulation and analysis, the limitations of these conclusions are unknown. However, we believe that the results considered above could attract the attention of nonlinear scientists to the natural systems because, despite their overwhelming complexity, these systems are governed by simple physical laws.

## Acknowledgements

I thank S. I. Abarzhi and K. R. Sreenivasan for suggestions, criticism and encouragement. The data were provided by the National Climatic Data Center at NOAA, by the Joint Institute for the Study of the Atmosphere and Ocean, and by the SIDC team, World Data Center for the Sunspot Index, Royal Observatory of Belgium. I also acknowledge that a software tool provided by K. Yoshioka was used for the computations.

## Footnotes

One contribution of 13 to a Theme Issue ‘Turbulent mixing and beyond: non-equilibrium processes from atomistic to astrophysical scales I’.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.