Gravity field models derived from Swarm GPS data

It is of great interest to numerous geophysical studies that the time series of global gravity field models derived from Gravity Recovery and Climate Experiment (GRACE) data remains uninterrupted after the end of this mission. With this in mind, some institutes have been spending efforts to estimate gravity field models from alternative sources of gravimetric data. This study focuses on the gravity field solutions estimated from Swarm global positioning system (GPS) data, produced by the Astronomical Institute of the University of Bern, the Astronomical Institute (ASU, Czech Academy of Sciences) and Institute of Geodesy (IfG, Graz University of Technology). The three sets of solutions are based on different approaches, namely the celestial mechanics approach, the acceleration approach and the short-arc approach, respectively. We derive the maximum spatial resolution of the time-varying gravity signal in the Swarm gravity field models to be degree 12, in comparison with the more accurate models obtained from K-band ranging data of GRACE. We demonstrate that the combination of the GPS-driven models produced with the three different approaches improves the accuracy in all analysed monthly solutions, with respect to any of them. In other words, the combined gravity field model consistently benefits from the individual strengths of each separate solution. The improved accuracy of the combined model is expected to bring benefits to the geophysical studies during the period when no dedicated gravimetric mission is operational.


Introduction
The Earth's Magnetic Field and Environment Explorers (more commonly known as the Swarm satellites), launched in November 2013, aim at improving the knowledge of the geomagnetic field (Haagmans 2004;Friis-Christensen et al. 2006, 2008Olsen et al. 2013). Nevertheless, the data gathered by the GPS and star tracker instruments are also useful for geodetic applications, in particular to measure the temporal variations of Earth's gravity field.
The main source of highly accurate gravimetric data describing the temporal changes in Earth's gravity field is the GRACE mission (Tapley et al. 2004), launched in March 2002. From these data, monthly gravity field models describing the global mass variations at scales of 300 km at Earth's surface are routinely produced (e.g. Bettadpur 2012;Meyer et al. 2012;Watkins and Yuan 2012;Dahle et al. 2012;Ditmar et al. 2013;Lemoine et al. 2013). The GRACE satellites exploit the low-low satellite-to-satellite tracking (ll-SST) measurement principle and take advantage of the KBR instrument, which provides inter-satellite range (ISR, or range for short) with μm accuracy (Dunn et al. 2003;Frommknecht et al. 2006;Kim and Lee 2009). On the other hand, the data collected by the Swarm satellites are regarded as highlow satellite-to-satellite tracking (hl-SST) observations and are precise only at the mm level. These data make it possible to estimate kinematic orbits (KOs), which consist of epoch-wise geometric fits of the pseudo-ranges derived from the hl-SST data and describe the position of the satellites with cm precision. As a consequence, the gravity field models produced from these data describe only the largest gravitational features, at around 2000 km. This was not only predicted by simulations (Gerlach and Visser 2006;Wang and Rummel 2012), but also shown from preliminary studies based on actual data Jäggi et al. 2014;Bezděk et al. 2014a; Dahle et al. 2014).
The GRACE satellites are expected to stop gathering data in the near future, due to the natural orbit decay and the degradation of the on-board batteries, in spite of the efforts that have been taken to prolong the mission lifetime (Herman 2012). As a consequence of the requirement to continue monitoring the Earth system, the GRACE Follow-On (GRACE-FO) mission is set to replace GRACE, albeit no sooner than August 2017 . To keep the time series of gravity field models uninterrupted in the gap between the two gravimetric missions, alternative data must be used. Using the hl-SST data, such as from the Swarm satellites, is a good option as demonstrated by Weigelt et al. (2013), Baur (2013), Sośnica et al. (2014) and Weigelt et al. (2014) on the basis of simulated data, as well as by Jäggi et al. (2016) and Bezděk et al. (2016) using Swarm hl-SST data.
A number of institutes have produced gravity field models from Swarm data, with different approaches (refer to Table 1).
These models have been produced from different KOs; the solutions from AIUB were produced on the basis of KOs produced at the same institute and the solutions from ASU and IfG were generated from the KOs produced at IfG. Although the models in Table 1 are accurate, there is still some room for improvement, to the benefit of the geophysical studies that exploit these data. This study proposes to determine the quality of the monthly gravity field models estimated from Swarm data from all three satellites. In analogy to the activities of the European Gravity Service for Improved Emergency Management (EGSIEM) currently performed for GRACE monthly solutions , we combine the various solutions, so that the advantages of one method should compensate for the weaknesses of another. For example, the short-arc approach is known to be particularly sensitive to temporal aliasing, i.e. the insufficient temporal sampling of the fast variations of Earth's gravity field (steps to address this issue are discussed by Kurtenbach et al. 2012); on the other hand, the acceleration approach is very sensitive to errors in the orbits (Ditmar et al. 2012). We determine the maximum spatial resolution of these models on a monthly basis and identify the geographical regions where the temporal variations of Earth's gravity field are in good agreement with what is observed by GRACE. The degree 2 coefficients are those provided in the respective models, i.e. no SLR data are considered in this study.
The Swarm Alpha and Charlie satellites orbit the Earth at around 460 km altitude, with a east-west separation no larger than 1.5°, while the Swarm Bravo satellite is at 520 km. All three satellites are in a near-polar orbit. The different inclination between the lower pair and Swarm Bravo (84.7° and 88°, respectively) will make the orbital planes perpendicular, in order to better observe Earth's magnetic field.
Although the Swarm satellites are equipped with accelerometers, the data gathered by these instruments have revealed inadequate for geodetic purposes (Siemes et al. 2015). Therefore, alternative strategies for processing the non-gravitational acceleration are employed, as listed in Table 2. In the same table, the handling of tidal and nontidal Atmosphere and Ocean processes is also shown.

Methods
This section starts with a brief overview of the gravity field estimation methods relevant to the models considered in this study. Second, the combination of the individual solutions is addressed. Finally, we introduce the concept of degree correlation, which is used to determine the wavelengths at which the Swarm solutions agree with the GRACE models.

Gravity field estimation methodology
In the following paragraphs, we give a brief description of the gravity field estimation techniques used in producing the gravity field models considered in this study. The variational equations approach (Reigber 1989) connects the measured distances from the GPS satellites, in the form of either GPS observations or GPS-derived KO data, to a set of unknown parameters which may include Stokes coefficients, initial state vectors, empirical accelerations, drag coefficients, instrument calibration parameters (e.g. accelerometer or metrology system) and other parameters which play a role in the dynamic equations of motion of the satellite(s). This is accomplished by linearising the mathematical model describing the motion of the satellite when considering a priori reference gravity field model and remaining a priori assumptions (such as analytically derived initial state vectors, initial guess for drag coefficients, calibration parameters provided by instruments manufacturer). The linearisation is performed (usually numerically) around the reference model response, i.e. the values of the unknown parameters computed from the reference model when considering the a priori assumptions. The celestial mechanics approach (CMA), used by AIUB to produce gravity field models from Swarm data, is a variant of the variational equations approach (Beutler et al. 2010).
The boundary value problem for short arcs, or shortarc approach in short, considers the double integration of Newton's equation of motion, resulting in a boundary value problem in the time domain (Mayer-Gürr et al. 2005a, b;Mayer-Gürr 2006). The relation between the unknown parameters and the observations, i.e. the satellite trajectory in case of hl-SST data and the range in the case of ll-SST data, is not linear, requiring linearisation (much like the variational equations approach). The unknown parameters are the Stokes coefficients and the boundary state vectors of each orbit arc; the observations are defined by the orbit and additionally by the range if ll-SST data are also available.
The acceleration approach exploits Newton's equation of motion directly (Rummel 1979). The method connects the double-differentiated orbital positions (or doubledifferentiated range) to the (relative) forces acting on the satellite. As a consequence, the acceleration approach (1) avoids numerical integration operations characteristic of other approaches, using instead the numerically efficient differentiation operator, and (2) uses the force model directly. The observations are usually transformed to the (quasi-) inertial reference frame before differentiation to avoid frame accelerations. The differentiation of noisy observations leads to the amplification of the high-frequency noise; however, it is possible to handle the highfrequency noise with a decorrelation procedure, such as frequency-dependent data weighting (FDDW) (Klees and Ditmar 2004;Ditmar and Sluijs 2004;Bezděk et al. 2011).
There are two subcategories of the acceleration approach: the point-wise (Austen et al. 2002;Reubelt et al. 2003;Chen et al. 2008) and averaged (Ditmar and Sluijs 2004;Liu 2008;Bezděk et al. 2014b). The pointwise acceleration approach differentiates interpolating functions fitted to the observations instead of the observations themselves, while the averaged acceleration approach differentiates the observations directly.

Combination of the gravity field models
The individual solutions, as provided by AIUB, ASU and IfG, are produced from the data of all three Swarm satellites. Along with the Stokes coefficients, the data files also contain the error estimates. These error estimates are formal in case of AIUB and of empirical nature in case of ASU and IfG. This can be seen in Fig. 1, particularly at the low degrees, where the error estimates of AIUB are overly optimistic.
The combination of the individual models, done on a monthly basis, cannot take the different types of error estimates into account because the combined solution would be unrealistically biased towards the model with more optimistic errors. For this reason, we combine the models described by the Stokes coefficients C ℓm and the corresponding error estimates σ ℓm using simple arithmetic averaging, as described by Eqs. 1 and 2.
(1) C The combination is done up to the maximum common degree and order, i.e. 40.

Degree correlation
The purpose of the degree correlation is to have a perdegree metric of the Swarm solutions that provides an indication of which spatial lengths are in good agreement with the GRACE solutions. This is accomplished by computing the correlation coefficient of the Stokes coefficients with the same degree in the GRACE and Swarm solutions (Tapley et al. 2004, supporting online material): It is important to note that the models C (Swarm) ℓm and C (GRACE) ℓm in Eq. 3 represent the time-variable signal. The time-variable signal is computed from the Swarm and GRACE solutions by subtracting a reference model, which in the case of this study was chosen to be GGM05G (Tapley et al. 2013). As in every correlation coefficient, a scaling of the gravity field model coefficients is not relevant, but the way in which they change with order is. This analysis in the frequency domain reflects the correlation of the spatial domain; it is equivalent to computing the spatial correlation between maps produced from individual degrees.

Results and discussion
The results are discussed in the following order: first, examples of degree amplitude spectra of the Swarm models are shown; second, the practical spatial resolution of the Swarm models is determined; finally, a more detailed comparison with the GRACE models is conducted, determining the regions where there is good agreement between the gravity field models produced from Swarm and GRACE data. We consider 10 months of data, from September 2014 to September 2015, excluding December 2014, May and June 2015. These excluded months could not be considered because the corresponding GRACE models are not available (cf. Information System and Data Center website 2016).

Degree amplitude spectra
In this section, we show illustrative examples of the timevariable signal in the Swarm models and their residual relative to GRACE.

Time-variable signal
The monthly Swarm solutions with respect to the static field GGM05G (which represent the time-variable signal in their lower degrees) for March 2015 are shown in Fig. 2. Their combination and the German Research Centre for Geosciences (GFZ) release 05a (produced from KBR data from GRACE and henceforth referred to as GFZ RL05a for brevity Dahle et al. 2012) is also shown in the same figure. The remaining months are somewhat similar, depicting the hydrological cycle over the year. It can be seen that all Swarm models have degrees with amplitudes comparable to GRACE at degrees below 10. Within this degree range and on average, the amplitude of the coefficients decreases slightly with increasing degree, indicating that they may represent some geophysical signal. The geophysical signal retrieved by Swarm is mainly hydrology because the tides are modelled and the non-tidal atmospheric and ocean gravity field variations are either modelled (as is the case with the ASU and IfG solutions) or shown to have no significant influence in the solutions of AIUB (cf. Jäggi et al. 2014, Figure 4, for the case of GOCE). On the other hand, above degree 10, the amplitude clearly increases in higher degrees, suggesting that noise is dominant at these spatial frequencies. A detailed analysis of the agreement between the spatial frequencies of the Swarm and GRACE models is presented in the next section.
Additionally, the degree amplitudes of the Swarm models are consistently higher than those of GRACE. The combined Swarm model tends to lower these amplitudes, which suggests that the errors in the low degrees are not negligible, as illustrated below.

Error estimates
The difference between the GRACE model and the Swarm models gives an estimation of the errors in the latter models, as shown in Fig. 3. This is a safe assumption given the fact that the gravity field models produced from GRACE take advantage of the much more accurate KBR data, in comparison with the hl-SST data used in producing the Swarm models.
Comparing Figs. 2 and 3, the degree amplitude of the error (as given by the residual relative to GRACE) is lower than that of the signal at the low degrees. These results demonstrate that the Swarm gravity field models describe, at a limited spatial resolution, the temporal variations of Earth's gravity field.
We took March 2015 as an illustrative month. Therefore, it does not mean that the relative quality of the models produced by AIUB, ASU and IfG is always the same as what is shown in Fig. 3. This fact is clearly shown in Fig. 4 top, which shows the cumulative degree amplitude up to degree 12 of the difference between Swarm and GRACE, for all institutes. It is also evident that the combination of the models of the considered institutes, even using the simple arithmetic averaging procedure, yields lower differences relative to GRACE.
As depicted in Fig. 4 top (and also in Fig. 3, with isolated exceptions at the lowest degrees), it is noteworthy that the combined model is consistently closer to GRACE than any of the individual models. The strengths of the different gravity field model estimation strategies complement each other sufficiently to improve the combined model. This result illustrates that there is no immediate strong motivation to choose one particular strategy over the other, and the best approach is to consider the combination of the widest possible range of strategies. The error estimates shown in Fig. 4 top rely on GRACE data (henceforth called estimated errors). Once this mission ends, such estimates will no longer be possible. For that reason, we show the cumulative amplitude at degree 12 of the error provided with the Swarm models in Fig. 4 bottom (referred to as predicted errors). Although the number of months is limited, it is possible to say that, in general, the month-to-month variations of the predicted and estimated errors match (i.e. if one increases, so does the other), with some exceptions (e.g. April in case of the combined solution, February in case of IfG and April in case of ASU). Regarding the scale between the two error types, the ones provided in AIUB's solutions underpredict the estimated errors roughly by a factor of 10, a factor of 4-6 in case of ASU and the combined solution, and a factor of 3 in case of IfG, except for the last 3 months, which resembles closely ASU. It is expected that once all institutes produce empirically calibrated errors, their combinations will likely produce more accurate error predictions. Nevertheless, this result gives some confidence in evaluating the quality of the (combined) Swarm model in the absence of GRACE data.

Estimation of the practical spatial resolution of the Swarm models
To have a quick identification of the maximum degree that represents a geophysical signal in the Swarm models, we use the concept of degree correlation. Assuming that the GRACE models are much more accurate than the Swarm models, the degree correlation between the Swarm and GRACE solutions gives an indication of how much signal is in the Swarm solutions and at which degrees. Figure 5 shows the result of applying Eq. 3 to the combined Swarm gravity field solution and the GRACE GFZ RL05a model, for the month of March 2015.
Another way of illustrating the level of correlation between the gravity fields from Swarm and GRACE is to calculate the cumulative degree correlation. This quantity accumulates the values of the degree correlations as the degree increases, as shown in Fig. 6.
As long as the cumulative degree correlation keeps growing with increasing degree, it is safe to say that the two gravity field models correlate well. As soon as the slope of the cumulative degree correlation becomes shallow and horizontal, the accumulation of strong correlations no longer takes place and the two gravity field models may not represent the same geophysical process. For March 2015, Fig. 6 shows that the correlation grows strongly up until degree 18. This figure also illustrates that the combined model correlates better with the GRACE model than any individual model.
In the remaining months, the (cumulative) degree correlations are similar but some differences are observed. Table 3 illustrates the differences between the degree correlations for the ten months considered in this study. The point of showing this table is to illustrate how far up the degree range the Swarm gravity fields are consistent with the GRACE models. The first column lists the months under analysis. The cumulative correlation at degree 10 is listed in the second column. The third column intends to show which of the lowest consecutive degrees are well correlated with GRACE data. Since the lowest degrees should be well correlated, it is expected that the correlation coefficient gradually decreases with increasing degree. In some months, this is not the case, because certain degree is poorly correlated, breaking the consecutive high correlation at the low degrees. The fourth column lists the degrees above the one reported in the third column that correlate well with GRACE. It might be the case that one particular degree has a poor correlation with GRACE but the following ones correlate well; therefore, this column intends to give more context to the second column (which by itself might not be sufficiently descriptive for the agreement between Swarm and GRACE at a particular month). Finally, the fifth column lists the largest degree below which the cumulative degree correlation is clearly increasing. The values in this column were determined visually, by locating the degree after which the cumulative degree curve becomes horizontal or with a negative slope. Fig. 4 Time series of the cumulative amplitude of the error at degree 12 of the gravity field models considered in the study, for the considered months, estimated from GRACE data (top) and predicted by the Swarm models (bottom) In general terms, Table 3 illustrates that the agreement of the combined Swarm solutions with GRACE varies somewhat over the considered months. The strongest correlations (reported in the second column) takes place in later months, namely August and September 2015. The reason might be twofold: a higher quality of the GPS data due to improvement in the receiver settings and the low ionospheric activity that is characteristic for summer months. October and November 2014, and January 2015 are months with the lowest value of the minimum degree with correlation lower than 0.4 (degree 2, as shown in the third column); at higher degrees, the correlation is strong (as shown in column 4), indicating that not all geophysical signal is restricted to degree 2. From the fifth column, one could say that October 2014 and August 2015 are not particularly good solutions (because the largest degree with significant degree correlation increase is 13, the lowest value); however, the remaining columns, namely the second and third columns, clearly indicate that August 2015, unlike October 2014, is a very good solution: it has the second largest cumulative correlation at degree 10 (second column), as well as the second largest value of the minimum degree with correlation lower than 0.4 (third column). In contrast, the solution for November has the largest value in the fifth column (listing the largest degree with significant degree correlation increase) but the second worst cumulative correlation at degree 10 (second column). It seems that, for this month, the signal is spread over a wide degree range (up to degree 22). From Table 3, especially from the second and last column and without considering any month in particular, we predict that there might be some geophysical signal in the combined Swarm model up to degree 20. We take degree 20 as an optimistic prediction for the maximum spatial resolution of the Swarm models, to be assessed in the next section.

Detailed comparison with the GRACE models
To better determine the maximum practical spatial resolution and the regions well observed by the Swarm gravity field models, we analyse the models in the spatial domain, as well as the two-dimensional (2D) spatial correlation between the Swarm and GRACE models [refer to Ditmar et al. (2012, Section 4.1) for the description of how the 2D spatial correlation is computed].
The analysis is done after applying a Gaussian smoothing with radii 833, 625 and 500 km to both Swarm and GRACE models. These smoothing radii are roughly related to degrees 12, 16 and 20, respectively. These radii relate to wavelengths of 1666, 1250 and 1000 km, which divide the half of Earth circumference, roughly 20,000 km, to produce the reported spherical harmonic degrees. The chosen smoothing radii are the result of the analysis presented in the previous section, which demonstrated that there might be some geophysical signal in the Swarm models up to degrees 12-20.
Referring to Fig. 7, the time-variable signal represented by the combined Swarm model is shown in the top left, and for GRACE on the top right. The bottom figures represent relations between these two models, more specifically their difference on the bottom left and their 2D spatial correlation on the bottom right. The figures in the bottom panels only show the gravity field for the land areas in order to derive meaningful statistics, since it is not expected that the data collected by the Swarm satellites are able to capture the smaller gravity signal associated with ocean dynamic topography. Figures 8 and 9 represent the same quantities but for smoothing radii of 625 and 500 km (Fig. 7 is the result of a smoothing radius of 833 km).
In these figures, the colour scale of the grids of the models and their difference (i.e. the top left, top right and bottom left plots) is the same. Figure 7 illustrates that the March 2015 Swarm model with a smoothing radius of 833 km correlates well with the GRACE model in many regions of the world, particularly Greenland, northern Russia, Argentina and Chile, and South Africa, as indicated by the dark red colour of these regions in the bottom right plot. The regions where the Swarm model does not correlate well with GRACE is Australia, Southeast Asia, West Antarctica and Central America. In spite of this, the mean correlation over land is 0.66. The root mean square (RMS) difference between the Swarm and GRACE models is on average 3.8 mm geoid height in land areas. There seems to be an overestimation of the amplitude of the signal in the Swarm model, as seen by the stronger colours in the top left figure, relative to the top right; the reason for this is still under investigation. Figure 8, showing the same as Fig. 7 but after applying Gaussian smoothing with 625 km radius, illustrates that this smoothing radius is sufficiently small to allow some smaller-scale features of the Swarm models to appear. These features are clearly associated with errors and are not useful to geophysical applications. Consequently, the mean correlation between the Swarm and GRACE models over land decreases to 0.63 (in comparison with 0.66 when applying a smoothing radius of 833 km). The RMS difference increases to 4.3 mm geoid height (from 3.8 mm for 833 km smoothing radius).
Finally, Fig. 9 illustrates the effect of applying a Gaussian smoothing with 500 km radius, which is unsuitable to suppress the noise in the Swarm models. The noise suppression is deficient enough to allow the characteristic signature of the geomagnetic equator to be seen on the top left figure; this feature is further described by Jäggi et al. (2016). As a consequence of the higher noise, the spatial correlation further decreases to 0.58 (compared to 0.66 when applying a smoothing radius of 833 km) and the RMS difference increases to 5.2 mm geoid height (from 3.8 mm for 833 km smoothing radius).
On the basis of Figs. 7, 8 and 9, it is possible to say that a smoothing radius of 833 km is the minimum to have good suppression of noise in the Swarm models and a good correlation with GRACE models.
Most solutions (with 833 km smoothing) have a correlation coefficient with GRACE above 0.6 and below 0.7, as illustrated in the top plot of Fig. 10, showing how the correlation between the Swarm and GRACE models varies over the 10 monthly solutions considered in this study. Although this is indicative of the agreement between the spatial variations of the Swarm and GRACE solutions, their difference is lost in this representation. For this reason, the bottom plot of Fig. 10   between the Swarm and GRACE models representing the time-variable signal, after applying Gaussian smoothing of the considered radii, for the 10 months considered in this study. In this representation, the considered month of March 2015 is the third worse and has almost the same value as November 2014, the month which has the worst correlation with GRACE data. This discrepancy illustrates the need to consider both the spatial correlations, as well as the actual differences between Swarm and GRACE. In addition, it also suggests that March 2015 is particularly affected by the overestimation of the geophysical signal in the Swarm models, while this issue is of less importance in the remaining solutions.
From these results, we have shown that the Swarm gravity field models represent the same geophysical signals as GRACE with RMS differences ranging from 2 to 4 mm geoid height and a spatial resolution up to degree 12 (or wavelengths of 1666 km).

Conclusions
In this study, we have compared the gravity field models produced from Swarm hl-SST data and from GRACE KBR data. We have combined the solutions computed with three different gravity field estimation methods and confirmed that the combined gravity field model consistently benefits from the individual strengths of each separate solution. We have determined that the spatial resolution with which the combined Swarm models are able to describe the time-variable gravity field of the Earth, on a monthly basis, is 1666 km (or up to degree 12). Compared to GRACE, the Swarm solutions differ on average by 2-4 mm in the geoid height and usually have a spatial distribution with a correlation coefficient better than 0.6 (with exception of 2 out of the 10 months under analysis).
The quality of the gravimetric data collected by the Swarm satellites is expected to increase as the mission progresses. First and foremost, the altitude of the satellites naturally decays with time, which will produce a larger gravitational signal and, therefore, higher-quality hl-SST gravity field estimates. The larger non-gravitational accelerations can be handled through frequencydependent data weighting (FDDW) (as demonstrated by Ditmar et al. 2007) or absorbed in the empirical accelerations (which effectively act as FDDW). Second, there are improvements being continually implemented in the GPS receiver, with the purpose of increasing the accuracy of the collected data (van den IJssel et al. 2015). Finally, as more data are collected, we are able to understand their peculiarities and better compensate for its deficiencies, namely the ionospheric disturbances over the geomagnetic equator and polar regions (Arnold et al. 2015).
On the basis of these results, we conclude that the Swarm hl-SST data are a source of additional information, together with data from other satellite missions tracked by hl-SST and SLR, to maintain the continued monitoring of the time-variable gravity field of the Earth in case there is a gap between GRACE and GRACE-FO. Not only are the Swarm data global but also gathered by three satellites, each equipped with two cold-redundant GPS receivers (i.e. these receivers cannot operate concurrently, but the redundant one can be switched on if needed). Therefore, it is very likely that these data will continue to be gathered well into the future and, in doing so, help to monitor Earth's system.