# Estimation of the gravimetric pole tide by stacking long time-series of GGP superconducting gravimeters

Yann Ziegler, Jacques Hinderer, Yves Rogister, Séverine Rosat

► To cite this version:

Yann Ziegler, Jacques Hinderer, Yves Rogister, Séverine Rosat. Estimation of the gravimetric pole tide by stacking long time-series of GGP superconducting gravimeters. Geophysical Journal International, 2016, 205, pp.77-88. 10.1093/gji/ggw007. hal-01282694

HAL Id: hal-01282694

https://hal.science/hal-01282694

Submitted on 5 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Estimation of the gravimetric pole tide by stacking long time-series of GGP superconducting gravimeters

Yann Ziegler, Jacques Hinderer, Yves Rogister and Séverine Rosat

IPGS/EOST, Université de Strasbourg/CNRS 5 rue René Descartes, F-67084 Strasbourg Cedex, France. E-mail: yann.ziegler@unistra.fr

Accepted 2016 January 6. Received 2016 January 5; in original form 2015 July 15

## SUMMARY

We compute the gravimetric factor at the Chandler wobble (CW) frequency using time-series from superconducting gravimeters (SG) longer than a decade. We first individually process the polar motion and data at each individual gravity station to estimate the gravimetric factor amplitude and phase, then we make a global analysis by applying a stacking method to different subsets of up to seven SG stations. The stacking is an efficient way of getting rid of local effects and improving the signal-to-noise ratio of the combined data sets. Using the stacking method, we find a gravimetric factor amplitude and phase of $1.118 \pm 0.016$ and $-0.45 \pm 0.66$ deg, respectively, which is smaller in amplitude than expected. The sources of error are then carefully considered. For both local and global analyses, the uncertainties on our results are reliably constrained by computing the standard deviation of the estimates of the gravimetric factor amplitude and phase for increasing length of the time-series. Constraints on the CW anelastic dissipation can be set since any departure of the gravimetric factor from its elastic value may provide some insights into the dissipative processes that occur at the CW period. In particular, assuming given rheological models for the Earth's mantle enables us to make the link between the gravimetric factor phase and the CW quality factor.

**Key words:** Time variable gravity; Earth rotation variations; Elasticity and anelasticity.

## 1 INTRODUCTION

The movement of the Earth rotation axis is made of various components. The Chandler wobble (CW) is a rotational eigenmode with a period around 435 sidereal days, whereas all the other components are motions forced by external sources with periods ranging from a few hours to several years and more. Contrary to the frequencies of the forced polar motions, which are mainly related to astronomical parameters, such as the positions and masses of the Sun and the Moon, the CW eigenfrequency only depends on Earth parameters. Another obvious distinction between the CW and the forced motions comes from the excitation sources. Although still debated, the CW excitation is thought to originate mainly from the oceans and atmosphere (e.g. Wahr 1982, 1983; Gross 2000; Brzeziński & Nastula 2002; Aoyama et al. 2003) and, maybe, earthquakes (Smylie et al. 2015). Geomagnetic jerks have also been proposed to explain phase jumps (Bellanger et al. 2001; Gibert & Le Mouël 2008).

Another important feature of the polar motion is its attenuation with time due to dissipative processes. Dissipation occurs in different ways: anelastic deformation in the mantle, viscomagnetic coupling at the core–mantle boundary, friction at the bottom of the oceans, etc. The quality factor *Q* is a parameter that quantifies the dissipation. If dissipation occurring in the mantle is the main source of attenuation, *Q* may be related to the rheological parameters of

the mantle. Its determination at the CW frequency can then provide additional constraints and new insight on the Earth rheology.

Besides, the movement of the Earth rotation axis induces a perturbation of the surface gravity field through (1) the variation of the centrifugal pseudo-force and (2) surface deformation and mass redistribution. These changes have been successfully observed thanks to superconducting gravimeter (SG) measurements. Because calibration changes in SGs are negligible and because their drift rates are quite small (Hinderer et al. 2007), it has been shown that SGs are well suited for studies of long-term polar motion (Loyer et al. 1999; Xu et al. 2004; Ducarme et al. 2006; Hu et al. 2007; Chen et al. 2009). Over the last 15 yr, departures from purely elastic response to polar motion have been sought in gravimetric signals. Essentially, within the given uncertainty range, all of the results obtained so far failed in revealing any anelastic effects, although the possibility of their existence has not been excluded.

In this study, we similarly address this question of anelasticity by investigating an alternative method that makes profit of the geographic coverage of the SG records. As shown in Section 3, the main difficulties in SG time-series processing are due to local gravimetric perturbations, in a large sense that encompasses both instrumental artefacts and local geophysical effects, the key point being that the needed corrections are both thorny to handle and site-specific. We endeavour to circumvent this complexity by stacking the gravity
time-series, as explained in Section 4, to decrease the relative impact of local effects unrelated to the polar motion. Doing so, we aim at increasing both the accuracy and the precision needed to study anelasticity at the CW period.

Importantly, this analysis relies on time-series longer than a decade, as is the common time span used in the stacking. This long duration of the recordings is one pillar of this kind of work knowing that the Chandler movement and the annual signal are only theoretically separable in the Fourier domain for time-series longer than 6.5 yr approximately, when the frequency resolution (equal to $1/T$ where $T$ is the length of the time-series) is, at least, equal to the separation in frequency between the Chandlerian and annual peaks. Thus, only the longest gravity records will let us distinguish the two components and provide reliable results.

## 2 OBSERVATION OF GRAVITY VARIATIONS

The ground measurement of the time variations of the gravity provides a multicomponent signal from which we want to extract the effect of periodic polar motion. We write the total measured gravity variation $\Delta g$, either in the time domain or in the frequency domain, as the sum of the polar motion perturbation $\Delta g_m$, a geophysical signal $\Delta g_{geo}$, a term of instrumental origin $\Delta g_{SG}$ and random noise n:

$$ \Delta g = \Delta g_m + \Delta g_{geo} + \Delta g_{SG} + n. \quad (1) $$

$\Delta g_{geo}$ contains the gravity variations due to the solid Earth tides and to the atmospheric, hydrological and oceanic perturbations. $\Delta g_{SG}$ combines the drift, occasional instrumental defects and secular or periodic artefacts.

The correction for or removal of $\Delta g_{geo}$, $\Delta g_{SG}$ and n will be considered in Section 3.2. It is impossible, by definition, to remove random noise from a signal but we can decrease its relative amplitude by combining different time-series, as explained in Section 4.

Once we will be able to determine $\Delta g_m$ from the observations, say $\Delta g_m^{obs}$, we will compare it to the theoretical gravimetric perturbation $\Delta g_m^{rig}$ of a rigid Earth model undergoing the same polar motion (Wahr 1985):

$$ \Delta g_m^{rig} = \Omega_0^2 r [\sin 2\theta(m_1 \cos \lambda + m_2 \sin \lambda) - 2m_3 \sin^2 \theta], \quad (2) $$

where $r$, $\lambda$ and $\theta$ are the spherical coordinates of the station, $\Omega_0$ is the unperturbed angular velocity of the system of reference attached to the Earth and $m_1$, $m_2$ and $m_3$ are the usual parameters defining the perturbation of the instantaneous rotation vector of the Earth $\omega = \Omega_0(m_1, m_2, 1+m_3)$. It is often convenient to write the polar motion in complex notation:

$$ \tilde{m} = m_1 + im_2, \quad (3) $$

where $i^2 = -1$.

Any departure of the observed gravity perturbation from the theoretical value given by eq. (2) may be interpreted as a manifestation of the Earth deformation. In the frequency domain, the ratio of the two perturbations is the gravimetric factor $\tilde{\delta}$:

$$ \tilde{\delta}(\omega) = \frac{\Delta \tilde{g}_m^{\text{obs}}(\omega)}{\Delta \tilde{g}_m^{\text{rig}}(\omega)}, \quad (4) $$

where $\omega$ is the angular frequency. If we write $\tilde{\delta}$, $\Delta \tilde{g}_m^{\text{obs}}$ and $\Delta \tilde{g}_m^{\text{rig}}$ as

$$ \tilde{\delta} = \delta e^{i\kappa} \quad (5) $$

$$ \Delta \tilde{g}_m^{\text{obs}} = A_{\text{obs}} e^{i\psi_{\text{obs}}} \quad (6) $$

$$ \Delta \tilde{g}_m^{\text{rig}} = A_{\text{rig}} e^{i\psi_{\text{rig}}}, \quad (7) $$

we have

$$ \delta = \frac{A_{\text{obs}}}{A_{\text{rig}}} \quad (8) $$

and

$$ \kappa = \varphi_{\text{obs}} - \varphi_{\text{rig}}. \quad (9) $$

If the deformation was purely elastic, $\kappa$ would be zero for any Earth model. Whether the deformation is elastic or not, $\delta$ depends on the Earth model. For instance, for the elastic spherical PREM model (Dziewonski & Anderson 1981), $\delta = 1.16$ if $\omega$ is small and significantly different from any eigenfrequency of the model. Aelasticity would make $\kappa$ different from zero and modify $\delta$. If the polar motion is not corrected for the influence of the oceanic pole tide, its value should increase up to 1.185 in Western Europe (Boy et al. 2000) and induce a phase of a few tenths of a degree (Chen et al. 2008). In several studies devoted to the estimate of the gravimetric amplitude and phase either at individual stations or in a global analysis (Loyer et al. 1999; Xu et al. 2004; Ducarme et al. 2006; Chen et al. 2009), departures from the elastic values have been attributed among other things to either measurements artefacts or geophysical phenomena that were not properly corrected for (e.g. hydrology or oceanic pole tide).

It is reasonable to assume that the gravimetric factor is constant over the duration of the records because it depends on geophysical parameters that vary on geological timescales. However, both the amplitude and phase of $\tilde{m}$ at the Chandler frequency, that is the CW, vary with time (Chao & Chung 2012). For instance, several phase jumps have been observed during the last century (Gibert et al. 1998; Malkin & Miller 2010). The changes in the amplitude are due to attenuation and irregular excitation. Even so, these variations impact both the theoretical and observed gravimetric perturbations in the same way and should not cause any change in the gravimetric factor.

Moving to the time domain, the signal can be modeled by a sum of n sine functions with constant amplitudes $A_j$, phases $\varphi_j$ and real frequencies $f_j$, that is,

$$ \Delta g_m(t) = \sum_{j=1}^{n} A_j \sin(2\pi f_j t + \varphi_j). \quad (10) $$

We can then compute the gravimetric factor amplitude and phase for a given frequency $f_j$ using formulae equivalent to the ones in the frequency domain:

$$ \delta(f_j) = \frac{A_j^{\text{obs}}}{A_j^{\text{rig}}} \quad (11) $$

and

$$ \kappa(f_j) = \varphi_j^{\text{obs}} - \varphi_j^{\text{rig}}. \quad (12) $$

For harmonic motions, it is obvious that the two definitions (8) and (11) are identical, whereas for excited, damped, anharmonic or any slightly more complex motion, the equation in time domain gives a practical way of calculating the gravimetric factor that may differ from the first one. Albeit $A_j$ and $\varphi_j$ are considered constants to be determined, one knows that they vary with time. The observed Chandler period is unknown too, because it may be distinct from the
<figure><img src="image_2.png" /><figcaption>Figure 1. Pole coordinates (top) and corresponding gravity perturbation of a rigid Earth model given by eq. (2) (bottom) at Strasbourg J9 gravimetric station.</figcaption></figure>

eigenperiod defined by the Earth's geophysical properties (Bizouard et al. 2011). Then, the period should just be one of the parameters that have to be adjusted. Furthermore, doing so, we will fit a hypothetical linear variation of the phase with time, since, for a given wave, it is equivalent to setting a constant phase with a modified frequency as shown in the following equality:

$$ \sin[2\pi f_j t + (\varphi_j^1 t + \varphi_j^0)] = \sin[(2\pi f_j + \varphi_j^1)t + \varphi_j^0]. \quad (13) $$

It is more complicated to take the amplitude variation into account, even in the case of a simple linear one, for numerical reasons: it is both very slow and unstable to fit sinusoids and cosinusoids with varying amplitude. That is why we will assume that all the waves have constant amplitudes. Based on the obtained misfits, this assumption will be shown a posteriori to be reasonable. In addition, as previously noted, even if the amplitudes vary in the two time-series, they must have identical variations which should not impact the gravimetric factor estimation.

# 3 POLAR MOTION AND GRAVITY DATA

## 3.1 IERS data

The Earth rotation data are provided by the International Earth Rotation and Reference Systems Service (IERS) as a set of six parameters with their uncertainties: the five standard Earth orientation parameters (EOPs), namely the celestial pole offsets $dX$ and $dY$, the perturbation of the rotation $UT1 - UTC$ and the pole coordinates $x_p$ and $y_p$, plus the length of day variation $\Delta LOD$. Details about their physical meaning and their measurement can be found in McCarthy & Petit (2003) and Bizouard (2014). In this work, we use the EOP 08 C04 time-series, started in 1962 and given at a daily frequency.

To compute the gravimetric perturbation for a rigid Earth as given by eq. (2), we need to express the three parameters ($m_1, m_2, m_3$)

as a function of the EOPs. Actually, $m_3$ is related to the variation of the angular velocity of the rotation of the Earth and its coupling with $m_1$ or $m_2$ compared to the coupling between $m_1$ and $m_2$ makes it negligible in the study of polar motion, as $m_3$ is two orders of magnitude smaller as indicated by Loyer et al. (1999).

The relation between $\tilde{m}$ as defined in eq. (3) and the EOPs was first established by Brzeziński & Capitaine (1993):

$$ \tilde{m} = p - i \frac{\dot{p}}{\Omega_0} + i \frac{\dot{p}}{\Omega_0} e^{i \Omega_0 t} \quad (14) $$

using the complex number $p = x_p - iy_p$ where $x_p$ and $-y_p$ are the coordinates of the celestial intermediate pole (CIP) as seen from the International Terrestrial Reference Frame and $P = dX \sin \epsilon_0 + idY$ is the pole offset in obliquity and longitude as seen from the International Celestial Reference Frame, with $\epsilon_0$ the mean obliquity of reference. The CIP is always less than a few centimetres apart from the instantaneous pole of rotation we are interested in.

Eq. (14) fills the gap between a purely geometrical description of polar motion and dynamic theories of Earth rotation that rely on the instantaneous rotation vector. For high-frequency motions, it is necessary to use this rigorous relation as pointed out by Loyer et al. (1999), but for low-frequency phenomena such as the CW or annual oscillation, the error in the computation of $\Delta g_m^{rig}$ is below 0.5 per cent, that is, below the accuracy we can reasonably achieve. This assertion has been verified by comparison of our results with or without using eq. (14), the difference on $\Delta g_m^{rig}$ being below 0.01 $\mu$Gal at mid-latitude. Thus, we simply use the relation

$$ \tilde{m} = p. \quad (15) $$

Fig. 1 shows the pole coordinates ($x_p, y_p$) between 1996 August and the end of 2009. Also shown is the corresponding gravity perturbation of a rigid Earth model given by eq. (2) at the Strasbourg J9 station. A polar motion of a few hundred milliarcseconds causes a radial gravity variation of at most a few microGals at mid-latitudes.
The 6.5 yr beating between the annual component and the CW is
clearly visible.

## 3.2 GGP data

The gravimetric data used in this work come from the Global Geo-dynamics Project (GGP; Crossley *et al.* 1999), now disrupted to be integrated in the newly created International Geodynamics and Earth Tides Service. Established in 1997, the GGP network was made of about 35 gravimetric stations that have provided decades of gravity recordings. Because of the closeness of the CW and annual signals in the frequency domain, only the uninterrupted time-series longer than the beating period of the two oscillations are useful for this study.

The hourly data are prepared by the International Centre for Earth Tides (ICET) and provided by the Information Systems and Data Center of the Geo Forschungs Zentrum. The local barometric changes at each station are also included to allow for local atmospheric pressure correction. All of the gravimetric data sets contain gaps and offsets that must be carefully managed because of their strong impact on low-frequency components of the entire signal. Superimposed on these artefacts and gravity variations, a linear drift of a few microGals per year is visible on recent time-series, whereas an exponential drift impacts the oldest and longest recordings. The removal of these drifts is included in the data processing described in the next subsections.

### 3.2.1 Correction for tides, atmospheric pressure and surface loadings

The first correction applied to the gravimetric data is the removal of the tidal signal, which has an amplitude up to several hundreds of microGals. However, the estimation of the gravimetric factor is barely modified whether the solid tides are removed or not. Indeed, the largest tidal components are high-frequency waves that will be filtered out and the low-frequency waves (semi-annual, annual, etc.) will be included in the fitting by means of eq. (10).

Second, we correct the gravity data for the local pressure effects by converting the pressure variations into gravimetric variations using an admittance factor of $-3 \text{ nm s}^{-2} \text{ hPa}^{-1}$.

Third, we use the loading data provided by the GGP Loading Service to correct for various loading effects. The Loading Service corrections encompass non-local barometric effects (Boy *et al.* 2002), local and non-local hydrological loading (Boy & Hinderer 2006) and oceanic loading (Boy *et al.* 2009).

Fourth, we correct the gravity data for the oceanic pole tide, this time-series having been added in the Loading Service products meanwhile (2015 February). The effect of the ocean tide loading caused by polar motion is computed on the assumption that the ocean response is static, which is consistent with altimetric observations (Desai 2002). However, it has been shown that the pole tide is dynamic and significantly larger in shallow seas such as the North and Baltic Seas (Miller & Wunsch 1973; Dickman 1979, 1988). Several stations used in this study are located in Europe and are probably affected by the dynamic oceanic pole tide, leading to an underestimation of the pole tide correction. Even so, the computation of the dynamic loading due to the pole tide at each GGP station is far beyond the scope of this study. Rough estimates of this effect by an artificial multiplication of the pole tide by factors between 2 and 10 confirmed that the increase of the correction by a given factor may significantly decrease the value of the estimated

gravimetric factor but no precise computation can be made without further modeling.

Fig. 2 illustrates the gravimetric corrections that we apply to the GGP data.

### 3.2.2 Offsets removal

The next step in the data processing is the proper correction of the offsets. Their influence on the study of the long-period phenomena has been pointed out, for example, by Xu *et al.* (2004). The strategy we have adopted in this work is to minimize iteratively the misfit between a model of the main polar motion components and the gravity residues. In the following order, we

(i) list the offsets: they are supposed to occur when (1) at least one data point is missing or (2) there is no real interruption but an obvious offset is apparent,

(ii) make a first rough offset correction using the median value of a few points on each side of the offset or gap; depending on the stations, between 8 and 30 points (i.e. hours) are used,

(iii) loop until a chosen criterion is met (stabilization of the rms or maximum number of iterations is reached, either having been used here) using the following steps:

(a) Fit of sinusoids with the following periods to the signal (see Section 3.2.4): 1 yr, half a year and the estimated CW period, plus a linear trend.

(b) Use of Latin Hypercube Sampling (LHS) method from McKay *et al.* (1979) to randomly sample new sets of offsets corrections (the times are fixed and we only look for offsets amplitudes).

(c) Application of these offsets corrections and computation of the corresponding rms (misfit to the model), to both the data and the theoretical perturbation deduced from polar motion data.

(d) Addition of the two rms for each tested set and search for the best combined rms to find the best offset correction for this iteration.

(e) Addition of the best correction to the total offset correction only if the new best combined rms is smaller than the previous one.

(f) Exit of this loop if the chosen criterion is met or return to (a).

(iv) application of the best offset correction to the residues.

The fitting procedure being detailed in Section 3.2.4, for now we just fit a few sinusoids of intermediate frequencies only to avoid an overfitting of the data at long period: 1, 0.5 yr and the CW period estimated from the polar motion data. Indeed, if there are many offsets in the time-series and we try to fit a long-period wave, the different segments could be artificially shifted to match this low-frequency oscillation of potentially high amplitude. The tides at 1305 and 205 d were initially fitted too, following Loyer *et al.*'s (1999) procedure, but considering their extremely small amplitudes which did not affect the results, they are not considered anymore.

The LHS method (McKay *et al.* 1979) used to determine the successive offsets corrections consists in a random sampling in the space of parameters, following a given random distribution. In other words, several complete sets of offset are randomly generated at each iteration and only the best one—the definition of this ‘best’ is given in the next paragraph—is finally applied. Here we use either a uniform or a Gaussian distribution for the sampling, depending on the confidence we have in the first rough correction from the median computation around offsets. For offsets associated with gaps of a
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Gravimetric corrections applied to the raw gravity time-series at the Strasbourg station. The local atmospheric correction is computed from the local pressure time-series at each station with an admittance factor of $-3.0668 \text{ mm s}^{-2} \text{ hPa}^{-1}$. The atmospheric effects along with the local and non-local hydrological corrections and oceanic loading are computed from the European Centre for Medium Range Weather Forecasts (ECMWF) reanalysis (ERA interim) surface pressure fields, assuming for the atmospheric loading computation a barotropic ocean model forced by air pressure and winds (MOG2D) when provided at the time of gravity records, or an inverted barometer ocean response otherwise.</figcaption></figure>

few weeks or more, it is more efficient to consider that all offsets values are equiprobable within a specific range. Although alternative algorithms exist, the LHS method is a good option as it is simple to implement and a fast way of trying a lot of possible corrections without any prior assumptions.

To determine the best offset correction to be applied, we add the rms from both the data and theoretical gravity perturbation. The use of a combined rms is an efficient way of avoiding overfitting of anomalous signals in the gravity data or introducing artificial non-reasonable offsets. However, the observed and theoretical signals have different amplitudes. Thus, they cannot be directly compared to the same model. Doing so, we would introduce a systematic bias in the difference between the theoretical gravity time-series and the model, which is only fitted on the gravity data. To avoid this problem without making any assumption on the amplitudes difference, we simply normalize both signals by their maximum amplitude. Then, we can compute consistent misfits between the fitted model on the one hand and observed or theoretical gravity perturbation on the other hand.

### 3.2.3 Filtering and decimation

After the offset correction, we linearly interpolate the data to fill the gaps smaller than one month and we do not include in the study the stations providing data with larger gaps.

To avoid to estimate many high-frequency waves, to eliminate short-period artefacts and to decrease the number of data points, we filter the signal using a low-pass finite impulse response (FIR) filter with no phase shift and decimate it to one day. The absence

of phase shifting is essential to preserve the phase of the signal we are interested in. With this filter, the same number of data points both at the beginning and end of the signal are lost. A happy consequence of this loss is the systematic elimination of possible quick drift at the beginning of some time-series. An ideal filtering would eliminate any signal above the CW frequency but because most of the gravimetric time-series do not exceed two decades, it is impracticable. An empirical estimation indicates that a filter length of roughly 30 yr—inducing a loss of 15 yr both at the beginning and end of the signal—is necessary to fully eliminate the annual signal without affecting the CW. Synthetic tests have shown that even with an attenuation as weak as $-20 \text{ dB}$ at the annual period and no significant impact at the CW period, a 10-yr long filter is required with a rectangular tapering window. Therefore, we have chosen to follow the filtering procedure used by Loyer *et al.* (1999), with a cut-off frequency at 90 d. The length of our filter is 2 yr, which means that 1 yr of signal is lost at the beginning and end.

### 3.2.4 Fitting

The last step of the data processing is the estimation of the sinusoids that compose the residual and theoretical signals. Their amplitudes and phases will be used to compute the gravimetric factor given by eq. (11). To do so, we fit a sine and a cosine of different amplitudes for each frequency, from which we deduce the corresponding amplitude and phase of the equivalent shifted sinusoidal wave as expressed in eq. (10). The periods of the forced oscillations are 18.6, 9.3, 1 and 0.5 yr. We also include in this set the observed CW period. Apart from the annual, semi-annual and Chandlerian
period, the others correspond to the waves with largest amplitudes at long period as indicated in Loyer et al. (1999). As a linear time variation of the phase of the signal appears as a time variation of the Chandler frequency, we first estimate the apparent CW frequency by means of the rotation data of the IERS and use this value in the fitting procedure for the residual time-series.

Because mantle anelasticity can only produce a negative phase lag, as reminded by Ducarme et al. (2006), a constraint on the fitting must be added to force a negative phase for the gravimetric factor when needed. This constraint is similar to the positivity constraint of the quality factor introduced by Florsch & Hinderer (2000) in the case of the Free Core Nutation, which is another free rotational mode of the Earth. For some of the stations, the fitting actually yields a positive phase, which does not respect the causality principle and has no physical meaning. In those cases, we have first estimated the phase of the Chandlerian oscillation from the polar motion data, then use it as an upper limit for the phase of the Chandlerian signal in gravity residuals, in order to keep $κ$ as given by eq. (9) negative.

Another component to fit is the instrumental drift. Although the time-series of the old SG (Metsahovi (ME) station) had an exponential drift, recent SG recordings (all of the other series used in this study) have a linear drift, except for a transient initial drift eliminated by the FIR low-pass filter.

For the fitting, we used the routines included in the lmfit Python package (Newville et al. 2014) with the least-squares method. They provide uncertainties estimates for each parameter we use with other methods described in Section 3.2.5 to estimate the reliability of our results.

### 3.2.5 Convergence and uncertainties estimate

To estimate the uncertainties associated to our results, we apply the fitting method to time-series whose length is progressively increased until they contain the whole data sets and we plot the amplitude and phase of the gravimetric factor as a function of the time span of the data. They should theoretically converge toward constant values for time-series that are longer than the beating period of 6.5 yr, although most of them do not completely as discussed in Section 6. This method is qualitatively used by Loyer et al. (1999). We use the variance of the amplitude and phase values after the 6.5 yr limit to estimate the uncertainties of the gravimetric factor. We denote the standard deviation obtained from the convergence by $σ^{CV}$, as opposed to the uncertainty $σ^{fit}$ obtained from the fitting, which is the formal least-squares error.

For the stations where the phase had to be constrained to be negative, the formal error is frequently equal to the phase itself, which means that the fitting could not provide the best solution without exceeding the imposed limit. In other words, the best estimate for the gravimetric factor phase could be zero for these stations, considering the physical constraint. In Tables 2 and 3, the affected stations are indicated with a formal error between parentheses.

When we consider the variability in the results reported in Section 5, using the convergence of the computed values obviously seems to yield a better quantification of the uncertainty than the usual misfit. Indeed, this method takes advantage of the intrinsic variability (whatever its origin) of gravity records with time to put less tight but more realistic constraints on the gravimetric factor estimated from GGP data. Of course, this comment only holds when the phase is negative without the need for any constraint.

## 4 STACKING

In the previous section, we have presented a step-by-step processing to remove the $Δg_{geo}$ (gravimetric corrections) and $Δg_{SG}$ (drift and offsets) terms from eq. (1) for any individual GGP time-series in order to extract the polar motion effect. None of these steps have decreased the noise.

A standard approach to increase the signal-to-noise ratio in gravimetric data was first used by Cummins et al. (1991) in their search for oscillatory gravity signals associated with core undertones. The method consists in stacking the data from different stations, assigning them a weight in accordance to the spatial pattern of the sought gravimetric signal. The perturbation of the vertical component of the gravity field associated to the CW is a degree 2, order -1 spherical harmonic, as depicted in Fig. 3.

To adapt the stacking method to polar motion study, we select the GGP stations with time-series longer than 10 yr to ensure a sufficient separation between the CW and annual signals. Moreover, to apply optimal offsets correction, we reject stations or time-series with long interruptions. Last, we need stations with simultaneous records in order to stack them. We therefore exclude time-series which stopped before 2008 (see Table 1). Among the nearly 35 stations of the GGP network, only nine are first retained. They are listed in Table 1 with the start and end dates of the recordings. The new ST time-series actually covers the whole time period between 1997 and today but instrumental problems on the tilts limited the data usability for this study till 2009.

After having processed the data as described in Sections 3.2.1–3.2.3, the time-series $Δg_j(t)$ at station j is the sum of a signal $s_j(t)$ and some uncorrelated noise n:

$$y_j(t) = s_j(t) + n(t). \quad (16)$$

The stacking formula is:

$$ \begin{aligned} S_{\ell}^{m}(t) &= \frac{\sum_{j} \bar{Y}_{\ell}^{m}(\theta_{j}, \phi_{j}) y_{j}(t)}{\sum_{j} |Y_{\ell}^{m}(\theta_{j}, \phi_{j})|^2} \\ &= \frac{\sum_{j} \bar{Y}_{\ell}^{m}(\theta_{j}, \phi_{j}) s_{j}(t) + \sum_{j} \bar{Y}_{\ell}^{m}(\theta_{j}, \phi_{j}) n_{j}(t)}{\sum_{j} |Y_{\ell}^{m}(\theta_{j}, \phi_{j})|^2} \end{aligned} \quad (17) $$

where $Y_{\ell}^{m}(\theta_{j}, \phi_{j})$ is the spherical harmonics of degree $\ell = 2$ and order $m = -1$ computed at station j of coordinates $(\theta_{j}, \phi_{j})$.

In addition to this spatial weighting and to take into account the variability in the quality of the stations, we similarly weighted the individual time-series $y_j(t)$ using the standard deviation on the convergence of $\delta$, namely the $\sigma_{\delta}^{CV}$ appearing in Table 3 for each station. The weights $w_j$ (last column of Table 1) are defined in such a way that they sum up to 1

$$w_j = \left( \sigma_j \sum_k \frac{1}{\sigma_k} \right)^{-1}, \quad (18)$$

where $\sigma_{\delta}^{CV}$ at station j is simply written $\sigma_j$ for the sake of clarity.

Then, we may defined $y_j^{\sigma}(t)$ as the time-series weighted by the individual uncertainties:

$$y_j^{\sigma}(t) = w_j y_j(t) \quad (19)$$

and use it instead of $y_j(t)$ in eq. (17).

Xu et al. (2004) made a global estimation of the gravimetric factor with a different subset of five stations but they did not spatially weight the time-series. Their weighting was only based on the standard deviation of the adjusted amplitudes and phases in the time domain. In addition, they could only work on a common time span
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Polar motion induced elastic gravity perturbation at randomly chosen date. The stations used in this paper are symbolized by white stars. At this date, the pole position was $x_p = 0.182947$ arcsec and $y_p = 0.283920$ arcsec. In this picture, the degree 2, order −1 spherical harmonic pattern used to weight the time-series in the stacking clearly appears.</figcaption></figure>

<table><caption>**Table 1.** First selection of GGP stations exhibiting the longest duration available. Columns Long. and Lat. are the stations coordinates in decimal degrees in longitude and latitude, respectively. The last column provides the weights as defined in eq. (18). CA and ST are limited to the new SG time-series (series from older instruments are available but they do not add any value in the stacking because they are noisier). For the ST station, we have used a time-series already processed for previous studies, hence the absence of gaps and offsets.</caption><thead><tr><th>Location (name)</th><th>Long.</th><th>Lat.</th><th>Start</th><th>End</th><th>Duration (yr)</th><th>Offsets</th><th>Weights</th></tr></thead><tbody><tr><td>Bad Homburg (BH)</td><td>8.6113</td><td>50.2285</td><td>2001.1</td><td>2012.3</td><td>11.2</td><td>16</td><td>0.114</td></tr><tr><td>Cantley (CA)</td><td>284.1927</td><td>45.5850</td><td>1997.5</td><td>2012.0</td><td>14.5</td><td>9</td><td>0.019</td></tr><tr><td>Canberra (CB)</td><td>149.0077</td><td>-35.321</td><td>1997.0</td><td>2012.0</td><td>14.9</td><td>39</td><td>0.340</td></tr><tr><td>Matsuhiro (MA)</td><td>138.2032</td><td>36.5439</td><td>1997.3</td><td>2008.5</td><td>11.2</td><td>66</td><td>0.118</td></tr><tr><td>Membach (MB)</td><td>6.0066</td><td>50.6093</td><td>1995.5</td><td>2012.0</td><td>16.4</td><td>27</td><td>0.045</td></tr><tr><td>Medicina (MC)</td><td>11.645</td><td>44.5219</td><td>1998.0</td><td>2012.0</td><td>14.1</td><td>8</td><td>0.118</td></tr><tr><td>Metsahovi (ME)</td><td>24.3958</td><td>60.2172</td><td>1994.6</td><td>2012.5</td><td>18.0</td><td>71</td><td>0.070</td></tr><tr><td>Moxa (MO)</td><td>11.6156</td><td>50.6447</td><td>2000.0</td><td>2012.0</td><td>12.0</td><td>26</td><td>0.106</td></tr><tr><td>Strasbourg (ST)</td><td>7.6850</td><td>48.6217</td><td>1997.1</td><td>2009.0</td><td>12.4</td><td>-</td><td>0.070</td></tr></tbody></table>

of nearly six years for the selected stations. Using a different approach, Ducarme et al. (2006) provided another global estimation of the gravimetric factor using a larger subset of nine stations but their goal was not to directly address the question of the signal-to-noise ratio or to get rid of local effects. The combination of multiple stations was rather used to increase the length of the data set they inverted, with redundancies in time. Here, we have chosen to first exploit the spatial and temporal dependency between the data by weighting and stacking them to decrease the influence of spatially incoherent signals and noise before estimating the global parameters. Last but not least and as already noted, another important feature of this study is the long duration of the time-series used in the stacking, the time-series sharing a common time span longer than a decade.

## 5 RESULTS

The main aim of this study being to estimate a global gravimetric factor from multiple stations, we only focus on the stations listed in Table 1. The results are summarized in Table 2. The values of the uncertainties $\sigma_{\delta}^{\text{fit}}$ and $\sigma_{\kappa}^{\text{fit}}$ are deduced from the misfit between the

fitted model and the residual filtered gravity, whereas the values of $\sigma_{\delta}^{\text{CV}}$ and $\sigma_{\kappa}^{\text{CV}}$ are estimated from the convergence of the parameters values for increasing length of the time-series: it is the standard deviation of the estimated values for the period ranging from the 6.5 yr threshold until the end of the time-series. An easy way to estimate the quality of the convergence is thus to compare the two kinds of uncertainties in Table 2 for the amplitude $\delta$ and phase $\kappa$. As an example, Fig. 4 illustrates the convergence obtained for one of the best (MC) and for one of the worst (CA) stations in terms of the stability of the results. It is noteworthy that an increase of the time-series length is in any case associated with variations of the gravimetric factor amplitude and phase and we will discuss this point further in Section 6.

Because the quality of the hydrology correction is one of the most difficult to assess, it is only applied when the obtained amplitude is closer to the expected elastic values. It is notable that both the amplitude and phase are impacted by the hydrology correction when applied.

Moving to our multistation analysis, the results for the stacking method are given in Table 3 for two different subsets of, respectively, five and seven SG stations. The selection of the stations in these
<table><caption>**Table 2.** Gravimetric factors amplitude and phase for the selected stations. The asterisk indicates that the hydrology correction has been applied.</caption><thead><tr><th>Location (name)</th><th>δ</th><th>σ<sub>δ</sub><sup>fit</sup></th><th>σ<sub>δ</sub><sup>CV</sup></th><th>κ(deg)</th><th>σ<sub>κ</sub><sup>fit</sup></th><th>σ<sub>κ</sub><sup>CV</sup></th></tr></thead><tbody><tr><td>Bad Homburg (BH)*</td><td>1.165</td><td>0.020</td><td>0.026</td><td>-0.39</td><td>(0.39)</td><td>0.03</td></tr><tr><td>Cantley (CA)*</td><td>1.315</td><td>0.030</td><td>0.177</td><td>-21.70</td><td>1.09</td><td>4.25</td></tr><tr><td>Canberra (CB)</td><td>1.066</td><td>0.013</td><td>0.010</td><td>-4.83</td><td>0.55</td><td>1.12</td></tr><tr><td>Matsuhiro (MA)</td><td>1.042</td><td>0.020</td><td>0.028</td><td>-0.20</td><td>(0.20)</td><td>0.02</td></tr><tr><td>Membach (MB)</td><td>1.145</td><td>0.022</td><td>0.075</td><td>-0.66</td><td>0.92</td><td>1.87</td></tr><tr><td>Medicina (MC)</td><td>1.103</td><td>0.011</td><td>0.034</td><td>-0.29</td><td>(0.29)</td><td>0.05</td></tr><tr><td>Metsahovi (ME)</td><td>0.948</td><td>0.021</td><td>0.049</td><td>-0.26</td><td>(0.26)</td><td>7.00</td></tr><tr><td>Moxa (MO)</td><td>0.967</td><td>0.016</td><td>0.020</td><td>-0.38</td><td>(0.38)</td><td>2.78</td></tr><tr><td>Strasbourg (ST)*</td><td>1.225</td><td>0.017</td><td>0.049</td><td>-2.32</td><td>0.66</td><td>2.31</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Convergence of the gravimetric factor amplitude and phase for increasing length of the time-series. Part (a) is for Medicina station (MC) and part (b) is for Cantley (CA). The vertical line indicates the 6.5 yr limit theoretically needed to separate the annual and Chandlerian oscillations; the horizontal red dashed line is the mean of the values after the 6.5 yr limit has been exceeded.</figcaption></figure>

<table><caption>**Table 3.** Gravimetric factors amplitude and phase for the stacking and multistation analysis. Seven stations: BH, CA, CB, MB, MC, ME and MO; five stations: BH, CB, MB, MC and MO. Ducarme *et al.* (2006) estimated a global amplitude with phases either set to the individual values or to zero for all of the stations. Their estimated phase is -1 d, which converts to -0.83 deg for a Chandler period of 435 d.</caption><thead><tr><th>Analysis</th><th>δ</th><th>σ<sub>δ</sub><sup>fit</sup></th><th>σ<sub>δ</sub><sup>CV</sup></th><th>κ(deg)</th><th>σ<sub>κ</sub><sup>fit</sup></th><th>σ<sub>κ</sub><sup>CV</sup></th></tr></thead><tbody><tr><td><strong>This study</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Spatial and uncertainty weighting</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Stack 7 stations</td><td>1.092</td><td>± 0.010</td><td>± 0.017</td><td>-0.32</td><td>(± 0.32)</td><td>± 0.01</td></tr><tr><td>Stack 5 stations</td><td>1.118</td><td>± 0.009</td><td>± 0.016</td><td>-0.45</td><td>± 0.30</td><td>± 0.66</td></tr><tr><td>Spatial weighting only</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Stack 7 stations</td><td>1.090</td><td>± 0.004</td><td>± 0.019</td><td>-0.36</td><td>(± 0.36)</td><td>± 0.01</td></tr><tr><td>Stack 5 stations</td><td>1.090</td><td>± 0.010</td><td>± 0.023</td><td>-0.34</td><td>(± 0.34)</td><td>± 0.03</td></tr><tr><td>Mean—see eqs (20) and (21)</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Mean 7 stations</td><td>1.079</td><td>± 0.006</td><td>± 0.008</td><td>-1.04</td><td>± 0.15</td><td>± 0.026</td></tr><tr><td>Mean 5 stations</td><td>1.081</td><td>± 0.007</td><td>± 0.009</td><td>-0.85</td><td>± 0.18</td><td>± 0.026</td></tr><tr><td><strong>Xu et al. (2004)</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>STD weighted</td><td>1.1613</td><td>± 0.0275</td><td>—</td><td>-1.30</td><td>± 1.33</td><td>—</td></tr><tr><td><strong>Ducarme et al. (2006)</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Ocean pole tide not corrected</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Mean</td><td>1.1788</td><td>± 0.0040</td><td>—</td><td>-0.83</td><td>—</td><td>—</td></tr><tr><td>Global (local phase)</td><td>1.1816</td><td>± 0.0047</td><td>—</td><td>local</td><td>—</td><td>—</td></tr><tr><td>Global (zero phase)</td><td>1.1797</td><td>± 0.0047</td><td>—</td><td>0</td><td>—</td><td>—</td></tr><tr><td>Ocean pole tide corrected</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Mean</td><td>1.1605</td><td>± 0.0040</td><td>—</td><td>-0.83</td><td>—</td><td>—</td></tr><tr><td>Global (local phase)</td><td>1.1612</td><td>± 0.0047</td><td>—</td><td>local</td><td>—</td><td>—</td></tr><tr><td>Global (zero phase)</td><td>1.1593</td><td>± 0.0047</td><td>—</td><td>0</td><td>—</td><td>—</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>Figure 5. Correlation between the annual and CW sinusoidal fitted amplitudes for increasing length of the time-series at station MB. The vertical line indicates the 6.5 yr limit. We can clearly see the pseudo-oscillation in the correlation.</figcaption></figure>

subsets is discussed in Section 6. The weighted mean E is given by:

$$E = \sum_j w_j y_j \qquad (20)$$

where the weights $w_j$ are given by eq. (18) in which $\sigma$ is either $\sigma^{\text{fit}}$ or $\sigma^{\text{CV}}$ for the amplitude or phase. The error $e$ associated with the mean is

$$e = \left( \sum_j \frac{1}{\sigma_j^2} \right)^{-\frac{1}{2}}. \qquad (21)$$

For increasing time spans, we have plotted the variations of the correlation between the fitted parameters for the annual and Chandlerian oscillations (Fig. 5). As expected, this correlation quickly decreases toward zero for time-series longer than the minimum duration needed to theoretically separate the annual signal from the CW. Even so, this decrease is not monotonic and exhibits oscillations of pseudo-half-period roughly similar to the 6.5 yr beating. This value is certainly related to the frequency difference between the two main signals contained in the data but its emergence in this case has yet to be explained. Similar fittings on synthetic signals exhibit this behaviour too, which would therefore be caused by a numerical artefact.

The main results for the individual stations and the stacking are summarized in Fig. 6.

# 6 DISCUSSION

If we first focus on the individual estimates of the gravimetric factor, the most obvious observation is that none of the amplitude or phase is easier to constrain and in particular some stations have outliers, for the amplitude, phase or both. As aforementioned, one of the source of errors lies in the offsets, especially when they are associated with other instrumental defects that prevent a simple realignment of the data points on either side of the offsets. Synthetic tests have shown that although the phase is actually sensitive to the offset correction and any long-period phenomenon that affect the trend of the signal, the amplitude is equally impacted by offsets too. As an

example, when the phase variation due to synthetic offsets reaches a few tenths of a degree, the amplitude undergoes variations of a few per cent.

Considering the obvious outliers, the explanations are essentially assumptions because of the diversity of possible sources of errors. For ME, the abnormally small amplitude may be due to the large number of offsets we had to correct taking into account an exponential drift, along with the position of the station at relatively high latitude where the polar motion effect is less important. The proximity of the ocean is another hint considering the additional noise. This last suggestion could be an explanation for the abnormally high amplitude at Cantley (CA) station which is close to the sea too. For Moxa (MO), the hydrology is probably one of the main explanation for the small amplitude. This hypothesis is suggested by the opposite way the amplitude and phase vary when we change the admittance factor used to compute the local hydrological effect at this station. Indeed, depending on the correction, an increase (improvement in this case) in the amplitude is accompanied by an increase of the phase which was already too large when we did not constrain it to be negative. Of course, if the Chandlerian and annual signals could be completely decorrelated—which is not the case here—the influence of hydrology would not be so strong because the hydrological effects are quite small at the Chandler frequency.

Amongst the gravimetric corrections we made on the time-series, the local hydrological correction has a particularly strong impact on the estimation of both the gravimetric factor amplitude and phase and on the convergence of their values in our tests. This is a serious issue because this correction is hard to improve without appropriate local hydrological studies and modeling. Actually, a lot of parameters and phenomena have an influence on the effect of hydrology on gravity measurements, the position of the station relative to the ground, the distribution of underground water masses both in the vadose and saturated zones being only some of them. For only some of the stations, the convergence is improved when the local hydrological correction is not applied (MC, for instance) but more importantly, both the amplitude and phase of the gravimetric factor are not clearly improved, if not worsened when hydrology is taken into account. A systematic determination of the local hydrological admittances at each GGP station would be the first step in an attempt to better consider hydrology in the gravimetric correction.

Another parameter which influences the gravimetric factor phase is of course the Chandler frequency. Everything else being equal, an artificial decrease of the Chandler period induces an increase of the gravimetric factor amplitude, which could be qualitatively interpreted as a tendency for the fitting to maintain the overall energetic content of the CW signal. Anyway, contrary to the offsets correction or local hydrology, the influence of a Chandler frequency variation is quite predictable, at least at the first order and reasonably small compared to other sources of error. The observed variations in polar motion data are indeed smaller than the spectral resolution we have in gravimetry.

Focusing on the stacking now, it is important to be aware that combining different stations does not necessarily improve the overall analysis quality. Using only a station with very good corrections clearly provides better results than stacking poorly corrected gravity time-series. Out of the nine stations we had initially selected, we did not retain MA and ST because of their usable time-series ending too early. Thus, we only kept the seven stations BH, CA, CB, MB, MC, ME and MO. Out of those these stations, we then only kept BH, CB, MB MC and MO for their reasonably small phase values and good convergence (small $\sigma_k^{\text{CV}}$).
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Gravimetric factor amplitude and phase for the individual stations and the stacking with five time-series. The stations name preceded by an asterisk are included in the five stations stacking. The vertical red dashed lines are the elastic values, namely $\delta = 1.16$ and $\kappa = 0$ (see Section 6). The error bars are the uncertainties from the convergence. The red diamonds for the stackings are the results with weights depending on the quality of each data set.</figcaption></figure>

After the stacking of the five best stations we have in the set, albeit the estimated amplitude is not completely satisfying, its uncertainty is smaller than for all of the individual estimates with a better convergence, that is, less sensitivity to the time-series length. The phase is even better recovered with quite a low uncertainty, considering the phase disparity in the individual estimates. Amongst the local effects we had to correct, the stacking could be especially efficient as for local hydrology. An interesting observation is the improvement of the gravimetric factor estimate when we select the five best stations amongst the initial set. It tends to show that we can take profit of the number of GGP stations to only extract a smaller subset of valuable data depending on the processes we are interested in.

We have already underlined that for none of the station the convergence of both the amplitude and phase values is totally achieved, even with the longest time-series. Of course we expected some variations after the theoretical 6.5 yr limit but the observed variations are quite significant for most of the stations compared to the sought accuracy. As far as we know, none of the results provided in the literature for the gravimetric factor, but Loyer et al. (1999), explicitly consider the influence of the time span of the gravimetric time-series. All of the previous works implicitly assumed that the longest the time-series, the better the estimate. Alternatively and as mentioned in Section 3.2.5, we consider longer time-series as a better way of estimating the uncertainties on the gravimetric factor since one of the biggest source of variations both in the amplitude and phase lies in the selection of the time span, revealing the need for even better corrections.

As already stated, the main application of the CW gravimetric factor estimation is to put some constraints on mantle anelasticity, especially through the computation of the expected gravimetric factor phase for classical rheologies of the mantle. Although the actual

computation of this phase for realistic Earth's models is beyond the scope of this work, we briefly discuss the method.

First, the complex gravimetric factor as a function of frequency $\omega$ is given by:

$$ \delta(\omega) = 1 - \frac{3}{2} \tilde{k}_2(\omega) + \tilde{h}_2(\omega) \quad (22) $$

with $\tilde{k}_2(\omega)$ and $\tilde{h}_2(\omega)$ being the complex Love numbers. The frequency dependency of $\tilde{k}_2$ and $\tilde{h}_2$ arises from the mantle anelasticity, characterized by the viscosity. Then, because the Love numbers are viscosity-dependent, eq. (22) provides a link between the viscosity and the gravimetric factor phase. To link these parameters to the quality factor $Q$, we may defined the parameter $\varepsilon$ as:

$$ \tilde{k}_2 = k_2(1 + i\varepsilon) \quad (23) $$

and use the following relation (Bizouard 2014):

$$ \varepsilon = \frac{1 - k_s / k_2}{2Q} \quad (24) $$

where $k_s$ is the secular Love number for the chosen Earth's model.

Using eqs (22) and (24), it is then possible to plot the phase $\kappa$ as a function of viscosity (because the Love numbers are functions of the viscosity) and estimate the corresponding quality factor $Q$. Whatever parameter is considered to be known among the phase, viscosity and quality factor, we can determine the two others for a given rheological model. Alternatively, if we know both the phase and viscosity or phase and quality factor, we can investigate the ability of different rheological models to explain the observations. In concrete terms, using the phases estimates from Table 1, we may put some constraints on the rheological behaviour of the mantle, as Benjamin et al. (2006) did for different geodetic observations, including the CW.
Regarding the amplitude of the gravimetric factor, Dickman (2005) suggested the need to use rotationally consistent Love numbers, which are about 10 per cent larger than the accepted values when the lack of core–mantle coupling is taken into account in the elastic deformation. It would imply an increase of 1–2 per cent of the elastic gravimetric factor amplitude, which is related to the Love numbers through eq. (22). The new elastic value would be closer to 1.18 and the interpretation of our results and those from previous works should be adapted accordingly.

Using the stacking method, we have obtained smaller uncertainties than previous global estimates (Xu *et al.* 2004) of the amplitude and phase of the gravimetric pole tide. However, these uncertainties are still too large to bring some constraints on the mantle rheological behaviour as just explained. So we may conclude that the influence of anelastic effects at the Chandler frequency is still too small to be proven from such gravity time-series analysis. None of the departures of the gravimetric factor amplitude and phase from the theoretical elastic values can be considered relevant, taking into account reasonable uncertainties on our results. On the contrary, this work neither excludes the existence of anelasticity even if its impact is beyond our detection yet.

## 7 CONCLUSIONS

Our estimates of the gravimetric factor at the Chandler frequency using both local and global analyses have shown the actual interest of the stacking method. In particular, local hydrological perturbations that may be so harmful in local studies seem to be significantly reduced when weighting and combining the time-series from several stations. Another essential component of this work lies in the length of the data sets we have used: more than a decade of simultaneous gravimetric recordings were included in the stacking, far more than the theoretical 6.5 yr limit needed to separate the annual and Chandlerian peaks. We have taken advantage of these long time-series to better estimate the uncertainties which are higher but more reasonable than the usual formal estimates from the misfit between gravity residuals and the fitted model.

The gravimetric factor amplitudes and phases we have obtained, either for individual stations or after the stacking, differ notably from previous estimates. These differences mainly come from four specific factors whose influence seems to be underestimated in previous works. The first one is the presence of offsets and gaps which are difficult to correct without introducing a bias in the study. Their effect cannot be easily determined considering the way they may compensate or cumulate. Synthetic tests have shown that the offsets can impact both the amplitude and phase, the latter being more sensitive to them. The second factor is the lack of proper correction for local hydrology and dynamic pole tide, both processes being quite difficult to model at a good level of accuracy. Once again, they strongly affect both the amplitude and phase. Local hydrology might be better modeled in the future at each GGP station but for the moment, it is not well enough known to be routinely and reliably used without taking the risk of introducing an additional bias. The third factor is related to the duration of the time-series which was only considered by Loyer *et al.* (1999) before this work. Indeed, the biggest variations still come from the choice of the data subsets or from their duration as shown by the convergence tests. This comment clearly proves the need to use very long, uninterrupted time-series to infer the reliability of any method or processing in the estimation of the gravimetric factor. The last factor is the geographic distribution of the stations, in relation with the stacking

method. Most of the long time-series being currently recorded in Europe, the advantages of the stacking are partly canceled by this highly heterogeneous distribution. However, by ensuring continuous recordings and adding new stations in the GGP network, this specific problem could be patiently solved.

As for the link between anelastic dissipation and the gravimetric factor phase, we have shown that our estimates do not exclude anelastic effect at the CW period although further refinements are needed to definitely answer this question.

## ACKNOWLEDGEMENTS

We would like to express our thanks to GGP managers for sending their data, to ICET for preparing them, to S. Dickman for his fruitful review and to an anonymous reviewer. We gratefully acknowledge financial support from the CNRS-INSU Programme National de Planétologie and from IPGS.

## REFERENCES

* Aoyama, Y., Naito, I., Iwabuchi, T. & Yamazaki, N., 2003. Atmospheric quasi-14 month fluctuation and excitation of the Chandler wobble, *Earth Planets Space*, **55**, 25–28.

* Bellanger, E., Le Mouël, J.-L., Mandea, M. & Labrosse, S., 2001. Chandler wobble and geomagnetic jerks, *Phys. Earth planet. Inter.*, **124**, 95–103.

* Benjamin, D., Wahr, J.M., Ray, R.D., Egbert, G.D. & Desai, S.D., 2006. Constraints on mantle anelasticity from geodetic observations, and implications for the J 2 anomaly, *Geophys. J. Int.*, **165**(1), 3–16.

* Bizouard, C., 2014. Le mouvement du pôle de l'heure au siècle, *Habilitation thesis*, Presses Académiques Francophones, 284 pp.

* Bizouard, C., Remus, F., Lambert, S.B., Seoane, L. & Gambis, D., 2011. The Earth's variable Chandler wobble, *Astron. Astrophys.*, **526**, A106, doi:10.1051/0004-6361/201015894.

* Boy, J.-P., Gégout, P. & Hinderer, J., 2002. Reduction of surface gravity data from global atmospheric pressure loading, *Geophys. J. Int.*, **149**, 534–545.

* Boy, J.-P. & Hinderer, J., 2006. Study of the seasonal gravity signal in superconducting gravimeter data, *J. Geodyn.*, **41**(1–3), 227–233.

* Boy, J.-P., Hinderer, J., Amalvict, M. & Calais, É., 2000. On the use of long records of superconducting and absolute gravity observations with special application to the Strasbourg station, France, *Cah. Cent. Eur. Géod. Séismol.*, **17**, 67–83.

* Boy, J.-P., Longuevergne, L., Boudin, F., Jacob, T., Lyard, F., Llubes, M., Florsch, N. & Esnoult, M.-F., 2009. Modelling atmospheric and induced non-tidal oceanic loading contributions to surface gravity and tilt measurements, *J. Geodyn.*, **48**, 182–188.

* Brzeziński, A. & Capitaine, N., 1993. The use of the precise observations of the celestial ephemeris pole in the analysis of geophysical excitation of Earth rotation, *J. geophys. Res.*, **98**(B4), 6667–6675.

* Brzeziński, A. & Nastula, J., 2002. Oceanic excitation of the Chandler Wobble, *Adv. Space Res.*, **30**(2), 195–200.

* Chao, B.F. & Chung, W.-Y., 2012. Amplitude and phase variations of Earth's Chandler wobble under continual excitation, *J. Geodyn.*, **62**, 35–39.

* Chen, X., Ducarme, B., Sun, H. & Xu, J., 2008. Loading effect of a self-consistent equilibrium ocean pole tide on the gravimetric parameters of the gravity pole tides at superconducting gravimeter stations, *J. Geodyn.*, **45**(4–5), 201–207.

* Chen, X., Kroner, C., Sun, H., Abe, M., Zhou, J., Yan, H. & Wziontek, H., 2009. Determination of gravimetric parameters of the gravity pole tide using observations recorded with superconducting gravimeters, *J. Geodyn.*, **48**(3–5), 348–353.

* Crossley, D. *et al.*, 1999. Network of superconducting gravimeters benefits a number of disciplines, *EOS, Trans. Am. geophys. Un.*, **80**(11), 121–132.

* Cummins, P., Wahr, J.M., Agnew, D.C. & Tamura, Y., 1991. Constraining core undertones using stacked IDA gravity records, *Geophys. J. Int.*, **106**(1), 189–198.
Desai, S.D., 2002. Observing the pole tide with satellite altimetry, *J. geophys. Res.*, **107**(April), 1–13.

Dickman, S.R., 1979. Consequences of an enhanced pole tide, *J. geophys. Res.*, **84**(B10), 5447–5456.

Dickman, S.R., 1988. The self-consistent dynamic pole tide in non-global oceans, *Geophys. J. R. astr. Soc.*, **94**(3), 519–543.

Dickman, S.R., 2005. Rotationally consistent Love numbers, *Geophys. J. Int.*, **161**(1), 31–40.

Ducarme, B., Venedikov, A.P., Arnoso, J., Chen, X.D., Sun, H. & Vieira, R., 2006. Global analysis of the GGP superconducting gravimeters network for the estimation of the pole tide gravimetric amplitude factor, *J. Geodyn.*, **41**(1–3), 334–344.

Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary reference Earth model, *Phys. Earth planet. Inter.*, **25**, 297–356.

Florsch, N. & Hinderer, J., 2000. Bayesian estimation of the free core nutation parameters from the analysis of precise tidal gravity data, *Phys. Earth planet. Inter.*, **117**, 21–35.

Gibert, D., Holschneider, M. & Le Mouël, J.-L., 1998. Wavelet analysis of the Chandler wobble, *J. geophys. Res.*, **103**, 27 069–27 089.

Gibert, D. & Le Mouël, J.-L., 2008. Inversion of polar motion data: Chandler wobble, phase jumps, and geomagnetic jerks, *J. geophys. Res.*, **113**, 1–8.

Gross, R.S., 2000. The excitation of the Chandler wobble, *Geophys. Res. Lett.*, **27**(15), 2329–2332.

Hinderer, J., Crossley, D. & Warburton, R.J., 2007. Gravimetric methods – superconducting gravity meters, in *Treatise on Geophysics*, Vol. 3, pp. 66–115, Elsevier B.V.

Hu, X.G., Liu, L.T., Ducarme, B., Xu, H.J. & Sun, H., 2007. Estimation of the pole tide gravimetric factor at the chandler period through wavelet filtering, *Geophys. J. Int.*, **169**(3), 821–829.

Loyer, S., Hinderer, J. & Boy, J.-P., 1999. Determination of the gravimetric factor at the Chandler period from Earth orientation data and superconducting gravimetry observations, *Geophys. J. Int.*, **136**(1), 1–7.

Malkin, Z. & Miller, N., 2010. Chandler wobble: two more large phase jumps revealed, *Earth Planets Space*, **62**(12), 943–947.

McKay, M.D., Beckman, R.J. & Conover, W.J., 1979. A comparison of three methods for selecting values of input variables in

the analysis of output from a computer code, *Technometrics*, **21**(2), 239–245.

Mc Carthy, D.D. & Petit, G., 2003. IERS Conventions (2003). (IERS Technical Note; 32) Frankfurt am Main: Verlag des Bundesamts für Kartographie und Geodäsie. 127 pp.

Miller, S.P. & Wunsch, C., 1973. The pole tide, *Nature*, **246**(01), 98–102.

Newville, M., Stensitzki, T., Allen, D.B. & Ingargiola, A., 2014. LMFIT: non-linear least-square minimization and curve-fitting for python, doi:10.5281/zenodo.11813.

Smylie, D.E., Henderson, G.A. & Zuberi, M., 2015. Modern observations of the effect of earthquakes on the Chandler wobble, *J. Geodyn.*, **83**, 85–91.

Wahr, J.M., 1982. The effects of the atmosphere and oceans on the Earth's wobble—I. Theory, *Geophys. J. R. astr. Soc.*, **70**, 349–372.

Wahr, J.M., 1983. The effect of the atmosphere and oceans on the Earth's wobble and on the seasonal variations in the length of day—II. Results, *Geophys. J. R. astr. Soc.*, **74**, 451–487.

Wahr, J.M., 1985. Deformation induced by polar motion, *J. geophys. Res.*, **90**(B11), 9363–9368.

Xu, J.Q., Sun, H. & Yang, X.F., 2004. A study of gravity variations caused by polar motion using superconducting gravimeter data from the GGP network, *J. Geod.*, **78**, 201–209.

## SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this paper:

(http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggw007/-/DC1).

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.