An Improved Model of the Earth’s Static Gravity Field Solely Derived from Reprocessed GOCE Data

Abstract

After it was found that the gravity gradients observed by the Gravity field and steady-state Ocean Circulation Explorer (GOCE) satellite could be significantly improved by an advanced calibration, a reprocessing project for the entire mission data set was initiated by ESA and performed by the GOCE High-level processing facility (GOCE HPF). One part of the activity was delivering the gravity field solutions, where the improved level 1b and level 2 data serve as an input for global gravity field recovery. One well-established approach for the analysis of GOCE observations is the so-called time-wise approach. Basic characteristics of the GOCE time-wise solutions is that only GOCE observations are included to remain independent of any other gravity field observables and that emphasis is put on the stochastic modeling of the observations’ uncertainties. As a consequence, the time-wise solutions provide a GOCE-only model and a realistic uncertainty description of the model in terms of the full covariance matrix of the model coefficients. Within this contribution, we review the GOCE time-wise approach and discuss the impact of the improved data and modeling applied in the computation of the new GO_CONS_EGM_TIM_RL06 solution. The model reflects the Earth’s static gravity field as observed by the GOCE satellite during its operation. As nearly all global gravity field models, it is represented as a spherical harmonic expansion, with maximum degree 300. The characteristics of the model and the contributing data are presented, and the internal consistency is demonstrated. The updated solution nicely meets the official GOCE mission requirements with a global mean accuracy of about 2 cm in terms of geoid height and 0.6 mGal in terms of gravity anomalies at ESA’s target spatial resolution of 100 km. Compared to its RL05 predecessor, three kinds of improvements are shown, i.e., (1) the mean global accuracy increases by 10–25%, (2) a more realistic uncertainty description and (3) a local reduction of systematic errors in the order of centimeters.

1 Introduction

More than ten years ago, the dedicated gravity field mission Gravity field and steady-state Ocean Circulation Explorer (GOCE) was launched to its science orbit by the European Space Agency (ESA) on March 17, 2009. As the first Earth Explorer, the mission objective was to determine the static part of the Earth’s gravity field with unprecedented precision of 1–2 cm in terms of geoid heights and 1 mGal in terms of gravity anomalies (Floberghagen et al. 2011). With the new measurement principle of satellite gravitational gradiometry (for details on the measurement principle see for instance Rummel and Colombo 1985; Rummel 1986; Rummel et al. 2011), these accuracies can be reached globally for spatial resolutions down to 100 km. The mission was completed October 21, 2013, after running out of fuel and the spacecraft reentered the Earth’s atmosphere on November 11, 2013.

Already with the availability of 71 days of data, first gravity field solutions were computed from the level 1b data sets, which already indicated the errors and inaccuracies of the existing global gravity field models (e.g., Pail et al. 2011). For gravity field recovery, the essential data sets are the precise science orbits (Bock et al. 2007, 2011, 2014), the star camera observations which determine the orientation of the spacecraft and the accurate components of the gravity gradient tensor, which are by design $V_{xx}$, $V_{yy}$, $V_{zz}$ and $V_{xz}$ (see Cesare and Catastini 2008; Rummel et al. 2011; Stummer et al. 2011, 2012; Stummer 2013). First models based on these data (slightly more than one repeat cycle) were published in 2011 as the so-called RL01 generation (see Pail et al. 2011). The level 1b and level 2 data were generated and prepared by the GOCE high-level processing facility (GOCE-HPF, see EGG-C 2010). Especially for the gravity gradients, a complex calibration procedure of the accelerometers and a removal of the angular rates (angular rate reconstruction) are required to derive gravity gradients at a useful error level for gravity field recovery (Cesare and Catastini 2008; Stummer et al. 2011; Siemes et al. 2012, 2019b; Siemes 2018b). Typical approaches for gravity field recovery from GOCE start at the level of the calibrated level 1b gravity gradients.

Three official approaches are pursued by GOCE-HPF, and five releases have been published so far. On the one hand, solutions of RL01 to RL05 mainly differ in the amount of data used (see later Fig. 3a for an overview). This is roughly the same for all approaches. In addition to that, there were improvements in the processing from RL03 to RL04 and the calibration of the level 1b data was improved (see Stummer et al. 2011, 2012). The RL05 versions were generated with the entire mission data and published in 2014. They were the first versions which included the entire GOCE mission data set.

The three approaches pursued by the official ESA GOCE-HPF processing are (1) the space-wise approach, (2) the direct-approach and (3) the time-wise approach (Pail et al. 2011). Within the following, the main characteristics of the approaches are shortly summarized. More specific details can be found in the cited references.

Primary goal of the space-wise approach is to compute a global GOCE-only model based on a multi-step least squares collocation approach (Reguzzoni 2003; Migliaccio et al. 2011; Reguzzoni and Tselfes 2009). The approach uses the kinematic precise science orbits to obtain the long wavelength signal (SST-hl) using an adopted version of the energy balance approach (e.g., Visser et al. 2003). The potential is predicted at a mean satellite altitude grid and converted to a spherical harmonic expansion yielding a SST-hl solution. This is used to reduce the SGG observations, which are filtered along the orbit using a Wiener filter. From these filtered SGG observations, locally regularized patches are predicted and composed to a global grid (Gatti et al. 2016). This global grids can again be converted to spherical harmonics with a covariance matrix determined by Monte Carlo samples (Migliaccio et al. 2008). Key feature of the space-wise models and grids is the local regularization, which is tailored to the local signal resulting in an improved local signal to noise ratio.

The direct approach uses the gravity gradients in the gradiometer reference frame (GRF) and sets up the observation equations for the second derivatives of the gravity potential directly in the location where the observations were captured. To account for the temporal correlations in the gravity gradient observations, a band-pass filter is used (Pail et al. 2011; Bruinsma et al. 2013, 2014). This filter takes out all frequencies outside of the measurement band. Inside the measurement band, the noise characteristics are assumed as white. Instead of using GOCE SST-hl to increase the accuracy of the long wavelengths, GRACE and LAGEOS normal equations are used since RL03. Thus, the models computed with the direct approach actually are satellite-only combination models.

In contrast to that, as the space-wise models, the time-wise approach uses only GOCE observations to derive an independent model of the Earth’s gravity field as it was measured by the GOCE satellite. Using the short-arc integral equations approach (since RL04), normal equations in the spherical harmonic domain are determined from the precise science kinematic orbits (Mayer-Gürr et al. 2005; Mayer-Gürr 2006). Using empirical auto-covariance functions and variance component estimation in addition to the covariance matrix of the kinematic positions, a refined data adaptive stochastic model for the kinematic orbit positions is used. This results in a more realistic covariance matrix of the resulting SST-hl solution, which later controls the relative weighting when combined with the normal equations assembled from the gravity gradients. As for the direct approach, the observation equations for the gravity gradient observations are set up in the GRF and handled as a time series along the orbit. To model the highly correlated noise, data adaptive AutoRegressive Moving Average (ARMA) and AutoRegressive (AR) models are estimated from a residual time series. For different periods (i.e., segments) of the time series, tailored stochastic models are designed, individually for each component of the gravity gradients (Schuh 1996; Pail and Plank 2002; Pail et al. 2011; Brockmann et al. 2014). Further aspects of the time-wise model will be discussed in the manuscript.

In addition to the three official approaches followed by GOCE-HPF, various other approaches were investigated. Several methods which are applied for the processing of SST-hl from kinematic orbits, especially for GOCE, were inter-compared by Baur et al. (2012). Their basic conclusion is that, except for the energy balance approach where the three-dimensional observable is reduced to a scalar value, all approaches lead to a comparable performance. But, depending on the stochastic modeling of the observations, the accuracy estimates are of different quality, which is essential for a combination with the gravity gradients.

Several approaches exist which include GOCE gravity gradient observations for global and regional gravity field recovery. We restrict ourselves to global satellite-only models here. On the one hand, there are approaches following the rotational invariants, which have the theoretical advantage that the observation equations are independent of the satellites’ rotation, and the observation equations do not have to be rotated to the GRF (Baur 2007; Baur et al. 2008). This is an advantage, as the satellites orientation is an observed quantity as well and it is not known without errors. But, as the higher order invariants are nonlinear functions of all tensor elements—including the inaccurate ones—for real data analysis linearization and substitution with model information are required, which again requires the satellites’ orientation for synthesization. Models like IGGT_R1 (Lu et al. 2018) indicate that a strong bias toward the used ’fill-in’ models is introduced.

Several authors applied similar approaches compared to the time-wise and the direct approach. Basically, they set up the least-squares observation equations in the GRF, but the way how the correlations are modeled and relative weights are estimated differs (e.g., Yi et al. 2013; Schall et al. 2014; Wu 2016; Xu et al. 2017). For instance, Schall et al. (2014) applied a short-arc approach, where the gravity gradient data were analyzed in short arcs. For the short arcs, an empirical covariance matrix could be derived from the spectrum. Additional empirical parameters are used to model the correlations between the individual short arcs.

Some other studies, similar to the direct approach, directly estimate a combined GRACE and GOCE model. One independent model is for instance GOCE DGM-1S (Farahani et al. 2013), where a lot of effort is spent on the stochastic modeling of the involved data sets to optimally combine the observations to derive the global gravity field model. In contrast to this, models like those of the GOCO series (GOCO 2020; Pail et al. 2010; Kvas et al. 2019b) make use of published individual satellite-only models, which are combined on the level of normal equations. E.g., the GOCO05s (Mayer-Gürr et al. 2015) combines the time-wise RL05 model with the ITSG-Grace2014s GRACE-only model and SLR normal equations.

Since the RL05 models, which already included the entire mission data set, a reprocessing campaign was performed by the GOCE-HPF to derive a new level 1b data set. The reprocessing covered all data sets which are required for gravity field recovery, i.e., the precise science orbit, the satellites orientation and the gravity gradients. For gravity field recovery, the most visible update is related to the gravity gradient calibration. Siemes (2018b) identified a correlation of the cross-track gravity gradient residual with the squared common mode accelerations and developed an approach to reduce the correlation, re-calibrating the gravity gradient by an additional quadratic factor. This study led to the reprocessing campaign, where the calibration procedure was extended by a quadratic factor in the standard processing. The full applied calibration is described in detail in Siemes et al. (2019b, 2019a). The result of this reprocessing is a new gravity gradient data set (EGG_NOM_1b) in version 202 (ESA GOCE-ODS 2020).

Various studies indicated, e.g., Jäggi et al. (2015), that the kinematic GOCE orbits were degraded depending on ionospheric activities. Within GOCE-only gravity field solutions, it was visible in terms of a long wavelength error correlating with the magnetic equator. Within the reprocessing of the kinematic orbits, affected GPS phase observations were carefully detected and down-weighted to significantly reduce the systematic error.

Within this manuscript, the processing of the time-wise RL06 solution is presented (GO_CONS_EGM_TIM_RL06, Brockmann et al. 2019). Based on the reprocessed level 1b data set and an improved processing, the entire mission data set is analyzed to obtain the more accurate RL06 solution. The manuscript is organized as follows. In Sect. 2, the theoretical background on the time-wise solution is summarized. It is discussed how the different observation types are processed and how they are combined to a final solution. Section 3 concentrates on the used data and its characteristics. The different data sets are analyzed individually to demonstrate the characteristics of the data as well as the consistency of the final solution. In Sect. 4, the combined model and its performance compared to the RL05 solution are shown. This paper ends with a short summary and some final conclusions on the resulting GO_CONS_EGM_TIM_RL06 solution in Sect. 5.

2 The Time-Wise Approach: An Overview

An essential feature of the GOCE time-wise global gravity field solutions is that they only depend on GOCE data, avoiding the usage of any other gravity field observations. Thus, they are well suited for further combination and reflect the Earth gravity field as observed by GOCE.

As the data of a single mission is processed starting at observation level, consistency is guaranteed and rigorous error modeling from the observations to the final model estimate is realized. This results in the second request to the time wise-models, that is to provide a realistic error description of the model in terms of a full covariance matrix. This section discusses the methodological basics applied in the processing of the new time-wise RL06 GOCE-only model (GO_CONS_EGM_TIM_RL06) and provides references to some deeper insights.

For a GOCE-only model, to achieve an acceptable performance of the entire spherical harmonic spectrum, a combined gravity field estimated from the kinematic orbits and the gravity gradients of good quality is essential. In addition, due to the polar gap (e.g., Sneeuw and van Gelderen 1997; Rudolph et al. 2002), some mathematical regularization has to be applied to constrain the data gaps and to improve the signal to noise ratio for the higher degrees.

2.1 Spherical Harmonic Representation

Conventionally, global gravity field models are represented by a truncated series of spherical harmonics. Following (e.g., Hofmann-Wellenhof and Moritz 2005), the finite series for the Earth’s gravitational potential $V\left( r,\theta ,\lambda \right)$ reads

$$\begin{aligned} V\left( r,\theta ,\lambda \right) = \frac{\text {GM}}{a}\sum _{l=0}^{l_{\max }} \left( \frac{a}{r}\right) ^{l+1}\sum _{m=0}^{l}\left( c_{lm}\cos \left( m\lambda \right) +s_{lm}\sin \left( m\lambda \right) \right) P_{lm}\left( cos \theta \right) , \end{aligned}$$

(1)

where $l_{\max }$ denotes the maximum degree and r, $\theta$ and $\lambda$ the coordinates of the evaluation point. GM and a are the gravitational constant of the Earth and the equatorial radius. $P_{lm}\left( \cdot \right)$ denote the fully normalized associated Legendre base functions.

Given the set of $u = (l_{\max }+1)^2$ fully normalized spherical harmonic parameters $\left\{ c_{lm},\,s_{lm}\right\}$ (fully normalized Stokes coefficients) for $l \in 0\dots l_{\max }, m\in 0\dots l$ ($m\in 1\dots l$ for sine coefficients), various functions of the potential can be evaluated (see Barthelmes 2013). Gravity field recovery from GOCE is an inverse problem. Given an observation $\ell _i$ at some location, it is expressed as some function of the potential

$$\begin{aligned} \ell _i = f\left( V\left( \cdot \right) \right) =f\left( \left\{ c_{lm},\,s_{lm}\right\} \right) . \end{aligned}$$

(2)

The function $f\left( \cdot \right)$ on the one hand depends on the type of the observations, which for GOCE are either the kinematic orbits position in the case of SST-hl or a component of the gravitational tensor in the case of the gravity gradient analysis.

Having a suitable function $f\left( \cdot \right)$ relating the observations to the spherical harmonic parameters, they can be estimated in a least-squares adjustment, minimizing the squared sum of residuals and introducing the inverse covariance matrix of the observations $\varvec{\Sigma }\left\{ \varvec{\ell }\right\} ^{-1}$ as a metric.

2.2 Processing of SST-hl Observations

The SST-hl observations are processed and made available by the Institute of Geodesy of the Graz University of Technology. Based on the kinematic orbits as input observations, normal equations are assembled based on the short arc integral equation approach (cf. Mayer-Gürr et al. 2005; Mayer-Gürr 2006) in the same way as done for ITSG GRACE analysis (Kvas et al. 2019a). Input for the generation of the GO_CONS_EGM_TIM_RL06 model is a set of normal equations

$$\begin{aligned} \mathbf {N}_{\mathrm{sst}}\mathbf {x}&= \mathbf {n}_{\mathrm{sst}} \end{aligned}$$

(3)

$$\begin{aligned} \lambda _{c,s}&:= \varvec{\ell }_{\mathrm{sst}}^{T}\mathbf {P}_{\mathrm{sst}} \varvec{\ell }_{\mathrm{sst}}, \end{aligned}$$

(4)

which is assembled from the official reprocessed kinematic orbit product.

Arcs of a length of 20 min are used. To account for the relative weighting, the given $3\times 3$ epoch covariance matrix is refined by arc dependent empirically estimated covariance functions ($\mathbf {P}_{\mathrm{sst}}$ in (4)). Following the GRACE processing strategy (Kvas et al. 2019a), temporal variations are reduced before the assembly of normal equations, resulting in a reference epoch for the model as 2010-01-01 (MJD 55197). The applied correction models are summarized by Table 1 (Mayer-Guerr 2019). Normal equations up to degree and order 150 serve as input for the RL06 model.

Table 1 Models applied to account for temporal variations of the gravity field

2.3 Processing of the Observed Gravity Gradients

The level 1b processing with all its steps (cf. Siemes 2018a) derives the gravity gradients from the measurements of the three-axis gradiometer which is realized by three pairs of three-axis accelerometers. The individual gravity gradients are the entries of the symmetric Marussi-tensor in the gradiometer reference frame (GRF) at a location $\mathbf {p}$ (e.g., Rummel et al. 2011). The gradiometer reference frame is a local Cartesian coordinate system, which is realized by the three axes of the gradiometer arms. The x-axis is roughly in the flight direction, the y-axis orthogonal to the orbital plane, and the z-axis is pointing radial downwards (cf. Rummel et al. 2011). The tensor represents the second derivatives of the gravitational potential in the GRF. Consequently, it is required to relate (1) in the earth fixed reference frame to the observations in the GRF. For this, in addition to the Earth rotation, the star camera observations are required (EGG_IAQ product), describing the orientation of the satellite in the inertial frame.

2.3.1 Observation Equations for the Gravity Gradients

As it is well known (e.g., Rummel et al. 2011; Pail et al. 2011), two of the six gravity gradients are by design of lower quality. The $V_{xy}$ and $V_{yz}$ gravity gradients have a noise level, which is about 100 times higher than the accurately measured $V_{xx}$, $V_{xz}$, $V_{yy}$ and $V_{zz}$ gradients. To indicate this different error level, Fig. 1 shows the power spectral density (PSD) of least-squares residuals for a representative period (2013, August and September) without outliers and gross errors. All gravity gradient components are shown in a logarithmic scale.

Fig. 1

Power spectral density of the gradiometer noise estimate for the different gravity gradient components. It is computed by Welch’s method (e.g., Oppenheim et al. 1999, Sect. 10.6) from residuals with respect to the EGM_TIM_RL06 solution from the time series August and September 2013

Within the time-wise approach, the observations of an individual gradient component $c\in \{xx,\,yy,\,zz,\,xz\}$ are treated as a time-series along the orbit. Without loss of generality, we can assume that we can decompose the entire data set into S segments $s\in \{1,\ldots ,S\}$, where all observations in a segment are gapless and equidistantly sampled. For GOCE the sampling rate is 1 Hz.

As for the case of GOCE only four components of the tensor are usable for gravity field recovery, it is important to not transform the tensor of observations to the earth-fixed frame where (1) is defined, but to transform the observation equations of the second derivatives of (1) to the GRF. In case the tensor was rotated, a mixing of precise and inaccurate gradients takes place, which results in an imperfect signal-to-noise ratio.

The transformed observation equations are calculated as follows. For a specific measurement location $\mathbf {p}_i$ of a segment s, the six second partial derivatives of (1) according to $\lambda$ and $\theta$ can be computed. The observation equations can be set up as a $6\times u$ design matrix, which can then be rotated from the earth-fixed reference frame to the GRF (using the Earth rotation parameters and the satellite orientation as derived from the star trackers) applying tensor rotation of the design matrix. This procedure is discussed in more detail in Brockmann et al. (2010) and especially in Hausleitner (1995). The observation equation for a single gravity gradient component c can be extracted from the associated row of the rotated time-series. Consequently, repeating the procedure for each measurement location, we can arrange least squares observation equations individually for the gradient components, i.e.

$$\begin{aligned} \varvec{\ell }_{c,s} + \mathbf {v}_{c,s}&= \mathbf {A}_{c,s} \mathbf {x} \end{aligned}$$

(5)

where the vector $\mathbf {x}$ contains the unknown Stokes coefficients in a specific ordering. Here, we already assume that the measured gravity gradients are reduced by time-variable gravity signal derived from models (as for SST-hl, cf. Table 1). The least-squares design matrix $\mathbf {A}_{c,s},$ and the observation vector $\varvec{\ell }_{c,s}$ contains only observation of the single gradient component c in the segment s. Consequently the observations in $\varvec{\ell }_{c,s}$ are equidistant and sorted in time.

2.3.2 Stochastic Model for the Gravity Gradients

A key feature of the time-wise approach is the modeling of the stochastic characteristics of the gravity gradients. As the noise is highly correlated (cf. Fig. 1), a stochastic model is important to derive best possible estimates for the unknown gravity field parameters and reasonable estimates for its accuracy. For that reason, the observations $\varvec{\ell }_{c,s}$ are treated as a time series along the orbit. Due to its efficiency, the correlations in that time series are modeled with the filter approach (Schuh 2003; Siemes 2008; Schuh and Brockmann 2018). Whereas in RL01 to RL05 AutoRegressive Moving Average (ARMA) filters were used to model the correlations (e.g., Siemes 2008; Krasbutter et al. 2011, 2014), in RL06 a robustified procedure to estimate AR-processes from the residuals is introduced (cf. Schubert et al. 2019).

Starting with residuals of (5)

$$\begin{aligned} \mathbf {v}_{c,s}&= \mathbf {A}_{c,s} \mathbf {x}_{\mathrm{REF}} - \varvec{\ell }_{c,s} \end{aligned}$$

(6)

for a single component and a single segment with respect to some reference model $\mathbf {x}_{\mathrm{REF}}$ which is either initially the former release (EGM_TIM_RL05 in this case) or the output from a former iteration of adjusting the coefficients of an AR process in a robustified procedure.

The robustification approach implements a series of statistical tests that cover typical systematic errors in highly correlated observation time series in order to detect outliers or non-stationarities, summarized as suspicious data. We incorporate different indicators to detect suspicious data or changes in the underlying process characteristics, i.e., changes in the mean value, the variance or in the signs of the residuals. The identified flags for these effects are used to make the process estimation procedure robust and to reject detected suspicious data for gravity field recovery and the decorrelation filter estimation.

The time series of the residuals $\mathbf {v}_{c,s}$ are seen as an individual realization of a stochastic process ${\varvec{{\mathcal {V}}}}$ given as a sequential signal ${\mathcal {V}}_i$. The process is characterized by an autoregressive model of order p

$$\begin{aligned} {\mathcal {V}}_i+{\mathcal {R}}_i&= \sum _{j=1}^p \alpha _j {\mathcal {V}}_{i-j} \ . \end{aligned}$$

(7)

where ${\mathcal {R}}_i$ represents a random noise contribution. The AR-process parameters are collected component-wise in a coefficient vector $\varvec{\alpha }_{c,s}$ for each segment s. As it is assumed that the residuals $\mathbf {r}_{c,s}$ are realizations of the white noise process ${\varvec{{\mathcal {R}}}}_{c,s}$, $\mathbf {r}_{c,s}$ can be statistically tested for white noise characteristics.

For the approach, we apply a series of tests, namely, a test for

• 1. Simple outliers, i.e., gross errors,

• 2. The deviation from zero mean within a certain period of time,

• 3. A change in the variance over a period,

• 4. The imbalance of signs and

• 5. A suspicious occurrence of sign changes.

Tests 2–5 are applied in moving windows over the time series leading to patches of outliers. The window size ${n_{T_l}}$ is chosen as a fraction of orbit period (e.g., the 30th part) and deployed with three different window sizes ($l=1,2,3$) for each test.

For the processing of RL06, the hypothesis tests have experienced a modification compared to Schubert et al. (2019), by introducing the more heavy-tailed Student’s t-distribution with statistics derived from robust scale estimators for data with outliers. In detail, the degrees of freedom of Student’s t-distribution are put in relation to the ratio of mean error $\sigma$ and average error $\vartheta$ (cf. (27) and Table 5 in Schuh and Korte 2017). Thus, the hypothesis tests in Schubert et al. (2019) are changed to consistently follow Student’s t-distribution or Fisher distribution, respectively, i.e.

$$\begin{aligned} {\mathcal {T}}_{\mathrm{outlier}}&= \frac{{\mathcal {R}}_i}{\widetilde{s}_{\mathrm{MAD}}} \sim {t(\nu )} \ , \end{aligned}$$

(8)

$$\begin{aligned} {\mathcal {T}}_{\mathrm{mean}}&= \frac{1}{\sqrt{n_{T_l}} \widetilde{s}_{\mathrm{MAD}} } \sum _{j=1}^{n_{T_l}} {\mathcal {R}}_j \ \sim \ t(\nu ) \ , \end{aligned}$$

(9)

$$\begin{aligned} {\mathcal {T}}_{\mathrm{var}}&= \frac{\sum _{j=1}^{n_{T_l}} {\mathcal {R}}_i^2}{n_{T_l} \widetilde{s}_{\mathrm{MAD}}^{2}} \ \sim \ F({n_{T_l}},\nu ), \end{aligned}$$

(10)

where $t(\cdot )$ and $F(\cdot ,\cdot )$ stand for the Student’s t-distribution and F-distribution. $s_{\mathrm{MAD}}$ represents the robust variance estimator $\widetilde{s}_{\mathrm{MAD}}=1.4826\,\mathrm{MAD}$ (e.g., Schuh and Korte 2017, Tab 2.) MAD (Median Absolute Deviation) estimator (e.g., Huber 1981, p. 107).

Finally, for estimating the AR-coefficients $\alpha _j$ we take use of a $k\sigma$-rejection-estimator in an iteratively re-weighted least squares (IRLS) estimator in the sense of Kleiner et al. (1979) which enables the treatment of single outliers and outlier patches. It is embraced by additional iterations realizing the adjustment of the degrees of freedom. The permanent recomputation of the variance estimators as the scale-step of the well-known IRLS-scheme is essential for the unbiasedness of the scale parameter.

For the adjustment of AR(p)-processes in RL06, the process order is chosen as $p=800$. The high order is required to model the characteristics below the measurement band of 5ṁHz. To further account for strong drifts, which mainly occur in the $V_{yy}$ observations, a difference filter is applied in advance of the process estimation and integrated to the final set of process coefficients.

For additional details see Schubert et al. (2019). Output of the stochastic model tuning are the process coefficients $\varvec{\alpha }_{c,s}$ and binary suspicious data flags $\mathbf {f}_{c,s}$ (use, do not use) individually for all segments s and for each of the usable gravity gradient components c.

2.3.3 Assembly of Normal Equations Per Segment and Gradient Component

The computational most demanding step is the assembly of the normal equations for the gravity gradients. Due to the number of observations (i.e., about $110\times 10^{6}$ epochs with four gradient components each), the number of unknown parameters ($\approx 90,000$ for RL06) and the high correlations in the data, high-performance computing is required to decorrelate the observations and to assemble the system of normal equations (Pail and Plank 2002; Plank 2004; Brockmann 2014).

The challenges are slightly relaxed, as discussed above, the partitioning into segments s is performed and observation equations are set up for the different components c individually. It is assumed that the data after a gap, so with the start of a new segment, are uncorrelated from the previous one, i.e.

$$\begin{aligned} \varvec{\Sigma }\left\{ \varvec{\ell }_{c,s},\,\varvec{\ell }_{c,s-1}\right\} =\mathbf {0}. \end{aligned}$$

(11)

Furthermore, as indicated by empirical cross-correlation estimates, it is assumed that the correlations between different gravity gradients in the measurement band are small and can be neglected, consequently

$$\begin{aligned} \varvec{\Sigma }\left\{ \varvec{\ell }_{c_1,s},\,\varvec{\ell }_{c_2,s}\right\} =\mathbf {0} \quad \text { for } c_1\ne c_2. \end{aligned}$$

(12)

This reduces the numerical complexity for the decorrelation, but there are segments with up to $6.4\times 10^{6}$ epochs with still 90,000 parameters for a single segment and their length is limited by bimonthly calibration maneuvers. This would result in a fully populated design matrix which has memory requirements of about 4 TB. Special techniques in a high-performance computing environment are implemented to efficiently solve the required steps. The applied approaches are discussed in detail in the literature (Brockmann 2014, Sect. 6.3). This step comprises i) the decorrelation with the component and segment dependent AR process from Sect. 2.3.2 to derive fully decorrelated observation equations

$$\begin{aligned} \bar{\varvec{\ell }}_{c,s} + \bar{\mathbf {v}}_{c,s}&= \bar{\mathbf {A}}_{c,s} \mathbf {x}, \qquad \text { with } \varvec{\Sigma }\left\{ \varvec{\ell }\right\} = \mathbf {I}, \end{aligned}$$

(13)

followed by ii) the computation of the least squares normal equations again per component and segment

$$\begin{aligned} \mathbf {N}_{c,s}&:= \bar{\mathbf {A}}_{c,s}^{T} \bar{\mathbf {A}}_{c,s}^{T} \end{aligned}$$

(14)

$$\begin{aligned} \mathbf {n}_{c,s}&:= \bar{\mathbf {A}}_{c,s}^{T} \bar{\varvec{\ell }}_{c,s} \end{aligned}$$

(15)

$$\begin{aligned} \lambda _{c,s}&:= \bar{\varvec{\ell }}_{c,s}^{T}\bar{\varvec{\ell }}_{c,s}. \end{aligned}$$

(16)

which form the least squares normal equation system for component c and segment s

$$\begin{aligned} \mathbf {N}_{c,s} \mathbf {x}_{c,s}&= \mathbf {n}_{c,s}. \end{aligned}$$

(17)

To avoid a recomputation of the computationally expensive steps within the relative iterative weighting and combination procedure (cf. Sect. 2.5) and having individual sub-solutions available for the internal analysis (cf. Sect. 3.2), all $\mathbf {N}_{c,s}$, $\mathbf {n}_{c,s}$ and $\lambda _{c,s}$ are stored on disk.

The decorrelation process with digital filters can be represented as a matrix multiplication,

$$\begin{aligned} \bar{\mathbf {A}}_{c,s}&:=\mathbf {F}_{c,s}\mathbf {A}_{c,s} \end{aligned}$$

(18)

$$\begin{aligned} \bar{\varvec{\ell }}_{c,s}&:=\mathbf {F}_{c,s}\varvec{\ell }_{c,s} \end{aligned}$$

(19)

which can be efficiently computed (cf. Brockmann 2014, Sect. 6.3.3) using the fact that $\mathbf {F}_{c,s}$, as it results from a causal time-invariant filter, is a band matrix of order p with Toeplitz structure.

As this requires equidistant time-series in $\varvec{\ell }_{c,s}$, the suspicious data flags are applied after the decorrelation by setting corresponding rows of $\bar{\mathbf {A}}_{c,s}$ to zero, using the information provided by $\mathbf {f}_{c,s}$. Using a massive parallel implementation and a two-dimensional block cyclic representation of all matrices from (13), the computation of (14)-(16) is straightforward using numerical libraries like PBLAS and ScaLAPACK (Blackford et al. 1997; Brockmann 2014).

Compared to the assembly of (14), the solution of (17) is fast¹. Consequently, the implemented software package tries to solve (17) by default after its assembly, which is of course only possible if (17) is positive definite which is the case for longer segments (about longer than seven days). These segment-wise and component-wise solutions will help to analyze the internal consistency and the quality of the data (cf. Sect. 3.2).

Due to (11) and (12) component-wise full normal equations

$$\begin{aligned} \underbrace{\sum _{s=1}^{S} w_{c,s}\mathbf {N}_{c,s}}_{:= \mathbf {N}_{c}} \mathbf {x}_c&= \underbrace{\sum _{s=1}^{S} w_{c,s}\mathbf {n}_{c,s}}_{:= \mathbf {n}_{c}} \end{aligned}$$

(20)

$$\begin{aligned} \mathbf {N}_{c}\mathbf {x}_c&= \mathbf {n}_{c} \end{aligned}$$

(21)

and the full gravity gradient normal equations

$$\begin{aligned} \underbrace{\sum _{c} \mathbf {N}_{c}}_{:= \mathbf {N}_{\mathrm{sgg}}} \mathbf {x}_{\mathrm{sgg}}&= \underbrace{\sum _{c} \mathbf {n}_{c}}_{:= \mathbf {n}_{\mathrm{sgg}}} \quad c\in \{xx,yy,zz,xz\} \end{aligned}$$

(22)

$$\begin{aligned} \mathbf {N}_{\mathrm{sgg}}\mathbf {x}_{\mathrm{sgg}}&= \mathbf {n}_{\mathrm{sgg}} \end{aligned}$$

(23)

can be computed and solved as well. $w_{c,s}$ is a segment and component depending weighting factor, which is estimated from the data, cf. Sect. 2.5. As the filters are estimated reflecting a complete decorrelation, $w_{c,s}\approx 1.0$ has to hold so that the estimation of the weights is another quality control.

2.4 Prior Information to Stabilize the GOCE-Only Solution

Due to the polar gap and in general the downward continuation, the inverse problem of determining the gravity field from GOCE observations is ill-posed (e.g., Sneeuw and van Gelderen 1997; Rudolph et al. 2002). Due to the gaps at the poles, especially the near zonal coefficients are highly correlated and the approximation within the polar areas is poor. Meanwhile, it is state of the art to include some prior information to constrain the field in the polar areas and to damp the amplified noise in the higher degrees to decorrelate the coefficients (e.g., Koch and Kusche 2002; Ditmar et al. 2003).

Various kinds of prior information are possible, e.g., additional data sources from GRACE (e.g., Bruinsma et al. 2014) or terrestrial campaigns (e.g., Zingerle et al. 2019b; Lu et al. 2020). To remain independent of other gravity field data, the time-wise modeling applies simple mathematical smoothness conditions to force the full signal toward zero using regularization. To account for the different characteristics, different approaches are applied for i) the polar gap and ii) the high degree coefficients.

Polar gap Whereas a strong regularization of the near zonal coefficients influenced by the polar gap (cf. van Gelderen and Koop 1997; Sneeuw and van Gelderen 1997) toward zero was applied in the former releases of the time-wise models, RL06 uses pseudo zero gravity anomaly observations. For that purpose, gravity anomaly observation equations in the polar caps on a highly oversampled regular grid are set up for the empirically determined starting degree 11 to maximum spherical harmonic degree 300 and full normal equations with a zero right hand sides are computed numerically. These normal equations are generated individually for the south pole (sp) area (below $-83^{\circ }$ with $0.5^{\circ }$ sampling in longitude and latitude direction) and the north pole (np) area (above $83^{\circ }$ with $0.5^{\circ }$ sampling in longitude and latitude direction). The generated observations are assumed to be uncorrelated and have the same accuracy.

This results in the fully populated least squares normal equations

$$\begin{aligned} \mathbf {N}_{a}&:= \bar{\mathbf {A}}_{a}^{T} \bar{\mathbf {A}}_{a}^{T}, \quad \text {for } a \in \left\{ \text {np, sp}\right\} \end{aligned}$$

(24)

$$\begin{aligned} \mathbf {n}_{a}&:= \bar{\mathbf {A}}_{a}^{T} \mathbf {0} = \mathbf {0} \end{aligned}$$

(25)

$$\begin{aligned} \lambda _{a}&:= 0. \end{aligned}$$

(26)

Note that they are assembled in a different numbering, as the coefficients below degree 11 are not contained.

High degrees To improve the signal-to-noise ratio for the higher degrees, a regularization of the full signal toward zero is applied. For this purpose, a diagonal regularization matrix is set up for all spherical harmonic coefficients above degree 200. The choice is motivated by a contribution analysis as discussed in detail in the upcoming Sect. 4.2. Stochastic pseudo observations stating that the coefficients are zero are introduced, with a variance which is derived from the average coefficient signal from Kaula’s rule of thumb (Kaula 1966), again to remain independent of any gravity model prior information. The observation equations lead to a diagonal normal equation matrix and again a zero right hand side for the high degree regularization (hd), whose entries depend only on the degree $l_i$ of the i-th coefficient

$$\begin{aligned} \mathbf {N}_{\mathrm{HD}}&, \quad \text {with } \mathbf {N}_{\mathrm{HD}}(i,i) = 10^{10} l_i^4 \quad l_i=201,...,300 \end{aligned}$$

(27)

$$\begin{aligned} \mathbf {n}_{\mathrm{HD}}&= \mathbf {0} \end{aligned}$$

(28)

$$\begin{aligned} \lambda _{\mathrm{HD}}&= 0. \end{aligned}$$

(29)

The scaling of $\mathbf {N}_{\mathrm{HD}}$ is derived with variance component estimation later.

2.5 The Combination Procedure

Due to the assumption that all individual normal equations as set up in Sects. 2.2–2.4 are independent, they can be easily combined by simply adding the normal matrices, right hand sides and $\lambda _i$. But, one has to account for the fact that the parameter spaces differ, as not all normal equations are set up for all spherical harmonic coefficients. Consequently, the numbering scheme differs and reordering strategies have to be applied.

As the gradiometry normal equations are set up for the entire parameter space, it can be assumed that the combined normal equations are set up in the numbering scheme the SGG normal equations were assembled in. Defining a reordering operator $\mathrm {reo}\left( \cdot \right)$ which performs the necessary row (and column) permutations, the combined normal equation system reads

$$\begin{aligned} \mathbf {N}&:= w^{(\nu )}_{\mathrm{sst}} \mathrm {reo} \left( \begin{bmatrix} \mathbf {N}_{\mathrm{sst}} & 0\\ 0 & 0 \end{bmatrix} \right)\\ &\quad +\sum _{s=1}^{S} \sum _{c \in C} w^{(\nu )}_{c,s}\mathbf {N}_{c,s} \nonumber \\&\quad +\sum _{a \in \{\text {np},\,\text {sp}\}} w^{(\nu )}_{a}\,\mathrm {reo} \left( \begin{bmatrix} \mathbf {N}_{\mathrm{a}} & 0\\ 0 & 0 \end{bmatrix}\right) \nonumber \\&\quad +w^{(\nu )}_{\mathrm{hd}}\,\mathrm {reo} \left( \begin{bmatrix} \mathbf {N}_{\mathrm{hd}} & 0\\ 0 & 0 \end{bmatrix}\right)\end{aligned}$$

(30)

$$\begin{aligned} \mathbf {n}&:= w^{(\nu )}_{\mathrm{sst}}\,\mathrm {reo} \left( \begin{bmatrix} \mathbf {n}_{\mathrm{sst}} & 0\\ 0 & 0 \end{bmatrix}\right) + \sum _{s=1}^{S} \sum _{c \in C} w^{(\nu )}_{c,s}\mathbf {n}_{c,s} \end{aligned}$$

(31)

$$\begin{aligned} \lambda&:= w^{(\nu )}_{\mathrm{sst}}\lambda _{\mathrm{sst}} + \sum _{s=1}^{S} \sum _{c \in C} w^{(\nu )}_{c,s}\lambda _{c,s}. \end{aligned}$$

(32)

introducing the $4S+4$ weighting factors $w^{(\nu )}_*$ for all normal equations. With initial values $w^{(0)}_* = 1$, the normal equations

$$\begin{aligned} \mathbf {N} \mathbf {x}&= \mathbf {n} \end{aligned}$$

(33)

can be solved for the spherical harmonic coefficients arranged in $\mathbf {x}$ via the Cholesky decomposition followed by forward and backward substitution, which yields the least squares estimate $\hat{\mathbf {x}}^{(\nu )}$.

With $\hat{\mathbf {x}}^{(0)}$, $\mathbf {N}_{*}$, $\mathbf {n}_{*}$, $\lambda _{*}$ and the number of observations used to compute the normal equation $n_*$, the weights $w^{(\nu )}_*$ can be estimated iteratively as inverse variance components (Förstner 1979; Koch and Kusche 2002). The procedure of reassembling (30) and (32) and solving for $\hat{\mathbf {x}}^{(\nu )}$ is iterated until convergence of the weights (cf. Brockmann 2014, Chap. 4).

After modeling the stochastic characteristics of the observations when setting up the normal equations and conversion of the weights, a realistic covariance matrix reflecting the uncertainty of the spherical harmonic coefficients is obtained by

$$\begin{aligned} \varvec{\Sigma }\left\{ \hat{{\varvec{{\mathcal {X}}}}} \right\} = \mathbf {N}^{-1}. \end{aligned}$$

(34)

3 Processing Details and Intermediate Results

Within the following, special attention is paid to the analysis of the gravity gradient observations, as they mainly contribute to the final result and are the observations which make the GOCE mission unique. Characteristics of the reprocessed input data compared to the previous calibration and the influence of the data improvement on gravity field determination are shown, analyzing sub-solutions and comparing them to their RL05 equivalents. Here, the observation types are accessed individually, before the finally combined final model is presented in the upcoming section.

3.1 Improvement of the Longer Wavelength Gravity Field Derived from Kinematic Orbits

Compared to the GO_CONS_EGM_TIM_RL05 solution, the reprocessed kinematic orbit data set was used as the input data set for gravity field recovery (cf. Sect. 2.2). Within the orbit reprocessing, ionosphere-induced effects were mitigated from the GPS phase observations by dedicated observation weighting (Grombein et al. 2019), which significantly reduce the magnetic equator artifacts which were present in the RL05 solutions (e.g., Jäggi et al. 2015; Grombein et al. 2019).

Fig. 2

Geoid heights of 400 km Gaussian-filtered RL05 and RL06 SST-hl solutions truncated at spherical harmonic degree 120 with respect to the XGM2016 model

To demonstrate the progress, the SST-hl normal equations were solved without regularization to determine a pure SST-hl solution. Figure 2 shows 400 km Gaussian-filtered SST-hl solutions in terms of geoid heights up to degree 120 with respect to the XGM2016 model as reference. The RL05 solution in Fig. 2a shows the systematic error caused by the ionosphere-induced error in the RL05 solution as a strong and sharp signal around the magnetic equator. Whereas the root mean square error (outside the polar gap) is 0.3 cm, the minimum and maximum reach − 1 cm and 1.7 cm, respectively. In Fig. 2b, it is shown that the effect of the magnetic equator is less pronounced in the RL06 SST-hl solution. Although its extent is slightly increased, the effect is less pronounced and the maximum values are in the range of $\pm 0.7\,\mathrm{cm}$. In addition, the root mean square error slightly decreased to 0.2 cm. This indicates the main progress of the SST-hl contribution.

3.2 Results from the Gravity Gradient Analysis

3.2.1 Stochastic Characteristics of Reprocessed Gravity Gradients

As discussed in Sect. 2.3, the gravity gradients were initially partitioned into gapless segments. After the insertion of new segments, the procedure converged to 85 segments for the RL06 analysis. Compared to RL05, Fig. 3a indicates that longer sequences are available, especially for the lower orbit phases toward the end of the mission. This shows that the data are more stationary and the new robust filter strongly increases the stability of the decorrelation filter estimation. A comparison to the previous time-wise RL01 to RL04 is shown as well, indicating the amount of data used for their computation. In addition, Fig. 3b focuses on RL06, there are 32 segments of the time-series which are longer than a week (vertically centered) which contain the majority of the observations. Furthermore, there are 15 segments which are shorter (less than a week), but contain usable data and are used for gravity field recovery (shifted slightly up to make them visible). The remaining 38 segments are short as well, but were completely rejected in the processing (gray, slightly shifted down). They typically correspond to satellite maneuvers or special events (cf. GOCE Flight Control 2014).

Fig. 3

Schematic overview of the gravity gradient data used in the generation of the time-wise models

All segments were sequentially processed to detect the suspicious data and to estimate the AR process as a decorrelation filter. Table 4 in “Appendix” lists the begin, the end and the amount of available epochs for all 47 usable segments. Furthermore, the usable number of observations and the percentage of detected suspicious data are provided individually for the components. In total, the level of detected outliers is slightly higher for $V_{xx}$ and $V_{xz}$ (1.9%, 1.6% respectively) compared to $V_{yy}$ and $V_{zz}$ (0.6%). Within the RL05 processing, the amount of outliers was only slightly higher, i.e., in the range of 2% to 3%². Details on the sensitivity of the different statistical tests used to identify the suspicious epochs are not discussed here. Examples are discussed in Schubert et al. (2019).

As the decorrelation filters are an estimated model for the noise characteristics of the gravity gradient data, the improved quality of the reprocessed gravity gradients can be nicely demonstrated comparing the power spectral densities (PSD) of the filters from the RL06 to the RL05 solution, which was based on the old gravity gradient data set. They are summarized for all components in Fig. 4, for RL06, only the stable estimates for the long segments ($> 7\ \hbox {days}$, cf. Fig. 3b are shown). The color code for RL06 corresponds to the segment IDs from Table 4. Compared to RL05, the following conclusions for the noise characteristics can be drawn

• The noise characteristics are more stationary in the reprocessed gravity gradient data set. Over the entire mission, the spectra of the noise estimates are nearly constant within the same component: There are only two exceptions, a) there are some variations in the small frequencies (1e−4 Hz to 5e−4 Hz) for the $V_{yy}$ noise component. It is interesting to see that it decreases toward the end of the mission, thus with orbit lowering. b) The noise level for the $V_{zz}$ gradient is higher in the beginning of the mission. This is a known fact, that after an anomaly on the satellite, the $V_{zz}$ performance was improved. The reason for that is still unknown.

• The frequency range with a flat white noise spectrum was extended toward smaller frequencies: For all components, it was shifted below the measurement band toward 2e−3 Hz from about 5e−3 Hz for the old data set (except $V_{yy}$).

• The general level of the noise in the measurement band was not significantly decreased: Except $V_{yy}$, which was strongly distorted, the level of the flat spectrum has not significantly changed.

• Below the measurement band, the noise level was reduced up to one order of magnitude: The long wavelength errors were significantly reduced.

• The $V_{yy}$ component shows the most impressive improvement: The new data set has now flat noise characteristics in the measurement band as well. The large variations are gone and the level of the noise decreased.

• The low orbit data, even at the lowest orbit, have now the same if not better quality compared to the nominal orbit data.

Fig. 4

Comparison of the PSD of the decorrelation filter estimates for RL05 (left column) and RL06 (right column) for all used gravity gradient components. Shown are the estimates for the 32 long segments. The color code corresponds to the segment numbers. Note that due to the different segmentation, colors for RL05 and RL06 may correspond to different time periods

All together, the estimated filters highlight the improved quality of the reprocessed gravity gradients and the new robustified decorrelation filter estimation procedure.

3.2.2 Performance of Reprocessed Gravity Gradient Data for Gravity Field Recovery Compared to RL05

With the suspicious data flag information and the estimated decorrelation filter models at hand, least-squares normal equations are set up for all of the segments and for all of the gravity gradient components cf. (14) and (15). In case a segment contains observations of about seven to ten days length, the segment and component-wise normal equations

$$\begin{aligned} \mathbf {N}_{c,s}\mathbf {x}_{c,s}&= \mathbf {n}_{c,s} \end{aligned}$$

(35)

are typically solvable for the parameters $\mathbf {x}_{c,s}$ and the accuracies³,

$$\begin{aligned} \varvec{\Sigma }\left\{ \mathbf {x}_{c,s}\right\}&= \mathbf {N}_{c,s}^{-1}. \end{aligned}$$

(36)

The solution of (35) represents a gravity gradient only estimate based on a single gradient component c from the data period of segment s. Although such solutions are of low quality, they can be used to access the quality of the components, the data of the period as well as the consistency. In fact, after the assembly of (35), the analysis software tries to solve for $\mathbf {x}_{c,s}$ by default, such that, in case $\mathbf {N}_{c,s}$ is positive definite, all computable intermediate sub-solutions are available.

Figure 5 shows an example solution for segment $s=3$ (cf. Table 4) and all four components in terms of degree (error) variance excluding the near zonal coefficients with order $m(l) < 0.113 l + 3$ (based on the rule found by van Gelderen and Koop 1997) as no regularization is applied. The chosen time period corresponds to the RL01 time series in the very beginning of the mission. As a reference model, to highlight consistency, EGM_TIM_RL05 is chosen, which is assumed to be more accurate, as based on all components and the entire mission data set, but of course it is not independent. This is an example of a quite long segment and is representative for the data in the nominal mission at the nominal orbit. For comparison, a RL05-like solution is shown (but truncated at degree 250).

Fig. 5

Degree (error) variances excluding the near zonal coefficients of component-wise gravity gradient-only solutions for segment $s=3$, cf. Table 4 of the RL06 processing (orange). As reference model EGM_TIM_RL05 is chosen. As comparison, the RL05 equivalent is shown in blue. The formal errors of all three solutions are shown as dashed lines of same color. Note that due to the good correspondence of empirical and formal error estimates, the lines overlap to a large extent

Looking at the solutions shown in Fig. 5, their absolute quality is low for the segment-wise and single gradient component analysis. Most important is the relative comparison of the RL05 and RL06 solutions to the reference model. The first conclusion is that there is no obvious situation, where the RL05 solution is better than the RL06 solution. Above degree 80 they seem to be identical. Below, RL06 seems to be more consistent to EGM_TIM_RL05 than the original RL05 sub-solution. The improvement for the lower degrees is not large, but significant for some degrees, most obviously seen for the xx, xz and yy components. This figure indicates that all the (sub-)solutions improved.

More obvious is the change of the formal error characteristics from RL05 to RL06 for the component wise sub-solutions. Again, especially in the low degrees and for xx, xz and yy, the formal error estimates decrease significantly. Due to the agreement of the formal and empirical error estimates of RL06 (solid and dashed lines correspond quite well), it is an indication that the newly derived filters (cf. Sects. 2.3.2 and 3.2.1) better reflect the true error characteristics. Furthermore, it can be seen that the error estimates of the RL05 sub-solutions were generally to pessimistic. It is expected and demonstrated later that the improved formal error characteristics are important as soon as different segments and components are combined, as this controls the degree and order dependent relative weighting. An improved formal error means that the individual contribution matrices $\mathbf {N}_{c,s}$ are more compatible and will produce better combination where the strengths of the different components are better kept. This will result in a more accurate combined solution. For all components, an improved quality of the reprocessed gravity gradient observations, even in the beginning of the mission, can be observed on the level of gravity field sub-solutions.

Fig. 6

Degree (error) variances excluding the near zonal coefficients of component-wise gravity gradient-only solutions for segment $s=20$, cf. Table 4 of the RL06 processing (orange). As reference model EGM_TIM_RL05 is chosen. As comparison, the RL05 equivalent is shown in blue. The formal errors of all three solutions are shown as dashed lines of the same color. Note that due to the good correspondence of empirical and formal error estimates, the lines overlap to a large extent

A similar analysis is shown by Fig, 6, but now for segment $s=20$ which is observed at an altitude already lowered by 30 km, where the satellite operated in a rougher environment. These low orbit data are very important for the final model, as the observations are captured with the best signal-to-noise ratio. It can be seen that the conclusions for the low degrees of Fig. 5 now hold true for the entire spectrum. Not only the formal errors improved and are more realistic, the solutions improved as well. Especially for yy, the huge progress achieved by the advanced calibration becomes visible with a huge reduction of the formal as well as the empirical error. The improvements are visible up to the highest degrees and even for the other components, but not as large.

The two segments were selected to be representative for the entire mission, one quite early and one quite late segment. The analysis could have been done with the same conclusions for any segment and component which is solvable.

Furthermore, it can be easily verified that the combination of different components improves. It is straightforward to combine, e.g., for both test segments $s=3$ and $s=20$, the four components on normal equation level using either variance component estimation to derive a weight or just using $w_{c,s}=1.0$ which is a good guess (cf. Sect.  3.2.3). The presentation of the results is not done here, as the combination of the segments and components of the entire data set reflect the same on average for the entire mission.

3.2.3 Combined Gravity Gradient Solution

All gravity gradient normal equations for the different segments and components can be combined in different ways. For that purpose, variance component estimation (see, e.g., Koch and Kusche 2002) is used to derive weights for the gravity gradient normal equations $w_{c,s}$. Different approaches are possible, i.e., weights can be computed combining only the gravity gradient normal equations and iterating variance component estimation or, directly combining the regularization matrices and the SST normal equations to stabilize the solution.

Within the computation of RL06, the latter approach was chosen, already providing weights for the SST and the regularization matrices, which will be discussed later. A further test for consistency and the applied stochastic model for the gravity gradients are the values estimated for the weights. In the case of a proper decorrelation filter estimation—which includes an estimate for the variance—$w_{c,s} = 1.0$ is expected for all $c\in \{xx,\,yy,\,zz,\,xz\}$ and $s\in \{1,\ldots ,S\}$. Actually, from the 188 estimated weights for the gravity gradient normal equations, only three weights are outside the interval $\left[ 0.99, 1.01\right]$ (rounded to the second digit). But these three (i.e., 1.13, 0.91 and 1.08) are still quite close to 1.0. Two of them belong to a short segment for which the filter estimation is by nature not that stable, for the third, no special characteristics could be identified.

Only the weights $w_{c,s}$ are extracted after convergence to compute a combined solution for all components $c\in C$, i.e., combining $\mathbf {N}_{c,s}$ and $\mathbf {n}_{c,s}$ for all $s\in S$ to obtain component-wise normal equations

$$\begin{aligned} \mathbf {N}_{c}&:= \sum _s w_{c,s} \mathbf {N}_{c,s} \end{aligned}$$

(37)

$$\begin{aligned} \mathbf {n}_{c}&:= \sum _s w_{c,s} \mathbf {n}_{c,s} \end{aligned}$$

(38)

and the corresponding solution together with its uncertainty estimate

$$\begin{aligned} \mathbf {N}_{c}\mathbf {x}_c&= \mathbf {n}_{c} \end{aligned}$$

(39)

$$\begin{aligned}&\Rightarrow \tilde{\mathbf {x}}_c = \mathbf {N}_{c}^{-1}\mathbf {n}_{c} \end{aligned}$$

(40)

$$\begin{aligned}&\Rightarrow \varvec{\Sigma }\left\{ \tilde{\mathbf {x}}_c\right\} = \mathbf {N}_{c}^{-1}. \end{aligned}$$

(41)

These sub-solutions are useful to study the sensitivities of the different gradient types for the gravity field recovery. The four available (sub-)solutions are again shown in terms of degree (error) variances, to assess the solutions derived from the three main diagonal components and from $V_{xz}$. Here, the XGM2016 (Pail et al. 2018) model is used as reference, which is built upon GOCO05S as background model and thus strongly depends on EGM_TIM_RL05 (further supported by GRACE, satellite altimetry and terrestrial data) in the degree range relevant for GOCE. Figure 7 shows again the consistency of empirical error estimates with respect to the reference model (colored solid lines) and the corresponding formal error estimates computed from the coefficients’ variances extracted from the diagonal of (41). Furthermore, the relative accuracy level of the components can be seen. Degree error variances from Fig. 7 conclude that $V_{zz}$ is the most important component, which is not surprising as the nearly radial gradient contains the strongest signal. Surprisingly, Fig. 7 suggests that the in-flight direction $V_{xx}$ is—seen per degree—the least important observation as it produces the least accurate solution.

Fig. 7

Degree (error) variances excluding the near zonal coefficients strongly influenced by the polar gap. Although not independent, XGM2016 (black) is used as reference. Empirical error estimates with respect to the reference model are shown by solid lines, formal error estimates by dashed lines. The different colors represent the four individual component-wise gravity gradient solutions

To demonstrate that this is a kind of wrong impression triggered by the degree variance representation, Fig. 8 shows the formal standard deviations of the coefficients in the triangular scheme. It can be seen that the order-dependence varies for the different gradient components. Whereas the $V_{zz}$ gradient-based solution hardly shows any order dependence outside the near-zonals, the other three do. The $V_{xx}$ and $V_{xz}$ based solutions show an increase toward the sectorial coefficients, especially $V_{xx}$ having the maximum value there. The degree error variances are dominated by the largest standard deviations, yielding a wrong impression for the importance. Consequently, only a coefficient-wise analysis yields a reasonable conclusion. Please note as well that the analysis of the degree variances neglected all correlations in (41), so the results should not be overinterpreted.

Fig. 8

Standard deviation of the individual spherical harmonic coefficients for the different sub-solutions from the individual gravity gradients in the triangular scheme

Finally, the gravity gradient solution based on all four components is computed, which is then compared again to the RL05 counterpart to demonstrate the improvements. This solution from all epochs and all components follows from the assembly of

$$\begin{aligned} \mathbf {N}_{\mathrm{SGG}}&:= \sum _c \mathbf {N}_{c} \end{aligned}$$

(42)

$$\begin{aligned} \mathbf {n}_{\mathrm{SGG}}&:= \sum _c \mathbf {n}_{c} \end{aligned}$$

(43)

and the solution of

$$\begin{aligned} \mathbf {N}_{\mathrm{SGG}}\mathbf {x}_{\mathrm{SGG}}&= \mathbf {n}_{\mathrm{SGG}} \end{aligned}$$

(44)

$$\begin{aligned}&\Rightarrow \tilde{\mathbf {x}}_{\mathrm{SGG}} = \mathbf {N}_{\mathrm{SGG}}^{-1} \mathbf {n}_{\mathrm{SGG}} \end{aligned}$$

(45)

$$\begin{aligned}&\Rightarrow \varvec{\Sigma }\left\{ \tilde{\mathbf {x}}_{\mathrm{SGG}}\right\} = \mathbf {N}_{\mathrm{SGG}}^{-1}. \end{aligned}$$

(46)

Figure 9 again uses the degree (error) variance representation to display the relative performance from RL05 to RL06. Again, XGM2016 is chosen as a reference model. The formal improvement can be seen by the degree error variances computed from the formal errors. It can be seen that improvements are visible for the entire spectrum. Furthermore, the empirical error estimates with respect to the reference model indicate that the RL06 solution is more consistent to XGM2016 than the RL05 solution. This is a remarkable quality indicator, as due to the usage of GOCO05S in the generation of XGM2016 it is biased toward EGM_TIM_RL05. Consequently, also the empirical error estimates contain a significant amount of XGM2016 errors, especially in the range of degree 80 to 220. To quantify the improvements, cumulative degree variances can be used. The cumulative formal error of RL06 improves by about 50% from degree 16 to 100 and about 25% for degree 200. For degree 280 it is still about 20%.

Fig. 9

Degree (error) variances excluding the near zonal coefficients for the solutions based on gravity gradients from RL05 (blue) and RL06 (orange). Although not independent, XGM2016 (black) is used as reference. Empirical error estimates with respect to the reference model are shown by solid lines, and formal error estimates by dashed lines

3.3 Effect and Weights for the Zero Regularization

As discussed before, regularization of the signal toward zero is applied to constrain the polar areas and to smooth the high degree coefficients. The high degrees were constrained with zeros using Kaula’s rule to derive mean variances for the zero pseudo observations of coefficients above degree 200. Variance component estimation is used to estimate a weight, which is estimated as $w_{\mathrm{hd}}=0.7808$.

More interesting are the weights for the zero gravity anomaly observations. After convergence, they are $w_{\mathrm{np}}=1.2672\times 10^{8}$ and $w_{\mathrm{sp}}=2.36357\times 10^{7}$ for the north and south polar gap. These weights correspondent to a standard deviation of $\sigma _{np} = 8.88\,\mathrm{mGal}$ in the northern gap and $\sigma _{sp} = 20.57\,\mathrm{mGal}$ in the southern, respectively. The ratio 2:1 roughly corresponds to the expected signal content in the caps as reflected by high resolution combination models using polar terrestrial data.

Fig. 10

Geoid differences of different constrained GOCE models with respect to XGM2016 at degree 200 (m)

Figure 10 shows two combined GOCE solutions (including SST-hl, gravity gradients and regularization), using both possible options to constrain the caps. Whereas the left column applies a strong regularization of the near zonal coefficients toward zero (setting up (27) for the near zonal coefficients with weight $w_{\mathrm{nz}} = 10.0$), the right column shows the solution using the zero gravity anomaly pseudo observations. Shown are geoid differences at degree 200 with respect XGM2016. It is obvious that all GOCE-based models produce a poor geoid inside the polar caps, but at least the gravity anomaly-based approach seems to produce a little less artifacts outside them. Consequently, the approach is selected as it seems to act more locally. To overcome the general deficit, an extended version of the GOCE time-wise (TIM_RL06E) model was compiled, which inserts data from several polar terrestrial campaigns (for details see Zingerle et al. 2019b) replacing the zero pseudo observations.

4 Final Results: The Time-Wise RL06 Solution

This section is not indented to completely validate the final results but to demonstrate the characteristics and to indicate some improvement with respect to the RL05 predecessor model.

The final model is computed setting up (30) and (31) and then by solving

$$\begin{aligned} \mathbf {N} \mathbf {x} = \mathbf {n} \end{aligned}$$

(47)

yielding

$$\begin{aligned} \mathbf {x}&= \mathbf {N}^{-1} \mathbf {n}, \end{aligned}$$

(48)

$$\begin{aligned} \varvec{\Sigma }\left\{ \mathbf {x}\right\}&= \mathbf {N}^{-1} \end{aligned}$$

(49)

the combined spherical harmonic parameters as well as the uncertainty estimate in terms of the full error covariance matrix.

4.1 Spectral Characteristics of the Model

In terms of degree error variances, the spectral characteristics are shown in Fig. 11 using the XGM2016 model as a reference. Figure 11a shows the individual solutions as comparison, i.e., the SST solution cf. Sect. 3.1 and the gravity gradient solution cf. Sect. 3.2.3. Near-zonal coefficients are excluded. It can be concluded that the orbit data mainly determine the solution up to degree 12 and still positively contribute to degree 40. Between degree 40 and 220, the solution is more or less completely identical to the gravity gradient solution and, starting above degree 220, the high degree regularization starts to damp the signal toward zero. Furthermore, except for the first 20 degrees, the error estimates seem to provide a consistent estimate. For the degree range from degree 80 to 220, the empirical error is above the formal error estimate, which corresponds to the XGM2016 error which is included in the empirical error estimate.

Figure 11b shows the combined RL06 solution in comparison with the RL05 counterpart, fading out the highly correlated near zonal coefficients, influenced still by the polar gap. Improvement is seen for all degrees above 20, slightly decreasing toward the higher degrees. In terms of degree variances, hardly a change in the low degrees is visible. In contrast to that, Fig. 11c includes the badly determined and highly correlated zonal coefficients. Consequently, the figure shows an increase in the degree (error) variances. In addition, a mismatch for the formal and empirical error estimates occurs. This becomes visible because the degree variance computation for the formal errors neglects the very strong correlations between coefficients of the same order. Consequently, the degree variance plot is hardly interpretable in a reasonable way. Nevertheless, some improvement can be seen from RL05 to RL06, but this is mainly contributed to the fact of a changing regularization approach.

Fig. 11

Characteristics of the combined solution in terms of degree (error) variances. As reference model XGM2016 is used (black). Solid lines are empirical errors with respect to the reference model, and dashed lines the formal errors

4.2 Analyzing the Contribution

Solution vector and full covariance matrices in the combination procedure on normal equation level cf. (30) and (31) can be seen as direct pseudo observations of the parameters, resulting in the design matrix $\mathbf {A}_i=\mathbf {I}$ for any of the involved groups g. Seeing the individual solution vectors $\mathbf {x}_{g}$ as observations, the local redundancy number $r_{g,i}$ for any ’observation’ i from $\mathbf {x}_g$ can be computed, or instead, directly the local contribution number $h_{g,i} = 1 - r_{g,i} = \mathbf {H}_g(i,i)$ which follows from the diagonal of

$$\begin{aligned} \mathbf {H}_g = \varvec{\Sigma }\left\{ \tilde{\mathbf {x}}\right\} \varvec{\Sigma }\left\{ \tilde{\mathbf {x}}_{g}\right\} ^{-1} = \varvec{\Sigma }\left\{ \tilde{\mathbf {x}}\right\} \mathbf {N}_{g} \end{aligned}$$

(50)

as a measure of contribution of the observation i from group g to the adjusted observations (which are in this case the combined parameters $\mathbf {x}$). As

$$\begin{aligned} \sum _{g \in G} h_{g,i} = 1 \end{aligned}$$

(51)

holds, $h_{g,i}$ can be seen as the contribution of group g to determine the i-th parameter (see, e.g., Sneeuw (2000, Sect. 6.2), Yi et al. (2013), Brockmann (2014, Sect. 5.4)). Note that in the case of correlated observations—this is the case here as we use parameters from an adjustment—$h_{g,i}$ can get negative and can be larger than 1.

The contribution can be computed for every data set used in (30), i.e., SST, gravity gradient components, high-degree and polar gap regularizations individually. The contribution measures are shown per group in Fig. 12 in the triangular scheme. The analysis highlights the importance and the complementary of the individual gravity gradients. Furthermore, it is shown that the SST solution contributes to the sectorial coefficients even at degree 100 with still 10% to 20%. The regularization in the polar areas only contributes to the near-zonals. The percentage is quite low, as the regularization acts as a decorrelation of the zonal coefficients. Due to the strong correlations at the low orders, the absolute number is hard to interpret. Finally, the regularization for the higher degree starts to contribute from degree 220 onwards. Above degree 280, it is dominating the observation and smooths the signal toward zero.

Fig. 12

Contribution analysis for all input and pseudo observations types, i.e., the diagonal of the matrix (50) in the triangular scheme

4.3 Spatial Characteristics and Quantifying the Improvement

As a final result, the improvement from RL05 to RL06 is shown spatially evaluating the spherical harmonic series. For that purpose, geoid heights are synthesized on a regular global grid as a difference to the XGM2016 model. The same could have been done with respect to EGM2008, the pre-GOCE era state-of-the-art global gravity field model, but these differences are dominated by the EGM2008 errors and will not obviously show the progress from RL05 to RL06. Figure 13 shows the situation for a truncation degree of 200 for both models, i.e., RL05 and RL06. It can been seen again that RL06 is more consistent to XGM2016 in terms of systematic differences and noise, although XGM2016 is biased toward RL05. This again underlines a quality improvement of RL06 and the quality of the combined model XGM2016. Obviously, the artifact around the magnetic equator present in RL05 (cf. Sect. 3.1) is significantly reduced. Nevertheless, systematic differences are visible over land which occur mostly in regions where the additional gravity data sets used in the compilation of XGM2016 are of lower quality (see, e.g., Fecher et al. 2017, Fig. 5). In those regions, XGM2016 cannot serve as a reference model, and their statistics like the root mean square error are unreliable and cannot serve as a quality measure for the RL06 solution. To assess the global mean improvement, one has to rely on the formal error estimates which are provided in Table 2. It summarizes the minimal, maximal and mean standard deviation of geoid heights and gravity anomalies which results from propagating the full covariance matrix of the spherical harmonic coefficients. The numbers are restricted to latitudes outside the polar gap. Thus, $-80^{\circ }$ to $80^{\circ }$ is used for the computation. Different truncation degrees are shown in order to show the accuracy depending on the spherical harmonic degree.

Table 2 Results of the variance-covariance propagation of the RL05 and RL06 solutions to geoid heights and gravity anomalies for the area outside the polar gaps depending on the truncation degree

On the one hand, the formal estimates can be used to quantify the expected mean accuracy of the model. For GOCE target resolution of 100 km it was 2.3 cm for RL05 and decreased now to 1.8 cm, which shows that the mission requirement was achieved. For gravity anomalies, it is even more obvious, it was already well below 1 mGal for RL05 (0.65 mGal) and is now down to 0.5 mGal for RL06. The 1 mGal level is now achieved at a spatial resolution of 90 km.

On the other hand, it can be used to numerically quantify the mean global improvement of the accuracy from RL05 to RL06. The mean standard deviation for the geoid reduces more than 22% at degree 200, for the geoid as well as for the gravity anomalies. The extreme values reduce as well. This is systematically the case for all truncation degrees. But for the highest common degree which is 280, reduction is only about 7–8%, which shows a decreased improvement with higher truncation degree. This has to be expected, as the quality of the reprocessed gravity gradients improved most for the lower frequencies (cf. Fig. 4). Nevertheless, the formal improvement of more than 22% at degree 200 is remarkable. Note that these values do not capture the reduced systematic errors as, e.g., the magnetic equator (cf. Sect. 3.1), which is not modeled in the covariance matrix which is constructed assuming stationarity.

Fig. 13

Geoid height differences of the RL05 and RL06 model truncated at spherical harmonic degree 200 with respect to the XGM2016 model

Furthermore, in regions where the dynamic ocean topography is smooth and no strong currents occur, the high-frequency noise is visibly reduced (cf. Fig. 13). Within this area, the marine gravity can be determined precisely from satellite altimetry, as the separation of the sea surface height into geoid and dynamic ocean topography is reliable. One of these smooth areas is located in the South Pacific Ocean, which is highlighted in the zoom in the right column of Fig. 13. This area is chosen to support the visual impression with some statistics. Those carefully chosen local areas are expected to be of high accuracy, so that the assumption that the XGM2016 which incorporates the altimetry derived marine gravity can serve as a reference might be true. Note that this is definitely not the case for the entire ocean and all global combination models.

For that purpose, Table 3 shows empirical statistics computed in that area (i.e., − 150° to − 120° longitude and − 45° to − 20° in latitude direction) from the differences, but now for truncation degrees from 180 to 300 in steps of 20 degrees and for geoid heights as well as gravity anomalies. In addition, the formal accuracy derived from covariance propagation to the corresponding quantity is determined for the specific truncation degree. The provided value is the mean value for that specific region. Formal and empirical improvement from RL05 to RL06 can be determined from the obtained values.

Table 3 Statistics of the differences of RL05 and RL06 to the XGM2016 model in terms of geoid heights and gravity anomalies for various truncation degrees for a small area in the South Pacific Ocean (cf. Fig. 13b, d) where the dynamic ocean topography is smooth

First of all, a decrease of all formal error is obtained from RL05 to RL06 as for the global mean. In contrast to that, the empirical errors (root mean square error—rms) show a slight increase in the rms for degree 280 for the geoid and 260 and 280 for the gravity anomalies. This increase happens due to the different regularization applied for the high degrees, including the implicit regularization due to the different maximal expansion.

Secondly, the statistics can be used to validate the formal errors. For RL06, a very nice consistency can be seen for rms values and the formal counterparts. For degrees 180 to 240, there is an excellent match for geoid heights as well as gravity anomalies. For degrees 260 to 300, the formal error is always larger compared to the rms. This very likely depends on the regularization (global average) in combination with the local signal characteristics of the studied area.

RL05 shows very similar characteristics, but the formal errors are always slightly larger than the rms. This is attributed to two facts. Firstly, XGM2016 is biased toward RL05 as used in the computation, though smaller rms values can happen. Secondly, due to the decorrelation filter estimation procedure in RL05, it is known that the formal error estimates of RL05 are pessimistic. This test shows that the quality of the error covariance matrix for RL06 improved.

The mean global formal improvements (more than 22%) are confirmed by the local formal errors in this specific area. Furthermore, improvements can be derived from the empirical error estimates. For this area, it is slightly less, about 13%. From these tests, it can be concluded that the improvement of RL06 compared to RL05 is somewhere in between 10 and 25%. Again, this does not include the reduction of the systematic effects, which were at centimeter level for the magnetic equator, which is an improvement locally in the order of 50%.

5 Summary and Conclusions

Based on the reprocessed gravity gradients of the GOCE mission and using an advanced processing approach, the new global gravity field model GO_CONS_EGM_TIM_RL06 was computed as a successor of the RL05 model published in 2014. In general, it follows the philosophy of the previous GOCE time-wise models, with the basic idea that it is based on GOCE observations only. As typically for global gravity field models, it is provided as a spherical harmonic expansion, truncated at degree 300. As done for the previous time-wise models, it is tried to provide an accuracy information of the spherical harmonic coefficients that is as realistic as possible, in form of the full error covariance matrix.

This paper provides an overview of the processing strategy with important details of the processing (cf. Sect. 2). Special focus is put on the gravity gradient processing and estimation of the decorrelation filters which are used as a stochastic model of the gravity gradients, as this was improved with respect to the RL05 analysis. In addition, as it is essential to determine the long wavelengths in any GOCE-only model, the short arc approach to set up the normal equations from the kinematic orbit positions is summarized. The required regularization strategies were discussed, and in the end we presented the combination procedure which is applied to derive the combined GOCE model using the kinematic orbits, gravity gradients as well as the regularization to derive a stable solution.

Section 3 describes the processing as well as the intermediate results, which includes the description of the used data, the characteristics of the data, as well as the characteristics and consistency of the gravity gradient solutions. For that purpose, intermediate models were produced, using only sub-sets of the observations which are analyzed with respect to their characteristics. In addition, a comparison to RL05 counterparts is performed to highlight the gain of the reprocessing campaign.

In the end, the characteristics of the final solution and the associated error estimates are summarized in Sect. 4. Again, focus is on indicating the improvement with respect to RL05 and some general quality analysis. A specific validation with external data is expected in the near future. The intention is to put the focus on the general characteristics of the model.

To conclude, an improved version of the GOCE time-wise solution and its error covariance matrix is available as the GO_CONS_EGM_TIM_RL06 model. It is shown that the improvements with respect to the previous RL05 solution are:

• The mean global accuracy decreased by 10−25% compared to RL05.

• Due to the improved decorrelation filter estimation, the quality of the covariance matrix was increased, yielding a more realistic error description.

• Systematic errors at centimeter level were reduced (e.g., magnetic equator in SST contribution or artifacts in gravity gradients reduced in re-calibration).

• The mission objectives of an accuracy of 1 mGal in terms of gravity anomalies are reached for spherical harmonic degrees below degree 220, i.e., up to 90 km spatial resolution.

Pessimistic accuracy estimates for the GO_CONS_EGM_TIM_RL06 solution are 2 cm in terms of geoid and 0.5 mGal in terms of gravity anomalies at GOCE target resolution of 100 km. The solution is already available at ICGEM (ICGEM 2020; Brockmann et al. 2019) and via the ESA GOCE Online Dissemination Service (ESA GOCE-ODS 2020) together with its full error covariance matrix as the official ESA time-wise solution. Furthermore, as discussed in Sect. 3, gravity gradient-based (sub-)solutions, normal equations of the contributors or unconstrained versions will be made available upon request. The formal equations were already used to construct the combined satellite only model GOCO06S (Kvas et al. 2019b), the polar extended model TIM_R6E (Zingerle et al. 2019b) and the high-resolution combined model XGM2019e (Zingerle et al. 2019a).

Appendix: Segmentation Used for Gravity Gradients

See Table 4.

Table 4 Final segmentation used for gravity gradient processing