Coupling a weather model directly to GNSS orbit determination

Neutral gas atmosphere bends and delays propagation of microwave signals in satellite-based navigation. Weather prediction models can be used to estimate these effects by providing 3-dimensional refraction fields to estimate signal delay in the zenith direction and determine a low-dimensional mapping of this delay to desired azimuth and elevation angles. In this study, a global numerical weather prediction model (Open Integrated Forecasting System (OpenIFS) licensed for Academic use by the 5 European Centre for Medium-Range Weather Forecast) is used to generate the refraction fields. The ray-traced slant delays are supplied as such – in contrast to mapping – for an orbit solver (GROOPS (Gravity Recovery Object Oriented Programming System) software toolkit of the Technical University of Graz) which applies the raw observation method. Here we show that such a close coupling is possible without need for major additional modifications in the solver codes. The main finding here is that the adopted approach provides a very good a priori model for the atmospheric effects on navigation signals, as 10 measured with the midnight discontinuity of Global Navigation Satellite System (GNSS) satellite orbits. Our interpretation is that removal of the intermediate mapping step allows to take advantage of the local refraction field asymmetries in the GNSS signal processing. Moreover, the direct coupling helps in identifying deficiencies in the slant delay computation because the modelling errors are not convoluted in the precision-reducing mapping. These conclusions appear robust, despite the relatively small data set of raw code and phase observations covering the core network of 66 ground-based stations of the International 15 GNSS Service over one-month periods in December 2016 and June 2017. More generally, the new configuration enhances our control of geodetic and meteorological aspects of the orbit problem. This is pleasant because we can, for instance, regulate at will the weather model output frequency and increase coverage of spatio-temporal aspects of weather variations. The direct coupling of a weather model in precise GNSS orbit determination presented in this paper provides a unique framework for benefiting even more widely than previously the apparent synergies in space geodesy and meteorology. 20


Introduction
Refraction in the neutral gas atmosphere bend and delay global navigation system satellite (GNSS) signals (Bevis et al., 1992). These atmospheric effects cannot be deduced based on GNSS measurements alone because the signal propagation is identical for all frequencies typically applied in GNSS. Some auxiliary information is therefore needed to correct these numerically integrated over 24 hours based on force models, such as Earth's gravity field, tidal forces, and radiation pressure (cf. Strasser et al., 2019, for a complete list). These dynamic orbits are then fitted to the observations in the least-squares adjustment by estimating their initial position and velocity as well as a set of solar radiation pressure parameters (Arnold et al., 2015), as this force cannot be modeled adequately in advance. In addition, small instantaneous velocity changes, pseudostochastic pulses (e.g., Hugentobler and Montenbruck, 2017), are estimated at the center of each 24-hour orbit arc to consider 70 orbit modeling deficiencies.
Next to the orbit parameters, various other parameters are estimated in the least-squares adjustment. These comprise static station positions, Earth orientation parameters, epoch-wise clock errors and constant signal biases at each receiver and satellite, the ionospheric slant total electron content (STEC) per group of observations between a receiver and satellite at one epoch, and phase ambiguities. The ambiguities are resolved to integer values during processing. 75 Following Petit and Luzum (2010), the tropospheric slant delay for a line-of-sight observation between a satellite and receiver is T SD = m h (e)D zh + m w (e)D zw + m g (e)[G N cos a + G E sin a] . (1) Here, e and a are the elevation and azimuth angles at the receiver antenna, D zh is the zenith hydrostatic delay, D zw is the zenith wet delay, and G N and G E are the horizontal delay gradients in north-south and east-west directions (the delays are where a, b, and c are mapping coefficients (Herring, 1992;Niell, 1996). The gradient mapping function is m g (e) = 1 sin e tan e + C with C = 0.0031 for the hydrostatic effect and C = 0.0007 for the wet effect (Chen and Herring, 1997).
In the default system, a priori tropospheric corrections are based on VMF3 operational data (Landskron and Böhm, 2018).
The model provides zenith hydrostatic, zenith wet, hydrostatic gradient, and wet gradient delays as well as corresponding mapping function coefficients (a, b, c) on a global 1x1°grid at a 6-hour sampling period. Values for a specific station location and point in time are bilinearly interpolated from the surrounding grid points in the space domain and linearly interpolated in 90 time. A height correction from orographic height to station height is applied as well (cf. Strasser et al., 2019).
In the experiment, on the other hand, OpenIFS-based slant delay sky-views (see Section 2.3) are used to determine a priori tropospheric delays. These sky-views comprise slant delays at a regular 1x1°azimuth-elevation grid ray-traced directly at the station coordinates with an hourly sampling period. The slant delay for a specific observation is then computed by interpolating the gridded slant delays to the azimuth and elevation of the observation and linearly interpolating in time between the hourly 95 delay values. The interpolation in space domain is done linearly in azimuth direction and using a degree-5 polynomial in elevation direction to properly cover the rapid delay changes at low elevations. Since these slant delays are based directly on a weather model, they include all hydrostatic, wet, and gradient effects and no mapping functions need to be used.
For individual measurements, there is a discrepancy ∆T SD between the a priori tropospheric slant delay and the actual delay affecting the measurement since the models are imperfect. If we assume that ∆T SD is solely due to troposphere, we can write 100 the tropospheric slant delay discrepancy by using Eq. 1 as Assuming that the hydrostatic delay is well modelled (∆D zh = 0), the residual zenith wet (∆D zw ) and gradient (∆G N , ∆G E ) delay parameters can be set up to consider the discrepancy between modeled and measured tropospheric influence. These station-wise parameters are set up or omitted in this study depending on the analysis. If included, ∆D zw is parameterized 105 as a degree-1 spline with 2-hourly nodes, while ∆G N and ∆G E are parameterized linearly over the 24 hours. The mapping functions required for this parametrization are determined using the VMF3 mapping coefficients. The same parametrization is used both in case of VMF3-based or OpenIFS-based a priori slant delays since the parametrization is independent of the a priori model.

Weather model 110
OpenIFS is a portable version of the Integrated Forecasting System (IFS) of the European Centre for Medium-Range Weather Forecasts (ECMWF). Essentially, it is a global weather prediction model with identical forecast skill as in the full IFS. Data assimilation is not included in OpenIFS, and thus external initial conditions of atmospheric state (temperature, wind, humidity, and surface pressure) are required. In this study, we use operational atmospheric analyses of ECMWF (Ollinaho et al., 2021) and OpenIFS version 43r3v1 that was part of the operational forecasting system at ECMWF from July 2017  between the twice-a-day analysis times.

Slant delays
The tropospheric slant delays are computed with ray tracing based on a least traveltime (LTT) operator (Eresmaa et al., 2008b). LTT operates on a 2-dimensional plane defined by the satellite and receiver positions and the centre of the Earth.
A 2-dimensional refractivity field is obtained by converting the 3-dimensional OpenIFS fields of atmospheric dry mass and 130 moisture content to refractivity and interpolating to a desired 2-dimensional plane. Interpolation is bi-linear in horizontal and assumes exponential refractivity profile in vertical between the model levels.
The LTT algorithm performs ray tracing to construct the path, expressed in polar coordinates, which satisfies the following system of differential equations: Here ds is a path element, r is distance from the Earth's centre, ψ is the counterpart for polar coordinates, θ is the zenith angle, and n the refractive index. Equations 5 are integrated with the fourth-order Runge-Kutta method starting from the receiver position (r rec = R + h rec ; ψ rec = 0) in the initial direction stated as the geometrical zenith angle of the satellite 140 (θ rec = θ geom ). The integration ends when the satellite altitude is reached (r end = r sat ). The integration yields a set of points (r; ψ) which satisfy these equations.
The total delay on a slanted path results from slowing down of the signal due to refractive index n > 1 and increasing signal path length due to signal bending. The path bending is due to the second term in Eq. 5c. Since geometrical zenith angle is used as the initial direction, an angular separation appears between the end point of the ray and the satellite location, i.e., r end = r sat 145 but ψ end = ψ sat . Therefore the calculation of the path and slant delay is repeated by using an updated θ rec as follows: In our implementation (Eresmaa et al., 2008b), the final slant delay is a linear combination of these two LTT calculations.
Additional accuracy of the starting zenith angle could be obtained with additional iterations.

150
The effect of Earth flattening was evaluated for the slant delay computation since the LTT algorithm accounts for the Earth oblateness applying the concept of Euler radius. The magnitude of this effect was evaluated by degrading the Earth ellipsoid to a sphere thus changing the 2-dimensional plane projection onto the Earth surface. The effect was found to vary from −0.5 cm to +1.5 cm depending on the receiver latitude and antenna azimuth angle.
Here the LTT solver is applied such that instead of computing ray paths exactly corresponding individual GNSS measurement 155 geometries, so-called sky-views are generated. They are constructed so that slant delays are calculated for zenith angles from 1 to 89 degrees and azimuth angles from 0 to 359 degrees with one degree increments in both directions, added with one calculation for the zenith delay. The sky-views are computed using the OpenIFS data at every full hour in December 2016 and June 2017 for the 66 selected IGS stations constituting the core network.

Performance metrics 160
In order to compare the experiments with the default system, the following performance metrics are defined. First, GNSS satellite orbits are determined for 24-hour periods. These 24-hour orbit arcs overlap at midnight (i.e., 00 UTC) between consecutive days, providing two independent positions r (expressed in meters) at a single epoch for each satellite. The differences between these positions are called orbit midnight discontinuities, as they represent jumps between two smooth orbit arcs. The discontinuity δr in the orbit position between consecutive days has components along (δr a ) and cross (δr c ) the orbit direction and a 165 radial (δr r ) component in the local vertical. Running GROOPS using only a priori models (i.e., not estimating any corrections through fitted parameters according to Eq. 4) and analyzing the discontinuities we can compare the goodness of the modelsthe smaller the discontinuities are, the better is the a priori modelling. The associated performance metric δr is defined as: δr = (δr a ) 2 + (δr c ) 2 + (δr r ) 2 for each satellite at each day boundary.

170
Second, as a part of the orbit determination, the a priori model for the troposphere can be corrected by fitting parameters in the least squares process (Eq. 4). The fitted residual zenith wet delay (∆D zw ) and gradient delays (∆G N and ∆G E ) are indicative of the goodness of the a priori model used -the smaller the corrections are, the better is the accuracy of the a priori model.

175
The following results are aimed to demonstrate how the newly assembled experimental configuration (OpenIFS) performs in relation to a well-established default system (VMF3). According to the performance metrics explained in section 2.4, an analysis over the midnight discontinuities has been carried out, along with a study of the fitted parameters for the month of

Orbit midnight discontinuities
The orbit midnight discontinuities are a good metric to analyze how any factor, for instance tropospheric modelling, affects the quality of GNSS products. Figure 1a shows a color map where each cell represents the difference in the midnight discontinuity between the a priori default and experiment systems for a satellite over two consecutive days, as measured with δr (see eq. 7 185 and expressed in units of cm. Here, blue (red) means that discontinuities in the experimental system are smaller (larger) than in the default system. The differences in Fig. 1a are generally very small, of the order of a few cm. The mean covering all satellites over the whole period is 0.025 cm, pointing to slightly smaller discontinuities in the experimental system, although statistically negligible.
Next is shown a time-series of the mean values of the norm δr for all satellites (Fig. 1b), where the default system is in blue 190 and the experiment in orange. Again, the results indicate that the experimental system has slightly smaller discontinuities. The mean values for the entire month are 1.961 cm for the experimental system and 1.986 cm for the default system. This difference, albeit small, is a rather systematic feature in Fig. 1b. Standard deviation around the mean value is 0.903 cm in the experimental system and 0.907 cm in the default system, showing a similar behaviour for both systems. The results presented in this section have been calculated not using Satellite G04 data due to malfunctions in the satellite or satellite-specific processing issues 195 which are unrelated to tropospheric modelling, that were deteriorating the results for both systems used in the experiment.

Fitted tropospheric parameters
Next, the magnitude of the fitted parameters (∆D zw , ∆G N and ∆G E ) is analyzed. Figure 2a shows The mean fitted zenith wet delay (Figure 2a) is 1.220 cm in the experimental system and 0.016 cm in the default system. This is indicative of a severe and consistent positive bias in the experimental system while it is negligible for the default system.
The positive bias can be interpreted such that in the experimental system the total slant delay is too small (i.e., the increased zenith wet delay compensates this deficiency in the least squares process). A test was carried out by segregating the receiver stations to three zones: Northern Hemisphere, Southern Hemisphere and the Tropics. This test does not reveal the sources of

215
In order to add confidence on the results presented here, the experiment was repeated for the month of June 2017. This period covers a different phase of the Earth on its orbit around the sun, and thus generally different circumstances for the orbit problem.
The central results are presented in Appendix A, and can be summarized as being very similar compared to the main study period. The conclusion is thus that the main study period of December 2016 seems to be representative for this type of studies.
4 Discussion 220 The default system of this study applies VMF3 delays and mapping function coefficients, which are computed from the ECMWF weather forecasts. The operational model version in December 2016 was very high-resolution, about 8 km grid spacing at 137 model levels. The disseminated forecast fields for VMF3 computation are however interpolated from the full model resolution to one by one degree horizontal resolution at 25 standard pressure levels, which are available at six-hourly intervals. Lastly, VMF3 coefficients can be interpreted as a low-order representation of the tropospheric delay which does not 225 represent, by design, azimuthal asymmetries in a receiver station sky-view. Despite these processing features, which all imply some loss of atmospheric information, VMF3 has proven to provide a hard-to-beat benchmark also in this study.
The experimental system applies a lower-resolution OpenIFS model version with about 31 km grid spacing at 91 levels, which rather closely corresponds to the system applied in ensemble prediction at ECMWF. In this respect, the default system has a clear advantage in terms of simulation accuracy of atmospheric dynamics and physical processes. On the other hand, the 230 subsequent processing in the experimental system is almost loss-less. We do interpolate the total tropospheric delay from a one by one degree sky-view to azimuth and elevation angles of an individual measurement, and interpolate in time from one-hourly model output fields to the measurement time. Other than that, the native OpenIFS refraction field is represented as such through the ray tracing step to the orbit solver. Technically, the ray tracing could be solved separately and directly corresponding to the azimuth and elevation angles of each measurement -in this case, time interpolation should also be omitted by gathering the 235 measurements into time-slots closest to the model output time. Note that the OpenIFS model time step is 15 minutes in the experimental system and atmospheric fields could be output at the same frequency to reduce time interpolation effects. Figure 2a is indicative of a sizable bias in the tropospheric a priori correction based on OpenIFS, where ∆D zw is at the level of 0.73 cm throughout the period. This result prompted us to investigate more closely the total delay mean differences between VMF3 and OpenIFS. Indeed, at very low elevation angles the total delay in the experimental system differs systematically from 240 the default system. This is very likely due to an approximation made in the LTT code regarding the geometrical effect on the signal path and the term cos θ nr ( ∂n ∂ψ ) r in Eq. 5c which is missing. Therefore, the ray tracer applied in the experimental system does not fully follow Rodgers (2000). Work is ongoing to update the code. Importantly, the experimental system provides us with a solid reference where the weather model-based data is directly interfaced with GNSS measurements and orbit solutions with no intermediate and possibly obscuring processing steps. Thanks to this property, it led us to detect the problem with the 245 implemented ray tracing code.
The following question remains. Figure 2 tends to suggest that the tropospheric a priori model in the experimental system is biased. How is it then possible that midnight discontinuities (Fig. 1) are comparable or even slightly smaller in the experimental system compared to the default system? The most likely explanation is that the azimuthal asymmetries of the tropospheric delay that are present in the experimental system but are missing from the default system do matter and contribute to the orbit 250 solutions. Essentially, the asymmetries are not systematic but they can systematically impact the system performance and the metrics presented here.
Finally, we have to keep in mind that only a small network of 66 stations are included in this study over a limited time period.
A more comprehensive experiment will be prepared later including better GNSS station coverage and enhanced tropospheric modelling. Our hypothesis is that in such a system, more atmospheric information is introduced in near-native format to the 255 orbit problem and it improves the solution accuracy.

Conclusions
Troposphere and stratosphere delay the propagation of navigation satellite signals and can lead to large errors in GNSS satellite orbit determination if not properly accounted for. In this article, a global numerical weather prediction model -OpenIFS of the European Centre for Medium-Range Weather Forecasts -is applied to generate atmospheric data of pressure, temperature, 260 and humidity. These are converted to 3-dimensional atmospheric refraction fields and used as an input for ray-tracing to solve the least travel time signal paths from satellites to receiver stations and to compute the associated delays. These delays are then directly used as a priori corrections of the atmospheric effects to solve GNSS satellite orbits with the GROOPS software toolkit of Graz University of Technology, Austria.
This new configuration to solve the GNSS orbit problem contains two novel aspects. First, the direct use of tropospheric 265 delays on slanted signal paths allows to fully account for the azimuthal asymmetries in the atmospheric refraction field, in contrast to traditional mapping functions which regularize the refraction field. Second, the intimate coupling of the numerical weather prediction model with the orbit solver allows to control the information flow between the two modules, including output frequency of the weather model, for instance.
The main finding here is that the new configuration provides good consistency of GNSS satellite orbits as measured with the 270 so-called orbit midnight discontinuities, i.e., how much satellite orbit initial positions need to be corrected between subsequent 24-hour orbit solutions. Another important finding is that the new configuration provides a solid reference where weather model-based data is directly interfaced with GNSS measurements and orbit solutions. Since the new configuration has fewer intermediate processing steps than the state-of-the-art methods, more direct diagnosis of the system performance is possible. In particular, this led us to detect a bias in our ray tracing computation that was undetected until now. Finally, the results indicate 275 that azimuthal asymmetries in the tropospheric delays, which are well-preserved in the new configuration, contribute to the accuracy of orbit solutions. These asymmetries are not systematic but they have a systematic impact on the orbit solutions.
At a more general level, the new configuration provides an enhanced control of geodetic and meteorological aspects of the orbit problem. This will allow us to widely benefit the apparent synergies in space geodesy and meteorology.