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

. In the framework of the COmbination Service for Time-variable Gravity ﬁelds (COST-G) gravity ﬁeld 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 cru-cial role in gravity ﬁeld determination. Therefore, we publish a set of data of various commonly used forces in orbit and gravity ﬁeld modelling (Earth’s gravity ﬁeld, 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


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 op-erationally 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.

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.
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 1 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.

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 (1) 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 twiceper-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 d max = max(|r − s|). (2)

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.

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 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 with ∇ being the gradient operator. The accelerations in cartesian coordinates can be obtained by a transformation of Earth-fixed (r, λ, ϑ)-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.

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, 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.

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): 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 Earthfixed frame.
Step 2: Corrections on c 20 -corresponds to Eq. IERS (6.8a) and uses Table IERS 6.5b: Corrections on c 21 , c 22 , s 21 , s 22 -corresponds to Eq. IERS (6.8b) and uses Table IERS 6.5a and Table IERS 6.5c: The amplitudes for the frequency dependent corrections for A ip f and A op f are listed in Tables IERS 6.5a, 6.5b and 6.5c. The Doodson angle argument reads as Finally, the Doodson arguments read as 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 n om1 = [0, 5, 5, 5, 6, 5]. The first element of the vector tells the periodicity of the tide.
The final coefficients c tot 2m and s tot 2m are obtained by adding together step 1 and step 2 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.

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 cos nm,f , s cos nm,f ) and retrograde (c sin nm,f , s sin nm,f ) coefficients in a sum over all tidal frequencies f using 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 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 −11 m s −2 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).
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 0, 0, 9.8 × 10 8 T m 2 s −1 (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 a relativistic = a Schwarzschild + a Lense Thirring + a de Sitter .
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 1 and t 2 is performed on the level of spherical harmonic coefficients to obtain the spherical harmonic coefficients at time

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 21 and s 21 are concerned and the contribution is computed from the polar motion parameters and the mean pole definition via x P , y P ,x P andȳ P form the wobble parameters (Eq. IERS (7.24)) via The pole tide and ocean pole tide computed in the benchmark data make use of the most recent secular pole (IERS, 2018) x P = 55.0 × 10 −3 + 1.677 × 10 −3 t, y P = 320.5 × 10 −3 + 3.460 × 10 −3 t, The relation between the spherical harmonic coefficients and the accelerations is again given by Eqs. 3 and 4.

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 where the coefficients c cos nm,f , c sin nm,f , s cos nm,f , s sin nm,f 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).

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).
where A Re nm , A Im nm , B Re nm , B Im nm are the coefficients from the model, m 1 and m 2 are the wobble parameters (Eq. 25), γ Re 2 = 0.6870, γ Im 2 = 0.0036 and the factor R n is given by 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 the load Love numbers. To obtain accelerations, the spherical harmonic coefficients are evaluated using Eqs. (3) and (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.

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.

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.

EPOS
The Earth Parameter and Orbit System (EPOS 2 ) 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).

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 8 M. Lasser et al.: Benchmark data for verifying background model implementations 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.

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).

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.

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 , 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  and operational processing of GRACE-FO data.

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).

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 −11 m s −2 , 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 −11 m s −2 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 −11 m s −2 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.

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 −11 m s −2 . 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.
Data availability. The data is available via ftp from ftp://ftp.tugraz.at/outgoing/ITSG/COST-G/softwareComparison/ (Mayer-Gürr and . It is freely accessible. A readme-file (00README_simulation.txt) guides the user through the dataset.
Author contributions. TMG and AK compiled and provide the dataset, ML, KHN, JML, IK and MW performed the computational work of the comparison with the respective software package. Concept, framework and interpretation of the study was accomplished by all authors.
Competing interests. The authors declare that they have no conflict of interest.
Special issue statement. This article is part of the special issue "European Geosciences Union General Assembly 2020, EGU Geodesy Division". It is a result of the EGU General Assembly 2020, 4-8 May 2020.