The complex dynamics of the seasonal component of USA’s surface temperature

The dynamics of the climate system has been investigated by analyzing the complex seasonal oscillation of monthly averaged temperatures recorded at 1167 stations covering the whole USA. We found the presence of an orbitclimate relationship on time scales remarkably shorter than the Milankovitch period related to the nutational forcing. The relationship manifests itself through occasional destabilization of the phase of the seasonal component due to the local changing of balance between direct insolation and the net energy received by the Earth. Quite surprisingly, we found that the local intermittent dynamics is modulated by a periodic component of about 18.6 yr due to the nutation of the Earth, which represents the main modulation of the Earth’s precession. The global effect in the last century results in a cumulative phase-shift of about 1.74 days towards earlier seasons, in agreement with the phase shift expected from the Earth’s precession. The climate dynamics of the seasonal cycle can be described through a nonlinear circle-map, indicating that the destabilization process can be associated to intermittent transitions from quasi-periodicity to chaos.


Introduction
The annual rotation of the Earth around the Sun provides a quasi-periodic solar forcing that continuously synchronizes the terrestrial climate.The resulting seasons observed on Earth represent the complex nonlinear response of atmosphere, land and oceans to this forcing and are one of the most important source of variability for the climate system.
Correspondence to: A. Vecchio (antonio.vecchio@fis.unical.it)In fact, even if climate changes are usually referred to trends in the average temperature records (Alley et al., 2003;Karl and Trenberth, 2003), many studies have shown that the analysis on the seasonal cycle of the Earth temperature could improve our knowledge about climate changes (Thomson, 1995(Thomson, , 1997;;Huybers and Curry, 2006;Pezzulli et al., 2005;Stine et al., 2009).Unlike the external solar forcing, which is almost constant from year to year, there is no guarantee that climatic seasons have to be the same each year.Actually seasonal variations cause more than 90% of the variance of a temperature record and represent one of the basic examples of the complex atmospheric response to external forcing (Thomson, 1995(Thomson, , 1997;;Mann and Park, 1996;Wallace and Osborn, 2002;Jones et al., 2003;Stine et al., 2009).
At a given latitude the seasonal cycle is generated by the balance between two main components: i) the direct insolation component, which varies with the tropical year (the time from equinox to equinox) τ tr 365.2422 days; ii) a mean energy transport effect related to the net radiation received by the Earth, which follows the anomalistic year (the time from perihelion to perihelion) τ an 365.2596 days (Thomson, 1995).The seasonal cycle can be characterized by two quantities, namely amplitude and phase.While changes in amplitude are commonly interpreted as caused by stochastic climate fluctuations (Huybers and Curry, 2006), the phase of the seasonal cycle, namely the synchronization of each season during time, is yet poorly investigated.The variability of the seasons is not something new in climate research.For example, Cook et al. (2000) pointed out that El Niño-Southern Oscillation (ENSO) has exhibited large changes in the amplitude and phase of the annual cycle as well as in the frequency/intensity of warm/cold events in the past century leading to inter-annual oscillations with multiple periods (Jin et al., 1994;Jiang et al., 1995).
A. Vecchio et al.: Season dynamics of Earth's surface temperature Previous studies about the phase of the global seasonality have underlined the presence of a phase shift.Thompson showed that, after 1940, the phase behavior started to change at an unprecedented rate with respect to the past 300 years (Thomson, 1995).Other authors indicated a global phaseshift towards earlier seasons (Mann and Park, 1996;Wallace and Osborn, 2002;Jones et al., 2003).A recent estimate (Stine et al., 2009) yields a global phase shift of the annual cycle of surface temperatures, over extratropical land between 1954 and 2007, of about 1.7 days.Even in this case the result seems to be highly anomalous with respect to the phase-shift of earlier periods.
The phase shift phenomenon has been attributed either to an increased temperature sensitivity to the anomalistic year forcing relative to the tropical year forcing due to Earth's precession (Thomson, 1995), or changes in albedo, soil moisture and short-wave forcing, also involved in changing modes of atmospheric circulation (Stine et al., 2009).The anomalous phase shift recorded after 1940 has been attributed to the increasing presence of CO 2 in atmosphere (Thomson, 1995;Stine et al., 2009).The phase shift is not predicted by the current Intergovernmental Panel on Climate Change models (Stine et al., 2009), thus representing a yet rather obscure physical effect of the climate system.The role and exact extent of natural and anthropogenic forcing for the climate evolution has been under much debate (Thomson, 1997;Keeling et al., 1996;Hasselmann, 1997;Rind, 2002;Alley et al., 2003;Karl and Trenberth, 2003;Scafetta and West, 2008), and any attempt to point out hidden aspects of climate dynamics is of considerable interest.

Data analysis
In the present paper we focused on the phase of the seasonal oscillation by using an ensemble of temperature time series T (t) from the United State Historical Climatology Network1 .The data set covers 110 yr from 1898 up to 2008 and is recorded by N=1167 different stations distributed on the whole USA.The seasonal component has been identified through the Empirical Mode Decomposition (EMD), a technique developed to process non-stationary data (Huang et al., 1998) and successfully applied in many different fields (Cummings et al., 2004;McDonald et al., 2007;Terradas et al., 2004;Vecchio et al., 2010).EMD decomposes the variance of each temperature record into a finite number of intrinsic mode functions (IMFs) and a residue, describing the mean trend, by using an adaptive basis derived from each data set (Huang et al., 1998), namely: (1) For a given time series, each IMF θ j (t) has its very own timescale and represents an oscillation with both amplitude and frequency modulations.This kind of decomposition is local, complete and in fact orthogonal (Huang et al., 1998;Cummings et al., 2004).The intrinsic mode functions yield instantaneous phases j (t) and, after a time derivative, the instantaneous frequencies ω j (t).For all records, the highest frequency IMF (j =0) represents the inter-seasonal stochastic component of the signal.High-order modes describe modulations of increasingly long periods, while the residue r m (t) represents the monotonically increasing local trend of temperature, commonly attributed to large scale warming since the urbanization contribution is smaller (see Peterson, 2003;Parker, 2006;Jones et al., 2008).The seasonal oscillation is described by the mode j =1 with a typical timescale of about τ 1 1 year.The statistical significance of information content for each IMF with respect to a white noise has been checked by applying the test developed by Wu and Huang (2004).
The 66% of the stations show an anomalous seasonal oscillation characterized by intermittent local decreases of the amplitude of the j =1 mode.For these stations the seasonal oscillation is spread over two EMD modes, namely the regular season is rediscovered when θ 1 (t) and θ 2 (t) are summed up.The remaining 34% of the stations shows a regular seasonal oscillation just isolated in the j =1 mode.An example from Evanston WY and Smithfield NC stations, characterized by regular and anomalous seasonal oscillation, respectively, is shown in Fig. 1 where the temperature records (panel a, b), the dynamics of θ 1 (t) (panels c, d), the instantaneous 1 (panels e, f) and unwrapped (panels g, h) phase, and the instantaneous frequency (panel i, j) are reported in a time interval of 30 yr.Each anomaly is associated with a strong variation of the instantaneous frequency and with an increasing phase-shift of the seasons, originating a destabilization of the phase 1 of seasonal oscillation.This phase shift is marked by steps in the unwrapped phase plot (panel h).All these features are not observed for the "regular" station (left column of Fig. 1).We have to remark that the amplitudes of the modes θ 1 are all above the 99% confidence level with respect to a white noise (Wu and Huang, 2004).
In Fig. 2 the time evolution of the EMD modes j =1 and j =2 (panel a, b) and their sum (panel c) for the Smithfield NC temperature record is shown.The time behavior of θ 1 (t) is typical of the mode mixing effect, which makes a single IMF consist of signals of widely separated scales (Huang et al., 1998).This feature of EMD is troublesome in signal processing where the main purpose is the signal cleanliness.In this framework, to effectively separate IMFs without mixed scales, the noise-assisted method named Ensemble EMD (Wu and Huang, 2009) has been developed.This approach consists in sifting an ensemble of white noise-added signal and treats the ensemble mean as the final true result.White noise series should cancel out in the averaging process, when a sufficient number of different realizations is  (a, b) refer to the rough data sets.The anomaly is underlined, in the right column, by the time evolution of: the IMF mode θ 1 (t) (d), the phase 1 (t) (f), the unwrapped phase (h) and the instantaneous frequency (j).These quantities have been calculated by using the first EMD mode j =1.
used, thus reducing the chance of mode mixing.By using the EEMD Wu et al. (2008) showed that the seasonal component of the analyzed temperature record is also spread over two modes.Our analysis indicates a correspondence between phase anomalies and spreading of seasonal oscillations over two modes.Hence it seems very interesting to investigate the role of the phase intermittency in the temperature records and to relate it to the physics of the system.To this purpose the EEMD approach does not seem suitable, since it represents an ad hoc mechanism to cancel the intermittency from IMFs.On the contrary, if the goal of the research is to study the intermittency phenomenon in the analyzed signals, then we need to analyze the IMFs as they are obtained by EMD.In this way an IMF, although affected by mode mixing, can be useful to identify when intermittency is present.Moreover, because of their orthogonality, IMFs, differently from EEMD modes, can be given some more direct physical interpretation.The EMD is a technique developed to be highly sensitive to the local frequency (or phase, being ω = dφ/dt) fluctuations.For records showing irregular seasonal oscillations the frequency is slightly different form the expected one during an anomaly.In these cases, the detection of two IMFs, necessary to describe the full contribution of the season, results from the properties of the EMD decomposition for which each mode is associated to a well defined time scale.If a given time scale is present only during a small portion of the signal, namely t * , the IMF describing this oscillation will be significantly different from zero only during t * .The mode j = 2 simply provides the value of the "anomalous" frequency and the time intervals in which it occurs.The meaningful quantity is the sum θ 1 + θ 2 describing the full contribution of the seasonal cycle to the temperature record.In this application, the usefulness of EMD resides in its ability to highlight the periods in which the frequency of the season is anomalous.We have to remark that the EMD represents a powerful tool to deseasonalize the temperature record (see e.g.Vecchio and Carbone, 2010) by subtracting the reconstructed seasonality θ 1 + θ 2 from the raw record.This kind of approach is more efficient than the classical deseasonalization procedured involving time averages, since the temperature records are far to be stationary.
The analysis of the statistical properties of phase-shift events represents an important tool to characterize the nature of this outstanding phenomenon.By considering the

sin fit Fourier
Method A 18.8 ± 0.4 20 ± 5 Method B 18.7 ± 0.2 18.5 ± 3.5 anomaly occurrence as a pointlike process, that is, each of them is supposed to occur at its starting time t i , we calculate the waiting time distribution (WTD), namely the probability density of the time intervals between two consecutive events P ( t).Quite interestingly the detected phase-shift events are strongly correlated in time.
Two independent procedures have been developed to recognize the time of occurrence of a season anomaly: 1. since an anomaly corresponds to a strong local variation of the instantaneous frequency of the seasonal IMF, it can be identified with the formation of a spike in the local frequency (see Fig. 1 panel j).Since in principle positive and negative excursions in frequency can occur, the time of occurrence of each season anomaly corresponds to the local maximum in the range where the absolute value of instantaneous frequency is greater than two standard deviations of its average 2. anomaly occurrence is detected from the j = 2 IMF characterized by a temporal behavior like those shown in Fig. 2 panel b.For this mode, the amplitude increases in correspondence of the season anomalies.The points of each interval where the absolute value of the amplitude exceeds two standard deviations of the chosen IMF are identified.For each interval the distance between extreme points, satisfying the previous threshold, defines the duration of the anomaly and identifies it.
The WTD calculated with the two different methods to identify is shown in Fig. 3.An exponential shape of WTD corresponds to Poisson processes.In our case roughly equispaced peaks, superimposed to the exponential decay, can be recognized, thus indicating that the WTD is not associated to a stochastic process but there is a dominant periodicity.This is clear by looking at the binned occurrence of seasonal anomalies as a function of time (Fig. 4): phaseshift events show an oscillating behavior.The main period P of the anomaly occurrence has been calculated through a sinusoidal fit over the red and black curve of Fig. 4, and by identifying the dominant peak in the Fourier periodogram (Fig. 5).The results are shown in Table 1.The uncertainties have been calculated from the fitting procedure and from the Fourier period resolution at the found peak.
The value of the period P is remarkably close to the 18.6year period of the Earth's nutation, which represents the Fig. 3. Waiting time distribution for the anomalies of seasonality for all the 1167 stations.Red and black curves refer to the different methods 1. and 2., described in the text, used to identify the anomalies.

Fig. 4.
Anomaly occurrence detected in our dataset.Red and black curves refer to the methods 1. and 2., described in the text, used to identify the anomalies.The period of modulation, calculated for both curves, is reported in Table 1.The blue curve refers to the cycle of the changing inclination of the Moon's orbit, with respect to the equatorial plane due to nutation, modeled by ψ=23 • 27 +5 • 09 sin( N t+π), where N =2π/18.6134yr (cf.Yndestad, 2006).stronger modulation to the Earth's precession.The same periodicity can be found in different geophysical phenomena connected with the Earth's nutation (Imbrie and Imbrie, 1980;Yndestad, 2006).Figure 4 shows that, even if weak differences (e.g. the double peaked maxima around 1955 and 1990), the general trend of the periodic modulation of the anomaly's occurrence is remarkably in phase with the inclination of the Moon's orbit with respect to the Earth's equatorial plane, all over the observation period, that is over at least six cycles.The linear Pearson's correlation coefficient Atmos.Chem.Phys., 10, 9657-9665, 2010 www.atmos-chem-phys.net/10/9657/2010/ between occurrence of seasonal anomalies (black curve in Fig. 4) and the inclination of the Moon's orbit (blue curve) assumes the value 0.57.Since the phase-shift anomalies cannot be considered as purely stochastic events, we are going to discuss the origin of the above periodicity.We think this is a strong evidence of nutational forcing on temperature records probably due to an influence of the anomalistic year variability on seasonal timing variability (see Thomson, 1995).The amplitude and phase variability of climate, which we call "season", is determined by the competing action of the direct insolation, having its maximum at the perihelion, and the net radiation received by the Earth, having its maximum at the summer solstice.In the Northern Hemisphere, summer solstice and perihelion are about 180 • out of phase.Consequently, slight perturbation on either component can destabilize the system by changing the resultant proportionality.We can conjecture that the motion of the Earth's axis due to the nutation, by affecting the insolation, can continuously perturb the climate system.Due to strong nonlinearities in the atmospheric system, the climate response to the annual cycle of the solar forcing can be surprisingly abrupt, for example as it happens in the case of the rapid onset of the Asian monsoon (Pezzulli et al., 2005).This coupling can generate impulsive destabilization of the phase which results in a global phase shift modulated by the nutation component.Being the temperature records strongly dependent on the local conditions, the nutation signal is detected only when a statistical analysis, over a significant number of stations, is performed.The above orbit-climate relation is amplified by the EMD, which is a technique extremely sensitive to the signal's phase shifts.
Three peaks, at low energy with respect to P , can be identified from the Fourier spectra of Fig. 5.They correspond to the periods P 1 = 13.9 ± 0.6 yr, P 2 = 10.1 ± 0.8 yr, P 3 = 8.5 ± 2.1 yr.Similar periodicities have been found in previous works.In detail, P 1 is consistent with the ∼ 15 yr periodicity in coastal surface air temperature in the Gulf of Alaska (GOA) (Wilson et al., 2007) attributed to large-scale coherent Pacific climate variability.P 2 could be related to the ∼ 11 yr periodicity in ice core sequences (Royer, 1993) attributed to solar cycle effects.P 3 might be attributed to changing tidal current speeds due to interannual variability of the lunar orbit, in particular to the 8.85 yr period of rotation of the lunar perigee around the Earth (McKinnell and Crawford, 2007).It must be remarked that a periodicity of about 7.8 yr has been also found in drought data (Cook et al., 1997).

A simple model that reproduces complex seasonality
To investigate the system response to the perturbation due to the nutation, the dynamics of the seasonal cycle has been described by the Phase Transition Curve (PTC), that is a basic tool to investigate recurrent dynamics as the time evolution of an oscillating system (e.g.Arnold, 1965;Glass and Mackey, 1979;Croisier et al., 2009;Glass, 2001, and references therein).The PTC is a map describing the dynamics of an angular variable α n , such that α n+1 =f (α n ) where f is a given function.In our case the variable α n is identified with 1 (t n ), namely the phase at a discrete time t n (n=0,1,...).A regular season is described by a linear map α n+1 =(α n +ω 0 )mod(2π ), where ω 0 is the intrinsic frequency of the cycle.In Fig. 6 empirical PTCs are reported for the two stations of Fig. 1 (Evanston WY in panel a and Smithfield NC in panel b).In both cases, the dynamics of the PTC locally deviate from the linear behavior.For both maps the points are spread around the straight line, corresponding to the linear PTC, because of weak stochastic fluctuations due to random vagaries of the weather at monthly scales.Moreover, large excursions away from the linear dynamics are observed in panel b) in correspondence of the anomalies.The graphic iteration of this PTC map around a detected anomaly, from 1972 to 1975, is shown in panel c).Nonlinearities of the climate system, that is the occurrence of fluctuations which occasionally destabilize the system, cause significant deviations from the linear dynamics.
The observed PTC, including the effect of anomalies, can be reproduced by a so-called circle map (Ott, 2002), that describes the dynamics of a system characterized by two competing forcing with different frequencies.According to the general theory of dynamical systems (Ott, 2002), the circle map, in presence of two forcing oscillations, can be expressed in a linear form where the winding number w corresponds to the ratio between the two competing frequencies.In our case w=τ tr /τ an is the ratio between tropic and anomalistic year.When w is irrational the dynamics of the map is quasi-periodic, namely the orbit obtained from the iteration densely fills the circle as time goes to infinity (Ott, 2002).A nonlinear coupling between the external periodic forces is described by adding a function f (φ n ) to the right hand side of Eq. ( 2).The most famous example is the sine-circle map investigated by Arnolds where k is a constant (Arnold, 1965;Ott, 2002).The presence of the nonlinear term destroyes the quasi-periodicity in favor of the frequency locking because of the coupling thus generating an interesting transition to a chaotic state.Roughly speaking, the set of periodic points, which is of zero measure for k=0, increases as k =0.For a fixed value of k, the rotation number as a function of w is an intricate sequence of periodic and quasi-periodic regions (Ott, 2002).Since in our case w changes because of the precession and nutation, the system moves into the net of periodic and quasi-periodic states and is continuously destabilized.The empirical PTC (panel b in Fig. 6) indicates that during the seasonal anomalies the regular dynamics of the seasonal oscillation is destabilized.This can be described by conjecturing that the frequency locking phenomenon depends on the nonlinear coupling with the atmosphere.In order to reproduce our result we modify the sine-circle map by adding a periodic perturbation of the winding number and a variable coupling parameter k, where the periodic term related to the nutation is R n =R 0 cos( N t n ) and the parameter k(R n ) is not kept constant and represents the response to this perturbation.In particular, we conjecture that the response of climate to the orbit perturbation is a threshold phenomenon, so that the behavior of k can be described by the following map where R th is a threshold value.When the inclination of the Earth's axis is greater than a critical value, the frequency locking of the seasonal cycle occurs abruptly.This can be reproduced for example by a simple on-off intermittent process (Platt et al., 1993).In this framework the parameter z n is suitably chosen to assume two values, say z n =1/2 or z n =4 with probability p and (1−p), respectively.When 1/3≤p≤0.47 the sequence of k n behaves as a typical on-off intermittency (Loreto et al., 1996).Results of the model Eqs.( 5) and ( 6), reported in Fig. 7, reproduce the observed behavior.

Conclusions
In this paper we have shown the results of EMD to analyze the temperature records of 1167 stations, covering 110 yr, over the whole USA.Our results confirm that a complete study of the phase of the seasonal component of temperature records is fundamental to investigate climate changes.These are not simply related to trends in temperature records since even phase-shifts of the seasonal component seem to play an important role.We identify intermittent periods in most temperature records causing local events of phase shift.We assume that the climate system, usually lying in a quasiperiodic state, can be occasionally destabilized, due to unbalances between direct insolation and radiation received by the Earth generated by the nutation, re-synchronizing itself in a relatively small time.In other words, an increased sensitivity to the anomalistic year variability influences the seasonal variability.The phase shift occurrence, in fact, is modulated by an oscillating component whose period and phase are remarkably close to the Earth's nutation.The occurrence of local phase shifts causes significant deviations of the system from the linear dynamics since the response of the system to this kind of perturbations is highly nonlinear.Our findings represent an indication that the phase of the season is influenced by the Earth's nutation.The global phase shift, underlined by different authors in the past (Thomson, 1995;Mann and Park, 1996;Wallace and Osborn, 2002;Jones et al., 2003;Stine et al., 2009), could be related to the global effect of the local dynamics of impulsive phaseshift events.EMD results indicate that each anomaly of the seasonal cycle of temperature corresponds to a local phase shift well identified by steps in the unwrapped phase plots.By making use of simple statistical arguments we investigate if the combined effect of the anomalies could result in a global phase shift of the temperatures as detected by some authors (Thomson, 1995;Mann and Park, 1996;Wallace and Osborn, 2002;Jones et al., 2003;Stine et al., 2009).This can be quantified by looking at the probability distribution P ( ) of =w/ω 0 (Fig. 8), where w=lim n→∞ [(φ n −φ 0 )/2π n] is the theoretical winding number.represents the phase shift, with respect to the initial phase, detected for each station The distribution is asymmetrical towards values greater than one, thus indicating an average shift towards earlier seasons.The maximum is found around the value 1.03 corresponding to a phase shift of about φ=|(1− )/ω 0 | 1.74 days with a standard deviation of about 0.08 days.This result is in close agreement with the global shift estimated by Stine et al. (2009).Thus the combined effect of the single phase shift on each station gives rise to a weak global phase shift of the seasonal cycle of surface temperatures of 1.74 days towards earlier seasons over 110 yr in agreement with past results.
We have to remark that the 18.6 yr periodicity has been found in other climatic time series like e.g, the air temperature and coastal sea surface temperature in the northeastern Pacific (Royer, 1993;McKinnell and Crawford, 2007), GOA coastal surface air temperatures (Wilson et al., 2007), drought occurrence (Currie, 1984;Cook et al., 1997), North Pacific Index and Pacific Decadal Oscillation Index (Yasuda et al., 2006) and has been commonly attributed to the external forcing of the lunar tides.Namely, the 18.6 yr lunisolar nodal cycle causes the long-term fluctuations of oceanic tides especially for the diurnal components.The high values of the tidal flow amplitude modulation, ∼ 20%, could affect the intermediate waters and large-scale oceanic circulation in the North Pacific.The coupling between bi-decadal variations in climatic data of northeast Pacific with the 18.6 yr nodal cycle has been modelled in terms of dissipation arising from local modulations of diapycnal mixing produced by the variable magnitude of diurnal tidal currents in the relatively shallow coastal ocean (Yasuda et al., 2006;McKinnell and Crawford, 2007).In other models, the bi-decadal fluctuation, interpreted by mid-latitude air-sea interactions with oceanic Rossby waves, has been obtained in a coupled ocean-atmosphere model and without the need for any external forcing (Yasuda et al., 2006).
It must be remarked that our model represents just a simple example to explain the observed behavior of the USA's temperature records, in terms of the variation of the insolation due to the nutation of the Earth.Previous papers, analyzing different climatic datasets, do not clearly indicate whether the bi-decadal periodicity is related to the variation of the insolation due to the nutation of the Earth and/or lunar tidal forcing.This interesting topic will be the object of future investigations.

Fig. 1 .
Fig. 1.Examples of EMD results from Evanston WY record (right column), characterized by regular seasonal oscillations, and from Smithfield NC (left column), showing season anomalies in a 30 yr time interval.(a,b) refer to the rough data sets.The anomaly is underlined, in the right column, by the time evolution of: the IMF mode θ 1 (t) (d), the phase 1 (t) (f), the unwrapped phase (h) and the instantaneous frequency (j).These quantities have been calculated by using the first EMD mode j =1.

Fig. 2 .
Fig. 2. Time evolution of the EMD modes j =1 and j =2 (a, b) and their sum (c) for the temperature record in panel (b) of Fig. 1.

Fig. 5 .
Fig.5.Fourier power spectra for the anomaly occurrence detected in our dataset.Red and black curves refer to the methods 1. and 2.

Fig. 6 .
Fig. 6.Phase Transition Curve α n+1 = f (α n ) calculated using the data from Evanston WY (a) and from Smithfield NC (b) stations.An example of iteration of the panel (b) map around a given seasonal anomaly is reported in panel (c).

Table 1 .
Main period P of the anomalies occurrence calculated through a sinusoidal fit (first column) and from the dominant peak in the Fourier periodogram (second column).