An initial investigation of the non-isotropic feature of GNSS tropospheric delay

Tropospheric delay is a significant error source in Global Navigation Satellite Systems (GNSS) positioning. Slant Path Delay (SPD) is commonly derived by multiplying Zenith Tropospheric Delay (ZTD) with a mapping function. However, mapping functions, assuming atmospheric isotropy, restrict the accuracy of derived SPDs. To improve the accuracy, a horizontal gradient correction is introduced to account for azimuth-dependent SPD variations, treating the atmosphere as anisotropic. This study uncovers that, amidst atmospheric dynamics and spatiotemporal changes in moisture content, the SPD deviates from that based on traditional isotropy or anisotropy assumption. It innovatively introduces the concept that SPD exhibits non-isotropy with respect to azimuth angles. Hypothesis validation involves assessing SPD accuracy using three mapping functions at five International GNSS Service (IGS) stations, referencing the SPD with the ray-tracing method. It subsequently evaluates the SPD accuracy with horizontal gradient correction based on Vienna Mapping Function 3 (VMF3) estimation. Lastly, the non-isotropic of SPD is analyzed through the ray-tracing method. The results indicate the smallest residual (1.1–82.7 mm) between the SPDs with VMF3 and those with the ray-tracing. However, introducing horizontal gradient correction yields no significant improvement of SPD accuracy. Considering potential decimeter-level differences in SPD due to non-isotropic tropospheric delay across azimuth angles, a precise grasp and summary of these variations is pivotal for accurate tropospheric delay modeling. This finding provides vital support for future high-precision tropospheric delay modeling.


Introduction
Tropospheric delay affects the accuracy of Global Navigation Satellite Systems (GNSS) positioning.It is caused by the neutral atmosphere, specifically the dry air and water vapor in the troposphere can refract GNSS signals, causing them to take a longer path through the atmosphere than they would in a vacuum.This delay can lead to the errors in the calculated position, especially for the signals that have longer passage in the troposphere (Yang et al., 2021).According to the elevation angles of a satellite, the tropospheric delay can be classified into Zenith Tropospheric Delay (ZTD) and Slant Path Delay (SPD).For high-precision GNSS navigation, as SPD can reach 20 m or more at elevation angles below 10° while ZTD is only about 2-3 m (Fan et al., 2019a, b), it is essential to precisely estimate tropospheric delay.
The most efficient method for estimating tropospheric delay in GNSS data processing uses meteorological data through the Ray-Tracing (RT) SPD.However, the raytracing is complex and the meteorological data involves delays, making it less commonly employed.To produce SPD, the most prevalent method is to multiply ZTD by a Mapping Function (MF) SPD.Models like the Saastamoinen, Hopfield, and Black models can directly derive ZTD, however they have lower precision and are frequently employed for standard single-point positioning (Saastamoinen, 1972;Black & Harold, 1978;Hopfield, 1969;Zhang et al., 2012).In precise single-point positioning, tropospheric delay is often estimated as a parameter, and the validation of this estimation is dependent on the accuracy of the observations.The development of the ZTD estimation method is able to meet the needs of positioning.Since many nations and regions have constructed global and regional Numerical Weather Models (NWM) models, the use of external meteorological data, such as the reanalysis data from NWM, to retrieve ZTD has become a research hotspot in recent years (Jones et al., 2020;Leuenberger et al., 2020;Zhou et al., 2020).
The Neill Mapping Function (NMF) (Niell, 2004), Vienna Mapping Function (VMF) series (Boehm & Schuh, 2004), and Global Mapping Function (GMF) (Böehm et al., 2006) are currently the most widely adopted as accurate mapping functions, all employing a third-order continuous fraction structure.These models utilize either spherical harmonic functions or the latest European Centre for Medium-Range Weather Forecasts (ECMWF) fitting values to represent the coefficients.Boehm & Schuh (2004) proposed the VMF based on the NWM, which was improved to VMF1.To enhance the precision of VMF1 in estimating SPD at low elevation angles, Landskron et al. (2018) proposed the more accurate VMF3.Furthermore, GPT is an empirical model of meteorological parameters based on VMF, providing a global grid with two horizontal resolutions of 5° × 5° and 1° × 1°.The main components of GPT consist of the empirical coefficients of the hydrostatic and non-hydrostatic mapping functions, as well as other meteorological quantities, making it a complete tropospheric empirical model.Studies have demonstrated that the accuracy of the VMF series models is generally superior to that of the NMF and GMF models (Guo et al., 2015).However, SPD varies with azimuth angles, and the existing mapping functions are based on the isotropy of the troposphere, assuming that the atmosphere above the station is isotropic in the horizontal direction.This results in the mapping function only considering the influence of elevation angle, neglecting the horizontal variation of the atmosphere, which restricts the accuracy of estimating the SPD (Ifadis et al., 1986).
A horizontal gradient term is frequently included in high-precision positioning algorithms to minimize the errors brought on by the mapping function which neglects the changes in SPD at various azimuth angles.The horizontal gradient component expands SPD in the azimuth domain as a first-order Fourier series, neglecting higher-order terms (Fan et al., 2019a, b), and parameter estimation is then used to solve for the anisotropic components in the east-west and north-south directions.These anisotropic components are multiplied by the corresponding horizontal gradient mapping function to obtain the "anisotropic" correction term.Currently, the main horizontal gradient mapping function used is a first-order continuous fraction structure in terms of elevation angle (Chen et al., 1997).However, the horizontal gradient assumes that the atmosphere is anisotropic, while the atmosphere flows and water vapor also have high spatiotemporal variability.Therefore, SPD is not simply isotropic or anisotropic.
Ray-tracing involves dividing the atmosphere into multiple units, such as thin layers or grids, and assuming a constant refractive index within each unit.Recursive calculations are then used to connect adjacent units, enabling the delay along the entire signal propagation path to be determined.Ray-tracing can be further categorized into one-dimensional, two-dimensional, and three-dimensional methods (Aghajany et al., 2017).The one-dimensional method assumes that the atmosphere is spherically symmetric and calculates only one azimuth angle.The two-dimensional method calculates at different azimuth angles, assuming that the trajectory of the ray at a certain azimuth angle always remains in the same plane.The three-dimensional method accounts for the leap-out-of-plane effect that may be caused by horizontal gradient (Hofmeister, 2016).Meteorological data plays an important role in ray-tracing, with the fifth generation of European Centre for Medium-Range Weather Forecasts Reanalysis (ERA5).It includes meteorological parameters such as temperature, pressure, and humidity, with a time delay of about 5 d.Studies have demonstrated that the ray-tracing method, based on the ERA5 reanalysis product, can obtain high-precision tropospheric delays (Zhou et al., 2020).However, despite its high accuracy, ray-tracing is not widely applicable for precise positioning due to computational complexity, time delay, and dependence on meteorological parameters.
In summary, many scholars have delved deep into the methods of estimating tropospheric delay and modeling.However, further exploration is required to investigate the relationship between SPD and azimuth angle (Fan et al., 2019a, b).Firstly, the accuracy of SPD estimated by VMFs and GMF mapping functions with the incorporation of the horizontal gradient needs to be evaluated (Du et al., 2020).Secondly, research on the horizontal variation of SPD is scarce, and neither the mapping function nor the horizontal gradient can accurately simulate the real path of signals passing through the troposphere (Qiu et al., 2020).By neglecting the variation of tropospheric delay at different azimuth angles, the improvement of estimation accuracy is limited.This is especially problematic in high-precision positioning and long baseline solutions, where the impact of tropospheric delay on the solution results is more significant (Davis et al., 1985;Hou et al., 2021).
SPD exhibits neither isotropic nor anisotropic behavior.This article presents a hypothesis that the SPD in the horizontal direction exhibits non-isotropic behavior.Importantly, this perspective highlights that despite consistent elevation angles, the SPDs at certain azimuth angles may demonstrate isotropic properties, while the others may exhibit anisotropic tendency.To verify this hypothesis, the study assesses the accuracy of estimated SPD with two approaches: MF-SPD and MF-SPD with horizontal gradients.These evaluations are conducted by using RT-SPD as a reference point for comparison.
One key aspect for improving the precision of tropospheric delay modeling lies in accurately describing the characteristics of SPD.The introduction of this concept provides potential pathways for enhancing the accuracy of SPD estimation.
This paper is arranged as follows.In the introduction the definition and the issues in the estimation of SPD are described, followed by an overview of the experiments which are designed to verify the conjecture that SPD is non-isotropic.In the methodology section, the basic principles of using mapping functions, horizontal gradient, and ray-tracing to estimate SPD are introduced, and the data and processing strategies are described in the next section.In the results and discussion, ZTD and horizontal gradient provided by International GNSS Service (IGS) and ERA5 reanalysis meteorological data are utilized.The SPD estimated by VMF1, VMF3, and GMF is evaluated using RT-SPD as a reference.The MF-SPD with the highest precision is selected and corrected with the horizontal gradient, and the difference between the corrected SPD and RT-SPD is evaluated to explore the existence of non-isotropy.Finally, the conclusions of the experiments are presented.

Mapping function and horizontal gradient
The SPD at an elevation angle ε can typically be expressed as the product of ZTD and the mapping function.In practical operations, the wet and dry components of SPD are calculated separately: where L represents SPD, L z w is the zenith wet delay, and L z h is the zenith hydrostatic delay, f i (ε) is mapping function, where subscript i represents the hydrostatic (h) or wet (w) component.The mapping functions used in this study, namely VMF1, VMF3, and GMF, can all be expressed by a third order continuous fraction with (1) different coefficient values.The general formula can be written as: where a i , b i , and c i are the coefficients of the mapping function, and v i is the correction term.Equation ( 2) can be expressed concisely as: For VMF1, the coefficients b h , b w , and c w can be con- sidered as constants with b h = 0.0029, b w = 0.0014, and c w = 0.04391 .The coefficients a h and a w can be obtained through the GPT2w model or interpolated from the 6-hourly gridded data provided on the VMF website (https:// vmf.geo.tuwien.ac.at/).The dry component c h can be obtained using the following equation (Böehm et al., 2006): where d is DOY represents the day of year, ϕ is latitude, c 0 , c 10 , c 11 , and are the constants that are related to the location of the station.
VMF3 adjusts the coefficients b and c of VMF1 using the least squares to better approximate RT-SPD at low elevation angles.The fitting formula for b h is (Landskron et al., 2018): where A 0 denotes the mean value, A 1 and B 1 the yearly amplitudes of the coefficient, and A 2 and B 2 its semi- annual amplitudes.
Compared with other mapping functions, the results of GMF are less affected by the bias caused by elevation while solving the problem of low time resolution of VMF1 coefficients.GMF shares the same coefficients b and c as VMF1, while the value of coefficient a and mean value a 0 are determined by Eq. ( 6) (Boehm et al., 2006): (2) The global grid of mean value a 0 and amplitude A are expanded into spherical harmonics coefficients based on the least-squares method.Using Eqs. ( 6) and ( 7), a can be determined for any location and date.In addition, the other coefficients in VMF3 and GMF can be computed using the empirical model GPT3.In practical applications, altitude conversion is necessary to normalize the hydrostatic mapping function from sea level to the target altitude in GPT3 (Landskron et al. 2018).
Equations ( 1) and ( 2) show that MF-SPD is independent of azimuth and isotropic for the same elevation angle.However, the Earth's atmosphere is not a spherical symmetry model, so estimating SPD with a spherical symmetry model will inevitably lead to errors.In high-precision positioning, horizontal gradient correction is usually introduced to correct the error caused by neglecting the difference in SPD at different azimuth angles (Macmillan et al., 2013).The software GAMIT is widely used in GNSS data processing (Davis et al., 1985), which includes a module for horizontal gradient correction.The horizontal gradient correction model in GAMIT is defined as: where C is a constant with the value of 0.003, G N is the atmospheric horizontal gradient parameter in the northsouth direction, G E is the atmospheric horizontal gradi- ent parameter in the east-west direction, σ is the azimuth angle.As can be seen from the Eq. ( 9) that the horizontal gradient assumes that SPD is anisotropic.

Ray-tracing
Ray-tracing is currently the most accurate approach for estimating the SPD.In this paper, RT-SPD was used as a reference to investigate the changes in SPD at different azimuth angles.The method is based on the ray propagation model, which assumes that signals propagate along a ray path in the atmospheric layers.Thus, by tracking the signal's path and speed of propagation in the atmosphere, the signal delay in the atmosphere can be calculated (Guo et al., 2015).
This article adopts Two-Dimensional (2D) ray-tracing.The computational approach for 2D ray-tracing bears a striking resemblance to that of One-Dimensional (1D) ray-tracing.In the 1D ray-tracing the atmospheric environment is modelled under the assumption of spherical ( 7) symmetry.Conversely, 2D raytracing is executed at different azimuth angles.However, it is presupposed that the ray path remains confined to a single plane.The fundamental divergence from the 1D methodology lies in the selection of disparate piercing points, while the underlying principle remains unaltered.Figure 1 serves as an illustrative representation of the raytracing approach, wherein the atmospheric is demarcated into multiple layers: In Fig. 1, the Earth is depicted as an ostensibly regular sphere, centered at O, while the atmosphere is partitioned into m discrete layers.P 1 designates a terrestrial station with the vertical axis at P 1 serving as the z-axis, and the horizontal direction as the y-axis.The atmospheric refractive index ( n ) for both dry air and water vapor at each path point is derived through bilinear interpolation and exponential interpolation.Refractivity ( N ) serves as an alterna- tive representation of the refractive index ( n ) and has been adopted to streamline the computational process.The relation between N and n is: The N of each layer is given by: where p w is the partial pressure of water vapor, T is the temperature, and ρ denotes the density of dry air.Z w is the compressibility factor for water vapor, and k 1 , k ′ 2 , k 3 are the refractivity coefficients that are related to the water vapor.Here, k ′ 2 is defined as (Hofmeister et al., (10)  (Rüeger, 2002): where the molar masses of dry and wet air are M h and M w , respectively, and the universal gas constant is R.
The calculation of tropospheric delay at any radius can be expressed as follows (Zhang et al., 2016): where τ (r) signifies the ranging delay, S(r) represents the distance traveled by rays, and V (r) stands for the linear distance traversed by the ray.
The expressions for these variables are elucidated as follows: where the Earth's radius is denoted as r 0 , while r rep- resents the separation between the two layers.η j charac- terizes the geocentric angles relative to the ray points P j , while θ j designates the angle of incidence.The ray propa- gation path between each stratum is denoted as S j , and the straight-line distance bridging two piercing points is encapsulated by V j .Calculation of the angle of incidence at any given r necessitates the application of Snell's law of refraction (Steen, 2000). (13 • r

Data and processing strategy
To assess the accuracy of MF-SPD and investigate the non-isotropic nature of tropospheric delay, the following data sets were selected.

Tropospheric delay products provided by IGS
The tropospheric products provided by the IGS were utilized in this study, including ZTD and horizontal gradient.The distribution of 75 IGS tracking stations and their ZTD values at DOY001 UTC18 in 2018 are shown in Fig. 2, where station names are positioned either to the left or right of the points.It can be observed that the ZTD values of IGS stations in high-latitude regions are smaller than those in mid and low-latitude regions.This can be ascribed to the greater amount of water vapor in low-latitude regions, leading to a subsequent increase in tropospheric wet delay.

Station information
In order to select suitable stations for sample analysis, this study opted for the 320 IGS stations worldwide with available tropospheric data at DOY001 UTC18 in 2018.The derived SPDs using the raytracing and three different mapping functions at each station were statistically examined, as depicted in Fig. 3.The x-axis represents the number of stations, while the y-axis indicates the average standard deviation of RT-SPD and the three MF-SPDs at each station.The selected elevation angle was fixed to 35°, with azimuth angles ranging (18) sin(θ(r)) = n(r) n(r − �r) • sin(e(r)) from 0° to 360° in 10° intervals.From the graph, the standard deviations at different stations fluctuate between 0 and 2.6 mm.Most stations exhibit the standard deviations concentrated within 1 mm, albeit not reaching zero.This illustrates that the neglection of horizontal variations in tropospheric delay at any station results in a bias.
To account for the potential impact of geographical factors on the experimental outcomes, five IGS stations-ABMF, BJFS, LAMA, NYAL, and RIO2-with a uniform global distribution were chosen for subsequent experiments.The latitude, longitude, and elevation of these stations are detailed in Table 1.

Atmospheric reanalysis data
The meteorological data used for estimating RT-SPD was obtained from the ERA5 reanalysis dataset provided by ECMWF.This dataset offers horizontal resolutions of up to 0.25° for atmospheric parameters at various pressure levels, with a temporal resolution of 1 h.In this study, the meteorological parameters at each station were interpolated for SPD calculations.The grid data used included the atmospheric pressure, temperature, relative humidity, specific humidity, and geopotential height at 18:00 UTC of DOY 1 2018.

Processing strategy
To evaluate the accuracy of SPD estimation using mapping functions and the horizontal gradient, as well as the hypothesis of non-isotropic SPD, three distinct strategies were employed in this study to estimate the SPD for five IGS stations at UTC 18:00 on DOY 001, 2018.Strategy 1 involved the utilization of three mapping functions: VMF1, VMF3, and GMF, to calculate the corresponding MF-SPD.VMF3 and VMF1 belong to the same series of mapping functions.VMF3, compared to VMF1, enhances the accuracy of estimating SPD at lower elevation angles.As the elevation angle decreases, the instability of SPD increases.Therefore, the selection of these two mapping functions from the widely used VMF mapping function series aims to examine the non-isotropic nature of SPD.The experimental procedure is outlined in Fig. 4.
The data preprocessing step is represented by the black box, followed by the calculation of SPD using VMF1 (yellow box), VMF3 (blue box), and GMF (green box), respectively.Initially, the position information of the IGS stations is retrieved, and the meteorological parameters required for VMF1 are obtained through GPT2w.Similarly, GPT3 is employed to acquire the necessary meteorological parameters for VMF3.Subsequently, the mapping function factors (mf i ) are computed based on these parameters.By employing the Saastamoinen model, the ZHD (Saas-ZHD) is determined, and the ZWD is obtained from the provided ZTD data by the IGS.The subscripts indicate the sources of the input meteorological parameters.The MF-SPD is then calculated by multiplying the ZHD and ZWD with the respective mapping function factors and then summing them.
In Strategy 2, the highest-accuracy MF-SPD obtained from Strategy 1 is further enhanced by incorporating  horizontal gradient correction.This correction utilizes the horizontal gradient parameters provided by the IGS.Strategy 3 involves the estimation of RT-SPD using a raytracing approach.For this purpose, a 2D ray-tracing program is developed.The program requires the station position information as well as the meteorological parameters from the ERA-5 reanalysis.Notably, this program is capable to retrieve RT-SPD for any point.

Accuracy assessment of MF-SPD
This study aims to assess the accuracy of MF-SPD by comparing the estimations with Strategy 1 and Strategy 3, thereby demonstrating the non-isotropic characteristics of SPD.Using station BJFS as an example, Figs. 5, 6, 7 depict the SPD profiles for the elevation angles ranging from 10° to 80° at UTC18 DOY001 2018.The plots present RT-SPD in blue, MF-SPD estimated using VMF1 in yellow, MF-SPD estimated using GMF in green, and MF-SPD estimated using VMF3 in purple.
From the figures, it is evident that MF-SPD, unlike RT-SPD, does not consider the intricate atmospheric conditions, but assumes a uniform tropospheric delay for different azimuth angles.The discrepancies between SPD estimations using different mapping functions are typically around 1-2 mm.As the elevation angle decreases, the disparities between RT-SPD and MF-SPD tend to increase.At station BJFS, the differences between MF-SPD and RT-SPD range from − 78 to 15 mm.Notably, Figs. 6 and 7 highlight a more pronounced non-isotropic behavior of SPD at the elevation angles below 30°.
To provide a comprehensive overview of the disparities in the RT-SPD estimations using different mapping functions at different stations, Table 2 presents the average residuals between the various MF-SPD and RT-SPD values at different azimuth angles for five stations at varying elevation angles.From Table 2, it can be observed that the general trend is that as the elevation angle decreases, the differences tend to increase.At lower elevation angles, station LAMA exhibits the smallest average difference with a minimum value of − 0.8 mm.Conversely, at higher elevation angles, station ABMF shows the largest average difference, reaching a maximum value of 152.6 mm.When the cutoff elevation angle is set to 30°, the maximum is 73.9 mm.The average residual values for VMF1, GMF, and VMF3 are 16.92 mm, 11.06 mm, and 10.25 mm, respectively.It can be observed that the SPD estimated with VMF3 is closer to the RT-SPD.These outcomes underscore the critical importance in considering azimuth angle variations in SPD during GNSS data processing.Neglecting these variations of SPD at different azimuths can lead to the errors from millimeters to decimeters.

Accuracy assessment of MF-SPD with horizontal gradient correction in VMF3 estimation
This section evaluates the impact of horizontal gradient correction on the accuracy of MF-SPD by comparing it with RT-SPD and investigate the presence of non-isotropic behavior in SPD.Based on the analysis in Sect."Accuracy assessment of MF-SPD", the results with VMF3 have a good agreement with those with the raytracing.Therefore, this section focuses on assessing the Calculation of meteorological parameters using GPT2W of meteorological parameters using GPT3 VMF1-mf h Saas-ZHD 2w

VMF1-SWD ZWD 2w
Location of IGS station accuracy of VMF3 with the introduction of horizontal gradient correction.Figures 8, 9, 10 illustrate the SPD plots at station BJFS for elevation angles ranging from 10° to 80° at UTC18 DOY001 2018.The blue curve represents RT-SPD, while the green curve represents MF-SPD with horizontal gradient correction.From Figs. 8, 9, 10, it is evident that the inclusion of horizontal gradient correction affects significantly the shape and symmetry of the SPD, leading to a departure from the circular shape centered at the origin.There is a global shift in the SPD values, and the values at different azimuth angles are no longer equal.Table 3 provides a quantitative analysis of the differences between MF-SPD and RT-SPD for various azimuth angles at each station, indicating that as the elevation angle decreases, the differences increase.The station LAMA has the smallest average difference of 0.8 mm, while ABMF exhibits the largest difference of 15.1 cm.When the cutoff elevation angle is set to 30°, the minimum difference is 2.70 mm, and the maximum can reach 9.26 cm, with an overall average difference of 2.04 cm.Despite the incorporation of the horizontal gradient term, the SPD errors do not    significantly decrease, and the maximum error remains at the decimeter level.
The limited improvement in SPD when adding the horizontal gradient term can be attributed to the simplistic assumption of isotropy in the horizontal gradient model.This model fails to capture the true variations of SPD at different azimuth angles.Its linear representation of the refractivity profile does not adequately account for the complex and non-isotropic nature of the atmosphere.As a result, the effectiveness of the horizontal gradient correction in PPP is constrained, especially at low altitude angles.
In summary, the insignificant improvement in SPD when incorporating the horizontal gradient in PPP can be attributed to the limitations of the horizontal gradient

Non-isotropic of SPD
From the curves of RT-SPD in the previous two subsections, it can be observed that the differences in tropospheric delays across different azimuth angles are commonly present, indicating non-isotropy rather than anisotropy.Taking the right panel of Fig. 10 as an example, with an elevation angle of 30° the tropospheric delays for azimuth angles of 180-360° are approximately equal, while there is a difference of around 0-1 cm among azimuth angles of 30-150°.Similarly, with an elevation angle of 10° (Fig. 10, right panel), the tropospheric delays for azimuth angles of 90-150° are approximately equal, whereas there is significant variation in tropospheric delays for other azimuth angles.Moreover, as the elevation angle increases, the differences in SPD between different azimuth angles at the same elevation angle become smaller.
To further investigate the hypothesis of non-isotropy in the troposphere, this study estimated the SPD for all azimuth angles for elevation angles of 10-50° for station BJFS station at 18:00 UTC on DOY001 in 2018.The SPD above the elevation angle of 50° was not considered due to minimal differences of SPD between azimuth angles.Additionally, the differences (∆) between each SPD and the average SPD at the corresponding elevation angle were calculated.The distribution of ∆ across different elevation angles and azimuth angles is illustrated in Fig. 11.
In Fig. 11, the radial axis represents the elevation angle, the azimuth angle is represented by the angular axis, and the color represents the values of SPD.Based on Fig. 11, considering the elevation angle range of 10°-20°, it can be observed that the tropospheric delay values exhibit some variation within four azimuth angle intervals: 45°-135°, 225°-315°, 315°-45°, and 135°-225°.More specifically, the values within the 45°-135° and 225°-315° intervals are similar, while the 315°-45° and 135°-225° intervals display noticeable disparities compared to the other intervals.
The maximum difference between the highest and lowest SPD values across different azimuth angles can reach up to 15 cm, with more pronounced discrepancies observed at lower elevation angles.These azimuth angle-related variations in tropospheric delay contribute to differences (∆) at the decimeter level.Consequently, disregarding the fluctuations in SPD across different azimuth angles, particularly at low elevation angles and during calculations involving long baselines, can impose limitations on the precision of the obtained solutions.

Conclusion
This paper puts forth the notion that the SPD exhibits non-isotropy at different azimuth angles.To validate this hypothesis, three mapping functions, NWP data, and the tropospheric products provided by IGS at five stations are utilized.With RT-SPD serving as the reference, the accuracy of MF-SPD is evaluated at eight elevation angles.Additionally, the introduction of horizontal gradient based on VMF3 estimation of SPD is examined, and the accuracy is analyzed for eight elevation angles and 36 azimuth angles.The following key findings emerge: (1) Compared to the SPD estimated with GMF and VMF1, the SPD estimated with VMF3 has a good agreement with RT-SPD, indicating its superior ability to calculate SPD.The minimum residual between VMF3-SPD and RT-SPD is 0.8 mm, while the maximum residual can reach 152.6 mm.When the cutoff elevation angle is set to 30°, the maximum is 73.9 mm, and the average is 14.08 mm.(2) The incorporation of horizontal gradient into MF-SPD does not lead to a significant reduction in residual when compared to RT-SPD, and the highest residual can still reach the decimeter level.Due to the inherent stochastic nature of water vapor movement, the application of horizontal gradient models is limited in accurately capturing the variations of SPD at low elevation angles and over extended periods.
(3) The analysis of SPD estimated with the ray-tracing reveals that, at the same elevation angle, the SPD exhibits overall anisotropic behavior across different azimuth angles.However, within specific ranges, it demonstrates isotropic characteristics.This finding provides evidence in support of the hypothesis that SPD is non-isotropic.Therefore, to simulate the signal path more accurately through the troposphere and improve positioning accuracy, it is necessary to consider the non-isotropic of SPD with different azimuth angles.
By proposing and investigating the concept of nonisotropy, this study introduces a new perspective that challenges traditional assumptions of isotropy and anisotropy in the troposphere.This recognition of the nonisotropy of tropospheric delay represents a departure from the norm.This innovative viewpoint illuminates the spatial intricacies and variations of atmospheric conditions, providing theoretical support for enhancing modeling accuracy and precision.
As future work, we plan to conduct a temporal and spatial analysis of the non-isotropic characteristics of the tropospheric delay.This analysis will investigate the influence of factors such as time and geographical location on these characteristics.Additionally, we plan to formulate a mechanistic definition for non-isotropic tropospheric delays using piecewise functions.Our goal is to establish a high-precision SPD model that closely reflects real atmospheric conditions.

Fig. 1
Fig. 1 Schematic of the Ray-tracing approach

Fig. 3
Fig. 3 Statistical results of the derived RT-SPD with three MF-SPDs, neglecting horizontal variations in SPD at global IGS stations

Fig. 11
Fig. 11 Distribution of ∆ at various azimuth angles and elevation angles

Table 1
Station latitude, longitude, and height

Table 2
Residuals between MF-SPD and RT-SPD at five stations

Table 3
Residual in MF-SPD estimation with the introduction of horizontal gradient correction at five stations