the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Description of the multiapproach 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 MayerGürr
Ulrich Meyer
Josef Sebera
C. K. Shum
Chaoyang Zhang
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 timevariable component with the purpose of estimating the hydrological and glaciological yearly cycles and longterm trends. Currently, the GRACEFO 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 GRACEFO, while the Swarm satellites have collected gravimetric data with their GPS receivers since December 2013.
We present highquality 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 10month gap between GRACE and GRACEFO, as well as the short gaps in the existing GRACE and GRACEFO 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 Timevariable Gravity Fields (COSTG) 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 highestquality 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 nongravitational 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 nongravitational 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 mid2015 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 largescale 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 GRACEFO 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 https://swarmdiss.eo.esa.int/, last access: 5 June 2020, follow Level2longterm and then EGF) and at the International Centre for Global Earth Models (http://icgem.gfzpotsdam.de/series/02_COSTG/Swarm, last access: 5 June 2020), as well as identified with the DOI https://doi.org/10.5880/ICGEM.2019.006 (Encarnacao et al., 2019).
Swarm is the fifth Earth Explorer mission by the European Space Agency (ESA), launched on 22 November 2013 (Haagmans, 2004; FriisChristensen 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 nearpolar orbit. In addition to a sophisticated instrument suite for observing the geomagnetic and electric fields, the Swarm satellites are equipped with highprecision, dualfrequency 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 longwavelength temporal variations with highquality GPS receivers on board lowEarth orbit (LEO) satellites (Zehentner and MayerGürr, 2014; Bezdek et al., 2016; Dahle et al., 2017). For Swarm, Teixeira da Encarnação et al. (2016) successfully demonstrated the observation of longwavelength 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 GPSbased gravity field solutions is to guarantee longterm observation of mass transport in the Earth system. The geophysical community has identified the need for continued monitoring of timevariable gravity for estimating the hydrological and glaciological yearly cycles and longterm trends (Abdalati et al., 2018). The US–German GRACE mission (Tapley et al., 2004) was by far the most important spaceborne global provider of the needed data for the period from April 2002 until July 2017. GRACE FollowOn (GRACEFO) was launched in May 2018 and is expected to continue the highquality observation of Earth's timevariable gravity field for at least 5 years (Flechtner et al., 2016). Thus a time gap exists between the GRACE and GRACEFO missions, and, importantly, no missions have yet been selected for the postGRACEFO period. It can thus be claimed that the only guarantee for sustained observation of timevariable gravity from space is constituted by spaceborne GPS receivers on LEO satellites. Moreover, the associated data can be used to fill the gap between the GRACE and GRACEFO 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 sealevel rise and (ii) the opportunity to provide independent estimates of temporal variations in lowdegree coefficients, in particular related to C_{2,0} and C_{3,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 Swarmbased observation of longwavelength timevariable gravity in support of the operational delivery of Swarmbased 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 shortarc approach (SAA) (MayerGü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 satellitetosatellite tracking (hlSST) 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 formationflying spacecraft can typically be derived with much better precision than the absolute positions by making use of ambiguityfixing schemes and due to cancellation of common errors (Kroes, 2006; AllendeAlba 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 nongravitational 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 GRACEFO 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 Timevariable Gravity Fields (COSTG), 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 https://swarmdiss.eo.esa.int/, last access: 5 June 2020, follow Level2longterm and then EGF) and at the International Centre for Global Earth Models (http://icgem.gfzpotsdam.de/series/02_COSTG/Swarm, last access: 5 June 2020), as well as identified with the DOI https://doi.org/10.5880/ICGEM.2019.006 (Encarnacao et al., 2019).
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 KOdriven 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.
^{1} ftp://ftp.aiub.unibe.ch/LEO_ORBITS/SWARM (last access: 5 June 2020). ^{2} ftp://ftp.tugraz.at/outgoing/ITSG/tvgogo/orbits/Swarm (last access: 5 June 2020). ^{3} http://earth.esa.int/web/guest/swarm/dataaccess (last access: 5 June 2020).
2.2.1 KBs produced at AIUB
Kinematic and reduceddynamic 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 reduceddynamic solution generated from zerodifferenced (ZD) ionospherefree GPS carrier phase observations. Reduceddynamic orbit parameters of the other satellite (Swarm C) are estimated by processing doubledifferenced (DD) ionospherefree GPS carrier phase observations with DD ambiguities resolved to their integer values. First, the Melbourne–Wübbena linear combination is analysed to resolve the widelane ambiguities, which are subsequently introduced as known to resolve the narrowlane ambiguities together with the reduceddynamic 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 reduceddynamic 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 reduceddynamic 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 leastsquares 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 bidirectional 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 reduceddynamic baseline, a KB solution is produced using the leastsquares (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 reduceddynamic 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 2sigma 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 bidirectionally to compute two solutions that are averaged according to the epochwise 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 ITSGGrace2016 (Klinger et al., 2016).
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 nongravitational 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 ITSGGrace releases (MayerGü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 KBand ranging instrument. Our software is not prepared to handle the full threedimensional (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) hlSST solution from TUD KOs, (2) hlSST + low–low satellitetosatellite tracking (llSST) solution from TUD KOs and KBs, (3) hlSST solution from AIUB KOs, and (4) hlSST + llSST solution from AIUB KOs and KBs. The four solution variants were produced for all seven test months.
2.3 Nongravitational accelerations
We assessed the quality of the Swarm GFMs when the nongravitational 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 nongravitational 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 nongravitational 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 inhouse orbital propagator NUMINTSAT (Bezdek et al., 2009) for processing the satellite orbital data, computing the coordinate transformations and generating the modelled nongravitational accelerations of each Swarm satellite. The computation of the nongravitational 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 longterm 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).
For the TUD model, the Near RealTime 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 nongravitational 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 singlesided 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 corotation, 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 Sunsatellite distance.

The Sun–Earth eclipse model takes into account atmospheric absorption and refraction, according to the analysis of nongravitational 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 alongtrack 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 crosstrack 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 nongravitational acceleration measurements are restricted to the available alongtrack Swarm C data.
2.4 Gravity field model estimation approaches
The estimation of the hlSST 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 socalled 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 doubledifferentiated 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 highfrequency noise. However, it is possible to mitigate the highfrequency 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 llSST data by Shang et al. (2015), the term related to the Earth's rotation cannot be neglected in the processing of hlSST 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 reduceddynamic positions (integrated on the basis of the reference background force model), and the observation covariance matrix constructed from the epochwise variance–covariance information distributed in the KO data files. As a consequence of this orbit smoothing procedure, we discard the warmup–cooldown edges of the daily data arcs. We further remove one cycle per revolution (CPR) sinusoidal and 3hourly 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 Shortarc approach
The shortarc approach (MayerGü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 highfrequency noise.
2.5 Combination
The individual GFMs are combined in the frame of the Combination Service of Timevariable Gravity Fields (COSTG) 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 timedependent 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 GRACEFO 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 geocentre 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 timevariable 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 coefficientwise time series.
We consider the entirety of the Swarm GFM time series, irrespective of the epochwise 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 GRACEtoGRACEFO gap, where no interpolation is performed. We detrend the time series of models at the level of the Stokes coefficients when computing nonlinear statistics, notably the epochwise 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 C_{2,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 C_{2,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 C_{2,0} estimates over large gaps would lead to unrealistic mass variations.
For this reason, we selected the C_{2,0} weekly time series from Loomis et al. (2019), since the necessary interpolation introduces negligible deviations. We are not advocating that the considered C_{2,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 Deepocean 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 cutoff value of 0.9 to be in deepocean areas. The cutoff 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 GRACEFO data in this regression, in order to be able to verify the continuity of the GRACEFO 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.
We selected the firstorder polynomial to represent bias and trend in the GRACE data. For the periodic parameters, we choose the year and semiyear 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 cosine components. Each set of parametric Stokes coefficients has an implicit time dependence which is evaluated coefficientwise at the epochs of the Swarm GFMs. We illustrate the general agreement between the climatological model and the GRACE data for the case of C_{3,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).
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 onboard 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 ambiguityfixed and kinematic ambiguityfixed baselines may be as low as 10–15, 4–6, and 3–5 mm on average for the radial, alongtrack, and crosstrack 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 lowquality 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 3 shows the degree amplitudes with respect to the static part of the GOCO release 05 satelliteonly gravity field model (GOCO05S) (MayerGü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 GRACEonly model 2016 (ITSGGrace2016) 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 hlSSTonly and the llSST + hlSST 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 llSST and hlSST 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 longwavelength part of the solutions, we have compared the individual solutions to ITSGGrace2016 monthly solutions in the spatial domain. The solutions are evaluated on an equiangular grid ($\mathrm{1}{}^{\circ}\times \mathrm{1}{}^{\circ}$), reduced by the corresponding ITSGGrace2016 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 highfrequency 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. KOonly solutions are already of very similar quality when compared to KBaugmented 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 hlSST 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 lineofsight (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 nongravitational 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 hlSSTonly 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 hlSSTonly GFMs.
3.2 Nongravitational accelerations
In this section, we present the intercomparison of the three types of nongravitational accelerations described in Sect. 2.3. Figure 4 compares three singlesatellite gravity field solutions derived from Swarm C data, considering the three nongravitational accelerations, for January 2015.
The SH degree difference amplitudes illustrate that the measured nongravitational accelerations improve the agreement of the lowest degrees of the Swarm C monthly solution with respect to the GOCO05S model (MayerGürr, 2015), which includes a timevariable component. We tested this comparison relative to the ITSGGrace2016 monthly GFM (MayerGürr et al., 2016) and observed similar results (not shown). The improvement at the lowest degrees in the Swarm C model when using observed nongravitational acceleration data is in accordance with what was reported by Klinger and MayerGürr (2016), relative to GRACE gravity field recovery.
In view of the lack of reliable measured nongravitational accelerations in Swarm A and Swarm B, the threesatellite Swarm GFM considers the ASU modelled nongravitational accelerations for these satellites. For Swarm C, we consider three cases where the nongravitational accelerations are either measured or represented by TUD or ASU's model. In this way, we isolate the effect of the three types of nongravitational acceleration data. The results for January 2015 are shown in Fig. 5, using ASU and TUD models and calibrated accelerometer data.
The threesatellite solutions that use modelled nongravitational 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 nongravitational accelerations, we derive the following model difference D, similar to rms:
with the median absolute deviation (MAD) analogous to SD when the median is considered instead of the mean and Δh being the $\mathrm{1}{}^{\circ}\times \mathrm{1}{}^{\circ}$ geoid height difference between the 500 km Gaussian filtering threesatellite Swarm models and both ITSGGrace2016 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.
The 2015 results indicate that the observed nongravitational accelerations improve the agreement between the threesatellite Swarm models and ITSGGrace2016 and GOCO05S, while that is not the case for 2016 (except for January 2016 and GOCO05S, when the GFM derived from TUDmodelled nongravitational accelerations agrees equally well with the one derived considering observed nongravitational accelerations). The comparison with GOCO05S intends to predict how possible it would be to assess the added value of the different types of nongravitational accelerations during those periods when there are no GRACE data. Other timedependent models were tested, but those do not agree as closely with GRACE monthly models (not shown).
The statistics in Table 6 imply that observed nongravitational accelerations are only beneficial when the amplitude of the nongravitational 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 11year 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 signaltonoise 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 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 EBAtype 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.
Another way of analysing the agreement between the individual solutions and GRACE is to derive percoefficient statistics of their temporal variations. One such statistic is the coefficientwise 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 perdegree 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 perdegree 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.
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.
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.
In the case of 1500 km smoothing, the Swarm differences have comparable amplitudes to the GRACE data since mid2015. 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 deepocean signal
This section illustrates the differences in SNR characteristic of the Swarm GFMs by computing separate statistics for land and deepocean areas, the latter defined in Sect. 2.6.2.
In Fig. 11 the rms of the deepocean 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) longterm trends since the data were detrended 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).
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 C_{20} published in Cheng and Ries (2018).
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 mid2015, 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.
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.
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 nontropical regions, as illustrated in Fig. 15.
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 nontropical regions. We note that the deepocean 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 nontropical regions, which are of comparable amplitude between themselves and complementary.
We now focus on the necessary smoothing to retrieve any deepocean 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 largescale 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 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 largescale 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 perdegree 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 largescale Swarm deepocean mass changes since mid2015.
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 lowdegree 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 GRACEFO.
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 perdegree correlation coefficients of Swarm and GRACE relative to the climatological model. The degree 2 coefficients (except C_{2,0}), which are particularly important for sealevel studies, are subsequently presented. Finally, we show the selected case of ${C}_{\mathrm{5},\mathrm{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. ${C}_{\mathrm{2},\mathrm{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 C_{5,5} seem to be poorly captured by the climatological mode, since neither Swarm nor GRACE correlated well.
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 mid2015, 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 C_{2,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).
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 monthtomonth changes (i.e. larger errors) and with frequent overshootings 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 intersatellite 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 ${C}_{\mathrm{5},\mathrm{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.
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–westbanded 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 largescale 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.
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 C_{2,0} replacement. The trend (and bias) is coestimated with yearly and semiyearly 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 mid2015. Swarm also agrees with the GRACEFO observation that the Greenland ice mass loss slowed down during the winter of 2018–2019, since both Swarm and GRACEFO 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 GRACEFO months agree particularly well with Swarm.
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 mid2015 is ignored, these statistics improve substantially (not shown).
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 mid2015 is discarded.
The Swarm monthly models are distributed on a quarterly basis at ESA's Earth Swarm Data Access (at https://swarmdiss.eo.esa.int/, last access: 5 June 2020, follow “Level2longterm” and then “EGF”) and at the International Centre for Global Earth Models (http://icgem.gfzpotsdam.de/series/02_COSTG/Swarm, last access: 5 June 2020), as well as identified with the DOI https://doi.org/10.5880/ICGEM.2019.006 (Encarnacao et al., 2019).
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 shortarc 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 nongravitational accelerations. We compared the threesatellite Swarm solution produced considering the DAA and nongravitational accelerations acting on Swarm C represented by the TUD and ASU nongravitational 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 nongravitational 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 nongravitational accelerations in the processing of GPSdriven GFMs, we would like to comment on the results of Ditmar et al. (2006) and Ditmar et al. (2008), who demonstrated that nongravitational 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 piecewiseconstant empirical parametrization does not require any modelling of nongravitational accelerations. In the case of ASU, who exploit a dedicated decorrelation procedure (which is a frequencydependent noise whitening procedure), their solutions benefit from drag, EARP, and EIRP models. Essentially, the inclusion of frequencydependent 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 mid2015 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 monthtomonth variations. We attribute this particularity to the lower accuracy of the GPS observations as gravimetric data, compared to GRACE's Kband 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 GRACEFO data, notably the atypical deviation around early 2016 in the ${C}_{\mathrm{5},\mathrm{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 GRACEFO 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 highfrequency noise in the differentiation of the KO positions into velocities.
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:  Ionospherefree 
GPS observations:  Code and carrier phase 
Estimator:  Bayesianweighted 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 reduceddynamic 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  
reduceddynamic 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:  Ionospherefree 
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 reduceddynamic precise orbit determination (POD) of 
approx. 120 d, nine iterations, 1^{∘} binning  
Data screening:  2 cm s^{−1} or larger time differences of the geometryfree 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 coestimated) 
GPS observations:  Code and carrier phase 
Estimator:  LS 
Arc length:  24 h 
Data weighting:  Elevation and azimuthdependent, epochwise VCE 
Transmitter PCV:  Empirical, estimated from 5.5 years of data, including data from several LEO missions (GRACE, 
Jason 2 & 3, MetOpA & MetOpB, Sentinel 3A, Swarm, TanDEMX, TerraSARX) (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:  Coestimated 
Sampling rate:  10 or 1 s (depending on L1B GPS data) 
Elevation cutoff 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 
B1 Delft University of Technology
Software:  Multiplesatellite orbit determination using Kalman filtering (van Barneveld, 2012) 
Linear combination:  n/a (the ionospheric frequencydependent influence is modelled) 
Estimator:  Iterative EKF 
Carrier phase ambiguities:  Integer, using the leastsquares ambiguity decorrelation 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:  Ionospherefree 
Estimator:  LS 
Carrier phase ambiguities:  Widelane and narrowlane integer ambiguity fixing with the Melbourne–Wübbena combination 
and the ionospherefree linear combination, respectively  
Receiver PCV:  Empirical 
B3 Common
Differencing scheme:  Doubledifferenced 
GPS observations:  Code and carrier phase 
Carrier phase ambiguities:  Integer 
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 GRACEonly static model, version 3 (Jäggi et al., 2011) 
Empirical parameters:  Daily and 15 min piecewiseconstant (constrained) 
Drag model:  None 
EARP and EIRP models:  None 
Nontidal model:  Atmosphere and Ocean Dealiasing Level 1B (Flechtner, 2011) 
Ocean tidal model:  2011 Empirical Ocean Tide model (Savcenko and Bosch, 2012) 
Permanent tide system:  Tidefree 
C2 Astronomical Institute Ondřejov
Software:  Developed inhouse 
Approach:  Decorrelated acceleration approach (Bezdek et al., 2014) 
Reference GFM:  ITG GRACEonly static model, 2010 (MayerGürr et al., 2010) 
Empirical parameters:  Daily constantpiecewise 
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) 
Nontidal model:  Atmosphere and Ocean Dealiasing Level 1B (Dobslaw et al., 2017) 
Ocean tidal model:  2004 finite element solution (Lyard et al., 2006) 
Permanent tide system:  Tidefree 
C3 Institute of Geodesy Graz
Software:  GROOPS 
Approach:  Shortarc approach (MayerGürr, 2006) 
Reference GFM:  GOCO release 05 satelliteonly gravity field model (MayerGü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:  RodriguezSolano et al. (2012) 
Nontidal model:  Atmosphere and Ocean Dealiasing 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 inhouse 
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:  Secondorder polynomial every 3 h, 1CPR 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) 
Nontidal model:  Atmosphere and Ocean Dealiasing Level 1B (Flechtner, 2011) 
Ocean tidal model:  2011 empirical ocean tide model (Savcenko and Bosch, 2012) 
Permanent tide system:  Tidefree 
C5 Common
Atmospheric tidal model:  Biancale and Bode (2006) 
Solid Earth tidal model:  IERS2010 
Pole tidal model:  IERS2010 
Ocean pole tidal model:  IERS2010 
Thirdbody perturbations:  Sun, Moon, Mercury, Venus, Mars, Jupiter, and Saturn, following the JPL Planetary and 
Lunar Ephemerides (Folkner et al., 2014)  
C_{2,0} coefficient:  Estimated alongside other coefficients 
Figures D1–D3 illustrate the time series for the zonal coefficients of degrees 3 to 6, respectively.
The zonal coefficient of degree 3 is an interesting case because both Swarm and GRACEFO observe a phase shift during late 2018, relative to the climatological model, which is well in phase with GRACE for the nonGRACEFO period (2003–2017). Swarm already captures this phase shift possibly as early as mid2017, 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 C_{2,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 mid2015. These are features intrinsic to our Swarm solutions.
3D  Threedimensional 
AA  Acceleration approach 
DAA  Decorrelated acceleration approach 
AC  Analysis centre 
AIUB  Astronomical Institute of the University of Bern, Switzerland 
AIUBGRACE03S  AIUB GRACEonly static model, version 3 
AOD1B  Atmosphere and Ocean Dealiasing Level 1B 
AOD1BRL06  Atmosphere and Ocean Dealiasing Level 1B RL06 
ANGARA  Analysis of nongravitational 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 MiniSatellite Payload 
CODE  Centre for Orbit Determination in Europe 
CMA  Celestial mechanics approach 
CoM  Centre of mass 
COSTG  Combination Service of Timevariable 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  Doubledifferenced 
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  Frequencydependent 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 SteadyState Ocean Circulation Explorer 
GOCO  Gravity Observation Combination 
GOCO05S  GOCO release 05 satelliteonly gravity field model 
GPS  Global Positioning System 
GRACE  Gravity Recovery and Climate Experiment 
GRACEFO  GRACE FollowOn 
GROOPS  Gravity Recovery Object Oriented Programming System 
hlSST  Highlow satellitetosatellite 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  Intersatellite range 
ITG  Institut für Geodäsie und Geoinformation, Germany 
ITSG  Institute of Theoretical Geodesy and Satellite Geodesy 
ITGGrace2010s  ITG GRACEonly static model, 2010 
ITSGGrace2016  ITSG GRACEonly model, 2016 
JB2008  JacchiaBowman 2008 
JPL  Jet Propulsion Laboratory, USA 
JPLPLE  JPL Planetary and Lunar Ephemerides 
KB  Kinematic baseline 
KBR  Kband ranging 
KO  Kinematic orbit 
L1A  Level 1A data 
L1B  Level 1B data 
L2  Level 2 data 
L2PS  Level 2 processing system 
LAMBDA  Leastsquares ambiguity decorrelation adjustment 
LEO  LowEarth orbit 
llSST  Low–low satellitetosatellite tracking 
LoS  Line of sight 
LS  Least squares 
MAD  Median absolute deviation 
MODK  Multiplesatellite 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 RealTime Density Model 
OSU  Ohio State University 
PCV  Phase centre variation 
POD  Precise orbit determination 
PSO  Precise or postprocessed science orbit 
RINEX  Receiver Independent Exchange Format 
RL05  Release 5 
RL06  Release 6 
Rms  Root mean square 
SAA  Shortarc approach 
SC  Stokes coefficient 
SH  Spherical harmonic 
SLR  Satellite laser ranging 
SNR  Signaltonoise 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  Zerodifferenced 
AB, ED, JIJ, EI, and JS produced the modelled nongravitational data. ED, JIJ, and EI preprocessed 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 coauthors reviewed the manuscript.
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 SLRderived C_{20} time series. The authors acknowledge the reviewers for their valuable remarks that helped improve the quality of the original manuscript.
This research was funded by the European Space Agency (contracts SWCODTUGS111 and SWCNDTUGS027, part of contract 4000109587/13/INB), 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 2017YFA06031033), 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).
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, https://doi.org/10.17226/24938, 2018.
AllendeAlba, G., Montenbruck, O., Jäggi, A., Arnold, D., and Zangerl, F.: Reduceddynamic and kinematic baseline determination for the Swarm mission, GPS Solutions, 21, 1275–1284, https://doi.org/10.1007/s102910170611z, 2017.
Barkstrom, B. R. and Smith, G. L.: The Earth Radiation Budget Experiment: Science and implementation, Rev. Geophys., 24, 379–390, https://doi.org/10.1029/RG024i002p00379, 1986.
Bettadpur, S.: UTCSR Level2 Processing Standards Document For Level2 Product Release 0006, Tech. rep., Center for Space Research, Austin, USA, available at: https://podaactools.jpl.nasa.gov/drive/files/allData/grace/docs/TN11_C20_SLR.txt (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, https://doi.org/10.1007/s0019001004017, 2010.
Bezdek, A.: Calibration of accelerometers aboard GRACE satellites by comparison with PODbased nongravitational acceler ations, J. Geodynam., 50, 410–423, https://doi.org/10.1016/j.jog.2010.05.001, 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, https://doi.org/10.1016/j.jog.2009.01.007, 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, https://doi.org/10.1016/j.asr.2013.11.031, 2014.
Bezdek, A., Sebera, J., Teixeira da Encarnação, J., and Klokocník, J.: Timevariable gravity fields derived from GPS tracking of Swarm, Geophys. J. Int., 205, 1665–1669, https://doi.org/10.1093/gji/ggw094, 2016.
Bezdek, A., Sebera, J., and Klokocník, J.: Validation of Swarm accelerometer data by modelled nongravitational forces, Adv. Space Res., 59, 2512–2521, https://doi.org/10.1016/j.asr.2017.02.037, 2017.
Bezdek, A., Arnold, D., Doornbos, E., Jäggi, A., Mao, X., Zehentner, N., Teixeira da Encarnacao, J., and Visser, P. N. A. M.: TN02: Swarm Data PreProcessing, Kinematic Baselines And Accelerometer Data, Tech. rep., https://doi.org/10.13140/RG.2.2.19891.99361, 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, https://doi.org/10.1016/j.asr.2018.04.041, 2018b.
Biancale, R. and Bode, A.: Mean annual and seasonal atmospheric tide models based on 3hourly and 6hourly ECMWF surface pressure data, Tech. rep., Deutsches GeoForschungsZentrum GFZ, Potsdam, Germany, https://doi.org/10.2312/GFZ.b10306011, 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, https://doi.org/10.2514/6.20086438, 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, https://doi.org/10.1007/s0019001609955, 2017.
Cheng, M. and Ries, J.: GRACE Technical Note 11: Monthly estimates of C20 from 5 satellites based on GRACE RL06 models, available at: ftp://podaacftp.jpl.nasa.gov/allData/grace/docs/TN11_C20_SLR.txt (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, https://doi.org/10.1029/2010JB000850, 2011.
Dach, R., Lutz, S., Walser, P., and Fridez, P.: Bernese GNSS Software Version 5.2, Bern Open Publishing, Bern, https://doi.org/10.7892/boris.72297, 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, https://doi.org/10.7892/boris.75876.2, 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, https://doi.org/10.1016/j.asr.2017.03.003, 2017.
Ditmar, P., Klees, R., and Liu, X.: Frequencydependent data weighting in global gravity field modeling from satellite data contaminated by nonstationary noise, J. Geodesy, 81, 81–96, https://doi.org/10.1007/s0019000600744, 2006.
Ditmar, P., Bezdek, A., Liu, X., and Zhao, Q.: On a Feasibility of Modeling Temporal Gravity Field Variations from Orbits of Nondedicated 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, https://doi.org/10.1007/9783540854265_36, 2008.
Dobslaw, H., BergmannWolf, I., Dill, R., Poropat, L., Thomas, M., Dahle, C., Esselborn, S., König, R., and Flechtner, F.: A new highresolution model of nontidal atmosphere and ocean mass variability for dealiasing of satellite gravity observations: AOD1B RL06, Geophys. J. Int., 211, 263–269, https://doi.org/10.1093/gji/ggx302, 2017.
Doornbos, E.: Thermospheric Density and Wind Determination from Satellite Dynamics, Springer Theses, Springer Berlin, Heidelberg, https://doi.org/10.1007/9783642251290, 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 SA31D2371, 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, https://doi.org/10.1029/2008JA013668, 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 storminduced disturbance winds, J. Geophys. Res.Space, 113, A11319, https://doi.org/10.1029/2008JA013541, 2008.
Encarnacao, J., Visser, P., Jaeggi, A., Bezdek, A., MayerGü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.: Multiapproach Gravity Field Models from Swarm GPS data, https://doi.org/10.5880/ICGEM.2019.006, 2019.
Flechtner, F.: GRACE AOD1B RL04 Quality Assurance, available at: http://op.gfzpotsdam.de/grace/results/grav/g007_aod1b_rl04.html (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 GRACEFO Laser Ranging Interferometer for Earth Science Applications?, Surv. Geophys., 37, 453–470, https://doi.org/10.1007/s107120159338y, 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: https://ipnpr.jpl.nasa.gov/progress_report/42196/196C.pdf (last access: 5 June 2020), 2014.
FriisChristensen, E., Lühr, H., Knudsen, D., and Haagmans, R.: Swarm – An Earth Observation Mission investigating Geospace, Adv. Space Res., 41, 210–216, https://doi.org/10.1016/j.asr.2006.10.008, 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, https://doi.org/10.5880/ICGEM.2019.006, 1998.
Guo, J. Y., Shang, K., Jekeli, C., and Shum, C. K.: On the energy integral formulation of gravitational potential differences from satellitetosatellite tracking, Celest. Mech. Dyn. Astr., 121, 415–429, https://doi.org/10.1007/s105690159610y, 2015.
Guo, X. and Zhao, Q.: A New Approach to Earth's Gravity Field Modeling Using GPSDerived Kinematic Orbits and Baselines, Remote Sens., 11, 1728, https://doi.org/10.3390/rs11141728, 2019.
Haagmans, R.: Swarm – The Earth's Magnetic Field and Environment Explorers, vol. 1279, ESA Publications Division, Noordwijk, The Netherlands, sp1279, available at: http://esamultimedia.esa.int/docs/SP_1279_6_Swarm.pdf (last access: 5 June 2020), 2004.
Jäggi, A., Hugentobler, U., and Beutler, G.: PseudoStochastic Orbit Modeling Techniques for LowEarth Orbiters, J. Geodesy, 80, 47–60, https://doi.org/10.1007/s0019000600299, 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, https://doi.org/10.1016/j.asr.2007.03.012, 2007.
Jäggi, A., Beutler, G., Prange, L., Dach, R., and Mervart, L.: Assessment of GPSonly Observables for Gravity Field Recovery from GRACE, International Association of Geodesy Symposia, 133, 113–123, https://doi.org/10.1007/9783540854265_14, 2009.
Jäggi, A., Meyer, U., Beutler, G., Prange, L., Dach, R., and Mervart, L.: AIUBGRACE03S, available at: http://icgem.gfzpotsdam.de/ (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.: Interagency comparison of TanDEMX baseline solutions, Adv. Space Res., 50, 260–271, https://doi.org/10.1016/j.asr.2012.03.027, 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, https://doi.org/10.1016/j.asr.2015.10.035, 2016.
Jäggi, A., Weigelt, M., Flechtner, F., Güntner, A., MayerGürr, T., Martinis, S., Bruinsma, S., Flury, J., Bourgogne, S., Steffen, H., Meyer, U., Jean, Y., Sušnik, A., Grahsl, A., Arnold, D., CannGuthauser, 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, https://doi.org/10.1093/gji/ggz238, 2019.
Jäggi, A., Meyer, U., Lasser, M., Jenny, B., Lopez, T., Flechtner, F., Dahle, C., Förste, C., MayerGürr, T., Kvas, A., Lemoine, J.M., Bourgogne, S., Weigelt, M., and Groh, A.: International Combination Service for Timevariable Gravity Fields (COSTG) – Start of operational phase and future perspectives, in: IAG Symposia Series, https://doi.org/10.1007/1345_2020_109, 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, https://doi.org/10.1007/s0019001811235, 2018.
Jekeli, C.: The determination of gravitational potential differences from satellitetosatellite tracking, Celest. Mech. Dyn. Astr., 75, 85–101, https://doi.org/10.1023/A:1008313405488, 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, https://doi.org/10.1007/s102910180733y, 2018.
Klinger, B. and MayerGürr, T.: The role of accelerometer data calibration within GRACE gravity field recovery: Results from ITSGGrace2016, Adv. Space Res., 58, 1597–1609, https://doi.org/10.1016/j.asr.2016.08.007, 2016.
Klinger, B., MayerGürr, T., Behzadpour, S., Ellmer, M., Kvas, A., and Zehentner, N.: The new ITSGGrace2016 release, in: EGU General Assembly, Research Gate, Vienna, Austria, https://doi.org/10.13140/RG.2.1.1856.7280, 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, https://doi.org/10.2514/6.19884292, 1988.
Kroes, R.: Precise Relative Positioning of Formation Flying Spacecraft using GPS, PhD thesis, Delft University of Technology, available at: http://resolver.tudelft.nl/uuid:1a68ee943d5544b99d8b25fa44e96922, 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 MassDriven Sea Level Rise, Geophys. Res. Lett., 46, 6910–6917, https://doi.org/10.1029/2019GL082929, 2019.
Lück, C., Kusche, J., Rietbroek, R., and Löcher, A.: Timevariable gravity fields and ocean mass change from 37 months of kinematic Swarm orbits, Solid Earth, 9, 323–339, https://doi.org/10.5194/se93232018, 2018.
Lyard, F., Lefevre, F., Letellier, T., and Francis, O.: Modelling the global ocean tides: modern insights from FES2004, Ocean Dynam., 56, 394–415, https://doi.org/10.1007/s102360060086x, 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, https://doi.org/10.1016/j.asr.2017.03.019, 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, https://doi.org/10.1016/j.actaastro.2018.03.009, 2018.
MayerGürr, T.: Gravitationsfeldbestimmung aus der Analyse kurzer Bahnbögen am Beispiel der Satellitenmissionen CHAMP und GRACE, Phd thesis, Rheinischen FriedrichWilhelms Universität Bonn, available at: https://www.researchgate.net/publication/253819808_ITGGrace2010_the_new_GRACE_gravity_field_release_computed_in_Bonn(last access: 5 June 2020), 2006.
MayerGürr, T.: The Combined Satellite Gravity Field Model GOCO05s, in: EGU General Assembly, EGU201512364, Vienna, Austria, 2015.
MayerGürr, T., Kurtenbach, E., Eicker, A., MayerGürr, T., Kurtenbach, E., and Eicker, A.: ITGGrace2010: the new GRACE gravity field release computed in Bonn, in: EGU General Assembly, EGU20102446, Vienna, Austria, available at: http://www.igg.unibonn.de/apmg/index.php?id=itggrace2010 (last access: 5 June 2020), 2010.
MayerGürr, T., Behzadpour, S., Ellmer, M., Kvas, A., Klinger, B., and Zehentner, N.: ITSGGrace2016 – Monthly and Daily Gravity Field Solutions from GRACE, https://doi.org/10.5880/icgem.2016.007, 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, https://doi.org/10.1007/s00190019012746, 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, https://doi.org/10.3390/rs11080956, 2019b.
Meyer, U.: Combination of monthly Swarm gravity fields applying variance component estimation, in preparation, 2020.
Montenbruck, O. and Gill, E.: Satellite Orbits, SpringerVerlag Berlin And Heidelberg Gmbh, Berlin, Heidelberg, 1st edn., https://doi.org/10.1007/9783642583513, 2000.
Olsen, N., FriisChristensen, 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øffnerClausen, 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, https://doi.org/10.5047/eps.2013.07.001, 2013.
Petit, G. G. and Luzum, B.: IERS Conventions (2010), available at: http://www.iers.org/TN36/ (last access: 5 June 2020), 2010.
Picone, J. M., Hedin, A. E., Drob, D. P., and Aikin, A. C.: NRLMSISE00 empirical model of the atmosphere: Statistical comparisons and scientific issues, J. Geophys. Res.Space, 107, 1468, https://doi.org/10.1029/2002JA009430, 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, https://doi.org/10.1111/j.1365246X.2006.03229.x, 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, https://doi.org/10.1007/BFb0010546, 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. CSRTM1601, Center for Space Research, University of Texas at Austin, Austin, https://doi.org/10.26153/tsw/1461, 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: http://download.csr.utexas.edu/pub/grace/Proceedings/Presentations_GSTM2011.pdf (last access: 5 June 2020), 2011.
RodriguezSolano, C. J., Hugentobler, U., Steigenberger, P., and Lutz, S.: Impact of Earth radiation pressure on GPS position estimates, J. Geodesy, 86, 309–317, https://doi.org/10.1007/s0019001105174, 2012.
Savcenko, R. and Bosch, W.: EOT11a – Empirical ocean tide model from multimission satellite altimetry, Tech. rep., Deutsches Geodätisches Forschungsinstitut, München, Germany, available at: https://epic.awi.de/36001/1/DGFI_Report_89.pdf (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, https://doi.org/10.5194/angeo371112019, 2019.
Seidelmann, P. K.: 1980 IAU Theory of Nutation: The final report of the IAU Working Group on Nutation, Celestial Mech., 27, 79–106, https://doi.org/10.1007/BF01228952, 1982.
Sentman, L. H.: Free molecule flow theory and its application to the determination of aerodynamic forces, LMSC448514, available at: http://www.dtic.mil/cgibin/GetTRDoc?Location=U2&doc=GetTRDoc.pdf&AD=AD0265409 (last access: 5 June 2020), 1961.
Shang, K., Guo, J., Shum, C., Dai, C., and Luo, J.: GRACE timevariable gravity field recovery using an improved energy balance approach, Geophys. J. Int., 203, 1773–1786, https://doi.org/10.1093/gji/ggv392, 2015.
Siemes, C.: Swarm satellite thermooptical properties and external geometry, Tech. rep., European Space Agency, available at: https://earth.esa.int/documents/10174/2563139/Swarm_thermooptical_properties_and_external_geometry.pdf (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, https://doi.org/10.1186/s4062301604745, 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, https://doi.org/10.1029/2004GL019920, 2004.
Teixeira da Encarnação, J. and Visser, P.: TN03: Swarm models validation, Tech. rep., TU Delft, https://doi.org/10.13140/RG.2.2.33313.76640, 2019.
Teixeira da Encarnação, J., Arnold, D., Bezdek, A., Dahle, C., Doornbos, E., van den IJssel, J., Jäggi, A., MayerGürr, T., Sebera, J., Visser, P., and Zehentner, N.: Gravity field models derived from Swarm GPS data, Earth Planets Space, 68, 127, https://doi.org/10.1186/s4062301604999, 2016.
Teunissen, P. J. G.: The leastsquares ambiguity decorrelation adjustment: a method for fast GPS integer ambiguity estimation, J. Geodesy, 70, 65–82, https://doi.org/10.1007/BF00863419, 1995.
van Barneveld, P. W. L.: Orbit determination of satellite formations, Phd thesis, Delft University of Technology, https://doi.org/10.4233/uuid:c5ac8599fca240ebadc6bbfeeec38fa, 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, https://doi.org/10.1016/j.asr.2015.06.002, 2015.
van den IJssel, J., Forte, B., and Montenbruck, O.: Impact of Swarm GPS receiver updates on POD performance, Earth Planets Space, 68, 85, https://doi.org/10.1186/s4062301604594, 2016.
van Helleputte, T.: GPS High Precision Orbit Determination Software Tools: User Manual (No. FDSSUM3110)., Deutsches Zentrum für Luft und Raumfahrt, Oberpfaffenhofen, available at: https://issfd.org/ISSFD_2007/73.pdf (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, https://doi.org/10.1007/s0019000303158, 2003.
Wahr, J., Nerem, R. S., and Bettadpur, S. V.: The pole tide and its effect on GRACE timevariable gravity measurements: Implications for estimates of surface mass variations, J. Geophys. Res.Sol. Ea., 120, 4597–4615, https://doi.org/10.1002/2015JB011986, 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 WPP308, Madrid, 2010.
Zehentner, N.: Kinematic orbit positioning applying the raw observation approach to observe time variable gravity, Doctoral dissertation, Graz University of Technology, https://doi.org/10.13140/RG.2.2.33916.33927, 2016.
Zehentner, N. and MayerGürr, T.: New Approach to Estimate Time Variable Gravity Fields from HighLow Satellite Tracking Data, in: International Association of Geodesy Symposia, Venice, Italy, Springer, Cham, 141, 111–116, https://doi.org/10.1007/9783319108377_14, 2014.
Zehentner, N. and MayerGürr, T.: Precise orbit determination based on raw GPS measurements, J. Geodesy, 90, 275–286, https://doi.org/10.1007/s0019001508727, 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, https://doi.org/10.1093/gji/ggv248, 2015.
 Abstract
 Introduction
 Methodology
 Results
 Data availability
 Conclusions
 Appendix A: Kinematic orbits
 Appendix B: Kinematic baselines
 Appendix C: Gravity field models
 Appendix D: Time series of zonal coefficients
 Appendix E: Storage basin time series
 Appendix F: Abbreviations
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Methodology
 Results
 Data availability
 Conclusions
 Appendix A: Kinematic orbits
 Appendix B: Kinematic baselines
 Appendix C: Gravity field models
 Appendix D: Time series of zonal coefficients
 Appendix E: Storage basin time series
 Appendix F: Abbreviations
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References