A New Approach to Earth’s Gravity Field Modeling Using GPS-Derived Kinematic Orbits and Baselines

Earth’s gravity field recovery from GPS observations collected by low earth orbiting (LEO) satellites is a well-established technique, and kinematic orbits are commonly used for that purpose. Nowadays, more and more satellites are flying in close formations. The GPS-derived kinematic baselines between them can reach millimeter precision, which is more precise than the centimeter-level kinematic orbits. Thus, it has long been expected that the more precise kinematic baselines can deliver better gravity field solutions. However, this expectation has not been met yet in practice. In this study, we propose a new approach to gravity field modeling, in which kinematic orbits of the reference satellite and baseline vectors between the reference satellite and its accompanying satellite are jointly inverted. To validate the added value, data from the Gravity Recovery and Climate Experiment (GRACE) satellite mission are used. We derive kinematic orbits and inter-satellite baselines of the twin GRACE satellites from the GPS data collected in the year of 2010. Then two sets of monthly gravity field solutions up to degree and order 60 are produced. One is derived from kinematic orbits of the twin GRACE satellites (‘orbit approach’). The other is derived from kinematic orbits of GRACE A and baseline vectors between GRACE A and B (‘baseline approach’). Analysis of observation postfit residuals shows that noise in the kinematic baselines is notably lower than the kinematic orbits by 50, 47 and 43% for the along-track, cross-track and radial components, respectively. Regarding the gravity field solutions, analysis in the spectral domain shows that noise of the gravity field solutions beyond degree 10 can be significantly reduced when the baseline approach is applied, with cumulative errors up to degree 60 being reduced by 34%, when compared to the orbit approach. In the spatial domain, the recovered mass changes with the baseline approach are more consistent with those inferred from the K-Band Ranging based solutions. Our results demonstrate that the proposed baseline approach is able to provide better gravity field solutions than the orbit approach. The findings may facilitate, among others, bridging the gap between GRACE and GRACE Follow-On satellite mission.


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. It has been mostly observed with the ultra-precise K-Band ranging (KBR) measurements from the Gravity Recovery and Climate Experiment (GRACE) mission [1] and its successor, i.e., the GRACE Follow-On (GFO) mission. However, there exists a gap between GRACE and GFO. 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 [2][3][4][5][6][7]. For that purpose, GPS-derived kinematic orbits are usually taken as pseudo-observations in gravity field modeling. Since GPS data are about three 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 [8][9][10]. As a result, the recovered gravity field solutions usually suffer from systematic errors, particularly at low degrees [11]. Therefore, efforts are still ongoing to further improve the GPS-based gravity field solutions. Among others, GPS-derived kinematic baselines between formation flying satellites have drawn much attention, due to their higher precision than the kinematic orbits.
To exploit the more precise space baselines, two approaches have been proposed in the context of the celestial mechanics approach to gravity field modeling [12]: The 'observation equation approach' and the 'GRACE-type approach' [13]. In the former, individual observation equations for both satellites are first established based on their positions; then the two individual observation equations are differenced to form the baseline observation equations; finally, the normal equations (NEQs) are established from the baseline observation equations. Since only the baselines are used, orbit parameters for one of the two satellites have to be fixed to the a priori values in this approach. Unfortunately, this operation will seriously degrade the gravity field solutions at very low degrees [13,14]. In the GRACE-type approach, individual observation and normal equations for both satellites are established based on their kinematic orbits in a first step; in a second step, the observation and normal equations for baseline vectors or baseline lengths are established based on kinematic baselines; the NEQs from the first and second step are combined to form the final NEQs in the third step. Unlike the observation equation approach, orbit parameters for both satellites are estimated in the GRACE-type approach. Recent studies based on this approach all use baseline lengths to perform gravity field recovery, the improvements are demonstrated to be insignificant [13] or negligible [15], when compared to the gravity field solutions based solely on kinematic orbits.
To make full use of the precise kinematic baseline vectors, we propose a new approach in this study, where kinematic orbits of the reference satellite and kinematic baseline vectors between the reference satellite and its accompanying satellite are taken as observations during gravity field recovery. Hereafter, we denote the new approach as the 'baseline approach'. Accordingly, the approach based solely on kinematic orbits is denoted as the 'orbit approach'. To demonstrate the added value of the new approach, data collected by the GRACE mission in 2010 are used. We first compute the kinematic orbits and baselines of the twin GRACE satellites. Then two sets of monthly gravity field solutions are produced: One is based on kinematic orbits of the twin satellites with the orbit approach (hereafter denoted as the 'orbit solution' where applicable), the other is based on kinematic orbits of GRACE A and kinematic baseline vectors between GRACE A and B with the baseline approach (hereafter denoted as the 'baseline solution' where applicable). To evaluate the quality of the GPS-based solutions, we use as the reference the WHU RL01 monthly gravity field solutions produced at the GNSS Research Center of Wuhan University [16,17]. Those solutions are mainly based on the GRACE KBR measurements and therefore are much more accurate than any GPS-based solutions.
This paper is organized as follows: Section 2 presents in detail the data and methods adopted for gravity field modeling. Results are shown in Section 3, and discussed in Section 4. Finally, Section 5 is left for conclusions.

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 [18,19]. Recently, the dynamic approach to gravity field modeling has been implemented in PANDA and successfully applied to produce GRACE monthly gravity field solutions [16,17,20,21]. In this approach, data processing consists of two steps when the 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 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.
Details about the data and models used in this study are also described in Table 1. The kinematic orbits of the twin GRACE satellites used in this study are produced in house based on a zero-difference data processing scheme. To achieve the best accuracy, both the single-difference ambiguities between GPS satellites at each GRACE satellite and double-difference ambiguities between the GRACE satellites and GPS satellites are resolved. Accuracy of the kinematic orbits obtained in this way is about 1.5 cm as inferred from satellite laser ranging residuals. More details can be found in Reference [22]. As to the kinematic baselines, they are produced based on a single-difference scheme data processing as described in Reference [23]. In this scheme, the positions of the reference satellite (which is GRACE A in this study), are not estimable, and have to be fixed to the a priori values. This requires accuracy of the a priori orbits to be better than 10 cm for millimeter precision baseline determination [24]. For that purpose, we exploit the abovementioned ambiguity-fixed kinematic orbits, which are fully qualified with an accuracy better than 2 cm. Accuracy of the kinematic baselines obtained in this way is demonstrated to be slightly better than 4 mm, as inferred from KBR residuals, which is consistent with other studies [13,25]. During the solution process, the accelerometer observations are used to model the non-gravitational forces. For each component in the accelerometer frame, the scale factors are estimated on a monthly basis, and the biases are modeled with a three-order polynomial and the polynomial coefficients are estimated on a daily basis to account for thermal variations. Finally, the arc length is set equal in our data processing scheme to 24 h.

The Orbit and Baseline Approach
In the orbit approach, only kinematic orbits of the twin GRACE satellites are taken as observations. The linearized observation equations are expressed as follows: where y is the satellite position vector, which is extracted from the kinematic orbits, and y 0 is its computed counterpart; x ng denotes the non-gravity parameter vector, which consists of satellite initial state vector, accelerometer scale and bias parameters (cf. Table 1); x g denotes the gravity parameters in the form of SHCs up to degree and order 60; d ng and d g denote the design matrices with respect to the non-gravity and gravity parameters, respectively; Finally, ε denotes the noise which includes noise in both the position and its computed counterpart. Data processing consists of three steps. First, individual observation and normal equations for the two satellites are established based on kinematic orbits. Second, the two individual NEQs are combined, and arc-specific parameters are pre-eliminated from the daily NEQs. Third, the daily NEQs are accumulated into monthly NEQs, which are eventually solved to obtain monthly gravity field solutions.
In the baseline approach, the observations consist of kinematic orbits of GRACE A and kinematic baseline vectors between GRACE A and B. The linearized observation equations are expressed as follows: where y AB and y 0 AB are the baseline vector and its computed counterpart, respectively; ε y AB denote the noise which consists of noise in both the baseline and its counterpart; other symbols are same to those in Equation (1). Data processing includes the following steps: First, kinematic orbits of GRACE B are derived from the kinematic orbits of GRACE A and the kinematic baseline vectors between GRACE A and B, i.e., y B = y A + y AB . Second, individual observation equations for both satellites are established based on the kinematic orbits. Third, the baseline observation equations are formed by differencing the two orbit observation equations, i.e., y AB = y B − y A . Fourth, individual NEQs are established from observation equations of kinematic orbits of GRACE A and observation equations of kinematic baseline vectors between GRACE A and B, respectively. Fifth, the individual NEQs are combined and arc-parameters are pre-eliminated from the daily NEQs. Finally, the daily NEQs are accumulated into monthly NEQs, which are eventually inverted to obtain monthly gravity field solutions.
As described above, six observations per epoch are used in both the orbit and baseline approach. While observations consist of the position vectors of GRACE A and B in the orbit approach, they are composed of the position vector of GRACE A and baseline vector between GRACE A and B in the baseline approach. It should be noted that the same kinematic orbits of GRACE A are used, and the same parameters are set up in both approaches (cf. Table 1). In the following calculations, we further replace the kinematic orbits of GRACE B with those constructed from kinematic orbits of GRACE A and kinematic baseline vectors between GRACE A and B in the context of the orbit approach. Our test experiments show that impacts of this replacement on the gravity field solutions are minor. Now, one may consider that the baseline approach is equivalent to the orbit approach, since the baseline vectors are just linear combinations of the two kinematic orbits used in the orbit approach and this will not change the final solutions in frame of the weighted least square adjustment if the stochastic models for different observations are properly accounted for [32]. However, this is not the case when the observations are contaminated by systematic errors, which has long been a problem for time-varying gravity field recovery from spaceborne GPS data. In that case, the baseline approach has two advantages at least over the orbit approach: First, the baseline vectors are free of common-mode errors from the GPS satellites, due to the adopted single-difference data processing scheme when deriving them from the GPS measurements. Thus, the baselines are more precise than kinematic Remote Sens. 2019, 11, 1728 5 of 14 orbits calculated with a zero-difference scheme. Second, according to the Hill's equation [33], errors in the initial state vectors and background force models would lead to constant or slowly varying perturbations in the computed a priori orbits which are used to establish the observation equations of kinematic orbits. While these systematic errors will enter into the O−C (observed minus computed) vectors and contaminate the final gravity field solutions in the orbit approach, they can be mitigated when forming the baseline observation equations through differencing the two orbit observation equations. This is particularly true for the GRACE mission, because errors for the two satellites are largely common, since the two spacecraft are identical and not far separated (about 220 km) in a coplanar orbit. Therefore, one may expect that the baseline approach will deliver better gravity field solutions than the orbit approach, which will also be demonstrated later in this article.

Data Weighting Scheme
Generally speaking, errors in kinematic orbits are often temporally correlated, although ambiguity-fixing can largely reduce the correlations [34]. On the other hand, deficiencies in the background force models also lead to correlated errors in the computed dynamic orbits [35]. As a result, noise in the residuals is usually frequency-dependent (i.e., colored). In this study, we adopt the frequency-dependent data weighting (FDDW) concept proposed by Reference [36] to account for the colored noise during the solution process. 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 [17]. To represent the dependence of noise on frequency, we consider noise power spectral density (PSD), which is estimated from postfit observation residuals. For that purpose, we start from a simple assumption of white noise in the observations during the first computation. Then the noise PSDs are estimated from the postfit residuals and are used for applying the FDDW scheme in the second computation. No further iteration is needed, due to the accurate background models. The reader is referred to Reference [36] for more details. Figure 1 shows the square-root PSD, hereafter denoted as PSD 1/2 , estimated from the postfit residuals of kinematic orbits of GRACE A and kinematic baselines between GRACE A and B for a typical month of January 2010. Since the PSD 1/2 s for kinematic orbits of both satellites and in both approaches show similar patterns, only those for GRACE A are displayed here. It can be observed that the PSD 1/2 s for both kinematic orbits and baselines exhibit clear frequency-dependent behaviors, which necessitates the FDDW scheme adopted in this study. A comparison between the two sets of PSD 1/2 s reveals that PSD 1/2 s for the kinematic baselines are generally lower than those for the kinematic orbits. On average, PSD 1/2 s for the kinematic baselines are reduced by 50, 47 and 44% at the along-track, cross-track and radial components, respectively, when compared to those for the kinematic orbits. Obviously, the kinematic baselines are of higher precision than the kinematic orbits.

Analysis in the Spectral Domain
As mentioned above, the KBR-based WHU RL01 monthly gravity field solutions are chosen as the ground truth. In fact, we also tried the official GRACE monthly solutions for that purpose. The results revealed that the differences were negligible. Therefore, to assess the computed GPS-based gravity field solutions, we subtract from them the WHU RL01 solutions and the obtained residual SHCs are interpreted as the 'true' errors of the GPS-based solutions. To facilitate analysis in the spectral domain, we calculate the degree errors (DEG) and cumulative errors (CUM) of different solutions in terms of geoid heights with the following formulas: where R is the reference radius (6,378,137.46 m), ∆C lm , ∆S lm are the residual SHCs with respect to WHU RL01. Figure 2 displays the degree errors of different GPS-based solutions and degree signals represented by WHU RL01 in January 2010. It can be seen that the degree errors of the baseline solutions are smaller than those of the orbit solutions, except for relatively low degree terms. To make the discussion more comprehensive, we further calculate the RMS (root mean square) for each SHC time series, as well as for the associated formal errors in 2010. The results are displayed in Figure 3, which in general are similar to those displayed in Figure 2. One can see that degree errors beyond degree 10 can be significantly reduced when the baseline approach is applied. One can also observe that the true errors clearly depart from the formal errors below degree 15 in both cases. This may likely be attributed to the presence of systematic errors in the observations, e.g., mismodeling of GPS antenna phase center, high-order ionosphere-induced errors as mentioned before. Table 2 further lists the cumulative geoid height errors up to degree 10, 20 and 60. It reveals that the cumulative errors up to degree 20 and 60 are reduced by 15 and 34%, respectively, when the baseline approach is applied. From these results, we can conclude that the proposed baseline approach can improve the gravity field beyond degree 10. As mentioned above, the KBR-based WHU RL01 monthly gravity field solutions are chosen as the ground truth. In fact, we also tried the official GRACE monthly solutions for that purpose. The results revealed that the differences were negligible. Therefore, to assess the computed GPS-based gravity field solutions, we subtract from them the WHU RL01 solutions and the obtained residual SHCs are interpreted as the 'true' errors of the GPS-based solutions. To facilitate analysis in the GPS antenna phase center, high-order ionosphere-induced errors as mentioned before. Table 2 further lists the cumulative geoid height errors up to degree 10, 20 and 60. It reveals that the cumulative errors up to degree 20 and 60 are reduced by 15 and 34%, respectively, when the baseline approach is applied. From these results, we can conclude that the proposed baseline approach can improve the gravity field beyond degree 10.   observe that the true errors clearly depart from the formal errors below degree 15 in both cases. This may likely be attributed to the presence of systematic errors in the observations, e.g., mismodeling of GPS antenna phase center, high-order ionosphere-induced errors as mentioned before. Table 2 further lists the cumulative geoid height errors up to degree 10, 20 and 60. It reveals that the cumulative errors up to degree 20 and 60 are reduced by 15 and 34%, respectively, when the baseline approach is applied. From these results, we can conclude that the proposed baseline approach can improve the gravity field beyond degree 10.

Analysis in the Spatial Domain
In this section, we compare the ability of the orbit and baseline approaches 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 [37] for that purpose. To account for the impact of the filter width on the results, we consider two radii-750 and 1000 km. Hereafter, we denote the corresponding solutions as the 'G750' and 'G1000' filtered solutions, respectively. As it was done in the spectral domain, we also compare the GPS-based solutions with the KBR-based WHU RL01 solutions. To keep consistency with the GPS-based solutions, the latter is also post-processed with the same Gaussian filters. In addition, the C 20 terms, to which GRACE is less sensitive, due to the polar orbits [38], are replaced in all gravity field solutions with the SLR-derived values [39]. As regards the degree 1 terms, which cannot be derived from GRACE data Remote Sens. 2019, 11, 1728 8 of 14 alone, we use the estimates obtained by combining GRACE data and geophysical models as described in References [40,41].
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 Reference [37]. The correction for the Earth's oblateness has been applied as proposed by Ditmar [42]. 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 WHU RL01 solutions allows us to assess the signals and noise of the mass anomalies inferred from the GPS-based solutions. Figure 4 displays geographical maps of the derived periodic annual signals in mass anomalies inferred from the G750 filtered solutions. It can be seen that after applying the G750 filter, the estimates are still dominated by strong high-frequency noise. In the case of the G1000 filtered solutions (Figure 5), high-frequency noise is largely suppressed, and the signal patterns inferred from the GPS-based solutions become similar to WHU RL01. Among others, annual signals over the Amazon river basin in South America can be clearly seen in the GPS-based solutions.      Figure 4, but for the G1000 filtered solutions.

Gridded Mass Anomalies
We have also computed the RMS differences of gridded mass anomalies between the GPS-based solutions and WHU RL01. The geographical distribution of the RMS differences is displayed in Figure 6. It can be observed that the RMS differences in the case of the baseline solutions are generally smaller than those in case of the orbit solutions when the G750 filter is applied. A further calculation reveals that the weighted mean of the gridded RMS differences (weighted by the cosine of latitude) in the case of the baseline solutions are reduced by 14 and 6% for the G750 and G1000 filtered solutions, respectively, when compared to those in case of the orbit solutions (cf. Table 3). These results demonstrate that the mass anomalies inferred from the baseline solutions are more consistent with WHU RL01 and, therefore, are more accurate, as compared to the orbit solutions. Table 3. Weighted mean of the gridded RMS differences of mass anomalies (cm) with respect to consistently filtered WHU RL01 solutions, for different GPS-based solutions.

Regional Mass Anomalies
To compare the performance of the orbit and baseline approaches in the context of mean regional signals, we select the Amazon river basin as the target regions. As shown in Figure 5, the Amazon river basin exhibits significant annual variations, which are mainly of hydrological origin [43,44]. Figure 7 displays the time series of mean mass anomalies over the Amazon river basin inferred from different solutions in 2010. To evaluate the consistency between the GPS-based solutions and WHU RL01, 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 baseline approach are more consistent with WHU RL01, as evidenced by larger correlation coefficients and smaller RMS differences (cf.

Regional Mass Anomalies
To compare the performance of the orbit and baseline approaches in the context of mean regional signals, we select the Amazon river basin as the target regions. As shown in Figure 5, the Amazon river basin exhibits significant annual variations, which are mainly of hydrological origin [43,44]. Figure 7 displays the time series of mean mass anomalies over the Amazon river basin inferred from different solutions in 2010. To evaluate the consistency between the GPS-based solutions and WHU RL01, 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 baseline approach are more consistent with WHU RL01, as evidenced by larger correlation coefficients and smaller RMS differences (cf. Table 4). Table 5 further lists the annual amplitudes and phases and 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 WHU RL01 does, particularly when the baseline approach is applied. The signals inferred from the baseline solutions are closer to those inferred from the WHU RL01 solutions when compared to the orbit solutions.
All these results lead us to conclude that the proposed baseline approach is able to provide better regional mass anomaly estimates, when compared to the orbit approach.     Table 5 further lists the annual amplitudes and phases and 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 WHU RL01 does, particularly when the baseline approach is applied. The signals inferred from the baseline solutions are closer to those inferred from the WHU RL01 solutions when compared to the orbit solutions. All these results lead us to conclude that the proposed baseline approach is able to provide better regional mass anomaly estimates, when compared to the orbit approach.

Discussion
As shown in the text, with the proposed baseline approach we have achieved more promising results than previous studies which are based on the GRACE-type approach [13,15], when exploiting GPS-derived kinematic baselines for Earth's gravity field recovery. While improvement of the gravity field solutions was insignificant or even negligible in previous investigations, it is shown to be significant in our case with cumulative errors up to degree 60 being reduced by 34%. We attribute this improvement mainly to the fact that we have made full use of the precise three-dimensional baseline vectors in our approach, whereas only the baseline lengths were used in previous studies.
Our results also show that the gravity field solutions still cannot benefit from the proposed baseline approach at the very low degree terms (below degree 10) at present. Because these low degree terms are very prone to systematic errors in the observations, further efforts are still needed to identify and reduce the possible systematic errors.
Another issue is that we have made use of the high precision accelerometer observations to account for the non-gravitational forces acting on the GRACE satellites, which may not be available for other formation flying satellites. Indeed, high precision accelerometer observations are indispensable to achieve the best possible gravity field solutions from hl-SST data. In fact, it is hardly possible to obtain realistic gravity field solutions if the non-gravitational forces cannot be adequately modeled either with accelerometer observations or with physical models. Thus, this is the prerequisite for both the orbit and baseline approach. Because the kinematic baselines are more precise than kinematic orbits (Figure 1), it is reasonable that the baseline approach can deliver better gravity field solutions, at least for formation flying satellites equipped with high precision accelerometers.

Conclusions
In this study, we have proposed a new approach to Earth's gravity field modeling based on GPS-derived kinematic orbits and baselines from formation flying satellites. In this approach, kinematic orbits of the reference satellite and kinematic baselines between the reference satellite and its accompanying satellite are taken as observations to perform gravity field recovery (denoted as the 'baseline approach' in the text). Compared to the orbit approach, which is based solely on kinematic orbits, systematic errors in observations and background force models can be suppressed in the proposed approach. To demonstrate the added value of the proposed approach, two sets of monthly gravity field solutions have been produced based on one year of kinematic orbits and baseline vectors of the GRACE satellites using both the baseline and orbit approaches. A comparison in the spectral domain shows that the gravity field solutions beyond degree 10 can be notably improved when the baseline approach is applied, with cumulative errors up to degree 20 and 60 being reduced by 15% and 34%, respectively, as compared to the orbit approach. Analysis in the spatial domain shows that the GPS-based gravity field solutions obtained with the baseline approach are more consistent with the KBR-based solutions. This is evidenced by smaller RMS differences and larger correlation coefficients with respect to the KBR-based solutions. These results, therefore, demonstrate that the proposed baseline approach can provide better time-varying gravity field solutions than the orbit approach. Our findings may facilitate, among others, bridging the gap between GRACE and GFO missions by using hl-SST GPS data from formation flying satellites.