On the nonstationarity of the decadal periodicities of the length of day

The Earth’s rotation rate is not constant, but changes on all observable timescales, from subdaily to decadal and longer. These variations are usually discussed in terms of variations in the length of the day (LoD) and are caused by processes acting within the interior, at the surface and outside of the Earth. Here, we investigate the presence of long-standing decadal variations in yearly LoD data covering the period from 1832 to 2009 by applying the Hilbert– Huang transform (HHT). The HHT has been slightly modified here to take into account the uncertainty of LoD values that has changed greatly in time due to the use of different LoD measurement techniques. The LoD time series has been completely decomposed into five intrinsic mode functions (IMF) and a residual trend. The estimation of instantaneous frequencies and related amplitudes of the obtained IMFs has allowed us to compute the Hilbert spectrum that has been used as the starting point for studying and discussing the stationarity of typical LoD timescale stationarity. The obtained results while showing the presence of multiple periodicities also indicate the absence of really stationary periodicities. Therefore, rather than considering the processes taking place in the Earth’s core as the result of a superposition of oscillations (i.e. stationary mechanisms) occurring on a discrete number of different timescales, it would be better to think of a superposition of fluctuations that are intermittent in both frequency and amplitude.


Introduction
The question of whether or not the rotation rate of the Earth is constant was investigated at least as early as the seventeenth century.However, it was only in the twentieth century that the observation of planetary motions confirmed its irregular behaviour.Nowadays, it is well known that the Earth's rate of rotation ( ), and hence the length of the day (LoD), is not constant in time, but exhibits changes of the order of a few parts in 10 8 , i.e. δ / ∼ 10 −8 .Though variations with the largest amplitude occur over the decadal timescale, observations of the evolving state of Earth rotation have revealed the occurrence of variability on many timescales ranging from days to centuries and longer.This wide range of timescales is the expression of the wide variety of processes influencing Earth rotation.These processes involve external tidal forces, superficial processes concerning the atmosphere, oceans, and hydrosphere, and internal processes taking place at the core-mantle boundary as well as within the solid Earth itself (Hide and Dickey, 1991).This is the reason why, for many years, scientists in the field of astronomy and Earth sciences have been focusing their attention on the interpretation of the Earth rotation fluctuations as a key to understanding fundamental terrestrial dynamical processes.
LoD variations occurring over the longest timescales are thought to be primarily a consequence of the tidal friction resulting from the gravitational attraction of the Sun and Moon on the rotating Earth.These interactions produce a secular slowing of rotation and consequently a secular increase in the LoD.The linear increase in the LoD is expected to be about 2.3 ms century −1 , but what is actually observed from the ancient astronomical record is an increase of about 1.7 ms century −1 .The difference between these two values, which is sometimes called the anomalous rate of increase in LoD, is believed to be caused by a linearly varying redistribution of mass within the Earth associated with the so-called post-glacial rebound.Indeed, although the ice began to melt about 18 000 yr ago and was mostly gone 5000 yr ago, the depressions of the Earth's surface have not completely disappeared yet.
Variations in the LoD occurring on relatively short timescales, from subannual to annual, characterized by an amplitude up to 2 ms have been explained primarily by the exchange of angular momentum between the solid Earth and its atmosphere and oceans, as deduced from many studies (e.g.Hide et al., 1980;Marcus et al., 1998;Gross et al., 2004).
Between these two kinds of variations are those occurring over timescales of several decades with an amplitude of a few milliseconds that cannot be easily explained by surface processes.These variations are generally believed to be caused by the transfer of rotational angular momentum between the Earth's liquid metallic core and the overlying solid mantle (Wahr, 1988;Jault, 2003).When the mantle gains angular momentum, its rotation rate increases, while the observed LoD decreases.Various mechanisms have been invoked to explain how the angular momentum is transferred across the core-mantle boundary.Commonly cited mechanisms involve viscous (Rochester, 1984), gravitational (Jault and Le Mouël, 1990;Buffet, 1996), electromagnetic (Holme, 1998) and topographic (Hide, 1969(Hide, , 1977) ) core-mantle coupling, but arguments for and against each of these mechanisms can be found as a consequence of the poor knowledge of the physical conditions in the core.
Since the Earth's magnetic field is produced by fluid flows occurring in the Earth's core, whatever the nature of the torque, the assumption that the decadal fluctuations are due to these flows suggests that there ought to be some correlation between the time variability of observed LoD and that of the observed Earth's magnetic field.Unfortunately, studying correlation is complicated because of the attenuation suffered by magnetic field variations travelling upward through the conducting mantle.When correlations have been identified, the time lag between changes in the LoD and similar changes in the observed magnetic field has been used to help constrain mantle electrical conductivity.
One of the keys to explaining the possible correlation between the observed LoD and Earth's magnetic field could be found in magneto-hydrodynamic torsional waves, i.e. in the rotation of core fluid on concentric cylinders coaxial with the Earth's rotation axis.It has been suggested that these torsional oscillations, with proper period and spatial structure, may transfer sufficient angular momentum across the coremantle boundary to produce the observed LoD variations on decadal timescales (Jault et al., 1988;Jackson et al., 1993).
Periodicities of torsional oscillations have not been assessed precisely yet.The reason for the large uncertainty in these periodicities arises from the poorly known physical conditions of the core and from the mathematical models constructed to describe the flow at the top of the core.Indeed, the periods of the fundamental modes of the torsional oscillations depend crucially on the assumed magnitude of the component of the magnetic field orthogonal to and away from the rotation axis (B s ), which is not well known and, partially, on the viscosity of the fluid core.For instance, according to Braginsky (1984), torsional oscillations can have a period of ∼ 60 yr if the amplitude of B s is around 0.2 mT, but if it is closer to 1-2 mT, the period of the free modes shifts to smaller values.Moreover, Gillet et al. (2010) have recently hypothesized the existence of torsional oscillations recurring in the core interior every 6 yr.These torsional oscillations, whose angular momentum accounts well for the change in the length of the day with a six-year period detected over the second half of twentieth century, seem to be carried by an internal field (B s ) of about 4 mT.
Besides Braginsky (1984) and Gillet et al. (2010), many other authors have estimated the characteristic periodicities of torsional oscillations and of the typical periodicities of LoD decadal variations.Currie (1973) analysed LoD data primarily from the 18th and 19th centuries using the maximum entropy method and found a periodicity of 57.5 yr.The same method was used by Jin and Thomas (1977), who found well-defined signals at 66 and 33 yr.Some years later, Zatman and Bloxham (1997) showed that the fluid flow at the surface of the core is consistent with the presence of two torsional oscillations travelling inside the core with periods of 76 and 53 yr.Hide et al. (2000) identified a dominant variability period of approximately 65 yr in LoD data.More recently, Roberts et al. (2007) confirmed the existence of a periodicity of around 60 yr by applying the empirical mode decomposition (EMD) analysis to LoD data and to the Earth's magnetic field by suggesting the existence of a correlation between the two signals with a well-determined lag.
Considering the nonlinearity and non-stationarity of LoD time series, we applied the complete Hilbert-Huang transform (HHT) here to try to gain a deeper knowledge of LoD typical periodicities.So, HHT was performed by decomposing the original signal into a series of monocomponents by means of EMD and then evaluating the associated Hilbert spectrum to study the properties of the original signal in the time-frequency-energy domain.Additionally, we have developed a procedure to estimate a Hilbert spectrum able to take into account the information on LoD actual measurement errors.In this way, the estimated LoD characteristic periodicities, as shown in the following, acquire a more reliable statistical meaning.Using the Hilbert spectrum it is possible to both study the temporal evolution of the instantaneous frequencies characterizing the signal and evaluate the actual presence and stationarity of the decadal periodicities, such as that of about 60 yr that have been found in the excess LoD data by other authors.
The paper is organized as follows: in Sect.2, we briefly describe the technique of the Hilbert-Huang transform introduced by Huang et al. (1998), in Sect. 3 we explain how HHT has been applied to test the existence and stationarity of the excess LoD decadal periodicities and, finally, in Sect. 4 we discuss and summarize our results.
2 The Hilbert-Huang transform: a brief introduction HHT provides a new method for analysing nonlinear and non-stationary time series, allowing the exploration of intermittent and amplitude-varying processes.HHT consists of the successive use of two mathematical techniques: EMD and Hilbert spectral analysis (Huang et al., 1998(Huang et al., , 2003;;Huang and Wu, 2008).The central idea of HHT is that a time series can be decomposed into simple oscillatory modes of significantly different average frequencies and a residue.These modes are called intrinsic mode functions (IMFs) and are directly obtained from the data with no a priori assumptions regarding their nature.IMFs are constructed to satisfy two specific conditions: 1. the number of extrema and of zero crossings must be either equal or differ at most by one; 2. the mean value of the envelope defined by the local maxima and of the envelope defined by the local minima is zero.
We will not linger over details of the complete procedure to perform EMD since it is reported in many scientific papers, for instance in Huang et al. (1998Huang et al. ( , 2003)), Huang and Wu (2008) and Flandrin et al. (2004).Here, we just emphasise that for the stoppage condition of the sifting process we use, as suggested by Huang et al. (1998), the Cauchy-type convergence criterion where the threshold value limiting the size of the standard deviation computed from two consecutive siftings has been set so as to have a diadic decomposition in the case of fractional Browian motion time series having the same length of actual time series (Flandrin et al., 2004).
After decomposing a signal into its monocomponents, i.e. its IMFs, it is possible to calculate the instantaneous frequencies and the associated energies using, for example, the traditional Hilbert transform.Here, however, we computed the Hilbert spectrum by estimating, for each IMF, the instantaneous frequency and associated energy through the so-called direct quadrature method (Huang et al., 2009).The estimation of instantaneous frequency from experimental data is not a trivial mathematical problem and many algorithms have been proposed.For instance, Huang et al. (2009) tested different methods for instantaneous frequency estimation, finding that the normalized Hilbert transform and direct quadrature gave the best results.One of the advantages of direct quadrature is that it is not affected by the occurrence of negative frequencies, a problem generally suffered by a traditional Hilbert transform.
The method of direct quadrature is based on the principle that a monocomponent signal, say the IMF c i (t) (where i = 1, . . ., k with k the number of IMFs into which the signal x(t) has been decomposed), can be written as the product of its envelope A i (t) and its carrier cos φ i (t) as follows: where φ i (t) is the phase function, and A i (t) and cos φ i (t) are the amplitude-modulated (AM) and the frequency-modulated (FM) parts of the signal, respectively.Thus, the instantaneous frequency f i (t) is given directly by differentiating the phase with respect to time, f i (t) = dφ i (t)/dt.The separation of a given IMF into its AM and FM parts has been achieved empirically following the iterative normalization scheme proposed by Huang et al. (2009).This scheme works as follows.Once all IMFs are found, all local maxima of the absolute value of the ith IMF, i.e. | c i (t) |, are connected in the first iteration by means of a cubic spline curve, say e i,1 (t).Then, the functions y i,1 (t), . . ., y i,n (t) are estimated via an iterative procedure as follows: where e i,k (t) is the cubic spline curve connecting all local maxima of the absolute value of the y i,(k−1) with k = 1, . . ., n. Normalization stops at iteration n when all the values of the function y i,n (t) are less than or equal to unity.
Having removed the amplitude modulation, y i,n (t) represents the FM part of the IMF c i (t) and, according to Eq. ( 1), y i,n (t) = cos φ i (t).Based on simple trigonometric relations, the instantaneous frequency for the ith IMF c i (t) is given by: while, according to Eq. ( 1), the AM part is given by: where A i (t) provides the instantaneous amplitude of the ith IMF c i (t) and reflects how the energy (defined as the amplitude squared), associated with the instantaneous frequency, changes with time.The Hilbert spectrum is obtained by plotting on the same graph the curves for the instantaneous frequency of each IMF as a function of time and the associated amplitude.The corresponding Hilbert spectrum H (f, t) (defined in terms of amplitudes, , is designed to represent the amplitude/energy in a time-frequency (or alternatively timeperiod) representation.

Data and analysis
In the present study we focused on a time series of excess length of day, LoD spanning May 1832-May 2009 at yearly intervals.In detail, we extended the smoothed length of day series LUNAR97, covering the period May 1832-May 1997 at yearly intervals (Gross, 2001), to the present using COMB2010 series.This consists of daily values Time [year] alues and uncertainties for the excess length of day, ∆LoD, spanning 1832.5-2009.5 at 1-year interhe Hilbert spectrum is obtained by plotting on the same graph the curves for the instantarequency of each IMF as a function of time and the associated amplitude.The corresponding spectrum H(f,t) (defined in terms of amplitudes, H(f,t) = A(f,t), or squared amplitudes, = A 2 (f,t)) is designed to represent the amplitude/energy in a time-frequency (or altername-period) representation.
a and Analysis resent study we focused on a time series of excess length of day, ∆LoD spanning 1832.5atyearly intervals.In detail, we extended the smoothed length of day series LUNAR97, g the period 1832.5-1997.5 at yearly intervals (Gross, 2001), to the present using COMB2010 his consists of daily values and uncertainties for the length of day from January 20, 1962 to , 2010 and can be downloaded from ftp://euler.jpl.nasa.gov/keof/combinations/2010/.So, it COMB2010 that yearly values of length of day and relative uncertainty for the last 12 years 009) have been evaluated.We notice that the LoD series of Gross (2001)  and uncertainties for the length of day from 20 January 1962 to 28 May 2010 and can be downloaded from ftp://euler.jpl.nasa.gov/keof/combinations/2010/.So, it is from COMB2010 that yearly values of length of day and relative uncertainty for the last 12 yr (1998-2009) have been evaluated.We notice that the LoD series of Gross (2001) has been obtained after many transformations.The LoD series is indeed derived from a Kalman Earth orientation filter based on a combination of independent Earth rotation measurements utilizing the techniques of optical astrometry, very long baseline interferometry (VLBI) and lunar laser ranging (LLR).The result is a smoothed and interpolated estimate of the length of day.This time series, as obtained, does not permit the analysis of periodicities of a few years because Gross (2001) has applied Gaussian filters to the original data of Jordi et al. (1994).Nevertheless, longer periodicities of the length of day, as decadal ones, can be appropriately investigated using this time series.
Figure 1 shows the LoD yearly values, i.e. values of the excess LoD in milliseconds with respect to the standard day of 86 400 s, for a time interval spanning 178 yr.Each LoD value is plotted together with its associated uncertainty (Gross, 2001) Since this time series is non-stationary and nonlinear, we studied it using HHT.For this reason, we pre-processed data using EMD and then estimated its Hilbert spectrum by direct quadrature as introduced by Huang et al. (1998Huang et al. ( , 2003) ) and Huang and Wu (2008).We belive that HHT allows us to better define values and stationarity of typical LoD periodicities and to gain a deeper insight into the underlying processes that are expected to influence Earth's rotation rate.So, we first applied EMD in the standard way without considering the effects of different error measurements and using as the stoppage criterion that discussed in Sect. 2.  value is plotted together with its associated uncertanity (Gross, 2001) Since this time series is non-stationary and nonlinear, we studied it using HHT.For we pre-processing data using EMD and then estimated its Hilbert spectrum by direct qu introduced by Huang et al. (1998Huang et al. ( , 2003) ) and Huang & Wu (2008).We belive that HHT   EMD yielded five prominent IMFs and a residue, as shown in Fig. 2. Figure 3 exhibits a comparison between the power spectral densities (PSDs) of the excess LoD ( LoD) and of the single IMFs into which the original signal has been decomposed.The sum of the IMF PSDs correctly reproduces the PSD of the LoD, confirming the completeness of our decomposition.Furthermore, we may notice that the shape of the LoD PSD is well described by a stretched Lorentzian shape, with α ∼ 2.8 and f 0 = [0.0137±0.0008]yr −1 .This suggests that the spectrum of LoD fluctuations is mainly broadband, thus supporting the presence of nonlinearities and also nonstationary features.We will return to this point in the next section.on among the power spectral densities (PSDs) of the excess LoD (∆LoD) and of the single which the original signal has been decomposed.The sum of the IMF PSDs correctly s the PSD of the ∆LoD confirming the completeness of our decomposition.Furthermore, otice that the shape of the ∆LoD PSD is well described by a stretched Lorentzian shape, 2.8 and f 0 = [0.0137±0.0008]yr −1 .This suggests that the spectrum of ∆LoD fluctuations broadband thus supporting the presence of nonlinearities and also non-stationary features.
eturn to this point in the next section.Figure 4 shows the average frequencies fi of the IMFs c i , plotted in Fig. 2, as a function of the IMF index i.Each average frequency has been computed by means of the Fourier PSD S i (f ) associated with each IMF, as follows: The error associated with each value of fi , as shown in Fig. 4, corresponds to the standard deviation, i.e. the root mean square of the second central moment that, similarly to the case of Eq. ( 6), has been obtained by applying the following relation: The observed scaling of the average IMF frequencies as a function of the IMF index is the classical where i is the IMF index and f c is a characteristic frequency.This means that the empirical mode decomposition acts as 0.01 The error associated with each value of fi , as shown in Figure 4, corresponds the standa i.e. the root mean square of the second central moment that, similarly to the case of Eq obtained by applying the following relation: The observed scaling of the average IMF frequencies as a function of the IMF-index is where i is the IMF-index and f c is a characteristic frequency.This means that the Em Decomposition acts as a nearly diadic filter bank suggesting that the nature of the a 210 series is similar to that of a fractional Brownian motion (Flandrin et al., 2004).a nearly diadic filter bank, suggesting that the nature of the analysed time series is similar to that of a fractional Brownian motion (Flandrin et al., 2004).From these average frequency values fi the corresponding average periodicities can be obtained.They are equal to: T1 = 5 yr; T2 = 10 yr, T3 = 21 yr, T4 = 59 yr and T5 = 90 yr.Similar periodicity values were obtained by Currie (1973), who analysed LoD data using the maximum entropy method and found periodicities of 10-11 yr, 21-22 yr and 57.5 yr.These values are nevertheless partially different from those previously obtained by Roberts et al. (2007), who found clear evidence for two periodicities of 30 yr and 65 yr.While we could partly confirm the longest periodicity of 65 yr having found an IMF, i.e. c 4 (t), with a mean periodicity of ∼ 59 yr, we did not find any evidence for the periodicity of 30 yr by using the Fourierbased approach on all the IMFs.Moreover, the Fourier spectrum of each IMF does not show any characteristic peak at the frequencies found by Roberts et al. (2007).In contrast, the PSDs shown in Fig. 3 are characterized by an energy distributed over a quite large interval of frequencies, suggesting a non-stationary character for the two periodicities indicated by Roberts et al. (2007).We believe that the observed periodicities should more reliably be interpreted under a probabilistic point of view in terms of average characteristic timescales.
To infer the origin of the differences between our results on LoD characteristic periodicities and those obtained by Roberts et al. (2007), we used EMD as already done by Roberts et al. (2007), but with the simple purpose of preprocessing the data.The real step forward in the analysis presented here is represented by the estimation of the Hilbert spectrum of excess LoD that, due to the way it is estimated, also takes into account the uncertainties associated with LoD estimation.By means of the Hilbert spectrum it is possible to move the investigation of excess LoD timescales from the time domain to the time-frequency-energy domain.direct quadrature method, explained in Sect.2, has been performed on each IMF.For each IMF the instantaneous frequency f i (t) and the associated amplitude A i (t) have been evaluated using Eqs.( 3) and ( 4).The results obtained are plotted in the time-frequency-energy space, with energy defined as the amplitude squared.The resulting representation is defined as the Hilbert spectrum H (f, t).Furthermore, to obtain a more robust estimate of the frequencies contained within the excess LoD data, we not only calculated the Hilbert spectrum starting from the excess LoD series, but we also developed a procedure to take into account the time-varying uncertainty associated with LoD estimation.So, we applied the HHT to a large ensemble of LoD i (t) time series with values defined as follows: where k is an index in the interval [1, N ], LoD(t) is the original excess LoD time series, k (t) is a random value from a Gaussian distribution such that the standard deviation of an infinite number of such values would be the measured uncertainty of the LoD(t), i.e. σ (t) = δ LoD(t) where σ (t) is the standard deviation of k (t).This corresponds to considering many of the possibly infinite LoD time series consistent with the estimated errors.This procedure ensures a more reliable estimation of the local frequencies, also providing a right evolution of their significance.Thus, for each time series, both EMD and direct quadrature have been performed to get the corresponding Hilbert spectrum H (f, t) on the time-frequency plane.This procedure has been iterated N = 100 times to explore the sample space of LoD(t) errors correctly.All resulting Hilbert spectra are successively averaged to get the ensemble averaged Hilbert spectrum H (f, t) .This way of applying the HHT must not be confused with the so-called ensemble empirical mode decompositions (EEMD) (Huang and Wu, 2008), where the average is done over the sets of IMFs before computing the Hilbert spectrum.Figure 5 shows the Hilbert spectrum H (f, t) , as obtained by averaging on the ensemble of N = 100 000 Hilbert spectra.The main feature of H (f, t) is its patchy character, which suggests the lack of stationary periodicities in the LoD time series.A more immediate interpretation of the characteristic periodicities contributing to the LoD fluctuations can be obtained by looking at Fig. 6, where the Hilbert spectrum is represented in the time-period-energy domain.
Here, the non-stationary feature of the contributing periodicities is more pronounced.In particular, there is no clear evidence of a rigorously constant periodicity in the range T ∈ [50, 70] yr.For this reason it would be more appropriate to refer to a range of periods within which characteristic periodicities are more likely to fluctuate.Indeed, the Hilbert spectra in Figs. 5 and 6 should be interpreted in terms of probability and statistical weight.This point is clearly shown in Fig. 7, where we have plotted the marginal spectrum H m (f ) computed using the ensemble averaged Hilbert arginal spectrum H m (T ) computed using the ensemble averaged Hilbert spectrum H(T,t) shown 6.The dashed lines correspond to the six Gaussian functions (see expression of Eq.11) whose supereconstructs the marginal spectrum H m (T ).
ely to fluctuate.Indeed, the Hilbert spectra in Figures 5 and 6 should be interpreted in terms bility and statistical weight.This point is clearly shown in Figure 7 where we have plotted inal spectrum H m (f ) computed using the ensemble averaged Hilbert spectrum H(f,t) , ned as follows: and t 2 are the starting and ending values of the time interval over which the analyzed signal d.This spectrum is very revealing about the mean energy distribution during the whole time in the frequency domain.We recall that, differently from the Fourier spectrum, the Hilbert has a statistical meaning.In detail, it is possible to identify 4 characteristic period bands yr, ∼ 48.5 yr, ∼ 68 yr and ∼ 84 yr, respectively.This result can be better visualized by g the marginal spectrum in the period domain using H(T,t) instead of H(f,t) .[15,150] yr.This has been done g the marginal spectrum using the following expression, le contributions are plotted in Figure 8, while in Table 1 are shown the main features of each ting Gaussian.spectrum H (f, t) , and defined as follows: where t 1 and t 2 are the starting and ending values of the time interval over which the analysed signal is defined.This spectrum is very revealing about the mean energy distribution during the whole time interval in the frequency domain.
We recall that, differently from the Fourier spectrum, the Hilbert spectrum has a statistical meaning.In detail, it is possible to identify four characteristic period bands at ∼ 24 yr, ∼ 48.5 yr, ∼ 68 yr and ∼ 84 yr, respectively.This result can be better visualized by evaluating the marginal spectrum in the period domain using H (T , t) instead of H (f, t) .Figure 8 shows the marginal spectrum H m (T ) computed in the period domain.This plot exhibits a more detailed structure than that shown in Fig. 7, showing at least two other possible main typical timescales around ∼ 80 yr and ∼ 90 yr, although most of the energy is concentrated in the range T ∈ [60, 80] yr.
To better identify the different main periodicities, the marginal spectrum H m (T ) has been decomposed into a superposition of single Gaussian bands in the interval T ∈ [15,150] yr.This has been done by fitting the marginal spectrum using the following expression: The single contributions are plotted in Fig. 8, while in Table 1 are shown the main features of each contributing Gaussian.To characterize the temporal behaviour of the identified typical period bands and to evaluate the stationarity of each period band, we have computed the so-called degree of stationarity DS(T ) of each period band.According to Huang et al. (1998), DS(T ) can be defined as follows: where t 1j and t 2j are the boundaries of the period band corresponding to T j ± σ j (as listed in Table 1) and t ∈ [t 1j , t 2j ].
For a stationary signal the degree of stationarity DS is expected to be zero.Departures from zero are a signature of non-stationarity.The last column in Table 1 gives DS values of the main characteristic periodicities.All values of DS are different from zero, indicating that non-stationarity is a feature common to all typical periodicities of excess LoD fluctuations.In particular, the periodicity of ∼ 79 yr characterized by the smallest value of DS seems to be more stationary than the other periodicities.However, we have to take into account the shortness of the LoD time series analysed here, which necessarily influences the value obtained for the degree of stationarity.A good calculation of this quantity must take into account the interval over which it is calculated.This interval should be much longer than the periodicity itself.

Summary and conclusions
The aim of this work is that of investigating the properties of decadal variations in the length of day with the purpose of gaining some new information that can be helpful in the modelling of the LoD and of the torsional oscillations in the fluid core.However, as will be discussed in what follows, our results can contribute to drawing a picture of these variations and of core motions which seem to be more complicated than expected.
In the first part of the paper we applied empirical mode decomposition to the time series of excess length of day.We obtained the same number of monocomponents as obtained by Roberts et al. (2007).The values of the mean frequency

P. De Michelis et al.: Nonstationarity of decadal LOD
of the different modes, which we have estimated by considering the (energy-weighted) mean frequency in the Fourier power spectrum, are also coincident to those obtained by Roberts et al. (2007), who used both the autocorrelation function and a visual method.The only exception is for the value of the mean frequency of the third mode where our results are slightly different.This difference is inherent in the empirical mode decomposition, which produces a set of components (IMFs) from the original data by selecting different parameters independently of each other.When using EMD it is possible to choose a number of variables such as the maximum number of sifting iterations, the stopping criteria as well as the relative thresholds.The obvious (yet critical) question is which set of the many possible choices of the sifting variables gives a meaningful result.For instance, Roberts et al. (2007) used two stopping criteria, each with two different thresholds, for pre-and post-1840 magnetic data to account for the different quality of data.Here, we preferred to use a single stopping condition with a single threshold and move the problem of frequency determinations to the second part of the HHT analysis.Indeed, the ultimate purpose of the EMD is to decompose a given data set into a finite and often small number of intrinsic mode functions that admit Hilbert transforms.We emphasise that the physical meaning of the decomposition into IMF of the original signal comes only in the totality of the decomposed components in the Hilbert spectrum (Huang et al., 1998).This is the reason why we have paid serious attention to the analysis of the Hilbert spectrum and have avoided giving a detailed physical meaning to each IMF component as in Roberts et al. (2007).Indeed, the broadband nature of LoD pointed out by traditional Fourier analysis clearly indicates the absence of a characteristic frequency/periodicity in the investigated range of scales, supporting the idea that the excess LoD fluctuations should not be considered the result of the superposition of a few simple linear waves/oscillations, but conversely are better described in terms of a random superposition of fluctuations characterized by different characteristic timescales.
The really innovative part of our study is contained in the second part of the paper, where traditional analyses are abandoned in favour of the so-called Hilbert-Huang transform (Huang et al., 1998), which we have applied in a completely original way.Taking advantage of the knowledge of errors in LoD measurements, we performed a sort of Monte Carlo approach to EMD Hilbert analysis (this simulation should not be confused with EEMD) to produce the most statistically likely Hilbert spectrum.
The patchy aspect of the Hilbert spectrum immediately suggested the lack of stationarity in the characteristic periodicities contained in the LoD time series.This lack of stationarity was confirmed by the estimation of the degree of stationarity DS for all the period bands found by decomposing the Hilbert marginal spectrum into a superposition of five Gaussian distributions.So, what we are finally capable of establishing is that LoD does not contain stationary periodic-ities, thus making it more appropriate to talk of period bands within which characteristic periodicities are more likely to fluctuate.
This result also implies that it would be more correct to introduce the concept of fluctuation in place of that of oscillation.An oscillation is an intrinsically stationary process.On the contrary, we have found that periodicities of the excess LoD are nonstationary and have a time-varying amplitude.Therefore, rather than considering the processes taking place in the Earth's core as the result of a superposition of oscillations (i.e.stationary mechanisms) occurring on a discrete number of different timescales, we should think of a superposition of fluctuations that are intermittent in both frequency and amplitude, as confirmed by the continuous power law spectrum of Fig. 3.Of course, this implies that processes occurring in the Earth's core are nonlinear, thus making the physical modelling of fluid flows in the core as well as of core-mantle coupling more complicated to perform.
We conclude with a purely speculative conjecture.The presence of a nonstationary and broadband spectrum for the excess LoD could be the counterpart of a wave turbulence process, perhaps a consequence of a turbulent convection phenomenon occurring in the internal fluid core.Clearly, at the present stage this is only a conjecture that requires more studies.
; this gets smaller with time.It is between 0.6 and 0.35 ms in the time interval May 1832-May 1955, between 0.1 and 0.02 ms in the period from May 1955 to May 1997, and around 0.02-0.01ms from May 1998 to May 2009.
, this gets smalle 180 It is between 0.6 and 0.35 ms in the time interval 1832.5-1955.5, between 0.1 and 0.0 period from 1955.5 to 1997.5, and around 0.02-0.01ms from 1998.5 to 2009.5.

185
better define values and stationarity of typical LoD periodicities and to gain a deeper ins underlying processes that are expected to influence Earth's rotation rate.So, we first ap in the standard way without considering the effects of different error measurements and stoppage criterion that discussed in Section 2.EMD yielded five prominent IMFs and a residue, as shown in Figure2.
power spectral density (PSD) of ∆LoD and the PSDs of each IMF.The dashed nlinear best fit of the ∆LoD PSD using a stretched Lorentzian shape (see Eq. 5).

4Fig. 3 .
Fig. 3. Comparison between the power spectral density (PSD) ofLoD and the PSDs of each IMF.The dashed line is a nonlinear best fit of the LoD PSD using a stretched Lorentzian shape (see Eq. 5).
From t frequency values fi the corresponding average periodicities can be obtained.They are = 5 yr; T2 = 10 y, T3 = 21 yr, T4 = 59 yr and T5 = 90 yr.Similar periodicity values were Currie (1973), who analyzed ∆LoD data using the maximum entropy method found of 10-11 yr, 21-22 yr and 57.5 years.These values are nevertheless partially differen 215 previously obtained by Roberts et al. (2007) who found clear evidence of two periodic and 65 yr.While we could partly confirm the longest periodicity of 65 yr having found c 4 (t), with a mean periodicity of ∼ 59 yr, we did not find any evidence of the periodi by using the Fourier based approach on all the IMFs.Moreover, the Fourier spectrum does not show any characteristic peak at the frequencies found by Roberts et al. (2007) 220 the PSDs shown in Figure 3 are characterized by an energy distributed over a quite larg frequencies suggesting a non-stationary character for the two periodicities indicated b al.(2007).We believe that the observed periodicities should more reliably be interpr probabilistic point of view in terms of average characteristic timescales.
ƒ, t)> semble averaged Hilbert Spectrum H(f,t) plotted in the time-frequency plane.Frequency resoluin the ln 2 (f ).The horizontal white line indicates the 65 yr characteristic periodicity.Spectrum H(T,t) reported in the time-period plane.The horizontal white ates the 65 yr characteristic periodicity.

11 Fig. 5 .
Fig. 5. Ensemble averaged Hilbert spectrum H (f, t) plotted on the time-frequency plane.Frequency resolution is 0.1 in the ln 2 (f ).The horizontal white line indicates the 65 yr characteristic periodicity.

Fig. 5 .
Fig. 5. Ensemble averaged Hilbert Spectrum H(f,t) plotted in the time-frequency plane.Frequen tion is 0.1 in the ln 2 (f The horizontal white line indicates the 65 yr characteristic periodicity.
Figure the marginal spectrum H m (T ) computed in the period domain.This plot exhibits a more structure than that shown in Figure 7, showing at least two other possible main typical es around ∼ 80 yr and ∼ 90 yr, although most of the energy is concentrated in the range 80] yr.tter identify the different main periodicities, the marginal spectrum H m (T ) has been decoma superposition of single Gaussian bands in the interval T ∈

Fig. 8 .
Fig.8.Marginal spectrum H m (T ) computed using the ensemble averaged Hilbert spectrum H (T , t) shown in Fig.6.The dashed lines correspond to the six Gaussian functions (see expression of Eq. 11) whose superposition reconstructs the marginal spectrum H m (T ).

Table 1 .
Second and third columns list the main parameters (average and standard deviation) of the Gaussian distributions used to decompose, by means of Eq. (11), the marginal spectrum H m (T ) shown in Fig.8.Fourth column gives the value of the degree of stationarity DS for each period T j .