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

# Improving GOCE cross-track gravity gradients

Christian Siemes<fnref n="1" />

Received: 12 January 2017 / Accepted: 25 May 2017 / Published online: 5 July 2017
© The Author(s) 2017. This article is an open access publication

**Abstract** The GOCE gravity gradiometer measured highly accurate gravity gradients along the orbit during GOCE’s mission lifetime from March 17, 2009, to November 11, 2013. These measurements contain unique information on the gravity field at a spatial resolution of 80 km half wavelength, which is not provided to the same accuracy level by any other satellite mission now and in the foreseeable future. Unfortunately, the gravity gradient in cross-track direction is heavily perturbed in the regions around the geomagnetic poles. We show in this paper that the perturbing effect can be modeled accurately as a quadratic function of the non-gravitational acceleration of the satellite in cross-track direction. Most importantly, we can remove the perturbation from the cross-track gravity gradient to a great extent, which significantly improves the accuracy of the latter and offers opportunities for better scientific exploitation of the GOCE gravity gradient data set.

**Keywords** Satellite gravity gradiometry · GOCE mission · Accelerometry · Calibration

## 1 Introduction

The Gravity field and Steady-state Ocean Circulation Explorer (GOCE) was launched on March 17, 2009, and re-entered on November 11, 2013. The objective of the GOCE mission is the determination of the mean gravity field with

an accuracy of 1–2 cm for the geoid and 1 mGal for the gravity anomaly field, and achieving both at a spatial resolution of 100 km half wavelength (Drinkwater et al. 2007). In order to reach this ambitious objective, GOCE carried two complementary measurement systems on board. The first is the satellite-to-satellite tracking (SST) instrument that was used to determine the position and velocity of the satellite, from which the gravity field can be retrieved. The second is the gravity gradiometer that measured the gravity gradients along GOCE’s orbit. The instruments contribute equally in the gravity field retrieval at a spatial resolution of approximately 1000 km, where the SST instrument is more accurate at lower and the gravity gradiometer is more accurate at higher spatial resolution.

Since the gravity gradiometer required a “quiet” environment for providing excellent performance, the satellite was equipped with a drag-free and attitude control system that kept the satellite Sun- and nadir-pointing and compensated non-gravitational forces in flight direction caused mainly by drag. The elements of that system included the gradiometer for measuring non-gravitational accelerations, three star sensors for measuring the attitude of the satellite, an ion propulsion system for counteracting non-gravitational forces in flight direction and three magnetic torquers for controlling the attitude (Floberghagen et al. 2011). In addition, the satellite was designed to be stable under varying thermal loads, had no moving parts in order to avoid microvibrations and its shape was optimized for minimizing the effect of aerodynamic drag.

The gravity gradiometer consisted of three pairs of ultra-sensitive three-axis accelerometers, where the accelerometers forming a pair were mounted at the ends of three orthogonal gradiometer arms, separated by half a meter. Each accelerometer measured the control voltages that were required to levitate a proof-mass at the center of a slightly

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

<footnote>Christian Siemes
christian.siemes@esa.int</footnote>

<footnote n="1">¹ RHEA for ESA - European Space Agency, ESTEC, Mission Science Division, Noordwijk, The Netherlands</footnote>
larger cage. The control voltages were applied to eight electrodes located on the inner walls of the cage of each accelerometer and were representative for the acceleration of the proof-mass relative to the cage. The proof-masses had square cuboid shape of 4 cm × 4 cm × 1 cm, which allowed for ground testing but made the accelerometer axis that was normal to 4 cm × 4 cm face less sensitive. This prevented the gravity gradiometer from measuring all elements of the gravity gradient tensor with high accuracy. Due to the orientation of the individual accelerometers, only the elements $V_{xx}$, $V_{yy}$, $V_{zz}$ and $V_{xz}$ of the gravitational tensor were measured with high accuracy (Rummel et al. 2011).

GOCE gravity gradients are used as input for many state-of-the-art gravity field models (Pail et al. 2011; Hashemi Farahani et al. 2013; Schall et al. 2014; Brockmann et al. 2014; Bruinsma et al. 2014), contributing significantly to those models' high accuracy and spatial resolution, with downstream applications in oceanography, geophysics and geodesy (Knudsen et al. 2011; Gruber et al. 2012; Fuchs et al. 2013; Becker et al. 2014; Hirt 2014). Beyond that, the direct use of gravity gradients in geophysical applications is currently evolving (Van der Meijde et al. 2015; Bouman et al. 2016). The GOCE mission provided a new and unique gravity gradient data set that is today the most accurate spaceborne measurement of the gravity field at high spatial resolution and is likely to remain so at least for the next decade.

The measured gravity gradient $\vec{V}_{yy}$, which points roughly into cross-track direction, is unfortunately subject to large perturbations in the regions around the geomagnetic poles, which was noted earlier, e.g., in Bouman et al. (2011). In particular, Petersheim et al. (2011) show that the perturbation cannot be removed by adjusting the scale factors for the accelerometer pair on the cross-track gradiometer arm. Though it was observed that magnitude of the perturbation is linked to the magnitude of non-gravitational acceleration (sum of drag, winds, and radiation pressure) in cross-track direction (Siemes et al. 2012), the root cause of the perturbation remained unclear until now. Ince and Pagiatakis (2016) find a correlation between the occurrence of the perturbation and more active space weather conditions. They suggest to perform an empirical correction based on a response analysis between the perturbation in $\vec{V}_{yy}$ and ionospheric dynamics represented by Poynting energy flux. Their empirical correction model is limited to regions where Poynting energy flux measurements are available, i.e., a global correction is not possible.

In this paper, we show that the perturbing effect can be accurately modeled by a quadratic function of the non-gravitational acceleration in cross-track direction, which implies that the root cause of the perturbation rests within the conversion from the control voltages to acceleration. As described in Lenoir et al. (2011), a quadratic factor occurs in the functional relation of the control voltage to the accel-

eration, when the proof mass position slightly deviates from the center of the accelerometer. For that reason, the method described in Frommknecht et al. (2011) was used to nullify the quadratic factor. This was achieved by measuring the quadratic factor and then adjusting the proof mass position accordingly. Therefore, it is surprising to find that the relation between acceleration and control voltages is best described by a quadratic function for at least one of the in total 18 accelerometer axes. More investigations are needed for understanding how a quadratic factor can still create a significant perturbation despite the efforts invested in physically nullifying it. Nevertheless, we demonstrate that the perturbation can be removed to a great extent from the measured gravity gradient $\vec{V}_{yy}$, which increases its accuracy significantly.

We specify the input data for the investigation presented in this paper in Sect. 2. The method is described in Sect. 3 and the results are shown in Sect. 4. We conclude the paper with a discussion of the results in Sect. 5, which includes recommendations for future investigations.

## 2 Input data

Our main input data are the gravity gradients and calibrated accelerations, both provided in the gradiometer reference frame and contained in the EGG_NOM_1b data products (Frommknecht et al. 2011). In addition, we use the gradiometer inertial attitude quaternions that describe the rotation from the gradiometer reference frame to the inertial reference frame and are also contained in the EGG_NOM_1b data products. The rotation from the inertial reference frame to the Earth-fixed reference frame is given in form of quaternions in the SST_PSO_2_ data products, which also contain the satellite positions of the reduced-dynamic orbit in the Earth-fixed reference frame. Furthermore, we use the GOCO03s gravity field model, which is in particular based on GOCE and GRACE mission data (Mayer-Gürr et al. 2012). It should be noted that we would obtain the same results with any other gravity field model that is based on GOCE gravity gradients.

From these input data, we calculate first the model gravity gradient $V_{yy}^{GOCO03s}$ from the GOCO03s model in the local north-oriented reference frame (The European GOCE Gravity Consortium EGG-C 2014). Then, we use the quaternions mentioned before to transform the model gravity gradients to the gradiometer reference frame as described in Rummel et al. (2011). The gradiometer reference frame is aligned with the satellite axes, i.e., the x, y and z-axis are pointing approximately into flight, nadir and cross-track direction, respectively. Finally, we calculate the residual $\Delta V_{yy}$ between the measured gravity gradient $\vec{V}_{yy}$ and the model gravity gradient according to
$$ \Delta V_{yy} = \bar{V}_{yy} - V_{yy}^{\text{GOCO03s}}. \qquad (1) $$

One could argue that the GOCO03s model might be affected by the perturbation in the measured gravity gradient $\bar{V}_{yy}$ since that was used in the gravity field retrieval. Therefore, we calculated the model gravity gradient $V_{yy}^{\text{ITG-Grace2010s}}$ in the same way also from the ITG-Grace2010s model (Mayer-Gürr et al. 2010) which is based only on GRACE mission data and therefore independent of the perturbing effects, and compared it to that calculated from the GOCO03s model. The RMS of the difference $V_{yy}^{\text{GOCO03s}} - V_{yy}^{\text{ITG-Grace2010s}}$ was only 1 mE, which we consider to be negligible for the investigation presented in this paper in view of the noise level of $\bar{V}_{yy}$, which is 10 mE/Hz$^{1/2}$ within the frequency range 5–100 mHz and larger outside (Rummel et al. 2011).

### 3 Method

It is instructive to discuss the signal content of the accelerometer data. A perfectly calibrated accelerometer measures the acceleration

$$ \mathbf{a}_i = -(\mathbf{V} - \mathbf{\Omega}^2 - \dot{\mathbf{\Omega}})\mathbf{r}_i + \mathbf{d} \qquad (2) $$

where $i = 1, \dots, 6$ indicates the accelerometer. Figure 1 shows how the accelerometers are arranged in the gradiometer reference frame. The acceleration $\mathbf{a}_i$ contains the signal of the gravity gradient $\mathbf{V}$, the centrifugal acceleration $\mathbf{\Omega}^2$, the angular acceleration of the satellite $\dot{\mathbf{\Omega}}$, and the non-gravitational acceleration of the satellite $\mathbf{d}$ where

<figure><img src="image_3.png" /><figcaption>Fig. 1 Arrangement of the six accelerometers in the gradiometer reference frame (GRF). Dotted lines indicate the less sensitive axis of the accelerometers</figcaption></figure>

$$ \mathbf{a}_i = \begin{bmatrix} a_{ix} \\ a_{iy} \\ a_{iz} \end{bmatrix}, \qquad (3) $$

$$ \mathbf{V} = \begin{bmatrix} V_{xx} & V_{xy} & V_{xz} \\ V_{xy} & V_{yy} & V_{yz} \\ V_{xz} & V_{yz} & V_{zz} \end{bmatrix}, \qquad (4) $$

$$ \mathbf{\Omega}^2 = \begin{bmatrix} -\omega_y^2 - \omega_z^2 & \omega_x \omega_y & \omega_x \omega_z \\ \omega_x \omega_y & -\omega_x^2 - \omega_z^2 & \omega_y \omega_z \\ \omega_x \omega_z & \omega_y \omega_z & -\omega_x^2 - \omega_y^2 \end{bmatrix}, \qquad (5) $$

$$ \dot{\mathbf{\Omega}} = \begin{bmatrix} 0 & -\dot{\omega}_z & \dot{\omega}_y \\ \dot{\omega}_z & 0 & -\dot{\omega}_x \\ -\dot{\omega}_y & \dot{\omega}_x & 0 \end{bmatrix} \qquad (6) $$

and

$$ \mathbf{d} = \begin{bmatrix} d_x \\ d_y \\ d_z \end{bmatrix}. \qquad (7) $$

$\mathbf{r}_i$ is the vector from the satellite center of mass to the center of the $i$-th accelerometer. It is practical to transform the acceleration $\mathbf{a}_i$ into differential mode accelerations

$$ \mathbf{a}_{dij} = (\mathbf{a}_i - \mathbf{a}_j)/2 = -\frac{1}{2}(\mathbf{V} - \mathbf{\Omega}^2 - \dot{\mathbf{\Omega}})(\mathbf{r}_i - \mathbf{r}_j) \qquad (8) $$

and common mode accelerations

$$ \mathbf{a}_{cij} = (\mathbf{a}_i + \mathbf{a}_j)/2 = \mathbf{d} \qquad (9) $$

for the accelerometer pairs $ij = 14, 25$ and $36$ where

$$ \mathbf{r}_1 - \mathbf{r}_4 = \begin{bmatrix} L_x \\ 0 \\ 0 \end{bmatrix}, \quad \mathbf{r}_2 - \mathbf{r}_5 = \begin{bmatrix} 0 \\ L_y \\ 0 \end{bmatrix} \quad \text{and} \quad \mathbf{r}_3 - \mathbf{r}_6 = \begin{bmatrix} 0 \\ 0 \\ L_z \end{bmatrix}. \qquad (10) $$

The quantities $L_x$, $L_y$ and $L_z$ are the distances between the accelerometers of pair 14, 25 and 36, respectively. Rearranging Eq. (8) gives the gravity gradients

$$ V_{xx} = -2a_{d14x}/L_x - \omega_y^2 - \omega_z^2, \qquad (11) $$

$$ V_{yy} = -2a_{d25y}/L_y - \omega_x^2 - \omega_z^2, \qquad (12) $$

$$ V_{zz} = -2a_{d36z}/L_z - \omega_x^2 - \omega_y^2, \qquad (13) $$

$$ V_{xy} = -a_{d25x}/L_y - a_{d14y}/L_x + \omega_x \omega_y, \qquad (14) $$

$$ V_{xz} = -a_{d14z}/L_x - a_{d36x}/L_z + \omega_x \omega_z \qquad (15) $$

and

$$ V_{yz} = -a_{d36y}/L_z - a_{d25z}/L_y + \omega_y \omega_z \qquad (16) $$
where $\omega_x$, $\omega_y$ and $\omega_z$ are the angular rates about the along-track, cross-track, and nadir axis of the gradiometer reference frame, respectively.

As mentioned before, Eq. (2) holds for the measurement $\mathbf{a}_i$ of a perfectly calibrated accelerometer. A real accelerometer is, however, subject to small imperfections and measurement noise. We relate the measurement $\hat{\mathbf{a}}_i$ of a real accelerometer to $\mathbf{a}_i$ by

$$ \hat{\mathbf{a}}_i = \hat{\mathbf{b}}_i + \mathbf{M}_i \mathbf{a}_i + \mathbf{K}_i \mathbf{a}_i^2 + \hat{\mathbf{n}}_i \quad (17) $$

where $\hat{\mathbf{b}}_i$ is the bias in the accelerometer measurements, $\mathbf{M}_i$ is the calibration matrix for a single accelerometer, $\mathbf{K}_i$ is a diagonal matrix holding the quadratic factors, and $\hat{\mathbf{n}}_i$ is measurement noise. The vectors and matrices are defined as

$$ \hat{\mathbf{a}}_i = \begin{bmatrix} \hat{a}_{ix} \\ \hat{a}_{iy} \\ \hat{a}_{iz} \end{bmatrix}, \quad (18) $$

$$ \hat{\mathbf{b}}_i = \begin{bmatrix} \hat{b}_{ix} \\ \hat{b}_{iy} \\ \hat{b}_{iz} \end{bmatrix}, \quad (19) $$

$$ \mathbf{M}_i = \begin{bmatrix} M_{ixx} & M_{ixy} & M_{ixz} \\ M_{iyx} & M_{iyy} & M_{iyz} \\ M_{izx} & M_{izy} & M_{izz} \end{bmatrix}, \quad (20) $$

$$ \mathbf{K}_i = \begin{bmatrix} k_{ix} & 0 & 0 \\ 0 & k_{iy} & 0 \\ 0 & 0 & k_{iz} \end{bmatrix}, \quad (21) $$

and

$$ \hat{\mathbf{n}}_i = \begin{bmatrix} \hat{n}_{ix} \\ \hat{n}_{iy} \\ \hat{n}_{iz} \end{bmatrix}. \quad (22) $$

For the square of vector $\mathbf{a}_i$, we use the convention

$$ \mathbf{a}_i^2 = \begin{bmatrix} a_{ix}^2 \\ a_{iy}^2 \\ a_{iz}^2 \end{bmatrix} \quad (23) $$

throughout this paper.

The calibration of the gravity gradiometer is described in detail in Frommknecht et al. (2011) and comprises two separate procedures. The first is the determination of the quadratic factors in $\mathbf{K}_i$ and the adjustment of the proof mass positions such that $\mathbf{K}_i \approx \mathbf{0}$. The quadratic factors were thus physically nullified and are consequently not considered in the Level 1b processing. The second procedure is the determination of the calibration matrices $\mathbf{M}_i$ from data collected during dedicated satellite maneuvers, which were performed typically every two months during mission lifetime. The calibration

matrices are taken into account in the Level 1b processing of the gradiometer data, whereas the biases $\hat{\mathbf{b}}_i$ are ignored. It should be noted that an extensive analysis presented in Siemes (2011) shows that the perturbing effect in $\tilde{V}_{yy}$ cannot be modeled by the calibration matrices $\mathbf{M}_i$. Therefore, we focus in this investigation on the quadratic factors $\mathbf{K}_i$.

In the same way that we calculated differential and common mode accelerations in Eqs. (8) and (9), we calculate now measured differential and common mode accelerations. Using Eq. (17) leads to

$$ \begin{bmatrix} \hat{\mathbf{a}}_{dij} \\ \hat{\mathbf{a}}_{cij} \end{bmatrix} = \frac{1}{2} \begin{bmatrix} \hat{\mathbf{a}}_i - \hat{\mathbf{a}}_j \\ \hat{\mathbf{a}}_i + \hat{\mathbf{a}}_j \end{bmatrix} \quad (24) $$

$$ \begin{aligned} &= \begin{bmatrix} \hat{\mathbf{b}}_{dij} \\ \hat{\mathbf{b}}_{cij} \end{bmatrix} + \begin{bmatrix} \mathbf{C}_{ij} & \mathbf{D}_{ij} \\ \mathbf{D}_{ij} & \mathbf{C}_{ij} \end{bmatrix} \begin{bmatrix} \mathbf{a}_{dij} \\ \mathbf{a}_{cij} \end{bmatrix} \\ &\quad + \frac{1}{2} \begin{bmatrix} \mathbf{K}_i - \mathbf{K}_j \\ \mathbf{K}_i & \mathbf{K}_j \end{bmatrix} \begin{bmatrix} \mathbf{a}_i^2 \\ \mathbf{a}_j^2 \end{bmatrix} + \begin{bmatrix} \hat{\mathbf{n}}_{dij} \\ \hat{\mathbf{n}}_{cij} \end{bmatrix}. \end{aligned} \quad (25) $$

The vectors $\hat{\mathbf{b}}_{cij}$ and $\hat{\mathbf{b}}_{dij}$ are the bias in the measured common and differential mode acceleration, respectively, the matrices $\mathbf{C}_{ij}$ and $\mathbf{D}_{ij}$ are common and differential calibration matrices, respectively, and the vectors $\hat{\mathbf{n}}_{cij}$ and $\hat{\mathbf{n}}_{dij}$ are noise in the measured common and differential mode acceleration, respectively, where

$$ \hat{\mathbf{b}}_{cij} = (\hat{\mathbf{b}}_i + \hat{\mathbf{b}}_j)/2, \quad (26) $$

$$ \hat{\mathbf{b}}_{dij} = (\hat{\mathbf{b}}_i - \hat{\mathbf{b}}_j)/2, \quad (27) $$

$$ \mathbf{C}_{ij} = (\mathbf{M}_i + \mathbf{M}_j)/2, \quad (28) $$

$$ \mathbf{D}_{ij} = (\mathbf{M}_i - \mathbf{M}_j)/2, \quad (29) $$

$$ \hat{\mathbf{n}}_{cij} = (\hat{\mathbf{n}}_i + \hat{\mathbf{n}}_j)/2, \quad (30) $$

and

$$ \hat{\mathbf{n}}_{dij} = (\hat{\mathbf{n}}_i - \hat{\mathbf{n}}_j)/2. \quad (31) $$

The calibrated common and differential mode acceleration $\tilde{\mathbf{a}}_{cij}$ and $\tilde{\mathbf{a}}_{dij}$, respectively, which are contained in the EGG_NOM_1b data products and used in the Level 1b processing for generating the gravity gradients, are related to the measured accelerations by

$$ \begin{bmatrix} \tilde{\mathbf{a}}_{dij} \\ \tilde{\mathbf{a}}_{cij} \end{bmatrix} = \mathbf{M}_{ij} \begin{bmatrix} \hat{\mathbf{a}}_{dij} \\ \hat{\mathbf{a}}_{cij} \end{bmatrix} \quad (32) $$

where

$$ \mathbf{M}_{ij} = \begin{bmatrix} \mathbf{C}_{ij} & \mathbf{D}_{ij} \\ \mathbf{D}_{ij} & \mathbf{C}_{ij} \end{bmatrix}^{-1} \quad (33) $$

is the inverse calibration matrix as defined in (Frommknecht et al. 2011). We transform Eq. (25) using the definitions in Eqs. (32) and (33) to
$$ \begin{bmatrix} \bar{a}_{dij} \\ \bar{a}_{cij} \end{bmatrix} = \begin{bmatrix} \bar{b}_{dij} \\ \bar{b}_{cij} \end{bmatrix} + \begin{bmatrix} a_{dij} \\ a_{cij} \end{bmatrix} + \frac{1}{2} M_{ij} \begin{bmatrix} K_i - K_j \\ K_i \quad K_j \end{bmatrix} \begin{bmatrix} a_j^2 \\ a_j^2 \end{bmatrix} + \begin{bmatrix} \bar{n}_{dij} \\ \bar{n}_{cij} \end{bmatrix} \quad (34) $$

where $\bar{b}_{cij}$ and $\bar{b}_{dij}$ are the bias in the calibrated common and differential mode acceleration, respectively, and $\bar{n}_{cij}$ and $\bar{n}_{dij}$ are the noise in the calibrated common and differential mode acceleration, respectively. Similar to Eq. (32), the bias and noise in the calibrated common and differential mode acceleration are related to those in the measured common and differential mode acceleration by

$$ \begin{bmatrix} \bar{b}_{dij} \\ \bar{b}_{cij} \end{bmatrix} = M_{ij} \begin{bmatrix} \hat{b}_{dij} \\ \hat{b}_{cij} \end{bmatrix} \quad (35) $$

and

$$ \begin{bmatrix} \bar{n}_{dij} \\ \bar{n}_{cij} \end{bmatrix} = M_{ij} \begin{bmatrix} \hat{n}_{dij} \\ \hat{n}_{cij} \end{bmatrix}, \quad (36) $$

respectively. As discussed before, we expect that the quadratic factors are close to zero. Noting that the inverse calibration matrices are close to an identity matrix, i.e., $M_{ij} \approx I$, we approximate

$$ M_{ij} \begin{bmatrix} K_i - K_j \\ K_i \quad K_j \end{bmatrix} \approx \begin{bmatrix} K_i - K_j \\ K_i \quad K_j \end{bmatrix} \quad (37) $$

without introducing significant approximation errors in Eq. (34), which then reads

$$ \begin{bmatrix} \bar{a}_{dij} \\ \bar{a}_{cij} \end{bmatrix} = \begin{bmatrix} \bar{b}_{dij} \\ \bar{b}_{cij} \end{bmatrix} + \begin{bmatrix} a_{dij} \\ a_{cij} \end{bmatrix} + \frac{1}{2} \begin{bmatrix} K_i - K_j \\ K_i \quad K_j \end{bmatrix} \begin{bmatrix} a_j^2 \\ a_j^2 \end{bmatrix} + \begin{bmatrix} \bar{n}_{dij} \\ \bar{n}_{cij} \end{bmatrix}. \quad (38) $$

The biases and quadratic factors in Eq. (38) enter the calculation of the gravity gradients as in Eqs. (11)–(16) directly by using the calibrated acceleration $\bar{a}_{dij}$ in place of $a_{dij}$ and indirectly through the angular rates, which are reconstructed from the calibrated acceleration $\bar{a}_{dij}$ together with star sensor data. For the angular rate reconstruction, the method described in Stummer et al. (2011) is deployed, which implicitly eliminates the biases $\bar{b}_{dij}$ by means of a high-pass filter. We assume in this paper that the quadratic factors that enter the angular rate reconstruction are zero, so that the reconstructed angular rates are not affected by biases or quadratic factors. Thus, we assume in this paper that the gravity gradients are only affected by biases and quadratic factors that enter Eqs. (11)–(16) directly through the calibrated acceleration $\bar{a}_{dij}$, though more investigations are needed to verify this assumption. In particular, we obtain the measured gravity gradient

$$ \bar{V}_{yy} = -2\bar{a}_{d25y}/L_y - \bar{\omega}_x^2 - \bar{\omega}_z^2 \quad (39) $$

where $\bar{\omega}_x$ and $\bar{\omega}_z$ are the angular rates resulting from the angular rate reconstruction. The reconstructed angular rates contain noise, which propagates to Eq. (39) through the centrifugal terms $\bar{\omega}_x^2$ and $\bar{\omega}_z^2$. We define $n_{\Omega_{yy}^2}$ as sum of the noise in the centrifugal terms such that

$$ \bar{\omega}_x^2 + \bar{\omega}_z^2 = \omega_x^2 + \omega_z^2 + n_{\Omega_{yy}^2} \quad (40) $$

and rewrite Eq. (39) to

$$ \bar{V}_{yy} = -2\bar{a}_{d25y}/L_y - \omega_x^2 - \omega_z^2 - n_{\Omega_{yy}^2}. \quad (41) $$

Next, we use Eq. (38) to replace the calibrated acceleration in Eq. (41) to obtain

$$ \bar{V}_{yy} = -2(\bar{b}_{d25y} + a_{d25y} + (k_{2y}a_{2y}^2 - k_{5y}a_{5y}^2))/2 + \bar{n}_{d25y}/L_y - \omega_x^2 - \omega_z^2 - n_{\Omega_{yy}^2}. \quad (42) $$

Subtracting the model gravity gradient $V_{yy}^{GOCO03s}$ as in Eq. (1) simplifies Eq. (42) to

$$ \Delta V_{yy} = -2\bar{b}_{d25y}/L_y - (k_{2y}a_{2y}^2 - k_{5y}a_{5y}^2)/L_y - 2\bar{n}_{d25y}/L_y - n_{\Omega_{yy}^2} - n_{V_{yy}^{GOCO03s}} \quad (43) $$

where $n_{V_{yy}^{GOCO03s}}$ is the noise in $V_{yy}^{GOCO03s}$, which is due to model errors and noise in the orbit as well as the attitude quaternions that were used to transform the gravity gradients to the gradiometer reference frame. In order to shorten Eq. (43), we introduce

$$ n_{yy} = -2\bar{n}_{d25y}/L_y - n_{\Omega_{yy}^2} - n_{V_{yy}^{GOCO03s}} \quad (44) $$

as the sum of all noise terms such that

$$ \Delta V_{yy} = -2\bar{b}_{d25y}/L_y - (k_{2y}a_{2y}^2 - k_{5y}a_{5y}^2)/L_y + n_{yy}. \quad (45) $$

Now it is important to know that $d_y$ is causing by far the largest signal dynamics (variation about mean) in the acceleration $a_{2y}$ and $a_{5y}$ as presented in Eq. (2). On April 1, 2012, for instance, the signal dynamics in $d_y$ are in the order of 1 µm/s², whereas $V_{yy}$, $\omega_x^2$, and $\omega_z^2$ cause acceleration signal dynamics of only 5 nm/s². We approximate therefore

$$ a_{2y} \approx a_{5y} \approx d_y = a_{c25y}. \quad (46) $$

Inserting Eq. (46) into (45) gives

$$ \Delta V_{yy} = -2\bar{b}_{d25y}/L_y - 2k_{d25y}a_{c25y}^2/L_y + n_{yy} \quad (47) $$

where $k_{d25y}$ is the differential quadratic factor defined as

$$ k_{d25y} = (k_{2y} - k_{5y})/2. \quad (48) $$
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Bandpass-filtered gravity gradient residuals $\Delta \tilde{V}_{yy}$ of present Level 1b gravity gradients (top panel), bandpass-filtered common mode accelerations $\tilde{a}_{c25y}$ (middle panel), and bandpass-filtered squared common mode accelerations $\tilde{a}_{c25y}^2$ (bottom panel)</figcaption></figure>

We find by rearranging Eq. (38) and using the approximation in Eq. (46)

$$a_{c25y}^2 = (\bar{a}_{c25y} - \bar{b}_{c25y} - (k_{2y} + k_{5y})a_{c25y}^2/2 - \bar{n}_{c25y})^2 \\ \approx \bar{a}_{c25y}^2 - 2\bar{b}_{c25y}\bar{a}_{c25y} + \bar{b}_{c25y}^2 \quad (49)$$

where we dropped all products that contain three or more factors or the noise term $\bar{n}_{c25y}$ since the acceleration and biases are small and the noise term is even smaller. In order to provide an indication of the size of $\bar{a}_{c25y}$ and $\bar{b}_{c25y}$, we refer again to the example of April 1, 2012, where $\bar{a}_{c25y}$ is in the order of 1-2 µm/s².

Inserting Eq. (49) into Eq. (47) gives

$$\Delta V_{yy} = -\frac{2}{L_y}\bar{b}_{d25y} - \frac{2}{L_y}k_{d25y}\bar{a}_{c25y}^2 + \frac{4}{L_y}k_{d25y}\bar{b}_{c25y}\bar{a}_{c25y} \\ -\frac{2}{L_y}k_{d25y}\bar{b}_{c25y}^2 + n_{yy}, \quad (50)$$

which represents our model for the perturbation in $\tilde{V}_{yy}$.

Equation (50) demonstrates that a nonzero quadratic factor $k_{d25y}$ causes a perturbation in $\tilde{V}_{yy}$ that is a quadratic function of the common mode accelerations $\bar{a}_{c25y}$. This means that strong dynamics in $\bar{a}_{c25y}$ in the regions near the geomagnetic poles give perturbations in $\tilde{V}_{yy}$ in the same region, which was noted earlier, e.g., in Bouman et al. (2011), Siemes et al. (2012) and Ince and Pagiatakis (2016). Since the quadratic factors are expected to be zero, their effects were not investigated until now. Instead, Bouman et al. (2011), Siemes (2011)

and Petersheim et al. (2011) focussed their investigations on effects of the scale factor $s_{d25y}$, which is the element in the second row and second column of the differential calibration matrix $D_{25}$.

Figure 2 compares the perturbation in $\Delta V_{yy}$ with the signal dynamics of common mode accelerations $\bar{a}_{c25y}$ and squared common mode accelerations $\bar{a}_{c25y}^2$ for a 2-day time window starting on 6 April 2016, 6:00 am. In order to highlight the perturbation, we bandpass filter all quantities to the frequency band 1–10 mHz. We signify bandpass-filtered quantities with a tilde in the following. We can clearly see that the perturbation in $\Delta \tilde{V}_{yy}$ is centered at 90° argument of latitude, which corresponds to the location along the orbit that is closest to the geographic North pole. In the same location along the orbit, we can also see strong dynamics in $\tilde{a}_{c25y}$, which show a similar, but not the same pattern as in $\Delta \tilde{V}_{yy}$. Also $\tilde{a}_{c25y}^2$ shows strong dynamics in the same location along the orbit and here the pattern matches almost perfectly the one in $\Delta \tilde{V}_{yy}$. We analyzed many more time windows, which always confirmed the almost perfect match between the patterns in $\Delta \tilde{V}_{yy}$ and $\tilde{a}_{c25y}^2$.

When estimating the parameters of the quadratic function, we bandpass-filtered Eq. (50) to the frequency band 1–10 mHz, which is the frequency range in which the perturbation is clearly visible as demonstrated in Fig. 2 and the gradiometer noise is expected to be low. The terms $-2\bar{b}_{d25y}/L_y$ and $-2k_{d25y}\bar{b}_{c25y}^2/L_y$ are much smaller than $n_{yy}$ after the bandpass filtering and are consequently neglected. The two remaining parameters $k_{d25y}$ and $cd25y \equiv -2k_{d25y}\bar{b}_{c25y}$
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Time series of estimated parameters $k_{d25y}$ (top panel) and $c_{d25y}$ (bottom panel). The dots indicate the limits of each time interval, for which the parameters were estimated</figcaption></figure>

were allowed to change linear in time and are estimated by minimization of $\tilde{n}_{yy}^2$. It should be noted that the effect of $c_{d25y}$ on $\Delta V_{yy}$ is similar to the effect of the element in the second row and fifth column of $M_{25}$.

We excluded the data from the estimation, when the satellite is not in drag-free mode either due to satellite anomalies or orbit maneuvers (GOCE Flight Control Team 2014). In practice, we estimate one pair of $k_{d25y}$ and $c_{d25y}$, i.e., four parameters due to the linear time dependency, for each time interval between two calibration maneuvers referred to as "satellite shakings" (Frommknecht et al. 2011) or any other satellite anomaly or orbit maneuver that causes the GOCE satellite to leave the drag-free mode. Satellite shakings were performed typically every 2 months.

# 4 Results

The estimated parameters $k_{d25y}$ and $c_{d25y}$ shown in Fig. 3 are both remarkably stable during mission lifetime. In this context, it should be noted that September 23, 2009, was the

last time the proof mass positions were adjusted in order to nullify quadratic factors (GOCE Flight Control Team 2014), which is before the period covered by the data shown in Fig. 3. A closer inspection reveals that $k_{d25y}$ seems to be approximately constant around -70.4 s$^2$/m before July 2010 and -67.6 s$^2$/m after December 2010. We can only speculate that this small change in mean might be related to the satellite anomaly that occurred in summer 2010 (GOCE Flight Control Team 2014), where it should be noted that the quadratic factors were determined several times afterward. They were found to be adequately nullified, so that the proof mass position was not adjusted. Such a change in mean is not visible for $c_{d25y}$, which appears to be constant around $55 \times 10^{-6}$.

As a representative example, we show in Fig. 4 the map of $\Delta \tilde{V}_{yy}$ for ascending orbital arcs in the period from March 16, 2012, to May 22, 2012, which is a period when the drag signal is already much larger compared to the beginning of the mission due to increased solar activity. Afterward, the drag signal increases further, also because of a sequence of orbit lowering maneuvers that started in November 2012.
**Fig. 4** Bandpass-filtered gravity gradient residuals $\Delta \tilde{V}_{yy}$ of present Level 1b gravity gradients (top panel), the bandpass-filtered quadratic factor correction for $\tilde{V}_{yy}$ (mid-panel), and the bandpass-filtered gravity gradient residuals $\Delta \tilde{V}_{yy}$ corrected for a quadratic factor (bottom panel). The green triangles indicate the location of the geomagnetic poles and the green line indicates the geomagnetic equator. In the top and mid-panel, the color scale is limited to ±30 mE, whereas peak values are larger than 400 mE

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

The top panel shows $\Delta \tilde{V}_{yy}$ as calculated from the EGG_NOM_1b data products. The perturbation in $\Delta \tilde{V}_{yy}$ occurs predominantly in the regions around the geomagnetic poles (green triangles) where its magnitude is at least 30 mE

as opposed to 5 mE at latitudes lower than 30°. For some of the orbits, the perturbation reaches peak values of more than ±400 mE, which is demonstrated in Fig. 5. In the example shown in Fig. 4, the perturbation is larger in the northern
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Bandpass-filtered gravity gradient residuals $\Delta \tilde{V}_{yy}$ of present Level 1b gravity gradients (blue line) and after correction for a quadratic factor (red line)</figcaption></figure>

hemisphere. It should be noted that there are also periods in which the opposite is the case.

The estimated quadratic function used for correcting the gravity gradients is illustrated in the mid-panel of Fig. 4 and matches almost perfectly the perturbation shown in the top panel. Even the small perturbing features at low latitudes are captured very well, which demonstrates that the correction improves the gravity gradient $\tilde{V}_{yy}$ globally, though the largest improvements are obviously in the regions around the geomagnetic poles.

The bottom panel of Fig. 4 shows $\Delta \tilde{V}_{yy}$ after removing the perturbation using the method introduced in Sect. 3. We clearly see that the perturbation in the regions around the geomagnetic poles is removed to a great extent. Due to reducing the range of the color scale from ±30 to ±5 mE, a much smaller perturbation that is parallel to the geomagnetic equator (green line) becomes visible. We suspect that it is due to small calibration errors that affect the angular rates $\omega_x$ and $\omega_z$ in Eq. (12). Since the attitude control of the satellite is performed by magnetic torquers, the angular rates correlate with the Earth's magnetic field. Therefore, small errors in, e.g., scale factors would result in perturbations that correlate with the magnetic field lines. Since the other gravity gradients are also affected by such small perturbations as shown in Siemes (2011), Figs. 26–30, we believe that it is possible to further improve the gravity gradients by revisiting the gravity gradiometer calibration (Rispens and Bouman 2011; Siemes et al. 2012), where the calibration model should be extended by quadratic factors. This is subject to ongoing investigations and is therefore not presented in this paper.

The reason for choosing linear time-dependent parameters is illustrated in Fig. 6. The top panel shows the bandpass-filtered gravity gradient residual $\Delta \tilde{V}_{yy}$ of present L1b data with respect to time and argument of latitude. We can clearly see that the residual is largest around 90° argument of latitude, which is near the north pole, and that there is a

significant variability over time. Most notably, the residual has the largest magnitude around April 23, 2012. The mid-panel shows $\Delta \tilde{V}_{yy}$ corrected for the effect of a quadratic factor where the parameters $k_{d25y}$ and $c_{d25y}$ are constant with respect to time. Though the residual is more than 10 times smaller than that without correction, we can clearly see remaining perturbing effects in March and May, whereas the residual is very small in the second half of April. The remaining perturbing effects are significantly reduced when the parameters $k_{d25y}$ and $c_{d25y}$ are allowed to change linear in time, which is illustrated in the bottom panel. Therefore, we favor the quadratic function where the parameters $k_{d25y}$ and $c_{d25y}$ are linearly time-dependent.

Figure 7 shows the amplitude spectral density (ASD) of the gravity gradient residual $\Delta V_{yy}$ before and after removing the perturbation (top and bottom panel, respectively) for different time intervals during mission lifetime. The ASD was calculated according to a variant of Welch's method (Welch 1967) where the mean was replaced by the median. As part of Welch's method, we used a Hanning window (Harris 1978), which is also known as Hann window. Before removing the perturbation, the ASD of $\Delta V_{yy}$ shows large errors in the frequency range 0.5–50 mHz. The magnitude of those errors shows a large variability over time, which we attribute to the variability over time of the non-gravitational acceleration $d_y$. At the beginning of the mission, the GOCE satellite is closer in time to the solar minimum and also flying at the initial, higher altitude. This means that the atmosphere density is still low and the drag environment is much more "quiet". As time advances, the solar maximum comes closer and the drag environment becomes less quiet. Starting in November 2012, the dynamics in the drag increase significantly not only due to the approaching solar maximum, but also due to the orbit lowering maneuvers since the atmosphere is denser at lower altitudes (Bruinsma et al. 2014).

It is striking to see that after removing the perturbation, the ASD of $\Delta V_{yy}$ is practically constant over time. In addition, it
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Bandpass-filtered gravity gradient residuals $\Delta \tilde{V}_{yy}$ of present Level 1b data (top panel) and the same after correction for a quadratic factor, where the parameters are constant (mid-panel) and linear time-dependent (bottom panel)</figcaption></figure>

also shows much lower errors than the one before removing the perturbation. For some of the time intervals toward the end of the GOCE mission, the improvement exceeds a factor of 10 in the frequency range 1–10 mHz. We conclude that removing the perturbation significantly improves the gravity gradient $V_{yy}$, which stresses the need for a more advanced calibration of the gradiometer that takes into account the quadratic factor $k_{d25y}$.

## 5 Discussion

The results show clearly that an unmodeled quadratic factor $k_{d25y}$ with a magnitude of approximately –70 s$^2$/m is responsible for the perturbation observed in the measured gravity gradient $V_{yy}$. Since we are only able to determine a differential quadratic factor, it is not possible to deduce the quadratic factors for the individual accelerometers. Never-

<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 7** Amplitude spectral density (ASD) of the gravity gradient residuals $\Delta V_{yy}$ of present Level 1b gravity gradients (top panel) and the same corrected for a quadratic factor (bottom panel)</figcaption></figure>

theless, it should be noted that the initial quadratic factors for the individual accelerometer reported in Frommknecht et al. (2011) have a much larger magnitude, reaching values up to -685.2 s²/m. Since the differential quadratic factor that we determined is smaller and because we have no evidence of perturbing effects due to other quadratic factors, we still consider the method for physically nullifying the quadratic factors to be very successful.

Smaller perturbations that correlate with Earth's magnetic field become visible in $\Delta V_{yy}$ after removing the large perturbing effect due to the quadratic factor $k_{d25y}$. Also the other gravity gradients are affected by such perturbations as shown in Siemes (2011) and we speculated in Sect. 3 that these are due to an imperfect gradiometer calibration. Since quadratic factors are not accounted for in the present calibration model, we recommend to extend the calibration model by quadratic factors, recalibrate the gradiometer data and reprocess the

gravity gradients contained in the EGG_NOM_1b data products.

An assessment of the impact on the gravity field model retrieval will reveal how the improvement described in this paper compares to earlier improvements of the gravity gradients described in Stummer et al. (2011). Even though this is not analyzed as part of this paper, we expect that the gravity field models and gravity gradient grids will benefit significantly from removing the perturbation. Regional studies that make direct use of the gravity gradients for geophysical modeling in the polar regions are expected to benefit most. For example, Bouman et al. (2014) demonstrated that GOCE gravity gradients combined with data of the GRACE mission are sensitive to ice mass changes of the West Antarctic ice sheet at basin scale, whereas a similar analysis for the Greenland ice sheet was hampered by the perturbation in
$V_{yy}$. The latter should be reassessed in view of the findings presented in this paper.

**Acknowledgements** The GOCE data were provided by the European Space Agency (earth.esa.int/goce). The GOCO03s model was downloaded from the Gravity Observation Combination (GOCO) project webpage www.goco.eu. Roger Haagmans is gratefully acknowledged for fruitful discussions and reviewing an early version of the manuscript.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

## References

Becker S, Brockmann JM, Schuh WD (2014) Mean dynamic topography estimates purely based on GOCE gravity field models and altimetry. Geophys Res Lett 41:2063–2069. doi:10.1002/2014GL059510

Bouman J, Fiorot S, Fuchs M, Gruber T, Schrama E, Tscherning C, Veichert s M, Visser P (2011) GOCE gravitational gradients along the orbit. J Geod 85:791–805. doi:10.1007/s00190-011-0464-0

Bouman J, Fuchs M, Ivins E, Van der Wal W, Schrama E, Visser P, Horwath M (2014) Antarctic outlet glacier mass change resolved at basin scale from satellite gravity gradiometry. Geophys Res Lett 41:5919–5926. doi:10.1002/2014GL060637

Bouman J, Ebbing J, Fuchs M, Sebera J, Lieb V, Szwillus W, Haagmans R, Novak P (2016) Satellite gravity gradient grids for geophysics. Sci Rep 6:21050. doi:10.1038/srep21050

Brockmann JM, Zehentner N, Höck E, Pail R, Loth I, Mayer-Gürr T, Schuh WD (2014) EGM TIM RL05: an independent geoid with centimeter accuracy purely based on the GOCE mission. Geophys Res Lett 41:8089–8099. doi:10.1002/2014GL061904

Bruinsma SL, Doornbos E, Bowman BR (2014) Validation of GOCE densities and evaluation of thermosphere models. Adv Space Res 54:576–585. doi:10.1016/j.asr.2014.04.008

Bruinsma S, Förste C, Abrikosov O, Lemoine JM, Marty JC, Mulet S, Rio MH, Bonvalot S (2014) ESA’s satellite-only gravity field model via the direct approach based on all GOCE data. Geophys Res Lett 41:7508–7514. doi:10.1002/2014GL062045

Drinkwater MR, Haagmans R, Muzi D, Popescu A, Floberghagen R, Kern M, Fehringer M (2007) The GOCE gravity mission: ESA’s first core Earth Explorer. In: Proceedings of 3rd International GOCE User Workshop, 6–8 November, Frascati, Italy, ESA SP-627, ISBN 92-9092-938-3, 1–8

Floberghagen R, Fehringer M, Lamarre D, Muzi D, Frommknecht B, Steiger C, Piñeiro J, da Costa A (2011) Mission design, operation and exploitation of the gravity field and steady-state ocean circulation explorer mission. J Geod 85:749–758. doi:10.1007/s00190-011-0498-3

Frommknecht B, Lamarre D, Meloni M, Bigazzi A, Floberghagen R (2011) GOCE level 1b data processing. J Geod 85:759–775. doi:10.1007/s00190-011-0497-4

Fuchs MJ, Bouman J, Broerse T, Visser P, Vermeersen B (2013) Observing coseismic gravity change from the Japan Tohoku-Oki 2011 earthquake with GOCE gravity gradiometry. J Geophys Res Solid Earth 118:1–10. doi:10.1002/jgrb.50381

GOCE Flight Control Team (2014) GOCE End-of-mission operations report. Technical report GO-RP-ESC-FS-6268, Issue 1, Revision 0. https://earth.esa.int/documents/10174/85857/2014-GOCE-Flight-Control-Team.pdf

The European GOCE Gravity Consortium EGG-C (2014) GOCE Level 2 Product Data Handbook. Technical report GO-MA-HPF-GS-0110, Issue 4, Revision 0. https://earth.esa.int/documents/10174/1650485/GOCE_Product_Data_Handbook-Level-2

Gruber Th, Gerlach C, Haagmans R (2012) Intercontinental height datum connection with GOCE and GPS-levelling data. J Geod Sci 2:270–280. doi:10.2478/v10156-012-0001-y

Harris F (1978) On the use of windows for harmonic analysis with the discrete Fourier transform. Proc IEEE 66(1):51–83

Hashemi Farahani H, Ditmar P, Klees R, Liu X, Zhao Q, Guo J (2013) The static gravity field model DGM-1S from GRACE and GOCE data: computation, validation and an analysis of GOCE missions added value. J Geod 87:843–867. doi:10.1007/s00190-013-0650-3

Hirt C (2014) GOCE’s view below the ice of Antarctica: satellite gravimetry confirms improvements in Bedmap2 bedrock knowledge. Geophys Res Lett 41:5021–5028. doi:10.1002/2014GL060636

Ince ES, Pagiatakis SD (2016) Effects of space weather on GOCE electrostatic gravity gradiometer measurements. J Geod 90:1389–1403. doi:10.1007/s00190-016-0931-8

Knudsen P, Bingham R, Andersen O, Rio MH (2011) A global mean dynamic topography and ocean circulation estimation using a preliminary GOCE gravity model. J Geod 85:861–879. doi:10.1007/s00190-011-0485-8

Lenoir B, Lévy A, Foulon B, Lamine B, Christophe B, Reynaud S (2011) Electrostatic accelerometer with bias rejection for gravitation and Solar System physics. Adv Space Res 48:1248–1257. doi:10.1016/j.asr.2011.06.005

Van der Meijde M, Pail R, Bingham R, Floberghagen R (2015) GOCE data, models, and applications: a review. Int J App Earth Observ Geoinform 35(A):4–15. doi:10.1016/j.jag.2013.10.001

Mayer-Gürr T, Kurtenbach E, Eicker A (2010) ITG-Grace2010. Institute of Geodesy and Geoinformation, Bonn. http://www.igg.uni-bonn.de/apmg/index.php?id=itg-grace2010. Accessed 13 Nov 2016

Mayer-Gürr T, GOCO consortium (2012) The new combined satellite only model GOCO03s. GGHS2012, Venice (Poster). http://www.goco.eu/

Pail R, Bruinsma S, Migliaccio F, Förste C, Goiginger H, Schuh WD, Höck E, Reguzzoni M, Brockmann JM, Abrikosov O, Veichert s M, Fecher T, Mayrhofer R, Krasbutter I, Sansö F, Tscherning CC (2011) First GOCE gravity field models derived by three different approaches. J Geod 85:819–843. doi:10.1007/s00190-011-0467-x

Petersheim N, Schlicht A, Stummer C, Yi W (2011) Impact of cross winds in polar regions on GOCE gradiometer data. In: Ouwehand L (eds) Proceedings of the 4th International GOCE User Workshop, 31 March–1 April 2011, Munich, Germany. ESA Publication SP-696, ISBN (Online) 978-92-9092-260-5, ISSN 1609-042X

Rispens SM, Bouman J (2011) External calibration of GOCE accelerations to improve derived gravitational gradients. J Geod Sci 1:114–126. doi:10.2478/v10156-010-0014-3

Rummel R, Yi W, Stummer C (2011) GOCE gravitational gradiometry. J Geod 85:777–790. doi:10.1007/s00190-011-0500-0

Schall J, Eicker A, Kusche J (2014) The ITG-Goce02 gravity field model from GOCE orbit and gradiometer data based on the short arc approach. J Geod 88:403–409. doi:10.1007/s00190-014-0691-2

Siemes C (2011) GOCE gradiometer calibration and L1b data processing. ESA Working Paper EWP-2384. earth.esa.int/goce

Siemes C, Haagmans R, Kern M, Plank G, Floberghagen R (2012) Monitoring GOCE gradiometer calibration parameters using accelerometer and star sensor data: methodology and first results. J Geod 86:629–645. doi:10.1007/s00190-012-0545-8

Stummer C, Fecher T, Pail R (2011) Alternative method for angular rate determination within the GOCE gradiometer processing. J Geod 85:585–596. doi:10.1007/s00190-011-0461-3
Stummer C, Siemes C, Pail R, Frommknecht B, Floberghagen R (2011)
Upgrade of the GOCE Level 1b gradiometer processor. Adv Space
Res 49:739-752. doi:10.1016/j.asr.2011.11.027

Welch PD (1967) The use of Fast Fourier Transform for the estimation of power spectra: a method based on time averaging over short modified periodograms. IEEE Trans Audio Electroacoust AU-15:70-73