Articles | Volume 12, issue 2
Data description paper
22 Jun 2020
Data description paper |  | 22 Jun 2020

Description of the multi-approach gravity field models from Swarm GPS data

João Teixeira da Encarnação, Pieter Visser, Daniel Arnold, Aleš Bezdek, Eelco Doornbos, Matthias Ellmer, Junyi Guo, Jose van den IJssel, Elisabetta Iorfida, Adrian Jäggi, Jaroslav Klokocník, Sandro Krauss, Xinyuan Mao, Torsten Mayer-Gürr, Ulrich Meyer, Josef Sebera, C. K. Shum, Chaoyang Zhang, Yu Zhang, and Christoph Dahle

Although the knowledge of the gravity of the Earth has improved considerably with CHAMP, GRACE, and GOCE (see appendices for a list of abbreviations) satellite missions, the geophysical community has identified the need for the continued monitoring of the time-variable component with the purpose of estimating the hydrological and glaciological yearly cycles and long-term trends. Currently, the GRACE-FO satellites are the sole dedicated provider of these data, while previously the GRACE mission fulfilled that role for 15 years. There is a data gap spanning from July 2017 to May 2018 between the end of the GRACE mission and start the of GRACE-FO, while the Swarm satellites have collected gravimetric data with their GPS receivers since December 2013.

We present high-quality gravity field models (GFMs) from Swarm data that constitute an alternative and independent source of gravimetric data, which could help alleviate the consequences of the 10-month gap between GRACE and GRACE-FO, as well as the short gaps in the existing GRACE and GRACE-FO monthly time series.

The geodetic community has realized that the combination of different gravity field solutions is superior to any individual model and set up the Combination Service of Time-variable Gravity Fields (COST-G) under the umbrella of the International Gravity Field Service (IGFS), part of the International Association of Geodesy (IAG). We exploit this fact and deliver the highest-quality monthly GFMs, resulting from the combination of four different gravity field estimation approaches. All solutions are unconstrained and estimated independently from month to month.

We tested the added value of including kinematic baselines (KBs) in our estimation of GFMs and conclude that there is no significant improvement. The non-gravitational accelerations measured by the accelerometer on board Swarm C were also included in our processing to determine if this would improve the quality of the GFMs, but we observed that is only the case when the amplitude of the non-gravitational accelerations is higher than during the current quiet period in solar activity.

Using GRACE data for comparison, we demonstrate that the geophysical signal in the Swarm GFMs is largely restricted to spherical harmonic degrees below 12. A 750 km smoothing radius is suitable to retrieve the temporal variations in Earth's gravity field over land areas since mid-2015 with roughly 4 cm equivalent water height (EWH) agreement with respect to GRACE. Over ocean areas, we illustrate that a more intense smoothing with 3000 km radius is necessary to resolve large-scale gravity variations, which agree with GRACE roughly at the level of 1 cm EWH, while at these spatial scales the GRACE observes variations with amplitudes between 0.3 and 1 cm EWH. The agreement with GRACE and GRACE-FO over nine selected large basins under analysis is 0.91 cm, 0.76 cm yr−1, and 0.79 in terms of temporal mean, trend, and correlation coefficient, respectively.

The Swarm monthly models are distributed on a quarterly basis at ESA's Earth Swarm Data Access (at, last access: 5 June 2020, follow Level2longterm and then EGF) and at the International Centre for Global Earth Models (, last access: 5 June 2020), as well as identified with the DOI (Encarnacao et al., 2019).

1 Introduction

Swarm is the fifth Earth Explorer mission by the European Space Agency (ESA), launched on 22 November 2013 (Haagmans, 2004; Friis-Christensen et al., 2008). Its primary objective is to provide the best ever survey of the Earth's magnetic field and its temporal variations as well as the electric field of the atmosphere (Olsen et al., 2013). Swarm consists of three identical satellites, two flying in a pendulum formation (side by side, converging near the poles) at an initial altitude of about 470 km and one at an altitude of about 520 km, all in near-polar orbit. In addition to a sophisticated instrument suite for observing the geomagnetic and electric fields, the Swarm satellites are equipped with high-precision, dual-frequency Global Positioning System (GPS) receivers, star trackers, and accelerometers. Many recent studies and activities have shown the feasibility of observing the Earth's gravity field and its long-wavelength temporal variations with high-quality GPS receivers on board low-Earth orbit (LEO) satellites (Zehentner and Mayer-Gürr, 2014; Bezdek et al., 2016; Dahle et al., 2017). For Swarm, Teixeira da Encarnação et al. (2016) successfully demonstrated the observation of long-wavelength temporal gravity. They produced solutions using three different approaches and showed that their combination resulted in improved observability of time variable gravity, a principle that has been suggested in the frame of the initiative of the European Gravity Service for Improved Emergency Management (EGSIEM) (Jäggi et al., 2019) and demonstrated for gravity field solutions based on the Gravity Recovery and Climate Experiment (GRACE) (Jean et al., 2018).

An important driver for producing LEO GPS-based gravity field solutions is to guarantee long-term observation of mass transport in the Earth system. The geophysical community has identified the need for continued monitoring of time-variable gravity for estimating the hydrological and glaciological yearly cycles and long-term trends (Abdalati et al., 2018). The US–German GRACE mission (Tapley et al., 2004) was by far the most important space-borne global provider of the needed data for the period from April 2002 until July 2017. GRACE Follow-On (GRACE-FO) was launched in May 2018 and is expected to continue the high-quality observation of Earth's time-variable gravity field for at least 5 years (Flechtner et al., 2016). Thus a time gap exists between the GRACE and GRACE-FO missions, and, importantly, no missions have yet been selected for the post-GRACE-FO period. It can thus be claimed that the only guarantee for sustained observation of time-variable gravity from space is constituted by space-borne GPS receivers on LEO satellites. Moreover, the associated data can be used to fill the gap between the GRACE and GRACE-FO missions (be it with a different quality in terms of spatial and temporal resolution). The measurement of Earth's gravitational changes with Swarm is further motivated by (i) the need to increase the accuracy of global mass estimates in order to properly quantify global sea-level rise and (ii) the opportunity to provide independent estimates of temporal variations in low-degree coefficients, in particular related to C2,0 and C3,0, which are weakly observed by GRACE. Under this motivation, this paper aggregates a series of studies and analyses that, respectively, motivate our processing choices and demonstrate the capabilities of the combined Swarm models to observe mass transport processes at the surface of the Earth on a monthly basis, in a way that is superior to any of its individual models.

The studies aim at improving Swarm-based observation of long-wavelength time-variable gravity in support of the operational delivery of Swarm-based gravity field solutions. It is a continuation of the activities described in Teixeira da Encarnação et al. (2016), which included the production of gravity field solutions using three different methods, referred to as the celestial mechanics approach (CMA) (Beutler et al., 2010), decorrelated acceleration approach (DAA) (Bezdek et al., 2014), and short-arc approach (SAA) (Mayer-Gürr, 2006). In this work, a fourth method, referred to as the improved energy balance approach (IEBA) (Shang et al., 2015), is added. The combination of the four gravity field solutions into combined models will be more advanced than in Teixeira da Encarnação et al. (2016), where a straightforward averaging was applied. In the results presented in this work, the weights are derived from variance component estimation (VCE) analogously with Jean et al. (2018), in order to arrive as close as possible at statistically optimal combined solutions, given the available combination strategies, as described in Sect. 2.5.

The nominal gravity field solutions will be based on kinematic orbit (KO) solutions, which consist of time series of position coordinates. These time series can be considered a condensed form of the original GPS high–low satellite-to-satellite tracking (hl-SST) observations, with no effect from dynamic models for the LEO satellites (the positions of the GPS satellite themselves are based on dynamic models, as usual). Three different KO solutions are produced by the Delft University of Technology (TUD), Astronomical Institute of the University of Bern (AIUB), and the Institute of Geodesy Graz (IfG) of the Graz University of Technology (TUG) (van den IJssel et al., 2015; Jäggi et al., 2016; Zehentner, 2016).

We also tested another potential innovation that could conceptually lead to improved gravity field solutions, that is the use of kinematically derived baselines for the two Swarm satellites flying in a pendulum formation. Kinematic baselines (KBs) between two LEO formation-flying spacecraft can typically be derived with much better precision than the absolute positions by making use of ambiguity-fixing schemes and due to cancellation of common errors (Kroes, 2006; Allende-Alba et al., 2017). The possible added value of KBs for the observation of temporal gravity field variations will be assessed making use of two different KB solutions by TUD (Mao et al., 2017) and the AIUB (Jäggi et al., 2007, 2009).

We also present a comparison of the quality of gravity field retrievals from Swarm C observational data making use of either the available accelerometer product for this satellite (Doornbos et al., 2015) or two different non-gravitational acceleration force models.

This paper is organized as follows. More details about the methodology are provided in Sect. 2. Results are included and discussed in Sect. 3. A summary, conclusions, and outlook are given in Sect. 4.

For the sake of brevity, we will refer to GRACE and GRACE-FO data simply as GRACE data, unless there is the need to be more specific. We also interchangeably use the terms solution (when relevant to a set of Stokes coefficients) and gravity field model (GFM).

The operational activities currently under way pertaining to the combined models described in this article are conducted in the frame of the Combination Service of Time-variable Gravity Fields (COST-G), under the umbrella of International Gravity Field Service (IGFS) International Association of Geodesy (IAG) (Jäggi et al., 2020), with additional support from the Swarm Data, Innovation and Science Cluster (DISC) and funded by ESA. The Swarm monthly models are distributed on a quarterly basis at ESA's Earth Swarm Data Access (at, last access: 5 June 2020, follow Level2longterm and then EGF) and at the International Centre for Global Earth Models (, last access: 5 June 2020), as well as identified with the DOI (Encarnacao et al., 2019).

2 Methodology

In this work, we mainly intend to present the capabilities of the Swarm GFMs, in terms of their particularities and data quality, and we typically refer to the relevant methodology in supporting literature. Nevertheless, this section discusses briefly some aspects of the various stages in the processing of the models, their combination, and, to better prepare the discussion of results in Sect. 3, the approach used in the analysis of the Swarm GFMs.

2.1 Kinematic orbits

The KOs are the observations from which the GFMs are estimated, since they are solely derived from the geometric distance relative to the GPS satellites. The different KO solutions are conceptually estimated in similar ways, but with the processing strategies described in detail in the references of Table 1. Furthermore, each analysis centre (AC) makes their own choices regarding the numerous assumptions and processing options for deriving their individual KO solutions, as listed in Appendix A. The reason for the different KO solutions is to provide various options for the ACs' individual GFM processing (see Sect. 2.4) and, in this way, reduce the impact of possible KO-driven systematic errors in the combined GFMs. It also enables the ACs to select which KO solution is more advantageous to the quality of their GFMs; consider that our gravity estimation approaches may be differently sensitive to the error spectra of the various KO solutions or have different requirements on the quality of the variance–covariance information provided with the kinematic positions. This selection is done at each AC and outside the scope of the current study.

2.2 Kinematic baselines

We investigate the added value of KBs in the quality of the Swarm GFMs, as presented in Sect. 3.1. The KB solutions, much in the same way as the KOs, are conceptually computed similarly, where fixing ambiguities is a necessary processing step to achieve the highest possible precision of the derived baselines. This constitutes the main motivation to include KBs in the estimation of the Swarm GFMs. The interested reader can find details in the references of Table 2; the main processing assumptions are listed in Appendix B and brief descriptions follow.

Table 1Overview of the kinematic orbits and the software packages used to estimate them.

1 (last access: 5 June 2020). 2 (last access: 5 June 2020). 3 (last access: 5 June 2020).

Download Print Version | Download XLSX

Table 2Overview of the kinematic baselines and the software packages used to estimate them.

Download Print Version | Download XLSX

2.2.1 KBs produced at AIUB

Kinematic and reduced-dynamic baselines are determined according to the procedures described by Jäggi et al. (2007, 2009, 2012). The positions of one satellite (Swarm A) are kept fixed to a reduced-dynamic solution generated from zero-differenced (ZD) ionosphere-free GPS carrier phase observations. Reduced-dynamic orbit parameters of the other satellite (Swarm C) are estimated by processing double-differenced (DD) ionosphere-free GPS carrier phase observations with DD ambiguities resolved to their integer values. First, the Melbourne–Wübbena linear combination is analysed to resolve the wide-lane ambiguities, which are subsequently introduced as known to resolve the narrow-lane ambiguities together with the reduced-dynamic baseline determination. For the KB estimation, the same procedure may be used but it turned out to be more robust to introduce the resolved ambiguities from the available reduced-dynamic baselines and not to make an attempt to independently fix carrier phase ambiguities in the KB processing. Exactly the same carrier phase ambiguities are therefore fixed in both the reduced-dynamic and the kinematic baseline determinations.

2.2.2 KBs produced at TUD

We take advantage of a forward and backward extended Kalman filter (EKF) that is run iteratively. The EKF initially runs from the first epoch to the last epoch of each 24 h orbit arc with 5 s step. The estimated float ambiguities and the corresponding covariance matrices (which are recorded for each epoch) are used by the least-squares ambiguity decorrelation adjustment (LAMBDA) algorithm in order to fix the maximum number of integer ambiguities (subset approach). The EKF smooths both solutions according to the bi-directional covariance matrices recorded at each epoch. In the next iteration, the smoothed orbit and fixed ambiguities are set as input and it is attempted to fix more ambiguities. The procedure is repeated until no new integer ambiguities are fixed.

After the convergence of the reduced-dynamic baseline, a KB solution is produced using the least-squares (LS) method. To this purpose, the same GPS observations and fixed integer ambiguities on the two frequencies are used, where one satellite (Swarm A) is kept fixed at the reduced-dynamic baseline solution. At least five observations are required on each frequency to form a good geometry. To minimize the influence of wrongly fixed ambiguities and residual outliers, a threshold of 2-sigma of the carrier phase residual standard deviation (SD) is set, which results in eliminating around 5 % of data, on average. A further screening of 3 cm is set to the rms of the kinematic baseline carrier phase observation residual. This makes it possible to screen out the epochs that are influenced by wrongly fixed ambiguities and bad phase observations. The kinematic baseline determination is also run bi-directionally to compute two solutions that are averaged according to the epoch-wise covariance matrices.

2.2.3 Inclusion of KBs in the estimation of Swarm GFMs

We exploit the variational equations approach (VEA) (Montenbruck and Gill, 2000) implemented at IfG in the inversion of gravity field considering both KOs and KBs. The VEA and its application to KOs and KBs corresponds to the processing scheme used for the production of the ITSG-Grace2016 (Klinger et al., 2016).

Table 3Periods considered in the analysis of the added value of different types of non-gravitational accelerations.

Download Print Version | Download XLSX

We selected a number of suitable test months with varying data quality, meeting the following criteria: GRACE monthly solutions are available for validation purposes, months with good GPS data quality are included as well as months with bad data quality, and some months should overlap with the test months selected in the non-gravitational acceleration study (Sect. 2.3) for the accelerometer data tests.

The descriptions good and bad data quality refer to several issues in the context of Swarm GPS data. Good means that an error found in the Receiver Independent Exchange Format (RINEX) converter is solved (fixed since 12 April 2016), the settings of the receiver tracking loop bandwidths are optimized (several changes during lifetime), and the ionospheric activity is at a low level. In contrast, the bad data hold for time periods for which these issues are not solved and the ionospheric activity is high. Finally, the intermediate data are during periods of lower ionospheric activity (relative to early 2015) but before the GPS receiver updates. In total we have selected seven test months: January and March 2015 refer to bad data quality, February and March 2016 refer to intermediate data quality, and June–August 2016 refer to good data quality.

The existing software exploiting VEA at IfG handles the Swarm KB data under the same processing scheme and handling of stochastic properties of the observations adopted for the generation of the ITSG-Grace releases (Mayer-Gürr et al., 2016). The observations derived from the Swarm KBs are introduced into the gravity inversion process as if they were collected by the K-Band ranging instrument. Our software is not prepared to handle the full three-dimensional (3D) information of the KBs, and the development of this capability is outside the scope of this study.

The KBs and KO solution are selected consistently from the same AC (i.e. TUD or AIUB) when producing the gravity field solution. In total four different GFM variants have been computed: (1) hl-SST solution from TUD KOs, (2) hl-SST + low–low satellite-to-satellite tracking (ll-SST) solution from TUD KOs and KBs, (3) hl-SST solution from AIUB KOs, and (4) hl-SST + ll-SST solution from AIUB KOs and KBs. The four solution variants were produced for all seven test months.

2.3 Non-gravitational accelerations

We assessed the quality of the Swarm GFMs when the non-gravitational accelerations are modelled following two distinct approaches and when they are represented by the Level 1B (L1B) accelerometer data from Swarm C (Siemes et al., 2016). One non-gravitational acceleration model was produced at the Astronomical Institute Ondřejov (ASU) and the other at the Delft University of Technology (TUD). We selected a number of periods for our tests (see Table 3), taking care to cover as much as possible different accelerometer data variability (arising from instrument artefacts) and signal amplitude, as well as ionosphere and geomagnetic activity, to cover different regimes of non-gravitational accelerations acting on the Swarm satellites. Moreover, we also chose months when GRACE gravity field solutions are available, to facilitate validation of the Swarm GFMs.

For the ASU model, we used the in-house orbital propagator NUMINTSAT (Bezdek et al., 2009) for processing the satellite orbital data, computing the coordinate transformations and generating the modelled non-gravitational accelerations of each Swarm satellite. The computation of the non-gravitational acceleration forces requires the knowledge of the physical properties of the satellite based on the information provided by the ESA: its mass, cross section in a specific direction, radiation properties of the satellite's surface, and a macromodel characterizing approximately the shape of the Swarm satellites. For neutral atmospheric density, we made use of the US Naval Research Laboratory Mass Spectrometer and Incoherent Scatter Radar (NRLMSISE) atmospheric model (Picone et al., 2002). We estimated the drag coefficient of each satellite by means of the long-term change in the orbital elements in order to consider realistic values. Further details of our approach can be found in Bezdek (2010) and Bezdek et al. (2014, 2016, 2017).

Table 4Overview of the gravity field estimation approaches.

Download Print Version | Download XLSX

For the TUD model, the Near Real-Time Density Model (NRTDM) software was employed (Doornbos et al., 2014). This software, as part of the “official” Swarm Level 2 Processing System (L2PS) infrastructure, is used in the L1B to Level 2 (L2) processing at TUD. A variety of models and parameters related to the non-gravitational forces are available in this software. For the current study, the following selection was made.

  • The Swarm panel model (macromodel) is based on Siemes (2019).

  • The panel orientation is dictated by Swarm quaternion data.

  • The satellite aerodynamics of single-sided flat panels are computed following Sentman's equations (Sentman, 1961), assuming diffuse reflection and energy flux accommodation set at 0.93.

  • The neutral densities are derived from the NRLMSISE thermosphere model, as well as temperature and composition dependence of Sentman's equations.

  • The velocity of the atmosphere with respect to the spacecraft is based on the orbit and attitude data, atmospheric co-rotation, and modelled thermospheric wind using the Horizontal Wind Model 07 (HWM07) (Drob et al., 2008) and the Disturbance Wind Model 07 (DWM07) (Emmert et al., 2008).

  • The solar radiation pressure (SRP) is computed taking into account absorption, diffuse reflection, and specular reflection, according to optical properties of the surface materials supplied by ESA and Astrium, and it considers the varying Sun-satellite distance.

  • The Sun–Earth eclipse model takes into account atmospheric absorption and refraction, according to the analysis of non-gravitational accelerations due to radiation pressure and aerodynamics (ANGARA) implementation (Fritsche et al., 1998).

  • The Earth infrared radiation pressure (EIRP) and Earth albedo radiation pressure (EARP) are based on the ANGARA implementation.

  • The monthly average albedo coefficients and infrared radiation (IR) irradiances are derived from the Earth Radiation Budget Experiment (ERBE) data (Barkstrom and Smith, 1986).

The equations for the algorithms and references for these models are available in Doornbos (2012) with updates specific to Swarm provided by Siemes et al. (2016).

For the Swarm C accelerometer data, we took advantage of the corrected L1B along-track accelerometer data (Siemes et al., 2016), which are distributed by ESA and processed in a single batch from July 2014 to April 2016. We applied a dedicated calibration method to the Level 1A (L1A) product ACCxSCI_1A for the cross-track and radial components (Bezdek et al., 2017, 2018b), but, as shown by Bezdek et al. (2018a), this approach was unable to recover the expected signal. For this reason, the non-gravitational acceleration measurements are restricted to the available along-track Swarm C data.

2.4 Gravity field model estimation approaches

The estimation of the hl-SST GFMs takes the KOs as observations, which describe the satellite's centre of mass (CoM) motion since in their production the processing of the L1B GPS measurements is corrected for location of the GPS antenna phase centre with the L1B Swarm attitude data. The KOs are suitable to gravimetric studies due to their purely geometric nature. Through a parameter estimation procedure, i.e. one of the strategies listed in Table 4, the gravity field parameters are derived from a functional relationship between the kinematic positions and gravity field parameters. Complementary to the KOs, numerous processing choices are made by the four gravity field ACs, as enumerated in Appendix C.

Each AC selects one KO solution to produce their so-called individual GFMs, as listed in Table 4. In contrast, the combined GFMs are derived from these individual solutions, as discussed in Sect. 2.5. The following subsections provide a brief recap of the selected methods. Elaborate details can be found in the cited literature.

2.4.1 Celestial mechanics approach

The celestial mechanics approach (Beutler et al., 2010), used at AIUB, is a variation in the traditional variational equations approach (Reigber, 1989), which linearizes the relation between the kinematic positions and the unknown Stokes coefficients as well as other unknown parameters that play a role in the dynamic model described by the equations of motion, such as initial state vectors, empirical accelerations, drag coefficients, and instrument calibration parameters (possibly) amongst others. Pseudostochastic pulses or accelerations are estimated to mitigate deficiencies of the a priori force model. The CMA has successfully been applied for gravity field determination from a number of LEO satellites, e.g. Meyer et al. (2019b).

2.4.2 Decorrelated acceleration approach

The decorrelated acceleration approach (DAA) (Bezdek et al., 2014, 2016) used at ASU connects the double-differentiated kinematic positions to the external forces acting on the satellite. This approach computes the geopotential harmonic coefficients from a linear (not linearized) system of equations. The observations are first transformed to the inertial reference frame before differentiation to avoid the computation of fictitious accelerations. The differentiation of noisy observations leads to the amplification of the high-frequency noise. However, it is possible to mitigate the high-frequency noise with a decorrelation procedure. We apply a second decorrelation based on a fitted autoregressive process to take into account the error correlations of the KOs.

2.4.3 Improved energy balance approach

The traditional energy balance approach (EBA) exploits the energy conservation principle to build a relation between the residual geopotential coefficients (relative to the reference background force model) and the deviations of the KO from the reference orbit on (Jekeli, 1999; Visser et al., 2003; Guo et al., 2015; Zeng et al., 2015). The main development of the IEBA, used at The Ohio State University (OSU), concerns the handling of the noise in the kinematic position and the weighting of the potential observations. Unlike the application of this approach to GRACE ll-SST data by Shang et al. (2015), the term related to the Earth's rotation cannot be neglected in the processing of hl-SST data. From the kinematic positions, the velocity is derived with 61 data points, a sliding window, and quadratic polynomial filter similar to Bezdek et al. (2014). The polynomial coefficients of the filter are estimated in a LS adjustment, with the observation vector being composed of position residuals between the kinematic positions and the corresponding reduced-dynamic positions (integrated on the basis of the reference background force model), and the observation covariance matrix constructed from the epoch-wise variance–covariance information distributed in the KO data files. As a consequence of this orbit smoothing procedure, we discard the warm-up–cool-down edges of the daily data arcs. We further remove one cycle per revolution (CPR) sinusoidal and 3-hourly quadratic polynomial signals from the potential observations derived from the smooth kinematic positions. We also take advantage of the observation covariance matrix to weight the filtered kinematic observations in the geopotential coefficient LS inversion. We do not apply any a priori constraints nor iterate the LS estimation since we take advantage of the linear relation between the potential observations and the geopotential coefficients.

2.4.4 Short-arc approach

The short-arc approach (Mayer-Gürr, 2006), used at IfG, formulates the relation between the geopotential coefficients and the kinematic positions as the boundary value problem resulting from the double integration of the equations of motion. This approach naturally defines the initial state vector as the boundary conditions of the integral equation, which are regarded as unknowns in the LS estimation along with the Stokes coefficients and other unknown parameters, such as empirical parameters. Additionally, the kinematic positions are treated with no explicit differentiation, thus circumventing the need to suppress the amplification of high-frequency noise.

2.5 Combination

The individual GFMs are combined in the frame of the Combination Service of Time-variable Gravity Fields (COST-G) of the IGFS, applying the methods developed during the EGSIEM project (Jäggi et al., 2019). We derive VCE weights in order to produce the combined GFMs from the individual GFMs produced at AIUB, ASU, IfG, and OSU. The VCE weights are derived at the solution level according to Jean et al. (2018), considering the individual models up to degree 20 only; if this is not done, the extremely high noise at the degrees close to 40 (the maximum degree of the individual solutions) dominates the estimation of the weights, which leads to a slightly worse agreement with GRACE (Teixeira da Encarnação and Visser, 2019). Irrespective of this, the maximum degree of the combined models is the same as that of the individual models (degree 40). We also tested the combination at the level of normal equations (NEQs) (Meyer et al., 2019a) but determined that the signal content was not in as good agreement with GRACE as the combination at the level of solutions with weights derived from VCE (Teixeira da Encarnação and Visser, 2019; Meyer, 2020). We attributed this result to the difficulty in calibrating the formal error types resulting from the different gravity field estimation techniques. There is the issue of the different types of error: some provide calibrated errors (e.g. DAA), while others provide the formal errors from the LS estimates (e.g. CMA). Another issue is the different error amplitude dependence with degree, thus preventing the errors from being calibrated with a simple bias. Finally, the time-dependent levels of errors in the individual models, which change their fidelity with time, and consequentially their optimum relative weights, were also a factor preventing us from successfully performing a combination at the NEQ level.

2.6 Assumptions in the gravity field model analyses

This section describes the set of assumptions considered in the analysis done in Sect. 3.3 and 3.4. Section 3.1 and 3.2 report parallel studies that were conducted with different background force models, better suited to their respective purposes.

We have chosen the release 6 (RL06) GRACE and GRACE-FO GFMs produced at the Center for Space Research (CSR) as comparison in our analysis of the Swarm GFMs. At the spatial scales relevant to Swarm, we have no reason to expect our results would change significantly if GRACE data produced at any other AC were used instead.

Unless otherwise noted, we apply a 750 km radius Gaussian smoothing, which we motivate in Sect. 3.4.1, to isolate the signal content in the Swarm models. The geo-centre motion has been ignored in our analysis, i.e. the degree 1 coefficients are always zero. The Combined GRACE Gravity Model 05 (GGM05C) static GFM (Ries et al., 2016) is subtracted from all Swarm and GRACE solutions in order to isolate the time-variable component of Earth's gravity field. The gravity field is presented in terms of equivalent water height (EWH), except for the statistics related to the correlation coefficient or when presenting coefficient-wise time series.

We consider the entirety of the Swarm GFM time series, irrespective of the epoch-wise quality because our objective is to give a complete overview of the quality and characteristics of our models. The analysis spans all available months during the Swarm mission, i.e. between December 2013 and September 2019. When comparing Swarm and GRACE directly, the Swarm time series is linearly interpolated to the time domain defined by the epoch of the GRACE solutions, except for the GRACE-to-GRACE-FO gap, where no interpolation is performed. We detrend the time series of models at the level of the Stokes coefficients when computing non-linear statistics, notably the epoch-wise spatial rms in Figs. 6, 9, 10, 11, 13, and 16.

2.6.1 Earth's oblateness

In our analysis, the proper handling of Earth's oblateness is not a trivial problem. In the case of GRACE, the mass estimates are improved if C2,0 is augmented with satellite laser ranging (SLR) data, which are provided in the form of the time series produced by Cheng and Ries (2018). Therefore, any comparison with mass variations derived from Swarm must also have the C2,0 coefficient replaced by the same time series. One could argue that simply discarding this coefficient would suffice for any comparison but we also intend to represent the actual mass changes observed by Swarm, notably in Sect. 3.5.3, where we show mass variations over large storage basins. Unfortunately, Earth's oblateness estimates provided by Cheng and Ries (2018) are exclusively available at those epochs when there are GRACE solutions. That essentially means that interpolating these GRACE–SLR C2,0 estimates over large gaps would lead to unrealistic mass variations.

For this reason, we selected the C2,0 weekly time series from Loomis et al. (2019), since the necessary interpolation introduces negligible deviations. We are not advocating that the considered C2,0 time series is in any way superior to other solutions, e.g. those of Cheng et al. (2011) (which is only available at the middle of calendar months) or Cheng and Ries (2018) (which is only available for epochs compatible with the GRACE monthly solutions); we have selected it purely under the consideration it was the most technically convenient option for our needs.

2.6.2 Deep-ocean areas

We consider the ocean mask of the areas away from continental masses illustrated in Fig. 1. To produce this mask, we start with a grid with a unit value over land areas, convert it to the spherical harmonic (SH) domain, apply Gaussian smoothing with a radius of 1000 km, convert it back to the spatial domain, and define those grid points with values below the cut-off value of 0.9 to be in deep-ocean areas. The cut-off value was selected on the basis of trial and error with the objective of generating an ocean mask with the desired and arbitrary buffer lengths, which for the results reported here remained equal to 1000 km. This procedure pushes the boundary of an ocean mask away from continental coastal areas and ignores islands. For the spatial scales relevant to the Swarm GFMs, we propose that this procedure is adequate.

2.6.3 GRACE climatological model

In the analyses conducted in Sect. 3.5.1, where we present the time series of selected Stokes coefficients, we use a parametric representation of Earth's temporal mass changes as observed by GRACE, which we refer to as the climatological model since it captures mass variations that are present in all 15 years of GRACE data. We do not use any GRACE-FO data in this regression, in order to be able to verify the continuity of the GRACE-FO data, relative to GRACE and to substantiate any deviation that is also observed by Swarm. This parametric regression is performed on the original CSR RL06 models, i.e. before any smoothing or masking.

Figure 1Deep-ocean mask shown as dark areas.

Figure 2Agreement between the GRACE climatological model and the GRACE data, exemplified by the C3,0 coefficient.


We selected the first-order polynomial to represent bias and trend in the GRACE data. For the periodic parameters, we choose the year and semi-year periods since these are dominant signals in the GRACE and Swarm data. We also modelled the S2, K2, and K1 tidal periods, with durations of 0.44, 3.83, and 7.67 years, respectively. These periods are driven by the orbital inclination of the GRACE satellites and produce strong aliasing in the GFM time series (Ray and Luthcke, 2006; Cheng and Ries, 2017). The linear regression of the 12 parameters is done independently for each SH coefficient, up to degree 40 (in agreement with the maximum degree of the Swarm models). This results in 12 sets of Stokes coefficients, one for each of the model parameters: bias, trend, and five periods represented by their sine and co-sine components. Each set of parametric Stokes coefficients has an implicit time dependence which is evaluated coefficient-wise at the epochs of the Swarm GFMs. We illustrate the general agreement between the climatological model and the GRACE data for the case of C3,0 in Fig. 2.

We regard this model as a good representation of the Earth system; it is by definition inferior to the original GRACE time series because it truncates the signal bandwidth to discrete frequencies. In spite of this, the assumed climatological model provides a measure to which both GRACE and Swarm can be compared. The differences between GRACE and this model should be regarded as the signal augmentation that GRACE brings, not as an error. We also regard the vastly different spatial sensitivity of Swarm compared to GRACE as an additional argument that the climatological model is able to represent the Earth system in a much more accurate way than Swarm, with the exception of large atypical mass variations (which are uniquely revealed by Swarm).

3 Results

Our results are shown in the following section, where we analyse the added value of KBs in Sect. 3.1, look into the effect of including accelerometer measurements of Swarm C in Sect. 3.2, provide an overview of the quality of our individual solutions in Sect. 3.3, quantify the quality of the combined solutions in Sect. 3.4, and illustrate their signal content in Sect. 3.5.

3.1 Kinematic baselines

This section is dedicated to quantifying the benefit of exploiting KBs in the quality of the GFMs derived from Swarm data, following the motivation and procedures described in Sect. 2.2.

Due to the decreasing ionospheric activity and the changes made to the Swarm on-board GPS receivers between 2015 and 2016 (van den IJssel et al., 2016), the consistency of the KB solutions has improved. Especially in summer 2016, the overall daily SD of the difference between the reduced dynamic ambiguity-fixed and kinematic ambiguity-fixed baselines may be as low as 10–15, 4–6, and 3–5 mm on average for the radial, along-track, and cross-track directions, respectively, while it is as high as 1–3 cm for 2015 in all three directions. It should be noted, however, that daily SD is always dominated by the low-quality kinematic positions over the polar regions. Eliminating such problematic data, the difference SD is consistently under 5 mm; therefore, the internal precision of the Swarm GPS data is of very good quality.

Figure 3Difference SH degree amplitudes of all four test solutions with respect to GOCO05S, for March 2015 (a), February 2016 (b), and June 2016 (c), regarded as representative of bad, intermediate, and good data quality, respectively.


Table 5Rms of geoid height differences in millimetres for different hl-SST-only results and the ll-SST and hl-SST Swarm solutions with respect to the corresponding ITSG-Grace2016 monthly solution.

Download Print Version | Download XLSX

Figure 3 shows the degree amplitudes with respect to the static part of the GOCO release 05 satellite-only gravity field model (GOCO05S) (Mayer-Gürr, 2015) in terms of geoid heights, representative of the results for bad, intermediate, and good data quality. For comparison the corresponding month from the ITSG GRACE-only model 2016 (ITSG-Grace2016) time series is also shown. For all months it can be seen that the solutions do not differ significantly. There are small differences between the two ACs (AIUB and TUD) as well as between the hl-SST-only and the ll-SST + hl-SST solutions. Differences are larger for those months with bad data quality (2015) and at the SH degree regions dominated by noise (above degree 15), with the ll-SST and hl-SST solutions showing larger degree amplitudes. For months with good data quality (June 2016) all four solutions display much smaller differences.

To quantify the impact on the long-wavelength part of the solutions, we have compared the individual solutions to ITSG-Grace2016 monthly solutions in the spatial domain. The solutions are evaluated on an equiangular grid (1×1), reduced by the corresponding ITSG-Grace2016 monthly solution, and filtered with a 500 km Gaussian filter, and finally the rms over all grid cells is computed. The filter width was selected so as to avoid suppressing all of the signal at degrees above 20 in order to assess the impact of KBs on the high-frequency noise as well. These results are summarized in Table 5.

Table 5 confirms what is depicted in Fig. 3; i.e. the inclusion of KBs in the gravity field estimation has no significant impact on the quality of the resulting GFM. KO-only solutions are already of very similar quality when compared to KB-augmented solutions, with small differences visible in the degree amplitude plots or the spatial rms having no discernible correlation with the data period (and therefore quality). In general, this confirms the findings of Jäggi et al. (2009), in that there are some small benefits for higher degrees when using KB; this was attributed to the elimination of errors common to both satellites by using DD observations. Our results suggest that common errors are already mostly absent in the computation of the Swarm KOs. Thus we found no added value in including KBs for the quality of Swarm GFMs.

Our results contrast with those of Guo and Zhao (2019), who demonstrated a noticeable improvement when KBs are used in conjunction with KOs to derive GFMs from hl-SST GRACE data. As the authors mention, their approach benefits from the 3D KB information, thus essentially increasing by a factor of 3 the number of observations. Although these components are most likely not completely independent, they provide observations with crucial information that is not available along the line-of-sight (LoS) component, in particular along the radial direction. We also note that the geometry of the GRACE formation provides a much more stable amplitude and attitude of resulting KBs, which may benefit the ambiguity fixing and, consequently, their overall quality. In the case of Swarm, the KBs are close to zero and flip their orientation by 180 at the poles. Additionally, GRACE accelerometer data were used to represent the non-gravitational accelerations, which is less straightforward for the Swarm satellites. These differences, i.e. 3D baselines, stable baseline length, and inclusion of accelerometer data, suggest that they may be necessary conditions for a positive added value of KBs to the quality of hl-SST-only GFMs. Finally, we also point out that the improvements reported in Guo and Zhao (2019) are only above SH degree 10, where the errors start to become dominant, thus reducing the practical added value of including baselines in the estimation of hl-SST-only GFMs.

3.2 Non-gravitational accelerations

In this section, we present the inter-comparison of the three types of non-gravitational accelerations described in Sect. 2.3. Figure 4 compares three single-satellite gravity field solutions derived from Swarm C data, considering the three non-gravitational accelerations, for January 2015.

Figure 4Swarm C gravity field solutions using TUD and ASU modelled non-gravitational accelerations, as well as measured nongravitational accelerations (January 2015).


The SH degree difference amplitudes illustrate that the measured non-gravitational accelerations improve the agreement of the lowest degrees of the Swarm C monthly solution with respect to the GOCO05S model (Mayer-Gürr, 2015), which includes a time-variable component. We tested this comparison relative to the ITSG-Grace2016 monthly GFM (Mayer-Gürr et al., 2016) and observed similar results (not shown). The improvement at the lowest degrees in the Swarm C model when using observed non-gravitational acceleration data is in accordance with what was reported by Klinger and Mayer-Gürr (2016), relative to GRACE gravity field recovery.

In view of the lack of reliable measured non-gravitational accelerations in Swarm A and Swarm B, the three-satellite Swarm GFM considers the ASU modelled non-gravitational accelerations for these satellites. For Swarm C, we consider three cases where the non-gravitational accelerations are either measured or represented by TUD or ASU's model. In this way, we isolate the effect of the three types of non-gravitational acceleration data. The results for January 2015 are shown in Fig. 5, using ASU and TUD models and calibrated accelerometer data.

Table 6Geoid height difference in millimetres between Swarm and GRACE GFMs.

Download Print Version | Download XLSX

The three-satellite solutions that use modelled non-gravitational accelerations in Swarm C are remarkably similar (see Fig. 5). In spite of this, note that using accelerometer data improved the agreement to GOCO05S for degrees 2 and 4.

To gather a better overview of the added value of the three types of non-gravitational accelerations, we derive the following model difference D, similar to rms:

(1) D = median Δ h 2 + MAD Δ h 2 ,

with the median absolute deviation (MAD) analogous to SD when the median is considered instead of the mean and Δh being the 1×1 geoid height difference between the 500 km Gaussian filtering three-satellite Swarm models and both ITSG-Grace2016 and GOCO05S, in the latitude band 85 from the Equator. We note that similar results were obtained using the CSR RL05 GRACE monthly solutions (not shown). The resulting differences are shown in Table 6.

Figure 5Three-satellite Swarm gravity field solutions using TUD- and ASU-modelled non-gravitational accelerations, as well as measured non-gravitational accelerations for Swarm C and ASU-modelled non-gravitational accelerations for Swarm A and B (January 2015).


The 2015 results indicate that the observed non-gravitational accelerations improve the agreement between the three-satellite Swarm models and ITSG-Grace2016 and GOCO05S, while that is not the case for 2016 (except for January 2016 and GOCO05S, when the GFM derived from TUD-modelled non-gravitational accelerations agrees equally well with the one derived considering observed non-gravitational accelerations). The comparison with GOCO05S intends to predict how possible it would be to assess the added value of the different types of non-gravitational accelerations during those periods when there are no GRACE data. Other time-dependent models were tested, but those do not agree as closely with GRACE monthly models (not shown).

The statistics in Table 6 imply that observed non-gravitational accelerations are only beneficial when the amplitude of the non-gravitational accelerations is larger than what was observed in 2016. This is likely related to the decreasing level of solar activity, which is approaching the minimum of its 11-year cycle (expected to reach the minimum in 2019). Through the influence of the solar radiation on the atmospheric density and resulting atmospheric drag, the low level of solar activity has a direct impact on the accelerometer measurements. The closer to the solar cycle minimum, the lower magnitude and variability of the accelerometer signal are. Another factor may be a potential worse performance of the accelerometer calibration procedure under low levels of solar activity, resulting from the lower signal-to-noise ratio (SNR) in the accelerometer data. In other words, the noise and (potentially) uncorrected artefacts in the accelerometer data of Swarm C are substantial enough to limit the usefulness of these data to gravimetric studies, except when the solar activity is high (as was the case in 2015) or when the satellites' altitude decays in the future. Given these characteristics and the continuing solar minimum, our Swarm models are not processed considering Swarm C accelerometer observations, but we plan to revisit this issue once the solar activity increases.

3.3 Individual Swarm models

In this section we illustrate the quality of the individual Swarm solutions. As described in Sect. 2.6, we directly compare Swarm at the epochs defined by GRACE under 750 km radius Gaussian smoothing.

Figure 6Time dependence between December 2013 and September 2019 of the epoch-wise cumulative degree amplitude (or global spatial rms) of the individual Swarm solutions, CSR RL06 GRACE, and their difference, considering 750 km smoothing.


Figure 6 shows a measure of the evolution of the quality of the individual Swarm solutions over the complete Swarm data period. We also plot the cumulative degree amplitude of GRACE, to illustrate the global spatial amplitude of the geophysical processes represented by these data. There is a clear improvement in the agreement of Swarm with GRACE, from rms differences as high as 40 cm geoid height in early 2014, down to 10 cm and below since 2016. We attribute this increase in quality to the decrease in solar activity and to the upgrades in the Swarm GPS receivers between 2015 and 2016 (van den IJssel et al., 2016; Dahle et al., 2017). As demonstrated in Sect. 3.4, the Swarm models contain large errors in the ocean areas, which dominate the global spatial rms difference; over land areas, the agreement with GRACE is much better.

The various individual solutions show different levels of quality. Generally speaking, the solutions from AIUB, ASU, and IfG cluster together as agreeing better with GRACE, with their dispersion narrowing down after 2016. This suggests that these approaches suffer differently in conditions of high solar activity, with ASU's models being the least sensitive overall. Possibly, ASU's efforts to minimize the amplification of the high frequencies when performing the double differentiation of the kinematic positions has the side effect of suppressing the negative effects of the high solar activity in the quality of the kinematic orbits. In contrast, OSU's solution consistently has lower agreement with GRACE. The velocity measurements, which are needed for IEBA (as well as any EBA-type approach), are to be derived from the kinematic positions by differentiation (then squared to obtain kinetic energy). The tedious data filtering and processing to approximate velocity errors is still imperfect, particularly in light of the spurious jumps in most of the kinematic orbits even in the cases without the GPS tracking signal degradation, e.g. from the southern Atlantic anomaly.

Figure 7Spectral agreement over the period between December 2013 and September 2019 in terms of the degree mean of the per-coefficient temporal rms difference, of the individual Swarm solutions, CSR RL06 GRACE, and their difference, considering 750 km smoothing.


Another way of analysing the agreement between the individual solutions and GRACE is to derive per-coefficient statistics of their temporal variations. One such statistic is the coefficient-wise temporal rms of the difference between the Swarm individual solutions and GRACE, thus producing a set of Stokes coefficients that describes the variability of that difference; from this set we compute the mean over each degree to represent the general agreement at the corresponding spatial wavelengths. The results are summarized in Fig. 7, which quantifies the agreement of Swarm and the GRACE climatological model in the spectral domain. Note that for most individual solutions, the rms difference decreases with degree as a result of the Gaussian smoothing, without which the curves would have a strong overall positive slope.

The ranking of quality of the individual solutions changes with spatial wavelength; for example, although OSU's solutions are consistently worse than IfG's as shown in Fig. 6, their degrees 2 and 3 are on average in equal or better agreement with GRACE. This diversity in the particularities of the various solutions is the main motivation for our practice to combine solutions derived from multiple gravity field estimation approaches. Unfortunately, as explained in Sect. 2.5, our combination is done at the solution level with weights derived from VCE, which means we loose the ability to weight the individual solutions differently in the degree domain and we cannot fully take advantage of the per-degree variations in quality of the individual solutions. Nevertheless, the VCE weights produce combined solutions with better agreement with GRACE than those combined at the NEQ level (Teixeira da Encarnação and Visser, 2019). From this we interpret that the benefits from per-degree weighting may not be as significant as the disadvantages of the combination at NEQ level, namely the different types of formal and calibrated errors, their different temporal evolution, and the difficulty in finding adequate empirical weights.

Figure 8Spectral correlation over the period between December 2013 and September 2019 in terms of the degree mean of the per-coefficient temporal correlation coefficient, of the individual Swarm solutions, CSR RL06 GRACE, and their difference, considering 750 km smoothing.


Finally, Fig. 8 shows the correlation of the Swarm time series with GRACE for the relevant spatial wavelengths. This figure is complementary to Fig. 7, since it does not illustrate the overall agreement (which is a measure of error) but the level that Swarm observes the same temporal evolution as GRACE (i.e. if Swarm sees the same proportional mass increase or decrease as GRACE). Understandably, the highest correlations correspond to the lowest degrees, not only because those are the signals with the highest amplitude (and therefore better observed by Swarm and GRACE) but also because of the smoothing. There is no obvious individual solution that stands out as being better correlated with GRACE, although ASU has the highest correlation coefficient for degrees 2, 4, and 7 to 9, while for AIUB that is the case for degrees 3 and 4. OSU's solution tends to correlate the least, except for degrees 4 and 8; this again indicates that a solution that may at first seem to be of consistently inferior quality may still provide a positive contribution to the combination. Also note that the correlations drop below 0.1 above degree 14 and remain relatively constant for higher degrees, indicating there is very little signal in the individual solutions that represent the same temporal variations as GRACE.

Figure 9Time dependence between December 2013 and September 2019 of the epoch-wise cumulative degree amplitude (or global spatial rms) of the combined Swarm solutions, CSR RL06 GRACE, and their difference, considering 750 km Gaussian smoothing.


3.4 Combined Swarm models

Having presented the individual Swarm GFMs in the previous section, we dedicate the current section to the analysis of the combined solutions. For more details about the combination strategy, refer to Teixeira da Encarnação and Visser (2019) and Meyer (2020). We determine the necessary intensity of smoothing of the Swarm models (Sect. 3.4.1) and illustrate the different sensitivity of the Swarm data to observe mass transport processes over land and ocean areas (Sect. 3.4.2).

3.4.1 Smoothing of the Swarm solutions

As demonstrated by Teixeira da Encarnação et al. (2016), the Swarm models do not seem to be sensitive to full wavelengths shorter than roughly 1500 km. We now update this assessment in light of the much longer time series and improved combination strategy than was the case in earlier publications. We compute the cumulative degree amplitude (which is proportional to the global spatial rms) of the difference between the Swarm and GRACE models and the unsmoothed GRACE climatological model, for two levels of smoothing radii: 750 km (Fig. 9) and 1500 km (Fig. 10).

For the 750 km case, the Swarm difference nearly always has the same amplitude as the Swarm signal itself. We will demonstrate in Sect. 3.4.2 that a significant portion of the amplitude of the Swarm difference is located over ocean areas and the agreement over land is significantly better. In spite of this, the lower amplitude of Swarm relative to the GRACE data suggests this smoothing intensity is inadequate to isolate the geophysical signal in the Swarm time series at the global scale.

Figure 10Time dependence between December 2013 and September 2019 of the epoch-wise cumulative degree amplitude (or global spatial rms) of the combined Swarm solutions, CSR RL06 GRACE, and their difference, considering 1500 km Gaussian smoothing.


In the case of 1500 km smoothing, the Swarm differences have comparable amplitudes to the GRACE data since mid-2015. We interpret this observation, given the conservative nature of the Swarm global rms difference, as indication that there is unnecessary suppression of the signal at spatial wavelengths from the two smoothing intensities considered in Figs. 9 and 10 (roughly 1500 to 3000 km, since we report smoothing radii).

We also repeated this exercise for the cases of no smoothing and 300 km smoothing radius. Those results indicated that the errors above degree 12 dominate the solution and produce monthly differences of negligible difference relative to the full Swarm spatial variability (not shown).

3.4.2 Land and deep-ocean signal

This section illustrates the differences in SNR characteristic of the Swarm GFMs by computing separate statistics for land and deep-ocean areas, the latter defined in Sect. 2.6.2.

In Fig. 11 the rms of the deep-ocean areas is shown in terms of the difference between the Swarm and GRACE solutions. As expected, the GRACE GFMs over the oceans have a relatively small amplitude, well under 2 cm EWH. Additionally, the Swarm GFMs show differences which are of much higher amplitude than the ocean signal represented by the GRACE data and barely different than the magnitude of the Swarm signal itself. In other words, Swarm is unable to resolve any monthly ocean signal with the spatial scales that GRACE can observe; however, the same cannot be said about (i) long-term trends since the data were de-trended prior to computing the statistics in Fig. 11 or (ii) aggregate measures such as the global ocean mass reported by Lück et al. (2018).

Figure 11Time dependence between December 2013 and September 2019 of the epoch-wise cumulative degree amplitude (or global spatial rms) of the combined Swarm solutions, CSR RL06 GRACE, and their difference, for deep-ocean areas considering 750 km smoothing.


We illustrate the agreement of Swarm and GRACE solutions in the spectral domain in Fig. 12 (which is produced in a similar way as Fig. 7). Similar to the evolution of the temporal agreement represented in Fig. 11, the spectral analysis illustrates that Swarm differs from the climatological model with amplitudes that surpass the signal, across all the spatial wavelengths, over the oceanic areas. The only exception refers to degree 2 but that is mainly driven by the consistent use of C20 published in Cheng and Ries (2018).

Figure 12Spectral agreement over the period between December 2013 and September 2019 in terms of the degree mean of the per-coefficient temporal rms difference, of the combined Swarm solutions, CSR RL06 GRACE, and their difference, for deep-ocean areas considering 750 km smoothing.


When it comes to land areas, the Swarm solutions agree with the climatological model much better than in the oceans. Figure 13 shows that since 2016, the Swarm difference with respect to the GRACE data has comparable amplitude. This means that Swarm is generally able to observe the majority of mass transport processes described by the climatological model (under Gaussian smoothing with 750 km radius), in particular after 2016. Prior to mid-2015, this is on average not the case, although we will demonstrate in Sect. 3.5.3 that regions where the mass transport signal is of substantial amplitude are reasonably well observed.

Figure 13Time dependence between December 2013 and September 2019 of the epoch-wise cumulative degree amplitude (or global spatial rms) of the combined Swarm solutions, CSR RL06 GRACE, and their difference, for land areas considering 750 km smoothing.


The analysis in the spectral domain summarized in Fig. 14 confirms that the difference with respect to the climatological model is of smaller amplitude than the signal therein represented up to degree 12. This result further confirms the result of Sect. 3.4.1 regarding de adequacy of smoothing the Swarm solutions with a Gaussian filter with 750 km radius.

Figure 14Spectral agreement over the period between December 2013 and September 2019 in terms of the degree mean of the per-coefficient temporal rms difference, of the combined Swarm solutions, CSR RL06 GRACE, and their difference, for land areas considering 750 km smoothing.


The results presented in Figs. 11–14 illustrate that the Swarm GFMs are unable to resolve the gravity signal in the oceanic regions at spatial lengths comparable to land areas. We observe that the discrepancy with respect to GRACE over the ocean is roughly 25 % larger than over land. We do not have a definitive explanation for this, other than the ionospheric activity may corrupt the estimated gravity field parameters over the oceans more significantly since away from land areas there is very little gravity signal to capture. In other words, the natural gravity variations over land are of sufficient amplitude to dominate the errors, at least enough to drive our statistics.

The higher accuracy over land could be explained by the ionospheric activity affecting mainly ocean areas, since those are mostly located along the Equator (e.g. the Pacific ocean). Masking the land areas could therefore remove the large land signals associated with hydrology and leave mostly the errors in the equatorial oceans. To test this hypothesis, we masked the Swarm and GRACE residual along tropical and non-tropical regions, as illustrated in Fig. 15.

Figure 15Time dependence between December 2013 and September 2019 of the epoch-wise cumulative degree amplitude (or global spatial rms) of the combined Swarm solutions masked over different regions, considering 750 km Gaussian smoothing.


It is clear that Swarm observes the tropical regions, which include regions with strong gravitational variations such as the Amazon basin and vast ocean areas in the Pacific, in as good agreement as the non-tropical regions. We note that the deep-ocean regions are not complementary of the land regions (i.e. the two domains do not cover the whole Earth; see Sect. 2.6.2) and it should not be expected that their spatial rms is proportionally larger than the tropical and non-tropical regions, which are of comparable amplitude between themselves and complementary.

We now focus on the necessary smoothing to retrieve any deep-ocean signal from the monthly Swarm models. We increased the smoothing intensity relative to what is discussed in Sect. 3.4.2 to demonstrate the capabilities of Swarm to contribute to ocean studies, in particular those related to large-scale mean dynamic ocean topography. For the case of global ocean mass Lück et al. (2018) already demonstrated an agreement with GRACE of less than 5 mm in terms of EWH. We tested smoothing radii of 1000, 1500, 3000, and 5000 km; the results for 3000 km are presented in Figs. 16 and 17.

Figure 16Time dependence between December 2013 and September 2019 of the epoch-wise cumulative degree amplitude (or global spatial rms) of the combined Swarm solutions, CSR RL06 GRACE, and their difference, for deep-ocean areas considering 3000 km smoothing.


Figure 17Spectral agreement over the period between December 2013 and September 2019 in terms of the degree mean of the per-coefficient temporal rms difference, of the combined Swarm solutions, CSR RL06 GRACE, and their difference, for deep-ocean areas considering 3000 km smoothing.


Figure 16 demonstrates that a smoothing radius of 3000 km is enough to reduce the spatial rms of the Swarm residual to amplitudes comparable to the signal in GRACE, particularly after 2016. This means that since 2016 Swarm has been observing ocean mass changes at the extremely coarse spatial scale of roughly 6000 km.

We further demonstrate Swarm's ability to resolve large-scale ocean mass changes in the spectral domain, Fig. 17. As illustrated, the smoothing radius of 3000 km is barely enough to, on average, decrease the degree average of the per-degree rms difference below the GRACE signal amplitude. Note that the spectral measure represented by the degree average considered the complete Swarm period, including the start of the mission, when the quality of the solutions was the lowest. Therefore, the smoothing radius of 3000 km is adequate to resolve large-scale Swarm deep-ocean mass changes since mid-2015.

Figure 18Per-coefficient correlation coefficient between the GRACE climatological model and Swarm.


3.5 Signal content

This section describes the geophysical signal represented by the Swarm models. We start by illustrating the time series of a few low-degree coefficients in Sect. 3.5.1. The variability of the Swarm model, and the patterns therein, is discussed in Sect. 3.5.2. We end with Sect. 3.5.3, where we give an overview of the capabilities of the Swarm models to observe large basin storage variations and how they compare to GRACE and GRACE-FO.

3.5.1 Low degrees

We now present the time series of a selection of low degree coefficients, without any smoothing applied. This section aims at illustrating the noise characteristics of the Swarm models in the time domain and how they compare to GRACE.

We give an overview of the per-degree correlation coefficients of Swarm and GRACE relative to the climatological model. The degree 2 coefficients (except C2,0), which are particularly important for sea-level studies, are subsequently presented. Finally, we show the selected case of C5,-1 that has an interesting temporal evolution and how Swarm and GRACE capture those signals. The time series of the zonal coefficients from degrees 3 to 5 are presented in Appendix D. Note that we represent the sine Stokes coefficients with negative order, e.g. C2,-1.

Figures 18 and 19 represent the correlation coefficient of the time series of Swarm and GRACE relative to the climatological model, including the early period of the mission when the quality of the Swarm models was lower. As expected, GRACE's coefficients correlated much more closely to the climatological model, as represented by the numerous dark red pixels in the triangular plot of Fig. 19. The overview of Swarm's correlation with the climatological model (Fig. 18) is dominated by values of around 0.2 (represented by a yellow colour), with some regions with average correlations of roughly 0.6 (represented by the red colour), notably for orders −5 to −3 and degrees 9 to 4. Furthermore, we observe some interesting common features in both Swarm and GRACE correlation plots, namely order −6 and C5,5 seem to be poorly captured by the climatological mode, since neither Swarm nor GRACE correlated well.

Figure 19Per-coefficient correlation coefficient between the GRACE climatological model and GRACE.


Figure 20Per-coefficient rms of the difference between the GRACE climatological model and Swarm.


Figure 20 illustrates one particularity of the Swarm models. The rms of the difference relative to the GRACE climatological model is heavily order dependent, with the even orders showing a larger rms than the odd orders (for degrees 4 and above); this effect is particularly striking for orders 6 and −6, as well as for 5 and −5. This feature is also present in the individual models (not shown), in spite of the fact that no order dependency is present in their formal errors. We cannot find an explanation for the discrepancy between the rms difference in even and odd orders.

Figures 21–24 show the time series of the degree 2 coefficients. They illustrate the general characteristics of Swarm coefficient time series: large signal amplitudes, in particular before mid-2015, as well as a general agreement in the average value, if one could imagine a heavy temporal smoothing operation. The last characteristic, which is extremely common for all the coefficients we have analysed (up to degree 6, not all shown here), finds a rare exception in C2,1, particularly before 2017. A possible explanation is related to the mean pole model (Wahr et al., 2015), which differs between our Swarm solutions (Appendix C) and CSR RL06 (Bettadpur, 2018).

Figure 21Coefficient C2,2 as observed by GRACE and Swarm, as well as represented by the GRACE climatological model.


Figure 22Coefficient C2,1 as observed by GRACE and Swarm, as well as represented by the GRACE climatological model.


Figure 23Coefficient C2,-1 as observed by GRACE and Swarm, as well as represented by the GRACE climatological model.


Figure 24Coefficient C2,-2 as observed by GRACE and Swarm, as well as represented by the GRACE climatological model.


Regarding the agreement of the temporal signal captured by Swarm and that captured by GRACE, it is generally possible to observe that Swarm tends to follow roughly in the same direction, albeit with large month-to-month changes (i.e. larger errors) and with frequent over-shootings before 2016. The large errors are the result of the Swarm solutions exploiting the less accurate kinematic positions as gravimetric observations (in comparison to the much more accurate inter-satellite ranges (ISRs) of GRACE). The errors tend to be larger before 2016, during the period of higher solar and ionospheric activity as well as prior to the GPS receiver tracking loop updates (van den IJssel et al., 2016; Dahle et al., 2017).

Figure 25 shows a representative case of a good agreement between Swarm and GRACE. The overall trend of the C5,-1 coefficient is well represented in the climatological model but fails to capture the abnormal deviation around early 2016, which is observed in a consistent way by GRACE and Swarm.

Figure 25Coefficient C5,-1 as observed by GRACE and Swarm, as well as represented by the GRACE climatological model.


3.5.2 Signal variability

The current section is devoted to presenting the signal variability in the Swarm solutions, shown in Fig. 26. The most striking features in the Swarm variability concern the strong geomagnetic Equator signature and the artefacts near the south magnetic pole (which is located due south of Tasmania, on the coast of Antarctica). Interestingly, there is no obvious signature close to the north magnetic pole (located north of Hudson Bay, west of Greenland). The geomagnetic Equator signature extends over land and ocean areas, notably the Sahara (somewhat less intensively), although it is possible to distinguish the signature of the strong geophysical signal over the Amazon basin. This artefact is also characterized by an obvious east–west-banded structure, which is well delineated over the central Atlantic, North Africa, mainland Southeast Asia, and western Pacific regions. In spite of these artefacts, we will demonstrate that Swarm is able to resolve monthly large-scale mass transport processes. For that purpose we look at the regions circumscribed by the red dashed rectangles in Fig. 26. We choose these regions because they are located at various geographical locations, are of different sizes, and are under the influence of different types of geophysical signals.

Figure 26Signal variability for Swarm during the period between December 2013 and March 2019, under 750 km Gaussian smoothing.

Figure 27Signal variability for GRACE during the Swarm period December 2013 to March 2019, including the earliest GRACE-FO solutions, under 750 km Gaussian smoothing.

Looking at the variability in the GRACE models over the same periods, Fig. 27 (produced in a way consistent with Fig. 26), there is no obvious signature of geomagnetic effects. Additionally, the variability over the oceans is very small, in comparison to land areas.

3.5.3 Large storage basins

In this section, we present time series of Swarm and GRACE average EWH over the areas highlighted in Figs. 26 and 27. Unlike the previous sections, the GRACE signal (relative to GGM05C) is calculated from the monthly RL06 CSR solutions, after 750 km smoothing and the usual C2,0 replacement. The trend (and bias) is co-estimated with yearly and semi-yearly sine and cosine periods, in order to be insensitive to phase differences at the beginning and end of the period under analysis. Instead of disclosing the constant term in the polynomial and sinusoidal regression, we prefer to report the average over the period under analysis as a measure of a constant bias.

We illustrate these time series with the example of Greenland and the Amazon, in Figs. 28 and 29, respectively. The remaining time series can be found in Appendix E. As was the case with the analysis of the low degrees, the time series are less smooth than GRACE, as a result of the increased influence of errors. In spite of this, the Swarm time series follows GRACE closely, with a correlation coefficient of 0.83 and 0.95 for Greenland and the Amazon, respectively. The trend is overestimated by 0.6 and 1.12 cm yr−1, respectively, mainly as a result of the higher errors before mid-2015. Swarm also agrees with the GRACE-FO observation that the Greenland ice mass loss slowed down during the winter of 2018–2019, since both Swarm and GRACE-FO lines are above the linear interpolation. During the summer of 2019, the ice mass loss in Greenland accelerated, which is consistently observed by both Swarm and GRACE. In the case of the Amazon basin, the GRACE-FO months agree particularly well with Swarm.

Figure 28Time series of EWH for the western Greenland region (latitude 60 to 85, longitude −60 to −37).


Table 7Bias and trend agreement, as well as correlation coefficient between the GRACE and Swarm time series for the selected basins, over the complete Swarm period (December 2013 to September 2019).

Download Print Version | Download XLSX

Table 7 provides an overview of the statistics derived from the time series of all analysed basins. The Swarm and GRACE time series agree on their average values between 1.50 cm (Amazon) and 0.78 cm (Orinoco), on their trend between 1.16 cm yr−1 (Orinoco) and 0.36 cm yr−1 (Congo and Zambezi) and on their correlation coefficient between 0.65 (Volga) and 0.95 (Amazon). All regions show a variety of values in their statistics, thus making it difficult to immediately identify which one is best observed. For example, although the Amazon time series shows the largest trend difference (in absolute value), it also has the highest correlation coefficient. If the period before mid-2015 is ignored, these statistics improve substantially (not shown).

Figure 29Time series of EWH for the Amazon basin (latitude −17 to 3, longitude −76 to −47).


Over the nine basins presented in this section and in Appendix E, the Swarm rms difference with respect to GRACE is 0.91 cm in terms of temporal mean and 0.76 cm yr−1 in terms of trend and shows an average correlation coefficient of 0.79 (bottom of Table 7). Note that the complete Swarm period was considered in deriving these statistics and represents a conservative estimate of the accuracy of Swarm if the early period before mid-2015 is discarded.

4 Data availability

The Swarm monthly models are distributed on a quarterly basis at ESA's Earth Swarm Data Access (at, last access: 5 June 2020, follow “Level2longterm” and then “EGF”) and at the International Centre for Global Earth Models (, last access: 5 June 2020), as well as identified with the DOI (Encarnacao et al., 2019).

5 Conclusions

We present Swarm GFMs resulting from the combination of four individual solutions computed from different gravity field estimation approaches: celestial mechanics approach (CMA), decorrelated acceleration approach (DAA), improved energy balance approach (IEBA), and short-arc approach (SAA). Two approaches (CMA and IEBA) exploit the KO solutions produced at AIUB and the other two (DAA and SAA) the KOs produced at IfG. The combination is done at the solution level, weighted by VCE; for the sake of brevity, we refer to Teixeira da Encarnação and Visser (2019) to demonstrate that our combination produces Swarm models in better agreement with GRACE than if the combination is done at the NEQ level.

We test the added value of KB in the quality of the Swarm GFM, when compared to the long wavelength signal recovered by GRACE, by computing seven GFMs during periods of different data quality. We demonstrate that the largest changes in the results appear during early 2015 (high ionospheric activity, before improvements in Swarm's GPS receivers) that translate into a slight deterioration of the quality of the Swarm solutions; see Table 5. For the 5 months analysed in 2016, considering two KO solutions, any improvement is minimal (0.1 to 0.2 mm geoid height, in five cases), negligible (in two cases), or slightly worse (0.1 mm in three cases). We conclude that any common errors that would be eliminated in the KB solutions are already (mostly) corrected in the KOs. For this reason, our Swarm GFMs do not consider KBs.

Another test regarding the added value of additional data took the form of including Swarm C non-gravitational accelerations. We compared the three-satellite Swarm solution produced considering the DAA and non-gravitational accelerations acting on Swarm C represented by the TUD and ASU non-gravitational acceleration models, in addition to exploiting the accelerometer measurements. Since the Swarm A and B satellites do not produce usable accelerometer readings, they are represented by the ASU model exclusively. The results indicate that the accelerometer observations are only beneficial in those cases when the amplitude of the non-gravitational accelerations acting on Swarm C are of higher amplitude than in quiet periods in solar activity, such as is the case since 2016. This may be the result of the potentially lower quality of the calibrated accelerations, caused by the lower SNR in the accelerometer observations.

Regarding the topic of non-gravitational accelerations in the processing of GPS-driven GFMs, we would like to comment on the results of Ditmar et al. (2006) and Ditmar et al. (2008), who demonstrated that non-gravitational accelerations are not needed for gravity field estimation and the quality of the GPS observations (and the resulting KOs) are the main drivers of the quality of the GFMs. Within our project, each AC is free to elect whatever processing strategy they deem to be most beneficial to their individual solutions, which is assessed internally. For example, AIUB has determined that the use of daily and 15 min piecewise-constant empirical parametrization does not require any modelling of non-gravitational accelerations. In the case of ASU, who exploit a dedicated decorrelation procedure (which is a frequency-dependent noise whitening procedure), their solutions benefit from drag, EARP, and EIRP models. Essentially, the inclusion of frequency-dependent data weighting (FDDW) is not within immediate reach to all ACs, in which case other processing strategies seem to produce comparable solution quality. In summary, we do not wish to contest previous results on this topic but clarify the differences in our processing choices.

We quantify the different quality of the various individual solutions and demonstrate that all have the potential to contribute positively to the quality of the combined Swarm time series of GFMs. We additionally explain that our approach to combine the individual GFMs at the solution level considering VCE weights is an effective way of overcoming the difficulty in combining solutions at the NEQ level when the corresponding normal matrices represent errors of different type, formal and calibrated in our case (Teixeira da Encarnação and Visser, 2019).

Masking the Swarm data separately over ocean and land, we demonstrate Swarm's combined models' ability to measure land mass transport processes, with a comparable spatial variability of the Swarm–GRACE residual to that of the GRACE signal, for the post mid-2015 period. Over the ocean areas, the spatial rms of the difference between Swarm and GRACE is always larger than the spatial rms of the latter. To resolve the oceanic signal, the Swarm data required a more intense Gaussian smoothing, with a radius of 3000 km.

We analyse the signal content of the Swarm models in terms of time series of the low degrees, spatial patterns of the temporal signal variability, and time series of large storage basins. Comparing the time series of isolated SH coefficients, we show that the Swarm data generally have a more erratic temporal evolution with sudden month-to-month variations. We attribute this particularity to the lower accuracy of the GPS observations as gravimetric data, compared to GRACE's K-band ranging (KBR) data. We also illustrate features in the Swarm data that are not captured by the GRACE climatological model, but confirmed by the GRACE and GRACE-FO data, notably the atypical deviation around early 2016 in the C5,-1 coefficient. By plotting the spatial patterns in the temporal variability of the Swarm data, we bring into evidence the strong signature over the geomagnetic Equator, showing strong meridional stripes, and over the south magnetic pole (but not on the north magnetic pole). In spite of this artefact, the strong mass variability over the Amazon basin is clearly visible. With regard to the time series of mass changes over large storage basins, Swarm agrees on average with GRACE (the climatological model was not relevant to this analysis) at 0.91  cm in terms of temporal mean, 0.76 cm yr−1 in terms of trend, and 0.79 correlation coefficient over the nine basins we considered. We show that Swarm agrees with the observation of GRACE-FO that the ice mass loss over Greenland seems to have slowed down during late 2018 and accelerated in the summer of 2019, in spite of the heavy signal dilution caused by the necessary smoothing to reduce the errors in the Swarm models. Although our Swarm models are already in production mode, we are considering several options to improve their quality. Given the high sensitivity of the KOs to ionospheric activity, we plan to focus our efforts to improve the weighting of the GPS observations (Dahle et al., 2017; Kermarrec et al., 2018; Schreiter et al., 2019). We also plan to decrease the disagreement between the individual solution produced at OSU and those at other ACs by including advanced algorithms for reducing the effects of jumps and the amplification of high-frequency noise in the differentiation of the KO positions into velocities.

Appendix A: Kinematic orbits

A1 Delft University of Technology

Software: GPS High precision Orbit determination
Software Tool (van Helleputte, 2004; Wermuth et al., 2010)
Differencing scheme: Undifferenced
Linear combination: Ionosphere-free
GPS observations: Code and carrier phase
Estimator: Bayesian-weighted LS
Arc length: 30 h
Data weighting: a priori weights equal to 1 m and 1 mm for code and phase observations (resp.)
Transmitter PCV: Official IGS08 ANTEX up to day 17/028, official IGS14 ANTEX from day 17/029 on
Receiver PCV: Empirically determined from stacking of reduced-dynamic POD residuals with 1 binning
Data screening: Minimum SNR of 10, minimum of six GPS satellites, code and phase outlier editing threshold of
2 m and 3.5 cm, respectively, 1 m or larger difference between estimated KO positions and with
reduced-dynamic precise science orbit (PSO)
Earth precession model: International Astronomical Union (IAU) 1976 (Lieske et al., 1977)
Earth nutation model: IAU 1980 (Seidelmann, 1982)
Earth orientation model: Centre for Orbit Determination in Europe (CODE) final Earth rotation parameters (ERP)

A2 Astronomical Institute of the University of Bern

Software: Bernese v5.3 (Dach et al., 2015; Jäggi et al., 2006)
Differencing scheme: Undifferenced
Linear combination: Ionosphere-free
GPS observations: Carrier phase
Estimator: Batch LS
Arc length: 24 h
Data weighting: Not applicable (n/a)
Transmitter PCV: Official IGS08 ANTEX up to day 17/028, official IGS14 ANTEX from day 17/029 on
Receiver PCV: Stacking of residuals from reduced-dynamic precise orbit determination (POD) of
approx. 120 d, nine iterations, 1 binning
Data screening: 2 cm s−1 or larger time differences of the geometry-free linear combination of
L1B GPS carrier phase observations
Earth precession model: International Earth Rotation Service (IERS) 2010 Conventions (Petit and Luzum, 2010)
Earth nutation model: IERS 2010 Conventions (Petit and Luzum, 2010)
Earth orientation model: CODE final ERP

A3 Institute of Geodesy Graz

Software: GROOPS
Differencing scheme: None
Linear combination: None (the ionospheric influence is co-estimated)
GPS observations: Code and carrier phase
Estimator: LS
Arc length: 24 h
Data weighting: Elevation and azimuth-dependent, epoch-wise VCE
Transmitter PCV: Empirical, estimated from 5.5 years of data, including data from several LEO missions (GRACE,
Jason 2 & 3, MetOp-A & MetOp-B, Sentinel 3A, Swarm, TanDEM-X, TerraSAR-X) (Zehentner, 2016)
Receiver PCV: Empirical, spherical harmonics (maximum D/O 60), derived from 38 months of data
Data screening: Implicit in VCE
Earth precession model: IAU 2006/2000A precession–nutation model (Petit and Luzum, 2010)
Earth nutation model: IAU 2006/2000A precession–nutation model (Petit and Luzum, 2010)
Earth orientation model: IERS Earth Orientation Parameter (EOP) 08 C04 (Petit and Luzum, 2010)

A4 Common

Carrier phase ambiguities: Float
Receiver clock corrections: Co-estimated
Sampling rate: 10 or 1 s (depending on L1B GPS data)
Elevation cut-off angle: 0
GPS orbits and clocks: Final orbits and 5 s clocks of Centre for Orbit Determination in Europe (Dach et al., 2017)
Swarm attitude: L1B attitude data
Appendix B: Kinematic baselines

B1 Delft University of Technology

Software: Multiple-satellite orbit determination using Kalman filtering (van Barneveld, 2012)
Linear combination: n/a (the ionospheric frequency-dependent influence is modelled)
Estimator: Iterative EKF
Carrier phase ambiguities: Integer, using the least-squares ambiguity de-correlation adjustment method (Teunissen, 1995)
Receiver PCV: Empirical phase centre variations (PCVs) and code residual variations (CRVs);
maps are estimated a priori for each GPS frequency

B2 Astronomical Institute of the University of Bern

Software: Bernese (Dach et al., 2015; Jäggi et al., 2006), development version
Linear combination: Ionosphere-free
Estimator: LS
Carrier phase ambiguities: Wide-lane and narrow-lane integer ambiguity fixing with the Melbourne–Wübbena combination
and the ionosphere-free linear combination, respectively
Receiver PCV: Empirical

B3 Common

Differencing scheme: Double-differenced
GPS observations: Code and carrier phase
Carrier phase ambiguities: Integer
Appendix C: Gravity field models

C1 Astronomical Institute of the University of Bern

Software: Bernese v5.3 (Dach et al., 2015; Jäggi et al., 2006)
Approach: celestial mechanics approach (Beutler et al., 2010)
Reference GFM: AIUB GRACE-only static model, version 3 (Jäggi et al., 2011)
Empirical parameters: Daily and 15 min piecewise-constant (constrained)
Drag model: None
EARP and EIRP models: None
Non-tidal model: Atmosphere and Ocean De-aliasing Level 1B (Flechtner, 2011)
Ocean tidal model: 2011 Empirical Ocean Tide model (Savcenko and Bosch, 2012)
Permanent tide system: Tide-free

C2 Astronomical Institute Ondřejov

Software: Developed in-house
Approach: Decorrelated acceleration approach (Bezdek et al., 2014)
Reference GFM: ITG GRACE-only static model, 2010 (Mayer-Gürr et al., 2010)
Empirical parameters: Daily constant-piecewise
Drag model: US Naval Research Laboratory Mass Spectrometer and Incoherent Scatter Radar atmospheric model
(Picone et al., 2002)
EARP and EIRP models: Knocke et al. (1988)
Non-tidal model: Atmosphere and Ocean De-aliasing Level 1B (Dobslaw et al., 2017)
Ocean tidal model: 2004 finite element solution (Lyard et al., 2006)
Permanent tide system: Tide-free

C3 Institute of Geodesy Graz

Software: GROOPS
Approach: Short-arc approach (Mayer-Gürr, 2006)
Reference GFM: GOCO release 05 satellite-only gravity field model (Mayer-Gürr, 2015)
Empirical parameters: Piecewise linear for each arc (ranging from 15 to 45 min)
Drag model: Jacchia–Bowman 2008 (Bowman et al., 2008)
EARP and EIRP models: Rodriguez-Solano et al. (2012)
Non-tidal model: Atmosphere and Ocean De-aliasing Level 1B RL06 (Dobslaw et al., 2017)
Ocean tidal model: 2014 finite element solution (Carrere et al., 2015)
Permanent tide system: Zero tide

C4 Ohio State University

Software Developed in-house
Approach: Improved energy balance approach (Shang et al., 2015)
Reference GFM: GRACE Intermediate Field 48 (Ries et al., 2011) up to degree and order (D/O) 200
Empirical parameters: Second-order polynomial every 3 h, 1-CPR sinusoidal every 24 h
Regularization: None
Drag model: US Naval Research Laboratory Mass Spectrometer and Incoherent Scatter Radar atmospheric model
(Picone et al., 2002)
EARP and EIRP models: Knocke et al. (1988)
Non-tidal model: Atmosphere and Ocean De-aliasing Level 1B (Flechtner, 2011)
Ocean tidal model: 2011 empirical ocean tide model (Savcenko and Bosch, 2012)
Permanent tide system: Tide-free

C5 Common

Atmospheric tidal model: Biancale and Bode (2006)
Solid Earth tidal model: IERS2010
Pole tidal model: IERS2010
Ocean pole tidal model: IERS2010
Third-body perturbations: Sun, Moon, Mercury, Venus, Mars, Jupiter, and Saturn, following the JPL Planetary and
Lunar Ephemerides (Folkner et al., 2014)
C2,0 coefficient: Estimated alongside other coefficients

Appendix D: Time series of zonal coefficients

Figures D1–D3 illustrate the time series for the zonal coefficients of degrees 3 to 6, respectively.

Figure D1Coefficient C3,0 as observed by GRACE and Swarm, as well as represented by the GRACE climatological model.


Figure D2Coefficient C4,0 as observed by GRACE and Swarm, as well as represented by the GRACE climatological model.


The zonal coefficient of degree 3 is an interesting case because both Swarm and GRACE-FO observe a phase shift during late 2018, relative to the climatological model, which is well in phase with GRACE for the non-GRACE-FO period (2003–2017). Swarm already captures this phase shift possibly as early as mid-2017, although the noisy character of the Swarm time series weakens this type of statement.

The zonal coefficient of degree 4 is one of the few examples where the Swarm time series shows a clear bias relative to GRACE and the climatological model, after 2017 in this case. As was the case for C2,1, we cannot explain such behaviour.

The zonal coefficient of degree 5 is an example of excellent agreement between all three time series. Swarm still shows the characteristic noise, as well as a higher overall disagreement before mid-2015. These are features intrinsic to our Swarm solutions.

Figure D3Coefficient C5,0 as observed by GRACE and Swarm, as well as represented by the GRACE climatological model.


Appendix E: Storage basin time series

Figure E1Time series of EWH for Alaska (latitude 56 to 65, longitude −151 to −129).


Figure E2Time series of EWH for the Congo and Zambezi basins (latitude −23 to −3, longitude 14 to 38).


Figure E3Time series of EWH for the Ganges–Brahmaputra basin (latitude 15 to 30, longitude 72 to 89).


Figure E4Time series of EWH for the northern Australia region (latitude −24 to −10, longitude 124 to 145).


Figure E5Time series of EWH for the Orinoco basin (latitude −3 to 12, longitude −72 to −59).


Figure E6Time series of EWH for the Volga basin (latitude 53 to 61, longitude 34 to 56).


Figure E7Time series of EWH for the western Antarctica region (latitude −80 to −70, longitude −140 to −85).


Appendix F: Abbreviations
3D Three-dimensional
AA Acceleration approach
DAA Decorrelated acceleration approach
AC Analysis centre
AIUB Astronomical Institute of the University of Bern, Switzerland
AIUB-GRACE03S AIUB GRACE-only static model, version 3
AOD1B Atmosphere and Ocean De-aliasing Level 1B
AOD1B-RL06 Atmosphere and Ocean De-aliasing Level 1B RL06
ANGARA Analysis of non-gravitational accelerations due to radiation pressure and aerodynamics
ASU Astronomical Institute (Astronomický ústav), AVCR, Ondřejov
AVCR Czech Academy of Sciences (Akademie ved Ceské Republiky), Czech Republic
CHAMP Challenging Mini-Satellite Payload
CODE Centre for Orbit Determination in Europe
CMA Celestial mechanics approach
CoM Centre of mass
COST-G Combination Service of Time-variable Gravity Fields
CPR Cycle per revolution
CRV Code residual variation
CSR Center for Space Research, UT Austin, USA
D/O Degree and order
DAA Decorrelated acceleration approach
DD Double-differenced
DISC Data, innovation and science cluster
DOI Digital object identifier
DWM07 Disturbance Wind Model 07
EARP Earth albedo radiation pressure
EGSIEM European Gravity Service for Improved Emergency Management, EU Horizon 2020
EIRP Earth infrared radiation pressure
EKF Extended Kalman filter
EBA Energy balance approach
EOT Empirical ocean tide model
EOT11a 2011 empirical ocean tide model
EWH Equivalent water height
EOP Earth orientation parameter
ERBE Earth Radiation Budget Experiment
ERP Earth rotation parameters
ESA European Space Agency
EU European Union
FDDW Frequency-dependent data weighting
FES Finite element solution global tide model
FES2004 2004 finite element solution
FES2014 2014 finite element solution
GFM Gravity field model
GHOST GPS High precision Orbit determination Software Tool
GGM05C Combined GRACE Gravity Model 05
GIF48 GRACE Intermediate Field 48
GOCE Gravity Field and Steady-State Ocean Circulation Explorer
GOCO Gravity Observation Combination
GOCO05S GOCO release 05 satellite-only gravity field model
GPS Global Positioning System
GRACE Gravity Recovery and Climate Experiment
GROOPS Gravity Recovery Object Oriented Programming System
hl-SST High-low satellite-to-satellite tracking
HWM07 Horizontal Wind Model 07
IAG International Association of Geodesy
IAU International Astronomical Union
ICGEM International Centre for Global Earth Models
IEBA Improved energy balance approach
IERS International Earth Rotation Service
IERS2010 IERS Conventions 2010
IfG Institute of Geodesy, TUG, Graz
IGFS International Gravity Field Service
IR Infrared radiation
ISR Inter-satellite range
ITG Institut für Geodäsie und Geoinformation, Germany
ITSG Institute of Theoretical Geodesy and Satellite Geodesy
ITG-Grace2010s ITG GRACE-only static model, 2010
ITSG-Grace2016 ITSG GRACE-only model, 2016
JB2008 Jacchia-Bowman 2008
JPL Jet Propulsion Laboratory, USA
JPL-PLE JPL Planetary and Lunar Ephemerides
KB Kinematic baseline
KBR K-band ranging
KO Kinematic orbit
L1A Level 1A data
L1B Level 1B data
L2 Level 2 data
L2PS Level 2 processing system
LAMBDA Least-squares ambiguity de-correlation adjustment
LEO Low-Earth orbit
ll-SST Low–low satellite-to-satellite tracking
LoS Line of sight
LS Least squares
MAD Median absolute deviation
MODK Multiple-satellite orbit determination using Kalman filtering
n/a Not applicable
NEQ Normal equation
NRLMSISE US Naval Research Laboratory Mass Spectrometer and Incoherent Scatter Radar atmospheric model
NRTDM Near Real-Time Density Model
OSU Ohio State University
PCV Phase centre variation
POD Precise orbit determination
PSO Precise or post-processed science orbit
RINEX Receiver Independent Exchange Format
RL05 Release 5
RL06 Release 6
Rms Root mean square
SAA Short-arc approach
SC Stokes coefficient
SH Spherical harmonic
SLR Satellite laser ranging
SNR Signal-to-noise ratio
SRP Solar radiation pressure
SD Standard deviation
TUD Delft University of Technology, Netherlands
TUG Graz University of Technology, Austria
UT Austin University of Texas at Austin
USA United States of America
VEA Variational equations approach
VCE Variance component estimation
ZD Zero-differenced
Author contributions

AB, ED, JIJ, EI, and JS produced the modelled non-gravitational data. ED, JIJ, and EI pre-processed the accelerometer data. PV, DA, ME, JIJ, EI, AJ, SK, XM, TMG, and CD produced and analysed the kinematic orbit and baselines. AB, ME, JG, AJ, JK, SK, TMG, UM, JS, CKS, CZ, YZ, and CD produced and analysed the gravity field models. JTE, DA, AB, ED, ME, EI, SK, XM, UM, JS, YZ, and CD took care of data curation. JTE, PV, DA, AB, XM, TMG, UM, and YZ analysed the results. JTE and PV coordinated the project activities and undertook funding acquisition. JTE, PV, DA, AB, ED, JIJ, AJ, XM, TMG, UM, CKS, and CD developed the relevant methodology for the production and analysis of orbits, baselines, accelerations, and gravity field models. JTE, PV, DA, AB, ED, XM, TMG, UM, and CKS wrote the article. All co-authors reviewed the manuscript.

Competing interests

The authors declare that they have no conflict of interest.


We thank the Swarm Data, Innovation and Science Cluster (DISC) for the support, flexibility, and guidance during the project activities. We also thank Mark Tamisiea for the fruitful discussions and critical viewpoints as well as Bryant Loomis for the weekly SLR-derived C20 time series. The authors acknowledge the reviewers for their valuable remarks that helped improve the quality of the original manuscript.

Financial support

This research was funded by the European Space Agency (contracts SW-CO-DTU-GS-111 and SW-CN-DTU-GS-027, part of contract 4000109587/13/I-NB), partially supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (grant XDA19070302), by the National Key Research & Development Program of China (grant 2017YFA0603103-3), by the National Natural Science Foundation of China (grant 41584016), and by the Ministry of Education, Youth and Sports of the Czech Republic (grant LTT18011).

Review statement

This paper was edited by Kirsten Elger and reviewed by two anonymous referees.


Abdalati, W., Gail, W. B., Busalacchi, A. J., Battel, S. J., Boland, S. W., Braun, R. D., Chen, S. S., Dietrich, W. E., Doney, S. C., Field, C. B., Fricker, H. A., Gille, S. T., Hartmann, D. L., Jacob, D. J., Janetos, A. C., Joseph, E., Macauley, M. K., Penner, J. E., Sorooshian, S., Stephens, G. L., Tapley, B. D., and Wilson, W. S.: Thriving on Our Changing Planet: A Decadal Strategy for Earth Observation from Space, National Academies Press, Washington, DC,, 2018. 

Allende-Alba, G., Montenbruck, O., Jäggi, A., Arnold, D., and Zangerl, F.: Reduced-dynamic and kinematic baseline determination for the Swarm mission, GPS Solutions, 21, 1275–1284,, 2017. 

Barkstrom, B. R. and Smith, G. L.: The Earth Radiation Budget Experiment: Science and implementation, Rev. Geophys., 24, 379–390,, 1986. 

Bettadpur, S.: UTCSR Level-2 Processing Standards Document For Level-2 Product Release 0006, Tech. rep., Center for Space Research, Austin, USA, available at: (last access: 5 June 2020), 2018. 

Beutler, G., Jäggi, A., Mervart, L., and Meyer, U.: The celestial mechanics approach: theoretical foundations, J. Geodesy, 84, 605–624,, 2010. 

Bezdek, A.: Calibration of accelerometers aboard GRACE satellites by comparison with POD-based nongravitational acceler- ations, J. Geodynam., 50, 410–423,, 2010. 

Bezdek, A., Klokocník, J., Kostelecký, J., Floberghagen, R., and Gruber, C.: Simulation of free fall and resonances in the GOCE mission, J. Geodynam., 48, 47–53,, 2009. 

Bezdek, A., Sebera, J., Klokocník, J., and Kostelecký, J.: Gravity field models from kinematic orbits of CHAMP, GRACE and GOCE satellites, Adv. Space Res., 53, 412–429,, 2014. 

Bezdek, A., Sebera, J., Teixeira da Encarnação, J., and Klokocník, J.: Time-variable gravity fields derived from GPS tracking of Swarm, Geophys. J. Int., 205, 1665–1669,, 2016. 

Bezdek, A., Sebera, J., and Klokocník, J.: Validation of Swarm accelerometer data by modelled nongravitational forces, Adv. Space Res., 59, 2512–2521,, 2017. 

Bezdek, A., Arnold, D., Doornbos, E., Jäggi, A., Mao, X., Zehentner, N., Teixeira da Encarnacao, J., and Visser, P. N. A. M.: TN-02: Swarm Data Pre-Processing, Kinematic Baselines And Accelerometer Data, Tech. rep.,, 2018a. 

Bezdek, A., Sebera, J., and Klokocník, J.: Calibration of Swarm accelerometer data by GPS positioning and linear temperature correction, Adv. Space Res., 62, 317–325,, 2018b. 

Biancale, R. and Bode, A.: Mean annual and seasonal atmospheric tide models based on 3-hourly and 6-hourly ECMWF surface pressure data, Tech. rep., Deutsches GeoForschungsZentrum GFZ, Potsdam, Germany,, 2006. 

Bowman, B., Tobiska, W. K., Marcos, F., Huang, C., Lin, C., and Burke, W.: A New Empirical Thermospheric Density Model JB2008 Using New Solar and Geomagnetic Indices, in: AIAA/AAS Astrodynamics Specialist Conference and Exhibit, August, American Institute of Aeronautics and Astronautics, Reston, Virigina,, 2008. 

Carrere, L., Lyard, F., Cancet, M., and Guillot, A.: FES 2014, a new tidal model on the global ocean with enhanced accuracy in shallow seas and in the Arctic region, in: EGU General Assembly, Vienna, Austria, 2015. 

Cheng, M. and Ries, J.: The unexpected signal in GRACE estimates of C20, J. Geodesy, 91, 897–914,, 2017. 

Cheng, M. and Ries, J.: GRACE Technical Note 11: Monthly estimates of C20 from 5 satellites based on GRACE RL06 models, available at: (last access: 5 June 2020), 2018. 

Cheng, M., Ries, J. C., and Tapley, B. D.: Variations of the Earth's figure axis from satellite laser ranging and GRACE, J. Geophys. Res.-Sol. Ea., 116, 1–14,, 2011. 

Dach, R., Lutz, S., Walser, P., and Fridez, P.: Bernese GNSS Software Version 5.2, Bern Open Publishing, Bern,, 2015. 

Dach, R., Schaer, S., Arnold, D., Prange, L., Sidorov, D., Susnik, A., Villiger, A., and Jäggi, A.: CODE final product series for the IGS,, 2017. 

Dahle, C., Arnold, D., and Jäggi, A.: Impact of tracking loop settings of the Swarm GPS receiver on gravity field recovery, Adv. Space Res., 59, 2843–2854,, 2017. 

Ditmar, P., Klees, R., and Liu, X.: Frequency-dependent data weighting in global gravity field modeling from satellite data contaminated by non-stationary noise, J. Geodesy, 81, 81–96,, 2006. 

Ditmar, P., Bezdek, A., Liu, X., and Zhao, Q.: On a Feasibility of Modeling Temporal Gravity Field Variations from Orbits of Non-dedicated Satellites, in: Observing our Changing Earth, edited by: Sideris, M., vol. 133, in:International Association of Geodesy Symposia, 307–313, Springer Berlin Heidelberg, Berlin, Heidelberg,, 2008. 

Dobslaw, H., Bergmann-Wolf, I., Dill, R., Poropat, L., Thomas, M., Dahle, C., Esselborn, S., König, R., and Flechtner, F.: A new highresolution model of non-tidal atmosphere and ocean mass variability for de-aliasing of satellite gravity observations: AOD1B RL06, Geophys. J. Int., 211, 263–269,, 2017. 

Doornbos, E.: Thermospheric Density and Wind Determination from Satellite Dynamics, Springer Theses, Springer Berlin, Heidelberg,, 2012. 

Doornbos, E., Bruinsma, S. L., Fritsche, B., Koppenwallner, G., Visser, P., Van Den IJssel, J., and Teixeira da Encarnação, J.: GOCE+ Theme 3: Air density and wind retrieval using GOCE data final report, Tech. rep., TU Delft, Delft, the Netherlands, 2014. 

Doornbos, E., Siemes, C., Teixeira da Encarnação, J., Perestý, R., Grunwaldt, L., Kraus, J., Holmdahl Olsen, P. E., van den IJssel, J., Flury, J., and Apelbaum, G.: Processing of Swarm Accelerometer Data into Thermospheric Neutral Densities, in: AGU Fall Meeting, Abstract SA31D-2371, 2015. 

Drob, D. P., Emmert, J. T., Crowley, G., Picone, J. M., Shepherd, G. G., Skinner, W., Hays, P., Niciejewski, R. J., Larsen, M., She, C. Y., Meriwether, J. W., Hernandez, G., Jarvis, M. J., Sipler, D. P., Tepley, C. A., O'Brien, M. S., Bowman, J. R., Wu, Q., Murayama, Y., Kawamura, S., Reid, I. M., and Vincent, R. A.: An empirical model of the Earth's horizontal wind fields: HWM07, J. Geophys. Res.-Space, 113, 1–18,, 2008. 

Emmert, J. T., Drob, D. P., Shepherd, G. G., Hernandez, G., Jarvis, M. J., Meriwether, J. W., Niciejewski, R. J., Sipler, D. P., and Tepley, C. A.: DWM07 global empirical model of upper thermospheric storm-induced disturbance winds, J. Geophys. Res.-Space, 113, A11319,, 2008. 

Encarnacao, J., Visser, P., Jaeggi, A., Bezdek, A., Mayer-Gürr, T., Shum, C., Arnold, D., Doornbos, E., Elmer, M., Guo, J., van den IJssel, J., Iorfida, E., Klokocnik, J., Krauss, S., Mao, X., Meyer, U., Sebera, J., Zhang, C., and Zhang, Y.: Multi-approach Gravity Field Models from Swarm GPS data,, 2019. 

Flechtner, F.: GRACE AOD1B RL04 Quality Assurance, available at: (last access: 5 June 2020), 2011. 

Flechtner, F., Neumayer, K.-H., Dahle, C., Dobslaw, H., Fagiolini, E., Raimondo, J.-C., and Güntner, A.: What Can be Expected from the GRACE-FO Laser Ranging Interferometer for Earth Science Applications?, Surv. Geophys., 37, 453–470,, 2016. 

Folkner, W. M., Williams, J. G., Boggs, D. H., Park, R. S., and Kuchynka, P.: The Planetary and Lunar Ephemerides DE430 and DE431, Interplanet. Netw. Prog. Rep, 42, available at: (last access: 5 June 2020), 2014. 

Friis-Christensen, E., Lühr, H., Knudsen, D., and Haagmans, R.: Swarm – An Earth Observation Mission investigating Geospace, Adv. Space Res., 41, 210–216,, 2008. 

Fritsche, B., Ivanov, M., Kashkovsky, A., Koppenwallner, G., Kudryavtsev, A., Voskoboinikov, U., and Zhukova, G.: Radiation pressure forces on complex spacecraft, Tech. rep., European Space Agency,, 1998. 

Guo, J. Y., Shang, K., Jekeli, C., and Shum, C. K.: On the energy integral formulation of gravitational potential differences from satelliteto-satellite tracking, Celest. Mech. Dyn. Astr., 121, 415–429,, 2015. 

Guo, X. and Zhao, Q.: A New Approach to Earth's Gravity Field Modeling Using GPS-Derived Kinematic Orbits and Baselines, Remote Sens., 11, 1728,, 2019. 

Haagmans, R.: Swarm – The Earth's Magnetic Field and Environment Explorers, vol. 1279, ESA Publications Division, Noordwijk, The Netherlands, sp-1279, available at: (last access: 5 June 2020), 2004. 

Jäggi, A., Hugentobler, U., and Beutler, G.: Pseudo-Stochastic Orbit Modeling Techniques for Low-Earth Orbiters, J. Geodesy, 80, 47–60,, 2006. 

Jäggi, A., Hugentobler, U., Bock, H., and Beutler, G.: Precise orbit determination for GRACE using undifferenced or doubly differenced GPS data, Adv. Space Res., 39, 1612–1619,, 2007. 

Jäggi, A., Beutler, G., Prange, L., Dach, R., and Mervart, L.: Assessment of GPS-only Observables for Gravity Field Recovery from GRACE, International Association of Geodesy Symposia, 133, 113–123,, 2009. 

Jäggi, A., Meyer, U., Beutler, G., Prange, L., Dach, R., and Mervart, L.: AIUB-GRACE03S, available at: (last access: 5 June 2020), 2011. 

Jäggi, A., Montenbruck, O., Moon, Y., Wermuth, M., König, R., Michalak, G., Bock, H., and Bodenmann, D.: Inter-agency comparison of TanDEM-X baseline solutions, Adv. Space Res., 50, 260–271,, 2012. 

Jäggi, A., Dahle, C., Arnold, D., Bock, H., Meyer, U., Beutler, G., and van den IJssel, J.: Swarm kinematic orbits and gravity fields from 18 months of GPS data, Adv. Space Res., 57, 218–233,, 2016. 

Jäggi, A., Weigelt, M., Flechtner, F., Güntner, A., Mayer-Gürr, T., Martinis, S., Bruinsma, S., Flury, J., Bourgogne, S., Steffen, H., Meyer, U., Jean, Y., Sušnik, A., Grahsl, A., Arnold, D., Cann-Guthauser, K., Dach, R., Li, Z., Chen, Q., van Dam, T., Gruber, C., Poropat, L., Gouweleeuw, B., Kvas, A., Klinger, B., Lemoine, J.-M., Biancale, R., Zwenzner, H., Bandikova, T., and Shabanloui, A.: European Gravity Service for Improved Emergency Management (EGSIEM) – from concept to implementation, Geophys. J. Int., 218, 1572–1590,, 2019. 

Jäggi, A., Meyer, U., Lasser, M., Jenny, B., Lopez, T., Flechtner, F., Dahle, C., Förste, C., Mayer-Gürr, T., Kvas, A., Lemoine, J.-M., Bourgogne, S., Weigelt, M., and Groh, A.: International Combination Service for Time-variable Gravity Fields (COST-G) – Start of operational phase and future perspectives, in: IAG Symposia Series,, in press, 2020. 

Jean, Y., Meyer, U., and Jäggi, A.: Combination of GRACE monthly gravity field solutions from different processing strategies, J. Geodesy, 92, 1313–1328,, 2018. 

Jekeli, C.: The determination of gravitational potential differences from satellite-to-satellite tracking, Celest. Mech. Dyn. Astr., 75, 85–101,, 1999. 

Kermarrec, G., Ren, L., and Schön, S.: On filtering ionospheric effects in GPS observations using the Matérn covariance family and its impact on orbit determination of Swarm satellites, GPS Solutions, 22, 66,, 2018. 

Klinger, B. and Mayer-Gürr, T.: The role of accelerometer data calibration within GRACE gravity field recovery: Results from ITSGGrace2016, Adv. Space Res., 58, 1597–1609,, 2016. 

Klinger, B., Mayer-Gürr, T., Behzadpour, S., Ellmer, M., Kvas, A., and Zehentner, N.: The new ITSG-Grace2016 release, in: EGU General Assembly, Research Gate, Vienna, Austria,, 2016. 

Knocke, P., Ries, J., and Tapley, B.: Earth radiation pressure effects on satellites, in: Astrodynamics Conference, American Institute of Aeronautics and Astronautics, Reston, Virigina,, 1988. 

Kroes, R.: Precise Relative Positioning of Formation Flying Spacecraft using GPS, PhD thesis, Delft University of Technology, available at:, 2006. 

Lieske, J. H., Lederle, T., Fricke, W., and Morando, B.: Expression for the precession quantities based upon the IAU (1976) system of astronomical constants, Astron. Astrophys., 58, 1–16, 1977. 

Loomis, B. D., Rachlin, K. E., and Luthcke, S. B.: Improved Earth Oblateness Rate Reveals Increased Ice Sheet Losses and Mass-Driven Sea Level Rise, Geophys. Res. Lett., 46, 6910–6917,, 2019. 

Lück, C., Kusche, J., Rietbroek, R., and Löcher, A.: Time-variable gravity fields and ocean mass change from 37 months of kinematic Swarm orbits, Solid Earth, 9, 323–339,, 2018. 

Lyard, F., Lefevre, F., Letellier, T., and Francis, O.: Modelling the global ocean tides: modern insights from FES2004, Ocean Dynam., 56, 394–415,, 2006. 

Mao, X., Visser, P. N., and van den IJssel, J.: Impact of GPS antenna phase center and code residual variation maps on orbit and baseline determination of GRACE, Adv. Space Res., 59, 2987–3002,, 2017. 

Mao, X., Visser, P., and van den IJssel, J.: The impact of GPS receiver modifications and ionospheric activity on Swarm baseline determination, Acta Astronaut., 146, 399–408,, 2018. 

Mayer-Gürr, T.: Gravitationsfeldbestimmung aus der Analyse kurzer Bahnbögen am Beispiel der Satellitenmissionen CHAMP und GRACE, Phd thesis, Rheinischen Friedrich-Wilhelms Universität Bonn, available at: access: 5 June 2020), 2006. 

Mayer-Gürr, T.: The Combined Satellite Gravity Field Model GOCO05s, in: EGU General Assembly, EGU2015-12364, Vienna, Austria, 2015. 

Mayer-Gürr, T., Kurtenbach, E., Eicker, A., Mayer-Gürr, T., Kurtenbach, E., and Eicker, A.: ITG-Grace2010: the new GRACE gravity field release computed in Bonn, in: EGU General Assembly, EGU2010-2446, Vienna, Austria, available at: (last access: 5 June 2020), 2010. 

Mayer-Gürr, T., Behzadpour, S., Ellmer, M., Kvas, A., Klinger, B., and Zehentner, N.: ITSG-Grace2016 – Monthly and Daily Gravity Field Solutions from GRACE,, 2016. 

Meyer, U., Jean, Y., Kvas, A., Dahle, C., Lemoine, J. M., and Jäggi, A.: Combination of GRACE monthly gravity fields on the normal equation level, J. Geodesy, 93, 1645–1658,, 2019a. 

Meyer, U., Sosnica, K., Arnold, D., Dahle, C., Thaller, D., Dach, R., Jäggi, A., Meyer, U., Sosnica, K., Arnold, D., Dahle, C., Thaller, D., Dach, R., and Jäggi, A.: SLR, GRACE and Swarm Gravity Field Determination and Combination, Remote Sens., 11, 956,, 2019b. 

Meyer, U.: Combination of monthly Swarm gravity fields applying variance component estimation, in preparation, 2020. 

Montenbruck, O. and Gill, E.: Satellite Orbits, Springer-Verlag Berlin And Heidelberg Gmbh, Berlin, Heidelberg, 1st edn.,, 2000. 

Olsen, N., Friis-Christensen, E., Floberghagen, R., Alken, P., Beggan, C. D., Chulliat, A., Doornbos, E., da Encarnação, J. T., Hamilton, B., Hulot, G., van den IJssel, J., Kuvshinov, A., Lesur, V., Lühr, H., Macmillan, S., Maus, S., Noja, M., Olsen, P. E. H., Park, J., Plank, G., Püthe, C., Rauberg, J., Ritter, P., Rother, M., Sabaka, T. J., Schachtschneider, R., Sirol, O., Stolle, C., Thébault, E., Thomson, A. W. P., Tøffner-Clausen, L., Velímský, J., Vigneron, P., and Visser, P. N.: The Swarm Satellite Constellation Application and Research Facility (SCARF) and Swarm data products, Earth, Planets Space, 65, 1189–1200,, 2013. 

Petit, G. G. and Luzum, B.: IERS Conventions (2010), available at: (last access: 5 June 2020), 2010. 

Picone, J. M., Hedin, A. E., Drob, D. P., and Aikin, A. C.: NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues, J. Geophys. Res.-Space, 107, 1468,, 2002. 

Ray, R. D. and Luthcke, S. B.: Tide model errors and GRACE gravimetry: towards a more realistic assessment, Geophys. J. Int., 167, 1055–1059,, 2006. 

Reigber, C.: Gravity field recovery from satellite tracking data, in: Theory of Satellite Geodesy and Gravity Field Determination, edited by: Sansò, F. and Rummel, R., vol. 25 in: Lecture Notes in Earth Sciences, Springer, Berlin, Heidelberg, 197–234,, 1989. 

Ries, J., Bettadpur, S., Eanes, R., Kang, Z., Ko, U., McCullough, C., Nagel, P., Pie, N., Poole, S., Richter, T., Save, H., and Tapley, B.: The Combined Gravity Model GGM05C, Tech. Rep. CSR-TM-16-01, Center for Space Research, University of Texas at Austin, Austin,, 2016. 

Ries, J. C., Bettadpur, S., Poole, S., and Richter, T.: Mean Background Gravity Fields for GRACE processing, GRACE Science Team Meeting, Austin, USA, 8–10 August 2010, available at: (last access: 5 June 2020), 2011. 

Rodriguez-Solano, C. J., Hugentobler, U., Steigenberger, P., and Lutz, S.: Impact of Earth radiation pressure on GPS position estimates, J. Geodesy, 86, 309–317,, 2012. 

Savcenko, R. and Bosch, W.: EOT11a – Empirical ocean tide model from multi-mission satellite altimetry, Tech. rep., Deutsches Geodätisches Forschungsinstitut, München, Germany, available at: (last access: 5 June 2020), 2012. 

Schreiter, L., Arnold, D., Sterken, V., and Jäggi, A.: Mitigation of ionospheric signatures in Swarm GPS gravity field estimation using weighting strategies, Ann. Geophys., 37, 111–127,, 2019. 

Seidelmann, P. K.: 1980 IAU Theory of Nutation: The final report of the IAU Working Group on Nutation, Celestial Mech., 27, 79–106,, 1982. 

Sentman, L. H.: Free molecule flow theory and its application to the determination of aerodynamic forces, LMSC-448514, available at: (last access: 5 June 2020), 1961. 

Shang, K., Guo, J., Shum, C., Dai, C., and Luo, J.: GRACE time-variable gravity field recovery using an improved energy balance approach, Geophys. J. Int., 203, 1773–1786,, 2015. 

Siemes, C.: Swarm satellite thermo-optical properties and external geometry, Tech. rep., European Space Agency, available at: (last access: 5 June 2020), 2019. 

Siemes, C., Teixeira da Encarnação, J., Doornbos, E., van den IJssel, J., Kraus, J., Pereštý, R., Grunwaldt, L., Apelbaum, G., Flury, J., and Holmdahl Olsen, P. E.: Swarm accelerometer data processing from raw accelerations to thermospheric neutral densities, Earth Planets Space, 68, 92,, 2016. 

Tapley, B. D., Bettadpur, S., Watkins, M., and Reigber, C.: The gravity recovery and climate experiment: Mission overview and early results, Geophys. Res. Lett., 31, L09607,, 2004. 

Teixeira da Encarnação, J. and Visser, P.: TN-03: Swarm models validation, Tech. rep., TU Delft,, 2019. 

Teixeira da Encarnação, J., Arnold, D., Bezdek, A., Dahle, C., Doornbos, E., van den IJssel, J., Jäggi, A., Mayer-Gürr, T., Sebera, J., Visser, P., and Zehentner, N.: Gravity field models derived from Swarm GPS data, Earth Planets Space, 68, 127,, 2016. 

Teunissen, P. J. G.: The least-squares ambiguity decorrelation adjustment: a method for fast GPS integer ambiguity estimation, J. Geodesy, 70, 65–82,, 1995. 

van Barneveld, P. W. L.: Orbit determination of satellite formations, Phd thesis, Delft University of Technology,, 2012. 

van den IJssel, J., Encarnação, J., Doornbos, E., and Visser, P.: Precise science orbits for the Swarm satellite constellation, Adv. Space Res., 56, 1042–1055,, 2015. 

van den IJssel, J., Forte, B., and Montenbruck, O.: Impact of Swarm GPS receiver updates on POD performance, Earth Planets Space, 68, 85,, 2016. 

van Helleputte, T.: GPS High Precision Orbit Determination Software Tools: User Manual (No. FDS-SUM-3110)., Deutsches Zentrum für Luft- und Raumfahrt, Oberpfaffenhofen, available at: (last access: 5 June 2020), 2004. 

Visser, P. N. A. M., Sneeuw, N., and Gerlach, C.: Energy integral method for gravity field determination from satellite orbit coordinates, J. Geodesy, 77, 207–216,, 2003. 

Wahr, J., Nerem, R. S., and Bettadpur, S. V.: The pole tide and its effect on GRACE time-variable gravity measurements: Implications for estimates of surface mass variations, J. Geophys. Res.-Sol. Ea., 120, 4597–4615,, 2015. 

Wermuth, M., Montenbruck, O., and Helleputte, T. V.: GPS high precision orbit determination software tools (GHOST), in: 4th International Conference on Astrodynamics Tools and Techniques, ESA WPP-308, Madrid, 2010. 

Zehentner, N.: Kinematic orbit positioning applying the raw observation approach to observe time variable gravity, Doctoral dissertation, Graz University of Technology,, 2016. 

Zehentner, N. and Mayer-Gürr, T.: New Approach to Estimate Time Variable Gravity Fields from High-Low Satellite Tracking Data, in: International Association of Geodesy Symposia, Venice, Italy, Springer, Cham, 141, 111–116,, 2014. 

Zehentner, N. and Mayer-Gürr, T.: Precise orbit determination based on raw GPS measurements, J. Geodesy, 90, 275–286,, 2016.  

Zeng, Y., Guo, J., Shang, K., Shum, C., and Yu, J.: On the formulation of gravitational potential difference between the GRACE satellites based on energy integral in Earth fixed frame, Geophys. J. Int., 202, 1792–1804,, 2015. 

Short summary
Although not the primary mission of the Swarm three-satellite constellation, the sensors on these satellites are accurate enough to measure the melting and accumulation of Earth’s ice reservoirs, precipitation cycles, floods, and droughts, amongst others. Swarm sees these changes well compared to the dedicated GRACE satellites at spatial scales of roughly 1500 km. Swarm confirms most GRACE observations, such as the large ice melting in Greenland and the wet and dry seasons in the Amazon.