1 Introduction

Knowledge of temporal variations of the Earth’s gravity field is of importance to understand large-scale mass transport at and below the Earth’s surface. In the past decade, it was mostly observed with the ultra-precise K-Band ranging (KBR) measurements from the gravity recovery and climate experiment (GRACE) mission (Tapley et al. 2004). However, the GRACE mission was completed in October 2017, so that there exists a gap between GRACE and its successor, GRACE Follow-On (GFO), which was launched in May 2018. Nowadays, there is a consensus that GPS-based high-low satellite-to-satellite tracking (hl-SST) can play an important role in bridging the gap (Bezděk et al. 2016; Guo et al. 2017a; Guo and Zhao 2019; Jäggi et al. 2014, 2016; Visser et al. 2014; Weigelt et al. 2013). For that purpose, kinematic orbits are usually derived from the hl-SST GPS data by a precise point positioning approach (Švehla and Rothacher 2005) and then taken as pseudo-observations in gravity field modeling. In the last decade, several techniques have been developed to convert kinematic orbits into a gravity field model (Baur et al. 2014). In this study, we will focus on the classical dynamic approach (CDA) (Reigber 1989). Our general goal is further improving the accuracy of gravity field modeling from kinematic orbit data.

Since GPS data are about 3 orders of magnitude less precise than the KBR data, they are only sensitive to temporal variations at low spherical harmonic degrees (< 20) of the Earth’s gravity field. Furthermore, it is well known that kinematic orbits are very sensitive to observation geometry and various systematic errors, e.g., mismodeling of GPS antenna phase center, high-order ionosphere-induced errors, near-field multipath, and uncalibrated hardware delays (Bock et al. 2014; Jäggi et al. 2009; Montenbruck et al. 2017). As a result, the recovered gravity field solutions usually suffer from systematic errors, particularly at low degrees (Jäggi et al. 2011a). Therefore, efforts are still ongoing to further improve the GPS-based gravity field solutions. For that purpose, previous authors mostly focused on developing dedicated algorithms to model systematic errors in GPS data, e.g., phase center variations (PCVs) (Bock et al. 2011; Jäggi et al. 2009) and high-order ionosphere-induced errors (Jäggi et al. 2014; Zehentner and Mayer-Gürr 2015). However, systematic errors may still exist in the kinematic orbits due to model deficiencies and/or unknown error sources. In this research, our focus is on the systematic errors that are constant or slowly varying in time, possibly showing a strongly non-stationary behavior. We propose to process the kinematic orbits through an epoch-difference (ED) scheme instead of the traditional undifferenced (UD) one in the context of the classical dynamic approach. When differences between neighboring positions are used as input instead of positions themselves, the systematic errors under consideration are largely eliminated, so that an improvement in the resulting gravity field solutions can be expected.

To demonstrate the added value of the proposed ED scheme, two sets of monthly gravity field solutions are produced using the classical dynamic approach from six years (2005–2010) of GRACE kinematic orbits. One is based on the traditional UD scheme, and the other on the ED one. To evaluate the quality of the GPS-based solutions, we use as the reference the latest ITSG-Grace2018 monthly gravity solutions produced at the Institute of Geodesy, Graz University of Technology (Mayer-Gürr et al. 2018). Those solutions are mainly based on the GRACE KBR measurements and therefore are much more accurate than any GPS-based solution.

This paper is organized as follows: Sect. 2 presents in detail the data and methods adopted for gravity field modeling. Results are shown in Sect. 3 and discussed in Sect. 4. Finally, Sect. 5 is left for conclusions.

2 Strategy

2.1 Data and models

Data processing in this study is performed with the Position And Navigation Data Analyst (PANDA) software, which is developed at the GNSS Research Center of Wuhan University and has been widely used in precise orbit determination for both GNSS satellites and low Earth orbiters (Liu and Ge 2003; Shi et al. 2008). Recently, the dynamic approach to gravity field modeling has been implemented in PANDA and successfully applied to produce GRACE monthly gravity field solutions (Guo et al. 2017b, 2018). In this approach, data processing consists of two steps when gravity field is estimated from kinematic orbits. In the first step, a priori dynamic orbits are computed by fitting to the kinematic orbits through numerically integrating the equation of motion defined by the a priori force models (cf. Table 1). During this process, orbits are expressed as truncated Taylor series with respect to the unknown parameters (cf. Table 1), about the computed a priori orbits. The partials with respect to those parameters are obtained by resolving the so-called variational equations. Only non-gravity parameters are adjusted at this step. In the second step, the gravity parameters are included. Based on the partial derivatives with respect to all unknown parameters, daily normal equations (NEQs) are set up for the purpose of the subsequent least-squares adjustment. Arc-specific parameters are then pre-eliminated, and the daily NEQs are accumulated into monthly NEQs. The monthly NEQ matrix is eventually inverted in order to obtain the corrections of the spherical harmonic coefficients (SHCs) with respect to the a priori gravity values.

Table 1 Data and models used for gravity field modeling
Full size table

Details about the data and models used in this study are also described in Table 1. For observing temporal variations of the Earth’s gravity field, the a priori static gravity field model should be as accurate as possible to reduce errors in the computed a priori dynamic orbits. For that purpose, we adopt the static part of the EIGEN-6C4 model, which is compiled from 25 years of satellite laser ranging observations of Lageos, 10 years of GRACE KBR data, 3.5 years of GOCE satellite gravity gradient observations, and terrestrial surface data (Förste et al. 2014). Importantly, old version of the products for the atmosphere and ocean de-aliasing and satellite attitudes have been used in this study, though more recent versions (RL06 and RL03, respectively) are currently available. This is in order to keep consistency with the other GPS-based solutions considered in this study. We compute our gravity field solutions only up to degree and order 60. Such a high maximum degree is definitely beyond the sensitivity range of GPS data. It was deliberately chosen to avoid spectral aliasing of high-frequency signals into low-degree coefficients.

The computed solutions are compared against those produced at the Institute of Geodesy, Graz University of Technology (Zehentner and Mayer-Gürr 2015), which are complete to degree and order 100. We also did a limited number of test computations up to degree 100. The results showed that the differences between the solutions complete to degrees 60 and 100 were in all cases negligible below degree 55. This implies that the chosen resolution has little influence on the long-wavelength gravity field determination from kinematic orbits. We remind that the same finding has also been reported by Baur et al. (2012). This justifies our decision to produce solutions up to a reduced maximum degree (i.e., 60), which allowed us to reduce the time of computations substantially. For the sake of further consistency, we use the orbits produced by Zehentner and Mayer-Gürr (2015) as input. Finally, the arc length is set equal in our data processing scheme to 24 h in order to keep consistency with the parameterization of the exploited kinematic orbits.

2.2 Epoch-difference scheme

The linearized undifferenced observation equation can be expressed as follows:

$$ {\mathbf{y}} = {\mathbf{A}}\left( {{\mathbf{x}}_{0} { + }\Delta {\mathbf{x}}} \right), $$
(1)

where y is the observation vector, A the design matrix, \( {\mathbf{x}}_{0} \) the vector composed of a priori values of parameters, \( \Delta {\mathbf{x}} \) the vector composed of parameter corrections to be estimated (cf. Table 1). By denoting the associated noise variance–covariance matrix as C, we can get the weighted least square (WLS) solution:

$$ \Delta {\mathbf{x}} = \left( {{\mathbf{A}}^{T} {\mathbf{C}}^{ - 1} {\mathbf{A}}} \right)^{ - 1} {\mathbf{A}}^{T} {\mathbf{C}}^{ - 1} {\mathbf{d}}, $$
(2)

where d = y − Ax0 is the ‘observed minus computed’ (O–C) vector. As mentioned above, the observations (the kinematic orbits in this study) may be contaminated by systematic errors. In addition, according to the Hill’s equation (Kaplan 1976), errors in the background force models would also lead to constant and time-varying perturbations in the computed a priori orbits. As soon as these perturbations cannot be explained by the adopted gravity field model, they form another source of non-stationary systematic errors in the data. All these systematic errors will enter into the O–C vector and contaminate the final solutions.

To mitigate the systematic errors, we introduce the epoch-difference matrix M. This matrix consists of blocks, one block per an uninterrupted fragment of the data set within an orbital arc:

$$ \left[ {\begin{array}{*{20}l} 1 \hfill &\quad { - 1} \hfill &\quad 0 \hfill &\quad 0 \hfill &\quad \cdots \hfill \\ 0 \hfill &\quad 1 \hfill &\quad { - 1} \hfill &\quad 0 \hfill &\quad \cdots \hfill \\ 0 \hfill &\quad {} \hfill &\quad 1 \hfill &\quad { - 1} \hfill &\quad \cdots \hfill \\ \vdots \hfill &\quad \vdots \hfill &\quad \vdots \hfill &\quad \vdots \hfill &\quad \vdots \hfill \\ 0 \hfill &\quad 0 \hfill &\quad \cdots \hfill &\quad 1 \hfill &\quad { - 1} \hfill \\ \end{array} } \right]. $$
(3)

Then, the linearized epoch-difference observation equation can be expressed as follows:

$$ {\tilde{\mathbf{d}}} = {\tilde{\mathbf{A}}}\Delta {\mathbf{x}} $$
(4)

where \( {\tilde{\mathbf{d}}} \) is the transformed data vector: \( {\tilde{\mathbf{d}}} = {\mathbf{Md}} \), and \( {\tilde{\mathbf{A}}} \) is the associated design matrix: \( {\tilde{\mathbf{A}}} = {\mathbf{MA}} \). The corresponding WLS solution is:

$$ \Delta {\mathbf{x}} = \left( {{\tilde{\mathbf{A}}}^{T} {\tilde{\mathbf{C}}}^{ - 1} {\tilde{\mathbf{A}}}} \right)^{ - 1} {\tilde{\mathbf{A}}}^{T} {\tilde{\mathbf{C}}}^{ - 1} {\tilde{\mathbf{d}}} $$
(5)

where \( {\tilde{\mathbf{C}}} \) is the error variance–covariance matrix associated with the transformed data vector. As soon as this matrix is close to matrix MCMT, solutions (2) and (5) are nearly equivalent (they would be fully equivalent if the matrix MCMT were invertible). This allowed previous authors to claim that a linear transformation of the data vector must not affect the final results (Ditmar and van Eck van der Sluijs 2004). In practice, however, the exploited error variance–covariance matrices are defined differently. They are either proportional to a unit one (white-noise assumption) or have the Toeplitz structure (color-noise assumption). In both cases, possible non-stationarity of data noise is ignored, since this would make an estimation of the error variance–covariance matrix problematic (if not impossible). As a result, different approaches may not be equivalent in practice. Better results could be expected in the case when non-stationary errors in the data are suppressed. We believe, therefore, that Eq. (5) is likely preferable in the presence of non-stationary systematic errors in the data, particularly if the errors are (piecewise) constant or show slow temporal variations. Indeed, let us assume that the O–C vector suffers from systematic errors and that these errors are partly formed by a bias \( {\varvec{\updelta}} \), which is constant per an uninterrupted data fragment (or per orbital arc). Then, the UD solution, according to Eq. (2), suffers from a bias given by:

$$ \left( {{\mathbf{A}}^{T} {\mathbf{C}}^{ - 1} {\mathbf{A}}} \right)^{ - 1} {\mathbf{A}}^{T} {\mathbf{C}}^{ - 1} {\varvec{\updelta}}. $$
(6)

On the other hand, the ED solution is free from such a bias because \( {\mathbf{M\delta }} = 0 \).

2.3 Data weighting scheme

Generally speaking, errors in kinematic orbits are often temporally correlated. This is mainly due to the float-estimated ambiguity parameters (Jäggi et al. 2011b; Montenbruck et al. 2018). This is the case, among others, for the kinematic orbits under consideration. On the other hand, deficiencies in the background force models lead to correlated errors in the computed dynamic orbits (Ditmar et al. 2012). A traditional way to handle such noise is based on the assumption that it is stationary (e.g., (Klees et al. 2003)). A correlated stationary noise is called colored (or frequency dependent). This implies that the statistically optimal data inversion can be realized by means of frequency-dependent data weighting (FDDW). Recently, the FDDW concept has successfully been applied to GRACE KBR data processing with the classical dynamic approach and has notably reduced noise in the WHU-RL01 monthly gravity solutions (Guo et al. 2018). As far as kinematic orbits are concerned, a successful application of FDDW was shown so far for sampling intervals not less than 30 s (Farahani et al. 2013; Guo et al. 2017a). To fully exploit the gravity information in the 10-s sampled kinematic orbits, we have applied an extended differentiation scheme similar to the one proposed by Baur et al. (2012). To be specific, the sampling points for epoch-differencing are chosen to be 30-s apart, whereas the differentiation moves along the original 10-s sampled orbit track. As such, there are three sets of observation equations, each set being shifted by 10 s. In fact, we also tried a larger differencing interval of 60-s and did some test experiments. The results revealed that the differences between the solutions produced with 30-s and 60-s intervals were negligible.

In this study, we adopt the FDDW variant proposed by Ditmar et al. (2007); we do so in the context of both the UD and ED schemes. To represent the dependence of noise on frequency, we consider noise power spectral density (PSD), which is estimated from the respective postfit observation residuals. For that purpose, we start from a simple assumption of white noise in the observations during the first iteration. Then, the noise PSDs are estimated from the postfit residuals and are used for applying the FDDW scheme in the second iteration. No further iterations are needed due to the accurate background models. The reader is referred to Ditmar et al. (2007) for more details. We would like to stress that the same procedure was applied to UD and ED residuals.

3 Results

3.1 Analysis in the spectral domain

As shown in Table 1, the accelerometer biases are estimated as piecewise constant parameters. It has been shown that accelerometer biases are rather effective to absorb model deficiencies (Guo et al. 2018). To investigate the reactions of the UD and ED schemes to the estimation interval of accelerometer biases, here we consider two intervals: 24 and 1.5 h. We did not try shorter intervals to avoid a risk of an over-parameterization and a signal absorption.

Figure 1 shows the square-root PSD, hereafter denoted as PSD1/2, estimated from the postfit residuals of the GRACE A satellite for a typical month of January 2010 (GRACE B behaves similarly and, therefore, is not addressed in the discussion below). It can be seen that the residuals obtained by different differentiation schemes exhibit clearly different noise patterns. In general, PSD1/2s in the case of the ED scheme are lower than those in the case of UD scheme up to about 8 mHz for both estimation intervals of accelerometer biases. This is not surprising, since the k-th order differentiation in the time domain corresponds to the multiplication of PSD1/2 with ωk in the frequency domain (where ω denotes the angular frequency). In other words, a differentiation reduces the low-frequency noise on the one hand and amplifies the high-frequency noise on the other hand.

Fig. 1
figure 1

PSD1/2s estimated from postfit residuals produced with different differentiation schemes and estimation intervals of accelerometer biases; the left and right columns are for the 24- and 1.5-h estimation intervals, respectively; the top, middle and bottom rows are for the along-track, cross-track and radial components, respectively

Full size image

In the case of the UD scheme, the low-frequency noise is notably reduced when an estimation interval of 1.5 h is applied, implying that the additionally estimated biases absorb errors in observations and background models. On the other hand, the PSD1/2s in the case of the ED scheme are rather insensitive to the choice of the estimation interval when compared to the UD scheme, as can be seen in Fig. 1. This is in line with our expectation that the ED scheme is more immune to slowly varying systematic errors than the UD scheme.

As mentioned above, the KBR-based ITSG-Grace2018 solutions (denoted here as ‘ITSG2018’ for brevity) are chosen as the ‘ground truth’. They consist of three sets of solutions, which have been computed up to degrees 60, 96, and 120, respectively. In this study, we use the solutions complete to degree 60. Therefore, to assess the computed GPS-based gravity field solutions, we subtract from them the ITSG2018 solutions. After that, we calculate the RMS (root mean square) for each SHC time series. Figure 2 displays the results in terms of geoid height errors per degree. It can be seen that the degree-errors in the case of the ED scheme are systematically smaller in the entire degree range than those in the case of the UD scheme for both estimation intervals of accelerometer biases. While a shorter estimation interval (1.5 h) improves the very low degree terms in the case of the UD scheme, the improvements are negligible when the ED scheme is applied. These results are well consistent with the analysis of residuals as shown in Fig. 1. Table 2 further lists the RMS values of cumulative geoid height errors up to degree 5, 10, 20 and 60. It reveals that the cumulative errors up to degree 20 can be reduced by at least 20% when the ED scheme is applied. From the above results, we can conclude that the proposed ED scheme is rather effective in mitigating non-stationary systematic errors in kinematic orbits and force models. It can notably improve the low-degree terms, which are of particular importance for recovery of temporal variations in the Earth’s gravity field. Finally, the obtained results allow us to choose the 1.5-hour estimation interval of accelerometer biases to perform the further calculations, since it can provide better gravity field solutions, particularly when the UD scheme is applied (cf. Table 2).

Fig. 2
figure 2

Geoid height errors per degree for the solutions produced with different differentiation schemes and estimation intervals of accelerometer biases (RMS values in 2010)

Full size image
Table 2 Cumulative geoid height errors (cm) up to degree 5, 10, 20 and 60 for solutions produced with different differentiation schemes and estimation intervals of accelerometer biases (RMS values in 2010)
Full size table

Hereafter, we also compare our GPS-based solutions with those produced at the Institute of Geodesy, Graz University of Technology (Zehentner and Mayer-Gürr 2015). Those solutions were computed to degree and order 100 based on the same kinematic orbits using the short-arc approach (SAA) (Mayer-Gürr 2006). During the solution process, no differentiation scheme was applied. In our study, we truncate those solutions to degree and order 60 and denote them as the ‘UD-SAA’ solutions. Similarly, our solutions produced with the UD and ED schemes are denoted as ‘UD-CDA’ and ‘ED-CDA’ solutions, respectively. To compare the ‘UD-SAA’ solutions with ours, we also subtract from them the ITSG2018 monthly solutions. After that, we calculate the RMS for each SHC time series spanning the same interval as before, i.e., from Jan. 2005 to Dec. 2010. Figure 3 displays the obtained RMS errors in terms of geoid heights per degree for different solutions. Apparently, the ED-CDA solution outperforms the other two in the entire degree range. We also note clear discrepancies at very low degrees (below degree 8) between the UD-SAA and UD-CDA solutions. This is in spite of the fact that both the SAA and CDA solutions are obtained on the basis of the Newton’s second law of motion and, therefore, theoretically must be equivalent. Probably, these discrepancies can be attributed to a different choice of some data processing parameters. A further investigation is needed to shed more light on this difference, but this is beyond the scope of this study. Table 3 lists the cumulative geoid height errors up to degree 5, 10, 20 and 60 for different solutions. Again, the errors at low-degree terms are considerably reduced (by at least 22%) in the case of the ED-CDA solution, when compared to the other two.

Fig. 3
figure 3

Geoid height errors per degree for different gravity field solutions (RMS values in 2005–2010)

Full size image
Table 3 Cumulative geoid height errors (cm) up to degree 5, 10, 20 and 60 of different solutions (RMS values in 2005–2010)
Full size table

3.2 Analysis in the spatial domain

In this section, we compare the ability of the ED and UD schemes to recover temporal gravity field variations and mass anomalies in the spatial domain. Since the unfiltered solutions are dominated by high-frequency noise, post-processing is required. In this study, we use the Gaussian filter (Wahr et al. 1998). To account for the impact of the filter width on the results, we consider two radii: 500 and 750 km. Hereafter, we denote the corresponding solutions as the ‘G500’ and ‘G750’ filtered solutions, respectively. As it was done in the spectral domain, we also compare the GPS-based solutions with the KBR-based ITSG2018 solutions. To keep consistency with the GPS-based solutions, the latter are post-processed with the same Gaussian filters. In addition, the C20 terms, to which GRACE is less sensitive due to the polar orbits (Cheng and Ries 2017), are replaced in all KBR- and GPS-based gravity field solutions by the SLR-derived values (Cheng et al. 2013). As regards the degree 1 terms, which cannot be derived from GRACE data alone, we use the estimates obtained by combining GRACE data and geophysical models as described in Sun et al. (2016, 2017).

To perform the analysis in the spatial domain, we first transform the monthly SHCs to mass anomalies in terms of equivalent water heights (EWHs) on a 1° × 1° grid as explained in Wahr et al. (1998). The correction for the Earth’s oblateness has been applied as proposed by Ditmar (2018). As mass variations (if exist) are primarily linear or/and seasonal, a deterministic model composed of an offset, trend, annual, and semi-annual terms is fitted to the time series of mass anomalies per grid node. Then, a comparison with the ITSG2018 solutions allows us to assess the signals and noise of the mass anomalies inferred from the GPS-based solutions.

To make the following discussion more comprehensive, we also quantify the level of random noise in the mass anomaly time series using an independent approach. The approach is based on the regularization concept and assumes that (1) the target signal is close to a combination of an arbitrary annual periodic function and a long-term linear trend; (2) noise in the time series is white. An estimation of the noise variance \( \sigma_{d}^{2} \) and signal variance \( \sigma_{x}^{2} \) is a part of the regularization procedure. For that purpose, the variance component estimation (VCE) technique (Koch and Kusche 2002) is used. The resulting value of \( \sigma_{d} \) is used as an estimate of standard deviation (SD) of random noise in the considered data time series. For a more extended presentation of this approach, the reader is referred to Ditmar et al. (2018).

3.2.1 Gridded mass anomalies

Figure 4 displays geographical maps of the derived trends and periodic annual signals in mass anomalies inferred from the G500 filtered solutions. It can be seen that after applying the G500 filter, the estimates are still dominated by strong high-frequency noise. In the case of the G750 filtered solutions (Fig. 5), high-frequency noise is largely suppressed and the signal patterns inferred from the GPS-based solutions become rather similar to ITSG2018. Many signals can be clearly seen in the GPS-based solutions in that case, e.g., linear signals over Greenland and West Antarctica; annual signals over the Amazon River basin (marked by the red polygons in the plots in the left column of Fig. 5) and South Africa; and many others.

Fig. 4
figure 4

Linear trends and periodic annual signals in mass anomalies in terms of equivalent water heights, inferred from different G500 filtered solutions. From top to bottom, the plots display the results based on the ITSG2018, UD-SAA, UD-CDA, and ED-CDA solutions, respectively

Full size image
Fig. 5
figure 5

Same as in Fig. 4, but for the G750 filtered solutions. Red polygons in the plots in the left column denote the Amazon River basin (marking it in the middle and right columns will degrade the visibility of the annual signals)

Full size image

We have also computed the RMSs over the time series of gridded mass anomaly differences between the GPS-based solutions and ITSG2018. The geographical distribution of the RMS differences is displayed in Fig. 6. It can be observed that the RMS differences in the case of the ED-CDA solutions are systematically smaller than those for the other two solutions. A further calculation reveals that the weighted mean of the gridded RMS differences (weighted by the cosine of latitude) in the case of the ED-CDA solutions is reduced by 19–23% (depending on the filter width), when compared to the UD-CDA solutions (cf. Table 4), and by 15–23% when compared to the UD-SAA solutions. These results demonstrate that the mass anomalies inferred from the ED-CDA solutions are more consistent with ITSG2018 and, therefore, are more accurate, as compared to the other two solutions.

Fig. 6
figure 6

Geographical distribution of the gridded RMS differences with respect to ITSG2018 for different GPS-based solutions. The top, middle, and bottom rows show the results for the UD-SAA, UD-CDA, and ED-CDA solutions, respectively

Full size image
Table 4 Weighted mean of the gridded RMS differences of mass anomalies (cm) with respect to consistently filtered ITSG2018 solutions, for different GPS-based solutions
Full size table

3.2.2 Regional mass anomalies

To compare the performance of the ED and UD schemes in the context of mean regional signals, we select the Amazon River basin (~ 6.22 × 106 km2) and Greenland (~ 2.26 × 106 km2) as the target regions. Again, Gaussian filters of 500-km and 750-km half-widths are applied.

As shown in Fig. 5, the Amazon River basin exhibits significant annual variations, which are mainly of hydrological origin (Chen et al. 2010; Xavier et al. 2010). As regards Greenland, it is known that mass anomalies there show seasonal variations superimposed by a strong negative long-term trend (Ran et al. 2018; Schrama et al. 2014; Shepherd et al. 2012; Siemes et al. 2013; Velicogna and Wahr 2013). Therefore, our selection of regions allows us to analyze the performance of the considered schemes over regions exhibiting different temporal behaviors of mass anomalies.

Figure 7 displays the time series of mean mass anomalies over the Amazon River basin and Greenland inferred from different solutions. To evaluate the consistency between the GPS-based solutions and ITSG2018, we calculate the correlation coefficients and RMS differences between the mass anomaly time series. The results confirm that the GPS-based solutions obtained with the ED scheme are more consistent with ITSG2018, as evidenced by larger correlation coefficients and smaller RMS differences (cf. Tables 5 and 6).

Fig. 7
figure 7

Time series of mean mass anomalies over the Greenland (top) and Amazon River basin (bottom) inferred from the G500 (left) and G750 (right) filtered solutions

Full size image
Table 5 Correlation coefficients between (1) mass anomaly time series over the Amazon River basin and Greenland inferred from the GPS-based solutions and (2) ITSG2018
Full size table
Table 6 RMS differences (cm) between (1) mass anomaly time series over the Amazon River basin and Greenland inferred from the GPS-based solutions and (2) ITSG2018
Full size table

Table 7 further lists the annual amplitudes and phases, as well as the associated formal errors over the Amazon River basin, inferred from different solutions. One can see that the GPS-based solutions (after applying a Gaussian filter) can reproduce the signals as ITSG2018 does, particularly when the ED scheme is applied. Regarding the mass anomaly trends over Greenland, all the GPS-based solutions are also consistent with ITSG2018, and the differences stay within the 2-σ error margins (Table 8). Remarkably, the signals inferred from the ED-CDA solutions are systematically closer to the ITSG2018 ones over both regions. This prompts not only that the ED-CDA solutions are more accurate but also that their higher accuracy is not compromised by an intrinsic signal damping.

Table 7 Annual amplitudes (cm) and phases (deg) and the associated formal errors of mass anomalies over the Amazon River basin derived from different solutions
Full size table
Table 8 Mass anomaly trends and the associated formal errors (Gt/yr) over Greenland inferred from different solutions. The time interval under consideration is 2005–2010
Full size table

In addition, one may note that the signal powers obtained over the two target regions are notably smaller than those reported by other studies. This is primarily due to the signal damping triggered by the truncated spherical harmonic expansion and Gaussian filtering (Chen et al. 2015). To achieve more realistic estimates of the signals, one has to correct for that signal damping. For that purpose, we up-scale the obtained estimates (Chen et al. 2007; Velicogna and Wahr 2006). In order to determine the scaling factors, we use as the truth the signals inferred from the GRACE mascon solutions. For that purpose, we adopt those generated by the Center of Space Research of the University of Texas at Austin, which are released in the form of SHCs up to degree and order 720 (Save et al. 2016). For the sake of consistency, we replace the C20 terms and restore the degree 1 terms in the mascon solutions as done in other solutions. The signal inferred from these mascon solutions is regarded as the true signal (St). Then, we truncate the mascon solutions to degree and order 60 and apply the Gaussian filters with radii of 500 and 750 km to them. The signal inferred from them is denoted as the filtered signal (Sf). After that, the scaling factor k is determined as: k = St/Sf. We consider the annual amplitudes to estimate the scaling factors for the Amazon River basin and the trends to estimate the scaling factors for Greenland. Finally, we use the scaling factors to up-scale the estimated annual amplitudes and trends over the Amazon River basin and Greenland, respectively.

The scaling factors over the Amazon River basin are determined to be 1.150 and 1.333 for the G500 and G750 solutions, respectively. The up-scaled annual amplitudes (Table 7) agree much better with other studies than the original ones. For example, Zhou et al. (2019) estimated the annual amplitude in 2005–2010 over the Amazon River basin as 19.5 ± 1.5 cm. Interestingly, one may note that all the annual amplitudes inferred from the GPS-based solutions are somewhat larger than those from the KBR-based solutions. A similar observation has also been reported in a recent study of Meyer et al. (2019). There, they reveal that the annual signals inferred from the SWARM GPS-based solutions are larger than the KBR-based estimates over the Amazon River basin. Specific reasons of this effect are not clear yet and require a further investigation.

Regarding the Greenland, the scaling factors are determined to be 2.077 and 2.621 for the G500 and G750 filtered solutions, respectively. The up-scaled trends are listed in Table 8. One can see that the scaled values agree well with those from literature. For instance, Shepherd et al. (2012) estimated the mean rate of Greenland Ice Sheet mass losses in 2005–2010 as − 263 ± 30 Gt/yr.

We also tried to estimate the scaling factors over the two target regions on a monthly basis. The results demonstrated that up-scaling in that way didn’t change our conclusions regarding the relative performance of the GPS-based solutions in general. However, we found that the estimated monthly scaling factors suffered from occasional large jumps. They occurred when the mean mass anomaly in a given region approached zero. Those jumps degraded the signal estimates in the months when they occurred, so that we ultimately decided to refrain from such an approach.

Noise SDs of the mass anomaly time series estimated with the VCE procedure are also considerably reduced when the ED scheme is applied (cf. Table 9). The observed error reduction is even more substantial than the one derived from the differences between GPS-based and ITSG2018 solutions (Table 6). Most probably, this can be explained by the presence of systematic errors in the compared time series, which cannot be taken into account by VCE.

Table 9 VCE-based noise standard deviations (cm) of the mass anomaly time series over the Amazon River basin and Greenland inferred from the GPS-based solutions
Full size table

All these results make us conclude that the proposed ED scheme is able to provide better regional mass anomaly estimates, when compared to the UD scheme.

4 Discussion

Since the ED scheme is applied to the residual quantity, i.e., the O–C vector, one may concern whether it may bias the estimates toward the a priori gravity field model. In our opinion, this is unlikely to happen. On the one hand, the ED scheme is applied not only to the O–C vector, but also to all the columns of the design matrix. As explained in Sect. 2.2, the ED scheme is even nearly equivalent to the UD one, provided that the error variance–covariance matrices are defined consistently. On the other hand, the adopted a priori gravity field model, i.e., EIGEN-6C4, only includes the static part of the Earth’s gravity field. If the ED scheme biased the estimates toward the a priori information, it would damp the temporal signals. However, this is not the case. In fact, the signals inferred from the ED-CDA solutions are closer to ITSG2018, showing no visible damping. The recovered mass loss trend over Greenland is even slightly larger than that inferred from the UD-CDA solutions (cf. Table 8).

To further clarify this issue, we have repeated the data processing, using an alternative static gravity field model, i.e., AIUB-CHAMP03S, as the a priori one. That model was compiled to degree and order 100 from 8 years of GPS data collected by the CHAMP mission from 2002 to 2009. As such, AIUB-CHAMP03S is definitely of a lower quality than EIGEN-6C4. To investigate the possible impacts of changing the a priori model, we produce the monthly solutions with both the UD and ED schemes over the period from January 2005 to December 2010. Then, we subtract from them their counterparts produced based on EIGEN-6C4 and calculate the RMS over the difference time series for each SHC. Figure 8a displays the results in terms of geoid height differences per degree (thicker lines), together with the degree-errors of the monthly solutions with respect to ITSG2018 (thinner lines). One can see that the degree-differences are 1–2 orders of magnitude smaller than the respective degree-errors for both monthly solutions [increased differences near degree 60 are mainly due to spectral aliasing caused by truncating the solutions to degree 60 (Baur et al. 2012)]. This indicates that a CHAMP-only static gravity field model is also qualified for temporal gravity field recovery from the GRACE GPS data and the impact of using it instead of a more accurate model is negligible for both data processing schemes. Remarkably, the degree-differences between the two sets of ED-CDA solutions are even smaller than those between the two sets of UD-CDA solutions in the entire degree range. This demonstrates that the ED scheme is probably somewhat more robust against changing the a priori static gravity field model, when compared to the UD scheme.

Fig. 8
figure 8

a Geoid height errors per degree (thinner lines) for the monthly solutions produced based on AIUB-CHAMP03S, as well as the degree-differences (thicker lines) with respect to their counterparts produced based on EIGEN-6C4 (RMS values in 2005–2010); the errors are produced as differences with respect to ITSG2018 model. b Geoid height errors per degree of the computed mean static gravity field solutions with respect to EIGEN-6C4

Full size image

Furthermore, we have also computed the mean static gravity field solutions over 2005–2010 using both the UD and ED schemes. Figure 8b shows the geoid height errors per degree for the two computed solutions with respect to EIGEN-6C4. In general, the ED solution again outperforms the UD one in the entire degree range. The cumulative errors up to degree 20 are reduced by 20% (from 1.45 cm to 1.16 cm) when the ED scheme is applied. These results are fully consistent with those obtained for monthly solutions (cf. Table 3).

To further demonstrate the superior performance of the ED scheme, we have also done some simulations. In a first step, we compute the dynamic orbits for the GRACE A satellite by integrating the equation of motion defined by the a priori force models as described in Table 1 with three exceptions. First, we use the latest version of the AOD products, i.e., AOD1B RL06, which are truncated to degree and order 100. Second, the accelerometer observations, which record the non-gravitational accelerations, don’t allow a direct use due to unknown scales and biases. Therefore, we calculate the non-gravitational accelerations using empirical models, which account for atmospheric drag, solar radiation, and Earth albedo. All the computed non-gravitational accelerations are transformed into the science reference frame using the satellite attitude data. As such, they can be used in the same way as the accelerometer observations. Third, we also incorporate the temporal gravity field signals, which are obtained by subtracting EIGEN-6C4 from the ITSG2018 monthly solutions. The orbits obtained in the first step are taken as noise-free observations in the subsequent temporal gravity field recovery. For that purpose, we set up two scenarios. In the first one, the nominal force models are error-free, i.e., they are identical to those used to produce the dynamic orbits. In the second one, the nominal force models suffer from errors. To be specific, the static gravity field model EIGEN-6C4 is replaced with a GPS-only model, i.e., AIUB-CHAMP03S. In addition, we further replace the ocean tide model EOT11a with FES2004 (80 × 80) (Lyard et al. 2006) and replace AOD1B RL06 with AOD1B RL05. For each scenario, the monthly solutions for January 2010 are produced using both the UD and ED schemes in exactly the same way as they have been produced from real data. Then, we compare the obtained solutions to the true model, i.e., ITSG2018. Figure 9a displays the results in terms of geoid height errors per degree for the first scenario. One can see that the errors of both solutions are orders of magnitude smaller than the signal. They likely represent a limited numerical accuracy of the computations. The differences between the UD and ED solutions are negligible: within 1% of their errors. Thus, the first scenario demonstrates that both UD and ED schemes do not introduce any bias into the estimates. Figure 9b shows the results for the second scenario. It can be observed that the degree-errors of the ED solution are generally smaller than those of the UD solution. This means that the ED solution is more accurate in the presence of force model errors.

Fig. 9
figure 9

Geoid height errors per degree of different solutions for a the first simulation scenario where the data are error-free and b the second scenario where the data suffer from errors in the background force models. The red and blue lines are for the UD and ED solutions, respectively. The green line denotes the temporal gravity field signals represented by ITSG2018

Full size image

To conclude, we hold the view that the ED scheme can indeed reduce the non-stationary systematic errors originated from kinematic orbits and background models. Thus, the gravity field recovery improves. This is in line with the fact that differencing is often utilized as an effective tool for reducing non-stationary errors in time series analysis (Box et al. 2015).

5 Conclusions

In this study, we have proposed to apply an epoch-difference scheme (ED) in the context of the classical dynamic approach to the Earth’s gravity field recovery from kinematic orbit data. We have shown that this leads to a substantial improvement of the estimated gravity field model parameters. We explain this by the presence of non-stationary systematic errors in kinematic orbits and background force models. We have shown that the ED scheme allows this noise to be largely suppressed, particularly when it concerns the disturbances that show slow temporal variations. As a result, the data noise becomes more stationary, so that a well-known frequency-dependent data weighting is capable of handling it properly.

The observed increase in the quality of data processing is particularly relevant in the estimation of low-degree terms, which contain most of the time-varying gravity signals. To demonstrate the added value of the proposed scheme, two sets of monthly gravity field solutions have been produced based on six years of kinematic orbits of the GRACE satellites using both the ED and undifferenced (UD) scheme. Furthermore, gravity field solutions produced with the UD scheme in the frame of the short-arc approach have also been considered (Zehentner and Mayer-Gürr 2015). As a reference, we have used the latest ITSG-Grace2018 monthly gravity solutions produced at the Institute of Geodesy, Graz University of Technology. A comparison in the spectral domain shows that the low-degree terms can be notably improved when the ED scheme is applied, with cumulative errors up to degree 20 being reduced by at least 20%. Analysis in the spatial domain also shows that gravity field solutions obtained with the ED scheme are more consistent with the ITSG-Grace2018 solutions. This is evidenced by larger correlation coefficients and smaller RMS differences when time series of mean mass anomalies over the Amazon River basin and Greenland are considered. Importantly, parameters of the estimated signals (more specifically, the amplitude of seasonal variations in the Amazon River basin and the linear trend in Greenland) show the best agreement with the ITSG-Grace2018-based estimates when the ED scheme is applied. This implies that the improved accuracy of the ED scheme is not compromised by a signal damping. We conclude that the proposed ED scheme is preferable for time-varying gravity field modeling, as compared to the traditional UD scheme.

Our findings may facilitate, among others, bridging the gap between GRACE and GFO missions by using hl-SST GPS data from other satellites.