# Constraints provided by ground gravity observations on geocentre motions

Y. Rogister, Anthony Mémin, S. Rosat, J. Hinderer, Marta Calvo

► To cite this version:

Y. Rogister, Anthony Mémin, S. Rosat, J. Hinderer, Marta Calvo. Constraints provided by ground gravity observations on geocentre motions. Geophysical Journal International, 2016, 10.1093/gji/ggw220. hal-01330819

HAL Id: hal-01330819

https://hal.science/hal-01330819

Submitted on 14 Oct 2021

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
<figure><img src="image_2.png" /></figure>

# Constraints provided by ground gravity observations on geocentre motions

Y. Rogister,<fnref n="1" /> A. Mémin,<fnref n="2" />,<fnref n="3" /> S. Rosat,<fnref n="4" /> J. Hinderer<fnref n="4" /> and M. Calvo<fnref n="4" />,<fnref n="5" />

Accepted 2016 June 7. Received 2016 June 6; in original form 2015 September 2

## SUMMARY

The geocentre motion is the motion of the centre of mass of the entire Earth, considered an isolated system, in a terrestrial system of reference. We first derive a formula relating the harmonic degree-1 Lagrangian variation of the gravity at a station to both the harmonic degree-1 vertical displacement of the station and the displacement of the whole Earth's centre of mass. The relationship is independent of the nature of the Earth deformation and is valid for any source of deformation. We impose no constraint on the system of reference, except that its origin must initially coincide with the centre of mass of the spherically symmetric Earth model. Next, we consider the geocentre motion caused by surface loading. In a system of reference whose origin is the centre of mass of the solid Earth, we obtain a specific relationship between the gravity variation at the surface, the geocentre displacement and the load Love number $h'_{1}$, which demands the Earth's structure and rheological behaviour be known. For various networks of real or fictitious stations, we invert synthetic signals of surface gravity variations caused by atmospheric loading to retrieve the degree-1 variation of gravity. We then select six well-distributed stations of the Global Geodynamics Project, which is a world network of superconducting gravimeters, to invert actual gravity data for the degree-1 variations and determine the geocentre displacement between the end of 2004 and the beginning of 2012, assuming it to be due to surface loading. We find annual and semi-annual displacements with amplitude 0.5–2.3 mm.

**Key words:** Time-series analysis; Reference systems; Time variable gravity.

## 1 INTRODUCTION

According to Hoffmann-Wellenhoff & Moritz (2005), Petit & Luzum (2010) or Wu et al. (2012), the geocentre is defined as the centre of mass CM of the system made up of the solid Earth and its fluid envelopes. Its motion in an inertial system of reference is determined by the forces exerted on the Earth by other celestial bodies. We will call geocentre motion the motion of CM with respect to a given terrestrial system of reference, which is therefore non-inertial.

The geocentre motion is associated with the long-wavelength deformation of the Earth due to surficial and internal processes. In a reference frame whose origin is either the centre of surface figure CF of the solid Earth surface or the centre of mass CE of the solid Earth, its magnitude reaches a few millimetres at the seasonal timescale and is mainly due to the redistribution of mass of the continental water, the oceans and the atmosphere (e.g. Dong et al.

1997; Blewitt et al. 2001; Feissel-Vernier et al. 2006; Rietbroek et al. 2012; Wu et al. 2012). The geocentre is also subject to a secular drift induced by both the present-day ice mass change in Greenland and Antarctica and the isostatic adjustment that followed the last Pleistocene deglaciation (Greff-Lefftz 2000; Métivier et al. 2010). Therefore, the accurate determination of the geocentre motion can provide significant insights into global change effects, which include mean sea level rise.

Moreover, the knowledge of the centre of mass position is essential to establish a Terrestrial Reference Frame (TRF). Indeed, the deformation of the Earth surface can be observed with geodetic observing systems such as Global Navigation Satellite Systems (GNSS), Very Long Baseline Interferometry, Satellite Laser Ranging (SLR) and Doppler Orbit determination and Radiopositioning Integrated by Satellite. It is described by the relative position of each point at the surface of the Earth in a Terrestrial Reference System (TRS). A TRF is a realization of a TRS obtained by providing a

<footnote n="1">¹École et Observatoire des Sciences de la Terre, Université de Strasbourg, France. E-mail: yves.rogister@unistra.fr</footnote>

<footnote n="2">²School of Physical Sciences, University of Tasmania, Hobart, Australia</footnote>

<footnote n="3">³Observatoire de la Côte d'Azur, Université Nice Sophia Antipolis, CNRS, IRD, Géoazur UMR 7329, France</footnote>

<footnote n="4">⁴Institut de Physique du Globe de Strasbourg, UMR 7516, Université de Strasbourg/EOST, CNRS, France</footnote>

<footnote n="5">⁵Observatorio Geofísico Central, Instituto Geográfico Nacional (IGN), Madrid, Spain</footnote>
set of globally distributed station coordinates (e.g. Altamimi et al. 2011). To accurately describe the motion of Earth’s surface, the TRF needs to be determined with a high level of accuracy (Plag et al. 2009).

A TRF is established by using positioning observations determined by one or several geodetic systems (e.g. Altamimi et al. 2007, 2011; Rülke et al. 2008) and corrected for known phenomena, such as the solid Earth tides, ocean tides and atmospheric loading (McCarthy & Petit 2004; Collilieux et al. 2010; Petit & Luzum 2010). The International Terrestrial Reference Frame (ITRF) is the TRF that results from combining all the geodetic techniques mentioned above. The ITRF origin is, however, currently determined with SLR data only (Altamimi et al. 2007, 2011). Because of the linear kinematic model employed to describe the time evolution of the coordinates of the ITRF stations, the ITRF origin is a long-term average of CM (Wu et al. 2012). As a result, the ITRF origin is not the instantaneous CM but moves with respect to CM. The methods used to estimate geocentre motion from GPS data are described by Dong et al. (2003). Wu et al. (2012) also described those methods, as well as inverse methods involving additional sets of data such as GRACE (Gravity Recovery And Climate Experiment) gravity data and data-assimilated ocean bottom pressure models (Rietbroek et al. 2012).

According to Altamimi et al. (2011), the ITRF origin accuracy achievable is at the level of or better than 1 cm. To improve that accuracy and, therefore, our understanding of global geodetic and geophysical phenomena, whose precise measurements depend on the position of the ITRF origin, we will develop a suggestion made by Plag et al. (2007) and investigate the constraints that ground gravity observations, which constitute a data set independent from the reference frame, can bring on the geocentre position. Two questions then arise: (1) what is the analytical relationship between geocentre displacement and surface gravity variation and (2) would the latter be detectable by superconducting gravimeters (SG), which are the most accurate relative gravimeters? The main aim of this study is to answer both questions.

The paper is organized as follows. In Section 2, we establish a general relationship, valid in any reference frame whose origin initially coincides with the entire Earth’s centre of mass, between the harmonic degree-1 gravity variation at a station, the displacement of the station and the geocentre motion. Next, in Section 3, we consider the usual treatment of surface loading with a Love number approach to derive a relationship between the gravity variation at the surface and the displacement of the whole Earth’s centre of mass in a system of reference whose origin is the centre of mass of the solid Earth. Section 4 is devoted to the computation of the power spectral density (PSD) of data acquired with an SG installed in Strasbourg, France. The PSD allows us to check that the annual gravity signal associated to the geocentre motion might be above the noise level. In Section 5, we invert ground gravity variation to retrieve its degree-1 components. The inversion of synthetic signals allows us to study the influence of the world distribution of the gravity stations on the quality of the retrieved degree-1 components. We then invert gravity data provided by six well-distributed stations of the Global Geodynamics Project (GGP; Hinderer & Crossley 2000), which is a world network of 38 SG stations, between 2004 and 2012. From the inverted degree-1 gravity variations, we derive time-series for the geocentre displacement by assuming that it is due to surface loading and estimate the annual and semi-annual displacements. We conclude in Section 6.

## 2 GRAVITY VARIATION DUE TO MOTION OF EARTH’S CENTRE OF MASS

We consider a spherical model of the whole Earth whose mass is $M$ and radius, $a$. The system of reference is chosen in a such a way that its origin initially coincides with the centre of mass CM. We denote by $r$ the radial coordinate of the spherical polar coordinates $(r, \theta, \varphi)$, $\theta$ being the colatitude and $\varphi$, the longitude. The gravity $g_0(r)$ at $r \ge a$ is

$$g_0(r) = -\frac{GM}{r^2}e_r, \qquad (1)$$

where $e_r$ is the radial unit vector pointing outward. $g_0(r)$ is the gradient of the gravitational potential $V_0(r)$

$$g_0(r) = -\nabla V_0(r), \qquad (2)$$

where

$$V_0(r) = -\frac{GM}{r}. \qquad (3)$$

Let us assume that the initial spherical configuration is slightly perturbed. A particle initially at $r$ is therefore displaced to $r + s$, where $s$ is the displacement vector in a given system of reference. We impose no specific requirements to the system of reference other than its origin must initially coincide with the initial position of Earth’s centre of mass CM. The change in gravity the particle experiences is the Lagrangian variation of gravity, which we shall denote by $g^\delta$. The gravity variation at $r$ is the Eulerian variation $g^\Delta$ that is related to the Eulerian variation of the gravitational potential $V^\Delta$ by

$$g^\Delta = -\nabla V^\Delta. \qquad (4)$$

If there are no other attracting masses than the Earth itself, a gravimeter moving with the ground would measure $g^\delta$ while an instrument fixed at $r$ would measure the gravity variation $g^\Delta$.

In a first-order approximation, the infinitesimal Lagrangian perturbation $g^\delta$ is related to $V^\Delta$ by (Dahlen & Tromp 1998, p. 64)

$$g^\delta = -\nabla V^\Delta + s \cdot \nabla g_0. \qquad (5)$$

For a spherical reference configuration, the radial component of this equation is

$$g^\delta = -\frac{\partial V^\Delta}{\partial r} - U \frac{dg_0}{dr}, \qquad (6)$$

where $U$ is the radial displacement. At $r = a$, this becomes

$$g^\delta(a) = -\left.\frac{\partial V^\Delta}{\partial r}\right|_{r=a} + 2\frac{U(a)}{a}g_0(a). \qquad (7)$$

When the total external gravitational potential $V = V_0 + V^\Delta$ is decomposed in a series of spherical harmonics, its degree-1 terms are directly related to the position $(x^{\text{CM}}, y^{\text{CM}}, z^{\text{CM}})$ of CM (Hoffmann-Wellenhoff & Moritz 2005, pp. 56–64):

$$V = {}_1^cV \sin\theta \cos\varphi + {}_1^0V \cos\theta + {}_1^sV \sin\theta \sin\varphi, \qquad (8)$$

where

$$({}_1^cV, {}_1^sV, {}_1^0V) = -g_0(x^{\text{CM}}, y^{\text{CM}}, z^{\text{CM}}) \qquad (9)$$

The negative sign comes from our definition of the gravitational potential which is the opposite of the definition by Hoffmann-Wellenhoff & Moritz (2005).
The three degree-1 coefficients $c_1^i V^\Delta$, $c_1^0 V^\Delta$ and $c_1^i V^\Delta$ of the Eulerian variation of the potential are therefore proportional to the components $s_x^{CM}$, $s_y^{CM}$ and $s_z^{CM}$ of the displacement of CM with respect to the system of reference:

$$ (c_1^i V^\Delta, c_1^0 V^\Delta, c_1^0 V^\Delta) = -g_0 (s_x^{CM}, s_y^{CM}, s_z^{CM}) \quad (10) $$

so that

$$ V^\Delta = -g_0 s_r^{CM}, \quad (11) $$

where $s_r^{CM}$ is the CM displacement in the direction of $\mathbf{r}$. From now on, $U(a)$ and $g_0(a)$ will be simply denoted by $U$ and $g_0$, respectively. The Lagrangian radial gravity variations of harmonic degree 1 at the surface are then given by

$$ (c_1^i g^\delta, c_1^0 g^\delta, c_1^0 g^\delta) = \frac{2}{a} g_0 (c_1^i U - s_x^{CM}, c_1^i U - s_y^{CM}, c_1^i U - s_z^{CM}) \quad (12) $$

where $c_1^i U$, $c_1^0 U$ and $c_1^i U$ are of course the degree-1 components of the radial surface displacement $U$. Thus, $g^\delta$ is

$$ g^\delta = \frac{2}{a} g_0 (U - s_r^{CM}). \quad (13) $$

We can quickly show from this formula that $g^\delta$ is independent of the system of reference, provided its origin initially coincides with the centre of mass of the spherically symmetric Earth model. Indeed, if we consider two systems of reference with two different origins in motion with respect to each other, the expressions of both the relative displacements $\mathbf{s}$ and $\mathbf{s}^{CM}$ in the two systems of reference will differ by a quantity which is the displacement of one origin with respect to the other. The difference between $U$ and $s_r^{CM}$ in the right-hand side of eq. (13) therefore makes the relative displacement of the two origins cancel out.

$g^\delta$ being measured by a gravimeter moving with the ground and $U$ being obtained by positioning measurements, the formula (12) provides a means to compute the displacement of the centre of mass with respect to the system of reference. For a difference $U - s_r^{CM}$ equal to 1 mm, the gravity variation given by eq. (13) is approximately $3 \times 10^{-10} g_0$ or $300 \text{ nGal} = 3 \text{ nm s}^{-2}$.

It should be noted that eq. (12) is valid for any kind of degree-1 deformation, whatever its cause and the Earth's rheological behaviour. In the next section, we derive a specific formula for the geocentre motion caused by surface loading.

## 3 SURFACE LOADING

Most often, studies of the motion of the centre of mass have their origin in surface loading problems. Love numbers are then frequently used to relate the gravity variations, surface displacement and centre of mass displacement to the load (e.g. Blewitt 2003). In this section, we connect this Love number approach to the formula (12) established in the previous section.

If we denote by $c_1^i \sigma$, $c_1^0 \sigma$ and $c_1^i \sigma$ the three degree-1 components of the surface mass load $\sigma$, similarly to the notation for the components of the Earth gravitational potential in eq. (8), the degree-1 gravitational acceleration due to the load $m_1^i g^{\text{load}}$ derives from the potential

$$ m_1^i \phi(r) = - \frac{4\pi G a^3}{3r^2} m_1^i \sigma, \quad (14) $$

where $m$ stands for $c$, 0 or $s$. The degree-1 load Love numbers $h'_i$ and $k'_i$ for the radial displacement and perturbation of the gravitational potential are defined, respectively, by

$$ m_1^i U = - \frac{m_i^i \phi(a)}{g_0} h'_i \quad (15) $$

and

$$ m_1^i V^\Delta = (1 + k'_i)^m \phi(a). \quad (16) $$

Relation (16) expresses that the total potential variation is the sum of the gravitational potential of the load and the perturbation of the Earth's self-gravitational potential, which is proportional to the former. Because of the spherical symmetry of the initial configuration, $h'_i$ and $k'_i$ do not depend on $m$. It is obvious from eqs (10) and (16) that $k'_i^{\text{CM}} = -1$ in a system of reference whose origin is CM. In such a system, $h'_i^{\text{CM}} = -1.288$ for the PREM elastic Earth model (Dziewonski & Anderson 1981). In a CE frame, $k'_i^{\text{CE}} = 0$ and $h'_i^{\text{CE}} = -0.288$. Further discussion of the load Love numbers in different reference frames can be found in Blewitt (2003).

Before we launch into the calculation of $m_1^i g^\delta$, we remind the reader how the analogue tidal problem is usually solved. The tidal gravimetric factor $\delta_\ell$ is a dimensionless number defined, for a given harmonic degree $\ell$, by the ratio between the Lagrangian gravity variation measured at the surface and the tidal attraction $m_\ell^i g^{\text{tide}}$ (Hinderer et al. 1991):

$$ m_\ell^i g^\delta = m_\ell^i g^{\text{tide}} \delta_\ell, \quad (17) $$

where $m_\ell$ is the harmonic order. Since we assume that the Earth deformation is linearly proportional to the perturbing potential, the gravimetric factor, like the tidal Love numbers $h_\ell$ and $k_\ell$ defined similarly to $h'_\ell$ and $k'_\ell$ with the load potential replaced by the tidal potential, depends on the Earth model and its rheological behaviour but is independent of the perturbation itself. The gravity variation at the surface is made up of three terms: the direct attraction of the perturbing body $g^{\text{tide}}$, the gravity variation due to the displacement of the gravimeter in the surrounding gravity field and the gravity variation due to mass redistribution inside the Earth. $\delta_\ell$ is therefore given by

$$ \delta_\ell = 1 + \frac{2}{\ell} h_\ell - \frac{\ell + 1}{\ell} k_\ell. \quad (18) $$

In the case of a loading problem, Longman (1963) showed that, right under the load, a relation similar to eq. (17) also holds:

$$ m_\ell^i g^\delta = m_\ell^i g^{\text{load}} \left( 1 + \frac{2}{\ell} h'_\ell - \frac{\ell + 1}{\ell} k'_\ell \right). \quad (19) $$

Right above the load, the relation between $m_\ell^i g^\delta$ and $m_\ell^i g^{\text{load}}$ is slightly different because $m_\ell^i \phi$ is proportional to $r^\ell$ inside the Earth and proportional to $r^{-\ell-1}$ outside. It is

$$ m_\ell^i g^\delta = m_\ell^i g^{\text{load}} \left( 1 - \frac{2}{\ell + 1} h'_\ell + k'_\ell \right), \quad (20) $$

which gives

$$ m_1^i g^\delta = m_1^i g^{\text{load}} (1 - h'_1 + k'_1), \quad (21) $$

for $\ell = 1$.

To establish a correspondence between eqs (12) and (21), we must relate $m_1^i g^{\text{load}}$ to the displacement of the centre of mass CM. That we will do in the CE frame. First, we note that the position vector $\mathbf{r}^{\text{load}}$ of the centre of mass of the load has components

$$ \mathbf{r}^{\text{load}} = (x^{\text{load}}, y^{\text{load}}, z^{\text{load}}) = \frac{4\pi a^3}{3 M^{\text{load}}} (c_1^i \sigma, c_1^0 \sigma, c_1^i \sigma), \quad (22) $$
where $M^{\text{load}}$ is the mass of the load. Consequently, above the load, the harmonic components of its gravitational potential (14) are

$$ (\dot{\phi}_1, \dot{\phi}_1, \dot{\phi}_0, \dot{\phi}_0) = -\frac{GM^{\text{load}}}{r^2}(x^{\text{load}}, y^{\text{load}}, z^{\text{load}}). \qquad (23) $$

Second, in the system of reference whose origin is the centre of mass CM of the system made up of the solid Earth and the load, $\mathbf{r}^{\text{load}}$ and the position $\mathbf{r}^{\text{CE}}$ of the centre of mass CE of the solid Earth obey the relation

$$ M^{\text{load}}\mathbf{r}^{\text{load}} + M^{\text{CE}}\mathbf{r}^{\text{CE}} = 0, \qquad (24) $$

where $M^{\text{CE}}$ is the mass of the solid Earth. Eq. (24) is the very definition of CM. If, initially, CM and CE coincide, the displacement of CE in the CM frame is $\mathbf{r}^{\text{CE}}$ and is the opposite of the displacement $\mathbf{s}^{\text{CM}}$ of CM in the CE frame: $\mathbf{r}^{\text{CE}} = -\mathbf{s}^{\text{CM}}$. Since

$$ \dot{\mathbf{g}}_1^{\text{load}} = -\frac{\partial \dot{\mathbf{\phi}}_1}{\partial r} \qquad (25) $$

and $\dot{\mathbf{k}}_1^{\text{CE}} = 0$, we therefore have at the surface

$$ (\dot{\mathbf{c}}_1^{\text{CE}}\mathbf{g}^{\delta}, \dot{\mathbf{s}}_1^{\text{CE}}\mathbf{g}^{\delta}, \dot{\mathbf{0}}_1^{\text{CE}}\mathbf{g}^{\delta}) = \frac{2}{a}\mathbf{g}_0(\dot{\mathbf{h}}_1^{\text{CE}} - 1)(\mathbf{s}_x^{\text{CM}}, \mathbf{s}_y^{\text{CM}}, \mathbf{s}_z^{\text{CM}}) \qquad (26) $$

in the CE frame. For the elastic PREM model, we have approximately

$$ \mathbf{g}^{\delta} = 3.966 \cdot 10^{-6} \mathbf{s}^{\text{CM}} \mathbf{m} \mathbf{s}^{-2} \qquad (27) $$

if $\mathbf{s}^{\text{CM}}$ is expressed in metre. A 1 mm displacement of CM with respect to CE therefore gives a gravity variation of approximately 400 nGal = 4 nm s⁻².

As a consistency check of our results, combining eq. (12) with eq. (26) yields

$$ (\dot{\mathbf{c}}_1^{\text{CE}}\mathbf{U}, \dot{\mathbf{s}}_1^{\text{CE}}\mathbf{U}, \dot{\mathbf{0}}_1^{\text{CE}}\mathbf{U}) = \dot{\mathbf{h}}_1^{\text{CE}}(\mathbf{s}_x^{\text{CM}}, \mathbf{s}_y^{\text{CM}}, \mathbf{s}_z^{\text{CM}}) \qquad (28) $$

in the CE frame, in agreement with eqs (13) and (14) of Blewitt (2003).

It is noteworthy that neither eq. (12) nor eq. (26) depend on the tangential displacement or Love number $\dot{l}'^{\text{CE}}$, which is defined by a relation similar to eq. (15) with $\dot{h}'_1^{\text{CE}}$ replaced by $\dot{l}'^{\text{CE}}$ and the radial displacement $\dot{\mathbf{u}}_1^{\text{CE}}$ replaced by the tangential scalar of the spheroidal displacement field.

In the next two sections, we explore how ground gravity data provided by the GGP network of 38 SG stations (Fig. 2), eqs (12) or (26) might offer valuable information on the geocentre motion complementary to the methods based on satellites-only observations.

## 4 POWER SPECTRAL DENSITY OF GRAVITY DATA FROM SUPERCONDUCTING METERS

Geophysical data and models of mass redistribution in the atmosphere, ocean, surface water and glaciers (Dong et al. 1997) or satellite observations of the displacement of geodetic stations (Blewitt et al. 2001; Wu et al. 2012) have revealed annual and seasonal geocentre and stations displacements with amplitude of order of a few millimetres either in a CE or CF frame. As mentioned in Sections 2 and 3, the gravity variation given by eq. (13) is approximately 300 nGal = 3 nm s⁻² for a difference $\mathbf{U} - \mathbf{s}_r^{\text{CM}}$ equal to 1 mm and the gravity variation given by eq. (26) is approximately 400 nGal = 4 nm s⁻² for a CM displacement of 1 mm caused by surface loading on an elastic Earth model. According to fig. 2 of van Camp et al. (2005), an annual gravity signal of a few microgals is possibly detectable with a single SG. The signal would contain

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Residual power spectral density of a 18-yr long series (1996–2014) of gravity data recorded with the superconducting gravimeter C026 at the Strasbourg station. 160 656 samples at a sampling rate of 1 hr were used. At a period of one year, the noise level is approximately −83 dB. Gravity variations larger than 620 nGal = 6.2 nm s⁻² would therefore be above the noise level. Also shown is the 300 nGal = 3 nm s⁻² level equivalent to a degree-1 gravity variation generated by a difference of 1 mm between the displacement of the Earth’s centre of mass CM and the radial displacement of the station (eq. 13).</figcaption></figure>

the total gravity variation, not only the degree-1 variation. That detectability threshold is largely confirmed by the residual PSD, shown in Fig. 1, of the total gravity variation recorded at the Strasbourg station between 1996 and 2014.

Before the PSD of Fig. 1 has been computed, the raw data had first been corrected for the oceanic and solid Earth tides by using the models NAO99 of Matsumoto et al. (2000) and DDW99 of Dehant et al. (1999), which include the main semi-annual, monthly, fortnightly, ter-monthly, diurnal, semi-diurnal and higher frequency tidal waves. For the annual wave, we use eq. (17) for an elastic Earth model with a real gravimetric factor $\delta_2$ equal to 1.16. The effects of the local atmospheric pressure were also removed, the local barometric admittance being −3 nm s⁻² hPa⁻¹, as well as the gravity effect of polar motion, the pole coordinates being publicly available on the International Earth Rotation Service website. The residual PSD has then been estimated using a modified Welch periodogram. The PSD of a Gaussian white noise is $\sigma_G^2 T_0$, where $\sigma_G$ is the standard deviation and $T_0$ is the sampling rate, and the PSD of an undamped harmonic signal of amplitude A is $A^2 NT_0/4$, where N is the number of samples. Consequently, the signal-to-noise ratio is $A\sigma_G\sqrt{N}/2$. The PSD of Fig. 1 has been computed with segments of N/4 points, an overlapping of N/8 points and the two parameters $T_0 = 3600$ s and N = 160 984, giving a total duration of 6694 d. It shows that at the annual period, the noise level is about −82 dB, which implies that an annual gravity variation larger than 620 nGal = 6.2 nm s⁻² can be detected. If it was a degree-1 variation, it would be equivalent to a difference $\mathbf{U} - \mathbf{s}_r^{\text{CM}}$ of 2 mm.

In the following, we will invert ground gravity data to retrieve the three degree-1 components $(\dot{\mathbf{c}}_1^{\text{CE}}\mathbf{g}^{\delta}, \dot{\mathbf{s}}_1^{\text{CE}}\mathbf{g}^{\delta}, \dot{\mathbf{0}}_1^{\text{CE}}\mathbf{g}^{\delta})$ of the gravity variation. We will use a stacking method that was initially proposed by Courtier et al. (2000) to search for the degree-1 Slichter modes, which are the three translation modes of the Earth’s inner core. It consists in inverting time variations of gravity recorded at the surface of the Earth, the data being weighted by the degree-1
spherical harmonics evaluated at the station coordinates. This stacking method has been successfully applied to the first detection of the seismic mode $_2S_1$ by using data from the GGP network (Rosat et al. 2003). Prior to the estimate of the geocentre motion from data of the GGP network, which is sparser than the global GPS network, we first need to perform tests with synthetic gravity signals to eliminate as much as possible any bias due to the network configuration.

# 5 INVERSION OF GROUND GRAVITY DATA

## 5.1 Tests with synthetic gravity signals

The tests we perform consist in inverting synthetic gravity variations for various sets of gravity stations, either real or fictitious, the number of stations ranging down from 6 up to 180. The total synthetic gravity variation at a station is obtained by computing the Earth's elastic response to the atmospheric loading quantified by pressure data provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) spanning the years 2002–2012.

To do so, we first decompose the ECMWF surface pressure fields into a sum of spherical harmonics with components $m_l^l P$. The system of reference is such that the z-axis is the pole axis and the x- and y-axes are in the equatorial plane, with x directed in the Greenwich meridian plane and y to the east. The components of the surface mass load are related to $m_l^l P$ by

$$m_l^l \sigma = \frac{m_l^l P}{g_0}. \qquad (29)$$

This surface mass distribution generates, above the load, the gravitational potential

$$m_l^l \phi(r) = - \frac{4\pi G a}{2l+1} \left(\frac{a}{r}\right)^{l+1} m_l^l \sigma. \qquad (30)$$

Eq. (14) is a specific case of eq. (30) for $l = 1$. $m_l^l g^{\text{load}}$ then derives from $m_l^l \phi(r)$.

Second, using eq. (20) with load Love numbers for PREM (Dziewonski & Anderson 1981), we obtain the gravity variation $m_l^l g^\delta$, in particular $m_l^l g^\delta$ given by eq. (21). We assume that the oceans respond as an inverted barometer to the atmospheric pressure variations. Indeed, the main periodic component of the geocentre motion is annual, which is long enough for the oceans to compensate for the atmospheric pressure changes by adjusting their surface level according to the inverted barometer hypothesis.

Third, the synthetic signals are the sum of $m_l^l g^\delta$ and Gaussian white noise with a standard deviation of $2 \text{ nm s}^{-2}$. This arbitrarily chosen but meaningful value can be compared to the standard deviation of $m_l^l g^\delta$ which ranges from 4 to $8 \text{ nm s}^{-2}$, depending on the location.

Fourth and final, we invert the synthetic signals by using the stacking method of Courtier et al. (2000) to obtain the degree-1 gravity variations $m_l^l g^{\text{inv}}$ that we compare to the degree-1 synthetic gravity variations $m_l^l g^\delta$ with the help of the relative error of the temporal residuals defined in the frequency domain by

$$m^l e = \sqrt{\frac{\sum (m_l^l g^{\text{inv}} - m_l^l g^\delta)^2}{\sum (m_l^l g^\delta)^2}}, \qquad (31)$$

the sums extending over gravity values computed at frequencies included in the range [$10^{-8}$, $10^{-7}$] Hz, which brackets the annual frequency. Table 1 sums up the relative errors estimated for eight networks shown in Figs 2 and 3:

<table><caption>**Table 1.** Relative error $m^l e$ (in per cent) defined by eq. (31) in the annual band [$10^{-8}$, $10^{-7}$] Hz for the inversion of synthetic data which are the Earth's elastic response to atmospheric loading provided by actual pressure data from the ECMWF between 2002 and 2012. The eight networks, which consist in sets of locations where the synthetic data are computed, are described in the text. Some locations correspond to stations of the GGP network, others correspond to fictitious stations. NY and SY are GGP stations at Ny-Ålesund, Svalbard and Syowa, Antarctica, respectively.</caption><thead><tr><th></th><th>Network</th><th><sup>c</sup><sub>e</sub></th><th><sup>s</sup><sub>e</sub></th><th><sup>0</sup><sub>e</sub></th></tr></thead><tbody><tr><td>1</td><td>Locations of 38 GGP stations</td><td>7.4</td><td>6.2</td><td>5.9</td></tr><tr><td>2</td><td>Locations of 6 GGP stations</td><td>6.3</td><td>6.4</td><td>6.3</td></tr><tr><td>3</td><td>Network 2 + Tahiti + Hawaii</td><td>6.4</td><td>6.3</td><td>6.3</td></tr><tr><td>4</td><td>Network 2 + NY + SY</td><td>7.3</td><td>5.9</td><td>6.2</td></tr><tr><td>5</td><td>Network 3 + NY + SY</td><td>7.2</td><td>6.0</td><td>5.9</td></tr><tr><td>6</td><td>180 gridpoints</td><td>6.9</td><td>6.0</td><td>6.1</td></tr><tr><td>7</td><td>62 continental points</td><td>7.4</td><td>6.5</td><td>6.2</td></tr><tr><td>8</td><td>118 oceanic points</td><td>5.5</td><td>6.1</td><td>5.8</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** GGP network, which includes 38 superconducting gravimeters (white diamonds, blue rectangles and green stars), and two fictitious stations in Tahiti and Hawaii. Network 2 of Table 1 is made of the six stations pictured by blue rectangles (CA = Cantley (Canada), BH = Bad Homburg (Germany), KA = Kamioka (Japan), TC = Tigo-Concepción (Chile), SU = Sutherland (South Africa), CB = Canberra (Australia)). The two green stars are polar stations in Ny-Ålesund (NY), Svalbard and Syowa (SY), Antarctica.</figcaption></figure>

(i) Locations of the 38 GGP stations pictured by white diamonds, blue rectangles and green stars in Fig. 2.

(ii) Locations of six GGP stations with a homogeneous world distribution, pictured by blue rectangles in Fig. 2.

(iii) Same six locations as in Network 1 plus two fictitious stations in Tahiti and Hawaii pictured by red dots in Fig. 2.

(iv) Same six locations as in Network 1 plus two locations corresponding to two polar GGP stations at Ny-Ålesund (NY), Svalbard and Syowa (SY), Antarctica, pictured by green stars in Fig. 2.

(v) Same eight locations as in Network 3 plus the two locations NY and SY.

(vi) Set of 180 fictitious stations regularly distributed at the surface of the Earth (Fig. 3).

(vii) Subset of Network 6 made of 62 continental stations pictured by black dots in Fig. 3.

(viii) Subset of Network 6 made of 118 oceanic stations pictured by blue diamonds in Fig. 3.

What could be called the minimum error of the inversion procedure is obtained by inverting, for any network, synthetic gravity data not containing additional white noise. The relative errors would be zero if the inversion procedure was exact but we find
<figure><img src="image_5.png" /><figcaption>**Figure 3.** 180 gridpoints forming Network 6 in Table 1. The 62 black dots are fictitious land stations and the 118 blue diamonds are fictitious ocean stations.</figcaption></figure>

$^m e = 6.1$ per cent, which is also the errors obtained with the synthetic data of Network 6, which include white noise. The reason that the error for this network, which is the densest network we consider, is unchanged with additional noise is that in the annual band the added noise is much smaller than the signal and does not deteriorate the results of the inversion when the density of the gravity stations over the world is high enough.

When discussing the results of Table 1, one should keep in mind
the relation (26) between the degree-1 gravity variations and the
geocentre displacement when a load deforms the surface of the
Earth. It shows how, in that case, the uncertainty in the estimate of
degree-1 gravity variations directly translates into uncertainty in the
geocentre displacement. The errors on the *c*, *s* and zero components
of the degree-1 gravity variations are actually errors on the *x*-, *y*-
and *z*-components of the geocentre displacement, respectively.

Although the number of stations varies a lot, from 6 to 180, the
errors are roughly the same for the eight networks, from 5.9 to
7.8 per cent. Of course, the best network, which performs equally
well for the three components of the geocentre displacement, is

Network 6, which is the densest network. The variations of the
relative errors between the networks can be explained by the ge-
ographical distributions of the stations. For instance, adding two
polar stations to either Network 2 or 3 diminishes a little the error
on the polar component of the geocentre displacement. Another
example is the better determination of the *x*-component of the geo-
centre displacement with a network of oceanic stations than with
a network of continental stations because there are more oceanic
stations close to both the Greenwich meridian in the Atlantic Ocean
and the antipodes in the Pacific Ocean.

Gravity variations of degree greater than 1 may induce regionally coherent signals that may contaminate estimate of degree-1 variations. It results that the degree-1 gravity variations may be difficult to isolate using our optimal Network 2, involving only six stations. To assess the effect of higher harmonic degrees on the estimate of the degree-1 gravity variations we perform a complementary test. The largest spatial variability of changes in atmospheric surface pressure is distributed between degrees 1 and 7, corresponding to gravity variations of a wavelength down to 2860 km, while that due to hydrology ranges between degrees 1 and 16, for wavelengths down to 1250 km (de Linage et al. 2009). Given that, we use for this test two years of the hydrology component of the updated European Spatial Agency (ESA) Earth System Model (ESM; Dobslaw et al. 2015), from 2005 January until 2006 December. The ESA ESM is provided as spherical harmonic coefficients or Stokes coefficients up to degree and order 180, which correspond to a wavelength of approximate 111 km. We first use the degree-1 Stokes coefficients to compute the degree-1 gravity variations $^m_1 g^\delta$ that result from the hydrological loading assuming the Earth's response is elastic. Then, we reconstruct hydrological fields by successively adding to the degree 1 degrees 2 to 180, denoting by $\ell_{max}$ the largest degree term of a specific field. We use Network 2 to invert for the degree-1 gravity variations induced by the hydrological loading from the reconstructed fields and we assess our inversions by computing the errors (31), plotted as a function of $\ell_{max}$ in Fig. 4. The errors are at the 6 per cent level when only considering degree 1 in the reconstructed fields. They are consistent with the values in Table 1.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Errors, given by eq. (31), for degree-1 gravity variations obtained at the six stations of Network 2 of Table 1 by inverting synthetic hydrological signals containing spherical harmonics terms up to degree $\ell_{max}$. For the three components of the degree-1 gravity variations, the errors are about 6 per cent for $\ell_{max} = 1$, as in Table 1. They fluctuate with increasing $\ell_{max}$ and reach almost constant values when $\ell_{max}$ is approximately 40.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) Spectra of the degree-1 gravity variations inverted from the gravity data acquired from late 2004 until early 2012 at the six GGP stations pictured by blue rectangles in Fig. 2 and forming Network 2 in Table 1. (b) Black curves: components of the geocentre displacement given by eq. (26), which relates the degree-1 gravity variations to the geocentre displacement supposedly caused by surface loading. Red curves: fits of the black curves, after removal of the trends, with semi-annual and annual sine functions whose amplitudes and phases are given in Table 2.</figcaption></figure>

As expected, gradually considering degrees up to 20 leads to larger errors, reaching up to almost 60 per cent for $^0e$. Adding signals of higher degrees leads to a convergence that is reached for degree 40. Errors converge to 18, 36 and 20 per cent for the *c*, *s* and zero components of the degree-1 gravity variations, respectively. This decrease of the errors by extending the hydrological fields to higher degrees shows that any regional or local signal, which is not globally coherent, will cancel out by stacking and hence helps in retrieving the gravity changes coherent to all stations to all stations, that is degree-1 variations.

To invert real data from the GGP network in the next section, we select Network 2, which at the same time contains fewer stations

and provides almost equal errors for the three components of the geocentre displacement because the geographical distribution of the stations is almost symmetrical with respect to both the equator and the Greenwich meridian.

## 5.2 Inversion of actual SG data

We thus pick the six stations forming Network 2 in Table 1. They are Bad-Homburg (Germany), Cantley (Canada), Canberra (Australia), Kamioka (Japan), Sutherland (South Africa) and Tigo-Concepción (Chile). Any European station other than Bad-Homburg would also
<table><caption>**Table 2.** Annual and semi-annual geocentre displacement. $s^{\text{CM}}$ is estimated in the CE frame, is assumedly caused by surface loading and is obtained from eq. (26) and measurements of surface gravity variations with Network 2 of Table 1 between 2004 and 2012. $\text{SLR}s^{\text{CM}}$ is estimated for the same period of time from geocentre motion time-series derived from SLR observations (Cheng et al. 2010). For each component, the fit is $A_a \cos 2\pi(t - t_0 + \gamma_a)/T + A_{sa} \cos 4\pi(t - t_0 + \gamma_{sa})/T$, where $T = 1$ yr, a and sa stand for annual and semi-annual, respectively, and $t_0$ is 2004 January 1.</caption><thead><tr><th rowspan="2">Component of CM displacement</th><th colspan="2">Annual</th><th colspan="2">Semi-annual</th></tr><tr><th>A<sub>a</sub> (mm)</th><th>γ<sub>a</sub> (day)</th><th>A<sub>sa</sub> (mm)</th><th>γ<sub>sa</sub> (day)</th></tr></thead><tbody><tr><td>γ<sub>x</sub><sup>CM</sup></td><td>0.5</td><td>106</td><td>1.8</td><td>171</td></tr><tr><td>γ<sub>y</sub><sup>CM</sup></td><td>1.2</td><td>271</td><td>1.5</td><td>176</td></tr><tr><td>γ<sub>z</sub><sup>CM</sup></td><td>1.6</td><td>88</td><td>2.3</td><td>1.7</td></tr><tr><td>SLR<sub>x</sub><sup>CM</sup></td><td>2.1 ± 0.5</td><td>315 ± 14</td><td>0.7 ± 0.3</td><td>93 ± 40</td></tr><tr><td>SLR<sub>y</sub><sup>CM</sup></td><td>2.6 ± 0.4</td><td>34 ± 11</td><td>0.2 ± 0.3</td><td>71 ± 68</td></tr><tr><td>SLR<sub>z</sub><sup>CM</sup></td><td>4.1 ± 1.3</td><td>4 ± 36</td><td>2.0 ± 1.3</td><td>186 ± 72</td></tr></tbody></table>

suit but Bad-Homburg is a good site in terms of noise levels (Rosat & Hinderer 2011) and length of time records. We correct the raw data from the end of 2004 until the beginning of 2012 the same way we corrected the Strasbourg data to compute the PSD in Section 4. Moreover, to remove the hydrological effects of high harmonic degrees and lessen their influence on the inversion for the degree-1 gravity variations, we correct the SG data from the GRACE time-series from degrees 2 to 50, provided by the CSR RL05 solutions.

The multistation inversion method provides the degree-1 gravity variations $\dot{m}_1^\delta g^\delta$ whose spectra are plotted in Fig. 5(a). If the geocentre displacement is entirely due to surface loading, it can be obtained in the CE frame from the three components $\dot{m}_1^\delta g^\delta$ by using eq. (26). It should be noted that, before invoking eq. (26), we have not assumed that the degree-1 gravity variations were caused by surface loading. In Fig. 5(b) are plotted the x-, y- and z-components of the geocentre displacement $s^{\text{CM}}$ as a function of time. High-frequency variations have been removed by sampling the gravity data at 15-d intervals, instead of the 1-hr sampling for the PSD in Fig. 1 and the spectra in Fig. 5(a). After removal of the trends from these time-series, they are fitted with both an annual and a semi-annual sine functions, whose amplitudes and phases are given in Table 2. The total fits are plotted in red in Fig. 5(b).

Table 2 of Wu et al. (2012) contains a collection of estimates of the annual geocentre displacement published between 1997 and 2010. Various methods based on various geodetic techniques, sometimes combined with GRACE observations or data-assimilated ocean bottom pressure model, spanning different time intervals included in the interval 1992–2010.6, provide a large scatter of amplitudes ranging between 0.0 and 4.2 mm, 1.6 and 4.9 mm and 1.5 and 7.6 mm for the x-, y- and z-components, respectively. Estimates based on GNSS are, however, dubious (Rebischung et al. 2014).

Our results given in Table 2 cannot be exactly compared to those measurements. The first reason is that the displacements reported by Wu et al. (2012) are measured with respect to CF whereas we used a CE frame. However, this is only a minor drawback since CF is a good approximation to CE (Dong et al. 1997; Blewitt 2003). The second and major reason is that the time span (2004–2012) of the ground gravity recordings we stacked does not exactly coincide with any of the time spans of the studies cited by Wu et al. (2012). Nevertheless, given the large scatter of the annual geocentre displacements provided by geodetic techniques, they are in a satisfactorily good agreement with our estimates using gravity data, especially if the errors between 18 and 36 per cent associated to the inversion of the gravity data are taken into account.

A better comparison can, perhaps, be made with the seasonal geocentre displacements $\text{SLR}s^{\text{CM}}$ between 2004 and 2012 that we

obtain from the time-series derived by Cheng et al. (2010) from SLR tracking to LAGEOS satellites (http://grace.jpl.nasa.gov/data/get-data/geocenter—last accessed 30 May 2016). The annual and semi-annual components of $\text{SLR}s^{\text{CM}}$ are also included in Table 2. For both $s^{\text{CM}}$ and $\text{SLR}s^{\text{CM}}$, the z-component is larger than the other two components, for either the annual or the semi-annual displacements. Again, if we take into account the errors, up to 36 per cent, associated to the inversion of the gravity data, most of the amplitudes of $s^{\text{CM}}$ and $\text{SLR}s^{\text{CM}}$ overlap.

Until today, all the estimates of the geocentre motion have relied on the assumption that it is caused by surface loading. One could free oneself from that hypothesis by combining precise positioning measurements with observations of the gravity variations and computing $s^{\text{CM}}$ in the CE frame by means of eq. (12), provided the degree-1 station displacement can also be determined without assuming that it is due to a surface load.

# 6 CONCLUSIONS

We have investigated the constraints that can be brought on the geocentre motion by ground observations with SGs.

Two key results of this paper are contained in the formulae (12) and (26). They are relations between the displacement of the centre of mass CM of the entire Earth and the degree-1 Lagrangian gravity variation at the surface. Eq. (12) is valid in any system of reference and requires that the radial ground displacement be known, either theoretically or observationally. Eq. (26) is valid in the CE frame when the Earth is subject to a surface load. It does not require any observation other than the gravity variation, but a spherical Earth model must be adopted in order to determine the load Love number $h'_1^{\text{CE}}$ which is contained in the formula. Interestingly, neither (12) nor (26) contains the horizontal displacement of the station.

We have stacked and inverted synthetic gravity signals to study the influence of the density and geographical location of the gravity stations on the estimate of the degree-1 gravity variation. We have then selected six stations of the GGP (Hinderer & Crossley 2000) that are symmetrically distributed about the equator and Greenwich meridian. We have inverted data from 2004 until 2012 to retrieve the degree-1 gravity variation and, using eq. (26) and the load Love number $h'_1^{\text{CE}} = -0.288$ for the elastic PREM (Dziewonski & Anderson 1981), we have obtained time-series for the three components of the displacement of CM in a CE frame. The fit of the series with annual and semi-annual signals has provided displacements from 0.5 to 2.3 mm, which are in good agreement with estimates from geodetic techniques such as SLR or GPS that also rely on the assumption that the geocentre motion is due to surface loading.
ACKNOWLEDGEMENTS

YR gratefully acknowledges financial support from the CNRS-
INSU *Programme National de Planétologie*. This work was ini-
tiated while AM was supported by an Australian Research Council
Super Science Fellowship (FS110200045).

REFERENCES

Altamimi, Z., Collilieux, X., Legrand, J., Garayt, B. & Boucher, C., 2007. ITRF2005: a new release of the International Terrestrial Reference Frame based on time series of station positions and Earth Orientation Parameters, *J. geophys. Res.*, **112**, B09401, doi:10.1029/2007JB004949.

Altamimi, Z., Collilieux, X. & Métivier, L., 2011. ITRF2008: an improved solution of the International Terrestrial Reference Frame, *J. Geod.*, **85**, 457–473.

Blewitt, G., Lavallée, D., Clarke, P. & Nurutdinov, K., 2001. A new global mode of Earth deformation: seasonal cycle detected, *Science*, **294**, 2342–2345.

Blewitt, G., 2003. Self-consistency in reference frames, geocenter defini-
tion, and surface loading of the solid Earth, *J. geophys. Res.*, **108**(B2),
doi:10.1029/2002JB002082.

Cheng, M.K., Tapley, B.D. & Ries, J.C., 2010. Geocenter variations from analysis of SLR data, in *IAG Commission 1 Symposium 2010*, pp. 19–25. Reference Frames for Applications in Geosciences (REFAG2010), Marne-La-Vallée, France, 4–8 October 2010.

Collilieux, X., Altamimi, Z., Coulot, D., van Dam, T. & Ray, J., 2010. Impact of loading effects on determination of the International Terrestrial Reference Frame, *Adv. Space Res.*, **45**, 144–154.

Courtier. N. et al., 2000. Global superconducting gravimeter observations and the search for the translational modes of the inner core, *Phys. Earth planet. Inter.*, **117**, 3–20.

Dahlen, F.A. & Tromp, J., 1998. *Theoretical Global Seismology*, Princeton Univ. Press.

Dehant, V., Defraigne, P. & Wahr, J., 1999. Tides for a convective Earth, *J. geophys. Res.*, **104**(B1), 1035–1058.

de Linage, C., Hinderer, J. & Boy, J.-P., 2009. Variability of the gravity-to-height ratio due to surface loads, *Pure appl. Geophys.*, **166**, 1217–1245.

Dobslaw, H., Bergmann-Wolrf, I., Dill, R., Forootan, E., Klemann, V., Kusche, J. & Sasgen, I., 2015. The updated ESA Earth System Model for future gravity mission simulation studies, *J. Geod.*, **89**, 505–513.

Dong, D., Dickey, J.O., Chao, Y. & Cheng, M.K., 1997. Geocenter variations caused by atmosphere, ocean and surface ground water, *Geophys. Res. Lett.*, **24**(15), 1867–1870.

Dong, D., Yunck, T. & Heflin, M., 2003. Origin of the interna-
tional terrestrial reference frame, *J. geophys. Res.*, **108**(B4), 2200,
doi:10.1029/2002JB002035.

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

Feissel-Vernier, M., Le Bail, K., Berio, P., Coulot, D., Ramillien, G. & Valette, J.-J., 2006. Geocenter motion measured with DORIS and SLR, and predicted by geophysical models, *J. Geod.*, **80**(8–11), 637–648.

Greff-Lefftz, M., 2000. Secular variation of the geocenter, *J. geophys. Res.*, **105**(B11), 25 685–25 692.

Hinderer, J. & Crossley, D., 2000. Time variations in gravity and inferences on the Earth's structure and dynamics, *Surv. Geophys.*, **21**, 1–45.

Hinderer, J., Legros, H. & Crossley, D., 1991. Global Earth dynamics and induced gravity changes, *J. geophys. Res.*, **96**, 20 257–20 265.

Hoffmann-Wellenhoff, B. & Moritz, H., 2005. *Physical Geodesy*, Springer-Verlag.

Longman, I.M., 1963. A Green's function for determining the deformation of the Earth under surface mass loads: 1. Computations and numerical results, *J. geophys. Res.*, **68**, 485–496.

Matsumoto, K., Takanezawa, T. & Ooe, M., 2000. Ocean tide models developed by assimilating TOPEX/POSEIDON altimeter data into hydrodynamical model: a global model and a regional model around, *Japan. J. Oceanogr.*, **56**, 567–581.

McCarthy, D.D. & Petit, G., 2004. IERS Conventions (2003), IERS Technical Note 32, Verlag des Bundesamts für Kartographie und Geodäsie Frankfurt am Main, p. 179.

Métivier, L., Greff-Lefftz, M. & Altamimi, Z., 2010. On secular geocenter motion: the impact of climate changes, *Earth planet. Sci. Lett.*, **296**, 360–366.

Petit, G. & Luzum, B., 2010. IERS Conventions (2010), IERS Technical Note 36, Verlag des Bundesamts für Kartographie und Geodäsie Frankfurt am Main, p. 179.

Plag, H.P. et al., 2009. The goals, achievements, and tools of modern geodesy, in *Global Geodetic Observing System: Meeting the Requirements of a Global Society on a Changing Planet in 2020*, eds Plag, H.-P. & Pearlman, M., Springer-Verlag, Berlin Heidelberg.

Plag, H.-P., Kreemer, C. & Hammond, W., 2007. Combination of GPS-observed vertical motion with absolute gravity changes constrain the tie between reference frame origin and Earth center of mass, Poster presented at the *National EarthScope Meeting*, Monterey, CA, 26–30 March 2007.

Rebischung, P., Altamimi, Z. & Springer, T., 2014. A collinearity diagnosis of the GNSS geocenter determination, *J. Geod.*, **88**, 65–85.

Rietbroek, R., Fritsche, M., Brunnabend, S.-E., Daras, I., Kusche, J., Schröter, J., Flechtner, F. & Dietrich, R., 2012. Global surface mass from a new combination of GRACE, modelled OBP and reprocessed GPS data, *J. Geodyn.*, **59–60**, 64–71.

Rosat, S. & Hinderer, J., 2011. Noise levels of superconducting gravimeters: updated comparison and time stability, *Bull. seism. Soc. Am.*, **101**, doi:10.1785/0120100217.

Rosat, S., Hinderer, J., Crossley, D. & Rivera, L., 2003. The search for the slichter mode: comparison of noise levels of superconducting gravimeters and investigation of a stacking, *Phys. Earth planet. Int.*, **140**, 183–202.

Rülke, A., Dietrich, R., Fritsche, M., Rothacher, M. & Steigenberger, P., 2008. Realization of the Terrestrial Reference System by a re-processed global GPS network, *J. geophys. Res.*, **113**, B08409, doi:10.1029/2007JB005231.

van Camp, M., Williams, S.D.P. & Francis, O., 2005. Uncertainty of absolute gravity measurements, *J. geophys. Res.*, **58**, B05406, doi:10.1029/2004JB003497.

Wu, X., Ray, J. & van Dam, T., 2012. Geocenter motion and its geodetic and geophysical implications, *J. Geodyn.*, **58**, 44–61.