## Abstract

The virtual geomagnetic pole (VGP) trajectories during some geomagnetic polarity reversals of different ages are marked by anisotropic behaviour. This recurrent phenomenon may be reflected in the paleomagnetic data, even if the transitional field was not completely recorded. As the long-scale geomagnetic variations have a confined oscillatory character, the VGP paths from stratigraphically controlled sequences may be described on the basis of sine and cosine functions, even if time is not the independent variable. Here we considered longitude (or space) as the independent variable which had to be ‘unrolled’ to overcome the 360° repetitions as the VGPs moved around the geographic pole.

Sixteen VGP series from the Early Cretaceous Serra Geral lava flows of southern Brazil were analysed using a modified version of the periodogram for uneven data series, and a combination of information approach. The combination of all the spectra, as in a stacking procedure, reduces noise and results in a smooth curve highlighting features of interest. We found a set of highest correlation wavelengths of approximately 167, 190, 209, 257, 277 and 368°. Phase analyses using two different methods revealed strikingly good coherence for some of these wavelengths, indicating that they are not only artefacts of the spectral analysis. Similar analysis of magnetostratigraphic data from the Icelandic Magmatic Province indicated that the two datasets may have wavelengths of approximately 165 and 270° in common. These results suggest quasi-periodic behaviour, possibly with sub-harmonic instabilities owing to the modulating effect of inner Earth's anisotropies influencing the pole trajectory.

Understanding of the long-term variations of the Earth's magnetic field requires a full description of the spatial–temporal characteristics of the field. Analyses of secular variation have been mainly based on directional dispersion as a function of latitude in recent times (Tauxe & Kent 2004 and references therein), and very few analyses refer to older fields (e.g. Smirnov & Tarduno 2004; Biggin *et al.* 2008). Long-term geomagnetic variations are better observed in sedimentary rocks, as they may continuously record the behaviour of the paleomagnetic field, although knowledge of the time of the remanence acquisition is usually imprecise. Furthermore, depending on depositional rate, the geomagnetic record may be naturally smoothed. Lava flow sequences tend to produce a more faithful magnetic record. However, the intermittent character of volcanism does not allow continuous records, and the time gap between two consecutive flows is difficult to ascertain precisely. The character of the magmatic activity is such that adjacent volcanic sections may exhibit distinct paleomagnetic records of the same event (Valet & Herrero-Bervera 2003). This behaviour may prevent the identification of some features of the geomagnetic field. However, multiple magnetostratigraphic records may display complementary information about certain geomagnetic characteristics that can be assessed if the information is adequately combined.

Magnetostratigraphic records allow the observation of virtual geomagnetic pole (VGP) trajectories throughout the time interval spanned by a given geological formation. One interesting issue resulting from the comparison of VGP trajectories during recent geomagnetic reversals is that transitional VGPs may follow preferred paths along some longitudinal bands (Clement 1991; Laj *et al.* 1991; Hoffman 1992; Gubbins & Love 1998; Love 2000). This type of observation made on sedimentary sequences has been viewed with scepticism by some authors (Barton & McFadden 1995; Quidelleur *et al.* 1995), who claim that the phenomenon is an artefact of the magnetization process in those rocks. However, some volcanic rocks also show geographical VGP confinement, as pointed out by Love (1998), Valet & Herrero-Bervera (2003), Coe & Glen (2004), and others.

Most of the observations of transitional VGP confinement are for reversals from the last 5 Ma, but older records (Late Permian through Early Cretaceous; Vizán *et al.* 1994) may show the same characteristics, at least on a statistical basis. Simulated reversals using the Glatzmaier–Roberts geodynamo model may also support (Coe *et al.* 2000) the preferred band hypothesis. Reversals may differ from one another, and the same reversal observed in different places may show distinctive records owing to the geomagnetic non-dipolar components that prevail at the sampling sites. However, the repetition of some observed patterns in long-term geomagnetic variations suggests the action of some regulatory phenomena (Clement 1991; Laj *et al.* 1991; Hoffman 1992; Gubbins & Love 1998; Coe *et al.* 2000; Costin & Buffet 2004; Leonhardt & Fabian 2007; Vizán & Van Zele 2008). Models for VGP paths (Kuznetsov 1999; Leonhardt & Fabian 2007) suggest that the trajectories may be conditioned by field anomalies superposing the dipolar field during reversals, leading to longitude confinements.

Unfortunately, many magnetostratigraphic datasets from igneous provinces are unsuitable for the investigation of transitional VGP confinement owing to their fragmentary character. They may show various polarity intervals but few or practically no transitional directions. However, they may carry valuable information on the behaviour of the magnetic field that may help better constrain geodynamo modelling. Volcanic magnetostratigraphic sequences are still underexplored, but advances may occur if new methods of analysis are developed and the problem is examined from a new perspective. This is the main objective of the present study.

In this study we apply a spectral analysis to the magnetostratigraphic sequences of the Lower Cretaceous Serra Geral Formation (Paraná Magmatic Province, southern Brazil). VGP sequences from several Serra Geral sections (Fig. 1) have already been obtained (Ernesto *et al.* 1990; Alva-Valdivia *et al.* 2003; Mena *et al.* 2006), but no complete reversal has been described despite the various polarity intervals recorded in the series (Fig. 2). However, if the VGP trajectories during polarity reversals or excursions are marked by anisotropic behaviour, it is possible that the datasets contain some recurrent features. In this way, if we approximate the VGP trajectory as a sum of oscillatory modes around a line (equator) in a 2D (latitude v. longitude) plane, then spectral properties (or correlations with sines and cosines) can be investigated. Because geochronologic control of the geomagnetic events is unavailable, we must proceed in a way that does not require an assumption about elapsed time. Here we address the problem by considering the spatial rather than the time behaviour of the VGP paths, taking longitudes as the independent variable.

## The dataset

We analyse the Serra Geral Formation data published by Ernesto *et al.* (1990). These data come from 16 stratigraphically ordered lava piles from the southern part of the Paraná Magmatic Province (Fig. 1) where the best rock exposures are concentrated. In a recent re-evaluation of the available ^{40}Ar/^{39}Ar data, Thiede & Vasconcelos (2010) proposed an age interval of 134±1 to 133±1 Ma for the magmatic activity in that area, with younger ages towards the north. This short time interval of about 1 Ma for the emplacement of the magma flows, and the small differences in latitude (less than 4°) between lava flow sequences, are good indicators that the rocks were magnetized by a relatively homogeneous paleomagnetic field. Most of the sequences record more than one polarity interval, although transitional directions are rare, as can be seen on the declination/inclination plots of Figure 2.

For the purpose of present analysis, the VGPs coordinates must be corrected of the effects of the South American plate displacement that has occurred since the Early Cretaceous. Hence, a correction of 5° to the whole dataset was applied taking into account the reference paleomagnetic pole given by Font *et al.* (2009) for the Early Cretaceous. The VGPs describe irregular loops around the geographic poles (north or south depending on the polarity state), causing longitudes to repeat. Therefore, the VGP trajectories must be ‘unrolled’, eliminating the longitude repetitions by adding 360° each time that data values moved backwards, as seen on the examples of Figure 3. This procedure does not affect the distribution of VGPs, and does not modify the spectral content in the series. This unwrapped sequence of longitude–latitude pairs can now be treated as an unevenly sampled ‘time’ series. Figure 4 shows the plots of the 16 unwrapped VGP sequences used in this work.

## Methods

We seek to detect the long-term features in the VGP trajectories. Therefore, we have to model a smooth curve that fits the main path described during the secular variation cycles or complete polarity changes, regardless of the short-wavelength features that are responsible for the stationary features causing the VGPs to cluster at some latitudes. Once the VGPs in each series are arranged in order of increasing longitudes (Fig. 4), the various data sequences can be treated as multiple independent samples of the same trajectory of a particle over an infinite 2D plane (latitude v. longitude). We want to recover the original function that gave rise to those samples. In a simplified picture, this function describes the trajectory of a particle with a general trend to the east/west, oscillating around the 90° parallels most of the time, and eventually crossing the equator. We have only randomly spaced sampled series for defining this function, but we can draw attention to three main constraints of the original function: (1) it is smooth (infinitely differentiable) – we are only interested in long period behaviour; (2) it is limited – all data must be contained in a limited space delimited by a finite number of parallels; and (3) it oscillates most of the time. The third condition suggests using Fourier basis to model the data, and the first two conditions allow us to do so mathematically.

The problem of estimating the original series from which we have only the unevenly sampled series may be solved by invoking the Popper–Bayes paradigm (Tarantola 2006): the data should not be used to choose a ‘best’ solution but instead to ‘falsify’ solutions over the parametric space. By ‘falsifying’ the author means setting a probability distribution over a subset of the parameter and improving that distribution by adding (combining) new experimental (or theoretical) information. In this framework, a Bayesian inversion scheme can be thought of as a combination of theoretical and experimental information, or as a combination of independent experimental information (Tarantola & Valette 1982). Therefore, the inversion of the spectral state of information functions we are dealing with can be the combination of the multiple records of the same events. These state of information functions are freely normalizable probability distributions over a finite discrete set of points in the frequency domain. For each frequency we associate a probability value given by its spectral power estimate – the periodogram ordinate. For this, we made use of the optimal least-squares fitting properties of the Lomb–Scargle (LS) periodogram (Lomb 1976; Scargle 1982) for decomposing the data series variance using sinusoidal bases.

The periodogram is a suitable tool for spectral analysis of unevenly sampled time series, although it produces very noisy estimates even when the sample data are only slightly noisy (Scargle 1982; Hernandez 1999; Zechmeister & KürSter 2009). A solution could be to smooth the periodogram estimates by enlarging the main lobe or applying filters (windows), therefore reducing resolution, and consequently reducing the main lobe power (and detection efficiency). As a result, the power leakage from an existing spectral component to neighbour frequencies is also reduced. This is the basis of the so-called ‘smoothed periodogram method’.

As in the case of evenly spaced data, we need to define the frequency set for which power will be estimated. Usually the maximum frequency corresponds to the ‘Nyquist frequency’, which is the inverse of the half the minimum distance between any two consecutive points, the minimum frequency being the inverse of the total time interval covered by the data series. The frequency increment is taken as one and a half of the inverse of the whole series length. When the investigated phenomenon is sampled in an almost regular grid containing only small gaps or a few missing points, a characteristic sampling interval may be estimated by the average interval or by more elaborate estimates of the number of independent points in frequency domain (Horne & Baliunas 1986).

In the present case the longitude difference between two consecutive points is extremely variable, and a new procedure is required for the estimate of the characteristic sampling interval as well as the Nyquist frequency. First, we set the frequency limits of the estimated spectra numerically in order to combine all information available in the series. The maximum available frequency range was computed as follows: from the inverse of the largest series length to the inverse of the shortest spacing between two consecutive points in any of the investigated data series (a kind of Nyquist frequency). Next, we reduced those limits in order to eliminate regions with highly incoherent signal from the estimated spectra (by comparison of all series), and regions with no resolution in the parameters estimation. With respect to spacing, we started from the smallest frequency spacing (the inverse of the greatest length), and progressively enlarged the spacing until we achieved smoother estimates for all individual series.

The challenge is to smooth those estimates, and this can be achieved through the combination of the information contained in the various data sequences. The obtained spectra can be compared and combined after being properly normalized. The adopted normalization parameter is the apparent total power, which corresponds to the power over the total investigated bandwidth. Each spectrum is normalized by the total area in such a way that the area under the curve is equal to unity; the most important aspect in this procedure is that the signal/noise ratio is preserved.

In the Bayesian approach we model our knowledge about a parameter rather than the parameter itself. Therefore, the determination of a restricted space solution does not necessarily mean that we know the parameter value with absolute precision (Tarantola 2006). We might interpret our obtained values (probability over frequencies) based on the data sample uncertainties and on the quantity (and quality) of the series. In this case, each sampled series represents an experiment of the spectral measurement or of the ‘roughness’ of the original series. In other words, if the distribution of the combined variable (estimation) is concentrated, it could be just because the combining distributions have small intersections or few common models, despite its low information content. In this approach we do not extract any central estimator (mean or variance), which according to Tarantola & Valette (1982) ‘would certainly lead to a loss of information’. These authors indicate the use of the discrete probability distribution function (state of information function) itself as an elementary datum.

The combination of the sampling functions furnishes two envelopes in the parametric space – one is larger and inclusive, and the other is smaller and exclusive. These envelopes respectively represent the union and intersection of the various states of information (as in usual set theory), and they correspond to the ‘OR’ and ‘AND’ operations (Tarantola & Mosegaard 2000). Mathematically, these operations can be performed as the sum (OR) and the product (AND) of the different estimates for each frequency. For the ‘AND’ operation we divide the resulting product by the null state of information function as defined by Tarantola & Valette (1982). We chose this function as a constant function over the interval, which is satisfactory in most cases as stated by those authors. The resulting combined estimate (combination of the multiple spectra) is also normalized to the total unit area.

## Results

### Periodogram analysis

In order to calculate the individual LS periodograms for the Serra Geral VGP sequences we had to define the sampling density as discussed above. The appropriate sampling density was then chosen empirically observing the following aspects: (1) as the number of points in the frequency domain increases more noise is introduced; (2) fewer points in the frequency domain imply more ‘scalloping error’, that is, the error that arises when the real signal is between two observing points in the frequency domain; and (3) the chosen spacing interval should allow observation of the common spectral features in most of the spectra.

The LS periodograms for the VGP sequences are given in Figure 5. The upper part of the figure (a–c) shows the spectrum for each of the investigated sequences. As expected, the spectra are extremely noisy, especially in the shorter wavelength range, and few sections give clear indications of real wavelengths. However, it is easily seen on the individual curves that some common features appear in all spectra, although sometimes as a weak signal (below 20% of the maximum peak height).

### Normalization and information states

The mean of all information in the individual spectra after normalization corresponds to the ‘OR’ curve (Fig. 5). Although still noisy, this is a smoother spectrum, mainly for long wavelengths, where probable real ‘models’ can be envisaged. The ‘AND’ curve (bottom of Fig. 5) represents a non-linear cascade of amplifiers, one for each data series furnishing one gain factor with encapsulated signal/noise ratio information. This curve combines only the common features in all individual spectra, and shows a surprisingly small number of discrete peaks, considering the large noise level and very low density of points relative to the length of the individual series. As a result we obtained a measurement of the probability over a subset of the frequency space that indicates which models fit better. In short, we used the information content in each series to falsify possible solutions in a subset of the parametric space. Therefore, we combined the general available information for all series in relation to the probability that a particular frequency can be a true Fourier component of the original trajectory.

Any spectral analysis is ultimately a correlation analysis of a finite, discretely sampled series with an infinite function basis (as in the LS periodogram), and so is subject to power leakage and aliasing. Therefore, we are interested not only in the maximum ordinate values of the resulting ‘spectra’, but we also have to investigate the local behaviour of the correlation maxima (neighbour frequencies), for which the continuity and the smoothness of the final distribution is important.

We found a set of highest correlation wavelengths of approximately 167, 190, 209, 257, 277 and 368°. The frequencies on the peaks λ_{1} = *c.* 167°, λ_{2} = *c.* 257° are respectively 4 and 5 times greater than that at λ_{3} = *c.* 368°. Then we infer a sub-harmonic relationship between λ_{1} and λ_{2}, namely λ_{1}/λ_{2} = *c.* 2/3, that is, λ_{1} and λ_{2} seem to be the second and third harmonic frequencies of a process for which λ_{f} = *c.* 514° would be the fundamental one.

## Performance of the method with simulated data

Synthetic data series possessing known harmonic composition, and with similar sampling characteristics to the real data studied in this work, were used to test the method. We considered a time series resulting from the superposition of sinusoidal components with known periods of 350, 480, 1100 and 4000 and amplitudes equal to 4, 6, 7 and 6, respectively. Noise proportional (5%) to the amplitude was added to each point. The resulting time series was then randomly sampled and uniformly distributed producing six different series as follows: three samplings of the whole series (21 cycles) with 25, 20 and 18 points (datasets 1–3); one sampling from the first two-thirds of the series with 15 points; one sampling from the last two-thirds of the series with 15 points; and one sampling with eight points from the middle third of the series. Figure 6 shows the resulting time series and sampled datasets.

The combined spectra for the six synthetic series are shown in Figure 7. The upper window (a) shows the spectrum for each series separately; the other windows show the result of the union and intersections of the state of information functions (operations ‘OR’ and ‘AND’). The spectral combinations are able to filter most of the noise or spurious signals, and they highlight the signal in all the series (‘AND’ operation) for which the probability is higher. Despite the poor quality of the generated data series, the procedure presented here shows excellent results considering the limited possible solutions in the parametric space. There are high-correlation regions in the OR and AND spectra, and the imposed periodicities (480, 1100 and 4000) in the original time series were satisfactorily recovered (487, 863 and 4207). In the investigated samples, no solution was found for frequency 350 based on the states of information used here.

The above results demonstrate that the method performs quite well for unevenly spaced time series affected by noise and with low density of points.

## Discussion

The wavelength correlation suggested in the Serra Geral periodogram of Figure 5 needs confirmation; it is not possible to distinguish between actual components and spurious spectral effects based only on this analysis. Alternatively, we may suspect a ‘coupling’ effect of the spatial propagation modes or a spatially periodic variability in the media. These results seem to reinforce the initial guess that there could be periodic or quasi-periodic components in the VGP trajectories, and although these trajectories are extremely complex, their long-term behaviour can be described by a relatively small number of sinusoidal components.

### Spectra stability

In the case of data series sampled on a regular grid, the application of spectral windows to reduce lateral lobes owing to spectral leakage is a common procedure. This practice is not recommended in the irregular case (Hernandez 1999) because the spectral windows will be asymmetrical and non-homogeneous owing to irregularities of the original sampling. Furthermore, adapting the windows in the case of data distribution with strong asymmetries in spacing, as is the present case, is not trivial. Therefore, we applied two variations of the rectangular window: (1) we removed the extreme points of each series; and (2) we gradually attenuated the three first and last points of the series to levels of 20, 30 and 70% of the original values from the extremity inward. In both cases, the objective was to modify the structure of the lateral lobes in each spectrum, altering their combination. We intended to attenuate side lobes relative to the main lobes in each case; in doing so, we tried to separate those features generated exclusively by spectral leakage from the existing components in the original series. The results are given in Figure 8. The two procedures preserved some characteristics of the former spectrum, but considering the shortness of the original series, the first adopted procedure apparently implied a considerable loss of information. Even so, peaks *c.* 255–270° still persisted, and peaks *c.* 450°–496°, which are close to the λ_{f} = *c.* 514° discussed above, are now evident.

The stability of the spectrum was also tested by disturbing the points of the original data series with uniformly distributed noise of the form (*x _{i} = x_{i} + K*rand*), with

*K*= 5, 10, 15, 20 and 30°. This kind of filter is widely used in image processing and acts as a low-pass or ‘anti-aliasing’ filter. In Figure 9 the stability of the spectrum is noticeable as the main features were maintained. This indicates that the spectral characteristics can really represent long-term components of the trajectory, although some power transfer among peaks still occurred for increasing values of

*K*, probably owing to aliasing.

### Phase determination

Phases were determined by two methods: the classical Fourier phase estimator and the Hocke's (1998) estimator based on the LS periodogram. Despite the poorness of the sampling series, there is good agreement among the estimated phases for each data series, as seen in the rose diagrams of Figure 10. Apparently, the results for the wavelengths of 167, 277 and 368° cluster better, and this can be checked in Figure 11, where the results for each data sequence can be visualized. The fact that one or more sinusoidal components are in phase cannot be seen as a chance result, indicating that at least one of the components might be real. Furthermore, in the tests of spectrum stability, a large amount of power was transferred to the *c.* 167 and *c.* 277° peaks as the original data were becoming noisier (Fig. 9), indicating that probably they are not spurious effects.

### Comparison with Icelandic data

At this point, it is instructive to investigate whether the characteristic wavelengths found in the Serra Geral spectra also represent a geometric feature in the VGP paths of other datasets. The Icelandic Magmatic Province is suitable for a comparative analysis, as the available paleomagnetic data (Kristjansson *et al.* 2004) show similar characteristics in the record of secular variation and reversals. The data come from sections exposed in the northern part of Iceland, and the ages are within the range of 9–5 Ma. The individual sections (coded GR, KG, AF, AS, GL, SG and VE; Kristjansson *et al.* 2004) are partly superposed, and the composite section displays 17 reversals and excursions, but each one appears in one or two flows. To perform the spectral analysis in the same bandwidth as the Serra Geral spectra, the longer profiles GR and KG were subdivided, and only the upper parts were analysed. The VGP longitudes were ‘continued’ (Fig. 12) according to the same criteria used in the Serra Geral data.

The ‘AND’ spectrum of the Iceland data (Fig. 13) is smooth, with pronounced peaks at 169.1 and 265.2°. The latter is predominant, which means that it is the most coherent feature in all the Iceland series. In the Serra Geral spectra, the highest peak is at *c.* 168°, although the doublet 257.6–277.3 also corresponds to a considerable part of the spectrum's energy (Fig. 5). This similarity between Iceland and the Serra Geral is somewhat surprising, as the two datasets refer to very different ages and geographic locations, and continuity in the geomagnetic field's behaviour is not necessarily expected.

## Spatial behaviour of the VGP paths

The Serra Geral magnetostratigraphic sequences are appropriate for this type of analysis despite the low density sampling of the long-term geomagnetic variations; the high number of available sections allows a considerable smoothing of the individual curves. Considering the performance with synthetic data and the quality of the real data spectra, it is reasonable to say that we have successfully detected some characteristic wavelengths in the VGP trajectories. Moreover, independent datasets like those from the Icelandic magmatism have indications of similar geometric characteristics in the VGP paths. It is important to stress that both the Serra Geral and Icelandic data record several geomagnetic polarity reversals, in both the N→R and R→N senses, and the record of the transitions themselves are scarce. Therefore, as a conclusion we may say that, although not fully observed, the secular variation cycles and/or polarity reversals show common features that are persistent for the time interval covered by the datasets, and that the general geomagnetic field behaviour may be similar during the Lower Cretaceous and younger ages.

The fact that the VGP trajectories may be modelled by the combination of some sinusoidal components does not mean that periodic phenomena were identified. We may only say that we have detected some possible wavelengths that are recurrent in all the investigated data sequences, and they clearly point to the existence of anisotropy in the VGP paths during secular variation cycles or polarity changes. Models for VGP paths during recent reversals were proposed by Leonhardt & Fabian (2007) and Kuznetsov (1999), among others, for which the trajectories may be conditioned by field anomalies superposing the dipolar field during polarity reversals, and leading to longitude confinements.

In their analysis of Jurassic to Early Cretaceous intermediate VGPs, Vizán & Van Zele (2008) found data concentrations mainly at longitude bands 300–330°E and 120–150°E, indicating a strong asymmetry of the deep-earth anomalies causing the deviation of the VGP trajectories. The angular differences between these two bands are in the interval 150–210°, and may be easily associated with the *c.* 168–270° wavelengths found in this work. One possible interpretation is that this is the range of longitudes over which VGPs undergo changes in latitude as they describe loops of secular variation or move from one polarity state to the other. Therefore, these results might represent some of the components of the oscillatory VGP trajectories. The recovery of the complete trajectory should be possible in principle by the summation of all components if the respective amplitudes and phase distribution are known.

## Conclusions

We have attempted to identify sinusoidal components in VGP trajectories recorded in magnetostratigraphic data of volcanic sequences using a combination of information approach. The starting point was the LS periodogram for unevenly sampled data series, although we do not consider the usual power spectrum, that is, in this case the independent variable is space and not time.

Without any characteristic sampling interval, and because data are too sparse in face of the phenomena of interest, the individual spectra did not give reliable information (periodicity or wavelength), or they exhibited a very weak signal. In the latter case, the problem could be by-passed, making use of some properties of the spectral estimator to suppress the spectral leakage. Power spreading across the (wavelength) frequency domain helped in the signal identification as the addition of more serial information (other data series) enhanced the common signal.

As in seismic signal stacking, the combination of multiple spectra reduces noise and highlights the common features even if the signal is weak. The solution in the parametric space is delimited in terms of possible existing highly correlated trigonometric components. As a result, we obtained a smooth curve highlighting the features of interest.

The great advantage is that this method does not require any data pre-processing (e.g. de-trending, filtering or pre-whitening) that would modify the power spectrum and possibly cause the displacement of spectral peaks or modification in power (Hernandez 1999).

We found a set of highest correlation wavelengths of approximately 167, 190, 209, 257, 277 and 368°. Phase analyses using two different methods revealed strikingly good coherence for some of the wavelengths, indicating that they might not be artefacts of the spectral analysis. Similar analysis of magnetostratigraphic data from the Icelandic Magmatic Province indicates that the two datasets have two wavelengths of approximately 165 and 270° in common. These results suggest quasi-periodic behaviour, possibly with sub-harmonic instabilities owing to the modulating effect of inner Earth's anisotropies influencing the pole trajectory.

## Acknowledgments

This work was funded by the Brazilian funding agency FAPESP (2004/05363-5). Thanks are due to W. Shukowsky and L.A. Diogo for valuable comments on the methodology. The paper was greatly improved by constructive comments and suggestions from X. Zhao, another unknown reviewer, and the special editor for this volume, L. Hinnov.

- © The Geological Society of London 2013