Lorito.pdf 573 ANNALS OF GEOPHYSICS, VOL. 46, N. 3, June 2003 Wavelet analysis on paleomagnetic (and computer simulated) VGP time series Stefano Lorito (1) (2), Grazia Giberti (1), Agata Siniscalchi (3) and Marina Iorio (4) (1) Dipartimento di Scienze Fisiche, Università «Federico II», Napoli, Italy (2) Istituto Nazionale di Geofi sica e Vulcanologia, Roma, Italy (3) Dipartimento di Geologia e Geofi sica,Università di Bari, Italy (4) Istituto IAMC, Geomare Sud, CNR, Napoli, Italy Abstract We present Continuous Wavelet Transform (CWT) data analysis of Virtual Geomagnetic Pole (VGP) latitude time series. The analyzed time series are sedimentary paleomagnetic and geodynamo simulated data. Two mother wavelets (the Morlet function and the first derivative of a Gaussian function) are used in order to detect features related to the spectral content as well as polarity excursions and reversals. By means of the Morlet wavelet, we estimate both the global spectrum and the time evolution of the spectral content of the paleomagnetic data series. Some peaks corresponding to the orbital components are revealed by the spectra and the local analysis helped disclose their statistical significance. Even if this feature could be an indication of orbital influence on geodynamo, other interpretations are possible. In particular, we note a correspondence of local spectral peaks with the appearance of the excursions in the series. The comparison among the paleomagnetic and simulated spectra shows a similarity in the high frequency region indicating that their degree of regularity is analogous. By means of Gaussian first derivative wavelet, reversals and excursions of polarity were sought. The analysis was performed first on the simulated data, to have a guide in understanding the features present in the more complex paleomagnetic data. Various excursions and reversals have been identified, despite of the prevalent normality of the series and its inherent noise. The found relative chronology of the paleomagnetic data reversals was compared with a coeval global polarity time scale (Channel et al., 1995). The relative lengths of polarity stability intervals are found similar, but a general shift appears between the two scales, that could be due to the datation uncertainties of the Hauterivian/Barremian boundary. Mailing address: Dr. Stefano Lorito, Dipartimento di Scienze Fisiche, Università «Federico II», Complesso Universitario Monte Sant’Angelo, Via Cintia, 80126 Napoli, Italy; e-mail: lorito@na.infn.it Key words CWT – VGP – geodynamo simulations – excursions and reversals – Earth’s orbital parameters 1. Introduction The long term phenomenology of the Earth’s magnetic field can be inferred from paleomagnetic data as they record the field behaviour on time scales which go from hundreds to millions of years. In particular, sedimentary data can achieve the time resolution to allow the extraction of details of the past geomagnetic field when the estimated sedimentation rate is stable and sufficiently high to affirm that samples pre- serve the imprinting of their first magnetization (Merrill and McFadden, 1999). A drawback of these data is the complexity of the formation/ magnetization process of sedimentary rocks and of their subsequent history. The measured magnetic signal could be thus have been fil- tered in a not well understood way (Merrill and McFadden, 1999). The need for accurate and non-standard analyses of the available data with all refined techniques therefore emerges (Dormy et al., 2000; McMillan et al., 574 Stefano Lorito, Grazia Giberti, Agata Siniscalchi and Marina Iorio 2001). Here we utilize the Continuous Wavelet Transform (hereafter CWT) together with the Discrete Fourier Transform (DFT) to analyze paleomagnetic sedimentary data. CWT is, in fact, a versatile analysis’ technique (Mallat, 1998) which, in principle, could investigate and discriminate among regular/singular and stationary/non-stationary data features. In this work, both paleomagnetic sedimentary (Iorio, 1995) and computer simulated (Glatzmaier et al., 1999) Virtual Geomagnetic Poles (VGP) latitude time series are analyzed. One aim is to explore how useful the CWT can be in estimating both the power spectrum of the whole data series than the time evolution of the spectral content. Similar studies, using CWT, have been performed by some authors on geomagnetic paleointen- sity data (Guyodo et al., 2000; Yokoyama and Yamazaky, 2000). This tool is recognized as sui- table to analyze non-stationary signals (Mallat, 1998; Torrence and Compo, 1998 and references therein), the periodicity over the whole series not being imposed a priori, as in every DFT analysis. Moreover, local spectral features result better characterized. The CWT uses windows that are shorter in time (and longer in frequency) the higher the frequency analyzed, allowing for a higher joint time-frequency resolution with respect to the windowed DFT (Mallat, 1998). To assess the reliability of the CWT analysis, the usual DFT is performed on all the data. In fact, the accuracy in estimating the decay of the power spectrum at high frequencies by means of CWT (Perrier et al., 1995) depends on the kind of Fourier spectrum and in turn on the nature of the signal. The CWT is also utilized to look for local regularity versus sharp variations in the latitude time series with the aim of localizing and cha- racterizing reversals and polarity excursions. A similar approach was followed by Alexandrescu et al. (1995) in studying geomagnetic jerks of the present geomagnetic field, there defined as singularities in the second derivative of the signal. The definitions of polarity reversals and excursions are a matter of debate (Merrill et al., 1996; Gubbins, 1999) both with regard to the amplitude and duration/stability of these processes. It is, however, generally accepted that reversals are fast transitions between oppo- site polarities characterized by a long polarity stability before and after the transition whereas the excursions would present instead a sudden return to the original polarity. We will detail later the definitions we use to quantify local regularities with the appropriate CWT. The efficacy of wavelet analysis depends on the particular mother wavelet employed, so the choice is usually made on the basis of the specific features one wants to investigate and quantify. We utilize two kinds of mother wavelets: i) the Morlet (Morlet et al., 1982), to ana- lyze local and global spectral content of the series; ii) the first derivative of Gaussian (1DOG), to localize and study reversals and polarity ex- cursions. This paleomagnetic data-set has already been analyzed by standard spectral and non spectral methods (Iorio et al., 1995, 1998, 2001; Tarling et al., 1999). In the framework of the time- frequency analysis (Mallat, 1998), this CWT study addresses a quantitative verification of some previously found results as: i) the influence of the main Milankovitch periods (Berger et al., 1992) on the paleomagnetic signal; ii) similitudes between observed and simu- lated data. 2. Data description 2.1. Paleomagnetic data series The paleomagnetic data come from bore- holes carried out at Mt. Raggeto (41.1°N, 14.2°E), Italian Southern Apennines. Rock samples were extracted from a carbonate sequence formed in the Early Cretaceous by deposition of shallow water sediments (Iorio et al., 1995, 1998). They constitute a composite core, with an overall length of 88 m. The (paleo)magnetic remanence intensity and directions were measured every 2 cm along the cores, so obtaining a sequence that spans an interval of approximately 3.2 Myr, when assuming that, on average, 2 cm were deposited in 720 yr (Tarling et al., 1999). As the sequence contains the Hauterivian-Barremian boundary (127 ± 1.6 Myr before present, Iorio et al., 1998), it was inferred that the sediments analyzed 575 Wavelet analysis on paleomagnetic (and computer simulated) VGP time series were deposited in the time interval included in the range 124-128 Myr before present. This geological age collocates a few million years before the so-called Cretaceous Superchron, which is generally considered a very long (about 35 Myr) period of normal polarity (Merrill et al., 1996). Here and in the following, normal pola- rity is defined as the present polarity of the geo- magnetic field. The samples were analyzed from a lithologic, petrological, textural and paleontological point of view (Iorio, 1995; D’Argenio et al., 1999). After these analyses, some intervals were excluded as corrupted by recent karst phenomena. Then, the longest evenly spaced data set is 2.3 Myr long. The estimated velocity of sedimentation of 1/360 cm yr-1 renders these data useful to study polarity reversals, if it is accepted that a transition occurs in a time interval ranging from 1 to 10 kyr. The stability of the sedimentation process, together with the relatively fast diagenesis of the samples (D’Argenio et al., 1999) and the above mentioned value of the sedimentation rate allow us to infer that the magnetization characteristics can be reliable for this kind of analysis (Merrill and McFadden, 1999). The VGP’s positions with respect to the mean paleomagnetic pole were calculated by means of standard paleomagnetic techniques (Tarling et al., 1999). Figure 1 presents the latitude of the longest evenly spaced sequence versus relative time. It corresponds to the data segment analyzed here. In this plot, when the latitude value is 90°, the VGP coincides with the mean paleomagnetic pole. The sequence shows a very high rate of variability and a prevailing normal polarity is evident. Nevertheless, several reversed chrones and excursions were identified and correlated with those in other pelagic Lower Cretaceous magnetostratigraphic sequences (Iorio et al., 1998). 2.2. Simulated data series The simulated data here analyzed (fig. 2) were produced by the geodynamo numerical simulation to which the authors refer as the tomographic case (Glatzmaier et al., 1999). This indicates that, as a condition for the heat flow at Core Mantle Boundary (CMB), a pattern was used that includes lateral variations as inferred by the interpretation of seismic tomographic results. The total heat flow through CMB was chosen equal to 7.2 TW. The equations were integrated using 3D spherical geometry and the simulation was carried ahead for a total time of 320 000 «model» years. This total time is estimated by the authors (Coe et al., 2000) long enough in order to assert that the field is self-sustained, when it is compared with the free dipole decay time (roughly 20 kyr). The results of the simulation are in the form of Schmidt quasi normalized Fig. 1. Mt. Raggeto VGP latitude versus relative time (Myr). The instant t = 0 Myr is assumed coincident with the beginning of the analyzed data segment. 576 Stefano Lorito, Grazia Giberti, Agata Siniscalchi and Marina Iorio Gauss coefficients, up to degree 21, of the scalar magnetic potential. The obtained magnetic field space spectrum is in good agreement in shape, but with a lower intensity, with respect to the spectrum based on Magsat data of the modern Earth’s magnetic field (Roberts and Glatzmaier, 2000). The derived VGP latitudes here used have a time step of 100 «model» years, considered in this paper as real years. It appears that the field reversed twice, the first reversal being faster, whereas in the second reversal the VGP latitude goes up to the normal polarity with a few minor excursions (Coe et al., 2000). Moreover, before of the first inversion two major excursions occur and, in between the two reversals, some major excursions or «rebounds» appear. 3. Wavelet transform The basics of the CWT are described in a number of papers and books (Mallat, 1998; Torrence and Compo, 1998 and references there- in). Here only the aspects of the CWT that are relevant to the analyses are detailed. The CWT of a function f Œ L2 ( ¬ ), is defined as (3.1) where (3.2) is the family obtained by scaling by s and trans- lating by u the argument of the mother wavelet y (t) Œ L1 (¬) I L2 (¬).The values Wf (u,s) of the CWT for a given pair (u, s) are called the wavelet coeffi cients. The mother wavelet must satisfy the admis- sibility condition (3.3) where ŷ (w) is the Fourier transform of y (t). This condition ensures the signal energy is conserved in the time-scale domain, i.e. that . (3.4) Different wavelet functions can be used and the choice is usually made on the basis of the specific features in the data set that one wants to investigate and quantify (Mallat, 1998). Here two kinds of mother wavelets are used: i) The Morlet wavelet, to estimate the local and global frequency content of the series; this choice is classical in harmonic time-frequency Wf u s f t t du s, ,( ) = ( ) ( )−∞ +∞ ∗∫ ψ ψ ψu s t s t u s u s s, , ,( ) = −    ∀ ∈ℜ > 1 0 C dψ ψ ω ω ω= ( ) < + ∞ +∞ ∫ ˆ 2 0 f t dt C Wf u s du ds s ( ) = ( ) −∞ +∞ −∞ +∞+∞ ∫ ∫∫ 2 2 20 1 ψ , Fig. 2. VGP latitude of the numerical geodynamo simulation (Glatzmaier et al., 1999) versus «model» time (kyr). 577 Wavelet analysis on paleomagnetic (and computer simulated) VGP time series sented by the terms . (3.8) The index i runs on the possible different sin- gularities present at each time t i and a i and b i are real constants. c(t) models a superposition of harmonic components and n(t) a Gaussian normal distributed noise. The chosen analyzing wavelet can detect singularities for which a i £ 1 (Mallat, 1998), the number of its vanishing moments being equal to 1. The quantities a i (Lipschitz exponents) and the t i values (corresponding to translations u i in the CWT) for the singular terms defined by eq. (3.8) are found by identification of the connected lines of local extrema (ridges) of the wavelets coefficients as functions of the scale s. It can be shown (Alexandrescu et al., 1995; Mallat, 1998) that the instants to which the ridges converge, as the scale decreases, give the localization in time t i of the terms in eq. (3.8). The local degree of regularity a i of these terms is found by analyzing the behaviour of the absolute values of the ridge functions, i.e. the local maxima of |Wf (u, s)|, versus the scale s in a neighborhood of the translations u centered around t i . It can be shown, in fact (Mallat, 1998) that, at fine scales, log , ~ log .Wf u s si( ) +( )α 1 2 As the above technique is applied here to localize and study reversals and polarity excursions in the VGP latitude time series and in view of the above mentioned debate on their definitions, to infer the analytical behaviour of these processes we adopt a working hypothesis. A reversal is considered a transition between the opposite polarities occurring in a few thousands years and characterized by a polarity stability before and after the transition for at least 20-50 kyr. For excursions, we assume a change in the latitude of the pole greater than 40-50 degrees, followed by a return to the original value in a time interval shorter than that occurring between a pair of reversals (Guyodo and Valet, 1999). As a consequence, we expect that when a reversal occurs the 1DOG will localize a ridge function with a Lipschitz coefficient ranging from 0 (which characterize a noiseless Heaviside function) to 1 (characterizing a noiseless ramp function). When analyses as the Morlet is sine-like shaped and limited in time (Morlet et al., 1982; Daubechies, 1992; Torrence and Compo, 1998). ii) The first derivative of a Gaussian (1DOG) to localize regular and singular behaviour in the signal (Mallat and Hwang, 1992), as could be considered reversals and polarity excursions. The normalized Morlet wavelet is defined as (3.5) where w0 is the central frequency. The presence of the correction term is necessary so that the admissibility condition holds (Mallat, 1998) and it is not specified in eq. (3.5) as numerically negligible for the choice w0 = 6 (Torrence and Compo, 1998) here adopted. The Morlet wavelet is employed to evaluate: i) The local wavelet power spectrum, i.e. the pattern of the Morlet squared coefficients |Wf (u, s)|2 in the time-scale plane (u, s). ii) The global wavelet power spectrum, i.e. the sum of the Morlet squared coefficients over the translations, that defines a consistent and unbiased estimator of the power spectrum (Torrence and Compo, 1998; Percival, 1995). iii) The time evolution of the squared coef- ficients over a chosen limited interval of scales. The last item allows us to analyze the energy content versus time of a given period in the signal as a definite factor of proportionality between the scale s and the Fourier equivalent period T exists for each wavelet (Torrence and Com- po, 1998 and references therein). In our case T = 1.03 s Dt, where Dt is the sampling rate of the data series. To estimate the local degree of regularity of a signal, we utilize the normalized 1DOG wavelet, i.e. the function (3.6) Following Alexandrescu et al. (1995), we assume here that the signal can be modeled as (3.7) where the searched singular behaviour is repre- ψ π ωt e ei t t( ) = +− −1 4 20 2 / correction term ψ π t d dt e t ( ) =     −2 1 2 2 2 . f t t t c t n ti i i i( ) = −( ) + ( ) + ( )+∑ β α t t t t t t t t i i i i i i −( ) = −( ) ≥ <    + α α for for 0 578 Stefano Lorito, Grazia Giberti, Agata Siniscalchi and Marina Iorio an excursion occurs, the 1DOG should localize two nearby (in time) ridges and the ridge function should decrease sharply at high scales because of the possible smaller amplitude of the process and/or shorter duration with respect to reversals. Ideally, an excursion in the signal should be seen as a bell function of small width which can possibly collapse into a d -function. 4. Morlet results The local wavelet spectrum obtained with the Morlet wavelet is shown in fig. 3, plotted as a 2D color image versus the scale s (converted in periods) and the translation u (times in our plot). The series has been padded with the necessary zeroes to limit, in the time-scale plane, the region of influence of the edges (delimited by the bold black line in the plot). This procedure also allows to decrease the computational time of the CWT by using a fast Fourier transform algorithm (Torrence and Compo, 1998). Local energy maxima appear at periods corresponding to the main orbital periods estimated for the Cretaceous, i.e. 22 (precession), 40 (obliquity), 100 (short eccentricity) and 400 (long eccentricity) kyr (Berger et al., 1992). Figure 4 shows the obtained global wavelet power spectrum in comparison with the Scargle periodogram (Scargle, 1982). Both are plotted in units of the latitude series variance s 2 = 1442 deg2 (standard deviation s � 38 deg). The Scargle periodogram is obtained padding the data series with zeroes as before. The white noise level at 90% for the Scargle periodogram is also plotted. The two curves give similar results, the wavelet spectrum being less scattered than its Fourier counterpart, due to the smoothing introduced by the Gaussian exponential present in the Morlet wavelet. It appears that the possible bias at high frequencies that the smoothing by the Gaussian exponential could introduce (Perrier et al., 1995), does not affect the calculated Morlet spectrum. In fig. 5, the red noise theoretical spectrum and the confidence levels at 50% and 90% are superimposed to the Morlet global periodogram. The red noise confidence level is calculated as- suming as noise-model a univariate lag-1 auto- regressive process (Torrence and Compo, 1998). The first order autocorrelation, after the sub- traction of the mean value, is equal to 0.88. With the exception of the peak related to orbital ob- liquity (� 40 kyr), the periods corresponding to the Cretaceous main orbitals emerge only over the 50% confidence level, and so they all are statistically not well resolved against the red noise background. To address the issue of distinguishing the peaks in the power spectrum from the oscillations inherent to a stochastic process, we analyze the time changes of the spectral content in bands around the equivalent periods where peaks Fig. 3. Local Morlet spectrum of Mt. Raggeto latitude (in units of time series variance) in the time-periods plane (see text for more comments). 579 Wavelet analysis on paleomagnetic (and computer simulated) VGP time series Fig. 4. Global Morlet spectrum and Scargle periodograms (in units of variance) of the Mt. Raggeto latitude. The 90% confidence level for the Scargle periodogram is represented by a dashed line. Fig. 5. Global Morlet spectrum of the Mt. Raggeto latitude, along with 90% and 50% red noise confidence levels. Spectral peaks corresponding to Cretaceous main orbitals are evidenced (see text for more comments). 580 Stefano Lorito, Grazia Giberti, Agata Siniscalchi and Marina Iorio emerge. To do that we follow Torrence and Com- po (1998), performing the weighted sum of the wavelets coefficients over a limited range of scales (equivalent periods) around three chosen periods (100, 40 and 22 kyr). The results for the three intervals selected are shown in fig. 6a-c. The per-interval red noise 50% and 90% confidence levels are plotted as horizontal dashed lines. For the scale interval including the short eccentricity orbital period (�100 kyr), we observe (fig. 6a) that the spectral energy emerges twice above the 90% confidence level, whereas all the other peaks are well above the 50% level. The long eccentricity orbital period (� 400 kyr) seems to modulate the whole pattern. A different and less regular behaviour appears in the scale intervals of data which refer to smaller orbital periods (fig. 6b,c). In both cases, the highest maxima dominate over the noise level and the spectral energies seem to have the same modulation found in the previous band. The global wavelet spectrum of the simulated data has been calculated, following the same procedure used for paleomagnetic data. This spectrum is plotted in fig. 7 together with that of Mt. Raggeto data. The two series exhibit comparable energy content. Both have been fitted with piece-linear shaped functions, each of which presents two different slopes. A roughly similar slope is found in the high-frequency regions, whereas in the low-frequency region paleomagnetic data exhibit a less steep spectrum. Fig. 6a-c. Scale averaged local Morlet spectrum, versus time (Myr). Intervals averaged are: a) 120-80 kyr; b) 50-30 kyr, and c) 30-15 kyr. The per-interval 90% and 50% red noise confidence levels are indicated by dashed lines. a b c 581 Wavelet analysis on paleomagnetic (and computer simulated) VGP time series This could be expected as a consequence of the low intensity of magnetization of these data, that could result in a higher amount of random noise and as a consequence in a flattening of the spectrum. 5. 1DOG results To analyze excursion and reversals, preli- minary tests using, beyond the 1DOG, the 2DOG and 3DOG (second and third derivatives of a Gaussian, respectively) have been carried out both on simulated and paleomagnetic data sets. The tests were undertaken to check out the possible different features that different wavelets can reveal in the data. In fact, increasing the order of derivatives, the number of zeroes and vanishing moments of the analyzing wavelet increases. As a consequence, it is possible to resolve local features with greater Lipschitz coefficients and the multiplicity of the ridges increases (Mallat and Hwang, 1992; Alexandrescu et al., 1995), so that they can be compared among them. The different results we obtain exhibit only a progressive loss of resolution as the order of derivatives increases, both in terms of coefficients images and ridges in them. As the processes we are exploring do not need, in our working hypothesis, more than one Fig. 7. Comparison between global wavelet spectra of paleomagnetic and simulated data (see text for more de- tails). 582 Stefano Lorito, Grazia Giberti, Agata Siniscalchi and Marina Iorio vanishing moment, the value of the maximum positive Lipschitz exponent we expect being one, we choose to work with the 1DOG function. This choice was made also taking into account that a more complex structure of the wavelet, when associated with the inherent noise and high variability of the paleomagnetic data, would induce at any scale correlations which render difficult the localization and interpretation of singular behaviours. We begin reporting the analysis of the simu- lated data. As they have a simpler structure, some of the results features will guide the analysis of the more complex paleomagnetic data results. The obtained wavelet coefficients are represented in fig. 8 as an image in the time-scale plane. The excursions and reversals, whose localization and characteristics can be easily identified by visual inspection of the data (fig. 2), also show up clearly as extrema (ridges) of the CWT. Only some among the found ridges are selected and plotted in fig. 8, i.e. the most energetic (highest values of the CWT modulus) and those connected from finest to coarsest scales. Figure 9 shows the corresponding ridges functions versus the scale s in a log-log plot. The reversals (last two plots in fig. 9) are char- acterized by an almost linear trend up to the biggest scales, whereas the excursions ridge functions die out at relatively smaller scales. All the ridges, with the sole exception of the fourth, have at the finest scales the positive slope which characterizes the transition, in the time series, between the regular behavior and the sudden change we are looking for. The Lipschitz coefficients are listed in table I, which also reports the total duration of the processes as estimated from the scale (converted in time) at which the ridge functions die out. This time interval roughly corresponds to the scale at which the full process has been scanned by the wavelet, as far as the excursions regard and coincide for the first two ridges almost exactly with the time interval during which the first Fig. 8. Image of 1DOG wavelet coefficients of the simulated data. Red lines indicates the selected ridges (see text for details). Blue lines indicate the boundary cone of influence. The coefficients are scaled with respect to their maximum. 583 Wavelet analysis on paleomagnetic (and computer simulated) VGP time series excursion fully develops (cfr. fig. 2). As far as the second excursion is concerned, worse results in terms of ridge shape, exponents and possibility of estimating durations, are obtained and attri- buted to the greater (in comparison with the first excursion) variability of the data around it. As regards the reversals, the difference between their exponents a (0.2 and 0.5, respec- tively) reflects quite accurately the character expected when looking at the data. The first exponent corresponds to a (slightly) noisy Heavi- side step and the second to that of a noisy ramp. Both ridge functions conserve their almost linear trend within the boundary cone of influence, a limit imposed by the finiteness of the data set. It has to be noted that the estimated duration of 45 kyr of the first reversal (table I) corresponds to the time interval of its polarity stability, suddenly interrupted by a energetic excursion after the reversal. The wavelet coefficients obtained for pa- leomagnetic data are represented in fig. 10 as an image in the time-scale plane, and the corresponding selected ridge functions in fig. 11. The selected ridges are the most energetic and those connected from finest to coarsest scale, as above. The values of t i and a i and of the total duration of the identified processes are listed in table II. In the table we indicate which among all the selected ridges we interpret as reversals or excursions. This interpretation needs specific comments. The variability of the data series reflects in a less clear behaviour of the ridges functions. Linear trends are interrupted by oscillations, some randomly distributed, some probably due to a positive interference with the orbital periods which envelop the entire data series. This feature does not allow us to decide, only looking at the plotted ridges functions as done with the simulated data, the character of the polarity change. To interpret the ridge characteristics it has been necessary to zoom on segments and look if data trends around ridge roots would confirm the indications coming from ridge functions. In any case, couples of nearby ridges were expected to indicate ex- cursions, whereas a prevalent linear trend up to the coarsest scales along with higher absolute values to indicate reversals. With this in mind, only the ridges numbered 4, 10, 13, 14 and 21 could be identified as reversals. To these, with the exception of ridge 4, are also associated the longest estimated intervals of stability (about 100 L o g 2 ˜W /W m a xÁ Log2s Fig. 9. Simulated data ridge functions versus the scale in log-log plots. Table I. Results of the analysis of the ridge functions (fig. 9) of simulated data (see text for details). Excursions ti (kyr) a Dt (kyr) Ridge 1 27.360 0.1 � 10 Ridge 2 38.960 0.2 � 12 Ridge 3 116.640 0.7 � 15 Ridge 4 123.840 - 0.3 > 8 Reversals ti (kyr) a t (kyr) Ridge 5 190.800 0.2 > 45 Ridge 6 295.920 0.5 > 20 1 2 3 4 5 6 584 Stefano Lorito, Grazia Giberti, Agata Siniscalchi and Marina Iorio Fig. 10. Image of 1DOG wavelet coefficients of Mt. Raggeto latitude (see caption of fig. 8 and text for more details). Fig. 11. Mt. Raggeto latitude ridge functions versus the scale in log-log plots. L o g 2 ˜W /W m a x Á Log2s 585 Wavelet analysis on paleomagnetic (and computer simulated) VGP time series kyr). The identification of ridges numbered 3, 7, 8 and 9 as reversal has been made looking at the data, their ridge functions being less definite and resembling those of the remaining excursions. Actually, the stability in polarity around the corresponding transitions is found scarce but they do not look, for amplitude and duration, like a sequence of nearby excursions. As no nearby ridges were found to confirm them as excursions, we classify them as unstable reversals and their way to die out is attributed to the correlations at high scales with the neighboring data. The overall not well definite character of the ridge functions for the paleomagnetic data is reflected in the given ranges of the a exponents. Their values are on average around zero, with the exception for the couple of ridges numbers 5 and 6, both corresponding to a «d -like» excursion, with a duration of around 2 kyr. Finally, the negative exponents of the ridges 1, 2 and 12 indicate a sudden variation, in this being similar to ridges 5-6. The overall estimated duration of all excursions is always much less than 30 kyr, shorter than those estimated for the reversals; Table II. Results of the analysis of the ridge functions (fig. 11) of Mt. Raggeto latitude. Legend for the fifth column: rev. and exc. mean, respectively, reversal and excursion. « + Æ -» indicates a transition from normal to reversed polarity, « - Æ +» the opposite transition. n. t i (kyr) a Dt (kyr) 1 56.16 - 0.5 < a < 0 > 20 exc. 2 66.24 a < 0 > 20 exc. 3 306.72 – 0.1 < a < 0.1 > 30 rev. + Æ – 4 331.20 – 0.1 < a < 0.6 > 30 rev. – Æ + 5 422.64 – 1 < a < – 0.5 ~ 2 exc. 6 424.80 a � – 0.7 ~ 2 exc. 7 686.88 a £ 0 > 30 rev. + Æ – 8 715.68 – 0.5 < a < 0.3 ~ 30 rev. – Æ + 9 830.16 – 0.5 < a < 0.3 > 30 rev. + Æ – 10 966.24 a � 0 > 90 rev. – Æ + 11 1141.20 a < 0.2 ? exc. 12 1157.76 a � – 0.6 > 15 exc. 13 1299.60 – 0.5 < a < 0.5 ~ 100 rev. + Æ – 14 1530.00 0 < a < 0.5 ~ 100 rev. – Æ + 15 1675.44 a < 0.2 ~ 25 exc. 16 1686.96 a � 0.3 ~ 20 exc. 17 1972.80 a � 0.2 ~ 20 exc. 18 1985.04 a � 0.2 ~ 20 exc. 19 2067.12 – 0.3 < a < 0.4 ~ 20 exc. 20 2083.68 – 0.2 < a < 0.4 ~ 20 exc. 21 2126.88 – 0.5 < a < 0.2 > 100 rev. + Æ – 586 Stefano Lorito, Grazia Giberti, Agata Siniscalchi and Marina Iorio ridge 11 constitutes the only exception as it gives an anomalously long timescale. All the MR2-MR12 Mt. Raggeto polarity changes as detected in Iorio et al. (1998) result related to some of the identified reversals, with the exception of MR10 and MR11 that emerge as excursions in this analysis. An indication of the nature of the process can be retrieved also by comparison of the localization of ridge functions with those of the reversals in the coeval interval of the Global Polarity Time Scale (Channell et al., 1995), shown in fig. 12. If the Hauterivian- Barremian boundary (H/B) is well located in both the series, no correlation appears but, looking at the tentative correlations among chrone boundaries and taking into account the uncertainties in the H/B datation (± 1.6 Myr), a better correspondence would be achieved by shifting the data. 6. Conclusions The global Morlet CWT spectrum of paleo- magnetic data reveals some peaks related to orbital frequencies with a feeble statistical significance. The same periods are again evi- denced in the local CWT analysis of three selected ranges (respectively around 22, 40 and 100 kyr), with a better significance. The long eccentricity orbital period (� 400 kyr) seems to modulate the whole pattern. These results also indicate the non-stationary character of the analyzed time series. Records of geomagnetic paleointensity by sediments, spanning the time interval 0 -1 Myr, have been analyzed by a number of authors (Channell et al., 1998; Guyodo and Valet, 1999; Guyodo et al., 2000; Yokoyama and Yamazaki, 2000) to address the relevance of orbital con- tributions. Controversial results emerge. The debate is mainly about the possibility of climatic contamination that could affect the paleointensity of sedimentary records. More recently, Yamazaki and Oda (2002) reported evidence of eccentricity signature in the inclination of a sedimentary sequence (0-2.25 Myr). Here we have tested with the CWT technique previous results (Iorio et al., 1995) from the older (124-128 Myr) Mt. Raggeto sequence. These two last results could be indications of orbital influence on the geody- namo, as a mechanism of climatic contamination of paleomagnetic directions seems less obvious than for paleointensity. On the other hand, we noticed that the local energy maxima in the Morlet analysis (especially at 22 kyr) occur at roughly the same times as most of the found excursions and some reversals. This could indicate two opposite situations: the polarity transitions are influenced by orbitals or, vice versa, the energy maxima reflect the average durations of the transitions. Discussions on the same problem is found also in Worm (1997), Yokoyama and Yamazaki (2000) and Kent and Carlut (2001). Fig. 12. Comparison among the positions of found reversals by means of 1DOG CWT (gray and white) and the global polarity time scale (Channell et al., 1995). White corresponds to reversed polarity while black and gray to normal polarity. 587 Wavelet analysis on paleomagnetic (and computer simulated) VGP time series The global wavelet spectra of simulated and paleomagnetic data exhibit roughly similar slopes in the high-frequency region and this can be due to the presence in both series of the same maximum degree of singularity (Mallat, 1998). The observed differences could be explained when considering the age of Mt. Raggeto data and that simulated data are based on the present day CMB seismic tomography. By using 1DOG we localized reversals or excursions in the paleomagnetic data, with some difficulties due to the data inherent noise. The analysis of the simulated data was a useful guide in this connection. The analysis of the CWT results could be more robust if applied also to other paleomagnetic data sets with the possible widest geographical and temporal distribution and if also magnetized continental and oceanic lavas with the proper time resolution were available. Such a detailed study would serve to confirm the debated features of paleomagnetic data and render them as data on which theoretical models could be validated. This, in particular, as far as the possible role of orbital periods on the geodynamo is concerned (Malkus, 1968; Buffett, 2000; Roberts and Glatz- maier, 2000). All the results of this work have been ob- tained by analysis of latitude data only. The longitude and intensity of remanence were not used at this first stage as they are more difficult to treat with respect to the latitudes. In view of the results obtained so far and the methodological experience acquired, their analysis, along with the analyses of some other paleomagnetic se- quences, is now in progress. Acknowledgements We wish to acknowledge Sergio Cavaliere for suggestions about signal analysis and inter- pretation. As far as CWT is concerned, special thanks are due to Umberto Amato and the IAM-CNR people. We are indebted with Gary Glatzmaier and Lionel Hongre for kindly provi- ding access to their computer simulated data. The paper was greatly improved during its re- vision, thanks to the criticism of Maurizio Fedi and of an anonymous referee. NOTES The numerical codes used in the Morlet wavelet analysis have been developed adapting the Wavelet software provided by C. Torrence and G.C. Compo, which is available at URL: http://paos.colorado.edu/research/wavelets/. The numerical codes for the DOG wavelet analysis have been developed starting from those included in the freeware tool Wavelab 802, available at URL: http://www- stat.stanford.edu/~wavelab/. REFERENCES ALEXANDRESCU, M., D. GIBERT, G. HULOT, J.-L. LE MOUËL and G. SARACCO (1995): Detection of geomagnetic jerks using wavelet analysis, J. Geophys. Res., 100 (B7), 12,557-12,572. BERGER, A., M.F. LOUTRE and J. LASKAR (1992): Stability of the astronomical frequencies over the Earth’s history for paleoclimate studies, Science, 255, 560-566. BUFFETT, A.B. (2000): Dynamics of the Earth’s core, Geophys. Monogr., 177, 37-61. CHANNELL, J.E.T., E. ERBA, M. NAKANISHI and K. TAMAKI (1995): Late Jurassic - Early Cretaceous time scales and oceanic magnetic anomaly block models, in Geochronology, Time Scales and Global Stratigraphic Correlation, edited by W.A. BERGGREN, D.V. KENT, C. SWISHER III and M.P. AUBRY, Soc. Econ. Paleontol. Mineral, Spec. Publ., 54, 51-63. CHANNELL, J.E.T., D.A. HODELL, J. MCMANUS and B. LEHMAN (1998): Orbital modulation of the Earth’s magnetic field intensity, Nature, 394, 464-468. COE, R.S., L. HONGRE and G.A. GLATZMAIER (2000): An examination of simulated geomagnetic reversals from a palaeomagnetic perspective, Phil. Trans. R. Soc. London Ser. A, 358, 1141-1170. D’ARGENIO, B., V. FERRERI, M. IORIO, A. RASPINI and D.H. TARLING (1999): Diagenesis and remanence acquisition in the Cretaceous carbonate sediments of Monte Raggeto, Southern Italy, in Palaeomagnetism and Diagenesis in Sediments, edited by D.H. TARLING and P. TURNER, Geol. Soc. London, Spec. Publ., 151, 147-156. DAUBECHIES, I. (1992): Ten Lectures on Wavelets (SIAM, Philadelphia), pp. 376. DORMY, E., J.P. VALET and V. COURTILLOT (2000): Numerical models of the geodynamo and observational constraints, Geochem. Geophys. Geosyst., 1, 2000GC00062. GLATZMAIER, G.A., R.S. COE, L. HONGRE and P.H. ROBERTS (1999): The role of the Earth’s mantle in controlling the frequency of geomagnetic reversals, Nature, 401, 885-890. GUBBINS, D. (1999): The distinction between geomagnetic excursions and reversals, Geophys. J. Int., 137, F1-F3. GUYODO, Y. and J.P. VALET (1999): Global changes in intensity of the Earth’s magnetic field during the past 800 kyr, Nature, 399, 249-252. GUYODO, Y., P. GAILLOT and J.E.T. CHANNELL (2000): Wavelet analysis of relative geomagnetic paleointen- 588 Stefano Lorito, Grazia Giberti, Agata Siniscalchi and Marina Iorio sity at ODP Site 983, Earth Planet. Sci. Lett., 184, 109-123. IORIO, M. (1995): High resolution palaeomagnetic analysis of Cretaceous shallow water Carbonates, Ph.D. Thesis, University of Plymouth, U.K. IORIO, M., D.H. TARLING, B. D’ARGENIO, G. NARDI and A.E. HAILWOOD (1995): Milankovitch cyclicity of mag- netic directions in Cretaceous shallow-water carbonate rocks, Southern Italy, Boll. Geofis. Teor. Appl., 46, 109-118. IORIO, M., D.H. TARLING and B. D’ARGENIO (1998): The magnetic polarity stratigraphy of a Hauterivian/ Barremian carbonate sequence from Southern Italy, Geophys. J. Int., 134, 13-24. IORIO, M., D.H. TARLING and B. D’ARGENIO (2001): Early Cretaceous geomagnetic behaviour compared to computer model, J. Geodyn., 32, 445-465. KENT, D.V. and J. CARLUT (2001): A negative test of orbital control of geomagnetic reversals and excursions, Geophys. Res. Lett., 28, 3561-3564. MALKUS,W.V.R. (1968): Precession of the Earth as the cause of geomagnetism, Science, 160, 259-264. MALLAT, S. (1998): A Wavelet Tour of Signal Processing (Academic Press, San Diego, CA), pp. 577. MALLAT, S. AND W.L. HWANG (1992): Singularity Detection and Processing withWavelets, IEEE Trans. Inf. Theory, 38, 617-643. MC MILLAN, D.G., C.G. CONSTABLE, R.L. PARKER and G.A. GLATZMAIER (2001): A statistical analysis of magnetic fields from some geodynamo simulations, Geochem. Geophys. Geosyst., 2, 2000GC000130. MERRILL, R.T. and P.L. MCFADDEN (1999): Geomagnetic polarity transitions, Rev. Geophys., 37, 201-226. MERRILL, R.T., M.W. MCELHINNY and P.L. MCFADDEN (1996): The Magnetic Field of the Earth: Paleo- magnetism, the Core, and the Deep Mantle (Academic Press, San Diego, CA), pp. 531. MORLET, J., G. ARENS, E. FOURGEAU and D. GIRARD (1982): Wave propagation and and sampling theory, 2. Sampling theory and complex waves, Geophysics, 47, 203-221. PERCIVAL, D.P. (1995): On estimation of the wavelet variance, Biometrika, 82, 619-631. PERRIER, V., T. PHILIPOVITCH and C. BASDEVANT (1995): Wavelet spectra compared to Fourier spectra, J. Math. Phys., 36, 1506-1519. ROBERTS, P.H. and G.A. GLATZMAIER (2000): Geodynamo theory and simulations, Rev. Mod. Phys., 72, 1081-1123. SCARGLE, J.D. (1982): Studies in astronomical time series analysis. II. Statistical aspects of spectral analysis of unevenly spaced data, Astrophys. J., 263, 835-853. TARLING, D.H., M. IORIO and B. D’ARGENIO (1999): Geomagnetic long-term secular variations in Italian Lower Cretaceous shallow-water carbonates, Geophys. J. Int., 137, 713-722. TORRENCE, C. and G.C. COMPO (1998): A practical guide to wavelet analysis, Bull. Am. Meteorol. Soc., 79, 61-78. WORM, H.U. (1997): A link between geomagnetic reversals and events and glaciations, Earth Planet. Sci. Lett., 147, 55-67. YAMAZAKI, T. and H. ODA (2002): Orbital influence on Earth’s magnetic field: 100 000-year periodicity in inclination, Science, 295, 2435-2438. YOKOYAMA, Y. and T. YAMAZAKI (2000): Geomagnetic paleointensity variation with a 100 kyr quasi-period, Earth Planet. Sci. Lett., 181, 7-14.