Detectability of deterministic non-linear processes in Earth rotation time-series-II. Dynamics (1999)

Geophys. J. Int. (1999) 137, 565^579 Detectability of deterministic non-linear processes in Earth rotation time-seriesöII. Dynamics V. Frede1 and P. Mazzega2 1 2 Observatoire de Paris, 61, avenue de l'Observatoire, 75014 Paris, France. E-mail: fred@hpvlbi.obspm.fr UMR 5566-CNRS-CNES,18 av. E. Belin, 31401 Toulouse cedex 4, France. E-mail: ciamp@boreal-ci.cst.cnes.fr Accepted 1998 November 27. Received 1998 November 9; in original form 1998 January 26 S U M M A RY We investigate the possibility of detecting non-linear low-dimensional deterministic processes in the time-series of the length of day (LOD) and polar motion components (PMX, PMY), ¢ltered to keep the period range [&1 day^100 days]. After each timeseries has been embedded in a pseudo-phase space with dimension DE ~5 or 6 (see Frede & Mazzega 1999, hereafter referred to as Paper I) we extract the geometric and dynamical characteristics of the reconstructed orbit. Using a local false neighbours algorithm and an analysis of the data local covariance matrix eigenspectrum, we ¢nd a local dimension DL ~5 for the three EOP series. The principal Lyapunov exponents averaged over the &27 years of observation (1970^1997) are positive. This result unambiguously indicates the chaotic nature of the Earth's rotational dynamical regime in this period range of £uctuations. As a consequence, some theoretical prediction horizons cannot be exceeded by any tentative forecast of the EOP evolution. Horizons of 11.3 days for LOD, 8.7 days for PMX and 8.1 days for PMY are found, beyond which prediction errors will be of the order of the rms of the ¢ltered EOP series, say 0.12 ms, 2.30 mas (milliarcsecond) and 1.57 mas respectively. From the Lyapunov spectra we estimate the Lyapunov dimension DLyap , which is an upper bound for the corresponding attractor dimension DA . We ¢nd DLyap (LOD)~4.48, DLyap (PMX)~4.90 and DLyap (PMY)~4.97. These determinations are in broad agreement with those of the attractor dimensions obtained from correlation integrals, i.e. DA (LOD)~4.5^5.5, DA (PMX)~3.5^4.5, DA (PMY)~4^5. We ¢nally show that the Earth's rotational state experiences large changes in stability. Indeed, the local prediction horizons, as deduced from the local Lyapunov exponents, occasionally drop to about 3.3 days for LOD in the years 1982^1984, 2.6 days for PMX in 1972^1973 and 2.6 days for PMY in 1996^1997. Some of these momentary stability perturbations of the Earth's rotation are clearly related to El Nin¬o events, although others lack an obvious source mechanism. Key words: Earth rotation, El Nin¬o, Lyapunov exponents, non-linear analysis, stability. 1 IN T ROD U C T I O N In this study we investigate the possibility of detecting a low-dimensional deterministic process in data time-series of the earth orientation parameters (EOP) which may come from non-linear atmospheric and oceanic excitations. Daily time-series of the polar motion components (PMX, PMY) and length of day (LOD) from 1970 to 1997 (10 080 samples in each series) are analysed. They result from a combination of individual series determined by various geodetic techniques (VLBI, SLR, GPS, etc.; Gambis 1996). We focus on EOP £uctuations with periods between a few days and 100 days. In this second part of our study, our goal is to ¢nd various ß 1999 RAS dynamical characteristics of the time-dependent Earth orientation from the data series alone. We will compare them with the characteristics (attractor dimensions, eventually horizon of prediction, see below) obtained by several authors for various timescales of the atmospheric (see e.g. Nicolis & Nicolis 1984; Essex et al. 1987; Fraedrich 1987; Tsonis & Elsner 1988) and oceanic (see e.g. Vallis 1986; Jin et al. 1994; Tziperman et al. 1994) £ows. In Paper I we obtained dimensional characterizations of the short-term variations in Earth rotation. First we estimated optimal time delays required to build multidimensional vectors y(n) from the scalar time-series of EOP in a system of statistically independent coordinates. From the averaged 565 566 V. Frede and P. Mazzega mutual information function (Fraser & Swinney 1986) we found time delays of 10, 15 and 18 days for the LOD, PMX and PMY series respectively. In order to apply the embedding theory to our data sets (Takens 1981) we then needed to retrieve the dimension DE of a pseudo-phase space where any reconstructed orbits from the y(n) vectors will be unfolded. Two approaches were tested for this task: we investigated changes in the global percentage of false neighbours (Abarbanel 1996) and changes in the slopes of the correlation integrals (Grassberger & Procaccia 1983) when the trial embedding dimension was incremented. From both methods we determined the embedding dimension DE to be between 5 and 7. This result agrees with the hypothesis of a EOP £uctuation source of low dimension and deterministic nature. However, we must bear in mind that it was impossible to give a precise value for DE because the EOP series are quite short and noisy. Then, from 3-D plots of 6-D reconstructed orbits (Paper I, Fig. 10), we pointed out that the general behaviour of a high-dimensional process such as a white or coloured noise is fundamentally di¡erent from those we obtained with the EOP time-series. While a stochastic noise ful¢ls the 3-D subspace in an erratic way or traces a fractal curve, the embedded EOP series trace a smooth trajectory which might originate from a di¡erentiable dynamical system. The goal of this second paper is to characterize the dynamics of the Earth rotation £uctuations through determinations of the dynamical (or local) dimension DL , the spectrum of the global Lyapunov exponents ji (with i~1, DL ), the timedependent principal local Lyapunov exponent j1 (t) and the attractor dimension DA . DE is the global dimension of the embedding space in which any reconstructed orbit from vectors of time-delayed coordinates is unfolded. The dimension DL gives the e¡ective number of dynamically active variables. It is obtained from the local distribution of the data vectors y(n) on the attractor. While the analysis of two time-series originating from the same source may lead to di¡erent embedding dimensions DE , the dynamical dimension DL will be the same. This property is another consequence of the embedding theorem (Takens 1981; Sauer et al. 1991). In any case, DL ¦DE gives the exact number of ordinary di¡erential equations one needs to describe the underlying system and capture its full dynamics. The spectrum of Lyapunov exponents ji provides useful information on the global stability of the system. Each exponent expresses how a small perturbation of the orbit in an eigendirection of the dynamics will evolve in time. Directions associated with null exponents are neutral in the sense that any perturbation applied along them will remain constant. Those directions associated with negative ji are stable: any perturbation component lying on them will asymptotically die out. Positive exponents characterize unstable directions along which a perturbation will grow exponentially during a ¢nite time. The rate of growth is precisely related to the exponent. Note that these directions change when going over the orbit. DL exponents can be computed from a reconstructed orbit in the embedding pseudo-phase space E*. A system which possesses only negative or null exponents is stable, whereas a chaotic system (characterized by its sensitivity to initial conditions) has at least one positive exponent. The Lyapunov exponents are invariants of the dynamics, so their determination does not depend on the particular orbit used to compute them. Of course, this statement should be treated with caution when dealing with orbits reconstructed with ¢nite and noisy data time-series or, from a dynamical point of view, when several attractors coexist in phase space, as will be discussed later. If the sum of the Lyapunov exponents is negative, the dynamical system is dissipative. The ability to predict the evolution of a time-series is intimately related to its principal Lyapunov exponent j1 . If j1 > 0, any error made on the present state of the system will grow quickly and reach unacceptable bounds. However, by scanning in E* the divergence of nearby orbit segments at various places along the attractor, it can be shown that the principal exponent experiences local variations. In other words, the stability of the system and its predictability depend on the system state or equivalently on time. The interest in computing local exponents j1 (t) is thus twofold: ¢rst it sets a theoretical horizon of prediction that cannot be exceeded for a given state of the observed system; second, for a nonautonomous system, it may hopefully provide some information on the time-varying external forcings which prompt system instabilities. If the system possesses an attractor in phase space, the determination of its dimension is also useful information. The attractor is a geometric ¢gure with dimension DA < DE on which all the orbits in phase space converge rapidly after transients have died out. For example, the attractor of a twoperiodic system is a 2-torus which captures the asymptotic behaviour of the system evolution. When the attractor dimension is an integer, we are dealing with periodic orbits, while a fractal dimension denotes the attractor `geometric' strangeness. Dynamical chaos is typi¢ed by the sensitivity of the system state evolution to initial conditions. At this stage non-linear analysis of the data series must be conducted very carefully. Indeed, we cannot conclude from the non-integer determination of DA the chaotic nature of the dynamical regime for two reasons: ¢rst, given the uncertainties in the evaluation of DA , the recovery of the integer dimension of, say, an n-torus is practically impossible from a ¢nite data set; second, there exist some strange non-chaotic attractors (with negative principal Lyapunov exponent; Ding et al. 1989). In any case, there is a fundamental di¡erence between chaotic regime and fractal attractor dimension. A fractal dimension DA is only relevant for characterizing the self-similarity of the orbit (for example of a random walk) or the way it ¢lls the phase space (for example on the Lorenz attractor). The notion of a chaotic regime adds the idea of sensitivity to initial conditions (Eckmann & Ruelle 1985; Osborne & Provenzale 1989). A complete dynamical characterization of any system is clearly required before giving any conclusion on the nature of its dynamical regime. All these dynamical characteristics belonging to the dynamical system theory are asymptotic quantities and assume in¢nite time spans to apply rigorously to a problem. We will thus meet several di¤culties whilst working with ¢nite and noisy time-series. As in Paper I, the robustness of the results produced by the diverse algorithms will be checked with the help of various tests concerning their sensitivity to data contamination and length of the time-series or sampling intervals. The dynamical analysis of well-known series will help us in our investigation. In Section 2, we will estimate the local dimension DL from the three EOP data series. Section 3 is devoted to the computation of the global Lyapunov exponents and constitutes a ß 1999 RAS, GJI 137, 565^579 Non-linear processes in Earth rotation time-seriesöII major step in this study. In Section 4 we estimate the time variations of the principal Lyapunov exponent j1 (t) over the 27 years of observations for the polar motion components and length of day. Section 5 deals with the determination of the attractor dimension. Besides answering the question of the detectability of a non-linear low-dimensional process in Earth rotation, several results of this study have fundamental implications concerning the nature and geophysical origin of observed £uctuations with periods shorter than 100 days; these are discussed in Section 6. 2 D ET E RM I NAT I O N OF TH E DY NA M ICA L D I M E N S I O N DL The embedding operation gave us the global dimension DE needed to unfold properly the orbits reconstructed from data vectors in delayed coordinates. Now we want to pinpoint the exact number of dynamical variables which are really active during the evolution of the system. This dynamical or local dimension DL is lower than or equal to DE. In order to illustrate the di¡erence between the embedding and local dimensions, let us consider a quasi-periodic system composed of two incommensurable frequencies. The trajectories lie on a 2-torus. We need a 3-D phase space to unfold this torus (DE ~3), whilst the e¡ective number of active dynamical variables is DL ~2. This is seen geometrically because a local proxy for the shape of the attractor is its tangent plane. The local dimension DL can also be estimated from observations of a single variable. It can be conceived as the number of autonomous ordinary di¡erential equations required to describe the source dynamics (and therefore the number of Lyapunov exponents to be determined). The method employed to compute the local dimension is the local false nearest neighbours (LFN) algorithm (Abarbanel & Kennel 1993), which looks locally at the dimension required to represent unambiguously the dynamical evolution. The method relies on the determination of local maps from neighbourhood to neighbourhood along the reconstructed orbit. These maps are then used to make local predictions of the orbit evolution. The predictions are statistically reliable only when the trial dimension used in the mapping is the local dimension. As will be shown below, we also found it very useful to look at the eigenspectrum of the local data covariance, so we add this option to the basic LFN algorithm. In practice, the LFN algorithm selects the local principal directions in phase space where most of the data lie. For this purpose we embed the data in a working space of dimension Dw higher than DE (here we choose Dw ~10). Taking a point of the pseudo-orbit y(i) and its No nearest neighbours yo (i), we form the data covariance matrix, i.e. C(i)~ No 1 X [yo (i){ym (i)] [yo (i){ym (i)]T , No o~1 (1) with the local mean vector ym (i)~ No 1 X yo (i) . No o~1 (2) We then calculate the eigenvalues of C(i) arranged in a non-increasing order: j1 §j2 §


§jn , ß 1999 RAS, GJI 137, 565^579 (3) 567 where n is the value of the trial local dimension. This operation is performed for each point of the trajectory. We ¢nally calculate the mean eigenvalues over the whole orbit and normalize them by the largest value to get the mean relative eigenvalues (MRE). In this approach the local dimension DL is given by the number of signi¢cant eigenvalues. They are associated with directions of phase space where data points y(i) really accumulate. Eigenvalues that are too small are not relevant and re£ect numerical approximations, covariance matrix illconditioning, truncation errors, data quantization or observational noise. For example, the local dimension of the Lorenz system is 3, whilst the attractor dimension is 2.06. If we calculate the MRE spectrum, we select two large eigenvalues and a smaller one which corresponds to the fractal direction of the Cantor set structure of the attractor. This third direction is much less populated as the Lorenz attractor locally lies in a quasi-plane. Therefore, its detection will be more sensitive to noise contamination. There is no general criterion to decide about the irrelevance of some part of the eigenvalue spectrum, so we shall ¢nd a threshold value for the signi¢cance of the MRE from the analysis of our contamination series. The next step in the basic LFN algorithm consists of computing the projected vectors yo' (i) of the No neighbours of y(i) on the eigenvectors. After k time steps g taken along the £ow they become the projected vectors yo' (izkg). We determine by a least-squares method the local linear map which transforms the yo' (i) into the yo' (izkg). The time evolution of the centre point y(i) is then predicted from the local mapping and compared to the true vector y(izkg). We look at the reliability of the local maps computed for every y(i) by examining the rms prediction error. We ¢x a threshold over which any prediction is labelled a bad one. Increasing the number of eigenvectors selected, we reach the local dimension DL when the percentage of bad predictions becomes independent of the number of local eigenvectors used in the mapping process. The choice of a number of neighbours No adapted to our problem was made by testing it on the Lorenz series lor (see Paper I, Section 2) of 10 080 samples. The expected local dimension DL ~3 was retrieved for No ~10 or 25. As we are working locally, a lower number of neighbours will lack statistical signi¢cance and will provide results much more sensitive to any noise contamination. A higher number would require polynomial mappings of higher degree to account for the curvature e¡ects on the attractor. The test performed with the Lorenz series also showed that the determination of DL based on the eigenvalues of the covariance matrix gives better results than the percentage of bad predictions, which has a marked tendency to increase the estimates of the local dimension. This is due to the limited number of data points available (10 080). Indeed, the same test performed on a series of 30 000 points gives the correct result. Moreover, we applied the LFN algorithm to a white noise series of 30 000 points and found no local dimension, the percentage of bad predictions remaining close to 100 per cent whatever the number of retained eigendirections. This is related to the fact that points are randomly distributed without preferred directions. Like for the Lorenz series, working with a short series ( f 0, 10 080 samples) gives spurious results by showing a ¢nite local dimension for a white noise. This drawback is signi¢cantly reduced when considering the MRE spectrum. The algorithm was also applied to the raw data of the random walk time-series f {2 (10 080 samples) and to the 568 V. Frede and P. Mazzega same data ¢ltered so as to keep the 1^100 day spectral band. The percentage of bad predictions for the raw coloured noise starts at only 3 and decreases slowly to zero when the trial local dimension is 6. Such behaviour is due to the fractal structure of the curve, which implies self-similarities characteristic of the Cantor set structure of an attractor. Nevertheless, there is no low-dimensional process here. This spurious result is suppressed when working with the ¢ltered time-series. Indeed, the fractal structure is lost and the ¢ltered f {2 series behaves more like a white noise. The percentage of bad predictions remains at 100 per cent for a trial local dimension of 9 before dropping suddenly to zero (probably as a consequence of the limited number of samples). Such behaviour is not characteristic of chaotic systems. Here a striking di¡erence occurs when performing the same test with a chaotic series. The percentage of bad predictions is nearly invariant when dealing with the Lorenz series lor or with its ¢ltered version. In both cases the correct local dimension DL ~3 is found, so we must be very careful if the percentage of bad predictions we obtain is very low in dimension 1 and decreases rapidly to zero, because this can be either a fractal curve or a chaotic process. The way to distinguish between the two possibilities is to analyse the high-pass-¢ltered time-series also. We now apply the LFN algorithm to the EOP series and to their companion series contaminated by power-law noise or low-dimensional processes. The three series have the same kind of behaviour (Fig. 1a). The percentage of bad predictions is quite high at the beginning and decreases progressively down to zero with increasing trial local dimension. DL as estimated from this analysis alone would be around 6^8 for the PMX and PMY series and 7^9 for LOD. The e¡ect of adding 10 per cent noise is, not surprisingly, enormous. Indeed, the local dimension determination is very sensitive to noise contamination because we are working in very con¢ned neighbourhoods (No ~25). Locally the noise has a large e¡ect by perturbing the optimal mapping recovery. The algorithm partly tries to predict the noise. It also becomes very di¤cult to choose a valid prediction horizon (kg in the explanation above). If the horizon is taken too short, we predict the noise only. If it is too long, we are faced with curvature e¡ects and orbit segment divergence that can be tackled with higher polynomial maps only, which in turn require a huge number of data samples in order to be constructed reliably. Adding 10 per cent white or coloured noise to the LOD time-series thus totally changes the behaviour of the percentage of bad predictions, as shown in Fig. 1(b). The noise is dominant and we cannot see the signal any more. This important remark shows us that if our data series are even slightly contaminated by noise, the DL estimated via the LFN algorithm will be signi¢cantly overestimated. On the other hand, the addition of 10 per cent Lorenz chaos to the LOD data does not have any signi¢cant e¡ect on the percentage of bad predictions. The same results are obtained with the PMX and PMY series. This test has shown that the relative behaviour of the EOP time-series and that of the high-dimensional timeseries are completely di¡erent. Moreover, the previous values for the EOP local dimensions are probably overestimated. Inspection of the mean relative eigenvalues (MRE) spectrum o¡ers a geometrical way to determine the local dimension DL . It consists of ¢nding the local directions along which the data are preferentially distributed. This test was performed on the original and contaminated EOP time-series as well as on the synthetic noise and Lorenz chaos. The main limitations of DL Figure 1. Percentage of bad predictions as a function of the trial local dimension DL . (a) For the LOD, PMX and PMY series. (b) For LOD (solid line) and its contaminated equivalents: zf 0 (dotted line), zf {2 (dashed line), zlor (long-dashed line). this algorithm come from the great local variations observed for the eigenvalues. A robust estimate of the MRE spectrum would require the use of time-series long enough to ¢nd the asymptotic characteristics of the data distribution in phase space. Moreover, the determination of very small eigenvalues is very di¤cult and the associated eigenvectors are given mostly by local noise. Note that this weakness of the MRE approach is shared with the LFN algorithm. Nevertheless, this method seems to give better results than the prediction method with our limited series. The local dimension estimate was more relevant for known series (e.g. for the Lorenz series lor and noises f 0 and f {2). We also checked that this geometrical method was more robust in the presence of noise; this is probably because no least-squares ¢t to the data for the determination of mapping coe¤cients is involved, nor any prediction skill. The MRE spectra for the three EOP series are shown in Fig. 2(a). The local dimensions obtained from this analysis are around 4^5 for the PMX and PMY series and 5^6 for LOD. They are systematically lower by two units than the determinations found with the percentage of bad predictions. We are more con¢dent in the results of the MRE analysis as it appears less biased by noise contamination. In Fig. 2(b) the MRE spectra for the LOD and for its contaminated companions are shown. The local dimension is overestimated in ß 1999 RAS, GJI 137, 565^579 Non-linear processes in Earth rotation time-seriesöII 569 Table 1. Summary of the dimensional and dynamical results for the ¢ltered time-series of the LOD and polar motion components PMX and PMY: time delay q in days; embedding dimension DE; local dimension DL ; principal Lyapunov exponent j1 ; sum of the Lyapunov exponents ji ; Lyapunov dimension DLyap ; attractor dimension DA ; prediction horizon HP in days. q (days) DE DL j1 & ji DLyap DA HP (days) LOD PMX PMY 10 5=6 5 0:35 {0:41 4:48 4:5=5:5 11:3 15 5=6 5 0:45 {0:08 4:90 3:5=4:5 8:7 18 5=6 5 0:48 {0:02 4:97 4=5 8:1 perturbation lying on the corresponding stable manifold. On the other hand, if the system possesses at least one positive Lyapunov exponent ji , its time evolution is sensitive to the choice of initial conditions in a particular direction. Any small perturbation d(0) exerted on the orbit at t~0 along that direction grows exponentially with time, so that d(t)~d(0) eji t . DL Figure 2. Mean relative eigenvalues of the local data covariances as a function of the trial local dimension DL . (a) For LOD, PMX and PMY. (b) For LOD (solid line) and its contaminated equivalents: zf 0 (dotted line), zf {2 (dashed line), zlor (long-dashed line). the presence of noise, although this artefact is less marked that with the LFN algorithm (compare with Fig. 1b). Once again the chaotic contamination has no signi¢cant e¡ect on the determination of DL . This possibly indicates the presence of low-dimensional chaos in the EOP series. In summary, the most reliable estimates are as follows: (1) DL ~5 or 6 for LOD; (2) DL ~4 or 5 for PMX; (3) DL ~4 or 5 for PMY. As a common local dimension we report the value of DL ~5 in Table 1, which summarizes the dimensional and dynamical characteristics we obtain for the three EOP series. 3 G L O BA L LYA P U NOV E X P O N E N TS The global Lyapunov exponents characterize the global stability of the regime of a dynamical system (see e.g. Thompson & Stewart 1987; Tu¢llaro et al. 1992). Each exponent indicates the rate of contraction or expansion of an orbit perturbation taken on a local eigendirection of the system dynamics (see below). A negative exponent implies a contraction of the orbit ß 1999 RAS, GJI 137, 565^579 (4) A zero Lyapunov exponent is always associated with a continuous dynamical system (a £ow). It corresponds to the absence of both contractions and expansions if the perturbation is exerted along the orbit itself (neutral manifold). Systems whose sum of Lyapunov exponents is negative are dissipative. Any volume of the n-dimensional phase space will shrink with time, although some directions may remain unchanged (along neutral manifolds) or even stretch (along unstable manifolds). For example, k-periodic systems have k zero Lyapunov exponent and n{k negative exponents. Initial conditions taken on a stable manifold will die out at an exponential rate and the orbit will `stick' to the attracting k-torus. The existence of positive exponents is typical of the chaotic nature of the system regime. Perturbations taken `outside' the attractor will be damped but two orbits will rapidly separate from each other on the attractor because of the instability. We want to ¢nd the global Lyapunov spectra of the Earth rotation dynamics from the EOP time-series of polar motion and LOD. We expect to ¢nd at least one zero exponent because we are dealing with a £ow underlying our data. Indeed, the atmospheric and oceanic excitations are ruled by £ows and Liouville's equations describing the Earth's rotational response are themselves a set of di¡erential equations. The accuracy of the determination of the zero exponent will give us an idea about the reliability of the Lyapunov spectrum recovery. The positive exponents are quite easy to ¢nd from a time-series because they correspond to dilatation directions of the nearby segments of the reconstructed orbit. On the other hand, it is di¤cult to determine the negative exponents because they correspond to the contraction of orbit perturbations, so their e¡ect disappears rapidly in time. As a consequence, they are associated with very small eigenvalues of the Jacobian matrix entering into their determination, which makes these matrices ill-conditioned. The global Lyapunov exponents are invariants of the motion evolving over the in¢nite time axis. 570 V. Frede and P. Mazzega This is a key point because working with a single time-series realization the global Lyapunov exponents we calculate are representative of the entire dynamical system regime. Several methods have been proposed to ¢nd the Lyapunov spectrum from a time-series (Eckmann et al. 1986; Bryant et al. 1990; Brown et al. 1991). The updated method we use is detailed in Abarbanel (1996). Consider a dynamical system with vector of state variables yi at a discrete time i described by a map yiz1 ~F ( yi ) , (5) where F is a di¡erentiable function. Using linearized dynamics, a small perturbation di applied to the orbit yi at time i will become at time iz1 diz1 ~DF ( yi )di , (6) where DF ( yi ) is the Jacobian matrix of F at yi . After L time steps we have dizL ~DF ( yizL{1 )DF ( yizL{2 )


DF ( yi )di :DF L ( yi )di . (7) The Lyapunov exponents are directly linked to the eigenvalues of the previous composition of Jacobian matrices. In order to calculate them, we form the symmetric Oseledec matrix: S( yi , L)~[DF L ( yi )T DF L ( yi )]1=2L , (8) where the superscript T is for transpose. In theory, the Oseledec matrix is built from an in¢nite sequence of matrix composition (L??) so that the choice of starting point yi is of no concern. In practice we start at the beginning of the embedded data series and compose the local matrix over the ¢nite series length. Therefore, we do not really obtain the global Lyapunov exponents but rather averaged exponents which are representative of the average stability of the system over the observation time span. When working with a time-series, we ¢rst have to determine the successive unknown Jacobian matrices. This is done by looking at a point in the pseudo-phase space and at its nearest neighbours. We then de¢ne a local map which models the evolution of the neighbours in time: yoi ?yo;iz1 , (9) where yo;iz1 is the oth neighbour of yi at time iz1. The local map is built as yo;iz1 ~ M X c(m, i)'m ( yoi ) , (10) m~1 where the 'm are a priori basis functions and the c(m, i) are local coe¤cients to be determined. The (a, b) component of the Jacobian matrix is de¢ned by DFab ( yi )~ M X m~1 ca (m, i) L'm ( yi ) . Lxb (11) ja ~ lim L?? 2L 1 X log[Raa (k)] , 2L k~1 (12) (13) where Raa (k) is the (a,a) component of the kth triangular matrix R(k). We see that this is an asymptotical limit, the matrix Q going to identity. Thus for ¢nite data sets, this condition cannot be achieved and we have to apply the recursive QR decomposition iteratively as explained in the next section. We have determined the averaged Lyapunov exponents for the EOP series and tested the result sensitivity to various contaminations. Dealing with a local dimension DL ~5, we evaluate both the forward and backward exponents in order to eliminate the possible false exponents. A true exponent changes sign when time is reversed (Abarbanel 1996). For PMX, PMY and LOD the ¢ve exponents are true and we can ¢nd two positive exponents. The principal exponents are 0.35, 0.45 and 0.48 for LOD, PMX and PMY respectively. We also obtain for the three series the expected zero exponent corresponding to the £ow direction of the data, with an accuracy of about 5 per cent. The accuracy of the positive exponent determinations is probably slightly better. We also obtain two negative exponents as illustrated in Fig. 3 for LOD. L The negative sum &D a~1 ja of the global exponents indicates the dissipative nature of the system (Table 1). The Lyapunov dimension DLyap (Kaplan & Yorke 1979) is directly evaluated from the exponents by K X ja DLyap ~Kz a~1 jjKz1 j , (14) Kz1 with &K a~1 ja > 0 and &a~1 ja < 0. The integer K is chosen so that K < DL and jKz1 < 0. DLyap is an upper bound of the attractor dimension DA . It is a dynamical characteristic of a chaotic motion. For its calculation, one needs to know the full Lyapunov spectrum, including the negative exponents. Working with noisy data and ¢nite series may be an obstacle to obtaining a precise determination of either of the negative exponents and then the Lyapunov dimension. For the Lorenz system, the Lyapunov exponent is exactly the same as that of the attractor (DA ~2.06), a value we recover with the lor series of &104 samples only. The Lyapunov dimensions obtained for the EOP time-series are as follows: (1) DLyap ~4:48 for LOD; (2) DLyap ~4:90 for PMX; (3) DLyap ~4:97 for PMY. Finally, horizons of prediction HP (in days) are evaluated for each EOP time-series from their respective principal exponents j1 by HP ~3:9/j1 . From this empirical Jacobian, we form the Oseledec matrix. Its ill-conditioning prevents the direct computation of its eigenspectrum, so we perform a recursive QR decomposition (Golub & van Loan 1989) leading to the following product of an orthogonal matrix Q and 2L upper-right triangular matrices R(:): S~QR(2L)R(2L{1)


R(1) . The Lyapunov exponents are computed from the diagonal terms of the R matrices as follows: (15) The factor 3.9 scales the predicting horizon to be the point where two trajectories, initially distant by 1 per cent of the mean attractor size RA , separate from each other by 0.5 RA (here we take the data series rms for RA ). HP is a crude estimate of the time until which it is possible to make a valid prediction of the system evolution. Over this time, the forecast error is of the order of the data rms. For the time-series ¢ltered at 100 ß 1999 RAS, GJI 137, 565^579 Non-linear processes in Earth rotation time-seriesöII 571 Global Lyapunov Exponents Figure 3. Values of the global Lyapunov exponents for LOD. Each global exponent is obtained as the average of 3|103 local exponents computed on the basis of the composition of L Jacobian matrices estimated at the rate of one per day (see text). L is given on the abscissa. days, the predicting horizon is about 11.3 days for LOD (recalling that the rms~0.12 ms), about 8.7 days for PMX (rms~2.3 mas) and 8.1 days for PMY (rms~1.6 mas). It has recently been shown that it is possible to ¢nd spurious positive exponents from random time-series (Tanaka et al. 1996). We have thus performed the same computations for the EOP time-series contaminated by noise and by the Lorenz series (see Paper I, Section 2). The addition of 10 per cent white noise greatly perturbs the determination of the Lyapunov exponents, but the results are clearly no longer valid because L &D a~1 ja > 0 (system neither conservative nor dissipative). This comes ¢rst from the fact that a white noise is a highdimensional process which cannot be embedded and second from the fact that the white noise component completely hides the underlying smooth EOP signal because the determination of the Jacobians is a local operation that is very sensitive to erratic point distribution over the attractor. The addition of a coloured noise f {2 or of a Lorenz signal does not signi¢cantly change the results. It thus seems that we could mistake a random walk for a chaotic signal. In fact, such a coloured noise might be detected by analysing the data series with reversed times. Indeed, in such a situation, after the normal sign correction accounting for the time reversal was applied, we obtained a positive sum of the exponents with the series f {2. Time reversal has been tested with the EOP series and never causes such an unphysical result. The Lorenz conß 1999 RAS, GJI 137, 565^579 taminated series not only successfully pass this test, but their principal Lyapunov exponents are changed by no more than 10 per cent. In summary, the estimates we obtained for at least the principal Lyapunov exponents j1 are quite robust with regard to noise contamination. They are positive for the three (¢ltered) EOP series. This result characterizes the dynamical chaos of the Earth rotation £uctuations in the period range [&1^100 days]. As a consequence, theoretical bounds (related to the system dynamics) exist beyond which no reliable forecast of the Earth's orientation variations at short periods can be given. In practice, an error of &0:001 ms on LOD will grow to 0:06 after 11 or 12 days. In the same way, errors of [0.02, 0.015] mas on the polar motion components [PMX, PMY] will reach about [1.15, 0.75] mas after &8.5 days. 4 L O CA L LYA P U NOV E X P O N EN T S Being invariants of the dynamics, the global Lyapunov exponents characterize the attractor as a whole. There can exist local variations of the system stability along an orbit. Those local variations are quanti¢ed by the local Lyapunov exponents (Abarbanel et al. 1991). They vary with the position of the state vector yi in phase space. Their knowledge is very useful for short-term predictions. If the principal exponent j1 is positive, its variations [j1 (t) hereafter] give equivalent time 572 V. Frede and P. Mazzega variations of the predicting horizon HP (t). System stability changes of this kind are well known, for example in large-scale atmospheric dynamics (see e.g. Palmer 1993; Palmer et al. 1994). In the previous section, we ¢nished the algorithmic part after the Oseledec matrix (eq. 8) was recursively decomposed (eq. 12). For convenience we call the `length' of the Oseledec matrix, S( yi , L)~[DF L ( yi )T DF L ( yi )]1=(2L) , the number L of local Jacobian matrices DF entering its composition, and yi its starting point on the attractor. The most reliable `global' Lyapunov exponents discussed in Section 3 are formed by averaging the local exponents of the Oseledec matrices of length L~512 estimated at 3000 starting locations. For short time spans (small L), the basic recursive QR decomposition is inadequate for providing accurate estimates of the Oseledec matrix eigenvalues. Indeed, the orthogonal Q matrix (eq. 12) does not ¢t the identity matrix. Applying the recursive QR decomposition iteratively overcomes this di¤culty, although preserving the eigenvalues of the decomposed matrix (Abarbanel 1996). In Section 3 we showed that S0 ~Q0 R0 (2L)R0 (2L{1)


R0 (1) , (16) where the zero subscript stands for the 0th iteration. We then de¢ne the matrix S1 ~Q0 T S0 Q0 (17) and perform a new recursive QR decomposition on S1 . We repeat this process until QK ¢ts the identity matrix to the desired accuracy. The sequence of Sk~1,K matrices have the same eigenvalues. The local Lyapunov exponents are obtained at the Kth iteration from the diagonal terms of the R matrices: ja ~ 2L 1 X log[RK[aa] ( j)] . 2L j~1 (18) With our time-series, the Q matrix converges to the identity within an accuracy of 10{6 after two iterations at most. Our aim is now to look at the variations of the (local) principal exponents j1 (t) over the years of EOP observations. All the results presented below are obtained with 28 years of observation (10 234 samples for each EOP series), ending in December 1997. This slightly longer EOP series results from minute di¡erences in data combination and ¢ltering routines (new combined IERS series). The previous dimensional and dynamical characterizations were recomputed for these slightly longer time-series without noticeably changing the results. We shall focus on those exponents obtained from local matrices of length L~8 and L~64 which are representative of stability £uctuations of the corresponding number of days. The 8 day stability is chosen so as to ¢t broadly the embedding time delay corresponding to the minimum interval required to have non-linearly independent data (see Paper I, Section 4). The 64 day stability is studied because we lose any mutual information between the successive measurements over 2 months and the corresponding local exponents will smooth out somewhat £uctuations in Earth's rotational stability that are too rapid. The time-series of local exponents j1 (t) for LOD is given in Fig. 4. We immediately note important variations with regard to the average exponent. The most striking feature is the peak for both the L~8 and the L~64 exponents in the year 1983. Throughout 1983, the local Lyapunov exponents are between two and three time their mean values, and consequently the corresponding horizon of prediction is divided by the same factor. Obviously, this transitory relative loss of stability in the Earth's rotational dynamics must be related to the strongest El Nin¬o on record. No other exceptional event is seen in this plot. Using the previous shorter LOD series of 10 080 samples, we obtain a second, although less marked, peak for the year 1977 Figure 4. Time variations of the LOD principal local Lyapunov exponent from 1970 to 1998 for a length L~8 days (solid line) and L~64 days (long-dashed line). ß 1999 RAS, GJI 137, 565^579 Non-linear processes in Earth rotation time-seriesöII 573 Figure 5. Time variations of the PMX principal local Lyapunov exponent from 1970 to 1998 for a length L~8 days (solid line) and L~64 days (long-dashed line). which might correspond to the mild El Nin¬o observed at that time in the sea surface temperature and in the sea level pressure. Local exponents for the PMX and PMY series are shown in Figs 5 and 6 respectively. The major event in the stability of the x-component of polar motion lasts over three successive years, from 1971 to 1974. The magni¢cation factor is larger than 3 (and the prediction horizon reduces to about 2.9 days only). A candidate excitation could be the two strong changes in polarity of the tropical ocean^atmosphere coupled system, with two La Nin¬a events £anking a well-marked El Nin¬o (1973). No other clear stability signal can be matched to known Figure 6. Time variations of the PMY principal local Lyapunov exponent from 1970 to 1998 for a length L~8 days (solid line) and L~64 days (long-dashed line). ß 1999 RAS, GJI 137, 565^579 574 V. Frede and P. Mazzega atmospheric or oceanic excitation except perhaps the 1992 El Nin¬o event where the prediction horizon seems to be reduced by a factor 2. More variability is seen in the local Lyapunov exponents of the PMY series (Fig. 6). Two relative losses of stability occur from mid-1975 to the end of 1976, and from mid-1982 to the end of 1983. If the second anomaly is induced by the famous El Nin¬o, the ¢rst could perhaps be related to La Nin¬a^El Nin¬o twins for the same period. No clear signature for the 1997 El Nin¬o appears in the EOP local exponents. An indistinct rise of both local exponents seems to start in Fig. 5 for PMX. This very unclear anomaly has yet to be con¢rmed during the ensuing few months, while El Nin¬o might be in partial retreat since December 1997. Remembering that we are dealing with EOP time-series high-pass ¢ltered at 100 days, the variability that we observe and link to external excitations is the result of geophysical events on that particular timescale. Even if the El Nin¬o event lasts several months, we here recover its signature on short periods in the Earth rotation. The predicting horizon we give is thus for the signal under 100 days only; it is obvious that working with raw data would give di¡erent values for the horizon of prediction and also a much larger attractor mean size, or equivalently a larger bound for any error growth. Stronger magni¢cation of the PMY local exponents occurs in 1996 (Fig. 6). No signi¢cant anomaly in the ocean^atmosphere circulation can be invoked as a candidate source for this particular perturbation in Earth rotation stability. In a recent catalogue of earthquake magnitudes (Engdahl et al. 1998) we found ¢ve very powerful events with magnitude >7.5 occurring in 1996 [01/01 with (north latitude, east longitude; magnitude) ({20:4, {174:2; M~7:9); 17/02 with ({00:9, z136:9; M~8:2); 10/06 with (z51:6, {177:6; M~7:9); 17/06 with ({07:1, z122:6; M~7:9); 12/11 with ({15:0, {075:7; M~7:7)]. This suggested relation between earthquake energy release and polar motion £uctuations has been studied in the past (see e.g. Chao & Gross 1987). Whilst no evidence for such a mechanism was found from the observations, the theoretical studies indicated that earthquake energy output larger by an order of magnitude or two would be required to pass the threshold of detectability with the current measurement techniques. The accuracy of the combined EOP series has been much improved during the last 10 years. Moreover, we do not correlate earthquakes with EOP measurements but with relative stability losses in rapid £uctuations of the Earth rotation. We are wary of drawing premature conclusions from this, but feel that there is a need to re-investigate this question. In Table 2, we report several such transitory relative losses of the Earth's rotational stability as detected from the local Lyapunov time-series j1 (t) for LOD, PMX and PMY. Most of them have an obvious signature on one or two of the EOP `stability' series only. We tested the robustness of these peaks in the j1 (t) by performing the same analysis with the contaminated EOP series (see Paper I, Section 2). All the events discussed above are still visible in the presence of 10 per cent rms of stochastic or deterministic contamination. The added white noise has the worst e¡ect by raising the background level of the j1 (t) series so as to partially merge the peaks. However, the signi¢cant peaks are still clearly seen with preserved values of the associated prediction horizons HP (t). The e¡ect of the random walk f {2 contamination is to increase spuriously some of the smaller peaks and reduce the Table 2. Date of events related to local variations of the principal Lyapunov exponent j1 (t) for LOD, PMX and PMY time-series. The upper number is the value of the local Lyapunov exponent and the lower one is the associated horizon of prediction HP in days. The last line gives the global principal Lyapunov exponent and the corresponding horizon of prediction for the 1970^1998 time interval. Series LOD PMX PMY 1972=1973 0:5 7:8 1:5 2:6 0:48 8:1 1976=1977 0:5 7:8 0:6 6:5 0:8 4:9 1982=1984 1:2 3:3 0:45 8:7 0:8 4:9 1991=1992 0:4 9:7 0:8 4:9 0:75 5:2 1996=1997 0:6 6:5 0:6 6:5 1:5 2:6 Global values 0:35 11:3 0:45 8:7 0:48 8:1 Event date higher ones, the events previously reported still being signalled by dominant peaks. The Lorenz contamination lor slightly tapers the peaks or equivalently induces a false increase of the local predicting horizons. In any case we do not attach too much importance to the exact values of the local Lyapunov exponents or prediction horizons. Their precise estimation is quite sensitive to data contamination. On the other hand, their relative change with time is robust and can be correlated with external excitations of Earth rotation. It should also be kept in mind that the time-dependent smoothing of the combined EOP series (with a cut-o¡ period as long as 30 days for the polar motion and 15 days for LOD before 1980; see Paper I, Section 2) is likely to have seriously damped those events a¡ecting the very short periods of Earth rotation £uctuations in the 1970s. The strength and weakness of our very tentative interpretations is discussed in Section 6. 5 D ET E RM I NAT I O N OF T H E AT T R AC TOR D I ME N S I O N DA. In Paper I, the Grassberger & Procaccia (1983) algorithm gave us a mean to evaluate the embedding dimension DE , de¢ned as the ¢rst trial dimension from which the slopes of the log[C(r)] versus log(r) diagram begin to saturate, C(r) being the correlation integral. Here we are interested in ¢nding an attractor dimension for the EOP time-series. An attractor is a geometric ¢gure lying in the phase space of a dissipative dynamical system which captures all the orbits belonging to a given basin of initial conditions after transients fade out. Several competing attractors can coexist with very intricate basin boundaries (see e.g. Sommerer & Ott 1993). Their attracting nature is intimately related to the existence of stable manifolds along which any orbit perturbation will shrink at an exponential rate (see Section 4). Therefore, DA < DE . The local DL and Lyapunov DLyap dimensions found in Sections 2 and 4, respectively, are also related to the attractor ß 1999 RAS, GJI 137, 565^579 Non-linear processes in Earth rotation time-seriesöII dimension: DA ¦DL < DA z1 , (19) DA ¦DLyap . (20) The dimension of the attractor DA is a good characterization of the dynamical regime of the non-linear system. For instance, a white noise process has no attractor since the embedded points are randomly distributed in an in¢nite-dimensional phase space without exhibiting any particular structure. The attractor of a k-quasi-periodic system is a k-torus (with k integer), a smooth well-de¢ned geometric ¢gure in phase space. On the other hand, a strange attractor is a complicated ¢gure whose dimension is not an integer. In such a case the dimension provides a measure of the Cantor set structure of the attractor. In the past, most dynamical systems having this kind of attractor were found to be chaotic, in particular through the analysis of their global Lyapunov exponents, but strange nonchaotic attractors were also discovered (Ding et al. 1989). We must thus be very careful before concluding from the fractal dimension of an attractor the chaoticity of the system. The chaotic nature of the rapid £uctuations in Earth rotation has already been established (Section 4), thus we have no such dilemma. It has also been shown that a coloured noise might exhibit a non-integer dimension as determined by the correlation integral (Osborne & Provenzale 1989) because its reconstructed curve in pseudo-phase space is fractal (self-similar and not di¡erentiable). From a projection in 3-D space of the EOP trajectories (Paper I, Fig. 10), we saw that the orbits were smooth and formed a ¢gure which did not ¢ll the space randomly like a white noise nor fractally like a random walk. As mentioned previously, the presence of positive Lyapunov exponents reveals the existence of a strange attractor in phase space. The sum of the global Lyapunov exponents being negative further indicates that dissipative mechanisms are active in the rapid £uctuations of Earth rotation. For the EOP time-series, the dimension of the attractor must be greater than 3 because the orbits are still self-crossing in a 3-D space. DA cannot be greater than 5 with regard to the determinations of the embedding dimension, local dimension and the Lyapunov dimension made previously. In order to determine the attractor dimension DA , we use the Grassberger & Procaccia (1983) algorithm based on the computation of the correlation integral C(r) (see Paper I, Section 6 for details). When the embedded data vectors y(i) populate an attractor, the correlation integral (`spatial correlation of points'), C(r), behaves like a scaling law. Consider that the y(i)s are homogeneously distributed on a surface: the number of pairs of points which are at a distance of less than r then grows proportionally to r2 as r is increased. For an attractor of dimension DA , we have the same kind of behaviour: C(r) behaves like rDA for r going to zero and the number of data points N going to in¢nity. We see that working with ¢nite sets of experimental data will be an obstacle for the convergence of the C(r) slope to the attractor dimension. Obviously r cannot be taken too small because of the truncation and limited accuracy of the data, nor can N go to in¢nity; that is, the embedded data vectors y(i) poorly and inhomogeneously sample the attractor. In theory the way the correlation integral is constructed gives a relevant measure for the attractor dimension because it is sensitive to the dynamical process ß 1999 RAS, GJI 137, 565^579 575 of coverage of the attractor and in that sense it is di¡erent from the usual Haussdorf dimension DH , even if they are very close (DA ¦DH ) (Eckmann & Ruelle 1985). The Haussdorf dimension results from a coverage of the attractor by cubes of length l. The number nb(l) of cubes needed to cover the attractor is related to the dimension by nb(l)&l {D. We see that DH is only sensitive to the geometric structure of the attractor as a whole and not to the density of points. Plots of log[C(r)] versus log(r) have been drawn for the three EOP time-series (see e.g. Paper I, Fig. 6 for LOD). We then select a range of log(r) values where the log[C(r)] curve is almost linear (note that this range depends on the trial embedding dimension). The mean slope over this curve segment is an estimate of the attractor dimension and the slope rms about its mean is reported as an error bar in Fig. 7. The higher the trial embedding dimension, the narrower the range of admissible log(r) values. As a consequence, error bars tend to be larger for higher embeddings, and our uncertainty estimates are mainly representative of the way we pick the slopes, not really of the e¡ect of various error sources on the slope determinations. No such theoretical accuracy bounds have been published. There is no obvious convergence of the slope to a constant value around an embedding dimension DE ~5, thus it is not straightforward to select a precise attractor dimension. The increase in the slope is much the same for LOD (Fig. 7a), PMX (Fig. 7b) and PMY (Fig. 7c). It increases until the trial embedding dimension reaches 5 or 6 and then grows more and more slowly, as if a random process, superimposed on the deterministic signal, was partly counterbalancing the saturation e¡ect of dynamical origin by pulling up the curve. The slope quasi-saturation is more evident for the PMX data series (Fig. 7b). For comparison, the slope of the correlation integral for a white noise process always equals the trial embedding dimension (it is shown as a dashed diagonal line in Fig. 7). Assuming an embedding dimension of 5 or 6, from Fig. 7 we obtain attractor dimensions of 4.5^5.5 for LOD, 3.5^4.5 for PMX and 4^5 for PMY. These results verify the inequality constraints relating DA to the local and Lyapunov dimensions (eqs 19^20). We are not very con¢dent in the way that these determinations of DA are obtained. Because of the limited number of data, their quantization and the presence of noise contamination, the correlation integral algorithm is unable to provide robust results (see also Roberts 1991). Indeed, in Paper I the saturation of the C(r) slopes was clearly detected (Paper I, Section 6 and Fig. 7) for an embedding dimension between 5 and 7. This behaviour is in fact mainly observed for small radii r corresponding to the part of the log[C(r)] curve that shows a `knee' (Paper I, Fig. 8); our DA estimates must be performed in a higher range of r values corresponding to linear segments of the plot. The knee is partly generated by `false' counting of pairs whose points are neighbours in phase space E* only because they are consecutive in time. Discarding these pairs (Theiler 1986) is of no real help for short time-series because the sparse coverage of the attractor provides a very restricted number of `true' neighbours. This knee might also be a genuine signature of the underlying dynamics. In fact, in the tests performed in Section 6 of Paper I we saw that the knee disappeared when the EOP data series were contaminated by noise, but remained unchanged when low-dimensional chaos was added. Eckmann & Ruelle (1985) also noted that such a knee might result from the observed dynamical variable 576 V. Frede and P. Mazzega Figure 7. Slope of the correlation integral as a function of the trial embedding dimension DE (for the error bars evaluation see Section 5): (a) for LOD (solid line) and corresponding slope for a white noise (dashed line); (b) for PMX (solid line) and corresponding slope for a white noise (dashed line); (b) for PMY (solid line) and corresponding slope for a white noise (dashed line). being the superimposition of two non-interacting dynamical system outputs. Moreover, data given at constant time steps cannot homogeneously cover the attractor. The corresponding steps in E* will be stretched when approaching the most external boundaries of the attractor. The C(r) estimate will be even more penalized by the data scarcity in these regions of E*. Noise contamination e¡ects will be most often observed for small radii. For large radii r, C(r) is only slightly representative of the attractor geometry. Indeed, when increasing r, orbit segments in E* with distance ¦r will progressively lose their dynamical connection. The slope of C(r) is more and more meaningless for increasing radii. We could hardly draw conclusions about the detectability of non-linear deterministic processes in EOP from the application of this algorithm alone. The range of attractor dimensions we give allows both integer and non-integer values. The computation of the Lyapunov exponents behaves better than the correlation interval (such a remark is made by Ruelle 1990). We thus take the Lyapunov dimension for a more reliable characterization of the attractor strangeness of the EOP series. 6 D IS CUS SI O N Dimensional and dynamical characteristics extracted from the ¢ltered EOP data series are summarized in Table 1. We have checked that most of the results reported are robust against data contamination. Such topics and more speci¢c issues are discussed in each section with regard to the basic ß 1999 RAS, GJI 137, 565^579 Non-linear processes in Earth rotation time-seriesöII geometrical or statistical principle used to design the corresponding algorithm. Globally speaking, however, we are faced with the same di¤culties throughout this study. The non-linear time-series analysis performed here relies on concepts of the dynamical system theory as well as on the possibility of obtaining information on the underlying source dynamics by embedding the data time-series in a reconstruction space. To operate these concepts rigorously, making them valid for characterizing asymptotic time evolution, one would need in¢nite uncontaminated time-series of measurements. Working with ¢nite sets of noisy data has been a recurrent problem throughout this study. Several algorithms designed and tested so as to give accurate results with long series of clean synthetic data provided mitigated results. For example, it appears that the dimensional characteristics (embedding, local and attractor dimensions) are overestimated when dealing with EOP series contaminated with highdimensional stochastic processes. Contamination by a lowdimensional chaotic process does not signi¢cantly change the determinations. Other assumptions are implicit in this approach. We mentioned that several attractors can coexist in a phase space. Due to external forcings or to noise, the system orbit may switch from one attractor to another in an apparently erratic manner. In the presence of a single attractor, the system can also experience bifurcations induced by time variations of some parameter controlling its dynamical regime. The characteristics obtained from the analysis of the corresponding time-series will probably result from some sort of averaging process. No systematic study has been conducted to ¢nd out the e¡ects of such qualitative change on the characteristics of the observed system. The ocean^atmosphere coupled system is complex enough to allow for the coexistence of competing attractors and certainly obeys bifurcation sequences. For example, the famous El Nin¬o is probably in subharmonic resonance with the atmospheric annual cycle (Jin et al. 1994). A possible scenario for its unpredictable evolution is a succession of jumps from one frequency-locked mode to another with possible transitions to chaos (Tziperman et al. 1994, 1995). Among the characteristics of the EOP we obtained, only the local Lyapunov exponents could detect such a time evolution of the dynamics. Relative losses in Earth rotation stability, with di¡erent strengths, are clearly related to El Nin¬o or La Nin¬a events (Section 4). We found two positive global Lyapunov exponents associated with each EOP time-series embedded in a DE ~5 phase space. These are indicators of the chaotic nature of the rapid £uctuation in Earth rotation in the years 1970^1997. Dissipation is also active, as is shown by the negative sum of the exponents. Because of the presence of stable manifolds in the pseudo-phase space E*, the orbit reconstructed from the EOP data series will not occupy a full 5-D volume, but will be con¢ned to an attractor of lower dimension DA . The Lyapunov dimension DLyap is an upper bound for DA. From the respective Lyapunov spectra we obtained DLyap (LOD)~4.48, DLyap (PMX)~4.90 and DLyap (PMY)~4.97. The attractor dimensions found from the correlation integral algorithm are broadly similar [DA (LOD)~4.5^5.5, DA (PMX)~3.5^4.5 and DA (PMY)~4^5]. Large uncertainties in their determination do not allow one to ascertain their non-integer nature from this computation only. Once again, these exponents and dimensions are only representative of the average instability of the ß 1999 RAS, GJI 137, 565^579 577 system, whilst several oceanic and atmospheric processes with di¡erent time- and space scales are superimposed to excite the observed £uctuations in EOP. Two studies have been undertaken to extract the dimension of an assumed attractor from LOD time-series. They both rely on the correlation integral approach as proposed by Grassberger & Procaccia (1983). From the long-period variations in LOD, Malinetskii et al. (1990) found a lower bound for the attractor dimension of 1:3+0:3. They used 6 month averaged data from 1656 to 1984. Here the underlying source is the magnetohydrodynamical processes in the Earth's innermost regions. Studying the rapid £uctuations, Tiwari et al. (1992) applied the approach to &5300 daily observations of LOD from 1976 to 1990. They found an attractor dimension between 5 and 7 and concluded chaoticity of the LOD variations. The reliability of this result obtained with the same algorithm, which appears unsatisfactory to us, is not discussed. Nevertheless, both the dimension and the associated large uncertainty are in agreement with our results in Section 5. Tiwari et al. (1994) also estimated the production rate of Kolmogorov entropy and found an intrinsic timescale of 7^8 days. This value is comparable to the horizon of prediction HP (LOD)~11.3 days that we estimated from the principal Lyapunov exponent j1 (LOD). Note that the way we convert j1 (LOD) into a prediction horizon is, to some extent, arbitrary and is intended to meet some operational requirement (HP is the time necessary for a 1 per cent rms error to grow to 50 per cent rms). Several similar studies were carried out for various timescales of the atmosphere dynamics. Atmospheric £ows are the dominant source of excitation for the rapid £uctuations of Earth rotation considered in this study. Attractors corresponding to hourly (Tsonis & Elsner 1988; DA ~7.3), seasonal (Essex et al. 1987; Fraedrich 1987; DA between 6 and 7) and climatological (Nicolis & Nicolis 1984; DA ~3.1) processes have been detected from the correlation integral analysis of data time-series (velocity winds, isotope records of deep sea cores, etc.). Some of these results were disputed because of the shortness of the data records (see e.g. Eckmann & Ruelle 1992). Nonetheless, chaos is probably ubiquitous over a wide range of time- and space scales in the £ows forcing Earth rotation, with di¡erent dimensional and dynamical characteristics portraying each mechanism at work. The Lyapunov dimension we obtain from the ¢ltered EOP series might be a signature of the averaged weather dynamical regime. Some of the relative losses of Earth's rotational stability are clearly related to La Nin¬a and El Nin¬o events. Unfortunately, no dimensional and dynamical characteristics of these major phenomena have been published. Stability analyses attempt to provide explanations for its pathological irregularity based on dynamical grounds (Vallis 1986; Jin et al. 1994; Tziperman et al. 1994, 1995). No long and dense enough data timeseries are available to undertake a non-linear analysis based on embeddings and compare the resulting El Nin¬o characterizations with those deduced from the EOP time-series. The impact of El Nin¬o on the LOD was observed soon after the strong 1982^1983 event (Rosen et al. 1984; see also Chao 1989). The associated atmospheric angular momentum as deduced from operational weather models accounts for more than 90 per cent of the LOD interannual variability during these years (Dickey et al. 1994). We should stress that our analysis sheds light on the causality relationship between 578 V. Frede and P. Mazzega ENSO (El Nin¬o/Southern Oscillation) events and £uctuations in Earth rotation from a completely di¡erent point of view. Indeed, the local Lyapunov exponents drawn in Fig. 4 for LOD and Figs 5 and 6 for PMX and PMY measure to some extent temporary chaoticity increases caused by an ENSO. A naive approach would be to try to ¢nd a systematic `response' in the local exponents to a known excitation. In fact, we see that an ENSO event is able to cause a deterioration in the stability of one or two EOP time-series without a¡ecting the other one (compare the e¡ect of the 1982^1983 El Nin¬o on the j1 (t)s for LOD and PMY with that for PMX). Moreover, the local stability of a system may depend on the system state so that the same external excitation will induce di¡erent perturbations in its time evolution. Such issues must be carefully investigated for the Earth's rotational system. All the EOP perturbations studied here concern the rapid £uctuations of Earth rotation. The signature of events such as El Nin¬o is seen in this period range, from a few days to 100 days, but lasts as long as the excitation is maintained. The data pre-processing is an important step that can partially exclude some signal components. In particular, the heavy smoothing performed in the early years of the IERS combined series has probably damped out some `local' events of the Earth's rotation dynamics, without compromising the global reconstruction, dimensional and stability analysis of the embedded orbits. Longer and less smoothed EOP time-series with a much higher sampling rate (not arti¢cially increased by data interpolations) would be highly valuable in order to obtain more accurate dimensional and dynamical characteristics of the reconstructed orbit. With a better coverage of the Earth rotation attractor we should be able to improve the resolving capability of the local Lyapunov exponents and to relate the corresponding stability variations to weaker or shorter excitation processes. 7 C O NCLUS IO N S For the ¢rst time, a complete set of concepts from the dynamical system theory has been exploited to characterize the dynamical behaviour of Earth rotation for periods under 100 days. Dimensional (Paper I) and dynamical (this paper) characteristics of the underlying data source have been extracted from three ¢ltered EOP time-series of daily observations spanning 27 years, using their embedding in a pseudophase space E*. The estimates of these characteristics are summarized in Table 1. We have checked that these values are quite robust against noise contamination of the time-series. However, some degree of uncertainty is associated with these determinations, mainly resulting from the fact that the data sets are limited (about 104 samples in each series). The orbits have been reconstructed from data vectors in time-delayed coordinates. The optimal delays (10 days for LOD, 15 days for PMX and 18 days for PMY) were obtained from the function of average mutual information (Paper I, Section 4). A 5- or 6-D space E* is shown to be adequate for unfolding the reconstructed orbits (Paper I, Sections 5 and 6). The three orbits corresponding to LOD and to the two components of polar motion traced smooth but complicated trajectories in E* that seem to have been generated from a di¡erentiable dynamical system. Their behaviour is in strong contrast with the erratic and fractal curves reconstructed from white noise and random walk time-series respectively (Paper I, Fig. 10). Geometric and statistical analyses of these orbits provided further characterization of the EOP time-series. We obtained several consistent clues establishing the chaotic nature of the rapid £uctuations in Earth rotation. A low-dimensional deterministic process is at work that can be captured with a set of DL ~5 active dynamical variables. Two positive Lyapunov exponents are associated with each EOP series. The corresponding intrinsic horizons of prediction are of 11.3 days for LOD, 8.7 days for PMX and 8:1 days for PMY. Beyond these time intervals, any initial error of 1 per cent of the series rms will exceed the 50 per cent rms level. Both this sensitivity of the system time evolution to its initial state and the dissipative nature of the underlying dynamical process (as demonstrated by the negative sum of the Lyapunov exponents) imply the presence of a strange attractor in the phase space of each series. The Lyapunov dimension DLyap is 4.48 for LOD, 4.90 for PMX and 4.97 for PMY. As an upper bound for the attractor dimension, it agrees with the dimensions DA obtained from the correlation integrals. Moreover, it was shown that major relative losses in Earth's rotation stability may be associated with La Nin¬a and El Nin¬o events. Such perturbations in the system stability have their signature in the time-series of the local Lyapunov exponents (Figs 4^6). The local horizons of prediction are correspondingly reduced and can sometimes reach values as low as 2.5^5 days (Table 2). Not all of the three EOP stabilities are equally degraded by external excitations. At the moment some observed strong stability losses are unexplained (see the 1996 peak in the local Lyapunov exponents of the PMY timeseries). A detailed stability analysis of the Earth's rotational system dynamics is required. Low-dimensional deterministic processes are detectable in the time-series of Earth rotation. AC K NOW L E D G ME N T S All the computations were performed on the CRAY YMP at the Centre National d'Etudes Spatiales (Toulouse). This study was supported by the Centre National de la Recherche Scienti¢que under contract 2E0025 (programme Mode©lisation et Simulation Numërique). R E F ER E NC E S Abarbanel, H.D.I., 1996. Analysis of Observed Chaotic Data, Springer Verlag, Berlin. Abarbanel, H.D.I. & Kennel, M.B., 1993. Local false nearest neighbors and dynamical dimensions from observed chaotic data, Phys. Rev. E, 47, 3057^3068. Abarbanel, H.D.I., Brown, R. & Kennel, M.B., 1991. Variation of Lyapunov exponents on a strange attractor, J. Nonlinear Sci., 1, 175^199. Brown, R., Bryant, P. & Abarbanel, H.D.I., 1991. Computing the Lyapunov spectrum of a dynamical system from an observed time series, Phys. Rev. A, 43, 2787^2806. Bryant, P., Brown, R. & Abarbanel, H.D.I., 1990. Lyapunov exponents from observed time series, Phys. Rev. Lett., 65, 1523^1526. Chao, B.F., 1989. Length of day variations caused by El Nin¬o -southern oscillation and quasi-biennal oscillation, Science, 243, 923^925. ß 1999 RAS, GJI 137, 565^579 Non-linear processes in Earth rotation time-seriesöII Chao, B.F. & Gross, R.S., 1987. Changes in the Earth's rotation and low degree gravitational ¢eld induced by earthquakes, Geophys. J. R. astr. Soc., 91, 569^596. Dickey, J.O., Marcus, S.L., Hide, R., Eubanks, T.M. & Boggs, D.H., 1994. Angular momentum exchange among the solid Earth, atmosphere and ocean: a case study of the 1982^1983 El Nin¬o event, J. geophys. Res., 99(B12), 23 921^23 937. Ding, M., Grebogi, C. & Ott, E., 1989. Evolution of attractors in quasiperiodically forced systems: from quasiperiodic to strange nonchaotic to chaotic, Phys. Rev. A, 39, 2593^2598. Eckmann, J.P. & Ruelle, D., 1985. Ergodic theory of chaos and strange attractors, Rev. Mod. Phys., 57, 617^656. Eckmann, J.P. & Ruelle, D., 1992. Fundamental limitations for estimating dimensions and Lyapunov exponents in dynamical systems, Physica D, 56, 185^187. Eckmann, J.P., Kamphorst, S.O., Ruelle, D. & Ciliberto, S., 1986. Liapunov exponents from time series, Phys. Rev. A, 34, 4971^4979. Engdahl, E.R., Van der Hilst, R. & Buland, R.P., 1998. Global teleseismic earthquake relocation with improved travel times and procedures for depth determinations, Bull. seism. Soc. Am., submitted. Essex, C., Lookman, T. & Nerenberg, M.A.H., 1987. The climate attractor over short timescales, Nature, 326, 64^66. Fraedrich, K., 1987. Estimating weather and climate predictability on attractors, J. atmos. Sci., 44, 722^728. Fraser, A.M. & Swinney, H.L., 1986. Independent coordinates for strange attractors from mutual information, Phys. Rev. A, 33, 1134^1140. Frede, V. & Mazzega, P., 1999. Detectability of deterministic non-linear processes in Earth rotation time-seriesöI. Embedding, Geophys. J. Int., 137, 551^564 (this issue). Gambis, D., 1996. Multi-technique EOP combinaisons, Proc. IGS Analysis Center Workshop, Silver Spring, Vol. 96, eds Van Scoy, P. & Neylan, R.E., Pasadena, JPL Publ. Golub, G.H. & van Loan, C.F., 1989. Matrix Computations, 2nd edn, John Hopkins University Press, London. Grassberger, P. & Procaccia, I., 1983. Characterization of strange attractors, Phys. Rev. Lett., 50, 346^349. Jin, F.F., Neelin, J.D. & Ghil, M., 1994. El Nin¬o on the Devil's staircase: annual subharmonic steps to chaos, Science, 264, 70^72. Kaplan, J.L. & Yorke, J.A., 1979. Chaotic behaviour in multidimensional di¡erence equations, in Lecture Notes in Mathematics, Vol. 730, pp. 228^237, Springer, Berlin. Malinetskii, G.G., Potatov, A.B., Gizzatulinan, S.M., Ruzmaikin, A.A. & Rukavishnikov, V.D., 1990. Dimension of geomagnetic attractor from data on length of day variations, Phys. Earth planet. Inter., 59, 170^181. Nicolis, C. & Nicolis, G., 1984. Is there a climatic attractor?, Nature, 311, 529^532. ß 1999 RAS, GJI 137, 565^579 579 Osborne, A.R. & Provenzale, P., 1989. Finite correlation dimension for stochastic systems with power-law spectra, Physica D, 35, 357^381. Palmer, T.N., 1993. Extended-range atmospheric prediction and the Lorenz model, Bull. Am. Met. Soc., 74, 49^65. Palmer, T.N., Buiza, R., Molteni, N., Chen, Y.C. & Corti, S., 1994. Singular vectors and the predictability of weather and climates, Phil. Trans. R. Soc. Lond., A348, 459^475. Roberts, D.A., 1991. Is there a strange attractor in the magnetosphere?, J. geophys. Res., 96(A9), 16 031^16 046. Rosen, R.D., Salstein, D.A., Eubanks, T.M., Dickey, J.O. & Steppe, J.A., 1984. An El Nin¬o signal in atmospheric angular momentum and Earth rotation, Science, 225, 411^414. Ruelle, D., 1990. Deterministic chaos: the science and the ¢ction, Proc. R. Soc. Lond., A427, 241^248. Sauer, T., Yorke, J.A. & Casdagli, M., 1991. Embedology, J. Stat. Phys., 65, 579^616. Sommerer, J.C. & Ott, E., 1993. A physical system with qualitatively uncertain dynamics, Nature, 365, 138^140. Takens, F., 1981. Detecting strange attractors in turbulence, in Dynamical Systems and Turbulence, pp. 366^381, eds Rand, D. & Young, L.S., Lecture Notes in Mathematics 898, Springer Verlag, Berlin. Tanaka, T., Aihara, K. & Taki, M., 1996. Lyapunov exponents of random time series, Phys. Rev. E, 54, 2122^2124. Theiler, J., 1986. Spurious dimension from correlation algorithms applied to limited time-series data, Phys. Rev. A, 34, 2427^2432. Thompson, J.M.T. & Stewart, H.B., 1987. Nonlinear Dynamics and Chaos, J. Wiley & Sons, New York. Tiwari, R.K., Negi, J.G. & Rao, K.N.N., 1992. Attractor dimensions in non-linear £uctuations of length of the day (LOD) variations, Geophys. Res. Lett., 19, 909^912. Tiwari, R.K., Negi, J.G. & Rao, K.N.N., 1994. Strange attractor in nonlinear £uctuations of length of the day (LOD) time series, in Nonlinear Dynamics and Predictability of Geophysical Phenomena, Geophys. Monogr., 83, IUGG Vol. 18, pp. 61^67. Tsonis, A.A. & Elsner, J.B., 1988. The weather attractor over very short timescales, Nature, 333, 545^547. Tu¢llaro, N.B., Abbott, T. & Reilly, J., 1992. An Experimental Approach to Nonlinear Dynamics and Chaos, Addison Wesley, New York.zannexes. Tziperman, E., Stone, L., Cane, M.A. & Jarosh, H., 1994. El Nin¬o chaos: overlapping of resonances between the seasonal cycle and the Paci¢c ocean-atmosphere oscillator, Science, 264, 72^74. Tziperman, E., Cane, M.A. & Zebiak, S.E., 1995. Irregularity and locking to the seasonal cycle in an ENSO prediction model as explained by the quasi-periodicity route to chaos, J. atmos. Sci., 52, 293^306. Vallis, G.K., 1986. El Nin¬o: a chaotic dynamical system?, Science, 232, 243^245.