Atmospheric bending effects in GNSS tomography

Abstract. In GNSS tomography, precise information about the tropospheric water vapor distribution is derived from integral measurements like ground-based GNSS slant wet delays (SWDs). Therefore, the functional relation between observations and unknowns, i.e. the signal paths through the atmosphere have to be accurately known for each station-satellite pair involved. For GNSS signals observed above 15 degrees elevation angle, the signal path is well approximated by a straight line. However, since electromagnetic waves are prone to atmospheric bending effects, this assumption is not sufficient anymore for lower elevation 5 angles. Thus, in the following, a mixed 2D piecewise linear ray-tracing approach is introduced and possible error sources in reconstruction of the bended signal paths are analyzed in more detail. Especially, if low elevation observations are considered, unmodeled bending effects can introduce a systematic error of up to 10−20ppm, on average of 1−2ppm into the tomography solution. Thereby, not only the ray-tracing method but also the quality of the a priori field can have a significant impact on the reconstructed signal paths, if not reduced by iterative processing. In order to keep the processing time within acceptable limits, 10 a bending model is applied for the upper part of the neutral atmosphere. It helps to reduce the number of processing steps by up to 85% without significant degradation in accuracy. Therewith, the developed mixed ray-tracing approach allows not only for a correct treatment of low elevation observations but is also fast and applicable for near real-time applications.

a bending model is applied for the upper part of the neutral atmosphere. It helps to reduce the number of processing steps by up to 85% without significant degradation in accuracy. Therewith, the developed mixed ray-tracing approach allows not only for a correct treatment of low elevation observations but is also fast and applicable for near real-time applications.

Introduction
For conversion of precise integral measurements into two-or three-dimensional structures, a technique called tomography has 15 been invented. In the field of GNSS meteorology, the principle of tomography became applicable with the increasing number of Global Navigation Satellite System (GNSS) satellites and the build-up of densified ground-based GNSS networks in the 1990s (Raymond et al., 1994;Flores Jimenez, 1999). Since then, a variety of tomography approaches based on raw GNSS phase measurements (Nilsson, 2005), double difference residuals (Kruse, 2001), slant delays (Flores Jimenez, 1999;Hirahara, 2000) or slant integrated water vapor (Champollion et al., 2005) has been developed for the accurate reconstruction of the 20 water vapor distribution in the lower atmosphere. An overview about the major developments within this field of research since Flores Jimenez (1999) is provided by Manning (2013).
While in most tomography approaches, observations gathered at low elevation angles are discarded (Bender et al., 2011;Champollion et al., 2005;Hirahara, 2000), straight line signal path reconstruction is sufficient for the determination of the path lengths. However, Bender and Raabe (2007) showed that especially low elevation observations can be a very useful source of strengthen the observation geometry and therewith, contribute to a more reliable tomography solution. However, a correct treatment of low elevation observations requires more advanced ray-tracing algorithms. The first paper which deals with bended ray path reconstruction in GNSS tomography was published by Zus et al. (2015), with main focus on the reconstruction of the signal paths for delay estimation but also for assimilation of GNSS slant delays into numerical weather prediction systems.
Most recently, Aghajany and Amerian (2017) published their results about 3D ray-tracing in water vapor tomography and 5 briefly analyzed its impact on the tomography solution.
Based on the existing studies, in the following, a more detailed discussion of possible error sources in signal path reconstruction is provided. Therefore, Sect. 2 describes the effect of atmospheric bending and its handling in GNSS signal processing.
Section 3 describes the principles of GNSS tomography and how the basic equation of tomography is solved for wet refractivity. Section 4 introduces the concept for reconstruction of signal paths using ray-tracing techniques. Hereby, the modified 10 piecewise linear ray-tracing approach is described -including its ability for reconstruction of the GNSS signal geometry. In Sect. 5, the defined ray-tracing approach is applied to real SWDs and its impact on the tomography solution is assessed and validated against radiosonde data. Section 6 concludes the major findings.

Atmospheric bending effects in GNSS signal processing
The effect of atmospheric bending on GNSS signals is related to the propagation properties of electromagnetic waves. In 15 vacuum, GNSS signals travel with the velocity of light. When entering into the atmosphere, the electromagnetic wave velocity changes, dependent on the electric permittivity ( ) and magnetic permeability (µ) of the atmospheric constituents and the frequency of the electromagnetic wave. The ratio between the velocity of light c in vacuum and the velocity ν in a medium defines the refractive index n.
For signals in the microwave frequency-band, n ranges from 0.9996 to 1.0004. Thus, n is usually replaced by refractivity N , expressed in mm/km (ppm).
The GNSS signal delay in the lower atmosphere, also known as slant total delay (ST D), is related to refractivity by the following equation (Bevis et al., 1992): The first term of Eq. (3) describes the change in travel time due to velocity changes along the true ray path R. The second term (about three orders of magnitude smaller than the first term) is related to the difference in geometrical path length between the true (R) and the chord signal path (S). According to Dalton's law, the refractivity of air can be split up into a hydrostatic and a wet component: N = N h + N w . Therewith, the GNSS signal delay reads: The slant wet delay (SW D) depends on the wet refractivity along the true ray path R.
The slant hydrostatic delay (SHD) results from the the hydrostatic refractivity along R and, by definition, from the additional 5 path length due to atmospheric bending.
While signal path S follows from the straight line geometry between satellite and receiver, the true signal path R depends in addition on the hydrostatic and the wet refractivity distribution along the signal path (see Sect. 3 for more details).
In GNSS signal processing, the integral along the signal path is usually replaced by the zenith delay and a mapping function.
10 Therefore, Eq. (4) is rewritten as follows: where ZHD is the zenith hydrostatic delay, ZW D is the zenith wet delay and mf h and mf w are the corresponding mapping functions, which describe the elevation (ε) dependency of the signal delay. The elevation and azimuth (α) dependent first-order horizontally asymmetric term G(ε, α) reflects local variations in the atmospheric conditions: see MacMillan (1995), Chen and 15 Herring (1997) or Landskron and Böhm (2018). In practice, e.g. when using VMF1 mapping function (Böhm et al., 2006) or similar mapping concepts, the tropospheric delay due to atmospheric bending is absorbed by the hydrostatic mapping function term mf h . Comparisons between ray-traced SHD(ε) and 'mapped' SHD(ε) = ZHD·mf h (ε) slant hydrostatic delays reveal that about 97 % of the atmospheric bending effect is compensated by the VMF1 hydrostatic mapping function (see Appx. A for further details).

The principles of GNSS tomography
According to Iyer and Hirahara (1993), the general principle of tomography is described as follows: where f s is the projection function, g s is the object property function and ds is a small element of the ray path R along which the integration takes place. In GNSS tomography, g s is usually replaced by wet refractivity N w , and integral measure f s by SW D (the prefactor of 10 6 vanishes if ds is provided in kilometer and SW D in millimeter).
A full non-linear solution of Eq. (9) for wet refractivity is not of practical relevance since according to Fermat´s principle, first order changes of the ray path lead to second order changes in travel time. In consequence, by ignoring the path dependency in the inversion of N w along ds and by assuming the ray path as a straight line, a linear tomography approach can be defined which is well applicable to SW Ds above 15 degrees elevation angle (Möller, 2017). However, with decreasing elevation angle, the true signal path deviates significantly from a straight line. In consequence, by ignoring atmospheric bending, a systematic 5 error is introduced in the tomography solution. In order to overcome this limitation, in the following an iterative tomography approach is defined in which the bended signal path is approximated by small line segments. Similar to the linear tomography approach, thereby the neutral atmosphere or parts of it are discretized in volume elements (voxels) in which the refractivity N w,k in each voxel k is assumed as constant. Consequently, Eq. (9) can be replaced by: where d k is the travelled distance in each voxel. Assuming l observations and m voxels, a linear equation system can be set up.
In matrix notation it reads: where SW D is the observation vector of size (l, 1), N w is the vector of unknowns of size (m, 1) and A is a matrix of size (l, m) which contains the partial derivatives of the slant wet delays with respect to the unknowns, i.e. the travelled distances d k 15 in each voxel.
Solving Eq. (11) for N w requires the inversion of matrix A.
The inverse A −1 exists if A is squared and if the determinant of A is non-zero, otherwise matrix A is called singular.

20
Unfortunately, singularity appears in GNSS tomography in most cases since the observation data is 'incomplete' and matrix A is not of full rank. Therewith, Eq. (13) becomes ill-posed, i.e. not uniquely solvable. In order to find a solution which preserves most properties of an inverse, in the following matrix A is replaced by pseudo inverse A + . According to Hansen (2000) the pseudo inverse is defined as follows: where U and V are orthogonal, normalized left and right singular vectors of A and matrix S is a diagonal matrix, which contains the singular values in descending order. In case a priori information (N w0 ) can be made available, it enters the tomography solution as first guess as follows: where matrix U , V and S are obtained by singular value decomposition of matrix A T · P · A + P c . The weighting matrices P and P c are defined as the inverse of the variance-covariance matrix C for the observations and C c for the first guess, respectively. Assuming, that the observations are uncorrelated, the non-diagonal elements of C and C c are zero and the diagonal elements are defined as follows: whereby σ ZW D = 2.5mm reflects the uncertainty of the ZW D. The values for σ T , σ q and σ p were taken from heightdependent error curves for pressure (p), temperature (T ) and specific humidity (q) as provided by Steiner et al. (2006) for the ECMWF (European Centre for Medium-Range Weather Forecasts) analysis data. For further details, the reader is referred 10 to Möller (2017).

Reconstruction of GNSS signal paths
Assuming that the geometrical optics approximation is valid and that the atmospheric conditions change only inappreciably within one wavelength, the signal path is well reconstructible by means of ray-tracing shooting techniques (Hofmeister, 2016;Nievinski, 2009). Thereby, the basic equation for ray-tracing, the so-called Eikonal equation, has to be solved for obtaining 15 optical path length L.
From Eq. (18), a number of 3D and 2D ray-tracing approaches have been derived for the reconstruction of ground-based and space-based GNSS measurements and of their signal paths through the atmosphere (Hobiger et al., 2008;Zou et al., 1999).
The main difference between both observation types is related to the observation geometry. While for space-based GNSS 20 observations derived in limb sounding, the bending angle is usually described as function of impact parameter a, for groundbased observations elevation and azimuth angle are used for characterizing the signal geometry. In consequence, the optimal ray-tracing approach will be significantly different for various observation geometries.
In order to find an optimal approach for operational analysis of ground-based measurements, Hofmeister (2016) carried out a number of exploratory comparisons. Based on the outcome, the 2D piecewise linear ray-tracer was defined as the optimal 25 reconstruction tool for the iterative reconstruction of the atmospheric signal delays including atmospheric bending. It is limited to positive elevation angles but it is fast and almost as accurate as the 3D ray-tracer. However, for the use in GNSS tomography, the ray-tracing approach had to be further modified. In the following, the developed ray-tracing approach but also its impact on the GNSS tomography solution are discussed in more detail.

Piecewise linear ray-tracer
The starting point for the 2D piecewise linear ray-tracer is the receiver position in ellipsoidal coordinates (ϕ 1 , λ 1 , h 1 ), the 'vacuum' elevation angle ε k (see Fig. 1) and the azimuth angle α under which the satellite is observed. In case of GNSS tomography, these parameters can be determined with sufficient accuracy from satellite ephemerides and the receiver position -assuming straight line geometry. Therewith, the initial parameters for ray-tracing (see Fig. 1), i.e. the geocentric coordinates 5 (y 1 , z 1 ) and the corresponding geocentric angles (η 1 , θ 1 ) read: where R G is the Gaussian radius, an adequate approximation of the Earth radius 15 with a and b as the semi-axes of the reference ellipsoid (e.g. GRS80). The z-axis connects the geocenter with the starting point, the y-axis is defined perpendicular to the z-axis in direction (azimuth angle) of the GNSS satellite in view. After setting the initial parameters, the 'true' ray path is reconstructed iteratively by making use of ray-tracing shooting techniques. Therefore, total refractivity derived from an a priori field is read in and pre-processed for ray-tracing. Hereby, the input data is interpolated vertically and horizontally to the vertical plane, spanned by the y-and z-axis.
In the ray-tracing loop, for each height layer h i+1 with i = 1 : (t − 1) whereby t defines the top layer of the voxel model, the geocentric coordinates and the corresponding angles are computed as follows: where d i is the reconstructed path length between height layer h i and h i+1 (h i+1 > h i ). It depends on the observation geometry but also on the atmospheric conditions (refractive indices n i and n i+1 ). By default, for our analysis, the spacing between two height layers h i and h i+1 was set to 5 m, which corresponds to a maximum path length d i of 100 m -assuming an elevation angle of 3 • (5m/sin3 • ).
The ray-tracing loop is repeated until ε t − ε k + g bend is smaller than a predefined threshold (e.g. 10 −6 degrees). While the elevation angle ε t is obtained by Eq. (29) for i = t − 1, the correction term g bend accounts for the additional bending above 30 the voxel model. Since atmosphere is almost in state of hydrostatic equilibrium, g bend can be well approximated by a bending model, like the one of Hobiger et al. (2008):  In order to reduce the impact of possible refractivity errors on the reconstructed ray paths and in further consequence on the tomography solution, ray-tracing was carried out iteratively. Therefore, the refractivity field obtained from the first tomography 15 solution replaces the initial refractivity field for ray-tracing for the next iteration and so on. The processing is repeated until N w converges.

The empirical ray-bending model
Besides the refractivity field, the quality of the reconstructed ray paths might be also affected by errors in the bending model as defined by Eq. (32). Comparisons of the bending model with ray-traced bending angles on a global 10 • x10 • grid over the 5 period of one year reveal that the error in bending is usually kept below 0.8 arcsec. Assuming a GNSS site near sea-level and an elevation angle of 5 • , an error in bending angle of ±0.8 arcsec causes an error in path length of up to ±10m, i.e. the reconstructed GNSS signal enters the voxel model slightly earlier or later than the observed GNSS signal. In Fig. 4, the bending error is visualized as pointing error at voxel model top. However, for the tomography solution this effect is too small to be significant. Thus, it can be concluded that the bending model of Hobiger et al. (2008) is well applicable for reconstruction

Ionospheric bending effects
Beyond, also the ionosphere influences GNSS signal propagation. In order to assess the impact of free electrons in the ionosphere (above 80 km altitude) on the signal path, the electron density model by Anderson et al. (1987) Figure 5 shows the obtained vertical profiles of ionospheric refractivity, exemplary for frequency f 1 . The higher the signal-frequency f the lower the phase velocity through the ionosphere and the less is its refraction. Following the approach by Wijaya (2010), the ray paths in the ionosphere were 10 reconstructed, separately for GPS L1 and L2. The analysis revealed significant path differences between the 'true' ray path and its chord line but also between the two signal-frequencies. Assuming a VTEC of 455 TECU and an elevation angle of 3 • , the maximum deviation from the straight line signal path is 800 m for L1 and 550 m for L2 respectively, at h = 400km, slightly below the layer of peak electron density. Fortunately, ray path deviation decreases significantly with decreasing VTEC and altitude to a few tens of meters at h = 13.6km (the upper rim of the troposphere at which the top of the voxel model was

Impact of atmospheric bending on the tomography solution
In the following, the differences between straight line and bended ray-tracing are further analyzed. For highest consistency, the ray-tracing approach defined in Sect. 4 was used for both, straight line and bended ray-tracing. The only difference is that in case of straight line ray-tracing the ratio n i /n i+1 in Eq. (27) was set to '1'. Thereby, it can be guaranteed that only the impact of atmospheric refraction is assessed.

Expected drying effect
In the beginning, the ray position is equal for both methods but diverges with increasing height. Thereby, the bended ray is travelling in most cases 'above' the straight ray, i.e. the straight ray enters the voxel model top 'earlier' than the bended ray.
This leads to the effect that the straight ray remains longer in the voxel model than the bended ray, i.e. the straight ray path within the voxel model (h < 13.6km) is systematically longer than the bended ray path. The differences between both ray 10 paths are plotted in Fig. 6 (top) as function of elevation angle. Therefore, ALARO model data was selected as input for the bended ray-tracer.
The additional ray path decreases rapidly with increasing elevation angle. Thus, a mixed ray-tracing approach can be defined, which considers ray bending only for ε ≤ 15 • . Beyond, the additional ray path is below 0.1 km, and straight line ray-tracing is sufficient for ray path reconstruction.
Figure 6 (top) shows also that in some cases even for low elevation angles the difference in path length is small (below 0.1 km). This appears when the ray enters the voxel model not through the top layer but through a later surface of the voxel model.

5
In this particular case, the difference in path length between both ray-tracing approaches is negligible (be aware that only the entire distance through all voxels is comparable for both ray-tracing approaches but not the individual distances in each voxel).
Figure 6 (bottom) shows the expected drying effects in the tomography solution caused by errors in the reconstructed signal paths assuming straight line geometry. Hereby, it is distinguished between the drying effect caused by the additional ray path (dN w 1 ) and the drying effect caused by the fact that the straight line travels through lower layers of the voxel model (dN w 2 ).

10
Both effects were assessed as follows: whereby SW D b and SW D s are the slant wet delays obtained by ray-tracing through ALARO model data along the bended 15 and the straight line ray path, respectively. The variables d k,b and d k,s are the corresponding path lengths within the voxel model. The sum of their differences along the ray paths are identical with the additional ray paths plotted in Figure 6 (top).
Both drying effects have to be considered as additive and are strongly connected to the current atmospheric conditions as well as to the parametrization applied for interpolation of the refractivity field. In our analysis we assumed an exponential decrease of refractivity between the vertical layers of the voxel model and applied a bi-linear interpolation method for horizontal 20 interpolation between the grid points.

Results from the Austrian GNSS tomography test case
In order to study the impact of bended ray-tracing on the tomography solution, a GNSS tomography test case was defined. The corresponding settings are summarized in Table 1.
Figure 7 (left) shows the differences in wet refractivity between Sol1 and Sol2 (as defined in Table 1). Even though on 25 average over all voxels no bias in wet refractivity is observed, specific voxels show differences in wet refractivity of up to 10ppm, particularly if due to bending different voxels than in the straight line solution are traversed.
Figure 7 (right) shows the differences in wet refractivity between the first two iterations of the mixed ray-tracing approach (Sol2). In this particular case, refractivity differences are smaller than 0.05ppm, which implies that the a priori model used for ray-tracing is already close to the 'true' atmospheric conditions, i.e. in this particular case no further iteration was necessary.

30
From all differences in wet refractivity over 248 epochs in May 2013, a maximum of 14.2ppm, a bias of 0.12ppm and a standard deviation of 0.24ppm were obtained. Although the bias and standard deviation over all voxels is small, differences of

Validation with radiosonde data
For validation of the mixed ray-tracing approach against straight line ray-tracing, the tomography derived wet refractivity fields were compared with radiosonde data at the airport of Innsbruck (ϕ i = 47.3 • , λ = 11.4 • , h = 579m). First, the radiosonde data 5 obtained once a day between 2 and 3 UTC were pre-processed, i.e. outliers in temperature were removed and dew point temperature was converted to water vapor pressure and further to wet refractivity. Finally, the radiosonde profiles were vertically interpolated to the height layers of the voxel model and the tomography derived wet refractivity fields were horizontally interpolated to the ground-position of the radiosonde launching site, respectively. Figure 8 shows the differences in wet refractivity as function of height above surface, exemplarily for two epochs in May, 2013. In both cases, the bended ray-tracing approach helps to reduce the tomography error by about 1−2ppm, especially in the lower 4 km of the atmosphere. Largest differences are visible when the bended ray traverses other voxels than its chord line. This appeared in about 2 % of the test cases, especially 5 if observations below 10 degrees elevation angle enter the tomography solution.

Conclusions
GNSS signals which enter the neutral atmosphere at low elevation angles (ε < 15 degrees) are significantly affected by atmospheric bending. In case the bending is neglected when setting up design matrix A, a systematic error of up to 10 − 20ppm, on average of 1 − 2ppm is introduced into the GNSS tomography solution. This error can be widely reduced if atmospheric 10 bending is considered in reconstruction of the signal paths. Therefore, a 2D piecewise linear ray-tracing approach was defined, which describes the bended GNSS signal path by small line segments. By limiting the length of the line segments to 100m in case of ε = 3 • or even shorter for higher elevation angles, the 'true' signal path can be widely reconstructed. However, the quality of the reconstructed signal paths depends primarily on the quality of the a priori refractivity field. Comparisons between refractivity fields derived from standard atmosphere and ALARO weather model data reveal that a refractivity error of 30ppm can cause a ray deviation of up to several hundred meters, i.e. the distance traveled in each voxel but also the number of traversed voxels is prone to misallocations. In consequence, reliable a priori data, e.g. derived from numerical weather model data, are recommended for GNSS tomography.
Nevertheless, if reliable a priori data are not available or if the quality is unknown, iterative ray-tracing helps for reducing 5 the impact of wet refractivity errors on the tomography solution. Therefore, the wet refractivity field obtained from an initial tomography solution is used for reconstruction of the signal paths for the next iteration. The processing is repeated until the tomography solution converges. This ensues usually after two iterations. Beyond, a bending model, like the one provided by Hobiger et al. (2008) helps to significantly reduce computational cost by describing the remaining bending in the higher atmosphere (above the voxel model). In consequence, the ray-tracer can be stopped right after the reconstructed signal leaves 10 the voxel model. In case of h t = 13.6km, the number of processing steps is reduced by 85%, which is a tremendous reduction in processing time without significant loss of accuracy.
In contrast, ionospheric bending effects have less impact on the GNSS tomography solution. Even during periods of solar maximum, ray path deviation caused by ionospheric bending is negligible for signals in L-band (1 − 2GHz). However, even if ionospheric bending has no impact on the tomography solution, first and higher order ionospheric effects should be taken into 15 account when processing GNSS phase observations.
Besides, comparisons with radiosonde data revealed that if atmospheric bending effects are considered in GNSS tomography, the quality of the tomography solution can be improved by 1 − 2ppm. Within the defined test case, especially voxels in the lower 4km of the atmosphere benefitted from the applied mixed ray-tracing approach. Due to significant optimization, the mixed ray-tracing approach ensures processing of large tomography test cases in adequate time. A test case with 72 GNSS 20 sites and 7 x 9 x 15 voxels can be processed in less than two minutes. Thus, the developed mixed ray-tracing approach is applicable also in near real-time and therefore well suited for operational purposes.
Code availability. The 2D piecewise linear ray-tracer for GNSS tomography as well as the RADIATE ray-tracer are part of the Vienna VLBI and Satellite Software (VieVS). The code of the RADIATE ray-tracer is available at https://github.com/TUW-VieVS/RADIATE. For more details to VieVS, the reader is referred to http://vievswiki.geo.tuwien.ac.at.

25
Appendix A: Unmodelled bending effects in the Vienna hydrostatic mapping function In case of VMF1 (Böhm et al., 2006) or similar mapping concepts, azimuthal asymmetry is not considered and for convenience, only a single hydrostatic mapping coefficient per site (a h ) is determined as follows: (A1) Figure A1. The unmodeled geometric bending effect in VMF1 hydrostatic mapping function (dg bend ), exemplary for VLBI sites Fortaleza, Brazil and Wettzell, Germany. Analyzed period: Jan-Feb 2014 where b h is 0.0029, c h depends on the day of year and latitude and mf h (ε) is defined as the ratio between SHD(3 • ) and ZHD, obtained by ray-tracing through numerical weather model data. For assessing the remaining unmodeled geometric bending dg bend (ε, α), ray-traced slant hydrostatic delays were compared with 'mapped' slant hydrostatic delays as follows: where ZHD is the zenith hydrostatic delay obtained by vertical integration, g bend (ε, α) is the geometric bending effect as 5 obtained by ray-tracing, mf h (ε) is the VMF1 hydrostatic mapping function determined by SHD(3 • )/ZHD and mf h0 (ε) is the hydrostatic mapping function determined by SHD 0 (3 • )/ZHD, whereby SHD(3 • ) and SHD 0 (3 • ) are the slant hydrostatic delays obtained by ray-tracing for an vacuum elevation angle ε k = 3 • with and without geometric bending, respectively. Figure A1 shows the remaining unmodeled geometric bending as obtained for six elevation angles (and 16 equidistant azimuth angles), exemplary for the two VLBI sites Fortaleza, Brazil (ϕ = −3.9 • , λ = 321.6 • , h = 23m) and Wettzell, Germany 10 (ϕ = 49.1 • , λ = 12.9 • , h = 669m). In case of ε = 3 • , almost no bending error is visible since mf h (ε) was tuned for this elevation angle. However, for other elevation angles, the unmodeled geometric bending is about 3 % of the slant hydrostatic delay, e.g. up to ±5mm at 5 • elevation angle. In case of Wettzell, dg bend (ε, α) is mostly negative, i.e. the 'mapped' SHD is smaller than the observed SHD and vice versa for Fortaleza. So far, these small variations are neglected when using VMF1 hydrostatic mapping function in GNSS signal processing.