The Performance of Different Mapping Functions and Gradient Models in the Determination of Slant Tropospheric Delay

Abstract

Global navigation satellite systems (GNSSs) have become an important tool for remotely sensing water vapor in the atmosphere. In GNSS data processing, mapping functions and gradient models are needed to map the zenith tropospheric delay (ZTD) to the slant total tropospheric delay (STD) along a signal path. Therefore, it is essential to investigate the spatial–temporal performance of various mapping functions and gradient models in the determination of STD. In this study, the STDs at nine elevations were first calculated by applying the ray-tracing method to the atmospheric European Reanalysis-Interim (ERA—Interim) dataset. These STDs were then used as the reference to study the accuracy of the STDs that determined the ZTD together with mapping functions and gradient models. The performance of three mapping functions (i.e., Niell mapping function (NMF), global mapping function (GMF), and Vienna mapping function (VMF1)) and three gradient models (i.e., Chen, MacMillan, and Meindl) in six regions (the temperate zone, Qinghai–Tibet Plateau, Equator, Sahara Desert, Amazon Rainforest, and North Pole) in determining slant tropospheric delay was investigated in this study. The results indicate that the three mapping functions have relatively similar performance above a 15° elevation, but below a 15° elevation, VMF1 clearly performed better than the GMF and NMF. The results also show that, if no gradient model is included, the root-mean-square (RMS) of the STD is smaller than 2 mm above the 30° elevation and smaller than 9 mm above the 15° elevation but shows a significant increase below the 15° elevation. For example, in the temperate zone, the RMS increases from approximately 35 mm at the 10° elevation to approximately 160 mm at the 3° elevation. The inclusion of gradient models can significantly improve the accuracy of STDs by 50%. All three gradient models performed similarly at all elevations and in all regions. The bending effect was also investigated, and the results indicate that the tropospheric delay caused by the bending effect is normally below 13 mm above a 15° elevation, but this delay increases dramatically from approximately 40 mm at a 10° elevation to approximately 200 mm at a 5° elevation, and even reaches 500–700 mm at a 3° elevation in most studied regions.

1. Introduction

Global navigation satellite systems (GNSSs) have become an important tool for remotely sensing water vapor in the atmosphere. The atmospheric information obtained from GNSS has provided a valuable complementary data source, not only for weather forecasting, but also for climate studies. The concept of global positioning system (GPS) meteorology was proposed by Bevis, et al. [1] to retrieve zenith wet delay (ZWD) and precipitable water vapor (PWV). ZWD represents the propagation delay in the zenith direction caused by the water vapor to the GNSS signal, while PWV is the depth of liquid water per unit area after precipitating all the water vapor over a given location at the zenith. PWV derived from GNSS has been proven to be accurate and reliable by comparing it to the value derived from water vapor radiometry (WVR) measurements [2], which is tuned to measure the microwave emissions of the vapor and liquid water molecules in the atmosphere at specific frequencies, with a root-mean-square (RMS) error of 1–2 mm [3,4]. Previous studies have shown that the GNSS-derived ZTD (the propagation delay of GNSS signal owing to the total effect of water vapor and dry air at the zenith), ZWD, and PWV [5,6] have a very high accuracy [7,8,9,10,11], as well as and the ability to capture the evolution of severe weather [12,13,14] and climate change [4,12,15,16,17,18], or can be used in operational weather forecasting by assimilating these variables into a numerical weather prediction (NWP) model [19,20]. However, ZWD and PWV are integrated values over a GNSS station and do not provide any information on the vertical distribution of water vapor, which means that severe weather evolution is not always captured by only monitoring ZWD or PWV [21].

Slant wet delay (SWD) or slant water vapor (SWV), which is the wet delay or the total precipitable water vapor in the slant direction, has the ability to capture the water vapor variation in both the horizontal and vertical directions and contains more information than that of the ZWD or PWV [21,22]. Ha, et al. [23] assimilated simulative GPS-derived SWD data, which excelled in the reconstruction of water vapor information in a hypothetical network of ground-based GNSS receivers and short-term rainfall prediction. Seko, et al. [24] compared water vapor’s vertical distribution from radial wind (RW) stemming from Doppler radar and SWV calculated by GPS, with RW stemming from Doppler radar and PWV calculated by GPS; the former improved the presentation of the vertical water vapor distribution. Bauer, et al. [25] demonstrated the remarkable improvement of nowcasting and short-range weather forecasting when assimilating GPS-derived slant total delays (STDs) into the numerical weather prediction model (NWM). Kawabata, et al. [26] indicated that assimilating STDs, rather than ZTDs or PWV, improves the recovery of the moisture and temperature fields and increases the accuracy of local heavy rainfall NWP forecasting.

Apart from assimilating STD, SWD, or SWV into NWP models, we can also reconstruct the 3D water variation with SWV using the tomography technique [27,28,29,30]. Flores, et al. [31] showed that the tomographic technique is a powerful tool for retrieving the 3D variation in tropospheric refractivity with the GPS observations from a regional dense network. Subsequently, numerous studies have been conducted to improve the performance of the tomographic technique in retrieving 3D water vapor fields [32,33,34,35,36,37,38].

The accuracy of the assimilated SWV determines the weather prediction performance and 3D tomography. Bender, et al. [39] also indicated that it is very important in GNSS tomography to use high-accuracy STDs with low elevations to reconstruct the boundary layer. In GNSS data processing, the STD cannot be directly estimated because too many unknowns would be involved. Usually, high-accuracy ZTDs can be derived with the double difference (DD) method, where some common error sources are removed. Additionally, ZTDs can be estimated with precise point positioning (PPP) [9,40,41,42,43,44,45,46,47].

In GNSS data processing, it is almost impossible to determine the tropospheric delay directly, because this will introduce too many unknowns into the equations. To reduce the unknowns of the tropospheric delay to be estimated, instead of estimating STD directly, the STD is usually represented as a function of ZTD (a piece-wise linear value), mapping functions, and gradient models. ZTD consists of zenith hydrostatic delay (ZHD) and ZWD, and the mapping function is divided into two components, hydrostatic and wet mapping functions [48]. Considering the effect of the asymmetry of the atmosphere on the ZTD estimation, various types of gradient models [49,50,51] have been introduced into GNSS data processing to further improve the accuracy of the STD. Therefore, slant hydrostatic delay (SHD) or SWD can be expressed as a function of ZHD or ZWD, respectively:

$$\mathrm{SHD}=\mathrm{ZHD}\times m_h\left(\epsilon \right)+G_h\left(G_{Nh},G_{Eh},\epsilon ,\alpha \right)+R_h,$$

(1)

$$\mathrm{SWD}=\mathrm{ZWD}\times m_w\left(\epsilon \right)+G_w\left(G_{Nw},G_{Ew},\epsilon ,\alpha \right)+R_w,$$

(2)

where $\alpha$ and $\epsilon$ are the azimuth and elevation, respectively; $m_h$ and $m_w$ are the hydrostatic and wet mapping functions, respectively; $G_h$ and $G_w$ represent the hydrostatic and wet gradient models, respectively; and $G_{Ni}$ and $G_{Ei}$ $\left(i=h,w\right)$ are the north–south and east–west gradient parameters, designated by indices of “h” (hydrostatic) or “w” (wet), respectively.

ZHD is acquired with a model for the surface pressure [52,53], and ZWD can be obtained through ZTD minus ZHD. In recent decades, based on the isotropy of the atmosphere, a mapping function has been developed that avoids directly ranking defects resulting from estimating STD with an arbitrary azimuthal angle and elevation angle [54]. The most commonly used mapping functions include CfA-2.2 [48], Ifadis [55], mapping temperature test (MTT) [56], NMF [57], GMF [58], and VMF1 [59], which were developed using either radiosonde observations or atmospheric reanalysis data.

To account for the effect of atmospheric anisotropy, many gradient models have been proposed for Very Long Baseline Interferometry (VLBI) and GNSS data processing. There have been few studies on the effect of gradient models during the past few decades. Furthermore, although hydrostatic and wet gradients are estimated separately in Equations (1) and (2), the total gradient, rather than the separate gradients, is currently used in postprocessing [60]. Therefore, three total gradient models are most commonly used today. These gradient models are applied in the GNSS post-processing software GAMIT [50], GIPSY [61], and Bernese [51].

As mentioned above, STDs derived from GNSS, especially those at low elevations, are very useful in retrieving 3D water vapor distributions or in weather forecasting by assimilating them into an NWP model. Therefore, it is important to investigate the accuracy of commonly used mapping functions and gradient models and their impact on the determination of STD. In this study, the performance of three mapping functions—NMF [57], GMF [58], and VMF1 [59]—and three gradient models—Macmillan [49], Chen and Herring [50], and Meindl, Schaer, Hugentobler and Beutler [51]—were studied using the atmospheric reanalysis dataset ERA-Interim provided by the European Centre For Medium-Range Weather Forecasts (ECMWF) [62]. A ray-tracing technique [63] was used to determine the STD, SHD, and SWD at nine elevations (i.e., 70°, 50°, 30°, 20°, 15°, 10°, 7°, 5°, and 3°) with the atmospheric information provided from ERA-Interim. The performance of the mapping functions and gradient models was assessed by comparing the SHD and SWD calculated using Equations (1) and (2) with the models obtained using the ray-tracing technique. Notably, the ZHD and ZWD used in Equations (1) and (2) were determined with the ERA-Interim dataset. To study the temporal and spatial properties of these mapping functions and gradient models, experiments were conducted in six regions, including the temperate zone, Qinghai–Tibet Plateau, Equator, Sahara Desert, Amazon Rainforest, and North Pole, for four seasons.

2. Data and Methodology

In this study, the ray-tracing technique was adopted to calculate the SHD and SWD at different elevations and azimuth angles with atmospheric information from the ERA-Interim dataset, which has a horizontal resolution of 0.75° × 0.75° and a vertical resolution of 37 levels, with a top level of 1 hPa, and is available at 0, 6, 12, and 18 Universal Time Coordinated (UTC) [64]. The pressure-level data are global grid data for geopotential, absolute temperature, specific humidity, and pressure (hPa). More details on the ray-tracing technique, mapping functions, and gradient models involved in this study are presented in the following subsections. Three mapping functions, NMF [57], GMF [58], and VMF1 [59], and three gradient models, including Macmillan [49], Chen and Herring [50], and Meindl, Schaer, Hugentobler and Beutler [51], are presented. The ray-tracing technique that was adopted to determine the propagation path of the electromagnetic signal considered as a ray is introduced in Section 2.3 [65].

2.1. Mapping Function

Marini [54] proposed a mapping function form with a continued fraction, which was then adopted by Niell [57], Böhm, Niell, Tregoning and Schuh [58], and Böhm, Werl and Schuh [59] by extending the form of the continued fraction to three terms. The mapping function is dependent on the elevation of the observation, which can be written as follows:

$$m_i\left(\epsilon \right)=\frac{1+\frac{a_i}{1+\frac{b_i}{1+c_i}}}{sin\epsilon +\frac{a_i}{sin\epsilon +\frac{b_i}{sin\epsilon +c_i}}}\left(i=h,w\right),$$

(3)

where $a_i$, $b_i$, and $c_i$ are different coefficients, designated by indices of “h” (hydrostatic) or “w” (wet). Each of the GMF, NMF, and VMF1 implementations follow this form by adopting different values for the coefficients.

For the NMF [57], the coefficients were determined using one-year meteorological data from 26 radiosonde stations mostly around the Northern Hemisphere, as well as the temperature and relative humidity profiles of the U.S. standard atmosphere for the northern latitudes of 15$°$, 30$°$, 45$°$, 60$°$, and 75$°$ in January and July. The coefficients of the NMF hydrostatic mapping function depend on the latitude, height of the site, and day of the year, while the coefficient of the wet mapping functions depends on the latitude.

Böhm (2006) published the GMF [58] considering the problem of the time delay of the VMF1 coefficients. Further, GMF has good agreement with VMF1, which is also an empirical function of the coordinates and day of the year. Moreover, compared with NMF, GMF has few annual errors and small regional height biases.

For the VMF [66], the coefficients $b_h$ and $c_h$ were determined based on the method of ray-tracing through the global grid data of ECMWF ERA40. VMF1 [59] is an update version of VMF [66]. Coefficient $b_h$ is assigned a new value, and coefficient $c_h$ is redefined as a function of the latitude and Julian day to correct the systematic errors of VMF related to climate zone and seasonality.

To investigate the effects of different mapping functions on the determination of STD, we first calculated STD using the following Equation (4):

$$STD\left(\epsilon \right)=m_h\left(\epsilon \right)·\mathrm{ZHD}+m_w\left(\epsilon \right)·\mathrm{ZWD}+\mathsf{\delta },$$

(4)

where ZHD and ZWD were calculated from the ERA-Interim dataset with Equation (28), and $\mathsf{\delta }$ is the residual, which is the unmodeled tropospheric delay.

2.2. Gradient Model

As mentioned above, the mapping function links the tropospheric delay in the zenith and slant directions in different directions without considering atmospheric asymmetry. The gradient models, which are a function of the azimuth, are an effective way to account for the atmospheric asymmetry effect.

Chen and Herring [50] and Herring [56] proposed a gradient model based on the “tilted” atmosphere assumption, which can be expressed as the first term of the gradient mapping function [67]:

$$G\left(G_N,G_E,\epsilon ,\alpha \right)=m_g\left(\epsilon \right)\left(G_{N}cos\alpha +G_{E}sin\alpha \right),$$

(5)

where G represents the gradient model; and $G_N$ and $G_E$ are the north–south and east–west gradient parameters, respectively.

$$m_g\left(\epsilon \right)=\frac{1}{sin\epsilon tan\epsilon +C^{\prime }}$$

(6)

where C is a constant of 0.0032 [50,56].

MacMillan (1995) [61] proposed a gradient model similar to the previous one by replacing $m_g\left(\epsilon \right)$ with $m\left(\epsilon \right)cot\epsilon$:

$$G\left(G_N,G_E,\epsilon ,\alpha \right)=m\left(\epsilon \right)cot\epsilon \left(G_{N}cos\alpha +G_{E}sin\alpha \right),$$

(7)

where $m\left(\epsilon \right)$ is a mapping function. Although the MTT dry mapping function was adopted as $m\left(\epsilon \right)$ in Macmillan’s study, Macmillan also indicated that no obvious changes were observed when a hydrostatic or wet mapping function was adopted.

Meindl, Schaer, Hugentobler and Beutler [51] proposed another gradient model using a zenith angle $z$ between the propagating path of satellite signals and the tropospheric zenith direction (i.e., the direction in which the tropospheric delay is at its minimum) to replace the angle $\tilde{z}$ between the propagating path of the satellite signals and the geometrical zenith direction to represent a “tilted” atmosphere:

$$G\left(G_N,G_E,\epsilon ,\alpha \right)=\frac{\partial f}{\partial z}\left(G_{N}cos\alpha +G_{E}sin\alpha \right),$$

(8)

where $z$ is the tropospheric zenith angle and $\frac{\partial f}{\partial z}$ is the derivative of the arbitrary mapping function with respect to the zenith angle $z$.

To investigate the effects of different gradient models on the determination of STD, we first calculated STD using the following Equation (9):

$$\mathrm{STD}\left(\epsilon ,\alpha \right)=m_h\left(\epsilon \right) \mathrm{ZHD}+m_w\left(\epsilon \right) \mathrm{ZWD}+G\left(G_N,G_E,\epsilon ,\alpha \right)+\mathsf{\Delta },$$

(9)

where G represents the gradient model; G_N and $G_E$ are the north–south and east–west gradient parameters, which need to be estimated together with the zenith delay in GNSS data processing; and $\mathsf{\Delta }$ is the residual, which is the unmodeled tropospheric delay.

2.3. Ray-Tracing

Ray-tracing is a technique for reconstructing the ray path with a three-dimensional refractivity field based on geometrical optics theory. Although the 3D ray-tracing method can provide a very accurate result for the retrieved ray path, it is time-consuming in an elliptical coordinate frame (i.e., WGS84) [68,69]. Normally, the 3D ray-tracing method can be reduced to a 2D method by fixing the azimuth of the ray path. Compared woth the 3D method, 2D ray-tracing can produce a comparable accuracy with much less computation complexity. Considering the complexity of the ray-tracing implemented in an elliptical coordinate, an approximation of the Euler radius can be carried out in a local spherical coordinate, which can produce a high ray-tracing result with much less calculation time [63,69,70].

According to the Euler formula, the Euler radius of the local spherical coordinate is written as

$$R_{\alpha }=\frac{MN}{M{sin}^2\alpha +N{cos}^2\alpha },$$

(10)

where $\alpha$, $M,$ and $N$ are the azimuth, the meridian, and the prime radii of curvature, respectively. $M$ and $N$ are defined by the latitude $\phi$ at the starting point $P_1$ (see Figure 1) and are given by

$$M=\frac{a\left(1-e^2\right)}{{\left(1-e^2{sin}^2\phi \right)}^{\frac{3}{2}}},$$

(11)

$$N=\frac{a}{{\left(1-e^2{sin}^2\phi \right)}^{\frac{1}{2}}},$$

(12)

where $a$ and $e$ are the semi-major axis and the first eccentricity of the reference ellipsoid (i.e., WGS84) [71], respectively.

The piecewise-linear (PWL) approach adopted in this study is a simple and effective 2D ray-tracing algorithm [72]. As shown in Figure 1, the $y–z$ plane, where the PWL approach is implemented, is defined by the geometric center of the Earth $O$, the starting point $P_1$, and the puncture point of the electromagnetic signal and the top troposphere $P_k$. The $z$-axis joins the geometric center, the origin point, of the Earth with the starting point $P_1$. The $y$-axis is orthogonal to the $z$-axis and restricted to the observation direction defined by the azimuth angle $\alpha$. The ${\theta }_i$, ${\eta }_i$, and $e_i$ for $i=1,2,\dots ,k$ denote the tangent angles defined by the ray trajectory and the height layer, the geometric angles to the ray points $P_i$, and the elevation angles at the height layer, respectively. The refractive index at the i-th layer is indicated as n_i. Here, the subscript $i$ of $i=1,2,\dots ,k-1$ is the layer index of the PWL approach and the height layer. In the loop of the PWL approach, based on the geometric rules and Snell’s Law, the calculations are as follows:

$$s_i=-r_{i}sin{\theta }_i+\sqrt{r_{i+1}^2-r_i^2{cos}^2{\theta }_i},$$

(13)

$$z_{i+1}=z_i+s_{i}sin{e}_i,$$

(14)

$$y_{i+1}=y_i+s_{i}cos{e}_i,$$

(15)

$${\eta }_{i+1}=arctan\frac{y_{i+1}}{z_{i+1}},$$

(16)

$${\theta }_{i+1}=arccos\left(\frac{n_i}{n_{i+1}}cos\left({\theta }_i+{\eta }_{i+1}-{\eta }_i\right)\right),$$

(17)

$$e_{i+1}={\theta }_{i+1}-{\eta }_{i+1},$$

(18)

with the initial conditions ${\theta }_1=e_1$, $z_1=r_1$, $y_1=0$, and ${\eta }_1=0$. In the first layer of ray-tracing, the elevation angle $e_1$ at starting point $P_1$ is determined based on the following [63]:

$$e_1=\epsilon +\mathsf{\Delta }{g}_{apbend},$$

(19)

where $\mathsf{\Delta }{g}_{apbend}$ is the a priori value of the bending angle, which can be calculated with the following [69]:

$$\mathsf{\Delta }{g}_{apbend}\left(\epsilon ,h_0\right)\left[deg\right]=\frac{C_{b}exp\left(-\frac{h_0}{6000}\right)}{tan\epsilon },$$

(20)

where $h_0$ is the ellipsoidal (geodetic) height of the starting point (in meters), and $C_b$ is the empirical constant set to 0.02.

Before starting the ray-tracing calculation, the refractive indices $n_i$ used in Equation (17) must first be defined from an atmospheric model (e.g., the ECMWF meteorological data). The refractive indices $n_i$ or the refractivity $N$ can be calculated with Equation (21), proposed by Thayer [73]:

$$N={10}^6\left(n-1\right)=k_1\frac{P_d}{T}{Z}_d^{-1}+k_2\frac{P_w}{T}{Z}_w^{-1}+k_3\frac{P_d}{T^2}{Z}_w^{-1},$$

(21)

where $P_d$ and $P_w$ are the pressure of dry air and water vapor (hPa), respectively; $T$ is the temperature in Kelvin; Z_d and Z_w are the compressibility factors for dry air and water vapor, respectively, which can be calculated with

$$Z_d^{-1}=1+P_d\left(57.97·{10}^{-8}\left(1+\frac{0.52}{T}\right)-9.4611·{10}^{-4}\frac{T-273.15}{T^2}\right),$$

(22)

$$\begin{array}{c}{Z}_w^{-1}=1+1650\frac{P_w}{T^3}(1-0.01317\left(T-273.15\right)+1.75\times {10}^{-4}{\left(T-273.15\right)}^2+1.44\\ \times {10}^{-6}{\left(T-273.15\right)}^3), \end{array}$$

(23)

respectively, as suggested by Owens [74].

The $k_1$,$k_2,$ and $k_3$ in Equation (21) are three empirical constants, which were set to the “best average” values, 77.6890 (K/hPa), 71.2952 (K/hPa), and 375,463 (K²/hPa), respectively, as calculated by Rüeger [75] and suggested by Nafisi, et al. [76].

According to Böhm, et al. [77] and Wallace and Hobbs [78], the water vapor pressure can be expressed as

$$P_w=\frac{qp}{\psi +\left(1-\psi \right)q},$$

(24)

where p is the atmospheric pressure; and $\psi$ is the ratio between the gas constant of dry air $R_d$ and the gas constant of water vapor $R_w$, defined as $\psi =R_d/R_w=0.622$.

As mentioned above, a geodetic reference system is needed for the ray-tracing algorithms adopted in this study; therefore, the geopotentials in the ERA-Interim dataset needed to be transformed to ellipsoidal heights following the procedures given in previous studies (see the works of [77,79,80]).

First, the geopotential height $h_d$ can be calculated with

$$h_d=\frac{Z}{g_n},$$

(25)

where $g_n$ is a constant gravity value and was set to 9.80665 $\mathrm{m}/{\mathrm{s}}^2$ [77].

Then, the orthometric (geometric) height can be calculated as follows:

$$\begin{array}{l}{h}_{orth}\\ =\frac{1}{2·1.57·{10}^{-7}}\\ -\sqrt{\frac{1}{{\left(2·1.57·{10}^{-7}\right)}^2}-\frac{1}{1.57·{10}^{-7}}\frac{h_d}{1-0.0026373cos\left(2\phi \right)+0.0000059{cos}^2\left(2\phi \right)}}.\end{array}$$

(26)

Finally, the ellipsoidal height $h_{ell}$ can be determined by Fotopoulos [81]

$$h_{ell}=h_{orth}+h_N,$$

(27)

where $h_N$ is the geoid undulation, which was obtained from the Earth Gravitational Model 2008 [80].

2.4. Methods

Two sets of STDs, SHDs, and SWDs are calculated in this study. One set is derived by applying the ray-tracing program to the ERA-Interim dataset to determine the slant delays at different elevations and azimuths. Using the equations given in Section 2.3, the ray-tracing-derived STD, SHD, and SWD can be calculated, according to the work of [72], as

$$\mathrm{STD}={\displaystyle \sum }_{i=1}^{k-1}\left[\left(n_i-1\right)s_i\right]+g_{bend},$$

(28)

$$\mathrm{SHD}={\displaystyle \sum }_{i=1}^{k-1}\left[\left(n_{h,i}-1\right)s_i\right],$$

(29)

$$\mathrm{SWD}={\displaystyle \sum }_{i=1}^{k-1}\left[\left(n_{w,i}-1\right)s_i\right],$$

(30)

where $g_{bend}$ is the tropospheric delay caused by the bending effect calculated with the following Equation (31) [72]:

$$g_{bend}={\displaystyle \sum }_{i=1}^{k-1}\left[s_i-cos\left(e_i-e_k\right)s_i\right],$$

(31)

where $g_{bend}$ is the tropospheric delay caused by the bending effect; and $n_h$ and $n_w$ are the hydrostatic and wet indices, respectively, which were determined, following the works of [48] and [72], as

$$n_h=1+\left(k_1\frac{R}{M_d}\right)\times {10}^{-6},$$

(32)

$$n_w=1+\left(k_2^{\prime }\frac{P_w}{T}+k_3\frac{P_w}{T^2}\right)\times {10}^{-6},$$

(33)

where $k_2^{\prime }$ is the refractivity empirical constant defined as $k_2^{\prime }=k_2-k_1\psi$; and $R$ and $M_d$ are the universal gas constant and molar mass of dry air, set to 8314.510 $\mathrm{J}/\mathrm{kmol}·\mathrm{K}$ and 28.9644 $\mathrm{g}/\mathrm{mol}$, respectively [82,83]. In addition to the slant delay, the zenith equivalent delay (ZTD, ZHD, and ZWD) can be calculated at a 90$°$ elevation using Equations (28)–(30).

The other set uses ray-tracing to calculate the ZHD and ZWD first and then calculates the SHD and SWD using Equations (28)–(30). The STDs, SHDs, and SWDs calculated by ray-tracing directly are regarded as the benchmark values in this study to assess the accuracy of the slant delays determined with Equations (4) or (9). To investigate the performance of these mapping functions and gradient models for different regions and seasons, a comparison was conducted in the six typical regions listed in

Table 1. The six regions and coordinates in the scheme (the notations “N”, “S”, “E”, and “W” indicate northern latitude, southern latitude, eastern longitude, and western longitude, respectively).

Region	Latitude Range	Longitude Range	Number of Grids	Height of Grids (m)
Temperature Zone	33° N–36° N	115° E–117° E	9	200
Qinghai–Tibet Plateau	29.5° N–32.5° N	88.5° E–91.5° E	9	4500
Equator	1.5° S–1.5° N	144° E–147° E	9	200
Sahara Desert	18° N–21°N	1.5° W–1.5° E	9	200
Amazon Rainforest	2.5° S–5.5° S	66° E–69° E	9	200
North Pole	85.5° N–90° N	1.5° W–1.5° E	10	200

for January, April, July, and October. For each region, the ray-tracing calculations were run at nine or ten selected grid points at eight azimuths (0$°$, 45$°$, 90$°$, 135$°$, 180$°$, 225$°$, 270$°$, and 315$°$) and 10 elevations (90$°$, 70$°$, 50$°$, 30$°$, 20$°$, 15$°$, 10$°$, 7$°$, 5$°$, and 3$°$). As the performance of mapping functions are closely related to regional climate characteristics, especially the variation of water vapor, the comparison studies were conducted in the abovementioned six regions with very different climate conditions. It should be noted that the height of the grids in the Qinghai–Tibet region is set to 4500 m. This is because, as the world’s largest and highest plateau, the Qinghai–Tibet region’s average elevations exceeds 4500 m.

3. Results

In this study, the effects of various mapping functions and gradient models on STD determinations were studied and discussed in detail. In Section 3.1, the variation of the water vapor in six regions was studied. In Section 3.2, the performances of the three mapping functions for the different seasons in the six selected regions were compared. In Section 3.3, the performances of the three gradient models were investigated, and the improvements in the STD estimates with or without the gradient model were studied. In Section 3.4, the bending effect of the GNSS signal propagation, which is normally ignored in GNSS remote sensing or positioning, was studied.

3.1. Variation of PWV

As the performance of mapping functions is closely related to water vapor variation, the temporal characteristics of PWV in these six regions are studied in this subsection. Figure 2 shows the variation of PWV during 2018. It can be clearly seen that the annual variation of PWV is much larger in the temperate zone and Sahara Desert than that in the Qinghai–Tibet plateau and the North Pole. Although the PWV value is very high in the Equator and the Amazon Rainforest region, the annual variation is not very strong. The black lines in Figure 2 represent the fitting model obtained with the following Equation (34):

$$\begin{array}{c}\mathrm{PWV}=a_0+a_{1}cos\left(2\pi \frac{doy}{365.25}\right)+a_{2}sin\left(2\pi \frac{doy}{365.25}\right)+a_{3}cos\left(4\pi \frac{doy}{365.25}\right)\\ +a_{4}sin\left(4\pi \frac{doy}{365.25}\right),\end{array}$$

(34)

where $doy$ is the day of the year; $a_0$ is the mean value; and ($a_1$, $a_2$) and ($a_3$,$a_4$) are the annual and semi-annual amplitude, respectively.

As shown in Figure 2 and

Table 2. The mean value, annual amplitude, and semi-annual amplitude of precipitable water vapor (PWV) in the six regions.

Region	Annual Mean PWV (mm)	Annual Amplitude (mm)	Semi-Annual Amplitude (mm)
Temperature Zone	23.17	19.98	5.62
Qinghai–Tibet Plateau	6.40	7.22	2.64
Equator	53.34	0.62	1.03
Sahara Desert	21.57	5.36	14.68
Amazon Rainforest	44.85	1.81	2.96
North Pole	5.23	4.05	1.38

, the largest mean value of PWV is found in the region of the Equator, with a value of 53.34 mm, and the smallest value is found in the North Pole, with a value of only 5.23 mm. The annual amplitude of PWV in the temperate zone is as large as 23.17 mm, while it is only 0.62 mm in the Equator region. However, the PWV in the Equator shows a very strong monthly variation.

3.2. Mapping Function

Three widely used mapping functions (i.e., NMF, GMF, or VMF1) were adopted in this study, and the STDs at the abovementioned nine elevations were calculated with Equation (4) and compared to the benchmark values determined using the ray-tracing program.

As shown in Figure 3, both the RMS and bias show a marked dependence on the elevations in all six regions. The RMS and bias are calculated with

$$\mathrm{RMS}=\sqrt{\frac{1}{N}{\displaystyle \sum }_{i=1}^N{\left(ST{D}_{ray-tracing}-ST{D}_{M+\left(G\right)}\right)}^2},$$

(35)

$$\mathrm{bias}=\frac{1}{N}{\displaystyle \sum }_{i=1}^N\left(ST{D}_{ray-tracing}-ST{D}_{M+\left(G\right)}\right),$$

(36)

where $ST{D}_{ray-tracing}$ and $ST{D}_{M+\left(G\right)}$ are the STD calculated by Equations (4), (9), and (28)–(30); the $RM{S}_M$ represents the RMS error of the STD determined with a mapping function only; and $RM{S}_{G+M}$ represents the RMS error of the STD determined with a mapping function and gradient model, which is studied in the following subsection.

For the three mapping functions, the RMSs of the STD are all smaller than 2 mm above the 30$°$ elevation and smaller than 9 mm above the 15$°$ elevation, but show a significant increase below 15$°$. For example, in the temperate zone, the RMS value increases from approximately 35 mm at the 10$°$ elevation to approximately 160 mm at the 3° elevation. The value of the RMS also shows an obvious dependence on the regions. In the Qinghai–Tibet Plateau region, the largest RMS value is approximately 50 mm, while in the other regions, the largest values are all greater than 100 mm. This result might be because the variation in the atmosphere in the Qinghai–Tibet Plateau region is not as great as that in the other regions, and thus can be modeled more accurately. A similar phenomenon can also be observed in the North Pole region, where the water vapor content is much lower than that in other regions, such as the Equator regions. In terms of RMS, these three mapping functions show similar performance above a 15$°$ elevation. However, for the elevations below 15$°$, compared with the GMF and NMF, the VMF1 mapping function performs better, especially below the 7$°$ elevation. The NMF has a much larger RMS at low elevations than at high elevations, and the RMS was two times greater than that of VMF1 in the North Pole region at 3$°$ elevations.

In terms of both the RMS and bias, above 15$°$, the mapping functions considered in this study perform relatively similarly in the selected six regions. However, at low elevations, the performances of these three mapping functions are substantially different, and VMF1 has a higher accuracy than those of the other two mapping functions.

Table 3. The biases of the slant total delays (STDs) from Vienna mapping function (VMF1) with respect to the STDs from ray tracing at different azimuths.

Azimuth	Temperate Zone	Qinghai–Tibet Plateau	Equator	Sahara Desert	Amazon Rainforest	North Pole
0°	101.1	72.4	60.8	118.1	61.8	−7.8
90°	1.2	17.1	18.9	34.6	−17.0	−9.7
180°	−49.1	−0.3	41.5	13.9	−10.1	−24.2
270°	−7.8	16.1	11.8	30.8	−11.6	−11.7

shows that the biases at different azimuths have clear discrepancies, especially in the north–south direction, with biases exceeding those in the east–west direction. Thus, at low elevations, the biases do not fluctuate around the value of zero. The biases of NMF and GMF for the Qinghai–Tibet Plateau and the biases at the North Pole are lower than zero, which may be because of the lack of enough meteorological data for the region in the ECMWF reanalysis data.

Figure 4 shows the RMS and bias of STD derived with VMF1 in January, April, July, and October. For the regions located in the Northern Hemisphere (i.e., the temperature zone, the Sahara Desert, and the Qinghai–Tibet Plateau), it is obvious that the RMS values and biases are larger during the summer than during the winter. For the Equator region, the RMS value is obviously greater in January and October than in April and July. For the North Pole, the RMS errors are larger in July and October than in January and April. This result confirms that the performance of VMF1 is highly correlated with both the season and location.

3.3. Gradient Model

Three widely used gradient models were discussed in Section 2.2, and their performance will be studied in this section.

In Equation (9), the STDs at the nine elevations and six azimuths were calculated using the ray-tracing program, after which the $G_N$ and $G_E$ were estimated using the least-square method. Using the ray-tracing-derived STD as a reference, the RMSs and biases of the STD derived with the mapping function plus the gradient model are shown in Figure 5. Notably, VMF1, which is more accurate than the other two mapping functions, was adopted in Equation (9).

The blue bar and yellow line are the RMS and bias, respectively, of the STD determined with VMF1 without the gradient model. All three gradient models studied show very similar performances in the six selected regions.

A statistic of the improvement in STD accuracy was calculated with the following Equation (37), and the results are shown in

Table 4. The accuracy improvement in VMF1 owing to the adoption of the Meindl gradient model.

Elevation	Temperate Zone (%)	Qinghai–Tibet Plateau (%)	Equator (%)	Sahara Desert (%)	Amazon Rainforest (%)	North Pole (%)
3°	59.11	36.64	42.78	39.92	66.73	51.65
5°	62.48	52.20	47.84	49.83	73.41	60.67
7°	62.84	57.70	50.46	54.00	75.97	62.17
10°	63.05	60.22	54.27	58.57	77.65	61.22
15°	62.87	64.86	56.85	61.18	77.91	60.53
20°	62.75	65.62	57.54	61.62	77.76	61.08
30°	62.21	64.44	57.68	61.51	78.18	61.30
50°	61.07	56.62	51.44	55.46	67.26	53.19
70°	57.54	38.89	51.59	54.17	68.41	36.17

.

$$\mathrm{Accuracy}\mathrm{Improvement}=\frac{RM{S}_M-RM{S}_{G+M}}{RM{S}_M},$$

(37)

where $RM{S}_M$ indicates the RMS of the STD determined with the VMF1 mapping function only, while $RM{S}_{G+M}$ is derived with the VMF1 mapping function together with the Meindl gradient model.

As the performances of the three mapping functions are similar, only the improvement results for the Meindl, Schaer, Hugentobler and Beutler [51] gradient model are shown in Table 4. From Table 4, we find that the accuracy improvement ranges from 37% to 79% owing to the adoption of the Meindl, Schaer, Hugentobler and Beutler [51] gradient model. Therefore, the inclusion of a gradient model is useful for reducing STD estimation errors. However, the STD error is still as large at 10–20 mm at a 5$°$ elevation owing to the unmodeled atmospheric inhomogeneity. Notably, as the ZHD and ZWD in this study are calculated from the atmospheric reanalysis data, the errors in the ZHD and ZWD can be neglected. The STD error discussed in this study can be regarded as an error caused by an error in the mapping function and gradient model. In terms of bias, the inclusion of a gradient model in the estimation of STD seems to have no obvious benefit in reducing the bias.

Figure 6 presents the RMSs and biases for the six regions in the four seasons and indicates that the RMS errors and biases have an obvious seasonal variation in the temperate zone, the Qinghai–Tibet Plateau, and the Sahara Desert. The seasonal variations in RMS errors are similar to the findings stated in Section 3.1; that is, both the RMSs and biases are larger in the summer than in the winter. This phenomenon is most obvious in the temperate zone, where the RMS error in July is approximately 2–3 times that in April and October and approximately 3–4 times that in January. Figure 6 also shows that the biases are all positive, except in the North Pole region. This result indicates that the ray-tracing-derived STD is normally larger than the value obtained with the ZHD and ZWD mapping functions and the gradient model. One possible reason for these positive biases is the bending effect, which will be discussed in the following subsection.

Figure 7 illustrates that, in most regions, the error of SWD is larger than that of SHD.

3.4. Bending Effect

As shown in Figure 1, a neutral atmosphere will affect the GNSS signal in two ways: (1) the signal will travel at a speed less than that in a vacuum, and (2) the signal will be bent owing to atmospheric refraction caused by the vertical changes in atmospheric density [84,85]. To investigate the effect of the bending on determining the STD, the delay caused by the bending effect in the six regions for the different months was calculated by Equation (31). As shown in Figure 1, i denotes the i-th height layer; $s_i$ is the distance of the ray point $P_i$; and the next ray point, $P_{i+1}$, $e_1$, is the actual starting elevation angle; while $e_k$ is the actual ending elevation angle. All these values can be determined during the ray-tracing procedure.

Figure 8 shows that the tropospheric delay caused by the bending effect is normally below 13 mm above a 15$°$ elevation. However, this value increases dramatically from approximately 40 mm at a 10$°$ elevation to approximately 200 mm at a 5$°$ elevation and even reaches 500–700 mm at a 3$°$ elevation in most regions, which means that this effect becomes nonnegligible at low elevations, especially tropospheric tomography. The bending effect on the Qinghai–Tibet Plateau is much smaller than that on the other regions because the height of the Qinghai–Tibet Plateau is normally above 4 km, and thus the tropospheric delay is smaller than that in the other regions because of a shorter signal propagation travel distance in the troposphere. The results in Figure 8 also indicate that the bending effect does not have a clear dependence on season.

4. Discussions

The performance of three mapping functions and three gradient models was studied comprehensively in this study. As suggested in the study, at low elevations, the performances of these three mapping functions are substantially different, and VMF1 has higher accuracy than that of the other two mapping functions. This means that, in the GNSS data processing, if the cutoff angle was set to 15$°$, then all three mapping functions would yield a similar result. However, if the cutoff angle was set to a value less than 15$°$, then the VMF1 mapping function would always be recommended in GNSS data processing. The performance of the mapping functions shows a very strong dependence on the location and the season. The RMS values and biases are obviously larger during the summer than during the winter. This is because the water vapor is greater and has more complex variation in the summer than in the winter, which makes it more difficult to accurately model the relationship between tropospheric delay in the zenith direction and slant direction. In addition, the RMS is much larger in the regions with a large PWV value (e.g., the Equator) than in regions with a small PWV value (e.g., the North Pole). The comparison shows that all three gradient models have very similar performance in the six selected regions. The adoption of the gradient models improved the accuracy of the STD. For example, in the temperate zone, the RMS error of the STD decreased from approximately 120 mm to approximately 50 mm at the 3$°$ elevation and from approximately 60 mm to approximately 20 mm at the 5$°$ elevation. Although the inclusion of a gradient model is useful for improving the accuracy of the STD, its error is still as large at 10~20 mm at a 5$°$ elevation owing to the unmodeled atmospheric inhomogeneity. In addition, although the coefficients of the VMF1 were determined using the observations at 3.2°, 5°, 7°, 10°, 15°, 20°, 30°, 50°, 70°, and 90° elevations, which will, to a certain degree, account for the bending effect, the nonlinear variation in the bending effect with the elevations still results in a large source of error in the determination of the STD using the mapping functions.

5. Conclusions

The accuracy of GNSS-derived STD is closely related to the adopted mapping functions and gradient models. In this study, the performances of various types of mapping functions (i.e., NMF, GMF, and VMF1) and gradient models (i.e., Chen, MacMillan, and Meindl) were assessed in six regions with typical climate conditions (the temperate zone, Qinghai–Tibet Plateau, Equator, Sahara Desert, Amazon Rainforest, and North Pole). The STD calculated with the ray-tracing technique from ERA-Interim is used as a reference to assess the STD calculated with ZHD, ZWD, and various mapping functions plus gradient models. Both the ZHD and ZWD are also calculated from the ERA-Interim. We first compared the STD derived with zenith delay and three types of mapping functions (i.e., NMF, GMF, and VMF1) to the STD derived with the ray-tracing method. Notably, in the first experiment, no gradient model was included in the determination of the STD. The results show that all three mapping functions perform similarly above 15° in terms of their RMS and bias, but VMF1 performs better than the other two mapping functions below a 15° elevation. This study also indicates that the performances of the mapping functions clearly depend on the season. For example, for the regions located in the Northern Hemisphere (i.e., the temperate zone, the Sahara Desert, and the Qinghai–Tibet Plateau), the RMS values and biases are larger during the summer than during the winter. This result occurs because the water vapor amount is greater and has more complex variation in the summer than in the winter, which makes it more difficult to accurately model the relationship between tropospheric delay in the zenith direction and slant direction.

In the second experiment, three types of gradient models and the VMF1 were included to determine the STD. The results show that all three types of gradient models perform similarly in the six selected regions. The inclusion of gradient models can significantly improve the accuracy of STD compared with the STD determined without gradient models. In most regions, this improvement is approximately 50%–60%, but in the Amazon rainforest, this improvement is approximately 74%. The results also show that the RMS of the STD has an obvious seasonal variation in the temperate zone, the Qinghai–Tibet Plateau, and the Sahara Desert regions. The SWD error is larger than that of the SHD in most regions.

This study shows that, although the coefficients of the VMF1 will, to a certain degree, account for the bending effect, the nonlinear variation in the bending effect with elevation still represents a large source of error in the determination of STD using the mapping functions.