Benchmark data for verifying background model implementations in orbit and gravity field determination software

In the framework of the COmbination Service for Time-variable Gravity fields (COST-G) gravity field solutions from different analysis centres are combined to provide a consolidated solution of improved quality and robustness to the user. As in many other satellite-related sciences, the correct application of background models plays a crucial role in gravity field determination. Therefore, we publish a set of data of various commonly used forces in orbit and gravity field modelling (Earth's gravity field, tides etc.) evaluated along a one day orbit arc of GRACE, together with auxiliary data to enable easy comparisons. The benchmark data is compiled with the GROOPS software by the Institute of Geodesy (IfG) at Graz University of Technology. It is intended to be used as a reference data set and provides the opportunity to test the implementation of these models at various institutions involved in orbit and gravity field determination from satellite tracking data. In view of the COST-G GRACE and GRACE Follow-On gravity field combinations, we document the outcome of the comparison of the background force models for the Bernese GNSS software from AIUB (Astronomical Institute, University of Bern), the EPOS software of the German Research Centre for Geosciences (GFZ), the GINS software, developed and maintained by the Groupe de Recherche de Géodésie Spatiale (GRGS), the GRACE-SIGMA software of the Leibniz University of Hannover (LUH) and the GRASP software also developed at LUH. We consider differences in the force modelling for GRACE (-FO) which are one order of magnitude smaller than the accelerometer noise of about 10⁻¹⁰ m s⁻² to be negligible and formulate this as a benchmark for new analysis centres, which are interested to contribute to the COST-G initiative.

1 Introduction

The correct application of background models plays a crucial role in many satellite-related sciences, such as orbit and gravity field determination. Nowadays not only one software is used to perform such tasks but various institutions around the world have set up their own processing schemes based on their in-house developed software packages. In order to minimise systematic differences caused by a diverse handling of well established models, for example, background force models, the effort of comparing software implementations is picked up regularly. Software comparisons are always a useful tool for detecting inconsistencies or even bugs in software implementations.

In the framework of the COST-G (Jäggi et al., 2020) gravity field solutions from different analysis centres (ACs) are combined to offer a consolidated solution of improved quality, robustness and reliability to the user. As the example of several services, like the International GNSS Service (Johnston et al., 2017) and the International Laser Ranging Service (Pearlman et al., 2002), showed, combining solutions from various institutions, which are computed by different and independent software packages, leads to improved results. The COST-G initiative was formally established in 2019 and operationally provides state-of-the-art monthly global gravity models from the Gravity Recovery And Climate Experiment (GRACE, Tapley et al., 2004), GRACE Follow-On (Landerer et al., 2020) and Swarm (Friis-Christensen et al., 2006). COST-G is a product centre of the International Gravity Field Service (IGFS) under the umbrella of the International Association of Geodesy (IAG).

Modelling background forces is vital to gravity field determination, especially when the time variable part is considered. Therefore, we publish a set of data of various commonly used forces in orbit and gravity field modelling evaluated along a one day orbit arc of GRACE, together with auxiliary data to enable easy comparisons. This data set is intended to be used as a reference data set and provides the opportunity to test the implementation of these models in various software packages. The COST-G consortium consists at the time of writing of ACs and further candidate ACs. Thus, we took the opportunity to test several software packages available within the COST-G consortium, which are the Bernese GNSS software (Dach et al., 2015) from AIUB (Astronomical Institute, University of Bern), EPOS (Zhu et al., 2004) of the German Research Centre for Geosciences (GFZ), the GINS software, developed and maintained by the Groupe de Recherche de Géodésie Spatiale (GRGS), the GRACE-SIGMA (Koch et al., 2020) software of the Leibniz University of Hannover (LUH) and the GRASP software (Weigelt et al., 2013) from LUH. Furthermore, we use this reference data to set a benchmark for new analysis centres, which are interested to contribute to the COST-G initiative.

The paper is organised in five sections, where Sect. 2 introduces the data set used and published in this study. Section 3 summarises the evaluation of each force model, intended to be easily comprehensible for programming. We discuss the result of a force comparison for one software package in more detail in Sect. 4 and give the summarised outcome of the comparisons for all COST-G ACs and candidate ACs. Finally, Sect. 5 concludes the results of this study and highlights the importance beyond the COST-G group.

2 Benchmark data set

The benchmark data set was compiled at the Institute of Geodesy (IfG) at Graz University of Technology. It consists of several accelerations a spacecraft experiences, evaluated along a given orbit, which are commonly used in orbit and gravity field determination. The models that describe the accelerations are evaluated along a one day GRACE orbit arc (integrated for 3 July 2008 using Encke's method, see Ellmer and Mayer-Gürr (2017) for more information on the integration) using the GROOPS software. The orbit has a 30 s sampling and is given in cartesian coordinates expressed in both Terrestrial Reference Frame (TRF), where ITRF2014 (Altamimi et al., 2016) is used, and the Celestial Reference Frame (CRF), where ICRF2 (Fey et al., 2009) is taken. The evaluations yield three dimensional accelerations for each sampling point at the location of the spacecraft. These accelerations are the main product of the benchmark data set. The underlying models are listed in Table 1.

Table 1 Accelerations modelled in the benchmark data set. The magnitude is a rough estimation for the given GRACE orbit, the list is ordered according to the magnitude of the forces.

¹ Förste et al. (2014). ² Folkner et al. (2009). ³ Petit and Luzum (2010). ⁴ Savcenko and Bosch (2011). ⁵ Carrere et al. (2016). ⁶ Dobslaw et al. (2017). ⁷ Desai (2002).

All resulting accelerations are expressed in the CRF, the gravity field is additionally given in TRF. For each acceleration one file is provided. Accelerations consisting of several components, such as ocean tides, are provided in separate files for the largest constituents. Each file provides time stamps in GPS time (expressed in Modified Julian Date [MJD]) and the respective x, y, z-accelerations. As the definition of frames may vary between different software packages (J2000, true system of epoch), the orbit is given in the terrestrial and the celestial reference frame, and additionally, the rotation between the two underlying frames is provided as quaternions and rotation matrices. In an attempt to also make the handling of the terrestrial and celestial reference frame comprehensible, as well as to enable an easier understanding of ocean and atmosphere tide models, several additional physical quantities are appended to the data set. The additional data is listed in Table 2. The complete data set can be found at the ftp server of IfG¹ including a description of the data and how the models are employed (see file 00README_simulation.txt). To illustrate the nature of the accelerations Fig. 1 shows the norm of each acceleration considered in the benchmark data in CRF. The mean orbital altitude is 466 km.

Figure 1 Norm of the benchmark accelerations.

Table 2 Additional data set

^* Bizouard et al. (2018).

2.1 Application of the benchmark data set

The main goal of the data is to provide a reference for basic software comparisons. It enables a comparison of the background force model implementations by evaluating the models at the given orbital positions in different software packages, and may serve as a reference for the handling of celestial and terrestrial reference frames. The most straight forward approach of comparison is to evaluate the force models with a given space geodesy software package and to print the resulting accelerations for each model. By subtracting the obtained accelerations from the benchmark data, differences may be revealed

$$\begin{array}{}\text{(1)}& \left[\begin{array}{c}{d}_x\\ d_y\\ d_z\end{array}\right]=\left[\begin{array}{c}{r}_x\\ r_y\\ r_z\end{array}\right]-\left[\begin{array}{c}{s}_x\\ s_y\\ s_z\end{array}\right].\end{array}$$

d denotes the difference between the reference and the software under evaluation, r the reference accelerations and s the accelerations resulting from the software's evaluation. Large differences to the reference accelerations may indicate potential implementation problems. It is very unlikely to obtain zero differences. Unless large systematic differences emerge, oscillating patterns around zero will usually be observed due to the orbital revolution, most commonly twice-per-revolution due to the Earth's oblateness at the poles and the almost polar reference orbit.

In the following we also use the maximum absolute deviation from the reference (d_max) to summarise the comparisons of several software packages

$$\begin{array}{}\text{(2)}& d_{\max }=\max \left(\right|\mathit{r}-\mathit{s}\left|\right).\end{array}$$

3 Summary of each force

In this section we give a few notes on each force in the benchmark data set. All formulae correspond with the IERS 2010 conventions (Petit and Luzum, 2010), however, they are provided in the way they are coded in the Bernese GNSS software. All constants can be found in the IERS 2010 conventions and the corresponding model descriptions. When referring to the equations from the IERS 2010 conventions we use the notation Eq.^IERS (n.n) to provide good reading flow.

3.1 Earth's Gravity field

The gravity field model used for the benchmark data set is EIGEN-6C4 (Förste et al., 2014), which contains a set of dimensionless, fully-normalised spherical harmonic coefficients representing the static part of the Earth's gravity field. In the benchmark accelerations it is evaluated from degree and order 2 to 180. In contrast to all other forces, the data set contains the accelerations expressed in both CRF and TRF. The Earth's potential and the coefficients are related by

$$\begin{array}{}\text{(3)}& \begin{array}{rl}V(r,\mathit{\lambda },\mathit{\vartheta })& ={\displaystyle \frac{{\mathrm{GM}}_{\mathrm{E}}}{r}}\sum _{n=\mathrm{0}}^N{\left({\displaystyle \frac{a_{\mathrm{E}}}{r}}\right)}^n\sum _{m=\mathrm{0}}^n\left[c_{nm}\cos \right(m\mathit{\lambda })\\ & +s_{nm}\sin \left(m\mathit{\lambda }\right)\left]P_{nm}\right(\cos \mathit{\vartheta }),\end{array}\end{array}$$

where c_nm and s_nm denote the normalised spherical harmonic coefficients of degree n and order m and P_nm are the normalised associated Legendre functions. The point of calculating the potential is defined by the spherical coordinates with the radius r, the co-latitude ϑ and the geographical longitude λ. a_E defines the reference radius of the Earth, GM_E the gravitational constant multiplied with the mass of the Earth. Both are associated with the gravity field model and can be found in the header of EIGEN-6C4 data file. Equation 3 corresponds to Eq.^IERS (6.1).

The accelerations and the potential are related by

$$\begin{array}{}\text{(4)}& {\mathit{a}}_{r,\mathit{\lambda },\mathit{\vartheta }}=\mathrm{\nabla }V,\end{array}$$

with ∇ being the gradient operator. The accelerations in cartesian coordinates can be obtained by a transformation of Earth-fixed ($r,\mathit{\lambda },\mathit{\vartheta }$)-frame to ($x,y,z$)-frame. We refer to Heiskanen and Moritz (1967) for a detailed derivation of the expression of the Earth's gravity field in spherical harmonics.

3.2 3rd body attractions

The influence of other celestial bodies (denoted by the subscript cb) than the Earth is computed with positions derived from the JPL DE421 (Folkner et al., 2009) ephemeris. All bodies are treated as point masses, thus,

$$\begin{array}{}\text{(5)}& {\mathit{a}}_{\mathrm{cb}}=-{\mathrm{GM}}_{\mathrm{cb}}\left({\displaystyle \frac{{\mathit{r}}_{\mathrm{sat}}-{\mathit{r}}_{\mathrm{cb}}}{|{\mathit{r}}_{\mathrm{sat}}-{\mathit{r}}_{\mathrm{cb}}{|}^{\mathrm{3}}}}+{\displaystyle \frac{{\mathit{r}}_{\mathrm{cb}}}{r_{\mathrm{cb}}^{\mathrm{3}}}}\right)\end{array}$$

expresses the accelerations caused by a celestial body. r_sat and r_cb denote the geocentric position vector of the satellite and the celestial body. The constants of GM_cb are given in the readme-file. Planets up to Saturn are taken into account, however, as tests showed, additionally introducing Neptune, Uranus and Pluto leads to negligible small differences. The positions of Sun and Moon, computed from the DE421 ephemeris, are also used for the Solid Earth tides and the relativistic corrections.

3.3 Solid Earth tides

Solid Earth tides are computed according to IERS 2010 conventions using the anelastic model. Besides fundamental quantities of the Earth, the solid Earth tides depend on the position of the Sun, the Moon and the load Love numbers. They affect the spherical harmonic spectrum up to degree four.

According to the IERS 2010 conventions the computation is divided into two steps. Step 1 computes the coefficients due to the tide generating potential for degree 2 and 3, as well as the effect of degree 2 on degree 4 coefficients. Step 1 is frequency independent, whereas step 2 states frequency dependent corrections for degree two.

Step 1 – corresponds with Eqs.^IERS (6.6) and (6.7):

$$\begin{array}{}\text{(6)}& {\displaystyle }\begin{array}{rl}{c}_{\mathrm{2}m}& =\mathit{\left\{}\right[k_{\mathrm{2}m}^{\mathrm{Re}}\cos \left(m{\mathit{\lambda }}_{\mathrm{M}}\right)+k_{\mathrm{2}m}^{\mathrm{Im}}\sin \left(m{\mathit{\lambda }}_{\mathrm{M}}\right)\left]P_{\mathrm{2}m}\right(\cos {\mathit{\vartheta }}_{\mathrm{M}})\\ & {\displaystyle \frac{{\mathrm{GM}}_{\mathrm{M}}}{\mathrm{5}{\mathrm{GM}}_{\mathrm{E}}}}{\left({\displaystyle \frac{a_{\mathrm{E}}}{r_{\mathrm{M}}}}\right)}^{\mathrm{3}}\mathit{\}}+\mathit{\{}\left[k_{\mathrm{2}m}^{\mathrm{Re}}\cos \right(m{\mathit{\lambda }}_{\mathrm{S}})+k_{\mathrm{2}m}^{\mathrm{Im}}\sin (m{\mathit{\lambda }}_{\mathrm{S}}\left)\right]\\ & P_{\mathrm{2}m}\left(\cos {\mathit{\vartheta }}_{\mathrm{S}}\right){\displaystyle \frac{{\mathrm{GM}}_{\mathrm{S}}}{\mathrm{5}{\mathrm{GM}}_{\mathrm{E}}}}{\left({\displaystyle \frac{a_{\mathrm{E}}}{r_{\mathrm{S}}}}\right)}^{\mathrm{3}}\mathit{\}},\\ & m\in \mathit{\{}\mathrm{0},\mathrm{1},\mathrm{2}\mathit{\}}\\ & s_{\mathrm{2}m}=\mathit{\left\{}\right[k_{\mathrm{2}m}^{\mathrm{Re}}\sin \left(m{\mathit{\lambda }}_{\mathrm{M}}\right)-k_{\mathrm{2}m}^{\mathrm{Im}}\cos \left(m{\mathit{\lambda }}_{\mathrm{M}}\right)\left]P_{\mathrm{2}m}\right(\cos {\mathit{\vartheta }}_{\mathrm{M}})\\ & {\displaystyle \frac{{\mathrm{GM}}_{\mathrm{M}}}{\mathrm{5}{\mathrm{GM}}_{\mathrm{E}}}}{\left({\displaystyle \frac{a_{\mathrm{E}}}{r_{\mathrm{M}}}}\right)}^{\mathrm{3}}\mathit{\}}+\mathit{\{}\left[k_{\mathrm{2}m}^{\mathrm{Re}}\sin \right(m{\mathit{\lambda }}_{\mathrm{S}})-k_{\mathrm{2}m}^{\mathrm{Im}}\cos (m{\mathit{\lambda }}_{\mathrm{S}}\left)\right]\\ & P_{\mathrm{2}m}\left(\cos {\mathit{\vartheta }}_{\mathrm{S}}\right){\displaystyle \frac{{\mathrm{GM}}_{\mathrm{S}}}{\mathrm{5}{\mathrm{GM}}_{\mathrm{E}}}}{\left({\displaystyle \frac{a_{\mathrm{E}}}{r_{\mathrm{S}}}}\right)}^{\mathrm{3}}\mathit{\}},\\ & m\in \mathit{\{}\mathrm{0},\mathrm{1},\mathrm{2}\mathit{\}}\end{array}\text{(7)}& {\displaystyle }\begin{array}{rl}{c}_{\mathrm{4}m}& =\mathit{\left\{}\right[k_{\mathrm{2}m}^{(+)}\cos \left(m{\mathit{\lambda }}_{\mathrm{M}}\right)\left]P_{\mathrm{2}m}\right(\cos {\mathit{\vartheta }}_{\mathrm{M}}\left){\displaystyle \frac{{\mathrm{GM}}_{\mathrm{M}}}{\mathrm{5}{\mathrm{GM}}_{\mathrm{E}}}}{\left({\displaystyle \frac{a_{\mathrm{E}}}{r_{\mathrm{M}}}}\right)}^{\mathrm{3}}\mathit{\right\}}\\ & +\mathit{\left\{}\right[k_{\mathrm{2}m}^{(+)}\cos \left(m{\mathit{\lambda }}_{\mathrm{S}}\right)\left]P_{\mathrm{2}m}\right(\cos {\mathit{\vartheta }}_{\mathrm{S}}\left){\displaystyle \frac{{\mathrm{GM}}_{\mathrm{S}}}{\mathrm{5}{\mathrm{GM}}_{\mathrm{E}}}}{\left({\displaystyle \frac{a_{\mathrm{E}}}{r_{\mathrm{S}}}}\right)}^{\mathrm{3}}\mathit{\right\}},\\ & m\in \mathit{\{}\mathrm{0},\mathrm{1},\mathrm{2}\mathit{\}}\\ & s_{\mathrm{4}m}=\mathit{\left\{}\right[k_{\mathrm{2}m}^{(+)}\sin \left(m{\mathit{\lambda }}_{\mathrm{M}}\right)\left]P_{\mathrm{2}m}\right(\cos {\mathit{\vartheta }}_{\mathrm{M}}\left){\displaystyle \frac{{\mathrm{GM}}_{\mathrm{M}}}{\mathrm{5}{\mathrm{GM}}_{\mathrm{E}}}}{\left({\displaystyle \frac{a_{\mathrm{E}}}{r_{\mathrm{M}}}}\right)}^{\mathrm{3}}\mathit{\right\}}\\ & +\mathit{\left\{}\right[k_{\mathrm{2}m}^{(+)}\sin \left(m{\mathit{\lambda }}_{\mathrm{S}}\right)\left]P_{\mathrm{2}m}\right(\cos {\mathit{\vartheta }}_{\mathrm{S}}\left){\displaystyle \frac{{\mathrm{GM}}_{\mathrm{S}}}{\mathrm{5}{\mathrm{GM}}_{\mathrm{E}}}}{\left({\displaystyle \frac{a_{\mathrm{E}}}{r_{\mathrm{S}}}}\right)}^{\mathrm{3}}\mathit{\right\}},\\ & m\in \mathit{\{}\mathrm{0},\mathrm{1},\mathrm{2}\mathit{\}}\end{array}\text{(8)}& {\displaystyle }\begin{array}{rl}{c}_{\mathrm{3}m}& =\mathit{\left\{}\right[k_{\mathrm{3}m}^{\mathrm{Re}}\cos \left(m{\mathit{\lambda }}_{\mathrm{M}}\right)+k_{\mathrm{3}m}^{\mathrm{Im}}\sin \left(m{\mathit{\lambda }}_{\mathrm{M}}\right)\left]P_{\mathrm{3}m}\right(\cos {\mathit{\vartheta }}_{\mathrm{M}})\\ & {\displaystyle \frac{{\mathrm{GM}}_{\mathrm{M}}}{\mathrm{7}{\mathrm{GM}}_{\mathrm{E}}}}{\left({\displaystyle \frac{a_{\mathrm{E}}}{r_{\mathrm{M}}}}\right)}^{\mathrm{4}}\mathit{\}}+\mathit{\{}\left[k_{\mathrm{3}m}^{\mathrm{Re}}\cos \right(m{\mathit{\lambda }}_{\mathrm{S}})+k_{\mathrm{3}m}^{\mathrm{Im}}\sin (m{\mathit{\lambda }}_{\mathrm{S}}\left)\right]\\ & P_{\mathrm{3}m}\left(\cos {\mathit{\vartheta }}_{\mathrm{S}}\right){\displaystyle \frac{{\mathrm{GM}}_{\mathrm{S}}}{\mathrm{7}{\mathrm{GM}}_{\mathrm{E}}}}{\left({\displaystyle \frac{a_{\mathrm{E}}}{r_{\mathrm{S}}}}\right)}^{\mathrm{4}}\mathit{\}},\\ & m\in \mathit{\{}\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3}\mathit{\}}\\ s_{\mathrm{3}m}& =\mathit{\left\{}\left[k_{\mathrm{3}m}^{\mathrm{Re}}\sin \left(m{\mathit{\lambda }}_{\mathrm{M}}\right)-k_{\mathrm{3}m}^{\mathrm{Im}}\cos \left(m{\mathit{\lambda }}_{\mathrm{M}}\right)\right]P_{\mathrm{3}m}\right(\cos {\mathit{\vartheta }}_{\mathrm{M}})\\ & {\displaystyle \frac{{\mathrm{GM}}_{\mathrm{M}}}{\mathrm{7}{\mathrm{GM}}_{\mathrm{E}}}}{\left({\displaystyle \frac{a_{\mathrm{E}}}{r_{\mathrm{M}}}}\right)}^{\mathrm{4}}\mathit{\}}+\mathit{\{}\left[k_{\mathrm{3}m}^{\mathrm{Re}}\sin \left(m{\mathit{\lambda }}_{\mathrm{S}}\right)-k_{\mathrm{3}m}^{\mathrm{Im}}\cos \left(m{\mathit{\lambda }}_{\mathrm{S}}\right)\right]\\ & P_{\mathrm{3}m}\left(\cos {\mathit{\vartheta }}_{\mathrm{S}}\right){\displaystyle \frac{{\mathrm{GM}}_{\mathrm{S}}}{\mathrm{7}{\mathrm{GM}}_{\mathrm{E}}}}{\left({\displaystyle \frac{a_{\mathrm{E}}}{r_{\mathrm{S}}}}\right)}^{\mathrm{4}}\mathit{\}},\\ & m\in \mathit{\{}\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3}\mathit{\}}\end{array}\end{array}$$

$k_{nm}^{\mathrm{Re}}$, $k_{nm}^{\mathrm{Im}}$, $k_{nm}^{(+)}$ denote the Love numbers (Table^IERS 6.3), r_S and r_M is the norm of the geocentric vector (r_S, r_M) to Sun and Moon, the Legendre functions P_nm depend on the cosine of the co-latitude ϑ_S,M of Sun and Moon in the TRF. λ_S,M is the geographical longitude of Sun and Moon in the Earth-fixed frame.

Step 2: Corrections on c₂₀ – corresponds to Eq.^IERS (6.8a) and uses Table^IERS 6.5b:

$$\begin{array}{}\text{(9)}& \mathrm{\Delta }{c}_{\mathrm{20}}^{\mathrm{corr}}=\sum _{f_{\text{long periodic}}}(A_{\mathrm{f}}^{\mathrm{ip}}\cos {\mathit{\theta }}_f-A_{\mathrm{f}}^{\mathrm{op}}\sin {\mathit{\theta }}_f)\end{array}$$

Corrections on c₂₁, c₂₂, s₂₁, s₂₂ – corresponds to Eq.^IERS (6.8b) and uses Table^IERS 6.5a and Table^IERS 6.5c:

$$\begin{array}{}\text{(10)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }{c}_{\mathrm{21}}^{\text{corr}}& =\sum _{f_{\text{diurnal}}}(A_{\mathrm{f}}^{\mathrm{ip}}\sin {\mathit{\theta }}_f+A_{\mathrm{f}}^{\mathrm{op}}\cos {\mathit{\theta }}_f)\\ & \mathrm{\Delta }{s}_{\mathrm{21}}^{\mathrm{corr}}=\sum _{f_{\text{diurnal}}}(A_{\mathrm{f}}^{\mathrm{ip}}\cos {\mathit{\theta }}_f-A_{\mathrm{f}}^{\mathrm{op}}\sin {\mathit{\theta }}_f)\end{array}\text{(11)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }{c}_{\mathrm{22}}^{\text{corr}}& =\sum _{f_{\text{semi diurnal}}}+A_{\mathrm{f}}^{\mathrm{ip}}\cos {\mathit{\theta }}_{\mathrm{f}}\\ & \mathrm{\Delta }{s}_{\mathrm{22}}^{\text{corr}}=\sum _{f_{\text{semi diurnal}}}-A_{\mathrm{f}}^{\mathrm{ip}}\sin {\mathit{\theta }}_{\mathrm{f}}\end{array}\end{array}$$

The amplitudes for the frequency dependent corrections for $A_{\mathrm{f}}^{\mathrm{ip}}$ and $A_{\mathrm{f}}^{\mathrm{op}}$ are listed in Tables^IERS 6.5a, 6.5b and 6.5c.

The Doodson angle argument reads as

$$\begin{array}{}\text{(12)}& {\mathit{\theta }}_{\mathrm{f}}=\langle {\mathit{n}}_{\mathrm{f}},\mathit{\beta }\rangle ,\end{array}$$

with θ_f being computed from the fundamental Doodson arguments β and the respective tidal frequency n_f. Both vectors have six elements. The computation of the fundamental Doodson arguments $\mathit{\beta }=[\mathit{\tau },s,h,p,N^{\prime },p_s]$ follows from the fundamental arguments of lunisolar nutation (Delaunay variables, cf.  Doodson and Lamb, 1921) $l,l^{\prime },F,D,\mathrm{\Omega }$ (Eq.^IERS (5.43)).

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& l={\mathrm{134.96340251}}^{\circ }+{\mathrm{1717915923.2178}}^{\prime \prime }t\\ & +{\mathrm{31.8792}}^{\prime \prime }{t}^{\mathrm{2}}+{\mathrm{0.051635}}^{\prime \prime }{t}^{\mathrm{3}}-{\mathrm{0.00024470}}^{\prime \prime }{t}^{\mathrm{4}}\\ & l^{\prime }={\mathrm{357.52910918}}^{\circ }+{\mathrm{129596581.0481}}^{\prime \prime }t-{\mathrm{0.5532}}^{\prime \prime }{t}^{\mathrm{2}}\\ & +{\mathrm{0.000136}}^{\prime \prime }{t}^{\mathrm{3}}-{\mathrm{0.00001149}}^{\prime \prime }{t}^{\mathrm{4}}\\ & F={\mathrm{93.27209062}}^{\circ }+{\mathrm{1739527262.8478}}^{\prime \prime }t-{\mathrm{12.7512}}^{\prime \prime }{t}^{\mathrm{2}}\\ & -{\mathrm{0.001037}}^{\prime \prime }{t}^{\mathrm{3}}+{\mathrm{0.00000417}}^{\prime \prime }{t}^{\mathrm{4}}\\ & D={\mathrm{297.85019547}}^{\circ }+{\mathrm{1602961601.2090}}^{\prime \prime }t-{\mathrm{6.3706}}^{\prime \prime }{t}^{\mathrm{2}}\\ & +{\mathrm{0.006593}}^{\prime \prime }{t}^{\mathrm{3}}-{\mathrm{0.00003169}}^{\prime \prime }{t}^{\mathrm{4}}\\ & \mathrm{\Omega }={\mathrm{125.04455501}}^{\circ }-{\mathrm{6962890.5431}}^{\prime \prime }t+{\mathrm{7.4722}}^{\prime \prime }{t}^{\mathrm{2}}\\ & +{\mathrm{0.007702}}^{\prime \prime }{t}^{\mathrm{3}}-{\mathrm{0.00005939}}^{\prime \prime }{t}^{\mathrm{4}}\end{array}\end{array}$$

t denotes the time interval between the current epoch and J2000.0 in Julian centuries, thus,

$$\begin{array}{}\text{(14)}& t={\displaystyle \frac{t_{\mathrm{MJD}}-t_{{\mathrm{MJD}}_{\mathrm{2000}}}}{\mathrm{36525}}},\text{with}{t}_{{\mathrm{MJD}}_{\mathrm{2000}}}=\mathrm{51544.5}.\end{array}$$

Finally, the Doodson arguments read as

$$\begin{array}{}\text{(15)}& \mathit{\beta }=\left[\begin{array}{c}\mathit{\tau }\\ s\\ h\\ p\\ N^{\prime }\\ p_s\end{array}\right]=\left[\begin{array}{c}{\mathit{\theta }}_{\mathrm{g}}+\mathit{\pi }-s\\ F+\mathrm{\Omega }\\ s-D\\ s-l\\ -\mathrm{\Omega }\\ s-D-l^{\prime }\end{array}\right],\end{array}$$

where θ_g denotes Greenwich Mean Sidereal Time (GMST) in angle units.

In the tidal frequency vector n_f each digit makes up one element of the multipliers n_i, the representation is hexadecimal to allow for frequencies greater than 999.999. An example for the long periodic om1-tide with a frequency of 55.565 would be ${\mathit{n}}_{om\mathrm{1}}=[\mathrm{0},\mathrm{5},\mathrm{5},\mathrm{5},\mathrm{6},\mathrm{5}]$. The first element of the vector tells the periodicity of the tide.

The final coefficients $c_{\mathrm{2}m}^{\text{tot}}$ and $s_{\mathrm{2}m}^{\text{tot}}$ are obtained by adding together step 1 and step 2

$$\begin{array}{}\text{(16)}& \begin{array}{r}{c}_{\mathrm{2}m}^{\text{tot}}=c_{\mathrm{2}m}+\mathrm{\Delta }{c}_{\mathrm{2}m}^{\text{corr}}\\ s_{\mathrm{2}m}^{\text{tot}}=s_{\mathrm{2}m}+\mathrm{\Delta }{s}_{\mathrm{2}m}^{\text{corr}}\end{array}.\end{array}$$

In order to provide all information needed for the evaluation, the Doodson arguments and fundamental arguments of nutation, computed for every position of the reference orbit, are also part of the data set. They are also needed for the computation of ocean and atmospheric tides. The accelerations caused by the solid Earth tides are eventually computed using Eqs. 3 and 4.

3.4 Ocean tides

Two different models are provided for the ocean tides: EOT11a (Savcenko and Bosch, 2011) and FES2014b (Carrere et al., 2016). Both are based on IfG's conversion of the corresponding grids to spherical harmonic coefficients. EOT11a is evaluated from degree 2 to 120, FES2014b from degree 2 to 180. EOT11a consists of 18 tidal frequencies, FES2014b of 34. Furthermore, the data set provides linear admittances, a short description and a MATLAB routine of their application. The spherical harmonics coefficients can be computed from the given prograde ($c_{nm,f}^{\text{cos}}$, $s_{nm,f}^{\text{cos}}$) and retrograde ($c_{nm,f}^{\text{sin}}$, $s_{nm,f}^{\text{sin}}$) coefficients in a sum over all tidal frequencies f using

$$\begin{array}{}\text{(17)}& \begin{array}{r}{c}_{nm}=\sum _f\left[\cos \left({\mathit{\theta }}_f\right)c_{nm,f}^{\text{cos}}+\sin \left({\mathit{\theta }}_f\right)c_{nm,f}^{\text{sin}}\right]\\ s_{nm}=\sum _f\left[\cos \left({\mathit{\theta }}_f\right)s_{nm,f}^{\text{cos}}+\sin \left({\mathit{\theta }}_f\right)s_{nm,f}^{\text{sin}}\right]\end{array},\end{array}$$

In this representation the Doodson-Warburg correction (see Table^IERS 6.6) is already applied.

To complete the tidal spectrum, admittances between the major tides can be computed using linear interpolation

$$\begin{array}{}\text{(18)}& \begin{array}{r}{c}_{nm,f}^{\text{cos}}=c_{nm,\mathrm{1}}^{\text{cos}}{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{2}}-t_f}{{\mathit{\theta }}_{\mathrm{2}}-t_{\mathrm{1}}}}{\displaystyle \frac{H_f}{H_{\mathrm{1}}}}+c_{nm,\mathrm{2}}^{\text{cos}}{\displaystyle \frac{{\mathit{\theta }}_f-{\mathit{\theta }}_{\mathrm{1}}}{{\mathit{\theta }}_{\mathrm{2}}-{\mathit{\theta }}_{\mathrm{1}}}}{\displaystyle \frac{H_f}{H_{\mathrm{2}}}}\\ s_{nm,f}^{\text{cos}}=s_{nm,\mathrm{1}}^{\text{cos}}{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{2}}-{\mathit{\theta }}_f}{{\mathit{\theta }}_{\mathrm{2}}-{\mathit{\theta }}_{\mathrm{1}}}}{\displaystyle \frac{H_f}{H_{\mathrm{1}}}}+s_{nm,\mathrm{2}}^{\text{cos}}{\displaystyle \frac{{\mathit{\theta }}_f-{\mathit{\theta }}_{\mathrm{1}}}{{\mathit{\theta }}_{\mathrm{2}}-{\mathit{\theta }}_{\mathrm{1}}}}{\displaystyle \frac{H_f}{H_{\mathrm{2}}}}\\ c_{nm,f}^{\text{sin}}=c_{nm,\mathrm{1}}^{\text{sin}}{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{2}}-t_f}{{\mathit{\theta }}_{\mathrm{2}}-t_{\mathrm{1}}}}{\displaystyle \frac{H_f}{H_{\mathrm{1}}}}+c_{nm,\mathrm{2}}^{\text{sin}}{\displaystyle \frac{{\mathit{\theta }}_f-{\mathit{\theta }}_{\mathrm{1}}}{{\mathit{\theta }}_{\mathrm{2}}-{\mathit{\theta }}_{\mathrm{1}}}}{\displaystyle \frac{H_f}{H_{\mathrm{2}}}}\\ s_{nm,f}^{\text{sin}}=s_{nm,\mathrm{1}}^{\text{sin}}{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{2}}-{\mathit{\theta }}_f}{{\mathit{\theta }}_{\mathrm{2}}-{\mathit{\theta }}_{\mathrm{1}}}}{\displaystyle \frac{H_f}{H_{\mathrm{1}}}}+s_{nm,\mathrm{2}}^{\text{sin}}{\displaystyle \frac{{\mathit{\theta }}_f-{\mathit{\theta }}_{\mathrm{1}}}{{\mathit{\theta }}_{\mathrm{2}}-{\mathit{\theta }}_{\mathrm{1}}}}{\displaystyle \frac{H_f}{H_{\mathrm{2}}}}\end{array}.\end{array}$$

The subscripts 1 and 2 denote the main waves, f the interpolated wave. H is the astronomic amplitude of the wave. The ocean tides are given in ICGEM-format (Barthelmes and Förste, 2011, “.gfc”) and IERS-format. The latter provides a smaller number of significant digits for the coefficients, hence, we expect differences at the level of 10⁻¹¹ m s⁻² when using the two different file formats. The tidal frequencies and the phase shift applied (Doodson-Warburg correction) are stated in the readme-file. For the evaluation using the IERS-format (potential and water height) we refer to Petit and Luzum (2010, chap. 6.3.1). Accelerations resulting from ocean tides are eventually obtained by Eqs. (3) and (4).

3.5 Relativistic corrections

Relativistic corrections are computed according to the IERS 2010 conventions for General Relativity using Eq.^IERS (10.12) and contain the Schwarzschild (Eq. 19), Lense-Thirring (Eq. 20) and de Sitter term (Eq. 21).

$$\begin{array}{}\text{(19)}& {\displaystyle }\begin{array}{rl}{\mathit{a}}_{\text{Schwarzschild}}& ={\displaystyle \frac{{\mathrm{GM}}_{\mathrm{E}}}{c^{\mathrm{2}}{r}_{\mathrm{sat}}^{\mathrm{3}}}}\left[\left(\mathrm{4}{\displaystyle \frac{{\mathrm{GM}}_{\mathrm{E}}}{r_{\mathrm{sat}}}}-{\mathit{v}}_{\mathrm{sat}}\cdot {\mathit{v}}_{\mathrm{sat}}\right){\mathit{r}}_{\mathrm{sat}}\right.\\ & \left.+\mathrm{4}\left({\mathit{r}}_{\mathrm{sat}}\cdot {\mathit{v}}_{\mathrm{sat}}\right){\mathit{v}}_{\mathrm{sat}}\right]\end{array}\text{(20)}& {\displaystyle }\begin{array}{rl}{\mathit{a}}_{\text{Lense Thirring}}& =\mathrm{2}{\displaystyle \frac{{\mathrm{GM}}_{\mathrm{E}}}{c^{\mathrm{2}}{r}_{\mathrm{sat}}^{\mathrm{3}}}}\left[{\displaystyle \frac{\mathrm{3}}{r_{\mathrm{sat}}^{\mathrm{2}}}}({\mathit{r}}_{\mathrm{sat}}\times {\mathit{v}}_{\mathrm{sat}})\left({\mathit{r}}_{\mathrm{sat}}\cdot \mathit{J}\right)\right.\\ & \left.+({\mathit{v}}_{\mathrm{sat}}\times \mathit{J})\right]\end{array}\text{(21)}& {\displaystyle }{\mathit{a}}_{\text{de Sitter}}=\mathrm{3}\left[-{\mathit{v}}_S\times \left({\displaystyle \frac{{\mathrm{GM}}_{\mathrm{S}}}{c^{\mathrm{2}}{r}_S^{\mathrm{3}}}}{\mathit{r}}_S\right)\right]\times {\mathit{v}}_{\mathrm{sat}}\end{array}$$

r_sat and v_sat are the geocentric position and the velocity of the satellite. c denotes the speed of light, J is the Earth's angular momentum per unit mass and can be set to ${\left[\begin{array}{c}\mathrm{0},\mathrm{0},\mathrm{9.8}\times {\mathrm{10}}^{\mathrm{8}}\end{array}\right]}^{\mathrm{T}}$ m² s⁻¹ (Petit and Luzum, 2010). The vectors r_S and v_S describe the geocentric position and velocity of the Sun in CRF. The relativistic accelerations in the benchmark data are the sum of the three components

$$\begin{array}{}\text{(22)}& {\mathit{a}}_{\text{relativistic}}={\mathit{a}}_{\text{Schwarzschild}}+{\mathit{a}}_{\text{Lense Thirring}}+{\mathit{a}}_{\text{de Sitter}}.\end{array}$$

3.6 Dealiasing

AOD1B RL06 (Dobslaw et al., 2017) is used as dealiasing product. The “glo” dataset, which is the sum of the atmospheric and oceanic contribution, is evaluated in the benchmark data. The degree one coefficients are not taken into account, it is evaluated from degree 2 to 180. The spherical harmonic synthesis of the dealiasing model follows Eqs. 3 and 4, additionally, as the data set is given in three hour sets, a linear interpolation between the neighbouring sets at time t₁ and t₂ is performed on the level of spherical harmonic coefficients to obtain the spherical harmonic coefficients at time t_i ($t_{\mathrm{1}}\le t_i\le t_{\mathrm{2}}$) (Eq. 23).

$$\begin{array}{}\text{(23)}& \begin{array}{r}{c}_{nm}=c_{nm}\left(t_{\mathrm{1}}\right){\displaystyle \frac{t_{\mathrm{2}}-t_i}{t_{\mathrm{2}}-t_{\mathrm{1}}}}+c_{nm}\left(t_{\mathrm{2}}\right){\displaystyle \frac{t_i-t_{\mathrm{1}}}{t_{\mathrm{2}}-t_{\mathrm{1}}}}\\ s_{nm}=s_{nm}\left(t_{\mathrm{1}}\right){\displaystyle \frac{t_{\mathrm{2}}-t_i}{t_{\mathrm{2}}-t_{\mathrm{1}}}}+s_{nm}\left(t_{\mathrm{2}}\right){\displaystyle \frac{t_i-t_{\mathrm{1}}}{t_{\mathrm{2}}-t_{\mathrm{1}}}}\end{array}\end{array}$$

3.7 Pole tide

“The pole tide of the solid Earth is generated by the centrifugal effect of polar motion” (Petit and Luzum, 2010, chap. 6.4). Only the coefficients c₂₁ and s₂₁ are concerned and the contribution is computed from the polar motion parameters and the mean pole definition via

$$\begin{array}{}\text{(24)}& \begin{array}{r}{c}_{\mathrm{21}}=-\mathrm{1.333}\times {\mathrm{10}}^{-\mathrm{9}}\left[\left(x_{\mathrm{P}}-{\overline{x}}_{\mathrm{P}}\right)-\mathrm{0.0115}\left(y_{\mathrm{P}}-{\overline{y}}_{\mathrm{P}}\right)\right]\\ s_{\mathrm{21}}=+\mathrm{1.333}\times {\mathrm{10}}^{-\mathrm{9}}\left[\left(y_{\mathrm{P}}-{\overline{y}}_{\mathrm{P}}\right)+\mathrm{0.0115}\left(x_{\mathrm{P}}-{\overline{x}}_{\mathrm{P}}\right)\right]\end{array}\end{array}$$

x_P, y_P, ${\overline{x}}_{\mathrm{P}}$ and ${\overline{y}}_{\mathrm{P}}$ form the wobble parameters (Eq.^IERS (7.24)) via

$$\begin{array}{}\text{(25)}& m_{\mathrm{1}}=(x_{\mathrm{P}}-{\overline{x}}_{\mathrm{P}}),m_{\mathrm{2}}=-(y_{\mathrm{P}}-{\overline{y}}_{\mathrm{P}}).\end{array}$$

The pole tide and ocean pole tide computed in the benchmark data make use of the most recent secular pole (IERS, 2018)

$$\begin{array}{}\text{(26)}& \begin{array}{rl}& {\overline{x}}_{\mathrm{P}}=\mathrm{55.0}\times {\mathrm{10}}^{-\mathrm{3}}+\mathrm{1.677}\times {\mathrm{10}}^{-\mathrm{3}}t,\\ & {\overline{y}}_{\mathrm{P}}=\mathrm{320.5}\times {\mathrm{10}}^{-\mathrm{3}}+\mathrm{3.460}\times {\mathrm{10}}^{-\mathrm{3}}t,\end{array}\end{array}$$

in units of arc seconds with t denoting the time interval between the current epoch and J2000.0 in Julian years, thus,

$$\begin{array}{}\text{(27)}& t={\displaystyle \frac{t_{\mathrm{MJD}}-t_{{\mathrm{MJD}}_{\mathrm{2000}}}}{\mathrm{365.25}}},\text{with}{t}_{{\mathrm{MJD}}_{\mathrm{2000}}}=\mathrm{51544.5}.\end{array}$$

The relation between the spherical harmonic coefficients and the accelerations is again given by Eqs. 3 and 4.

3.8 Atmospheric tides

Atmospheric tides are modelled using AOD1B RL06 product and contain all twelve tidal constituents from degree 2 to 180. Additionally, the accelerations caused by S1 tide only are given separately. The evaluation of the atmospheric tides follows

$$\begin{array}{}\text{(28)}& \begin{array}{r}{c}_{nm}=\sum _f\left[\cos ({\mathit{\theta }}_f+{\mathit{\chi }}_f)c_{nm,f}^{\text{cos}}+\sin ({\mathit{\theta }}_f+{\mathit{\chi }}_f)c_{nm,f}^{\text{sin}}\right]\\ s_{nm}=\sum _f\left[\cos ({\mathit{\theta }}_f+{\mathit{\chi }}_f)s_{nm,f}^{\text{cos}}+\sin ({\mathit{\theta }}_f+{\mathit{\chi }}_f)s_{nm,f}^{\text{sin}}\right]\end{array},\end{array}$$

where the coefficients $c_{nm,f}^{\text{cos}}$, $c_{nm,f}^{\text{sin}}$, $s_{nm,f}^{\text{cos}}$, $s_{nm,f}^{\text{sin}}$ are given by the model for each tidal frequency f. The angle argument θ_f follows from Eq. (12), the Doodson-Warburg phase correction χ_f depends on the tidal frequency. To obtain accelerations, the spherical harmonic coefficients can be evaluated by using Eqs. (3) and (4).

The model is given in ICGEM-format as well. Using this format, the application follows the method described for the ocean tides (see Eq. 17). For the evaluation using the IERS-format (potential) we refer to Petit and Luzum (2010, chap. 6.3.1).

3.9 Ocean pole tide

Similar to the solid Earth pole tide, the ocean pole tide is a result of the centrifugal effect of polar motion on the oceans. The implementation for the benchmark data set follows the IERS 2010 conventions employing the Desai model (Desai, 2002), which is given in spherical harmonic coefficients, representing a self-consistent equilibrium model. The benchmark accelerations make use of the degrees and orders 2 to 180. Load Love numbers are given externally as the conventions only state a few of the low degree Love numbers. The formulae follow Eqs.^IERS (6.23a) and (6.23b).

$$\begin{array}{}\text{(29)}& \begin{array}{r}{c}_{nm}=R_n\left[A_{nm}^{\mathrm{Re}}\left(m_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{Re}}+m_{\mathrm{2}}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{Im}}\right)+A_{nm}^{\mathrm{Im}}\left(m_{\mathrm{2}}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{Re}}-m_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{Im}}\right)\right]\\ s_{nm}=R_n\left[B_{nm}^{\mathrm{Re}}\left(m_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{Re}}+m_{\mathrm{2}}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{Im}}\right)+B_{nm}^{\mathrm{Im}}\left(m_{\mathrm{2}}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{Re}}-m_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{Im}}\right)\right]\end{array},\end{array}$$

where $A_{nm}^{\mathrm{Re}}$, $A_{nm}^{\mathrm{Im}}$, $B_{nm}^{\mathrm{Re}}$, $B_{nm}^{\mathrm{Im}}$ are the coefficients from the model, m₁ and m₂ are the wobble parameters (Eq. 25), ${\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{Re}}=\mathrm{0.6870}$, ${\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{Im}}=\mathrm{0.0036}$ and the factor R_n is given by

$$\begin{array}{}\text{(30)}& R_n={\displaystyle \frac{{\mathit{\omega }}_{\mathrm{E}}^{\mathrm{2}}{a}_{\mathrm{E}}^{\mathrm{4}}\mathrm{4}\mathit{\pi }G\mathit{\rho }}{{\mathrm{GM}}_{\mathrm{E}}{g}_{\mathrm{eq}}}}\left({\displaystyle \frac{\mathrm{1}+k_n^{\prime }}{\mathrm{2}n+\mathrm{1}}}\right),\end{array}$$

with ω_E being the nominal mean Earth's rotation velocity, G the gravitational constant, ρ the density of sea water, g_eq the gravity at the equator and $k_n^{\prime }$ the load Love numbers. To obtain accelerations, the spherical harmonic coefficients are evaluated using Eqs. (3) and (4).

4 Comparisons within COST-G

The benchmark data set was created, used and examined within the COST-G initiative, where AIUB, GFZ, GRGS and IfG are currently acting as COST-G ACs. LUH is a candidate AC. To augment the combination effort and in particular to rule out large systematic differences in the implementation of background force models, all contributing groups performed a comparison with the benchmark data using their own software packages. This includes the Bernese GNSS software (AIUB), EPOS (GFZ), GINS (GRGS), GROOPS (IfG), GRACE-SIGMA (LUH) and GRASP (LUH). The intention of this comparison is to test the implementation and handling of the background models in each software package. Every software tries to reproduce the reference accelerations as closely as possible by employing the provided rotation between TRF and CRF. The GROOPS software serves as a reference as it was used to compute the benchmark data set.

4.1 Software packages

The software packages follow different approaches of modelling gravity fields from satellite data. Even though data is treated differently, we expect a high level of agreement with the benchmark data for background model handling. The following sub-sections give a brief introduction to each package.

4.1.1 Bernese GNSS software

The Bernese GNSS software is a scientific software package, used by more than 700 institutions around the world. It features space geodetic applications, mainly high-precision, multi-GNSS data processing for ground networks (e.g., Prange et al., 2017) as well as precise satellite orbit determination and thereof deriving gravity field solutions (e.g., Meyer et al., 2016). The software is developed, maintained and distributed at AIUB. In its core it employs the Celestial Mechanics Approach of orbit determination, be it for high or low flying Earth orbiting satellites or planetary geodetic applications (Arnold et al., 2015). The software is written in Fortran.

4.1.2 EPOS

The Earth Parameter and Orbit System (EPOS²) is a software package designed and used for operational Precise Orbit Determination (POD). It consists of tools for data orbit analysis, orbit integration, orbit improvement, orbit predictions, normal equation handling and simulation of observations, as well as for gravity field computation. It served and serves numerous satellite missions, such as Mir, Space Shuttle, LAGEOS, CHAMP, GRACE, GOCE, GRACE-FO, Swarm, TerraSAR-X, TanDEM-X, Envisat or Jason, being able to deal with SLR, GPS, DORIS, radar altimeter or GRACE K-Band-Ranging data. EPOS is based on the dynamic approach of orbit modelling and EPOS is one of the software of the current COST-G ACs that is written in Fortran. EPOS is developed, maintained and used at GFZ, and has also been installed at a few other institutions around Europe. Recent applications are the computation of monthly GRACE and GRACE-FO gravity field solutions (Dahle et al., 2019) or investigations on local ties for the datum realisation of global terrestrial reference frames (Glaser et al., 2019).

4.1.3 GINS

The GINS (Géodésie par Intégrations Numériques Simultanées) software is developed and maintained by the GRGS of the French space agency. It is a multi-technique space geodetic software, capable of processing data from GNSS, SLR, VLBI, DORIS and inter-satellite ranging. It is used for operational processing of all space geodetic observation techniques.

4.1.4 GRACE-SIGMA

The GRACE-SIGMA (GRACE-Satellite orbit Integration and Gravity field analysis in MAtlab) software is a recent development specifically designed for the processing of GRACE and GRACE-FO data. The software is developed at Institut für Erdmessung (IfE) of LUH. It is written entirely in MATLAB and uses strongly vectorised modules for modelling of disturbing forces, orbit propagation and orbit improvement. The integration of satellite ephemerides and state and sensitivity matrices is performed using an efficient in-house developed numerical integration approach (Naeimi, 2018). The gravity field estimation is based on a generalized dynamic orbit determination using variational equations. The software is applied for the computation of monthly GRACE and GRACE-FO gravity field solutions (Koch et al., 2020).

4.1.5 GRASP

The GRAvity Satellite Processing engine (GRASP) software is dedicated to gravity field recovery from kinematic positions of satellites (Weigelt et al., 2013). It uses the acceleration approach of orbit modelling, thus, the kinematic positions are numerically differentiated twice in order to make them express a quantity directly related to the force field that is acting on the satellite. GRASP is developed at LUH and written in MATLAB.

4.1.6 GROOPS

The Gravity Recovery Object Oriented Programming System (GROOPS) used at IfG is a software suite for geodetic applications. Its feature set includes the determination of GNSS orbits, clocks and ground station networks (Strasser et al., 2019), static and time-variable gravity field solutions from satellite data, and regional gravity field modeling with terrestrial data. GROOPS is written in C++ and makes heavy use of low level Basic Linear Algebra Subprograms (BLAS) and LAPACK (Linear Algebra PACKage) subroutines. It uses the Message Passing Interface (MPI) communication protocol for parallelization and is therefore capable to run on large distributed systems. GROOPS software serves as a reference within these software comparisons as the benchmark data set was compiled using the capabilities of GROOPS. Recent applications are the computation of the ITSG‐Grace2018 time series (Kvas et al., 2019) and operational processing of GRACE-FO data.

4.2 Example results using the Bernese GNSS Software

To show the results of the software comparison for all provided force models, we evaluate all force models with the Bernese GNSS software and compute the norm of the difference vector between the benchmark accelerations and the accelerations from Bernese. We expect the magnitude of the differences to be small as the handling of the background forces is supposed to be generally compatible. Figure 2 shows that the largest differences result for the solid Earth tides, the smallest for the relativistic effects. Generally, the magnitude of differences is close to the numerical precision of 64 bit arithmetic of about 16 decimal places (IEEE-754, 1985).

Figure 2 Norm of the difference between the benchmark data and the evaluation with the Bernese GNSS software.

Figure 3 Maximum absolute deviation of the difference between the reference and the evaluation by the respective software of the COST-G ACs for each force considered in the benchmark data set.

4.3 Results from COST-G

Each of the above mentioned institutions performed the software comparison with the benchmark data. EPOS and GINS currently contribute only with a subset of forces, whereas all others accomplished the comparison for all accelerations provided in the reference data set. The limit for the difference in evaluating the models along the reference orbit was set to 10⁻¹¹ m s⁻², thus, at least one order of magnitude lower than the accelerometer noise in the high-precision axes of GRACE (Touboul et al., 1999). As an absolute threshold, this limit does not take into account that the different forces influence orbit and gravity field solutions in a different way. Thus, for instance, a relatively large difference does not necessarily map to a final solution. Furthermore, we formulate this data set and a threshold of 10⁻¹¹ m s⁻² as a benchmark for new analysis centres which are interested to contribute to the COST-G initiative. The performance is shown in Fig. 3, the dashed black line marks the threshold. A result below that line is considered to agree with the benchmark data sufficiently well. All COST-G ACs and candidate ACs fulfill the requirement of 10⁻¹¹ m s⁻² for the tested accelerations. Consequently, we consider systematic errors related to the application of background forces to be eliminated to the extent possible when compared to the GRACE observation precision. The solid Earth tides, due to its complexity, turned out to be the most challenging acceleration of the benchmark data set. Every software in the comparison delivers differences slightly above the numerical precision.

5 Conclusions

We test and publish benchmark data of forces commonly used for gravity field and orbit determination purposes. The data set consists of orbital positions and accelerations evaluated along a one-day GRACE orbital arc. It is intended to enable fundamental software comparisons and bug detection and shall serve as a valuable contribution to the community. It is potentially interesting to any user of space geodetic software but especially to those groups where a solution is created from individual software packages by combination. The benchmark data was examined with the software packages currently used in the COST-G service. These packages agree with each other in the usage of the background models at a level of less than 10⁻¹¹ m s⁻². Further influences, such as the use of a set of different Earth rotation parameters, on the accelerations are scheduled for further investigations. The data set will be used, among other criteria (see https://cost-g.org, last access: 1 December 2020), as a benchmark for groups who are interested to join the COST-G initiative as an analysis centre.