TMF: A GNSS Tropospheric Mapping Function for the Asymmetrical Neutral Atmosphere

: Tropospheric mapping function plays a vital role in the high precision Global Navigation Satellites Systems (GNSS) data processing for positioning. However, most mapping functions are derived under the assumption that atmospheric refractivity is spherically symmetric. In this paper, the pressure, temperature, and humidity ﬁelds of ERA5 data with the highest spatio-temporal resolution available from the European Centre for Medium-range Weather Forecast (ECMWF) were utilized to compute ray-traced delays by the software WHURT. Results reveal the universal asymmetry of the hydrostatic and wet tropospheric delays. To accurately represent these highly variable delays, a new mapping function that depends on elevation and azimuth angles—Tilting Mapping Function (TMF)—was applied. The basic idea is to assume an angle between the tropospheric zenith direction and the geometric zenith direction. Ray-traced delays served as the reference values. TMF coefﬁcients were ﬁtted by Levenberg–Marquardt nonlinear least-squares method. Comparisons demonstrate that the TMF can improve the MF-derived slant delay’s accuracy by 73%, 54% and 29% at the 5 ◦ elevation angle, against mapping functions based on the VMF3 concept, without, with a total and separate estimation of gradients, respectively. If all coefﬁcients of a symmetric mapping function are determined together with gradients by a least-square ﬁt at sufﬁcient elevation angles, the accuracy is only 6% lower than TMF. By adopting the b and c coefﬁcients of VMF3, TMF can keep its high accuracy with less computational cost, which could be meaningful for large-scale computing. Contributions: analysis, Funding acquisition, D.Z. and J.G.; Methodology, D.Z.; Project J.G.; Resources, W.M.; D.Z.; J.G.; Writing—original D.Z.;


Introduction
Tropospheric delay refers to the effect caused by the propagation of the radio signals among the neutral atmosphere, which can be divided into a hydrostatic part and a wet part [1]. Many regional or global tropospheric delay models have been built to reduce the tropospheric delay error, which can be divided into two categories, depending on whether meteorological factors are needed or not [2]. Models of the first category use pressure, temperature, and humidity as their input parameters, such as Hopfield [3], Saastamoinen [4], Davis [5], Baby [6], Ifadis [7], Askne, Nordius [8], and MSAAS [9]. If in-situ meteorological observations are not available, the standard atmosphere [10][11][12] or empirical meteorological models [13][14][15] may also be used in many GNSS data processing applications. The second category doesn't rely on meteorological measurements, such as UNB [16], MOPS [17], TropGrid [18], ITG [19], IGGTrop [20], and SHAtropE [21].
However, due to the irregular spatial and temporal distribution of water vapor, it is challenging to precisely model the wet part of the tropospheric delay. Thus it has been a commonly used strategy for Global Navigation Satellites Systems (GNSS) data processing to estimate the tropospheric zenith delay [22][23][24], especially for high precision applications [25,26]. The estimated Zenith Wet Delay (ZWD) can be converted into Precipitable Water Vapor (PWV) [27][28][29], and therefore GNSS meteorology has gradually become a fundamental and effective method for sounding the atmosphere under any weather condition. Barriot, et al. [30] proposed an approach based on perturbation theory, with the ability to separate eddy-scale variations of the wet refractivity.
The mapping function has been used to scale the slant delays from various elevation angles to the zenith direction. Consequently, the mapping function's accuracy has significant and direct impacts on the determination of the ZWD and station coordinate. Since Marini first proposed the continued fraction form [31], almost all modern mapping functions, including Ifadis [7], MTT [32], NMF [33], IMF [34], UNBabc [35], VMF [36], GMF [37], VMF1 [38], GPT2 [39], and VMF3/GPT3 [40], have taken it as their model expression. Each mapping function has two subtypes: the hydrostatic part and the wet part. The main difference among various mapping functions is the specific value of each coefficient.
However, the Marini concept mapping functions were built on the assumption of the neutral atmosphere's spherical symmetry [41][42][43], which can be clearly seen from the expression being independent on azimuth (will be discussed in Section 2). This assumption holds only approximately even for the troposphere's normal state, mainly due to the atmospheric bulge, high variation of tropospheric meteorological parameters such as water vapor and temperature. Therefore, such mapping functions would degrade the estimated ZWD and station height in the GNSS data processing. The tropospheric delay's horizontal gradients, including a North-South and an East-West component, have been used to model the tropospheric delay's anisotropy [32,[44][45][46]. The inclusion of gradient models can significantly improve the accuracy of slant delays [43], station positions [44,[47][48][49][50], zenith delays [51,52], and PWV [27,53]. Nevertheless, only total gradients [54][55][56] can be estimated in the GNSS data processing, since the hydrostatic and wet gradient mapping functions are very similar. Spherical harmonics were used by Zhang [57] (using ray-traced delays) and Zhang, et al. [58] (using GPS-derived delays) to replace the mapping function and gradients. However, many more unknown parameters have to be fitted for those approaches.
To overcome the shortcomings due to the assumption that atmospheric refractivity is spherically symmetric, we tested a new mapping function-TMF-where a concept of tilting the tropospheric zenith by an angle introduced by Gardner [59], Herring [32], Chen, et al. [44], Meindl, et al. [49] is utilized in this study. The TMF takes not only the elevation but also the azimuth as its input parameter. Ray-tracing [60] through Numerical Weather Model (NWM) is one of the most accurate approaches to obtaining tropospheric delays. Hence, ERA5 data [61] of the highest spatio-temporal resolution provided by the European Centre for Medium-range Weather Forecast (ECMWF) was adopted for computing ray-traced delays, using the software WHURT programmed in FORTRAN and developed by Zhang [2]. In the second part, we discuss some critical algorithms for ray-tracing. A detailed definition of the TMF is given. In the third part, we firstly investigate the asymmetry of the slant tropospheric delays. Then the coefficients of TMF are fitted by the Levenberg-Marquardt nonlinear least-squares method, using ray-traced tropospheric delays. Four fitting schemes were compared, with a different spatial resolution of NWM and different sampling on elevations and azimuths. The performance of TMF against mapping functions based on the VMF3 concept, without or with an estimation of gradient parameters, is presented in the results and discussion section. The summaries and conclusions are given in the last part.

Materials and Method
The flow chart of the study procedure is presented in Figure 1. Firstly, ERA5 data for several IGS stations containing the pressure, temperature, and humidity fields, were Remote Sens. 2021, 13, 2568 3 of 22 retrieved from the ECMWF. Then, ray-tracing was performed at various elevation and azimuth angles for each IGS station. Finally, the TMF was fitted using mapping factors calculated by ray-traced slant delays and zenith delays.
gave the relationship below: where d P is the partial pressure of dry gases (in hPa); w P is the water vapor pressure (in hPa); T is the temperature (in degree Kelvin); d Z is the compressibility factor for dry air; w Z is the water vapor compressibility factor; and i k (i = 1, 2, 3) are refractivity constants, of which different versions are suggested by various authors [63]. In this study, we set constants in accordance with most ray-tracing software [64] [9,65]. The water vapor pressure can be calculated using the following equation [66]: where ε is the ratio of the molar masses of water vapor and dry gases, q is the specific humidity, and P is the total pressure.
Assuming that air in the troposphere behaves as an ideal gas, Equation (2) becomes [4,64]: ; h N and w N are denominated as the hydrostatic and wet components of the refractivity, respectively [5]. The refractivity of the troposphere, N, is defined by: where n is the refractive index, which represents the ratio of the propagating speed of radio waves in a vacuum and in the troposphere, respectively. Smith and Weintraub [62] gave the relationship below: where P d is the partial pressure of dry gases (in hPa); P w is the water vapor pressure (in hPa); T is the temperature (in degree Kelvin); Z d is the compressibility factor for dry air; Z w is the water vapor compressibility factor; and k i (i = 1, 2, 3) are refractivity constants, of which different versions are suggested by various authors [63]. In this study, we set constants in accordance with most ray-tracing software [64], which are k 1 = 77.689 K/hPa , k 2 = 71.2952 K/hPa and k 3 = 375, 463 K 2 /hPa [9,65]. The water vapor pressure can be calculated using the following equation [66]: where ε is the ratio of the molar masses of water vapor and dry gases, q is the specific humidity, and P is the total pressure. Assuming that air in the troposphere behaves as an ideal gas, Equation (2) becomes [4,64]: and N w are denominated as the hydrostatic and wet components of the refractivity, respectively [5].

Tropospheric Delay
The tropospheric delay can be defined as the discrepancy between the propagation time of a GNSS signal in the troposphere and in a vacuum: where STD, SHD and SWD denote the slant total delay, slant hydrostatic delay, and slant wet delay, respectively; τ is the geometric delay.

ERA5 and Appropriate Land Cover Radius
The ERA5 reanalysis data produced by the ECMWF was used in this study. As the fifth generation of ECMWF atmospheric reanalysis of the global climate, ERA5 covers the period from 1979 to the present. It provides hourly estimates of atmospheric variables such as air pressure P, temperature T, specific humidity q, and geopotential Φ, at a horizontal resolution of 0.25 • (~31 km), on 137 model levels from the surface up to 0.01 hPa (~80 km). We downloaded some regional ERA5 data containing daily parameters {P, T, q, Φ} (in the format of grib) and found the size of a unit grid (0.25 • × 0.25 • ) is~15.5 Kb. There are 1,036,800 unit grids for a global file, and the file size can reach up to~15.5 Gb, which will lead to too significant a computational burden. Hence, it is necessary to find the appropriate longitude and latitude cover area of the ray-tracing for a specific station.
In Figure 2, according to the law of sines, we get: where R = M, M is the meridian radius of curvature, H max is the height of the neutral atmosphere top; θ is the cutoff angle. Therefore ϕ Δ can be calculated by: and λ Δ can be calculated by: Therefore ∆ϕ can be calculated by: and ∆λ can be calculated by: where r = N · cos ϕ is the radius of the parallel circle, N is the radius of curvature in prime vertical.
The data set chosen for our study is comprised of data from 12 IGS (International GNSS Service) stations distributed worldwide (see Figure 3). Therefore ϕ Δ can be calculated by: and λ Δ can be calculated by: is the radius of the parallel circle, N is the radius of curvature in prime vertical. The data set chosen for our study is comprised of data from 12 IGS (International GNSS Service) stations distributed worldwide (see Figure 3).

Eikonal Equation and 3D Ray-Tracing
With Maxwell's equations and considering a medium free of currents and charges, and a short vacuum wavelength, the Eikonal equation can be written as [70]: where S is the optical path length; → x is a position vector; S( → x ) is referred to as the Eikonal, which can be used to compute ∇S( → x ) to get the ray direction; n is the refractive index of a medium at the position → x . For a length L in a ray-based system, we can get: Which can be split into two first-order differential equations as: where → v denotes the tangent of the trajectory at the point → x . Firstly, the meteorological parameters at the point are interpolated using the adjacent grid points (see Section 2.1.5. Interpolation of meteorological parameters), and then the refractivity (see Equation (4)) is computed. The gradient of the refractivity can be obtained by two adjacent points along the ray in a spherical coordinate system (r, λ, ϕ) via Equation (13), and partial derivatives can be calculated analytically (detailed information can be found in Appendix B of Hobiger, et al. [70]).
where r is the radial distance, λ is the longitude, and ϕ is the latitude. Finally, the tropospheric delay can be integrated by Equation (6). For a given outgoing elevation angle, the ray-tracing will be performed three times by adding small changes on the initial elevation angle, and a quadratic function can fit delays at the given angle.

Interpolation of Meteorological Parameters
As shown in Figure 4, during the ray-tracing process, to obtain meteorological parameters of the point L with longitude λ, latitude ϕ and height H, we may first retrieve the parameters' values of its eight adjacent grid points from the ERA5 data. Since the heights of these grid points on the same model level are always different from each other, the interpolation can be performed vertically to the height H at first and then horizontally by a bilinear method.  In the vertical direction, temperatures and humidity can be interpolated by a linear method. For the total pressure P and water vapor pressure W P , the following formulas are utilized in the vertical interpolation [13,16,73]: In the vertical direction, temperatures and humidity can be interpolated by a linear method. For the total pressure P and water vapor pressure P W , the following formulas are utilized in the vertical interpolation [13,16,73]: Remote Sens. 2021, 13, 2568 7 of 22 where P 0 , T 0 , P W 0 are meteorological values at the start point; g 0 is the gravity acceleration for a point at latitude ϕ and height h, the formula for the normal gravity [74] is used here; R d is the mean specific gas constant for dry air; γ is the water vapor lapse rate parameter, and β is the temperature lapse rate.
Since ERA5 provides data hourly, and the GNSS receiver always logs observation every 1 or 30 s, linear time interpolation between two epochs can be carried out to calculate the meteorological parameters for the time of interest.

Height Transformation
ERA5 is based on a geopotential height system, whereas the ray-tracing calculation is performed on the geometrical (ellipsoidal) system. Therefore, it is necessary to transform heights of the ERA5 data into geometrical height using: where h is the ellipsoidal height; H is the orthometric height; N is the geoid undulation; Φ is the geopotential provided by ERA5; γ(ϕ, λ, h) is the mean gravity acceleration between the geoid and the point (ϕ, λ, h).

Asymmetsry Demonstration of Slant Delays
Tropospheric slant delays were ray-traced hourly with 1 • steps at elevations angles from 3 • to 90 • (the zenith direction), and 1 • steps on azimuths from 0 • to 359 • . Therefore there are 734,400 (86 × 360 × 24) grid points for a site each day. To show the spatial asymmetry of the SHDs and SWDs visually, we created skyplots in various directions by removing the average value over all azimuths for each elevation. We calculated the mean value, Root Mean Square (RMS), and range (max minus min) for ray-traced SHDs and SWDs at a given elevation angle θ according to the following formulas: where L i (θ) is the slant delay at the i-th azimuth of elevation angle θ, L(θ) represents the mean value of L i (θ), and n is the number of azimuths.

Tropospheric Mapping Function
In the space geodesy, we may use the general formula to represent the tropospheric delay as follows: where θ and ϕ are the elevation and azimuth angle, respectively; STD, SHD and SWD denote the slant total delay, slant hydrostatic delay, and slant wet delay, respectively; ZHD and ZWD stand for the zenith hydrostatic delay and zenith wet delay respectively; m h (θ, ϕ) and m w (θ, ϕ) are the hydrostatic and wet mapping functions, respectively. Almost all modern mapping functions, e.g., NMF [33], GMF [37], VMF1 [38] and VMF3 [40], take the continued fraction form [31] as their expressions: where a, b and c are coefficients. Nowadays, VMF1 is one of the most popular mapping functions in GNSS analysis. The "b" and "c" coefficients of VMF1 are given by empirical equations, whereas the "a" coefficient is determined epoch-wise (00, 06, 12, 18 UTC) from ray-traced delays at 3 • elevation. VMF3 is the upgrade version of VMF1, which gets an average "a" on eight equally spaced azimuth angles. The ray-tracing of VMF3 is performed through three different NWM types, resulting in three different versions of VMF3. ERA-Interim (EI) is only published in blocks every few months, not available in real-time; OP is made available at about 18 UTC for the previous day, and FC can be made available in real-time at about 09 UTC for the following day.
In Equation (22), mapping factors on all azimuths at a specific elevation angle are assumed to be equal. This assumption is not always valid, mainly due to the high spatial and temporal variation of the meteorological parameters, especially the water vapor. As a result, it may introduce non-negligible errors to the estimated ZWD parameters, the station height and receiver clock parameters.
Tilting the zenith direction in the mapping function by an angle is one possibility to represent a tropospheric gradient [32]. The representation given by Chen, et al. [44], which has been used by many scientific GNSS software such as GAMIT [26], and Bernese [25], can be written as: where G NS and G EW are North-South and East-West components of the tropospheric gradient parameters, respectively, and m f g = 1 sin θ·tan θ+C i is the gradient mapping function, C h = 0.0031 for the hydrostatic [44], C w = 0.0007 for the wet [44], and C t = 0.0032 for the total [32] gradient mapping function.
With Equation (23), gradients can be estimated. The gradients can be split into a hydrostatic and a wet part, which are the integration of the North and East directional derivatives of the "hydrostatic" and "wet" parts of the refractivity along with the altitude [30]. One should note that it is difficult to separate the hydrostatic and wet gradients in the GNSS analysis [75] due to the similarity of the two gradient mapping functions (see Figure 5).
Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 23 In Equation (22), mapping factors on all azimuths at a specific elevation angle are assumed to be equal. This assumption is not always valid, mainly due to the high spatial and temporal variation of the meteorological parameters, especially the water vapor. As a result, it may introduce non-negligible errors to the estimated ZWD parameters, the station height and receiver clock parameters.
Tilting the zenith direction in the mapping function by an angle is one possibility to represent a tropospheric gradient [32]. The representation given by Chen, et al. [44], which has been used by many scientific GNSS software such as GAMIT [26], and Bernese [25], can be written as:  (23), gradients can be estimated. The gradients can be split into a hydrostatic and a wet part, which are the integration of the North and East directional derivatives of the "hydrostatic" and "wet" parts of the refractivity along with the altitude [30]. One should note that it is difficult to separate the hydrostatic and wet gradients in the GNSS analysis [75] due to the similarity of the two gradient mapping functions (see Figure 5). However, since the spatial distribution of the hydrostatic delays and the wet delays are always different, it is reasonable to estimate gradient parameters respectively for the two delays, which could be performed based on ray-traced delays from NWM. However, since the spatial distribution of the hydrostatic delays and the wet delays are always different, it is reasonable to estimate gradient parameters respectively for the two delays, which could be performed based on ray-traced delays from NWM.

The Definition of TMF
The refractivity at sea level inevitably increases from the warmer regions toward the colder regions. In the lower troposphere, the surfaces of constant refractivity acquire a general slope [59]. One method that can be used to model azimuthal asymmetry is based on a "tilted" atmosphere assumption [32]. Meindl, et al. [49] gave the relations between the tropospheric and the geometric zenith direction to derive the gradient parameters' concept. In this study, we applied this relation to the mapping functions.
As shown in Figure 6, assuming that there is a small angle β between the tropospheric zenith direction Z and the geometric zenith direction z, we get: where z denotes the tropospheric zenith angle; z is the geometric zenith angle and z = 90 • − θ; ϕ 0 is the azimuth of Z with respect to z; ϕ represents the azimuth of the GNSS observation.  Introducing tropospheric elevation angle  θ as an argument for the mapp tions, and utilizing the Taylor series to approximate the slant delay linearly as:  Equation (25) only the first-order term of the Taylor series is kept, which m duce the linearization error. Therefore in the mapping function expression, we re β and 0 ϕ as two variables, which can be estimated together with coefficients i c , directly by mapping factors of ray-traced slant delays over zenith delay. Ther Introducing tropospheric elevation angle θ as an argument for the mapping functions, and utilizing the Taylor series to approximate the slant delay linearly as: where θ = 90 • − z; G NS and G EW are the tropospheric gradient parameters; and ∂m f ∂θ is the derivative of the arbitrary mapping function m f with respect to the elevation angle θ. Equation (23) may also be obtained from Equation (25) by using.m f = 1 sin θ . Equation (25) only the first-order term of the Taylor series is kept, which may introduce the linearization error. Therefore in the mapping function expression, we reserve the β and ϕ 0 as two variables, which can be estimated together with coefficients a i , b i and c i , Remote Sens. 2021, 13, 2568 10 of 22 directly by mapping factors of ray-traced slant delays over zenith delay. Therefore, the Equation (22) can be rewritten as: where , the unique value for θ when θ = 90 • is to avoid the multivalued problem at the geometric zenith direction, andm i θ are called TMF. With this new mapping function, which depends both on the elevation and azimuth angle, it would be unnecessary to estimate the gradient parameters.

Fitting of TMF
The mapping function can be expressed as: where RSD i (θ, ϕ) and RZD i are the ray-traced slant delays (i = h for SHD, and i = w for SWD) and zenith delays (i = h for ZHD and i = w for ZWD), respectively. According to Equation (26), the parameters to be estimated for the TMF are {β, ϕ 0 , a i , b i , c i }. The nonlinear fit was achieved by the Levenberg-Marquardt method.
Since the ERA5 data is very hard-disk-consuming and the ray-tracing is rather timeconsuming, several computation approaches were investigated. As listed in Table 1 , which were picked in such a way as to cover the whole elevation range [34], and 24 azimuths (15 • steps from 0 • to 345 • ). It should be noted that the ray-traced zenith delay RZD i was computed at the geometric zenith direction, not the tropospheric zenith direction. It would be rigorous to perform ray-tracing again at the fitted tropospheric zenith direction, computing the factor between the geometric zenith delay and the tropospheric zenith delay, and multiply it by the initial TMF. However, since the β is very small (with an observed mean value of 31" found by Meindl, et al. [49] and 25" in our experiments), the difference between the two delays is negligible. The mapping functions' coefficients were determined through a nonlinear least square adjustment over all selected elevation and azimuth angles. After that, mapping factors outputted from these fitted mapping functions were multiplied by the ray-traced zenith delay (ZHD or ZWD) to get the TMF-derived slant delays (SHD or SWD). Ray-traced slant delays by the first scheme (ERA5 with a horizontal resolution of 0.25 • × 0.25 • , at all 86 elevations and 360 azimuths) were used as truth values.
After the appropriate scheme was determined, a series of mapping functions that are based on the TMF, and the VMF3 concept, without or with gradient estimation, were designed and tested. As listed in Table 2, "WHU-" denotes that the ray-traced delays of WHURT serve as the basis for the calculation of mapping function coefficients, using the same model as VMF3 (see Equation (22)); LS means Least-Squares; eight equally spaced azimuth angles (0 • :45 • :345 • ) were chosen to be consistent with VMF3 [40]; VMF3 in the last three column means the value of the coefficient or gradients were retrieved from the VMF3_EI site-wise version.

Performance Test of TMF
In order to assess the errors of mapping functions effectively, their function values at some selected elevations and azimuths were multiplied by the ray-traced zenith delay to get the nominal slant delays; namely, the MF (mapping function)-derived slant delays. Then the RMS of differences between the ray-traced slant delays and the MF-derived slant delays were calculated by: where m = 86, n = 24 are the number of elevation angles and azimuths, respectively; SD ray ij is the ray-traced SHD or SWD, and SD MF ij is the MF-derived SHD or SWD. The improvement percentage of a mapping function against another mapping function was computed as: According to the rule of thumb [34,36], one-third of the mapping function delay error at the lowest elevation angle included in the analysis can be seen as station height error, which indicates the impact of the mapping functions on the position domain without processing GNSS observations.

Tropospheric Delay Asymmetry
The tropospheric delays' asymmetry can be assessed visually by skyplots with the removal of the average value over all azimuths on each elevation angle. Due to space limitation, only a few of them are present here exemplarily to demonstrate the spatiotemporal variability. Figures 7 and 8 Figure 7d shows much more significant anisotropy (please note the disparity in the bounds of the colour bars) than Figure 7c, which means that there was a quick variation of SWD from 0:00 to 5:00 UTC. This may be due to the fast-changing distribution of humidity, the typical summer weather conditions at Shanghai, where the SHAO station is located.
The situation is a little different on 26 December. As shown in Figure 8, although the elevation-dependent pattern is similar to Figure 7, SHD shows more spatial variability than SWD this time. The SHD range can reach up to ~10 cm at 5° elevation, while the SWD range is no more than several centimeters. This result may be caused by the fact that there is much less water vapor in winter than in summer in Shanghai. However, a comparison between Figures 8c,d shows that the temporal variations of SWD are still more complicated than SHD.  Figure 7d shows much more significant anisotropy (please note the disparity in the bounds of the colour bars) than Figure 7c, which means that there was a quick variation of SWD from 0:00 to 5:00 UTC. This may be due to the fast-changing distribution of humidity, the typical summer weather conditions at Shanghai, where the SHAO station is located.
The situation is a little different on 26 December. As shown in Figure 8, although the elevation-dependent pattern is similar to Figure 7, SHD shows more spatial variability than SWD this time. The SHD range can reach up to~10 cm at 5 • elevation, while the SWD range is no more than several centimeters. This result may be caused by the fact that there is much less water vapor in winter than in summer in Shanghai. However, a comparison In order to get a more quantitative investigation on the results at low elevations, we summarise the statistical result of slant delays at four specific elevations: 5°, 10°, 15° and 20° in Tables 3 and 4. There is much in common between the two tables. Firstly, the SHD and SWD and their range and RMS all tend to increase positively as elevation angle decreases. However, the SHD is always 6-15 times as large as the SWD. Secondly, the RMS of SHD and SWD at an elevation above 15° are mainly at the level of several millimetres. Furthermore, the RMS and range at the 5° elevation are always ten times as lager as that at the 20° elevation, both for SHD and SWD. Results above indicate that both the SHD and the SWD may present decametric asymmetry at low elevations.  In order to get a more quantitative investigation on the results at low elevations, we summarise the statistical result of slant delays at four specific elevations: 5 • , 10 • , 15 • and 20 • in Tables 3 and 4. There is much in common between the two tables. Firstly, the SHD and SWD and their range and RMS all tend to increase positively as elevation angle decreases. However, the SHD is always 6-15 times as large as the SWD. Secondly, the RMS of SHD and SWD at an elevation above 15 • are mainly at the level of several millimetres. Furthermore, the RMS and range at the 5 • elevation are always ten times as lager as that at the 20 • elevation, both for SHD and SWD. Results above indicate that both the SHD and the SWD may present decametric asymmetry at low elevations.

TMF Fitting
The results of the four fitting schemes introduced in Table 1 are listed in Table 5, in which elevations are divided into two bands: low elevation (3 • ≤ θ ≤ 15 • ) and high elevation (15 • < θ < 90 • ). As shown in Table 5, there is no apparent difference for bias and RMS between the two elevation and azimuth angle selection strategies (1 vs. 2, or 3 vs. 4). However, the horizontal resolution of ERA5 has a significant impact on the results. The RMS of the fitted SWDs based on ERA5 with 1 • × 1 • horizontal resolution is four and 7~9 times larger than that of the 0.25 • × 0.25 • , at low and high elevation angle bands, respectively. Results for SHD are similar but a little better. Hence, we use Scheme 2 in Table 1 (ERA5 with a horizontal resolution of 0.25 • × 0.25 • , at 18 selected elevations and 24 azimuths) to implement ray-tracing in the following research, which aims to keep a balance between the computational accuracy and the efficiency. Table 5. Statistic result of TMF-derived slant delays, which are the product of the TMF and the ray-traced zenith delays. The meaning of the postfix numbers is listed in Table 2.

Elevation
Angle

TMF Performance
The MF-derived slant delays are the production of the mapping factors and the raytraced zenith delay. The discrepancy compared with the ray-traced slant delays can directly reflect the accuracy of a mapping function. Figure 9 shows the RMS scatters of the discrepancies between MF-derived and raytraced delays at the 5 • elevation angle for the 12 globally distributed IGS stations listed in Table 1, at 96 epochs on four days (doy: 74, 202, 246, and 360) in 2018. Since each station's computation is independent of each other, such a graph would be an excellent way to reflect the global applicability of a mapping function. As shown in Figure 9, the TMFabc performs the best globally, with RMS of almost the same level for each station. Notably, the station FAIR and YAKT have the highest accuracy of wet delays. This may be due to the cold and dry weather on high latitudes in the northern hemisphere. reflect the global applicability of a mapping function. As shown in Figure 9, the TMFabc performs the best globally, with RMS of almost the same level for each station. Notably, the station FAIR and YAKT have the highest accuracy of wet delays. This may be due to the cold and dry weather on high latitudes in the northern hemisphere.
In contrast, there are much more variations with the symmetric mapping functions based on the VMF3 concept (WHU-VMF3a, WHU-VMF3a_gg, and WHU-VMF3a_g), especially for the wet delays. However, comparing WHU-VMF3a_gg with WHU-VMF3a shows that the estimation of gradients improves the result dramatically, although not as good as that of TMFs.
The performances of the original VMF3 (VMF3abc and VMF3abc_gg) are the worst but not surprising since the resolution of the NWM used by us (0.25° × 0.25°, 137 model levels, 1 h) is superior to that of VMF3 (1° × 1°, 25 pressure levels, 6 h [40]). VMF3abc_gg improves accuracy mainly on hydrostatic delays but very slightly on wet delays. We can conclude that the resolution of the NWM limits the accuracy of VMF3-derived wet gradients.  Tables 6 and 7 give the statistics results of the hydrostatic and wet part of all mapping functions (except for the WHU-VMF3a_g, of which the hydrostatic and wet delays cannot be separated), respectively. These two tables clearly show that the differences between various models are primarily in the RMS; the differences concerning the bias are minimal.
As shown in Table 6, TMFabc performs the best among all hydrostatic mapping functions, with the minimal RMSall (RMS for all elevation angles) and the minimal RMS5° (RMS for the 5° elevation angle). Compared with WHU-VMF3a, TMFabc improves the RMSall and the RMS5° by 68% and 74%, respectively. By the extra estimation of hydrostatic gradients, WHU-VMF3a_gg improves these by 62% and 65%, respectively. The accuracy of WHU-SMFabc_gg is only slightly lower than TMFabc. For VMF3, with gradients, the RMSall and the RMS5° can be improved by 35% and 39%, respectively. However, it seems that there are some slight systematic errors between the VMF3-derived and the WHURT- In contrast, there are much more variations with the symmetric mapping functions based on the VMF3 concept (WHU-VMF3a, WHU-VMF3a_gg, and WHU-VMF3a_g), especially for the wet delays. However, comparing WHU-VMF3a_gg with WHU-VMF3a shows that the estimation of gradients improves the result dramatically, although not as good as that of TMFs.
The performances of the original VMF3 (VMF3abc and VMF3abc_gg) are the worst but not surprising since the resolution of the NWM used by us (0.25 • × 0.25 • , 137 model levels, 1 h) is superior to that of VMF3 (1 • × 1 • , 25 pressure levels, 6 h [40]). VMF3abc_gg improves accuracy mainly on hydrostatic delays but very slightly on wet delays. We can conclude that the resolution of the NWM limits the accuracy of VMF3-derived wet gradients. Tables 6 and 7 give the statistics results of the hydrostatic and wet part of all mapping functions (except for the WHU-VMF3a_g, of which the hydrostatic and wet delays cannot be separated), respectively. These two tables clearly show that the differences between various models are primarily in the RMS; the differences concerning the bias are minimal. As shown in Table 6, TMFabc performs the best among all hydrostatic mapping functions, with the minimal RMS all (RMS for all elevation angles) and the minimal RMS 5 • (RMS for the 5 • elevation angle). Compared with WHU-VMF3a, TMFabc improves the RMS all and the RMS 5 • by 68% and 74%, respectively. By the extra estimation of hydrostatic gradients, WHU-VMF3a_gg improves these by 62% and 65%, respectively. The accuracy of WHU-SMFabc_gg is only slightly lower than TMFabc. For VMF3, with gradients, the RMS all and the RMS 5 • can be improved by 35% and 39%, respectively. However, it seems that there are some slight systematic errors between the VMF3-derived and the WHURT-derived hydrostatic slant delays, which would be mainly due to the difference in the NWMs. By adopting the b and c coefficients of VMF3, TMFa can achieve higher accuracy than WHU-VMF3a_gg, with less computational cost than TMFabc, which could be meaningful for large-scale computing.
In Table 7, all MF-derived wet delays' biases are closer to zero than the hydrostatic ones, but mostly with larger RMSs. TMFabc is still the most accurate, followed by the WHU-SMFabc_gg and TMFa. There is no significant difference between TMFabc and WHU-SMFabc_gg for most stations at most epochs, except for some particular conditions. Figure 10 shows this exemplarily for IGS station SHAO on the doy 202, 2018, when typhoon "Ampil" passed Shanghai. During 3:00-5:00 and 18:00-20:00 UTC, there are apparent improvements for TMFabc against WHU-SMF3abc_gg. WHU-VMF3a_gg improves the RMS all and the RMS 5 • of WHU-VMF3a by 59% and 62%, respectively.   Figure 11 shows the RMS scatters of all kinds of MF-derived total delays at the 5 • elevation angle, which are calculated by rms total = rms 2 h + rms 2 w , where rms h , and rms w are the RMS 5 • of the hydrostatic and wet MF-derived delays respectively. TMFabc performs the best. The improvement percentage of TMFabc against WHU-VMF3a, WHU-VMF3a_g, and WHU-VMF3a_gg can reach up to 73%, 54%, and 29%, respectively. Even by directly adopting the two coefficients b and c of VMF3 instead of estimating them, TMFa can still improve WHU-VMF3a, WHU-VMF3a_g WHU-VMF3a_gg by 68%, 47%, and 18%, respectively. It is worth noting that the improvement percentage of WHU-VMF3a_gg against WHU-VMF3a_g is 35%, which indicates that a separate estimation of the hydrostatic and wet gradients is preferable to a coarse estimation of the total gradients. Compared with WHU-VMF3a, the reduced RMS 5 • ratio between WHU-VMF3a_gg and TMFabc is 2.8/3.3, which is close to the value found by Landskron, et al. [46] that two-thirds of the azimuthal asymmetry can be described by the first-order gradients based on the VMF3 mapping function. VMF3a_g, and WHU-VMF3a_gg can reach up to 73%, 54%, and 29%, respectively. Even by directly adopting the two coefficients b and c of VMF3 instead of estimating them, TMFa can still improve WHU-VMF3a, WHU-VMF3a_g WHU-VMF3a_gg by 68%, 47%, and 18%, respectively. It is worth noting that the improvement percentage of WHU-VMF3a_gg against WHU-VMF3a_g is 35%, which indicates that a separate estimation of the hydrostatic and wet gradients is preferable to a coarse estimation of the total gradients. Compared with WHU-VMF3a, the reduced RMS5° ratio between WHU-VMF3a_gg and TMFabc is 2.8/3.3, which is close to the value found by Landskron, et al. [46] that twothirds of the azimuthal asymmetry can be described by the first-order gradients based on the VMF3 mapping function.

Discussion
Investigation of 360-degree ray-traced delays based on NWM revealed that both the SHD and SWD show azimuthal asymmetry, reaching up to decimeter level at low elevation angles below 15°. When a symmetric mapping function is used, the accuracy would not be as bad as expected, since the cutoff elevation angle is always set to 10° or 15° for most geodetic applications, which are typically based on double-differenced solutions. However, significant correlations between troposphere zenith delay parameters, station height, and receiver clock parameters are found [48]. The situation may be considerably improved by lowering the elevation cutoff angle. According to the rule of thumb, the delay path error at a 5° elevation will map to the station height at a ratio of 1/3 [34,36]. There-

Discussion
Investigation of 360-degree ray-traced delays based on NWM revealed that both the SHD and SWD show azimuthal asymmetry, reaching up to decimeter level at low elevation angles below 15 • . When a symmetric mapping function is used, the accuracy would not be as bad as expected, since the cutoff elevation angle is always set to 10 • or 15 • for most geodetic applications, which are typically based on double-differenced solutions. However, significant correlations between troposphere zenith delay parameters, station height, and receiver clock parameters are found [48]. The situation may be considerably improved by lowering the elevation cutoff angle. According to the rule of thumb, the delay path error at a 5 • elevation will map to the station height at a ratio of 1/3 [34,36]. Therefore, it is essential to model variations of the slant delays on the elevation and the azimuth.
The gradient parameters have been used as a supplement for a symmetric mapping function to model the azimuthal asymmetry for decades. However, in former GNSS studies [43][44][45]51], only the total gradients can be estimated, due to the difficulty of distinguishing the hydrostatic and wet gradient mapping functions. The WHU-VMF3a_g can be deemed a simulation for such tropospheric delay strategy, of which accuracy is not as good as that of separate estimation of hydrostatic and wet gradients (WHU-VMF3a_gg). In contrast, TMF can achieve better accuracy than a VMF3-like symmetric mapping function enhanced by estimation of respective gradients, of which coefficients b and c are not estimated but obtained from climatology data, and a is determined at a single elevation angle. The result of SMFabc_gg indicates that: (1) If a symmetric mapping function is not accurate enough, it would be difficult to precisely model the anisotropy part of the tropospheric delays by gradients. Therefore, estimating coefficients by least-squares on sufficient widespread elevation angles is necessary. (2) Assuming a proper mapping function is applied, gradients can describe the bulk of the tilting troposphere in most situations. However, there could still be some non-negligible high-order residuals [46], especially under some particular conditions such as severe weather scenarios [76].
The assumption of a tilted troposphere can explain most of the actual distribution of the Earth's troposphere. However, the tropospheric delay error of TMF is about ±1.2 cm at the 5 • elevation angle, corresponding to a station height error of ±4 mm. A tilted troposphere may not fully explain the variability of the tropospheric delays, which is also affected by other factors, such as the highly variable water vapor. Therefore, further study should be investigated on how to model the residual delays more precisely based on TMF, if higher accuracy is required. Besides, various NWMs featuring high spatio-temporal resolution produced by different agencies could be utilized together to validate each other and achieve a more robust mapping function.

Conclusions
In this research, ERA5 data retrieved from ECMWF with the highest spatial-temporal resolution was applied to compute the tropospheric slant delays by the ray-tracing method. It is found that azimuthal variation of tropospheric delay at low elevation angles can reach up to several decimeters. Traditional three order continued fraction mapping functions have been developed based on the assumption of atmospheric spherical symmetry, in which the azimuthal variation of tropospheric delays is neglected. To overcome shortcomings caused by such mapping functions, TMF is built by tilting the zenith direction of the mapping function, of which coefficients were fitted by Levenberg-Marquardt nonlinear least square method. We found that the horizontal resolution of NWM has a significant impact on mapping functions. However, it is not necessary to choose all elevations and azimuths. Fitting at 18 selected elevations and 24 azimuths are proved to be enough to produce comparable results with complete sampling.
The performances of several mapping functions based on the TMF and VMF3 concept were assessed quantitatively by using the ray-traced slant delays as reference values. Results show that TMF can improve 54% against the VMF3-like symmetric mapping function enhanced by estimating total gradients (WHU-VMF3a_g). According to the rule of thumb, a similar improvement in height parameter estimation is expected in GNSS analysis, which needs to be verified in our further study. A global grid-wise TMF can be developed for an arbitrary station's interpolation. To balance the accuracy and the computational cost, TMFa may be a good choice. Moreover, it would be necessary for the real-time application to derive TMF based on high-quality operational or forecast versions of NWM, including in-situ meteorological observations if possible.