Topical Review
Tropospheric and ionospheric media calibrations based on global navigation satellite system observation data
^{1}
Telespazio VEGA Deutschland GmbH c/o European Space Operations Centre (ESA/ESOC), Navigation Support Office,
RobertBoschStraße 5,
64293
Darmstadt, Germany
^{2}
DEIMOS Space c/o European Space Operations Centre (ESA/ESOC), Flight Dynamics,
RobertBoschStraße 5,
64293
Darmstadt, Germany
^{3}
PosiTim UG c/o European Space Operations Centre (ESA/ESOC), Navigation Support Office,
RobertBoschStraße 5,
64293
Darmstadt, Germany
^{4}
VisionSpace Technologies c/o European Space Operations Centre (ESA/ESOC), Navigation Support Office,
RobertBoschStraße 5,
64293
Darmstadt, Germany
^{5}
European Space Operations Centre (ESA/ESOC), Navigation Support Office,
RobertBoschStraße 5,
64293
Darmstadt, Germany
^{6}
European Space Operations Centre (ESA/ESOC), Flight Dynamics,
RobertBoschStraße 5,
64293
Darmstadt, Germany
^{*} Corresponding author: Joachim.Feltens@esa.int
Received:
3
July
2017
Accepted:
6
March
2018
Context: Calibration of radiometric tracking data for effects in the Earth atmosphere is a crucial element in the field of deepspace orbit determination (OD). The troposphere can induce propagation delays in the order of several meters, the ionosphere up to the meter level for Xband signals and up to tens of meters, in extreme cases, for Lband ones. The use of media calibrations based on Global Navigation Satellite Systems (GNSS) measurement data can improve the accuracy of the radiometric observations modelling and, as a consequence, the quality of orbit determination solutions.
Aims: ESOC Flight Dynamics employs ranging, Doppler and deltaDOR (DeltaDifferential OneWay Ranging) data for the orbit determination of interplanetary spacecraft. Currently, the media calibrations for troposphere and ionosphere are either computed based on empirical models or, under mission specific agreements, provided by external parties such as the Jet Propulsion Laboratory (JPL) in Pasadena, California. In order to become independent from external models and sources, decision fell to establish a new inhouse internal service to create these media calibrations based on GNSS measurements recorded at the ESA tracking sites and processed inhouse by the ESOC Navigation Support Office with comparable accuracy and quality.
Methods: For its concept, the new service was designed to be as much as possible depending on own data and resources and as less as possible depending on external models and data. Dedicated robust and simple algorithms, well suited for operational use, were worked out for that task. This paper describes the approach built up to realize this new inhouse internal media calibration service.
Results: Test results collected during three months of running the new media calibrations in quasioperational mode indicate that GNSSbased tropospheric corrections can remove systematic signatures from the Doppler observations and biases from the range ones. For the ionosphere, a direct way of verification was not possible due to nonavailability of independent third party data for comparison. Nevertheless, the tests for ionospheric corrections showed also slight improvements in the tracking data modelling, but not to an extent as seen for the tropospheric corrections.
Conclusions: The validation results confirmed that the new approach meets the requirements upon accuracy and operational use for the tropospheric part, while some improvement is still ongoing for the ionospheric one. Based on these test results, green light was given to put the new inhouse service for media calibrations into full operational mode in April 2017.
Key words: Troposphere / ionosphere / media calibrations / fitting/interpolation techniques / vector methods / CSPformat / orbit determination / radiometric tracking
© J. Feltens et al., Published by EDP Sciences 2018
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Atmospheric delays are significant contributors to the error budget of satellite tracking. The total atmospheric effect on satellite signals comprises two major parts: 1) the contribution from the troposphere and 2) the contribution from the ionosphere (e.g. Seeber, 2003). The troposphere is the neutral lower part of the Earth atmosphere extending from ground up to 50 km (stratosphere included). For frequencies below 30 GHz, the troposphere is a nondispersive medium, i.e. the tropospheric range delays do not depend on satellite signal frequency. For mathematical modelling, the troposphere can be described as a mixture of two ideal gases: dry air and water vapour, whereby the dry part contributes about 90% of the total tropospheric refraction and the wet part the remaining 10%. Typical total tropospheric range delays are between 2.3 m in zenith direction and up to almost 30 m at low elevations (e.g. Seeber, 2003). Above the troposphere resides the ionosphere, extending from 50 km up to about 2000 km, where it transits into the plasmasphere, which smoothly transfers into the interplanetary plasma at several 10000 km altitude (e.g. Davies, 1990). There exist no sharp transitions between ionosphere, plasmasphere and interplanetary plasma. Ratcliffe (1972) uses the following formulation to describe nicely the ionosphere as “The part of the atmosphere above about 60 km, where free electrons exist in numbers sufficient to influence the travel of radio waves …”. The ionization in the ionosphere is primarily triggered by solar radiation. In addition, the composition of atmospheric gases at different altitudes, the geomagnetic field, transport and recombination processes play an important role to establish the typical ionospheric structures in terms of different layers, equatorial anomaly, auroral oval, etc.; the interested reader is referred at this point to the pertinent literature (e.g. Davies, 1990; Ratcliffe, 1972; Komjathy, 1997; Kelley, 1989). In contrast to the troposphere, the ionosphere is a dispersive medium, i.e. the delay a radio signal suffers by the ionosphere depends on the signal's frequency. Typical values of ionospheric range delays (for singlefrequency measurements) range between 10 cm and up to 500 m. Ionospheric range delays are inversely proportional to the square of the satellite carrier frequency and depend on elevation (e.g. Seeber, 2003).
ESOC Flight Dynamics has to apply tropospheric and ionospheric corrections to the tracking data of ESA spacecraft (range, Doppler, deltaDOR) to achieve accurate orbit determinations (OD), especially for the interplanetary space probes. In addition for deltaDOR (http://www.esa.int/Our_Activities/Operations/About_delta_DOR), these media calibrations are also needed for quasar tracking. Concerning frequency bands used, there is currently a change ongoing from SBand tracking to XBand and KaBand tracking; SBand is actually no longer used at ESA for tracking deepspace missions. Thus the troposphere is considered to be the dominant term while the ionosphere becomes less significant.
So far, ESOC Flight Dynamics has computed its tropospheric and ionospheric calibrations inhouse based on empirical background models (Saastamoinen, 1973, Klobuchar, 1975, Jakowski et al., 2011) or, under specific agreements, received support from the Jet Propulsion Laboratory (JPL) in Pasadena, California, U.S.A, in the form of timenormalized polynomials. These polynomials will be described in detail in Section 3. On the other hand, the ESOC Navigation Support Office generates as part of its routine Global Navigation Satellite Systems (GNSS) data processing for the International GNSS Service (IGS, http://www.igs.org/) GNSS orbit and clocks data, tropospheric zenith delays and ionospheric slant Total Electron Content (sTEC) observables, among others, on regular basis (http://navigationoffice.esa.int/about/activities/). In addition, the ESOC Navigation Support Office operates and maintains its ESA's GNSS Observation Network (EGON) permanently recording multifrequency data (http://navigationoffice.esa.int/ESA's_GNSS_Observation_Network_(EGON).html), and some of these stations are also equipped with WEather STation (WEST) meteorological data recorders. These WEST data instruments belong to the ESOC Ground Systems Engineering Department. − So an inhouse generation of media calibrations by using this data is an obvious task, with the target to be as independent as possible from external data and models.
The ground stations to be supported are primarily the three ESA Deep Space Network (DSN) sites Cebreros/Spain (CEBR), Malargüe/Argentina (MGUE) and New Norcia/Australia (NNOR) (http://www.esa.int/Our_Activities/Operations/Estrack). Further sites may be included in future.
The rest of this paper is structured as follows: Section 2 points out general aspects of the algorithms design and concepts. Section 3 explains how the tropospheric calibrations are computed. Section 4 explains how the ionospheric calibrations are computed. Section 5 presents results obtained with the newly established media calibration service provided by the Navigation Support Office to Flight Dynamics. Section 6 concludes the paper by summarizing the new media calibrations service implementation at ESOC. Finally Appendix A provides some supplementing algorithms for ionospheric intersection point (IPP) position vector computation, which is needed for the ionospheric media calibrations computation.
2 Principles for new approach to compute media calibrations
Three key points were important when working out the concept for this new inhouse media calibrations service:

independency from external (especially nonESA) sources in terms of data, models and software as far as possible. − Nevertheless, global GNSS data collected with the IGS network is needed to get highprecision orbits and clocks, which are required to obtain reliable tropospheric zenith delay solutions, Section 3;

well suited for operational needs: Algorithms are simple and robust, optimized for efficient coding;

consideration of alternatives to retain the service in cases when ESA observation data (GNSS, meteo data from the ground sites) recording is interrupted.
With regard to the above three key points, the algorithms to compute the tropospheric calibrations, Section 3, and the ionospheric calibrations, Section 4, were designed in a way that they rely on ESA observation data as far as possible and on models as less as possible. It is assumed that an empirical approach based on observation data better reflects, apart from measurement errors and noise, the actual status of the troposphere/ionosphere rather than a model typically describing some mean or smoothed troposphere/ionosphere. Therefore the usage of models is foreseen in the algorithms presented in Sections 3 and 4 only in cases of nonavailability of observation data.
3 Tropospheric media calibrations
3.1 Tropospheric calibrations provided by JPL
Historically, tropospheric media calibrations processing at ESOC Flight Dynamics is aligned to JPL data formats, which may change in future. So the essentials of tropospheric calibrations provided by JPL (JPL, 2008) are summarized in the following for a better understanding of the own new approach, Section 3.2:

assumption of azimuthal symmetry, i.e. no gradients;

neither spacecraft nor frequencyspecific;

tropospheric calibrations are represented per station and in 6^{h} time intervals by one timenormalized polynomial for the wet (ZWD) − and one timenormalized polynomial for the dry zenith delay (ZHD), Eq. (2);

JPL employs the Niell Mapping Function (NMF) (Niell, 1996, 2000) to compute slant range tropospheric corrections of these ZWDs and ZHDs. − The NMF is a function of season, station latitude and height, but not of surface meteorological data;

the 1^{st} order (or background) tropospheric calibration is computed with socalled seasonal models for the ZWD and ZHD component (JPL, 2008; Chao, 1973) which do not depend on realtime data. The seasonal model parameters are delivered to the user in form of a timenormalized Fourier series representation, separately for the ZWD and ZHD:
A_{i}, B_{i} – fitted Fourier coefficients;
N – degree of Fourier series expansion;
T – calibration time [sec];
S – start time [sec];
P – period of fundamental mode [sec].
The seasonal models, as well as the JPL media calibrations in general, are delivered to the users in the Control Statement Processor (CSP) format (JPL, 2008);

on top of these ZWD and ZHD background calibrations, time series of deltacorrections of tropospheric zenith delays, wet and dry, are derived from nearrealtime weather and GPS data and delivered to the users in form of normalizedtime polynomial representations:
C_{i} – fitted polynomial coefficients;
N – degree of polynomial representation;
T – calibration time [sec];
S – start time [sec];
E – end time [sec].
The Eq. (2) polynomial representation is used to compute both the wet and dry deltacorrections. The user has then to calculate the actual tropospheric zenith delay, separately for wet and dry, as sum of the value obtained from the seasonal model (background) plus the deltacorrection computed with the actual polynomial.
3.2 Computation of tropospheric calibrations by the ESOC Navigation Support Office
Tropospheric media calibrations are independent of spacecraft and, below 30 GHz, also of frequency. Thus per day, one ACSII file is produced in CSPformat by the ESOC Navigation Support Office, listing for ESA's three DSN sites CEBR, MGUE, NNOR polynomial representations of wet and dry deltazenith delays in 6 hour time intervals. Such a CSPfile is generated once per day for the day before in rapid mode. In final mode the same process is repeated with 4 days latency and more GNSS observation data. Each polynomial representation, Eq. (2), describes then wet and dry deltazenith delays w.r.t. the background ZWD and ZHDs computed with the seasonal model.
Before entering into the polynomial fit, Eqs. (3)–(5), the tropospheric zenith delay values are preprocessed as follows:

a SINEX.tro file (ftp://igs.org/pub/data/format/sinex_tropo.txt) obtained from daily routine GNSS processing at the Navigation Support Office provides 48^{h} time series of total tropospheric zenith delays (ZTDs) and ZWDs for all stations considered (Fig. 1);

in addition, at the ESA DSN sites CEBR, MGUE, NNOR, WEather STation (WEST) data is permanently recorded and enters into the media calibrations process as 1 minute samples. The WEST data comprise, among others, temperature, pressure and humidity;

detection and removal of large outliers from the tropospheric zenith delay series;

the Saastamoinen (1973) Model is employed to compute the ZHD component of the ZTD with the WEST data. In the case that pressure data from the WEST stations is not available, corresponding weather parameters are computed with the Global Pressure and Temperature model (GPT) (Böhm, 2007; Steigenberger et al., 2009; Lagler et al., 2013);

the ZHDs are corrected for the height difference between the GNSS/IGS antenna on DSN site, w.r.t. which the GNSS measurements and parameter estimation were made, and the reference point of the DSN antenna (point of antenna azimuth axis with elevation axis crossing), using a formula of Estefan and Sovers (1994);

the ZTDs and ZWDs provided in the SINEX.tro file were determined by using temperature and pressure values from the GPT model; this is basically a setup for the routine IGS processing. On the other hand, Flight Dynamics prefers, due to its systems setup, to get its tropospheric calibrations referred to the Saastamoinen Model, with the dry part computed from WEST data pressure measurements. Thus the ZHDs computed with the Saastamoinen Model from WEST data, are subtracted from the ZTDs of the SINEX.tro file, leading to ZWDs referred to the Saastamoinen Model. − The ZWDs provided directly by the SINEX.tro file are in this way ignored for tropospheric media calibrations processing;

the JPL wet and dry seasonal model is subtracted from the above ZWD and ZHD series, resulting in corresponding deltaZWD and deltaZHD (δZWD and δZHD) time series into which polynomial expressions, Eq. (2), are then fitted, Eqs. (3)–(5).
For the tropospheric media calibrations service, the Navigation Support Office produces once per day a SINEX.tro file containing 48^{h} time series of ZWDs and ZTDs for all DSN sites considered. Processing is done in rapid and in final mode. In rapid mode, this computer run is submitted via an automatic job scheduler (crontab) every day at 14:30 UT and covers the time interval from 12:00 UT of today (observation day) 48^{h} backwards to 12:00 UT of the day before yesterday. In a subsequent computer run, this SINEX.tro file enters as input into the polynomials fitting and CSPfile generation. In final mode, the same process is repeated at 02:02 UT, 4 days after the end of the observation day, with more observation data.
Figure 1 shows how these computer runs of 48^{h} SINEX data generation are scheduled for a sequence of several days i, i + 1, i + 2, i + 3, … Also indicated are the rapid mode launch times at 14:30 UT. The 48^{h} run of each day is processed independently of the other days. Of this 48^{h} series of ZWDs and ZTDs, δZWD and δZHD series are then established as described above. Into the δZWD and δZHD series thus obtained, separately for wet and dry, a sequence of eight 6^{h} polynomial representations is fitted. In order to avoid jumps at the interfaces between one polynomial representation and the next, the eight 6^{h} polynomials sequence is fitted to the complete 48^{h} deltazenith delays series within one fitting process, as indicated in Eq. (5). In order to suppress jumps at the interfaces, constraints are introduced between successive polynomial representations in terms of offset, slope and slope rate, depending on the user's setting. Constraints are introduced between two successive polynomial representations in terms of conditional Eqs. (3.3a–c): (3a) (3b) (3c) where – function, 1^{st} and 2^{nd} derivative of preceding polynomial; – function, 1^{st} and 2^{nd} derivative of subsequent polynomial.
The polynomial arguments are computed according to Eq. (2): (4) where
X – argument of preceding polynomial;
Υ – argument of subsequent polynomial;
S_{P}, E_{P} – start/end time of preceding polynomial;
S_{Q}, E_{Q} – start/end time of subsequent polynomial.
The interface epoch between the two successive polynomial representations is defined by the end time E_{P} of the preceding polynomial P(X), which is equal to the start time S_{Q} of the subsequent polynomial Q(Y), i.e. E_{P} = S_{Q}. Putting these values for T into Eq. (4) leads to X_{E} = +1 and Y_{S} = −1 as arguments for the constraining Eq. (3a–c):
There is an agreement with Flight Dynamics that the polynomial representation for the tropospheric media calibrations should be of polynomial degree Four. This choice was a consequence of a trial and error process, during which it was observed that in some cases high degree polynomials were following the noise of the zenith delay estimates, which are uncorrelated and not smoothed, thus creating artificial slopes in the calibrations that were worsening the results.
Equation (5) displays the observation equations structure of two successive polynomial representations P(X) and Q(Y) with the interface constraints between them. These constraining equations are multiplied with weights W_{offs}, W_{slpe}, W_{slrt}, which can be set by the user. Optimal values for these weights were found out by trying: W_{offs} = 100, W_{slpe} = 100, W_{slrt} = 0. An introduction of a Zeroweight leads to a nonconstraining of the respective parameter, i.e. the slope rate is with the actual setting not constrained. But principally the slope rate could be constrained too and is thus marked with “eventually” in Eq. (5).
In the ideal case, the interface epochs of the eight polynomial representations coincide with measurement epochs of the tropospheric zenith delays. In this way, the last deltatropospheric zenith delay (δTZD, standing here for δZWD or δZHD) value of the preceding polynomial is also the first deltazenith delay value of the subsequent polynomial, i.e. this deltazenith delay value appears twice in the observation equations (δTZD_{m} in Eq. (5)). For the sake of simplicity and better readability, Eq. (5) is set up here for 3^{rd} degree polynomial representations only and only for a sequence of two successive polynomial representations. (5) Figure 2 displays the principal structure of the complete Amatrix for a 8 × 6^{h} polynomials fit, also for the ideal case of no gaps in the tropospheric zenith delay data and tropospheric zenith delay epochs coinciding with the interface epochs of the polynomial representation time spans. For simplicity, only three tropospheric zenith delay epochs are shown for each polynomial representation.
So, in summary, one Amatrix is set up and solved in a least squares fit to obtain the estimated coefficients for all polynomials.
Equation (5) and Figure 2 represent the ideal case in terms of full tropospheric zenith delay data availability and tropospheric zenith delay epochs coinciding with the interface epochs of the eight polynomial representations. However, deviations in tropospheric zenith delay values availability and distribution from these ideal conditions can happen in daily routine operational processing. Thus the above described approach had to be enabled to handle also such cases:

in the case that an interface epoch does not coincide with one of the tropospheric zenith delay epochs, the fit of the preceding polynomial takes all δTZD values falling with their epochs into its respective representation time interval and the fit of the subsequent polynomial takes all δTZD values falling with their epochs into its respective representation time interval, i.e. no δTZD value will, due to nonexistence, be introduced twice at the interface epoch. The constraining equations will, nevertheless, be established at the interface epoch;

if in a given polynomial representation time span the number of available δTZD values is due to gaps less or equal to the specified polynomial degree, the polynomial degree is for that time span set to One (linear fit) in the case that at least two δTZD observables are available and that the affected 6^{h} interval is not the first or the last interval of the 8 × 6^{h} polynomials fit. If only one observable is available and/or the affected interval is not the first or the last interval of the 8 × 6^{h} polynomials fit, the polynomial degree is set to Zero (constant fit). This processing scheme is retained as long as at least one observable falls into the respective time interval. This approach of setting the polynomial degree generally to One or Zero in case of sparse data is to avoid oscillations that may propagate into the neighbouring polynomials;

a constraining is only made between successive polynomial representations. In the case that due to a longer observation data gap there are no δTZD values available for one or more complete polynomial representations / 6^{h} time intervals at all, no constraining is introduced between the last polynomial before − and the first polynomial after this gap.
Figure 3 shows two examples of sequences of polynomial 8 × 6^{h} fitting in the case of gaps. Obviously the constraining scheme, Eqs. (3) and (5), in combination with linear or constant polynomial fitting over the gaps, is able to strongly mitigate the effects of these data gaps. The neighbouring polynomials enclosing the affected 6^{h} time interval remain almost unaltered.
Expressed in CSPformat, the four central 6^{h} polynomials (blue in Fig. 1) are finally written into the daily CSPfile being delivered to Flight Dynamics. At Flight Dynamics, the polynomials are evaluated and added onto the JPL Seasonal Model values to compute absolute ZWD and ZHD corrections for an ESA station. Then with the NMF, the computed tropospheric zenith delays are mapped into any direction for the correction of ESA spacecraft/quasar tracking data recorded at that station.
In the case of nonavailability of GNSS observation data and/or systems, Flight Dynamics will compute tropospheric corrections based on the Saastamoinen Model evaluated with meteo data recorded locally on the affected ground sites.
Fig. 1 Time schedule for the daily computer runs at the Navigation Support Office to establish 48^{h} time series of ZWDs and ZTDs, and tropospheric media calibrations. 
Fig. 2 Ideal case Amatrix structure. 
Fig. 3 Examples of polynomial fitting in the case of data gaps, marked in red. 
4 Ionospheric media calibrations
4.1 Ionospheric calibrations provided by JPL
Historically, ionospheric media calibrations processing at ESOC Flight Dynamics is aligned to JPL data formats, which may change in future. So the essentials of ionospheric calibrations provided by JPL (JPL, 2008) shall be summarized in the following for a better understanding of the own new approach, Section 4.2:

specific for each station/spacecraft/frequency. − Upon need, ionospheric calibrations can be rescaled to other frequencies very easily;

ionospheric calibrations are computed per station for each tracking pass (from satellite rise to satellite set);

per pass, the ionospheric calibrations are provided as lineofsight corrections w.r.t. a reference frequency of 2295 MHz, which must be applied to the satellite ranging data;

ionospheric calibrations are actually represented per pass in form of a polynomial as function of normalized time, Eq. (2). − Nevertheless, JPL has also the capability to provide multiple piecewisecontinuous polynomials per pass;

no background model, like the JPL Seasonal Model for the troposphere, is applied, i.e. the ionospheric polynomials provide absolute slant range delays;

the user has to provide satellite trajectory information, giving the direction to the satellite as function of time, as input for the ionospheric calibrations computation. For deltaDOR observations, the quasar direction as function of time is required too.
4.2 Computation of ionospheric calibrations by the ESOC Navigation Support Office
As with the tropospheric calibrations, also for the ionospheric calibrations an approach has been worked out, which is based on own GNSS tracking data recorded and evaluated by the Navigation Support Office. This approach is simple, robust and easy to handle, relying as far as possible on inhouse data and software and being as far as possible independent from external sources. − The only exception from this rule is the usage of the selfcontained Neustrelitz TEC Model − Global (NTCMGL) (Jakowski et al., 2011) for the computation of ionospheric background values.
The following input is required from the Navigation Support Office:

GNSS data (RINEX) from the sites of the ESA DSN as well as from the global IGS network, including its preprocessing;

GNSS satellite orbits and clocks from routine products generation for the IGS;

sTEC and Differential Code Bias (DCB) values from routine ionosphere products generation for the IGS.
The following input is required from Flight Dynamics:

ESA spacecraft/quasar trajectory information for the specified 48 hour time interval.
Interplanetary spacecraft/quasar passes never exceed 15 hours significantly. Thus per day an ionospheric media calibration is computed covering the 48hours timespan from 0^{h} UT of today (observation day) back to 0^{h} UT of the day before yesterday. As with the troposphere, also ionosphere processing is conducted in rapid and in final mode. In rapid mode, this job is submitted via an automatic job scheduler (crontab) every day at 07:33 UT. In order to avoid short, truncated passes at the beginning and at the end of this 48hours timespan, only those passes are considered which fall completely into the 48hours interval. All other passes are considered in the jobs of the day before or after. The pass criterion is for a spacecraft/quasar to have an elevation angle of ≥5° at the considered station. Per day and ESA spacecraft/quasar, one CSPfile is generated covering all complete passes of this spacecraft/quasar at all considered stations for the last 48 hours from 0^{h} UT to 0^{h} UT. In final mode, the same processing is repeated at 07:03 UT, 4 days after the end of the observation day.
The developed method is based on the evaluation of dualfrequency GNSS data which is permanently recorded with receivers owned by the Navigation Support Office at numerous ESA tracking sites, among them the DSN sites CEBR, MGUE and NNOR. After preprocessing, series of uncalibrated sTEC values at each ESA DSN to all visible GNSS satellites over the specified 48hours time interval are produced from the RINEX data. These sTEC values have to be corrected for the GNSS receiver and satellite DCBs (Schaer et al., 1998). The GNSS receiver and satellite DCBs are obtained from the global ionosphere model estimates conducted operationally each day at the Navigation Support Office, in which these DCBs are part of the estimation process. It should be noted that these DCBs cannot be estimated in an absolute sense so a certain reference (datum) has to be defined. This is done by specifying that the sum of the satellite DCBs is Zero. To avoid significant daytoday variations in the DCB estimates, e.g. due to a missing satellite on a specific day changing the average, the daily DCB estimates of the last 7 days are combined. This smooths out the daytoday variations and avoids jumps due to satellites missing on individual days, e.g. due to manoeuvres. It does, however, not resolve the change in the mean DCB level when a new satellite is added to the GNSS constellation.
Then the algorithm works as follows per DSN station and ESA spacecraft/quasar: Per epoch, the unit vectors pointing from the station to all GNSS satellites and to the ESA spacecraft/quasar are simply computed from their position vectors (Fig. 4): (6) where_{} – unit vector pointing from DSN site to spacecraft i;_{} – geocentric earthfixed position vector of spacecraft i (i = GNSS_{1}, GNSS_{2}, … GNSS_{n}, ESAspacecraft/quasar);
– geocentric position vector of DSN site.
Next, the cosines of the angles between all G_{j} = 1, …, n GNSS satellites and the ESA interplanetary spacecraft (ESA)/quasar (Q), as viewed from the DSN site, are computed: (7)
The smaller the angle between a GNSS satellite and the ESA spacecraft/quasar is, the closer its cosine value comes to +1, i.e. GNSS satellites viewed close to the ESA spacecraft/quasar will display cosine values near to +1, while far distant GNSS satellites will get cosine values close to −1, e.g. when both GNSS and ESA spacecraft/quasar are close to horizon but diametrical in direction (e.g. Fig. 6). Thus the are well suited to establish a weighting when combining the GNSS sTEC values to a sTEC value into the direction of the ESA spacecraft/quasar. The following weighting function was selected for the ionospheric media calibrations: (8)
i.e. for . With the power κ in Eq. (8) the distribution of weights might be further tuned. − Nevertheless, κ = 1 is used for the ionospheric calibrations.
By introducing these weights, a weighted mean of the GNSS sTEC values, depending on the angular distances of the GNSS satellites from the ESA spacecraft/quasar, could principally be computed as: (9)
Tests, however, showed that Eq. (9) does not provide satisfactory sTEC values into ESA spacecraft/quasar direction. The above principal approach, Eqs. (6)–(9), had thus to be improved and supplemented by the following measures in order to meet anticipated accuracy margins:
Fig. 4 Considered vector geometry between DSN station, GNSS satellites and ESA spacecraft/quasar. 
Fig. 6 Maximal possible geocentric angle ξ_{max} for given mean earth radius R_{M} and ionospheric shell height h_{I}_{.} 
4.2.1 Measure 1
An alternative to directly processing the sTEC values, is to map them with a mapping function, Eq. (10), into the vertical in order to compute the weighted mean vTEC and then to map this weighted mean vTEC back into the slant range direction of the ESA spacecraft/quasar: (10) whereζ = 90^{0} − el – zenith distance at ionospheric shell of maximum ionization (IPP, Appendix A). Standard Mapping Function, IONEX and NTCMGL; Modified Single Layer Model Mapping Function.
Depending on parameters (R_{M}, h_{I}, α) setting, Eq. (10) can either be evaluated as the Standard Mapping Function, or as the Modified Single Layer Model Mapping Function (MSLM, ftp://ftp.unibe.ch/aiub/users/schaer/igsiono/doc/mslm.pdf) (Schaer, 1999). The MSLM is an optimal fit to the extended slab model mapping function used by JPL (Mannucci et al., 1998; Coster et al., 1992).
The relation between sTEC and vTEC is then: (11)
With Eq. (11) the weighted mean of the vTEC can be expressed in terms of sTEC and mapping function, by modifying Eq. (9): (12)
By solving Eq. (12) for sTEC_{ESA/Q} and introducing Eq. (10) for the mapping functions, one obtains after the some algebra the following formulation: (13) where
C = h_{I} ⋅ (2 ⋅ R_{m} + h_{I}) – auxiliary constant.
The zenith distances ζ_{ESA/Q,Gj} required in Eq. (13) have to be computed at the IPP positions, where the slant ranges to the ESA spacecraft/quasar and to the GNSS satellites intersect the ionospheric shell of maximum ionization h_{I}. The computation of the IPP position vectors is described in Appendix A, whereby here the unit vector pointing into range direction is already known, Eq. (6), and must not be first computed from azimuth and elevation, i.e. the calculations can enter directly at Eq. (A.4) into the Appendix A algorithm. With these IPP position vectors, the required cosines of the zenith distances are evaluated: (14) where R_{m} – mean earth radius, for values see Eq. (10);h_{I} – ionospheric shell height, for values see Eq. (10).
4.2.2 Measure 2
When inspecting Figure 5 it becomes obvious that the angular distances between the ESA spacecraft/quasar and the GNSS satellites, as viewed from the station, can provide a rather distorted picture of the local ionosphere above, whose principal layered structure is indicated in Figure 5 by its shell of maximum ionization at altitude h_{I}. So appears GNSS_{2} seen from the station much closer to the ESA spacecraft/quasar than the angular distance between the IPP position vectors of both would justify. Thus a probably higher weight is appointed to GNSS_{2} as it in reality would deserve. In order to achieve a more objective weighting, the angular distances ξ between the IPP position vectors (green in Fig. 5) should be computed and enter into the weight computation, rather than the angular distances ψ as viewed from the station (red in Fig. 5). So Eq. (7) is modified to: (15)
Before the cosine values obtained with Eq. (15) are entered into Eq. (8) for weight computation, their value range should be increased. While the angular distances ψ, as viewed from the station (red in Figs. 5 and 6), can achieve values between 0° ≤ ψ ≤ 180° (180° when both GNSS satellite and ESA spacecraft/quasar appear at the tracking site close to horizon, but almost diametrical in direction), the angular distances ξ between the IPP position vectors (green in Figs. 5 and 6) can get maximum values of 0° ≤ ξ ≤ ξ_{max} ≈ 40° only (Fig. 6). This angular range limitation of ξ will also reduce the range of cosine values for weight computation by Eq. (8) considerably. In order to achieve a wider range of angular arguments for the cosine function and thus a better weighting, the ξ are artificially enlarged as follows: (16)
Instead of the cosine values obtained from Eq. (15), the following cosine values should then enter into Eq. (8) for weight computation, based on Eq. (16): (17) where (Fig. 6, green triangles).
Fig. 5 Computing weights from the angles enclosed by the IPP position vectors. 
4.2.3 Measure 3
The weighted mean sTEC value computed with Eq. (13) (considering Eqs. (15) and (17)) does represent the sTEC in the direction to the ESA spacecraft/quasar only with limited accuracy and is “corrected” with the aid of the NTCMGL background model as follows: (18a)
The right hand side of Eq. (18a) comprises the following parameters:

Weighted mean TEC ratio between actual level of ionization and NTCMGL background level:

TEC ratio between actual level of ionization and NTCMGL background level at the jth GNSS satellite:

… mapping function value, Eq. (10), evaluated at ESA spacecraft/quasar IPP location;

… mapping function value, Eq. (10), evaluated at jth GNSS IPP location;

… sTEC measured by jth GNSS satellite, appearing also in Eq. (9);

… vTEC corresponding to measured ;

… NTCMGL vTEC computed at jth GNSS IPP location;

… NTCMGL vTEC computed at ESA spacecraft/quasar IPP location;

… NTCMGL vTEC computed at weighted mean IPP location, Eq. (18d);

Unit vector pointing from ground station into the direction of weighted mean sTEC, Eq. (13), componentwise mean of GNSS range unit vectors:
The unit vectors pointing from ground station to the j GNSS satellites were computed with Eq. (6). is needed to compute the IPP location , Appendix A, corresponding to the weighted mean sTEC. At this IPP, a NTCMGL vTEC value has to be computed too for Eq. (18a).
In summary, NTCMGL vTEC values need to be computed for Eq. (18a) at the IPPs of the ESA spacecraft/quasar, all GNSS satellites and the weighted mean sTEC direction.
Equation (18a) can be summarized as follows: The weighted mean sTEC, as obtained from Eq. (13) (considering Eqs. (15) and (17)), is corrected for a sTEC offset between the IPP location of the ESA interplanetary spacecraft/quasar direction and the IPP corresponding to the direction of the weighted mean over all GNSS satellite directions involved. This offset is computed as NTCMGL vTEC difference between the ESA spacecraft/quasar IPP location and the IPP location of the weighted mean sTEC direction, mapped into the slant range at the ESA spacecraft/quasar IPP location. By the difference making, potential biases between the NTCMGL and measured GNSS sTEC observables cancel. To account for the actual ionospheric activity level, corresponding to the actual ionospheric overall charge w.r.t. the NTCMGL background level, at the IPP position of each GNSS satellite a background TEC value is computed with the NTCMGL too and set in ratio to the actually measured sTEC at that IPP, Eq. (18c). The vTEC difference, ESA spacecraft/quasar IPP minus weighted mean sTEC IPP, is then scaled by the weighted mean ratio, Eq. (18b), to account for the actual ionization level w.r.t. the ionization level represented by the NTCMGL. As closer as the IPP positions into ESA interplanetary spacecraft/quasar direction and of the weighted mean sTEC coincide, as closer this difference will approach Zero, and finally vanish in case of coincidence. In addition, by applying this correction, also the situation of unfavourable geometrical conditions will be mitigated, e.g. in case that an ESA spacecraft/quasar is not surrounded by directions to GNSS satellites. − Normally this should not happen, since interplanetary missions are typically located in − or close to the ecliptic plane, a region which is well covered by the GNSS constellations. Measure 3 included _ see Figure 11a
Figure 7 shall illustrate how the weighted mean sTEC values (red & violet curves) approach closer and closer the sTEC curve of an excluded GNSS satellite (blue), which served as reference, by introducing in sequence Measure 1–3, whereby the final curve (Measure 3) is shown in Figure 11a. A description of the figure details and all curves displayed and their meaning is also provided beneath Figure 11.
Theoretically, Eqs. (6)–(18) could also work with only one GNSS satellite. Of course the algorithm will in practice not be applied in this way, but this indicates that it is very robust and stable, and thus well suited for operational use. As an indicator of quality for the interpolated sTEC value into ESA interplanetary spacecraft/quasar direction, the following additional parameter can be evaluated, in parallel to the weighted mean: (19)
Ξ is something like a “mean” of the weights of all GNSS satellites, based on Eq. (8). Ξ ≈ 1 would indicate a close angular distance of most of the GNSS satellites to the ESA spacecraft/quasar, Ξ ≈ 0.5 medium angular distance and Ξ ≈ 0 far angular distance, i.e. the closer Ξ comes to One, the better the sTEC value for the ESA spacecraft/quasar could be determined.
The angle between ESA spacecraft/quasar and weighted mean sTEC direction, either in terms of range from station (Eq. (20) top) or in terms of IPPs (Eq. (20) bottom), can serve as additional accuracy measure: (20)
The angle enclosed by ESA spacecraft/quasar direction and sTEC weighted mean direction should be preferably close to Zero and its cosine value should approach One.
In the case of nonavailability of GNSS observation data and/or systems, Flight Dynamics will compute ionospheric corrections with the NTCMGL.
For the ionospheric CSPfiles writing, the computed sTEC values along an ESA interplanetary spacecraft/quasar pass are then expressed in terms of group delays (e.g. Davies, 1990; Komjathy, 1997) w.r.t. to a reference frequency of 2295 MHz. This frequency was selected to be consistent with the current processing system set up at Flight Dynamics. (21) whereg_{ion} – ionospheric code group delay [m];40.3 – constant originating from the AppletonHartree refractive index formula [m^{3}/s^{2}];f – satellite signal frequency [Hz], here set to 2295 × 10^{6} Hz (rescaling to group delays at another frequency f by multiplication with {2295 × 10^{6} Hz/f}^{2});sTEC – total electron content along the signal's slant range path through the ionosphere [m^{−2}];
Finally, a normalizedtime polynomial of degree Five, Eq. (2), is fitted to this sTEC series of ionospheric slant range delays. One polynomial is fitted per pass. Per day and ESA spacecraft/quasar, one CSPfile is produced covering all complete passes of this spacecraft/quasar at all considered stations for the last 48 hours from 0^{h} UT to 0^{h} UT.
An implementation of media calibration extrapolations has so far not been considered and is not planned for the time being.
Fig. 11 (a) Example from Test Case 1; (b) example from Test Case 1 with sTEC simulated by NTCMGL; (c) example from Test Case 2; (d) example from Test Case 3. 
Fig. 7 (a) Computing sTEC_{ESA/Q} only with Eqs. (6)–(9); (b) Measure 1 included; (c) Measure 2 included; (d) Measure 3 included → see Figure 11a. 
5 Results
5.1 Tropospheric media calibrations
The tropospheric media calibrations were validated by analysing orbit determination statistics obtained when using the media calibrations derived from the Saastamoinen (1973) Model on one hand and when using the media calibrations computed with the approach described in Section 3 on the other. The Saastamoinen Model employed for the comparison makes use of real weather data collected at the stations, and includes both the wet and dry components of the tropospheric zenith delay. The Niell Mapping Function (Niell, 1996, 2000) is used in both cases to map the zenith delays to the desired elevation.
Figure 8 shows an example of a fitted 8 × 6^{h} polynomial sequence for the dry part (Fig. 8a) and for the wet part (Fig. 8b). The station in Fig. 8 exemplary case is MGUE. The blue noisy curves display the δTZDs derived from GNSS observables, the red curves the sequence of fitted polynomials. Only the four central polynomials of the 8 × 6^{h} sequence enter then into the CSP file of the respective day (blue in Fig. 1). Table 1 summarizes the RMS values obtained on the example day of Figure 8 (12 April 2016, 12:00 – 14 April 2016, 12:00; 16.104–16.105 in the plot headers indicate the time range for the central 24^{h} in YY.DDD) for all three ESA DSN sites. The RMS values were computed from the differences of GNSSderived δTZD minus computed polynomial value. A 15^{m} time resolution in the SINEX.tro file lead to 97 tropospheric zenith delay samples in the considered 48^{h} time interval.
In the following, the validation results for the new corrections used within the orbit determination of deepspace probes are presented. Statistics of orbit determinations (OD) made by using the tropospheric media calibrations computed with the approach described in Section 3 (called TRCal in the following) were compared with corresponding statistics when the tropospheric media calibrations from the legacy Saastamoinen Model were employed. For all the test cases presented for the validation of tropospheric corrections, the ionospheric delay is treated by the means of the NTCMGL model.
Two test scenarios, based on GAIA and LISA pathfinder routine orbit determination, have been used. Both satellites are orbiting in a very quiet dynamic environment, better suited for observation models comparison compared to planetary orbiters. The configuration of the OD is briefly summarised in the following paragraphs. For more details on deep space OD terminology and ESOC OD software refer to Moyer (1971, 2000) and Budnik et al. (2004).

GAIA: The spacecraft is orbiting in proximity of the SunEarth libration point L2, thus is mainly tracked during local evening, night or early morning. GAIA is extensively tracked from all of the ESA DSN stations, providing a very good test case for atmospheric delays validation. A data arc of three weeks duration (07 – 28 March 2017) has been selected for the orbit determination. Up to three Doppler passes per day are available from all ESA deep space stations, while, due to mission constraints, only a single ranging session per day is possible, split into two 10minutes periods at the beginning and the end of all NNOR passes. No station keeping manoeuvres were present in the chosen period, however some dynamic parameters needed to be estimated such as the corrections to the Solar Radiation Pressure (SRP) model and to the propulsion model, given that GAIA attitude is controlled with lowthrust propulsion. GAIA is always tracked at elevations higher than 10°. The Doppler measurements are compressed to oneminute count intervals, and the frequency of range measurements, when available, is one data point per minute. More details about GAIA orbit determination can be found in Budnik et al. (2012);

LISA Pathfinder: The satellite is orbiting in proximity of the SunEarth libration point L1, and is tracked by a single ESA station at Cebreros, mostly during local morning and daytime for a total duration of 8 hours per day, during which both Doppler and range measurements are continuously acquired. For this scenario a data arc of 18 days duration was used (13 Feb. − 03 Mar. 2017). LISA pathfinder is tracked with elevation down to 6°, which poses a major challenge in terms of tropospheric corrections. Due to this reason, for every type of calibration employed two runs were made, the first including all data and the second including only data at elevations higher than 10°. Also for this test case, the Doppler measurements are compressed with oneminute count intervals, while range measurements are sampled every 10 minutes. The dynamic parameters estimated in the orbit determination include only the satellite state vector and the corrections for the internal gravity between the satellite and the test masses used to control the forcefree motion of the satellite (Armano et al., 2016).
For both scenarios, the orbit determination has been performed with the ESOC operational orbit determination software AMFIN (Budnik et al., 2004). A range bias per pass is estimated in the fit to compensate for time dependent unmodelled effects. The aposteriori dispersion of the range biases is usually a good indicator of possible systematic errors in both the dynamics and the observations modelling, and of possible radial errors in the reconstructed trajectory.
The reported statistics include the overall RMS of all Doppler and range residuals, intended as the measured minus computed values of the observations based on the estimated solution (trajectory plus estimated parameters). The dispersion, in terms of RMS, of the estimated range biases is also included.
Figure 9 gives an overview of the evolution of the Doppler residuals over the entire data arc used for the GAIA OD, in this case with tropospheric corrections derived from TRCal. In Table 2 the statistics of the observation residuals are compared between the two calibration methods. It can be seen that all the figures present a noticeable improvement when passing from the Saastamoinenbased corrections to the TRCal ones. Even though the Doppler residuals RMS is dominated by the presence of some very noisy passes, likely due to high solar/geomagnetic activity, the benefit of GNSSbased corrections can be appreciated especially when the satellite is tracked at low elevation. In Figure 10, some tails at the beginning and the end of the passes, i.e. at low elevation, are visible in the Doppler residuals when the Saastamoinen Model is used (top plot), while they disappear when the TRCal are used (bottom plot). These tails are typical for a systematic error in the tropospheric zenith delays over relatively long time periods (tens of minutes to hours) mapped down to low elevation, and they can have a larger impact on the OD accuracy than a plain uncoloured noise of even larger amplitude like the one visible in some passes in Figure 9.
Table 3 reports the comparison of LISA Pathfinder orbit determination statistics. Also in this case the statistics on the range biases and Doppler residuals improve when the new GNSS based calibrations are employed. In particular, the Doppler residuals show a significant improvement (more than 10%) when the data at low elevation is included in the OD. The RMS of the range residuals does not show a noticeable improvement when changing the calibration method. This could be a consequence of the larger background noise, compared to the GAIA case, possibly masking small variations in the residuals.
On top of the formal validation of the operational tropospheric calibrations presented here, the corrections generated in preoperational mode were used successfully for over one year for GAIA routine orbit determination.
Fig. 8 (a) Example of drypart polynomials; (b) example of wetpart polynomials. 
RMS summary on tropospheric polynomial fitting, 12 April 2016, 12:00–14 April 2016, 12:00.
Fig. 9 Overview of Doppler residuals over the full arc used in the GAIA OD; tropospheric corrections based on TRCal. 
GAIA OD statistics comparison for tropospheric corrections.
Fig. 10 Comparison of 2way Doppler residuals with Saastamoinen (top plot) and TRCal (bottom plot) corrections. 
LISA Pathfinder OD statistics comparison for tropospheric corrections.
5.2 Ionospheric media calibrations
The ionospheric media calibrations were validated by two basic methods:

The ESA spacecraft/quasar was replaced by a GNSS satellite, whose sTEC values are known from its dualfrequency measurements. These measured sTEC values were then compared with the weighted mean sTEC derived from the dualfrequency data of all the other GNSS satellites, Section 4. In addition, corresponding sTEC values were interpolated from JPL IONEX files along the excluded GNSS satellite's trajectory, as well as computed with the NTCMGL;

An analysis of orbit determination statistics obtained by using media calibrations based on the legacy NTCMGL model on one hand and by using the media calibrations computed with the new approach, as described in Section 4, on the other.
5.2.1 Validation method 1
Figure 11 stands exemplarily for the results of a systematic validation of the ionospheric calibrations, that were conducted by replacing the ESA spacecraft/quasar with a GNSS satellite excluded from the polynomials computation and whose measured sTEC values served then as reference for the weighted mean sTEC values computed from the measurement data of all the other GNSS satellites, Eqs. (6)–(18). Three test cases were considered:

Test Case 1: 19 July 2015, 00:00 − 21 July 2015, 00:00, normal ionospheric conditions;

Test Case 2: 8 May 2016, 00:00 − 10 May 2016, 00:00, strong magnetic storm with Kp values up to 7 (the planetary Kindex, Kp, is a measure for the irregular variations of the Cartesian components of the geomagnetic field and is used to characterize geomagnetic storms, e.g. Komjathy (1997) and http://www.swpc.noaa.gov/products/planetarykindex);

Test Case 3: 8 March 2016, 00:00 − 10 March 2016, 00:00, solar eclipse over the Pacific (http://www.timeanddate.com/eclipse/solar/2016march9) which was at least partly visible from NNOR.
Per test case, five GNSS satellites, GPS as well as GLONASS, were arbitrarily chosen to serve in sequence as reference for the validation. For Test Case 3, GNSS satellites that were affected by the solar eclipse had been selected. Then per excluded GNSS satellite and pass/arc of that satellite, a CSPfile was computed from the data of all the other GNSS satellites and the weighted mean sTEC values were then compared with the sTEC series measured by the excluded GNSS satellite. Per test case, about 70 passes/arcs at the three ESA DSN sites CEBR, MGUE, NNOR were evaluated. Figure 11 shows examples of the coincidence between the curves resulting from different sources to compute the ionospheric media calibrations for the three test cases.
Figure 11 displays the slant range corrections [m] w.r.t. to a 2295 MHz reference frequency, Eq. (21), in the following colours:

Blue: Measured sTEC of the excluded GNSS satellite serving as reference;

Red: Weighted mean sTEC obtained from the other GNSS satellites (noisy curve), Eq. (18);

Violet: sTEC values computed from the CSP polynomial (this is always a smoothed version of the red weighted mean curve, the polynomial degree is 5);

Green: sTEC interpolated from a JPL IONEX file along the excluded GNSS satellite's trajectory;

Yellow: sTEC computed with the NTCMGL along the excluded GNSS satellite's trajectory;

Brown: cosθ, cosς, Eq. (20), should always be close to One;

Pale blue: Ξparameter, Eq. (19), should always be close to One.
In addition, Figure 11 presents the RMS values of differences:

RMSgw: measured sTEC of excluded GNSS satellite minus weighted mean sTEC;

RMSgj: measured sTEC of excluded GNSS satellite minus interpolated JPL IONEX sTEC;

RMSwj: weighted mean sTEC minus interpolated JPL IONEX sTEC.
The following aspects should be kept in mind when inspecting Figure 11:

the JPL IONEX files include always all available GNSS satellites, i.e. also the one selected in the analysis as reference;

what the blue curve (reference GNSS satellite) displays, are measurements too, i.e. also the reference sTEC values are subject to noise and inaccuracies;

the red weighted mean curve (then smoothed by the polynomial fit, violet CSP curve) and the pale blue Ξcurve display noisy behaviour. Since the weighting method is purely geometricbased, this noise may be induced by changing GNSS constellations (some satellites rise, while other set);

the analysis follows here the trajectory of the excluded GNSS satellite. While interplanetary missions are typically located in − and around the ecliptic, GNSS satellites can achieve considerably higher latitudes (GPS 55°, GLONASS 64°.8, Galileo 56°).
Figures 11a, c, d show the different curves for sTEC values derived from measured data; the only exception here is the NTCMGL curve (yellow). A systematic shift of the NTCMGL curve w.r.t. the other GNSS measurementsbased curves is clearly visible. This “bias” reflects the difference of the actual ionisation level w.r.t. a background ionisation level represented by the NTCMGL.
Figure 11b shows for the same pass/arc of Figure 11a the situation when all measured sTEC values were “simulated” with the NTCMGL. Only the JPL IONEX curve remains measurementsbased; this could not be changed. Now the measurementsbased green JPL curve displays nicely the offset between the actual level of ionospheric ionisation and the background ionisation level represented by the NTCMGL. Since here also the reference sTEC values were simulated with the NTCMGL, the blue curve is in this special case always identical with the yellow curve and thus not visible (the RMS values displayed in Figure 11b are meaningless). The analysis by taking NTCMGL simulated sTEC values as “observables” seems to confirm that the weighted mean sTEC method, as described in Section 4, works principally.
In order to get a more detailed picture about the overall performance of the new weighted mean sTEC method versus JPL, the three types of RMS values, displayed in Figure 11 (not Fig. 11b), were analysed for all considered passes/arcs and the three test cases in the form of histograms (Fig. 12). In Figure 12 the RMS types RMSgw, RMSgj and RMSwj (from left to right) of all passes/arcs (≈70 per test case) are shown in the form of histograms for each of the three test cases (from top to down). Since RMS values are naturally always positive, the resulting histograms comprise only the positive part of the abscissa. In the ideal case, the majority of RMS values should be small; thus the resulting histograms should display a peak close to the ordinate and then decay to the right along the abscissa. The strength of that decay and the height, location and width of the peak are an indicator on how good the media calibrations were able to reproduce the reference sTEC values of the excluded GNSS satellite under the different ionospheric conditions of the three test cases.
In addition to the histograms, Table 4 summarizes the mean, standard deviation, minimal and maximal RMS values of the different RMS types for the three test cases.
When inspecting Figure 12 and Table 4, one can come to the following conclusions:

under normal ionospheric conditions, Test Case 1, the overall distribution of the RMSgw values (weighted mean sTEC, red histogram) and the RMSgj values (JPL, green histogram) are of comparable shape. The RMSwj values (brown histogram) summarize then the differences between weighted mean sTEC and JPL IONEXderived sTEC. Smaller values for the mean (peak position) and the standard deviation (width) of the RMSgw values compared to the RMSgj values indicate a slightly better performance of the new weighted mean sTEC approach over the JPL IONEXderived sTEC numbers. The range of maxima and minima is of the same magnitude for both the RMSgw values and the RMSgj values;

under strong magnetic storm conditions, Test Case 2, mean and standard deviation are slightly smaller for the RMSgw values (weighted mean sTEC, red histogram) compared to the RMSgj values (JPL, green histogram) too. Also the range of outliers, indicated by the maxima and minima, is smaller for the RMSgw values than for the RMSgj values. Thus for Test Case 2, the RMSgw values are slightly favourite over the RMSgj values too;

during the solar eclipse, Test Case 3, mean and standard deviation of the RMSgw values (weighted mean sTEC, red histogram) are clearly smaller than for the RMSgj values (JPL, green histogram). On the other hand, Table 4 shows that the spread between maxima and minima is larger for the RMSgw values compared to the RMSgj values;

in all cases the RMSgj values (JPL, green histograms) seem to display more sensitivity to outliers (higher bars on the right part of the abscissa in Fig. 12) than the RMSgw values (weighted mean sTEC red histograms).
Fig. 12 Histogram analysis of RMS values. 
Mean, standard deviation, minimal and maximal RMS values of the three test cases.
5.2.1.1 Summary of validation method 1
The ionospheric media calibrations computed with the new weighted mean sTEC method, Section 4, seem to display a slightly better performance compared to the JPL IONEX data. As pointed out above, the reference GNSS satellite could not be excluded from the JPL data, in this way slightly favouring the JPL data. And as pointed out above too, the sTEC values from the reference GNSS satellite were also measurements subjected to noise and errors, i.e. they were not absolute. So in summary, the weighted mean sTEC media calibrations compare favourably to those computed from the JPL IONEX files.
5.2.2 Validation method 2
For the validation of the ionospheric calibrations in the context of a deepspace satellite orbit determination (OD), the same procedure and the same scenarios were used as the ones described in Section 5.1.
The legacy model for the ionospheric corrections is the NTCMGL (Jakowski et al., 2011), which was implemented in recent years and demonstrated to reliably improve the quality of the orbit determination compared to the previously used Klobuchar model (Klobuchar, 1975). This model is here compared in terms of OD statistics with the calibrations based on the approach described in Section 4, called IOCal in the following.
The effect of the ionospheric delay on range and Doppler observations is less pronounced than the tropospheric one for Xband satellites. The best available corrections for troposphere, i.e. the ones based on GNSS, need to be used in order to reduce the overall noise level of the tracking data and make the comparison of orbit determination statistics easier. For the same reason, for the LISA Pathfinder scenario only the test case with elevation larger than 10° was kept for this assessment.
In Table 5 the residuals and range bias statistics are reported for the GAIA scenario. It can be observed that no significant change is visible in terms of range and Doppler residuals; however, the range bias RMS increases by a factor two compared with the legacy model. This might indicate that some systematic source of error could still be present in the IOCal corrections. The reconstructed orbits differ by up to 450 m, which is in line with the formal uncertainty of the OD solution based on range and Doppler at declination lower than 10° (Budnik et al., 2012).
To better understand this somehow unexpected result, a dedicated investigation was carried out, which exposed some issues at the level of preprocessing the GNSS data, like an intersystem bias between the GPS and GLONASS constellation that was not properly accounted for. The results reported in this section already take into account the updated calibrations, that are in general more consistent with other sources of data and models compared to the initial runs, but still do not improve in a significant way the quality of the orbit determination. Further improvements at preprocessing level have already been identified, which require longerterm developments that cannot be incorporated in this paper results.
Table 6 shows the comparison of the OD statistics for the LISA Pathfinder test case. A minor but insignificant improvement can be observed for this scenario in the Doppler residuals as well as in the range biases, although its magnitude is essentially within the level of noise of the measurements. The two solutions can be considered equivalent both in terms of statistics and orbit, with a maximum difference of about 90 m, which is within the uncertainty region for this type of orbit tracked with range and Doppler only.
It has to be pointed out that due to time restrictions, an extensive validation of the ionospheric corrections generated with the new technique at this stage was still not possible, and different results could be obtained when adding more test scenarios. Also it has to be pointed out that range and Doppler observables acquired at Xband or higher frequencies are less sensitive to the ionospheric delay than deltaDOR measurements, thus further work needs to be carried out on Flight Dynamics side to validate the IOCal also with this last observation type.
GAIA OD statistics comparison for ionospheric corrections.
LISA Pathfinder OD statistics comparison for ionospheric corrections.
5.2.2.1 Summary of validation method 2 and open points
Although the ionospheric calibrations based on the described method displayed a relatively good agreement in most of the cases with the ones retrieved via IGS, the results of the validation in the OD environment did not show the expected improvement in the radiometric observation residuals, especially for the case of GAIA, which is tracked from all three ESA DSN sites, frequently during daynight transitions. This can be partly justified by the good performance of the NTCMGL model which, despite its simplicity, proved so far to be a good background model for all of those missions for which no other sources of ionospheric corrections were available. Further investigations are currently being carried out to improve the settings of the ionospheric calibration software and consequently the accuracy of the corrections.
6 Conclusions
The availability of accurate tropospheric and ionospheric delay information is essential for deepspace orbit determination. Until now, tropospheric and ionospheric media calibrations are computed by ESOC Flight Dynamics based on empirical background models. If arranged for certain missions, these media calibrations can be also provided by JPL. The media calibrations are processed in the form of CSPfiles. Both, tropospheric and ionospheric delays are represented in these CSPfiles in form of timenormalized polynomials directly providing corrections in terms of range delays. Flight Dynamics uses these media calibrations to correct the tracking data of ESA spacecraft (range, Doppler, deltaDOR), especially for the interplanetary missions. In the case of deltaDOR, these media calibrations are also required for quasar tracking. The ground stations considered are the three ESA DSN sites CEBR, MGUE, NNOR. The media calibrations system might be extended to further stations in future.
A new inhouse generation of media calibrations − troposphere and ionosphere − at ESA/ESOC, whereby the ESOC Navigation Support Office produces these media calibration products based on data accumulating from its routine products generation for the IGS anyway (GNSS orbits and clocks, ZTDs and ZWDs, sTEC data), replaces the old set up used by Flight Dynamics to generate its media calibrations. The target is to become as independent as possible from external data and models. The approach is primarily GNSSobservations based, recorded with own receivers at the ESA DSN sites and processed with own software and systems. Simple and robust algorithms were worked out, well suited for operational needs. In case of nonavailability of fresh GNSS data, the operational process is bypassed by the usage of simple background models.
The new service is now (March 2018) running in operational mode since 14 months. The results of the validation in the operational orbit determination context show that the use of GNSS based tropospheric corrections represents a significant improvement in the quality of radiometric observations modelling compared to the legacy Saastamoinen Model. The ionospheric calibrations, although partially validated with a dedicated procedure, did not demonstrate yet to improve significantly the orbit determination results, and further validation efforts are currently ongoing to improve their reliability and accuracy.
Acknowledgements
The authors are grateful to ESOC Flight Dynamics and the ESOC Navigation Support Office for providing the infrastructure, data and budget to realize this project. The authors thank the ESOC Ground Systems Engineering Department for the routine delivery of the WEST data needed for the tropospheric calibrations. The Navigation Support Office version of the NAvigation Package for Earth Observation Satellites (NAPEOS) software could be used to implement the media calibrations as new component. The authors are thankful for fruitful discussions with other ESOC colleagues improving the concept of the new modelling approach. Finally, the authors thank the journal Editor and the referees for their fruitful comments by evaluating the manuscript. The editor thanks two anonymous referees for their assistance in evaluating this paper.
Appendix A : IPP computation from station position vector and azimuth & elevation of satellite
In socalled single layer representations of the ionosphere (e.g. Schaer, 1999), as for instance used for the NTCMGL and the currently available IONEX files (www.igs.org), the whole ionospheric electron content is assumed to be condensed on an infinitesimal thin shell enveloping the Earth as hollow sphere at an altitude of typically 450 km. IPPs are the intersection points of satellite range vectors through that shell. This Appendix A describes the computation of IPPs from the given ground station position vector and azimuth/elevation of a satellite as observed at that station.
From the components of a given earthfixed Cartesian ground station position vector , unit vectors pointing into local North and East can be derived (Feltens, 2008): (A.1)
With the unit vectors and , a unit vector pointing at station location to the satellite can be established in two steps:
Fig. A.1 Computation of unit vector pointing to satellite from azimuth & elevation at station. 

Unit vector pointing into the direction of the projection of the satellite range vector onto the local horizon plane spanned up by and (Fig. A.1, greenred plane and dark blue vector):

Unit vector pointing into the direction of the satellite range (Fig. A.1, violet vector):
Knowing the unit vector pointing from station into satellite range direction, the position vector of the IPP can be computed (Fig. A.2):
Fig. A.2 Computation of the position vector of the IPP. 
The magnitude of the IPP position vector is (Fig. A.2, left): (A.4) where
R_{m} – mean earth radius, Eq. (10)
h_{I} – height of ionospheric shell, Eq. (10) Next, sine and cosine values of the auxiliary angles α and β have to be computed (Fig. A.2, right): (A.5)
Applying the cosine rule (e.g. Sigl, 1977): (A.6)
And considering: γ = 180^{0} − (α + β) and cos {180^{0} − (α + β)} = − cos(α + β), Eq. (A.6) can be modified to: (A.7)
And with cos(α + β) = cosα ⋅ cosβ − sinα ⋅ sinβ, Eq. (A.7) is modified to: (A.8)
Finally by setting and taking the square root of Eq. (A.8), one gets the formulae (A.9) for computing the IPP position vector: (A.9)
References
 Armano M. et al. 2016. SubFemtog Free Fall for SpaceBased Gravitational Wave Observatories: LISA Pathfinder Results. Phys Rev Lett 116: 231101. [CrossRef] [PubMed] [Google Scholar]
 Böhm J. 2007. Tropospheric Delay Modelling at Radio Wavelengths for Space Geodetic Techniques. Habilitation Thesis at the Vienna University of Technology, Institute of Geodesy and Geoinformation, Vienna, Austria, Geowissenschaftliche Mitteilungen, Iss. 80. [Google Scholar]
 Budnik F, Morley TA, Mackenzie RA. 2004. ESOĆs System for Interplanetary Orbit Determination: Implementation and Operational Experience, 18th International Symposium on Space Flight Dynamics, Munich (Germany), http://issfd.org/ISSFD_2004/papers/P0118.pdf. [Google Scholar]
 Budnik F, Morley T, Croon M. 2012. Orbit Reconstruction for the GAIA mission. 23rd Symposium on Space Flight Dynamics, Pasadena, California, http://issfd.org/ISSFD_2012/ISSFD23_OD2_5.pdf. [Google Scholar]
 Chao CC. 1973. A New Method to Predict Wet Zenith Range Correction from Surface Measurements. DSN Progress Report, Vol XIV, https://ipnpr.jpl.nasa.gov/progress_report2/XIV/XIVG.PDF. [Google Scholar]
 Coster AJ, Gaposchkin EM, Thornton LE. 1992. RealTime Ionospheric Monitoring System Using the GPS. Technical Report 954, ESCTR91179, ADA256 916, Lincoln Laboratory, Massachusetts Institute of Technology, Lexington, Massachusetts, https://pdfs.semanticscholar.org/fd9d/527c1747b72636159ac5ee6fcf1997e01e66.pdf. [Google Scholar]
 Davies K. 1990. Ionospheric Radio. IEE Electromagnetic Wave Series, 31, Peter Peregrinus Ltd., ISBN 0 86341 186 X. [Google Scholar]
 Estefan JA, Sovers OJ. 1994. A Comparative Survey of Current and Proposed Tropospheric RefractionDelay Models for DSN Radio Metric Data Calibration. JPL Publication 9424, National Aeronautics and Space Administration, Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA, https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19950013586.pdf. [Google Scholar]
 Feltens J. 2008. Vector methods to compute azimuth, elevation, ellipsoidal normal, and the Cartesian (X, Y, Z) to geodetic (φ, λ, h) transformation. J Geod 82: 493–504. [CrossRef] [Google Scholar]
 Jakowski N, Hoque MM, Mayer C. 2011. A new global TEC model for estimating transionospheric radio wave propagation errors. J Geod 85: 965–974, DOI: 10.1007/s0019001104551. [CrossRef] [Google Scholar]
 JPL. 2008. TRK 223 Media Calibration Interface. 820013, Deep Space Network (DSN), External Interface Specification, JPL D16765, Revision C (draft), http://dawndata.igpp.ucla.edu/download.jsp?file=documents/Gravity/DATA_SET_DESCRIPTION/TRK223_REVC_L5.PDF. [Google Scholar]
 Kelley MC. 1989. The Earth's Ionosphere, Plasma Physics and Electrodynamics, Academic Press Inc., Harcourt Brace Jovanovich, Publishers, San Diego New York Berkeley Boston London Sydney Tokyo Toronto. [Google Scholar]
 Klobuchar JA. 1975. A FirstOrder, Worldwide, Ionospheric, TimeDelay Algorithm. Bedford, Mass. United States Air Force, Air Force Systems Command, Air Force Cambridge Research Laboratories, Ionospheric Physics Laboratory. [Google Scholar]
 Komjathy A. 1997. Global Ionospheric Total Electron Content Mapping Using the Global Positioning System. Ph.D. dissertation, Department of Geodesy and Geomatics Engineering Technical Report No. 188, University of New Brunswick, Fredericton, New Brunswick, Canada. [Google Scholar]
 Lagler K, Schindelegger M, Böhm J, Krásná H, Nilsson T. 2013. GPT2: Empirical slant delay model for radio space geodetic techniques. Geophys Res Lett 40: 1069–1073, DOI:10.1002/grl.50288. [CrossRef] [Google Scholar]
 Mannucci AJ, Wilson BD, Yuan DN, Ho CH, Lindquister UJ, Runge TF. 1998. A global mapping technique for GPSderived ionospheric total electron content measurements. Radio Sci 33: 565–582. [CrossRef] [Google Scholar]
 Moyer TD. 1971. Mathematical formulation of the Double Precision Orbit Determination Program (DPODP), NASAJPL report JPLTR321527, Jet Propulsion Laboratory, https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19710017134.pdf. [Google Scholar]
 Moyer TD. 2000. Formulation for Observed and Computed Values of Deep Space Network Data Types for Navigation, Deep Space Communications and Navigation Systems, Jet Propulsion Laboratory, Publication 007, https://descanso.jpl.nasa.gov/monograph/series2/Descanso2_all.pdf. [Google Scholar]
 Niell AE. 1996. Global mapping functions for the atmosphere delay at radio wavelengths. J Geophys Res 101: 3227–3246. [NASA ADS] [CrossRef] [Google Scholar]
 Niell AE. 2000. Improved atmospheric mapping functions for VLBI and GPS. Earth Planets Space 52: 699–702. [CrossRef] [Google Scholar]
 Ratcliffe JA. 1972. An introduction to the ionosphere and magnetosphere, Cambridge University Press, Great Britain. [Google Scholar]
 Saastamoinen JJ. 1973. Contributions to the theory of atmospheric refraction. Bulletin Géodesique 107: 13–34. [CrossRef] [Google Scholar]
 Schaer S. 1999. Mapping and predicting the Earth's Ionosphere using the Global Positioning System, PhD Thesis, AIUB, University of Berne, Switzerland. [Google Scholar]
 Schaer S, Gurtner W, Feltens J. 1998. IONEX: The IONosphere Map Exchange Format Version 1. In: Proceedings of the IGS Analysis Centres Workshop. Darmstadt, Germany. [Google Scholar]
 Seeber G. 2003. Satellite Geodesy, 2nd edition, Walter de Gruyter, Berlin, New York. [CrossRef] [Google Scholar]
 Sigl R. 1977. Ebene und Sphärische Trigonometrie, Hermann Wichmann Verlag, Karlsruhe. [Google Scholar]
 Steigenberger P, Böhm J, Tesmer V. 2009. Comparison of GMF/GPT with VMF1/ECMWF and implications for atmospheric loading. J Geod 83: 943–951, DOI: 10.1007/s0019000903118. [CrossRef] [Google Scholar]
Cite this article as: Feltens J, Bellei G, Springer T, Kints MV, Zandbergen R, Budnik F, Schönemann E. 2018. Tropospheric and ionospheric media calibrations based on global navigation satellite system observation data. J. Space Weather Space Clim. 8: A30
All Tables
RMS summary on tropospheric polynomial fitting, 12 April 2016, 12:00–14 April 2016, 12:00.
Mean, standard deviation, minimal and maximal RMS values of the three test cases.
All Figures
Fig. 1 Time schedule for the daily computer runs at the Navigation Support Office to establish 48^{h} time series of ZWDs and ZTDs, and tropospheric media calibrations. 

In the text 
Fig. 2 Ideal case Amatrix structure. 

In the text 
Fig. 3 Examples of polynomial fitting in the case of data gaps, marked in red. 

In the text 
Fig. 4 Considered vector geometry between DSN station, GNSS satellites and ESA spacecraft/quasar. 

In the text 
Fig. 6 Maximal possible geocentric angle ξ_{max} for given mean earth radius R_{M} and ionospheric shell height h_{I}_{.} 

In the text 
Fig. 5 Computing weights from the angles enclosed by the IPP position vectors. 

In the text 
Fig. 11 (a) Example from Test Case 1; (b) example from Test Case 1 with sTEC simulated by NTCMGL; (c) example from Test Case 2; (d) example from Test Case 3. 

In the text 
Fig. 7 (a) Computing sTEC_{ESA/Q} only with Eqs. (6)–(9); (b) Measure 1 included; (c) Measure 2 included; (d) Measure 3 included → see Figure 11a. 

In the text 
Fig. 8 (a) Example of drypart polynomials; (b) example of wetpart polynomials. 

In the text 
Fig. 9 Overview of Doppler residuals over the full arc used in the GAIA OD; tropospheric corrections based on TRCal. 

In the text 
Fig. 10 Comparison of 2way Doppler residuals with Saastamoinen (top plot) and TRCal (bottom plot) corrections. 

In the text 
Fig. 12 Histogram analysis of RMS values. 

In the text 
Fig. A.1 Computation of unit vector pointing to satellite from azimuth & elevation at station. 

In the text 
Fig. A.2 Computation of the position vector of the IPP. 

In the text 