GNSS signal ray-tracing algorithm for the simulation of satellite-to-satellite excess phase in the neutral atmosphere

Traditionally, GNSS space-based and ground-based estimates of tropospheric conditions are performed separately. It leads to limitations in the horizontal (e.g., a single space-based radio occultation proﬁle covers a 300km slice of the troposphere) and vertical resolution (e.g., ground-based estimates of troposphere conditions have spacing equal to stations’ distribution) of the tropospheric products. The ﬁrst stage to achieve an integrated model is to create an effective 3D ray-tracing algorithm for the satellite-to-satellite (radio occultation) path reconstruction. We verify the consistency of the simulated data with the RO observations from the Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC-1) Data Analysis and Archive Center (CDAAC) in terms of excess phase and bending angle. The results show that our solution provides an effective RO excess phase, with a relative error varying from 35% at the height of 25–30km (1.0–1.5 m) to 0.5% at heights 5–10km (0.1–1m) and 14 to 2% at heights below 5km (2–14m). The bending angle retrieval on simulated data attained for high-resolution ray-tracing, bias lower than 2% with respect to the observed bending angle. The optimal solution takes about 1s for one transmitter–receiver pair with a tangent point below 5km altitude. The high-resolution processing solution takes 3 times longer.


Introduction
The effective use of observational information from satellite data is an important prerequisite for global weather prediction.As stated by Bauer et al. (2015) in their fundamental paper on the future of weather prediction, the provision of physically consistent initial conditions is one of the key challenges to attain a 1 km resolution global weather forecast.This request is directly linked to the reliability of observations and thus of interest to the meteorological remote sensing community and GNSS researchers in particular.
In this paper, we describe the development of a 3D raytracing (RT) GNSS Geometric Optics simulator, which in follow-up studies will become part of the 3D integrated GNSS tomography model.Traditionally, in GNSS tomography, the 3D refractivity field is estimated from integrated ground-based GNSS measurements, using inverse Radon transform (Flores et al. 2000).The technique's reliability heavily depends on the density of scanning rays (path delays) and the intersection angles between these rays (Kak and Slaney 2001).Hence, for ground-based GNSS tomography, the transmitter-receiver geometry does not favor the high vertical resolution and leads to solution instability (Trzcina and Rohm 2019;Liu et al. 2019).Conversely, Moeller et al. (2019) has shown that a tomographic "peeling" approach with space-based observations can significantly increase the vertical resolution if applied to radio occultation measurements.The scanning rays, in our solution, will come from ground-and space-based GNSS observations.This approach can benefit from the high density of ground-based observations and detailed vertical structures to capture small-scale atmospheric features.The growing density of ground-based networks (e.g., E-GVAP network contains now 3500 receivers) and an increasing number of satellites (e.g., 121 Spire Global radio occultations satellites producing 10000 measurements per day) make this research well aligned with current development in the satellite observation segment.
The simultaneous use of space-based and ground-based observations in the tomography concept requires a unified ray-tracing solution, which is demonstrated in this manuscript.As the ground-based ray-tracing solution has been studied to a great extent (e.g., Kačmařík et al. 2017;Hobiger et al. 2008;Nievinski and Santos 2010) and the signal propagation path above 15 • elevation angle is well represented as a straight line connecting transmitter and receiver (Möller and Landskron 2019) in this study, we focus only on space-based ray-tracing Høeg et al. (1995).Conversely, in RO observations, the measurable effect is directly linked to the bending of the signal in the atmosphere measured as an excess path, reaching 100 m at 5 km height (Steiner and Kirchengast 2005).
There are two major formulations used to simulate spacebased GNSS signal propagation: geometric optics (GO) and wave optics (WO).GO assumes that the signal frequency is very high and there are no sharp changes in the propagation environment Nievinski and Santos (2010).An excellent study employing tracing GNSS signal was presented by Beyerle et al. (2002).However, the focus is on reflected signals, and the majority of the study (3769 events) is based on GO 2D ray-tracing for both phase and amplitude to form differences between observed and ray-traced electromagnetic fields.Similarly, Healy (2001) developed the geometric optics simulator to verify the impact of spherical symmetric approximation and horizontal gradients, on the bending angle and tangent point height.Further studies also employed geometric optics forward operator (Syndergaard et al. 2005) to assimilate GNSS RO refractivity.In climate change monitoring simulation studies, using RO constellation focused on Upper Troposphere Lower Stratosphere also 2D RT (Foelsche et al. 2008) was applied.The GO raytracing was also used to assess the impact of severe weather on ground-based and space-based observations.In the study by Norman et al. (2014), researchers demonstrate that, if the strong horizontal gradients are present (4 ppm/km on the ground, overall increase of absolute value from 320 to 350 ppm/km), the tangent point for the lowermost ray will deviate up to 5 km and the excess phase will be 2 m larger.It scales down exponentially with height but could be a major impact on troposphere tomography as the signal may cross different model elements.
In more complex conditions, when the GNSS signal is subject to an amplitude variation and diffraction on the strong vertical gradient, the assumption of geometric optics may be violated; therefore, wave optics techniques need to be adopted.The signal is represented as a complex field propagating through the variable refraction landscape.The WO is used to effectively handle interfering signals in the presence of irregular water vapor structures; therefore, it is less vulnerable to multipath (Gorbunov 2002).The WO technique could be used to simulate spikes introduced by higher order terms of refractivity leading to a better understanding of water and ice impact on the radio occultation retrievals Hordyniec (2018).Mixed GO and WO methods were studied in one of the latest works on simulating RO with the RT technique by Gorbunov et al. (2019).The authors review the latest techniques of RO assimilation into numerical weather models (NWM).According to their results, both bias and standard deviation of bending angle calculated as a difference between simulations and Spire Global data are lower for 3D GO RT simulations than 1D RT simulations by around 0.03 rad in the bottom 5 km of the troposphere.These discrepancies are more visible in the tropics than in other locations on the globe.
In summary, the 2D GO ray-tracing is a well-established solution to simulate the propagation of the GNSS signal and the out-of-plane horizontal gradients.The spherical symmetry assumption will reduce the quality of the simulation; hence, 3D GO ray-tracing should be employed for severe weather signals.In more challenging conditions, with signal amplitude variation and diffraction on the strong vertical gradient, the wave optics technique needs to be employed; however, it is computationally expensive and difficult to apply in tomography models, as there is no "signal path" that could feed into inverse Radon transform.Therefore, the 3D GO RT method will be applied in this study to simulate the signal, whereas both GO and WO were used to retrieve bending angle.Geometrical optics (GO) ray-tracing simulations used in this study are not necessarily expected to produce results close to observation.The influence of unmodeled ionosphere in the upper part of the atmosphere, multipath or ducting conditions in lower part of the troposphere can be significant (Jensen et al. 2004); however, as shown in Beyerle et al. (2002), GO approximation represents well the majority of troposphere conditions (in Beyerle et al. (2002) only 14 out of 3783 cases required wave optics (WO) ray-tracing approach).
Our ray-tracer is tested on 18 different settings scenarios and is an extension of the work of Möller (2017).We selected 10 ROs as the test object from work by Healy (2012).Simulation results were compared with CDAAC's atmPhs files for excess phase, bending angle; we also tested the alignment of CDAAC's refractivity and vertical gradient of this parameter, with ERA5 ray-traced values.
The paper is structured as follows: after this Introduction Sect. 1, the RO data and the weather models are described in Sect. 2. 3D RT with the tested settings and the validation process are presented in Sect.3. Section 4 describes the simulation and validation results.The discussion and conclusions of the study are presented in Sects.5 and 6, respectively.

Radio occultation data
The RO data used in this study come from CDAAC (2022) data center and contain data from the COSMIC-1 mission managed by the National Oceanic and Atmospheric Administration (NOAA).The COSMIC-1 constellation consisted of six low Earth orbit (LEO) micro-satellites launched in 2006.For the conducted analysis in particular the atmPrf and atmPhs files were used.The first one stores information about full-resolution profiles of physical parameters such as dry pressure, dry temperature, refractivity, bending angle, and impact parameter.The second contains the atmospheric excess phase by setting or rising GNSS satellite as measured on the L1, L2, and LC (ionosphere-corrected) bands.In the first case, the data consist of ordered consecutive epochs of satellite-receiver pair coordinates sampled with a frequency of 50 Hz (CDAAC 2021), while the second is ordered by geometric or impact height variable.
In this study, we selected 10 RO cases reported by Healy (2012).The author divides ROs into four categories.The conditions are based on the vertical gradient of refractivitythe rate of change of refractivity with respect to changes in corresponding altitude (dN /dz).For this work, we slightly modified them for application in the ERA5 model, which is more smooth than observations.The first category means standard refractive conditions with a vertical gradient of refractivity at a selected height smaller than one-third of the gradient required for ducting (dN /dz ≤ (−)60 km −1 ).The second and third category profiles are getting more complex due to the increasing number of perturbations ((−)60 < dN /dz < (−)100 km −1 and (−)100 < dN /dz ≤ (−)140 km −1 ).The fourth category refers to the conditions where the maximum refractivity gradient is dN /dz > (−)140 km −1 .From each category, we selected at least one RO profile.All selected cases are described in Table 1.
To clearly present the variability of selected RO profiles, the differences in the vertical gradient of refractivity for each category were plotted.These values are shown in Fig. 1 as a function of height.The refractivity values were derived from atmPrf files.

ERA5 model reanalysis
As an input to the ray-tracing algorithm, we used in this study the fifth-generation European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalysis (ERA5) of the global climate.ERA5 delivers the most upto-date global reanalysis with data from 1950 to the present.The model assimilates weather data from observations (e.g., satellite imagery) and measurements (e.g., weather stations) with relatively high quality.In some cases even outperforms satellite-based precipitation products (Xu et al. 2022).More details of ERA5 can be found in C3S (2017).In this study, data from ERA5 were used for selected cases at epochs when the ROs occurred in a 1-hour time span with a horizontal resolution of 0.25 • .The area for which the data were analyzed covered the entire globe at all 37 available pressure levels (from 1 to 1000 hPa).As the maximum height of the model is limited by the 1 hPa pressure level, which corresponds to an altitude of approximately 80 km, for heights above 80 km the U.S. Standard Atmosphere 1976 (US76) was used.The transition between models is discontinuous due to low values of refractivity (< 0.1) at that height.US76 is a static Earth atmospheric model of density, temperature, and viscosity.It divides the atmosphere into layers with linearly distributed absolute temperature and geopotential altitude.More about US76 can be found in USSA (1976).

3D ray-tracing algorithm
The 3D ray-tracing algorithm is based on the reconstruction of the signal path within the grid model (3D refractivity field).For our studies, an extension of the original 3D ray-tracing module from ATom (Atmospheric Tomography) GNSS software package (Möller 2017) is used.The necessary modifications were inspired by Javaherian et al. (2020).The first step of ray-tracing in our approach (shooting technique) is based on the prior knowledge of transmitter and receiver positions, as well as information about weather parameters and thus, refraction at nodal points of the grid.The transmitter and receiver are the GPS and LEO satellites, respectively.In the current version, coordinates expressed in earth-centered inertial (ECI) coordinate system (Herman 1995) are delivered from atmPhs files and recalculated (Petit et al. 2010) to earth-centered-earth-fixed (ECEF) coordinate system.The total refractivity values N are calculated from pressure p, temperature T , and water vapor partial pressure e derived from ERA5.grib files (CDS 2022) using the N-apriori with refractivity coefficients k 1 = 77.689K/hPa, k 2 = 71.295K/hPa, and k 3 = 375463 K 2 /hPa.

Signal processing scheme
Our algorithm shoots a signal toward the LEO satellite with three declared parameters.The first of them is the step size (d s ) referring to the length of a single segment of the signal path.The second parameter is the unit vector ( ûv ), pointing from the GNSS satellite to the LEO satellite, i.e., the equivalent of the two shooting angles (e.g., azimuth and elevation angle).After reaching the destination (as defined below), its distance to the LEO is calculated; if the distance exceeds the value of the third declared parameter-distance difference (d d )-the program performs ( ûv ) correction and shoots again in the direction of the satellite.The entire algorithm is repeated until it hits the LEO satellite according to a specified value equal to or less than d d or is terminated after a specified number of iterations (k = k limit ).The scheme of ray-tracing is presented in Fig. 2 in four main steps: (1) shooting toward the satellite with initial ( ûv ) (Fig. 2A); (2) calculating the gradients (Fig. 2B); (3) calculating next ray position (Fig. 2C); and (4) verifying the d d difference (Fig. 2D).Then, the iter- Each iteration begins with the calculation of successive nodal points of the signal path according to the following equations: where X contains the coordinates of a nodal point on a straight line (determined in Cartesian coordinates), X b is the coordinate vector of a nodal point in a deflected line (determined in Cartesian coordinates) and i is the nodal point number (i = 0-transmitter coordinates) assuming X 0 = X b 0 .The X i and X b i coordinates are additionally converted from the Cartesian to ellipsoidal latitude, colatitude, as well as geodetic altitude or ellipsoidal altitude using the SPICE toolkit (Acton Jr 1996;Acton et al. 2018) (S-spacecraft, Pplanet, I-instrument, C-orientation, E-events).Those three nodal point coordinates are named G i for a straight line and G b i for a deflected line, respectively.The conversion is done to fit ERA5 apriori refractivity field original coordinate system (WGS84 for latitude and longitude; EGM96 for altitude).ECMWF makes a simplification assuming WGS84 coordinates to be a sphere.
The algorithm basically calculates two signal paths, straight and bent.In the case of a straight line, Eq. 2 is invariably used.In the case of bent lines, the formula becomes a function of height (X i H ). The new positions of the nodal points are calculated according to Eq. 2 until X i H reaches a defined height threshold (e.g., H lim = 100 km).Below this height, the algorithm additionally calculates the refractivity value (N i ) at the point G b i using the 3D spline interpolation and then converts it to refractive index (n) according to Eq. 4: where N i is the refractivity (Eq. 1) at the point G b i .

Interpolation of refractivities
The interpolation scheme begins by verifying if the point G b i lies within the boundaries of the grid model.If it lies outside of the model boundaries, there are three variants.The first is when a point is located below the minimum height of the grid model (G H min ).The refractivity value is then assigned to the point with the coordinates (Gr , where Gr and Gr λ i are the colatitude and longitude of the nearest to the G b,β i and G b,λ i grid nodes.The second variant is for the point located above the maximum height of the grid model (G H max ).The algorithm uses then the U.S. Standard Atmosphere 1976 USSA (1976) model, from which the refractivity values are calculated at specific levels up to a height of H lim .Then the refractivity is interpolated to the point (Gr ), where G b,H i is the height of the nodal point, utilizing spline interpolation (McKinley and Levine 1998).If the longitude or colatitude of the point is outside of the grid range but altitude is in the range, the refractivity values are interpolated using the same method for the point (Gr The third variant is when the point lies within the model boundaries.Then the refractive values are interpolated to the point G b,β i , G b,λ i , and G b,H i also by the spline interpolation method.The spline method is set to minimize the total curvature of the interpolated surface.It fits a function to a specified number of nearest input points while passing through the sample points (McKinley and Levine 1998).After calculating the value of the refractive index, the refractivity gradient (∇N ir ) of ray point with X b i coordinates is calculated according to the scheme below.

Calculations of gradients
The mean model value of the refractive index is calculated at each altitude level from G b i coordinates.The algorithm finds the boundary altitude levels for each ray point and calculates the vertical gradient of refractivity according to Eqs. 7 and 8: where N H max and N H min are mean refractivities at Gr H max and Gr H min altitude levels and dN vox dH is the vertical gradient of refractivity within the voxel containing the ray point.
For differentiation at each voxel, additional linear interpolation of dN vox dH is introduced to obtain where dN ivox dH is the dN vox dH value at G b,H i altitude.The linear interpolation is made with reference to G b,H i where the minimum value is dN vox dH calculated for the previous voxel (in reference to ray path direction), and the maximum value is dN vox dH from Eq. 7.
For horizontal gradients of refractivity-the rate of change of refractivity with respect to changes in the corresponding distance in colatitude/longitude direction-we calculate the refractivity gradient for colatitude β and longitude λ following Eqs.9 and 10 consecutively.
where N β min are obtained similarly to their counterparts from altitude direction according to Eq. 5.In the longitude direction, we have where N λ max and N λ min are refractivities at Gr Transformed to Cartesian coordinates (∇ N XY Z ir ) (Gates 2004), it reads: where the G b,βl i stands for latitude of G b,β i value.In Eq. 12, simplification is performed using spherical equations.Since the presented algorithm calculates gradients based on the relative differences of coordinates in a voxel, the distinctions between gradient calculations with ellipsoidal and spherical coordinates are negligible.

Calculations of the unit vector
The next step is to calculate the second derivative of the wave equation (h i ) at new position i and the new ray position according to Javaherian et al. (2020): With the knowledge of h i , we obtain the ray position update g i and its normalized version ĝi as follows: The new ray point position is: After receiving the new ray position, the new ray tangent vector for the next (in ray path direction, i + 1) ray position is calculated following Eq.17.
For X b and X, the process is repeated until X H i > H lim with Eq. 2 and 3.The destination is declared reached and the signal modeling process is terminated when X H i is greater than the receiver height (X H f ).Between H lim and X H i , Eq. 2 is used assuming a straight line.At this point, the length of the vector real distance r d is calculated as follows: where X f is the position of the receiver.If r d > d d , a correction to the initial ûv is made based on the correction factor ( ûc ): where , j is the point with the lowest r d value.The corrected unit vector reads: where k is the number of signal iteration.

Next iterations and excess phase estimation
After the unit vector ( ûv ) is calculated, the operator ends the first iteration.The next iteration starts if the condition r d > d d is fulfilled.In such a case, the entire process is repeated.Otherwise (i.e., if r d < d d ), the iteration process is terminated.
The resulting total delay (in meters) along the reconstructed signal paths is obtained as follows: where a = 1 and c are, respectively, the first and the last ray points with altitudes below G H max .

Calculation of bending angle
In order to comprehensively test our solution, we decided to calculate the bending angle using the RO processing package (ROPP) ROPP ( 2022).The first step of retrieving the bending angle was to smooth the simulated excess phase values using the Savitzky-Golay (S-G) function (Schafer 2011).The S-G function is a smoothing method based on a polynomial least squares approximation to noisy data.In our study, we used the third-order filter with a frame length of 101.The settings used should have only a small effect on the absolute values of the simulated phase excesses but should be sufficient to smooth the differences between successive excess phase values.Smoothed values were then used in the ROPP tool with geometric optics (GO) and wave optics settings (WO) to retrieve bending angles.Assuming the local symmetry, the ROPP determines bending angle α at a given impact parameter p from the geometry of satellite position and velocity data.In a simpler approach-GO-only one ray is observed at one time.The bending angle is calculated here from ray directions to LEO and GNSS satellite ( ûL and ûG , respectively) according to Eq. 22.
Then, based on the Doppler shift, excess phase, and center of the local curvature of the Earth at the occultation point the values of ûL and ûG are iteratively solved.
In the case of WO processing, the ROPP algorithm uses the canonical transform algorithm described by Gorbunov and Lauritsen (2004) and Fourier integral operator (FIO) to transform the wave field.The wave field is used here as a coordinate system in which each ray is a single-valued function of the new coordinate.The resulting FIO is a composition of a reference signal, Fourier transform, and an amplitude function.The last part of processing is an ionospheric correction made by a linear combination of bending angles at a common impact parameter.More processing details about the ROPP tool with the theoretical background can be found in ROPP (2022).

Experiment design
To validate the ray-tracing results, cumulative excess phase, bending angles, and refractivity profiles are derived from the LC (ionosphere-corrected) measurements in the atmPhs files.
Each of the 10 selected RO cases comprises 2000 to 4000 transmitter-receiver (GNSS satellite-RO satellite) location pairs.We ray-traced every 10th measurement in the atmPhs file, which leads to 200 to 400 rays per RO event.For each event, the processing time for ray-tracing was recorded.The major ray-tracing parameters are (see Fig. 3, for graphical explanation): a. step size (d s ): the incremental distance that the ray-tracer adds to the previous signal location (see Eqs. 2 and 3), b. threshold to satellite (d d ): the maximum distance accepted as convergence criterion for the ray-tracer (see Eq. 18), c. grid source as we alter the refractivity data source for ray-tracing (see Eq. 4), changing from a one-dimensional ERA5-based profile (1DERA5; vertical profile of refractivity calculated as the mean value over all 3DERA5 refractivities per altitude layer), three-dimensional ERA5based data (3DERA5; the refractivities N i and gradients of refractivity ∇ N XY Z ir were calculated using interpolated refractivities from a 3D ERA5 grid), and finally, the CDAAC retrieved profiles from the atmPrf files were used as an input for the ray-tracer (COSMIC-1), d. horizontal resolution as we verify whether the horizontal gradients of refractivity smoothing will have an impact on the excess phase, e. vertical resolution as we verify whether using different numbers of vertical layers will have an impact on the calculations of vertical gradients of refractivity, ranging from basic pressure levels in ERA5 up to all layers for which CDAAC refractivity products (atmPrf) are provided.
Overall we got 18 different simulation settings to find the optimal solution in terms of time and accuracy.Those are two different horizontal resolutions (1 • and 5 • ), three different values for the step size (1 km, 5 km, and 50 km), and three different grid sources (1DERA5, 3DERA5, COSMIC).In addition, four different numbers of vertical layers (33 and 642 for 1DERA5 and 3DERA5 as well as 40, 142, and 642 for COSMIC) and three different values for the thresholds to satellite (0.01 km, 0.1 km, and 0.3 km) were analyzed.Table 3 in the appendix summarizes all applied settings.One of them (number 7 with 3DERA5-based grid, 5 km d s , 0.01 km d d , 642 vertical layers, and 5 • horizontal resolution) was designed to cover the most precise ray-tracer settings with very high computational costs to check its capabilities.

Agreement between observed and simulated excess phase
The basic validation of our RT implementation, yet in some cases challenging, is direct comparison of excess phase simulated minus observed excess phase.Figure 4 presents relative differences in the excess phase from two selected simulation settings in comparison with the atmPhs excess phase on the LC band for each RO case.Labels A and B in Fig. 4 correspond to the 33 vertical layers 3DERA5 grid and 1DERA5, respectively.The heat maps present sample settings fast in terms of computation time, however with lower quality.Labels 'C' and 'D' are 40 vertical layers COSMIC-1 grid (where the COSMIC-1 grid model is a spherically symmetric field) with d s = 5 and d d = 0.3 km and very dense 3DERA5 grid (5th and 7th setting numbers, respectively, from Table 3 in Appendix A).The 'C' solution is based on observed values of refractivity, while the 'D' is the one with very time-consuming highest quality processing.
Differences are based on rays where the ray point's height was within a defined height range.The full results of the comparison for each setting option as a mean value for all cases can be found in Appendix A.
The results show the best agreement between the simulated and the CDAAC (atmPhs file) excess phase is obtained Note that the largest difference up to an altitude of 30 km in all cases was 35%.At the corresponding altitude of 25 to 30 km, this can be described as a 1-1.5 m difference in excess phase.At the altitudes of most interest, below 10 km, the average difference is 4% (4-6 m excess phase) for the COSMIC-1-based solutions and 8% (8-12 m excess phases) for the 3DERA5 profiles.However, if we take a look only at 'D' profiles with a very low threshold and a large number of vertical layers, significant improvement can be noticed with an average bias around 3%.
Comparison It is also significant that above 10 km altitude, the simulated excess phase values are higher than those provided by the CDAAC files, while below 10 km they are too low.Moreover, the results show no dependency between retrieval quality and profile complexity represented by the categories in Fig. 1.

Agreement between observed and simulated bending angle
Next to the absolute values of the simulated excess phase, we checked the quality of our results by recalculating them to the bending angle.Figure 5 presents differences between bending angles obtained with the ROPP tool from our simulations and observed profiles.For each case, we used geometric optics (GO) and wave optics (WO) processing schemes with  3 in Appendix A).The number in brackets indicates the RO category (see Table 1).Negative values mean that the simulated excess phase is lower than the one observed by CDAAC, while a positive that it is higher.Black values mean that the tangent point of the simulated signal was above 5 km altitude  3, respectively).
In the case of WO processing, the bias between observed and simulated bending angle for 'D' (higher resolution) profiles is lower than 2%, while for 'A' (lower resolution) profiles it is below 3% in 8 out of 10 cases.Cases 5 and 7 reached differences up to 3.5%.
GO processing produced significantly higher biases up to 3% for the 'D' (higher resolution) solution and up to 6% for the 'A' solution (lower resolution).Two very high spikes were noted in Case 4 with differences around 8%.

Parameters of refraction and vertical gradients of refractivity dN ivox /dH
Regarding the differences between successive program settings, we verified the differences between interpolated refractive values at given levels.These differences were calculated between the refractivity values from the CDAAC reference profiles and the 3DERA5-based as well as COSMIC-1based solutions.Figure 6 presents the percentage difference between the refractivity values obtained from 3DERA5 and COSMIC-1-based solutions with those provided by the CDAAC files.
As expected, the simulations based on the COSMIC-1-based solutions as input for the ray-tracer are almost  3 in Appendix A and option 'A' from Fig. 4) Fig. 6 Percentage difference in refractivity between simulations based on 1) 3DERA5 (red line) as well as 2) COSMIC-1 (black line) and atmPhs excess phase on LC for selected height ranges for all ten simulated cases.The numbers in brackets are the numbers of the category (see Table 1).The figures were based on one ray path from one arbitrarily selected moment (transmitter-receiver position) when all simulation cases reached similar values of excess phase (in terms of ±2%) consistent with data taken straight from reference atmPhs files.The mean differences are around ±4% for the profiles up to 40 km altitude and around ±1% below 10 km altitude.For simulations derived from a 3DERA5-based solution, a pronounced relative difference of up to 12% (absolute difference of up to 8 in refractivity) is seen.The mean difference of the 3DERA5-based solution is around 6%.The difference for both 3DERA5 and COSMIC-1-based solutions comes from differences in simulated and observed ray path geometry.In terms of 3DERA5, they are also influenced by disagreement in 3DERA5-based refractivity grid values and the reference atmPhs files.
In Fig. 7a similar analysis is provided for vertical gradients of refractivity dN ivox dH derived from the 3DERA5 and the COSMIC-1 data, respectively.Vertical gradient of refractivity profiles can be a more accurate tool to verify the variability of refractivity than its absolute values.Table 1) As with refraction, the COSMIC-1-based refractivity vertical gradient of refractivity values overlaps much more with the reference ones compared to the 3DERA5-based solution.In this case, however, the differences are much larger and amount to about 14% in case 6 and up to 10% in the rest of the cases for the COSMIC-1-based solution.For 3DERA5 they reach up to 20% in individual cases.The reason for the differences is the same as in the case of absolute refractivity values.Those come from (1) differences in simulated and observed ray paths geometry, thus different initial unit vectors, and (2) differences in refractivity using 3DERA5 and COSMIC-1-based grids as input.It is difficult to state explicitly the impact of these errors on the final excess phase values; thus, in Sects.4.4-4.9, the impact of different settings on the final excess phase results is presented.

Impact of step size
It is clearly visible that changing the step size from 1 to 5 km has almost no impact on the results (Fig. 8).However, the change in step size from 5 to 50 km seems to be more significant and probably linked to atmospheric structures, which cannot be resolved anymore by the ray-tracer.This is particularly visible within the height range from 20 to 40 km, where the relative differences are up to 10%.Comparing different grid sources, it is also worth noting that the change of d s from 5 to 50 km for the 3DERA5-based solution gives higher values of excess phase.The same change for the COSMIC-1-based 142 layers solutions yields slightly smaller relative differences in excess phase.In the case of 40 layers, the relative differences are larger and negative proving the importance of the grid source.

Impact of convergence criteria (d d )
The d d parameter has no clear effect on the simulated excess phase in the atmosphere at heights above 40 km as the relative differences (Fig. 9) are lower than ±0.5%.On the contrary, in the 20 and 40 km altitude range, the relative differences around 1 − 5% between the solutions with d d = 0.01 and d d = 0.3 km are visible and up to 1.5 % when d d is changed from 0.1 to 0.3 km.This effect is independent of the input grid, i.e., it appears in both the 3DERA5 and COSMIC-1based solutions.Below 5 km the effect increases to almost 5%.We can find two possible explanations for this effect.In the case of altitudes between 20 and 40 km it is connected with the iteration cycle.The 20 to 40 km altitude range is where the ray tracer needs to do the second/third iteration of signal processing with small but significant initial unit vector changes.Ray paths simulated with an d d lower than 0.3 km need additional iteration.Thus, it changes the ray path geometry and final excess phase.In the case of the lower troposphere, the accuracy of targeting the receiver seems to have a significant impact.A more detailed discussion of this effect is provided in the next section.

Impact of grid source
As expected, the grid source seems to have a very significant influence on the results.It is also clearly visible in Fig. 10.
The highest relative differences in excess phase can be seen above 70 km altitude, where COSMIC-1 and 3DERA5- based solutions yield larger excess phase values than the 1DERA5 solution by up to 20%.The reason for that can come from the fact of extrapolating the data to this height but it is worth noting that the relative differences in the upper parts of the analyzed atmosphere drop as the length of d s increases.In addition, it is interesting to note the difference between solutions at altitudes between 20 and 40 km.Once more, there are significant relative differences in excess phase values.At this level, a 10% difference is about 0.5 m of excess phase.Also in this height range, there is a slight decrease in relative error when step size d s increases to 50 km.It seems that the grid sources have a significant impact on the accuracy of the simulated excess phases, too, especially when using lower values for step size d s .

Impact of horizontal resolution
The results in Fig. 11 show that the final excess phase bias does not depend on the horizontal resolution of the grid model.However, the results are slightly different when the grid is denser (around −0.3% at heights below 50 km altitude).The largest difference is visible above 80 km and is up to −5.9% which can be a result of extrapolating the data to this level.The analyzed solutions do not differ significantly for different step sizes (d s equal 5 and 50 km).

Impact of vertical resolution
According to Fig. 12, we can see a significant impact of the vertical profile density on the COSMIC-1-based solution.The denser the profile, the more accurate the ray-traced excess phases, regardless of step size d s .
It is worth noting that extending the number of layers from 142 to 642 has almost no impact (lower than 1%) on the results.However, increasing the number of layers from 40 to 142 seems to have significant impact (about 10%), mainly in a height range between 20 and 40 km altitude.It provides positive feedback concerning the computation time discussed in Sect.4.9.

Impact of the ray-tracing settings on the processing time
The results summarized in Table 2 show that processing time differs significantly between various settings.The most timeconsuming simulation was the one which we assumed to have the best accuracy-solution number 7 from Table 3 (see Appendix A)-with on average 3 times longer computation time in comparison to second place on this list (see Fig. 15).The most significant difference in terms of processing time between used settings is seen for 1 • horizontal resolution,

Discussion
In the following section we inspect our results in light of previous research, in three major blocks: Sect.5.1 comparison to observations and RO processing results, Sect.5.2 impact of the ray-tracing algorithm fine-tuning, and Sect.5.3 where we identify limitations of our ray-tracing algorithm.

Excess phase data validation
There are many other studies considering the RO simulations using the geometric Høeg et al. (1995); Poli and Joiner (2004); Wee et al. (2010) or wave optics method (Hordyniec et al. 2019;Gorbunov and Kirchengast 2018;Beyerle et al. 2002;Benzon et al. 2015), but the direct comparison of sim-ulated to observed excess phase is rarely presented.We only found the study by Wee et al. (2010) who demonstrated the difference for single CHAMP profile on the level of 0.5 m between ray-traced (2D RT GO) and observed excess phase (5-10 km), and these are similar to the 1-4% reported in this study Fig. 4. Another paper, (Shao et al. 2009), to some extend similar to ours, with large number of observations (1158 RO events) but, simulate the signal in the 2D vertical plane only in the straight lines around the local tangent points, reports that there is up to 5% difference in relative difference between observed and ray-traced signal.These two results show similar performance to our solution (see, e.g., Figure 4).Figure 4 shows that the simulated excess phase values are similar to the ones provided in the atmPhs files (with relative differences ranging from 0 to 20% in excess phase).The results seem to be not dependent on the complexity of the profiles indicated in the paper by Healy (2012), as similar level of discrepancy is observed for profile Case 01 (leftmost graph in 4) and for profile Case 10 (rightmost graph in 4).
The bottom 5 km in our solution is not always resolved; in few cases (Fig. 4 black boxes) the ray-tracer stops as it cannot converge.We attribute this limitation to inhomogeneous water vapor distribution in the lower troposphere, lack of ionosphere model, and the multipath of signal propagation (Zou et al. 1999;Gorbunov et al. 1996;Gorbunov 2002;Jensen et al. 2004).
Related issue is that, in the most of the settings 'A', 'B', 'C' (Fig. 4 with large threshold of 0.3 km and low number of vertical layers 33 or 40), simulated excess phase is negatively  biased.This effect is less visible in 'D' setting Fig. 4 and setting number 7, in Table 3, with smaller threshold 0.01 km and 642 vertical layers.The average value of percentage difference is around 3% below 5 km altitude.The follow-up procedure that estimates the bending angle from excess phase rate was applied similarly to other studies, e.g., Healy (2001).We actually used only two simulation settings, namely 'A' (large threshold of 0.3km and 33 vertical layers) and 'D' (smaller threshold 0.01 km and 642 vertical layers).In the case of 'A' settings, the solution with a higher value of threshold (0.3 km) gives less accurate results with biases up to 3.5 %; however, it should be still acceptable for, e.g., GNSS tomography processing as the absolute values of excess phase do not differ too much above 5 km altitude (see Fig. 4).We demonstrated that for settings 'B' the bias of up to 2% between simulated and observed one with WO processing (Fig. 5).The value of 2% is comparable to 4% as reported by Healy (2001).
The quality of the ray-traced signal paths depends primarily on the quality of the refractivity field (discrepancy between the numerical field and observations) (Möller and Landskron 2019).Therefore, the source of refractive values at a given altitude will be the most critical factor.In the case of refraction (Fig. 6), it should be noted that there is a high agreement between the simulated refraction profiles from the COSMIC-1-based solutions with their original (atmPhs files) observations.The relative differences that arise from the simulations on the ERA5 (3DERA5 and 1DERA5) data indicate that a ± 10 % difference in refraction values at levels above 40 km altitude does not significantly affect the final excess phase value.This difference is expressed as a percentage in Fig. 6 and, in most cases, does not exceed 5-10%.The most significant relative differences appear above the 20 km height level.In the paper by Shao et al. (2009), the differences between GPS RO data and simulated data reach a standard deviation of 1-3 refraction units at heights below 14 km altitude.That can be considered as a 1-2% difference which is similar to our differences between COSMIC-1-based solutions and observations from atmPhs files at heights below 14 km.
The values of the vertical gradients of refractivity in Fig. 7 are in turn directly translated into the refraction values shown in Fig. 6.This can be seen especially in the 4th case profile in which the simulated COSMIC-1-based data above 20 km altitude deviate significantly from their originals.The reason for this may be the interpolation method used and the toosparse distribution of profile layers in some of the analyzed solution (33 or 40 layers in 'A', 'B' and 'C').
We would like to underline that the ray-traced satelliteto-satellite excess phase values obtained in this work will be further used as input for space-and ground-based GNSS tomography studies.The additional observations strengthen the observation geometry and therefore contribute to a more reliable tomographic solution (Möller and Landskron 2019).In the ground-based tomographic solutions, the slant total delay (STD) uncertainty is up to 8-9 cm at very low elevation angles (< 5 • ) (Kačmařík et al. 2017) with similar results for 2D ray-tracing study (Hordyniec et al. 2018).As reported by Brenot et al. (2019), the mean relative error of slant wet delay (SWD) is ± 9.1%.Therefore, we expect that the attained satellite-to-satellite excess phase quality, better than 10% in most cases (Fig. 4), is sufficient for tomography.
Moreover, as presented in Bossema et al. (2018), the tomographic reconstruction depends on the distribution of the projection angles.The RO excess path of ≈ 100 m in the lowermost section of the troposphere (below 5 km) is going to traverse the tomography model that spans at least 300 km, with at least 5 to 10 voxels (horizontal resolution).A large number of empty voxels can significantly increase the RMSE of the tomography solution (Yao et al. 2020).Filling them with additional RO simulated data can work similarly to, e.g., adding virtual observation (Zhang et al. 2021) and decreasing tomography uncertainties.Hence, the RO excess phase can provide the new projection angle, almost perpendicular to the ground-based signal delay observations, improving the tomographic retrievals of water vapor content.
Considering integrated GNSS tomography as the final purpose of the algorithm presented in this manuscript, it is worth highlighting the differences in the use of 3D RT instead of 2D RT on ground-based observation.The expected outcome is the reduction of differences between observed and simulated slant delays by 0.1 to 4 cm (Nafisi et al. 2011).

Influence of ray-tracing settings
According to different settings of the ray-tracer (Figs. 8,9,10,11,12), we found that changing the grid source (1DERA5, 3DERA5, COSMIC-1) as well as the number of vertical layers and threshold d d might have the largest influence on the results, whereas the influence of horizontal resolution and step size d d might be considered insignificant.The relative differences between 1 • and 5 • horizontal resolution (Fig. 11) are visible only in the highest simulated altitude (above 80 km) and might be connected to different nodes, which rays pass through after reaching a height of 90 km (the place where we start to calculate refractivity).Another viable explanation is that at these altitudes refractivities are extrapolated which limits the quality of retrieval.Additionally, the 1 • solution needs fewer iterations to hit the satellite, assuming a threshold to the satellite of d d = 0.3 km but still takes almost 10 times longer to receive the results.Thus, we consider the 5 • as sufficient.
In the case of d d , the largest relative differences between 0.1 and 0.3 km (around 3 to 4%) appear between 20 and 40 km (Fig. 9) and by around 5% when we change from 0.01 to 0.1km.The increase of d d changes the number of itera-tions and, in consequence, the ray path geometry as seen in the unit vector ûv .The atmospheric region between 20 and 40 km altitude is where the refractivity starts to increase and the bent ray path starts to differ significantly from the straight ray path.While d d changes from 0.01 to 0.3 km, the signal stops iterating earlier, after reaching a 0.3 km distance to the satellite (cases of 0.1 or 0.01 km need more iterations), but the place of crossing the 90 km altitude border remains close to each other.Furthermore, we have to note that the change of d d does not mean the change in tangent point location by the value of d d due to ray geometry.However, it significantly reduces the bias between simulated and observed excess phases.This can be explained by the fact that with 50Hz of sampling and an average speed of LEO satellite around 25.000 km/h, the change of the receiver position may reach up to 200 m during one 'shot'.Then, with high values of d d the algorithm may target the incorrect position of the receiver.Therefore, the change of d d makes the ray-tracer simulate a more consistent path of the signal.The relative differences in interpolated refraction at 20 to 40 km altitude due to differences in grid source (3DERA5 vs. 1DERA5 vs. COSMIC-1 solutions) or step size (1 km vs. 5 km vs. 50 km) can increase this impact.It is visible in Figs. 6,8,10 and 9.The discrepancy in presented excess phase differences between Fig. 9 and Appendix A (0 vs. 5%) comes from the fact that in Fig. 9, we only compare the results of simulated solutions, while in Appendix A, we compare our solutions with CDAAC observations.
In our study, based on the experiments, we consider the d d value of 0.3 km as the best option for simulations in terms of quality and computation time.It is reasonable to use a larger d d value to obtain the simulation results faster (e.g., in GNSS tomography) as the absolute bias between simulated excess phase is very low above 5 km altitude.However, as a stand-alone tool, the best option is to use lower d d with the value of 0.01 km.
Another variable analyzed was step size d s (Fig. 8).It may be interesting to note that a step size of 50 km provides an excess phase that is slightly closer (by 1 − −2%) to the reference solution than using a step size of 5 km below 5 km altitude, according to Appendix A. We might assume that the denser grid or denser ray point field gives higher values of the final excess phase due to deeper penetration of the signal into the atmosphere.We should also consider it as a more precise solution.The possible explanation for that is the same as in the case for d d -the more rays cross the border of 5 km altitude, the higher mean difference we obtain, as in general, the ray-tracer simulates slightly too low values of excess phase below 10 km altitude.In the case of high relative differences above 80 km, the explanation can be the same as for horizontal grid resolutions, i.e., the ray passes through the atmosphere in different places or it is a result of extrapolation of refractivity values.
The difference between the 1 and 5 km step size seems to be insignificant, and it takes much longer (5 times) to simulate the signal with d s = 1 km.Therefore we consider d s = 5 km as the best trade-off between good quality and calculation time.Because of the dynamically changing vertical gradients of refractivity, in this study, we did not look into height adaptive d s values as was stated in the paper of Liu et al. (2018).As for relative differences between 20 and 40 km, it is worth considering the grid source and the number of vertical layers (Fig. 12).It is not a surprise that the best-fitting solutions are COSMIC-1-based solutions with 142 and 642 layers, especially for altitudes between 20 and 40 km.Forty layers in a COSMIC-1-based solution is not enough; however, increasing the number above 142 seems to have no influence on the results, while it is highly computation time-consuming.It is worth noting that the more layers we have, the more significant problems with convergence are obtained (see Sect. 5.3).In the case of the grid source, by analyzing the results more carefully (Fig. 10), we can see that COSMIC-1-based rays passing through the troposphere need more iterations than the ones from 3DERA5 due to problems with convergence.Meanwhile, 1DERA5 and 3DERA5-based simulations need the same number of iterations to hit the receiver.Hence, the relative differences must come directly from the refraction values in the profiles.That is expected and in agreement with the findings in Möller and Landskron (2019).In general, 1DERA5-based solutions have higher values of excess phase than the ones from 3DERA5 and COSMIC-1.The relative differences between COSMIC and ERA5 may also come from different refractivity coefficients.The one used in our study (Eq. 1) differs slightly from the one used in CDAAC processing with k 1 = 77.60K/hPa (vs.77.689 K/hPa in our study) and k 3 = 373000 K 2 /hPa (vs.375463 K 2 /hPa in our study).Values of k 2 are the same (k 2 = 71.295).We do not consider those differences to make a significant impact on the 3DERA5-based results; however, in further studies, they should be set to same numbers.Regardless of the data source, the computation time remains the same (Table 2).Depending on d s , we can see that significant differences, similar to the other analyzed parameters, appear at levels of 20 to 40 km altitude.
The observed discrepancies are not related to the horizontal resolution, which showed no difference between 1 and 5 • (Fig. 11).This is the only case analyzed where the same d s , grid, and d d were applied.The reason, therefore, lies between these three variables.More specifically, it comes from the too low density of grid layers at heights of 20-40 km for the 33 (1DERA5 and 3DERA5-based) and 40 layers (COSMIC-1based) analyzed in most cases.Hence, the lack of differences between 142 and 642 grid layers is described below.Denser simulations (with d s = 1 km) and more precise solutions (with d d = 0.1 km) can cancel out this effect, but it becomes more pronounced for less precise solutions.Averaging the refrac-tion values in the 1DERA5-based simulations, as shown in Fig. 10, also brings up a significant difference.

Limitations
As mentioned in the Introduction Sect. 1, geometrical optics ray-tracing simulations used in this study are not necessarily expected to produce results close to observation due to multipath or ducting conditions in lower parts of the troposphere (Jensen et al. 2004).Geometrical optics is needed here for GNSS tomography purposes and therefore cannot be replaced with wave optics techniques.We are aware that possible problems connected with large uncertainties of excess phases and low levels of ray paths convergence may affect tomography solution when relatively large d d values of 0.1 and 0.3 km are used.However, these problems are connected strictly with chosen technique and will be tested in our further studies.Furthermore, when higher accuracy is needed, the solution with a very low d d of 0.01 km can be used.A good solution to balance the computation time and quality of the results may be also to mix the settings between high d d solution (e.g., 5 • , 5 km d s , 33-142 layers, 0.3 km d s ) for altitudes above 5 km and low d d (e.g 5 • , 5 km d s , 642 layers, 0.01 km d s ) below 5 km altitude.
In terms of low convergence, as visible in Fig. 4, all of the 3DERA5-based simulation cases went below 5 km altitude (most of them stopped converging at 3 km altitude.)However, for 5 out of 10 cases, the COSMIC-1-based simulations did not converge.Figure 13 shows that the ray passing through the COSMIC-1-based refraction grid starts to repeat its positions during the consecutive iterations at some height level and, thus, cannot reach the assumed threshold to the satellite.Looking more closely into Eq.21, we can see that ûv is divided by a factor of 2. We consider '2' to be the most effective.However, changing the factor makes a significant influence on the number of iterations and ray convergence.The solution for that issue might be a dynamic factor instead of 2 based on the number of iterations.On the other hand, it is possible to increase the value of d d ; however, it might be costly for the accuracy of final excess phase value simulations.
The difference between 3DERA5 and COSMIC-1-based simulations can as well come from the lack of horizontal gradients of refractivity.Ray-tracer was designed to work in 3D space while in COSMIC-1-based solutions, the grid is, in fact, one-dimensional copied to each horizontal node.Due to zero-valued horizontal gradients of refractivity, the algorithm stops iterating earlier as it moves only between sparsely distributed vertical layers and jumps between two ray paths closest to accurate results.
The second issue might be related to the interpolation of the refractions between grid nodes.The spline technique used in this paper appears to be effective at low altitudes below The numbers in legends are iteration numbers.The first iteration is omitted due to the initial unit vector leading the ray path through the Earth's surface 20 km; however, it deviates significantly at altitudes around 20 to 40 km.Changing the way of interpolation or taking denser vertical profiles into account may affect the results, however, it will be tested in our further studies.

Conclusions
Taking into account all collected data from 18 types of simulation settings and 10 RO cases, it should be concluded that the developed ray-tracer is a stable and fast way to simulate RO.The simulation accuracy and calculation time depend on the settings used.During the study, settings that have the most significant impact on the final excess phase value were determined.
Differences between non-grid source-based settings reach up to 3% of excess phase (1 m) value at altitudes below 10 km and 11% (up to 20 cm) above this height.For obvious reasons, in addition to the vertical profile, the most significant impact on the results is related to the grid data source as the difference between COSMIC-1 and 3DERA5 reaches up to 1.5 m excess phase.However, in most cases, it is below 0.5 m.Another important factor is the threshold to the satellite with relative differences reaching up to 5 % below 5 km altitude.Those changes significantly affect calculation time (up to 500%).The calculation time for the fastest settings is a bit over 1 s for one transmitter-receiver pair with a tangent point below 5 km altitude, which turns into 120 s to 180 s per one RO event with 50 Hz (or 12 s to 18 s for 5Hz), thus, up to 2000 (or 200) different excess phase values.For the initial settings with the longest processing time, the simulations take over 9 s.Therefore, we can conclude that in 3D geometrical optics, ray-tracing with the use of high-resolution refractivity models is applicable only when the highest precision of simulations is needed.Then, it is worth decreasing the threshold to satellite as much as possible.
The lack of simulations below 5 km altitude in some of the cases of the COSMIC -1-based solution is due to diffi-culties with signal convergence.The situation is worth further investigation by modifying Eq.21 so that it is more dependent on the signal iteration number.In other simulations below 5 km altitude, significant differences between simulations and observations are due to the influence of multipath effects.However, this cannot be readily resolved because of the use of the GO RT technique.GO RT is more dependent on irregular water vapor structures in the lower troposphere than WO.
Overall the difference between COSMIC-1-based solutions and excess phase LC observations reaches 5% of excess phase values at altitudes below 5 km ( 4.0 m), 2% at 5 km to 10 km ( 0.8 m), and around 11% ( 0.2 m) above 10 km.Results are comparable to other studies of this type in terms of mean relative error of simulations around 10%.The bending angle retrieved from ray-traced solution agrees with observations within 2% relative error for high-resolution ray-tracing and up to 3.5% for low-resolution ray-tracing.
Suggested settings combining speed, stability, and solution accuracy are step size d s = 5 km, the threshold to satellite d d = 0.3 km, approximately 140 layers in the vertical profile, and a 5 • horizontal resolution of the grid.Changing the number of vertical layers to 642 and decreasing d d to 0.01 km will increase the quality of the results, however, significantly increase the calculation time (by 461%).
The differences arising in the interpolation of vertical gradients of refractivity at altitudes above 20 km (up to 10% for COSMIC-1-based solutions) affect the obtained results and require further study, but they do not have a significant impact on rays reaching the lower troposphere.

A Appendix
Table of excess phase relative differences between simulations and atmPhs files (Figs. 14,15 and Table 3) Table of simulation time of one iteration of signal processing at selected height threshold on each setting.

Fig. 2
Fig. 2 Scheme of ray-tracer calculations.Panel A presents the initial ( ûv ) on a straight and bent line.Panel B presents the gradient of the refractivity calculation scheme.Panel C shows the transformation of gradients to the Cartesian coordinate system and the calculation of the

Fig. 3
Fig. 3 Ray-tracer validation parameters.The top row shows the signal discretization parameters: Panel a step size, Panel b. the threshold to the satellite.The bottom row shows the parameters related to the grid source: Panel c grid source, Panel d the horizontal resolution, Panel e between 'A' (d s = 5 km, d d = 0.3 km, 5 • horizontal resolution on 33 vertical layers 3DERA5 grid) and 'B' (d s = 5 km, d d = 0.3 km, 5 • horizontal resolution on 33 vertical layers 1DERA5 grid) solutions, which covers all the same settings except the difference in 1D vs 3DERA5 profile, shows slightly lower values of bias (around 3% lower on average) when the 'A' solution was used.In all of the analyzed cases, 3DERA5-based simulations were closer to the observations than 1DERA5-based.

Fig. 4
Fig. 4 Percentage difference between simulated and atmPhs excess phase on LC band in selected height ranges for all simulated cases.Label 'A' corresponds to a solution with d s = 5 km, d d = 0.3 km, 5 • horizontal resolution on 33 vertical layers 3DERA5 grid (5th setting in Table 3 in Appendix A).Label 'B' corresponds to a solution with d s = 5 km, d d = 0.3 km, 5 • horizontal resolution on 33 vertical layers 1DERA5 grid (9th setting in Table 3 in Appendix A).Label 'C' corresponds to a solution with d s = 5 km, d d = 0.3 km, 5 • horizontal

Fig. 5
Fig. 5 Differences between simulated and observed bending angle.Dashed lines were processed with the wave optics option of ROPP, and solid lines with geometric optics options of ROPP.Red lines are raytracing options with very high density of grid model and low threshold 'D'd s = 5 km, d d = 0.01 km, 5 • horizontal resolution on 642 ver-

Fig. 7
Fig. 7 Percentage difference between simulated and atmPhs refractivity-based vertical gradients of refractivity values for selected height ranges for all ten simulated cases.Red lines are the 3DERA5based, while black lines are the COSMIC-1-based simulations.The

Fig. 8 Fig. 9
Fig. 8 Relative differences in excess phase from simulations with varying step size (1 km, 5 km, 50 km) and grid sources (3DERA5 and COSMIC-1) as input.Negative values suggest that the first solution mentioned in the field 'step size' has lower values

Fig. 10
Fig. 10 Percentage difference in excess phase from simulations with varying GRID sources (1DERA5, 3DERA5, and COSMIC-1) and step sizes (d s = 1 km, 5 km, and 50 km) for all ten simulated cases.Negative values suggest that the first solution mentioned in the field 'GRID source' has lower values

Fig. 12
Fig. 12 Percentage difference in excess phase values from simulations with varying vertical resolutions (40, 142, and 642 layers) and step size (5 km, 50 km), for all ten simulated cases.Each column shows the differ-

Fig. 13
Fig. 13 Convergence of the same ray path (same epoch and coordinates of GPS and LEO).Panel a: 3DERA5 grid (33 vertical layers), and panel b: COSMIC-1 (40 vertical layers).The numbers in legends are iteration numbers.The first iteration is omitted due to the initial unit vector leading the ray path through the Earth's surface

Fig. 14
Fig.14Percentage difference between simulations and CDAAC file excess phase for every solution on selected height and band.X labels are sorted in the following way: horizontal grid resolutiond s -GRID

Table 2
Percentage difference in mean computation time for each setting.Values in the 'Difference' row mean how much longer (in percent) it takes to calculate the excess phase with the settings specified in the row values Parameter

Table 3
Ray-tracer settings used in the study Flowchart of the ray-tracing algorithm used in this study.'Eq.number' refers to the equations used in this study.Green dots represent the positions of the ray points where interpolation and gradient calculations are performed, while red dots denote ray points located outside the voxel model.H lim denotes the altitude of the receiver Fig.16