Least travel time ray tracer version 2 (LTT v2) adapted to the grid geometry of the OpenIFS atmospheric model

Electromagnetic signals commonly used in geodetic applications, such as the Global Navigation Satellite System (GNSS), undergo bending and delay in the neutral gas atmosphere of the Earth. The least travel time (LTT) concept is one of the approaches to model signal slant delays via a ray tracing (RT) procedure. In this study, we developed an LTT-based RT algorithm (LTT v2), where the three-dimensional refractivity field of the atmosphere is based on the atmospheric model data. This representation is complete in a sense that the domain of the RT conforms to the native grid geometry of the atmospheric model. In principle, the LTT-based RT algorithm is seen as an extension of an atmospheric model for signal delay evaluation. The atmospheric states are generated using a global numerical weather prediction model, the Open Integrated Forecast System of the European Centre for Medium-Range Weather Forecasts. In the LTT v2 model, some physical and numerical approximations are improved compared to the original implementation, called “LTT v1”. We compare the slant delays products of the two models. Additionally, a comparable modelling setup is created with the state-of-the-art VieVS Ray Tracer (RADIATE). The skill of slant delay estimation is assessed using metrics that are indicative of the quality of GNSS products derived using the GROOPS (Gravity Recovery Object Oriented Programming System) orbit solver software toolkit of the Graz University of Technology. The metrics used are GNSS orbit midnight discontinuities (MDs) and residuals of ground station precise positioning with respect to the IGS14 reference. Employment of slant delay products of the LTT RT algorithm in GNSS processing shows similar performance with v1 and v2. The GNSS orbit MDs are reduced by around 3 % when using the LTT v2 model, while root-mean-square residuals of ground station precise positioning are 5 % lower with LTT v1. The consistency of both metrics is improved slightly using LTT v2, as seen by the metrics' standard deviation values. Intercomparison with RADIATE indicates significantly better performance of LTT v2, which we attribute entirely to the much larger amount and lossless utilization of weather model data as input to LTT v2 versus RADIATE.

1 Introduction

The neutral gas atmosphere of the Earth delays propagation of electromagnetic waves. The delay develops due to wave group velocity being less than in the vacuum, i.e. the speed of light. In addition, the inhomogeneity of air invokes the bending of a wavefront from a straight line. Thus, the time delay of an electromagnetic signal in the Earth's atmosphere consists of delay along the path and excess of the path length due to bending, both of which are often measured in length units. In space geodesy, the measured distances are thus larger than the geometrical distances.

Signals of the Global Navigation Satellite System (GNSS) are affected by the neutral gas atmosphere (a.k.a. “tropospheric”) delay. The tropospheric delay is one of the major error sources in GNSS processing (Landskron and Böhm, 2018; Wilgan et al., 2017). The state-of-the-art modelling setup, described in the International Earth Rotation and Reference Systems Service Conventions (Petit and Luzum, 2010), represents the tropospheric delay between a surface site and a satellite as an empirical mapping function with adjustable parameters (Böhm et al., 2006): zenith wet and hydrostatic delays (ZWDs, ZHDs), mapping coefficients (providing elevation dependence) and so-called tropospheric gradients (providing azimuth dependence). The a priori mapping parameters are estimated using ray tracing (Landskron and Böhm, 2018; Hofmeister and Böhm, 2017). Ray tracing (RT) is a direct approach to computing the delay by solving the signal propagation through given atmospheric conditions in a specified direction. Quality of the RT procedure is therefore of great importance for obtaining accurate GNSS products. Various signal RT methods are available for geodetic applications. Hofmeister (2016) developed the RADIATE ray tracer for signals of the microwave and the optical frequency ranges, based on the Eikonal equation. Another example is the least travel time (LTT) ray tracer by Eresmaa et al. (2008). It is based on Snell's law equation expressed in a polar coordinate system, as explained by Rodgers (2000). These slant delay models are designed differently and are simplified in one way or another by assumptions in favour of efficiency and utility.

A signal RT algorithm must be supplied with a three-dimensional representation of the atmospheric refraction index. Generally, the microwave atmospheric refraction index is a function of pressure, temperature and relative humidity (Bevis et al., 1994). These crucial variables usually comprise the atmospheric states provided by numerical weather prediction (NWP) models. NWP systems, whose data can be employed for signal RT, are continuously improving their spatial and temporal representation of flow dynamics and physical processes as well as advancing the use of Earth observations in data assimilation (Bauer et al., 2015). From the data assimilation perspective, an RT algorithm can be seen as an operator from model state space (model variables) to observation space (delays of signals traversing the atmosphere). This, in fact, motivates implementing an RT algorithm that is structurally integrated into an NWP system, interconnecting the advances in both NWP and signal delay modelling. It has been previously shown that the use of signal delays in NWP data assimilation requires an accurate RT model (Järvinen et al., 2007). With this motivation, we revised the LTT RT model by ensuring lossless utilization of atmospheric information and lifting various unphysical constraints. We developed a new software with a reworked LTT ray tracer at its core: the LTT version 2 (LTT v2). This software implementation is designed to optimally utilize NWP data produced by the Open Integrated Forecasting System model (OpenIFS) of the European Centre for Medium-Range Weather Forecasts (ECMWF). The software can be modified to other NWP data input, if needed.

The new LTT v2 model (Vasiuta, 2024) is compared with two other models which are capable of using the OpenIFS input data. The first reference is the LTT code by Eresmaa et al. (2021), which is a modification of the original LTT implementation by Eresmaa et al. (2008). This model is referred to as “LTT v1”. The second reference is the RADIATE ray tracer. Given its distinct way of employment of weather data, a special configuration of the OpenIFS atmospheric states is necessary before input into RADIATE. To ensure a fair comparison, a specific configuration is established for RADIATE, providing the ray tracer with the same weather data and aligning its output structure to match that of LTT v1 and v2. This way, all three models produce 6 hourly slant delays snapshots for each ground station used in the GNSS processing.

We test the quality of the slant delay models by solving the GNSS precise orbit determination (POD) using a direct coupling of POD and NWP models (Navarro Trastoy et al., 2022). Estimation of tropospheric parameters is omitted, and a priori slant delays from ray tracers are used throughout the processing. We compute orbit and station position solutions based on observations of Global Positioning System (GPS) constellation in December 2016 using GPS multi-frequency measurements from 205 ground stations. The GNSS processing is performed with the GROOPS toolkit (Mayer-Guerr et al., 2021), which is modified to enable various representations of tropospheric delay estimates. The quality of GNSS products is assessed using two sets of metrics: orbit midnight discontinuity (MD) and precise point positioning (PPP) deviations. The MD metric evaluates the quality of batch orbit estimation, which inherently has a disconnection of two consecutive batch POD solutions both in position and velocity (Montenbruck and Gill, 2012). We assess positional MD, not the velocity. The PPP-related metrics are root-mean-square deviation and standard deviation from a reference station position. These PPP metrics inform of the accuracy of solutions and show the consistency of the PPP procedure, respectively. This analysis is utilized to determine the relative skill of the three ray tracers: LTT v1, LTT v2 and RADIATE.

In Sect. 2, we explain the state-of-the-art concept of slant delay modelling via RT. In Sect. 3, we describe our methodology of slant delay modelling, as well as the features of LTT v2. Section 4 explains the skill assessment of the three slant delay algorithms via GNSS processing. The slant delay product comparison and the skill assessment are reported in Sect. 5. Finally, Sect. 6 attempts to attribute the skill of slant delay modelling to different factors.

2 Background

2.1 Atmospheric ray tracing methods

The classical approach to ray-tracing (RT) a signal in a refractive medium, such as the Earth's atmosphere, is by using the so-called Eikonal equation. The Eikonal equation (Hamilton, 1828) describes an electromagnetic wave propagating through a slowly varying medium, connecting physical (wave) optics and geometric (ray) optics. It reads

$$\begin{array}{}\text{(1)}& \Vert \mathbf{\nabla }L{\Vert }^{\mathrm{2}}=n(\mathit{r}{)}^{\mathrm{2}},\end{array}$$

where L is the optical path length, ∇L the components of ray direction, n the known refractive index of the medium and r the position vector. The surfaces of equal optical length (L(r)=const.) are called geometrical wave surfaces or geometrical wavefronts.

This equation is usually described in the Hamiltonian canonical formalism (Born et al., 1999; Nilsson et al., 2013):

$$\begin{array}{}\text{(2a)}& {\displaystyle }& {\displaystyle }H(\mathit{r},\mathbf{\nabla }L)\equiv {\displaystyle \frac{\mathrm{1}}{a}}\left\{(\mathbf{\nabla }L\cdot \mathbf{\nabla }L{)}^{\frac{a}{\mathrm{2}}}-n(\mathit{r}{)}^a\right\}=\mathrm{0},\text{(2b)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\mathit{r}}{\text{ d}u}}={\displaystyle \frac{\partial H}{\partial \mathbf{\nabla }L}},\text{(2c)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\mathbf{\nabla }L}{\text{ d}u}}=-{\displaystyle \frac{\partial H}{\partial \mathit{r}}},\text{(2d)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}L}{\text{ d}u}}=\mathbf{\nabla }L\cdot {\displaystyle \frac{\partial H}{\partial \mathbf{\nabla }L}},\end{array}$$

where a is a value related to the parameter of interest u, and H(r,∇L) is the Hamiltonian. When applying RT to the computation of delay along the path in the atmosphere, the parameter u is the length s along the ray, a equals to 1 and a spherical coordinate system is used (Hofmeister, 2016; Nafisi et al., 2012; Böhm and Schuh, 2013). The formulation of the tropospheric delay problem is thus the following:

$$\begin{array}{}\text{(3a)}& {\displaystyle }& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }H(r,\mathit{\vartheta },\mathit{\lambda },L_r,L_{\mathit{\vartheta }},L_{\mathit{\lambda }})\\ {\displaystyle }& {\displaystyle }\equiv {\left(L_r^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}{L}_{\mathit{\vartheta }}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\vartheta }}}{L}_{\mathit{\lambda }}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}-n(r,\mathit{\vartheta },\mathit{\lambda },t)=\mathrm{0},\end{array}\text{(3b)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}r}{\mathrm{d}s}}={\displaystyle \frac{\mathrm{1}}{\mathit{\omega }}}{L}_r,\text{(3c)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\mathit{\vartheta }}{\mathrm{d}s}}={\displaystyle \frac{\mathrm{1}}{\mathit{\omega }}}{\displaystyle \frac{L_{\mathit{\vartheta }}}{r^{\mathrm{2}}}},\text{(3d)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\mathit{\lambda }}{\mathrm{d}s}}={\displaystyle \frac{\mathrm{1}}{\mathit{\omega }}}{\displaystyle \frac{L_{\mathit{\lambda }}}{r^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\vartheta }}},\text{(3e)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{L}_r}{\mathrm{d}s}}={\displaystyle \frac{\partial n(r,\mathit{\vartheta },\mathit{\lambda },t)}{\partial r}}+{\displaystyle \frac{\mathrm{1}}{\mathit{\omega }r}}\left({\displaystyle \frac{L_{\mathit{\vartheta }}^{\mathrm{2}}}{r^{\mathrm{2}}}}+{\displaystyle \frac{L_{\mathit{\lambda }}^{\mathrm{2}}}{r^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathit{\vartheta }}}\right),\text{(3f)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{L}_{\mathit{\vartheta }}}{\mathrm{d}s}}={\displaystyle \frac{\partial n(r,\mathit{\vartheta },\mathit{\lambda },t)}{\partial \mathit{\vartheta }}}+{\displaystyle \frac{\mathrm{1}}{\mathit{\omega }}}{\displaystyle \frac{L_{\mathit{\lambda }}^{\mathrm{2}}}{r^{\mathrm{2}}{\sin }^{\mathrm{3}}\mathit{\vartheta }}},\text{(3g)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{L}_{\mathit{\lambda }}}{\mathrm{d}s}}={\displaystyle \frac{\partial n(r,\mathit{\vartheta },\mathit{\lambda },t)}{\partial \mathit{\lambda }}},\text{(3h)}& {\displaystyle }& {\displaystyle }d=L-G=\underset{S}{\int }n(r,\mathit{\vartheta },\mathit{\lambda },t)\mathrm{d}s-G,\end{array}$$

where ds is the path element, r the radial distance, ϑ the co-latitude ($\mathrm{0}\le \mathit{\vartheta }\le \mathit{\pi }$), λ the longitude ($\mathrm{0}\le \mathit{\lambda }\le \mathrm{2}\mathit{\pi }$) and ω is equal to refractive index $n(r,\mathit{\vartheta },\mathit{\lambda },t)$. The ray propagation is solved by integrating Eqs. (3b)–(3g) with respect to (r,∇L) starting from some initial value (r₀,∇L₀). The delay d is obtained by subtracting geometrical distance G from optical path length L along the path S (Eq. 3h).

To solve Eqs. (3b)–(3g) in a general three-dimensional case, one must develop a rule to compute gradients of the refractive index n, which is a challenging task in practice. Furthermore, the computational complexity of a strict approach surges with an increase in zenith angle and atmospheric data resolution. Thus, many practical ray tracer implementations benefit in computation and design by simplifying refractive index gradient effects and adding other assumptions. For example, a three-dimensional RT can be reduced to a two-dimensional RT by setting horizontal gradient components of the refractive index to zero: $\partial n/\partial \mathit{\vartheta }=\mathrm{0}$ and $\partial n/\partial \mathit{\lambda }=\mathrm{0}$. In this way, the ray propagation is limited to a plane of fixed azimuth.

The analytical approach of Thayer (1967) further simplifies the computation by assuming the refractive index is only dependent on height and refractivity between two neighbouring layers to be approximated using a power-law relation. In this scenario, the refractive index has only one vertical profile; hence RT constrains slant delays to have azimuthal symmetry. Another 2D RT approach is the so-called “piecewise-linear ray tracing” explained by Hobiger et al. (2008) and Hofmeister (2016). In this approach, the ray is propagated according to simple Snell's law, where the ray bends only at the half-height layer. The horizontal refractivity gradient is also neglected here, implying azimuthal asymmetry since refractivities on the vertical layers are obtained by bi-linear interpolation from a horizontal grid. Both approaches are employed in the RT program RADIATE, with the “piecewise-linear” method being the default option.

An alternative design of atmospheric signal RT was proposed by Rodgers (2000). In their approach, Rodgers (2000) considered the Eikonal Eq. (1) in a two-dimensional coordinate system, instead of restricting the general 3D case. From the 2D Eikonal equation, the law of ray propagation was derived as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}\mathit{ϵ}}{\mathrm{d}s}}={\displaystyle \frac{\mathrm{1}}{n(s,t)}}{\displaystyle \frac{\partial n(s,t)}{\partial t}},\end{array}$$

where ϵ is the direction of propagation and s, t the coordinates along and perpendicular to the ray, respectively. By expressing Eq. (4) in polar coordinates and introducing angle θ between the ray direction and the radius vector, the system of differential equations is obtained as follows:

$$\begin{array}{}\text{(5a)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}r}{\mathrm{d}s}}=\cos \mathit{\theta },\text{(5b)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\mathit{\psi }}{\mathrm{d}s}}={\displaystyle \frac{\sin \mathit{\theta }}{r}},\text{(5c)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}s}}=-\sin \mathit{\theta }\left({\displaystyle \frac{\mathrm{1}}{r}}+{\displaystyle \frac{\partial n(r,\mathit{\psi })}{\partial r}}\right)+{\displaystyle \frac{\cos \mathit{\theta }}{n(r,\mathit{\psi })\cdot r}}{\displaystyle \frac{\partial n(r,\mathit{\psi })}{\partial \mathit{\psi }}}.\end{array}$$

Here ds is a path element, r the distance from the Earth's centre of curvature (a.k.a. the Euler radius of curvature), ψ the horizontal counterpart for polar coordinates, θ the local zenith angle and n the refractive index. Equations in (5) are integrated numerically starting from the ground receiver position ($r_{\text{rec}};{\mathit{\psi }}_{\text{rec}}=\mathrm{0}$) in the initial direction defined as θ₀. Since θ₀ is initially unknown, finding a signal path between ground receiver and satellite, satisfying Eq. (5), is an initial value problem. The resulting tropospheric slant delay d is expressed as follows:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}d& =\underset{S}{\int }n(r,\mathit{\psi })\mathrm{d}s-G\\ & =\underset{S}{\int }\left(n\right(r,\mathit{\psi })-\mathrm{1})\mathrm{d}s+\left(\underset{S}{\int }\mathrm{d}s-G\right)\end{array}.\end{array}$$

Equation (6) explicitly shows two components of the delay: delay along the path and the excess of the path length due to bending, so-called geometrical delay.

2.2 The least travel time algorithm version 1

The least travel time (LTT) algorithm by Eresmaa et al. (2021) is based on the theory of Rodgers (2000). In LTT, the starting direction is set as the geometrical zenith angle of the satellite (θ₀=θ_g). The integration ends at the satellite altitude (r_end=r_sat). Yet, due to bending, an angular separation appears between the endpoint of the ray and the satellite location, i.e. r_end=r_sat but ψ_end≠ψ_sat. Therefore, the ray propagation is repeated by using an updated ${\mathit{\theta }}_{\mathrm{0}}^{\prime }$ as follows:

$$\begin{array}{}\text{(7)}& {\mathit{\theta }}_{\mathrm{0}}^{\prime }={\mathit{\theta }}_{\mathrm{0}}-({\mathit{\psi }}_{\text{end}}-{\mathit{\psi }}_{\text{sat}})={\mathit{\theta }}_{\mathrm{0}}-\mathrm{\Delta }\mathit{\psi }.\end{array}$$

The final slant delay is approximated as a linear combination of the two delay estimates with different starting zenith angles following Eqs. (6) and (7).

Several assumptions are made in the LTT v1. First, the ambiguity of the initial zenith angle θ₀ is not solved strictly. Also, the geometrical delay term in Eq. (6) is neglected. Moreover, only the radial refractive index gradient is present ($\frac{\partial n(r,\mathit{\psi })}{\partial \mathit{\psi }}$ is zero in Eq. 5c). These assumptions make LTT v1 computationally efficient and convenient to use in data assimilation of GNSS signal delays in NWP models employing the LTT operator's adjoint counterpart (Eresmaa and Järvinen, 2006; Järvinen et al., 2007). However, they add nonphysical uncertainty to the slant delay products.

2.3 Coordinate system and gravity

The ray tracing (RT) of a signal in the atmosphere deals with geometrical variables (lengths, angles). Hence, for the RT procedure, a reference coordinate system, which is consistent with the geometry of the Earth, should be selected. On the other hand, atmospheric data from numerical weather models are typically based on a coordinate system consistent with the Earth's gravity. To use the information from NWP models in RT applications, it is necessary to bring station position data and atmospheric model grid to the same height system. The ground station height is usually expressed with respect to a reference ellipsoid, while the heights in NWP models are with respect to a geopotential. Transformations for a position at a geographic location (ϕ,λ) between geopotential height and orthometric height, as well as orthometric height and ellipsoidal height, are as follows (Hobiger et al., 2008):

$$\begin{array}{}\text{(8a)}& {\displaystyle }& {\displaystyle }H={\displaystyle \frac{\mathrm{\Phi }}{\stackrel{\mathrm{‾}}{\mathit{\gamma }}(\mathit{\varphi },\mathit{\lambda },H)}}={\displaystyle \frac{\mathit{\zeta }\cdot g_{\mathrm{0}}}{\stackrel{\mathrm{‾}}{\mathit{\gamma }}(\mathit{\varphi },\mathit{\lambda },H)}},\text{(8b)}& {\displaystyle }& {\displaystyle }h=H+N(\mathit{\varphi },\mathit{\lambda }),\end{array}$$

where Φ is the geopotential value, g₀ the standard acceleration due to gravity, ζ the geopotential height, $\stackrel{\mathrm{‾}}{\mathit{\gamma }}(\mathit{\varphi },\mathit{\lambda },H)$ the mean acceleration due to gravity, H the orthometric height, N the geoid undulation and h the ellipsoidal height.

2.4 Refractivity

The refractive index of air n is commonly expressed with refractivity value N as shown in Eq. (9a). Under the ideal gas assumption, refraction of microwave signal follows the relationship in Eq. (9b), being a function of dry air and water vapour pressures (p_d and e) and temperature T (Bevis et al., 1994). Modifying the relationship for N to be a function of pressure p, temperature T and specific humidity q, which are the common atmospheric variables in NWP models, the new formula was proposed by Eresmaa and Järvinen (2006) leading to Eq. (9c).

$$\begin{array}{}\text{(9a)}& {\displaystyle }& {\displaystyle }n=\mathrm{1}+{\mathrm{10}}^{-\mathrm{6}}\cdot N,\text{(9b)}& {\displaystyle }& {\displaystyle }N=k_{\mathrm{1}}{\displaystyle \frac{p_{\mathrm{d}}}{T}}+k_{\mathrm{2}}{\displaystyle \frac{e}{T}}+k_{\mathrm{3}}{\displaystyle \frac{e}{T^{\mathrm{2}}}},\text{(9c)}& {\displaystyle }& {\displaystyle }N={\displaystyle \frac{k_{\mathrm{1}}p}{T}}+{\displaystyle \frac{(k_{\mathrm{2}}-k_{\mathrm{1}})qp}{(\mathit{ϵ}+(\mathrm{1}-\mathit{ϵ}\left)q\right)T}}+{\displaystyle \frac{k_{\mathrm{3}}qp}{(\mathit{ϵ}+(\mathrm{1}-\mathit{ϵ}\left)q\right)T^{\mathrm{2}}}},\end{array}$$

where $\mathit{ϵ}={\mathrm{M}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}/{\mathrm{M}}_{\text{dry}}\approx \mathrm{0.622}$ is the ratio of molar masses of water vapour and dry air. In LTT v1 (and v2), it is assumed that refractivity decreases exponentially with an increase in height between adjacent vertical levels of an atmospheric model. In the LTT ray tracer, coefficients are set in accordance with Bevis et al. (1994): k₁=77.60 K hPa⁻¹, k₂=70.4 K hPa⁻¹ and $k_{\mathrm{3}}=\mathrm{3.739}\times {\mathrm{10}}^{\mathrm{5}}$ K² hPa⁻¹. The RADIATE ray tracer uses coefficients by Rüeger (2002) with k₁=77.695 K hPa⁻¹, k₂=71.295 K hPa⁻¹ and $k_{\mathrm{3}}=\mathrm{3.754}\times {\mathrm{10}}^{\mathrm{5}}$ K² hPa⁻¹ (Landskron, 2018). The differences in slant delays due to different coefficients are rather small, about 1 mm at 5° elevation angle (Nafisi et al., 2012). Microwave refractivity increases slightly due to the presence of atmospheric particulate matter, such as hydrometeors, dust and biogenic aerosols (Solheim et al., 1999). However, in this study, their impact in slant delay computation is ignored.

3 Methodology

3.1 Improvements contained in LTT v2

Before discussing physical assumptions of the new algorithm, the key concepts underpinning the LTT v2 ray tracer are reviewed. The domain of the LTT v2 ray tracer solver (or the LTT domain grid) is a two-dimensional grid of refractivity N(i,j) and radius r(i,j). The first index i corresponds to the vertical model level of the NWP model and j the position of a vertical profile (more conveniently, a “column”). Thus, i and j are the vertical and horizontal coordinates in the LTT domain grid. Additionally, the station height, satellite position and top-of-atmosphere conditions are specified. This arrangement is illustrated in Fig. 1.

The LTT domain is constructed by interpolating refractivity values from the regular NWP model grid onto the propagation plane, which intersects the Earth at a great circle defined by the station location and azimuth. The interpolated columns are equally separated by angular separation Δψ, which is equal to a longitude increment Δλ of the regular Gaussian grid in the NWP model. As a default configuration, the LTT algorithm represents a propagation plane with 60 columns for NWP data resolution 640×1280 (Δλ=0.28125°) and 120 columns for 1280×2650 (Δλ=0.140625°), which can be changed, if needed. Hence, the default limit of the great circle distance from the receiver to the point where the signal enters the atmosphere is around 1900 km, and the LTT ray tracer is capable of computing tropospheric delay for as low as around 2° elevation angle signals. The contribution of slant delay from above the NWP model top is modelled with the Saastamoinen formula, and the ray is a straight line (Saastamoinen, 1972).

Figure 1 Illustration of the key concepts in the LTT algorithm. The RT integration starts at the ground station R and terminates at the radius r_S (red points). The satellite is located at S (black point). The signal path estimate is denoted by black solid curve inside the model domain, while propagation is linear above the model atmosphere (dashed black line). Atmospheric refractivity is computed at the LTT domain grid points (blue dots), spanning from the ground (orography, purple solid) to the highest model level of the NWP model (model top, grey solid). Vertical coordinates of the ray and the LTT domain are bonded to the idealized geoid (geoid, blue solid curve).

The domain's radial separation is inhomogeneous. The radius values of the ground station, the signal path and the domain grid are bonded to the selected coordinate system. In the LTT approach (both in LTT v1 and LTT v2), the shape of the geoid in the vicinity of the ground station is approximated with a sphere of a radius of curvature R_Eu, which is the Euler radius of curvature in the WGS-84 reference ellipsoid (Rodgers, 2000). This way, an origin of polar coordinates for Eq. (5) is the centre of curvature at the station position. The radial coordinates of the domain's grid points are computed separately for each vertical profile as the model level height plus R_Eu, as illustrated in Fig. 1. Model level heights are computed by solving the hydrostatic equation (Eq. 13), as explained in Sect. 3.2. Station height is provided to the LTT domain as the mean sea level height. It is pre-computed from ellipsoid height with Eq. (8b) using undulation products of Earth's geopotential model EGM2008 (Pavlis et al., 2012).

The RT equations are restored to exactly follow the theory of Rodgers (2000). Hence the horizontal refractivity gradient term has now been included in Eq. (5c). The along-plane horizontal gradient of a refractive index $\frac{\partial n(r,\mathit{\psi })}{\partial \mathit{\psi }}$ at any location (r,ψ) is computed under the assumption of refractivity changing exponentially in the vertical direction and linearly in the horizontal direction in the sub-grid scale, as follows:

$$\begin{array}{}\text{(10a)}& {\displaystyle \frac{\partial n(r,\mathit{\psi })}{\partial \mathit{\psi }}}={\mathrm{10}}^{-\mathrm{6}}\left({\widehat{N}}_{j+\mathrm{1}}\left(r\right)-{\widehat{N}}_j\left(r\right)\right)/\mathrm{\Delta }\mathit{\psi },\end{array}$$

$$\begin{array}{}\text{(10b)}& \begin{array}{rl}{\widehat{N}}_j\left(r\right)& =N(i,j)\\ & \cdot \text{exp}\left(-\text{ln}{\displaystyle \frac{N(i,j)}{N(i+\mathrm{1},j)}}\cdot {\displaystyle \frac{r-r(i,j)}{r(i+\mathrm{1},j)-r(i,j)}}\right),\end{array}\end{array}$$

where $N(.,.)$ and $r(.,.)$ are the known discretized values of refractivity and radius, such that point (r,ψ) is inside the grid cell $(i,j),(i+\mathrm{1},j),(i+\mathrm{1},j+\mathrm{1}),(i,j+\mathrm{1})$. ${\widehat{N}}_j\left(r\right)$ is the interpolated refractivity value for the vertical profile j at radius r. Δψ is the constant horizontal separation between columns. Horizontal gradient in the ray tracer has a minor impact on slant delay, altering the 5° elevation slant delay by less than 1 mm. Therefore, this term can be neglected in the operational setup to reduce computational cost.

Figure 2 Illustration of the horizontal gradient computation at point A located at (r,ψ). The adjacent refractivity values (N) needed in the formula are indicated in boldface. The intermediate values $\widehat{N}\left(r\right)$ are vertically interpolated refractivities, and the gradient is computed by Eq. (10a). The blue dots indicate the LTT domain grid.

The geometrical delay term in Eq. (6) is included in the new LTT v2 algorithm. Direct estimation of a path length integral in Eq. (6) is prone to large numerical errors. Therefore, the extension of delay due to bending (i.e. geometrical delay) is calculated based on simple geometrical assumptions, which are illustrated in Fig. 3a. The geometrical delay Δ is the difference between the straight and the real path, approximated with the following formula:

$$\begin{array}{}\text{(11)}& \mathrm{\Delta }=\underset{S}{\int }\mathrm{d}s-G\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{R}_{\text{ToA}}({\mathit{\psi }}_{E^{\prime }}-{\mathit{\psi }}_E)\cdot \mathit{\delta }\mathit{\theta }\cdot \cos {\mathit{\theta }}_E,\end{array}$$

where R_ToA is the radius at the top of atmosphere, ${\mathit{\psi }}_{E^{\prime }}$ and ψ_E the ψ coordinates of E^′ and E, δθ the difference between the geometrical and real zenith angles at the receiver, and θ_E the ray's zenith angle at E. In Fig. 3a, ΔTEE^′ is much larger than ΔREE^′, and the atmospheric part of the signal (RE) is assumed to be a straight line. Hence, the geometrical delay is $RE+ET-RT$. Since the atmospheric part of the signal (RE) is assumed to be a straight line (Fig. 3a), the geometrical delay can be refined by calculating the length of the RE curve conventionally. This refinement, however, is not implemented in the LTT v2.

The impact of geometrical delay Δ is seen at 20° elevation and below, reaching 50 cm at 5° elevation. At the same time, the value of Δ can vary significantly at low elevations, as this approximation is sensitive to the station height and orography. In Fig. 3b, Δ less than 10 cm for 5° elevation occurs only for stations POL2 and UNSA (IGS, 2024b, c), which both have steep orography in some azimuth directions.

Figure 3 (a) An idealized setup for geometrical delay calculation. The solid bold black curve is the signal path from a transmitter T to the receiver R, entering the atmosphere at E. A straight signal path would reach the atmosphere at E^′. The geometrical delay is $RE+ET-RT$. (b) Two-dimensional histogram of geometrical delays. The geometrical delays are computed for 256 stations at 72 azimuth angles and 85 zenith angles on 5 December 2016. The total number of values is around 3.8×10⁷. Colour bar indicates the number of occurrences of geometrical delay Δ in each bin (40 equal zenith angle bins and 500 equal delay bins).

In LTT v2, the initial value problem of finding the zenith angle of a ray when it reaches the receiver (or initial zenith angle, θ₀) is solved iteratively. Similarly to Eq. (7), the correction of initial zenith angle estimate ${\mathit{\theta }}_{\mathrm{0}}^n$ yields the next value ${\mathit{\theta }}_{\mathrm{0}}^{n+\mathrm{1}}$:

$$\begin{array}{}\text{(12)}& {\mathit{\theta }}_{\mathrm{0}}^{n+\mathrm{1}}={\mathit{\theta }}_{\mathrm{0}}^n-\mathit{\alpha }\mathrm{\Delta }{\mathit{\psi }}^n,\end{array}$$

where α is an arbitrary scaling parameter, $\mathrm{0}<\mathit{\alpha }\le \mathrm{1}$ and Δψⁿ the deviation of angle ψ from the target value at the nth iteration.

Figure 4 provides an example of an initial value iterative search in a selected LTT domain (see Fig. 1) for a variety of satellite zenith angles. The initial-end conditions diagram is in the left plot, and produced slant delays are on the right, computed for the first 10 iterations by Eq. (12). As seen in Fig. 4 (left), the initial zenith angle adjustment is increasingly larger the larger satellite zenith angle is. Curiously, the convergence of the deviation Δψⁿ is almost perfectly exponential, which implies that the uncertainty of initial θ₀ is proportional to the error of the end condition. Also, after many iterations, the deviation Δψ weakly depends on satellite zenith angle, in this example, reaching around 10⁻⁵° by iteration #9 for all zenith angle values. The slant delay value is adjusted throughout the iterative process and converges after 5–6 iterations, as shown in Fig. 4 (right). In the LTT v2 algorithm, the number of iterations to produce the final delay value is fixed to 8 and $\mathit{\alpha }=\mathrm{1}/\sqrt{\mathrm{2}}$; no threshold criteria on convergence is applied.

Figure 4 Difference between the geometrical zenith angle of the satellite (θ_G) and the ray zenith angle at the receiver (θ₀) versus the angular separation between the endpoint of the ray (ψ_f) and the satellite location (ψ_sat) for the first 10 iterations. The values in grey shades are plotted for different satellite zenith angles, ranging from 10 to 86° with 1° steps. The coloured lines indicate the convergence trajectory for selected zenith angles (a). The resulting delay values at 5° elevation produced at each iteration (b). Scaling parameter α is $\mathrm{1}/\sqrt{\mathrm{2}}$. Data are produced using the LTT v2 algorithm for GNSS station JIUFENG (JFNG) on 5 December 2016 at 07:00 UTC for azimuth angle 45°.

3.2 Data input from the OpenIFS weather model

In this paper, all three RT algorithms (LTT v1, LTT v2 and RADIATE) are supplied with numerical weather data of the OpenIFS weather prediction model, which is the atmospheric component of the Integrated Forecasting System (IFS) of the European Centre for Medium-Range Weather Forecasts (ECMWF) (ECMWF, 2017a). OpenIFS is a portable version of IFS intended for academic use. It produces forecasts with a skill identical to the full IFS. However, data assimilation is not included in OpenIFS, and the initial conditions must be supplied separately with the optimal ECMWF state estimates (Lean et al., 2021; Rabier et al., 2000). The version of OpenIFS used here is 43r3v1, which was part of the operational forecasting system at ECMWF from July 2017 to June 2018.

The forecast setup is as follows. 12 h forecasts are produced using initial states at 00:00 and 12:00 UTC for each day during December 2016. The model's time integration proceeds with 10 min step, and the output interval of the state is 1 h, resulting in 744 atmospheric states per month. The output contains the following atmospheric variables: surface geopotential, logarithm of surface pressure, and temperature and relative humidity at model levels. The atmospheric simulations are performed at a resolution of T_L 1279, corresponding roughly to 18 km horizontal resolution at the Equator, at 137 vertical levels (+ a surface level) (ECMWF, 2022). Administration of the atmospheric simulations and output data is facilitated by the OpenEPS workflow manager (Ollinaho et al., 2021; Ollinaho, 2020).

The weather data are output at model levels, the heights of which are computed by the LTT program. In OpenIFS (ECMWF, 2017b), the construction of model-level orthometric heights h starts by defining the height of the surface level h_N+1 using the surface geopotential. The lowest model level h_N follows the model surface 10 m above it. The heights of model levels h_i are computed from the lowest model level N to the top of the atmosphere (level 1), as follows:

$$\begin{array}{}\text{(13a)}& {\displaystyle }& {\displaystyle }{h}_{N+\mathrm{1}}={\displaystyle \frac{z}{g}},\text{(13b)}& {\displaystyle }& {\displaystyle }{h}_N=h_{N+\mathrm{1}}+\mathrm{10}\mathrm{m},\text{(13c)}& {\displaystyle }& {\displaystyle }{\mathit{h}}_i={\mathit{h}}_{i+\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{g}}{R}_{\mathrm{d}}{T}_{\mathrm{1}/\mathrm{2}}^{*}\ln \left({\displaystyle \frac{p_{i+\mathrm{1}}}{{\mathit{p}}_i}}\right),\text{(13d)}& {\displaystyle }& {\displaystyle }{\mathit{p}}_i={\displaystyle \frac{{\mathit{a}}_i+{\mathit{a}}_{i+\mathrm{1}}}{\mathrm{2}}}+{\displaystyle \frac{{\mathit{b}}_i+{\mathit{b}}_{i+\mathrm{1}}}{\mathrm{2}}}{p}_{\mathrm{s}},\end{array}$$

where z is the surface geopotential (in m² s⁻²), g the constant mean sea level gravity (9.80263 m s⁻²), R_d the dry air constant (287.04 $\mathrm{J}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{kg}}^{-\mathrm{1}}$), $T_{\mathrm{1}/\mathrm{2}}^{*}$ (K) the half-layer virtual temperature, and p_i (Pa) the half-level pressure at level i computed from the surface pressure p_s using the sigma coefficients a and b. Before vertical column construction, the atmospheric variables are bi-linearly interpolated from the horizontal grid.

3.3 OpenIFS input data to the RADIATE ray tracer

To deploy RADIATE with OpenIFS data, conversion to a format compliant with RADIATE's rigid requirements is needed. RADIATE requires weather data in a hardcoded format, mandating a spatial representation of data in a regular grid from 90 to −90° in latitude and 0 to 359° in longitude. The horizontal resolution of weather data is scaled to comply with default 1° × 1° of RADIATE input. Despite being technically feasible to interpolate weather data with finer resolution, Hofmeister et al. (2015) found no strong evidence to suggest changing the default resolution. The vertical resolution must adhere to the 25 pressure levels used in ECMWF operational forecasts starting from 1 hPa downwards. The input file for RADIATE is a plain text file, featuring a hardcoded section containing information about resolution, extent, the list of pressure levels, the number of parameters and the date of weather data, followed by the actual data – geopotential height Z (m), specific humidity q and temperature T (K).

3.4 Production of slant delays

Information about the signal tropospheric slant delay is required in GNSS data analysis. To employ the modelled delays in GNSS precise orbit determination (POD), we utilize the direct slant delay coupling approach of Navarro Trastoy et al. (2022), where the delay estimates are provided to the POD solver as is, omitting the use of low-order parameterization by, for example, Vienna mapping functions (Landskron and Böhm, 2018). The three ray tracers (LTT v1, LTT v2 and RADIATE) are configured to produce these products.

The slant delay products are as follows: slant delay values are computed in a regular azimuth–elevation grid for each station. These are called skyviews. Their resolution is 5° × 1° in azimuth and elevation, respectively. Technically, the ray tracer algorithm is launched multiple times with a virtual GNSS satellite located at the grid points of a skyview. The delay estimate for an individual GNSS observation is obtained by spatially interpolating the skyview value for the azimuth–elevation angles, and temporally between the output epochs.

The RADIATE skyviews are produced as follows. The user is required to create an AZEL (azimuth–elevation) file, detailing the directions of observations. We are interested in a regular grid of azimuths and elevations; thus, we use the option of creating a uniform AZEL file. The file format is specified in a separate file, where elevations and the number of uniformly distributed azimuths (set at 72 in our case) are defined.

The frequency of the slant delay products is four times per day (every 6 h). The LTT v1 and LTT v2 algorithms digest slant delays with any time resolution, e.g. the hourly OpenIFS model output. However, we are limited by the design of the RADIATE software, which allows the desired number of output epochs per day to be set to a maximum of 4. Hence, a 6 h interval is selected to ensure a fair comparison.

3.5 LTT v2 as a software

We assembled the LTT v2 model into an openly available Fortran-language software with a modular structure. The weather data are uploaded to LTT v2 from a GRIB (GRIdded Binary) format file using the ecCodes (Najm, 2021) routines. Other required inputs are the station list containing a four-character identifier, position and height for each station, as well as the IFS sigma coefficients defining the model levels. The LTT v2 program is configurable via the so-called namelist files, which specify the format of delay products and station height definition.

The core of the new ray tracer is the refined LTT operator routine, which maps the LTT domain, the satellite and the station coordinates (Fig. 1) to the slant delay value via solving the Eq. (5). They are integrated as follows. Each NWP model level is divided into equal differential steps, and the signal path increments are computed in each step via a 4-degree Runge–Kutta scheme. The number of steps per level and the number of initial zenith angle iterations are adjustable in the program. The baseline computing setup has 10 steps per level and 8 initial zenith angle iterations. A single baseline LTT RT execution costs 1–10 CPU ms, depending on the station height and satellite elevation, and the 72×85 skyview takes 40 CPU s to compute.

The LTT v2 program is designed to produce slant delay skyviews; however, other modelling strategies are possible. A viable option is to perform computation per observation, given the location, time, azimuth and elevation of the event. Thus, no interpolation is needed from a skyview value, but time interpolation remains. The main issue of this strategy is increased computational overhead due to the construction of the LTT domain for each observation. The GNSS network produces tens of millions of observations per day, which makes this approach unrealistic. However, geodetic systems with higher accuracy demands and fewer observations can benefit from per-observation slant delay computation.

4 GNSS processing and performance metrics

Assessing the quality of ray tracing in a neutral gas atmosphere is complicated due to the lack of an absolute reference. The approach here is to apply a software to process GNSS constellations and ground station networks. This allows validation of the different ray tracers by comparing resulting GNSS products with all other inputs remaining unchanged. GROOPS (Gravity Recovery Object Oriented Programming System) is used here to produce the GNSS products. GROOPS (Mayer-Guerr et al., 2021) is an open-source software toolkit for gravity recovery that also includes a package for processing GNSS constellations and ground station networks. In GROOPS, the “GNSSProcessing” program is used, which solves orbits and ground network station locations simultaneously using least-squares adjustment.

The products of GNSS data processing are position and velocity states, as well as clock parameters, ionospheric electron content, integer ambiguities and phase biases. Since there is no absolute reference, the processing is evaluated by measuring the stability of the solution. Discussions by Strasser et al. (2019) and Zus et al. (2021) suggest that among all the most informative are orbit and ground station positions. Hence, we choose two metrics to validate ray tracing models: midnight discontinuities of satellite orbits for consecutive days, and precise point positioning error of ground stations. GNSS products are computed daily in observation batches of 24 h. The daily products are thus independent of one another. The midnight discontinuity (MD) means the distance between the last point of one day's orbit and the first point of the following day (Navarro Trastoy et al., 2022; Massarweh et al., 2021). For the precise point positioning (PPP), the daily station positions of the set of 205 stations from the International GNSS Service (IGS) are computed and compared with a long-time-series reference position from IGS14 (Rebischung and Schmid, 2016). We compute root mean square (RMS) and standard deviation (SD) of the differences to the reference positions. RMS indicates how close the PPP solutions are to the reference IGS14 solutions, whereas SD indicates the consistency of the PPP solutions.

The GNSS POD experiments are performed producing both metrics using the slant delay estimates from three ray tracers: RADIATE (Landskron, 2018), LTT v1 (Eresmaa et al., 2021) and LTT v2. Thirty days of GNSS data are processed in December 2016.

5 Results

5.1 Delay estimation and GNSS processing with LTT

Our processing setup using LTT-based ray tracers comprises production of slant delay estimates and creation of GNSS orbit solutions and ground station positions. The new version of LTT (v2) is compared to the previous version (v1). First, the whole monthly products of tropospheric slant delays are compared in a statistical representation. Figure 5 summarizes the differences between these models, showing the residual statistics for the slant delays in skyviews produced during 1 month for the set of 256 ground stations. The red areas on the histogram correspond to difference values of the overwhelming majority of slant delays, while the green and blue areas show the extent of rare extreme discrepancies between the models. Since the delay values are computed for a variety of weather and geographical conditions, Fig. 5 can be seen as a probability analysis of the gross errors in slant delay models, in analogy to Monte Carlo simulations.

Figure 5 Two-dimensional histogram of ray-traced slant delay difference: LTT v2 minus LTT v1 estimates. The data sample comprises around 1.9×10⁸ slant delays (124 epochs at 256 stations and 72 azimuth angles by 85 zenith angles). Colour bar indicates the number of occurrences of delay differences in each bin. There are 85 equal zenith angle bins and 500 equal delay bins.

The difference between the two LTT model versions is asymmetric. As seen in Fig. 5, slant delays of LTT v2 are typically smaller than those produced with LTT v1 at high elevations, whereas below about 75° zenith angle the opposite is true. Two effects might explain this behaviour: the refinement of the initial conditions in LTT v2 leads to a preferably smaller starting zenith angle θ₀ (Fig. 4), thus decreasing the delay, whereas the addition of the geometrical delay term to LTT v2 increases the delay at low elevations (Fig. 3b).

The second part of analysis is obtaining metrics for GNSS network solutions informed by each tropospheric delay model. The first metric used here is the midnight discontinuity (MD) for GPS satellites, specifically, the average positional MD of all satellite orbits (Δ) and the standard deviation related to this average (σ_Δ):

$$\begin{array}{}\text{(14a)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }={\displaystyle \frac{\mathrm{1}}{N_{\text{GPS}}}}\sum _{i=\mathrm{1}}^{N_{\text{GPS}}}\sqrt{\mathit{\delta }{r}_{\mathrm{a},i}^{\mathrm{2}}+\mathit{\delta }{r}_{\mathrm{c},i}^{\mathrm{2}}+\mathit{\delta }{r}_{\mathrm{r},i}^{\mathrm{2}}},\text{(14b)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{\Delta }}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N_{\text{GPS}}}}\sum _{i=\mathrm{1}}^{N_{\text{GPS}}}{\left(\sqrt{\mathit{\delta }{r}_{\mathrm{a},i}^{\mathrm{2}}+\mathit{\delta }{r}_{\mathrm{c},i}^{\mathrm{2}}+\mathit{\delta }{r}_{\mathrm{r},i}^{\mathrm{2}}}-\mathrm{\Delta }\right)}^{\mathrm{2}}},\end{array}$$

where N_GPS is the number of processed satellites in a GPS constellation (usually equal to 32), and $(\mathit{\delta }{r}_{\mathrm{a}},\mathit{\delta }{r}_{\mathrm{c}},\mathit{\delta }{r}_{\mathrm{r}})$ indicate the position difference in along-track, cross-track and radial direction between the last point of one day's orbit and initial point of the following day. The Δ and σ_Δ are computed daily for December 2016, resulting in 30 pairs of values. These values are shown in Fig. 6, with Δ on the left and σ_Δ on the right.

Figure 6 Comparison of daily midnight discontinuities (Δ, a) and their standard deviations (σ_Δ, b) with tropospheric slant delay estimates from LTT v1 and LTT v2 for December 2016. Points are the mean values over the entire satellite constellation, computed with Eq. (14). The star indicates the mean value, also reported in Table 1.

The second GNSS processing metric is related to the precise point positioning (PPP) of the ground stations involved in the processing. The position is defined by north, east and height components. We analyse the residuals x of daily positions of the set of 205 stations for 30 d in December 2016.

$$\begin{array}{}\text{(15)}& \mathit{x}=\mathit{r}-{\mathit{r}}_{\mathrm{IGS}\mathrm{14}},\end{array}$$

where r_IGS14 is the reference daily position and r the result of the conducted processing. The obtained positions are network solutions which comply with the net-zero rotation condition. In reality, the resulting reference frame is different each day due to noise and natural variations; however, position values are not transformed from this fluctuating frame towards the IGS14 reference frame. From the monthly time series of x, the metrics are computed for each ground station as follows:

$$\begin{array}{}\text{(16a)}& {\displaystyle }& {\displaystyle }{\mathit{x}}_{\text{RMS}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{30}}}\sum _{d=\mathrm{1}}^{\mathrm{30}}{\mathit{x}}_d^{\mathrm{2}}},\text{(16b)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_x=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{30}}}\sum _{d=\mathrm{1}}^{\mathrm{30}}({\mathit{x}}_d-\stackrel{\mathrm{‾}}{\mathit{x}}{)}^{\mathrm{2}}},\end{array}$$

where $\stackrel{\mathrm{‾}}{\mathit{x}}$ is the mean residual, x_RMS the root mean square and σ_x the standard deviation of residual in east, north and height components. Squaring of vectors is performed element-wise. The distributions of x_RMS and σ_x among different stations are shown in Figs. A1 and A2 of the Appendix, respectively. The mean values over all stations are shown in Table 1.

Table 1 The metrics for GNSS processing supplied with tropospheric slant delays from LTT v2 and LTT v1 ray tracers. Percentages show the improvement (positive) or deterioration (negative) of the LTT v2 algorithm compared to LTT v1. Values are computed as follows: a LTT v1 value minus LTT v2 value divided by the LTT v1 value; average x_RMS and σ_x are transformed into Euclidean distance.

Table 1 summarizes values for both MD and PPP metrics and their respective standard deviations over the orbit days (for MD) and over the set of ground stations (for PPP). LTT v1 and LTT v2 demonstrate comparable skill, with v1 performing better at precise point positioning and v2 better at orbit determination. Standard deviation, which is indicative of consistency of processing, is better when using the LTT v2 model. As seen from Fig. 6, the MD values vary greatly on a day-to-day basis; thus, the improvement from LTT v1 to v2 is seen as statistically insignificant. The same conclusions can be drawn from Figs. A1 and A2. This implies that modifications in ray tracing models from the old to the new version have a very minor or even negligible effect.

5.2 Results of using the special RADIATE setup

The VieVS Ray Tracer (RADIATE) is deployed with the same processing setup as the LTT model. Figure 7 shows the discrepancies between LTT v2 and RADIATE slant delay estimates accumulated over 1 month. According to Fig. 7, slant delays of LTT v2 deviate much more from RADIATE-based values than from LTT v1. For 1° zenith angle delays, the mean delay difference is −0.685 cm, and the standard deviation is ±0.645 cm. The corresponding values for discrepancies between LTT v2 and v1 are $-\mathrm{0.583}\pm \mathrm{0.287}$ cm (Fig. 5). The LTT–RADIATE delay difference is close to a normal distribution with a preference for RADIATE delays being slightly larger than those of LTT v2 at most zenith angle range. At below 10° elevation, however, LTT v2 usually produces larger delays.

Figure 7 Two-dimensional histogram of ray-traced slant delay difference: LTT v2 minus RADIATE estimates. The data sample comprises around 1.9×10⁸ slant delays (124 epochs at 256 stations and 72 azimuth angles by 85 zenith angles). Colour bar indicates the number of occurrences of delay differences in each bin. There are 85 equal zenith angle bins and 500 equal delay bins.

The midnight discontinuity metric shows that LTT v2 is an overall improvement over RADIATE, showing more consistent and accurate orbital products for most of daily solutions.

Figure 8 Comparison of daily midnight discontinuities (Δ, a) and their standard deviations (σ_Δ, b) with tropospheric slant delay estimates from RADIATE and LTT v2 for December 2016. Points are the mean values over the entire satellite constellation, computed with Eq. (14). The star indicates the mean value, also reported in Table 2.

The station precise point positioning residuals are generally smaller when using the LTT algorithms compared to RADIATE. In the RADIATE experiment, the height estimates are inconsistent for some stations, leading to high average RMS and σ values. These occasional high values do not appear in Figs. A1 and A2.

Table 2 The metrics for GNSS processing supplied with tropospheric slant delays from LTT v2 and RADIATE ray tracers. Percentages show the improvement (positive) or deterioration (negative) of the LTT v2 algorithm compared to RADIATE. Values are computed as follows: a RADIATE value minus LTT v2 value divided by the RADIATE value; average x_RMS and σ_x are transformed into Euclidean distance.

6 Discussion

The three RT algorithms are principally similar. They employ a two-dimensional regime of ray propagation, discrete path bending, and continuous refractivity change being linear in horizontal and exponential in vertical directions. The weather forecast data are the same for these RT algorithms. The delay estimates by the three algorithms are provided in the same format and resolution for the verification procedure of obtaining global navigation solutions. However, verification via the GNSS processing demonstrates different performances of these three RT models. Many details seem to be important.

First, the atmospheric state is input differently in RADIATE and LTT. RADIATE employs global NWP data at 1° × 1° resolution at 25 hardcoded atmospheric pressure levels, which are interpolated from the full OpenIFS model resolution definition. LTT v1 and LTT v2, in contrast, access the original NWP data at the native resolution (around 0.14° × 0.14°) at 137 model levels. Downscaling both the horizontal and vertical resolutions truncates small-scale atmospheric phenomena, leading to less informative delay in RADIATE than in LTT. It is worth noticing that the time resolution of slant delays is 6 h, and time-linear interpolation is applied. This hinders the influence of time evolution and movement of weather patterns both in RADIATE and LTT delay estimates.

Second, the RT domain grid is created differently in the three RT algorithms. Both LTT versions employ a sequence of model level heights to build up the grid vertically, while RADIATE has a different setup. LTT v1 constructs full-pressure model level heights, while LTT v2 makes use of half-level pressure values. In addition, the gravity acceleration g in the hydrostatic equation (Eq. 13) is changed to a constant value in LTT v2. While not physically correct, constant g with respect to height is technically consistent with the OpenIFS model. It is important to note that the LTT v1 model is a modification of the slant delay operator by Eresmaa et al. (2008). The latter was intended as an extension of the High-Resolution Limited Area Model (HIRLAM), whose model level definition is different from that of OpenIFS.

Although the main physical assumptions in the three RT algorithms are the same, the details in the implementation are different. In LTT v2 we include plausible physical improvements. Attribution of impacts of different improvements is not performed in full since they influence each other in the GNSS processing. Considering the importance of correctly utilizing the weather model data, improving the physical assumptions is a somewhat minor role. We emphasize that LTT v2 has a principal advantage over LTT v1 and RADIATE models since it is directly compatible with the Integrated Forecasting System (IFS) at any present or future resolution.

7 Conclusions

We have developed a signal delay model which is closely tied to the OpenIFS numerical weather model ensuring nearly lossless representation of the atmospheric state for solving a signal propagation. The LTT v2 algorithm has a more robust mathematical foundation of signal propagation than the previous version (LTT v1). The signal delay model validation via GNSS processing concludes that the LTT-based approach outperforms the RADIATE software. We attribute this mainly to the much larger amount and lossless utilization of weather model data as input to the delay modelling. The two LTT models, versions 1 and 2, on the other hand, have fairly similar skill, suggesting that the physical improvements of ray tracing have less impact on slant delay quality than lossless adaptation of numerical weather model data as input to the signal delay modelling. This supports the conclusion that piecewise-linear and Rodger's implementation of Eikonal ray tracing lead in practice to very similar performance. In comparison with earlier developments in the atmospheric signal delay of GNSS and other geodetic systems, the main innovation we introduce is the relaxation of all assumptions regarding the use of weather model data. It is now possible to fully enjoy the increasing quality and skill of weather model data.

From the perspective of the future of RT applications, the newly developed LTT v2 is well positioned for the emerging next-generation high-resolution global and regional models (nextGEMS, 2024; ECMWF, 2024a; Hohenegger et al., 2023; Rackow et al., 2025). Since LTT v2 can ingest atmospheric model output data at any horizontal and vertical resolution, only a minimal amount of post-processing of the model fields is required. This means that LTT v2 is ready to take full advantage of the fine-scale details of meteorological information provided by the next-generation high-resolution models.

Appendix A: Statistics of station position residuals

The process of GNSS ground station positioning yields one residual estimate for one station per day, producing 205 time series with 6150 data points in total in our experiment. We performed three experiments: LTT v1, LTT v2 and RADIATE. The distribution of the root mean square and the standard deviation of the residual time series are shown in this section.

A1 Root mean square

Figure A1 Distribution of root mean square of position residuals x_RMS for IGS 205 ground stations. In the three-by-three table of figures, columns distinguish between the components of x_RMS, and rows distinguish between different slant delay products.

A2 Standard deviation

Figure A2 Distribution of standard deviation of position residuals σ_x for IGS 205 ground stations. In the three-by-three table of figures, columns distinguish between the components of σ_x, and rows distinguish between different slant delay products.