Improved Empirical Coefficients for Estimating Water Vapor Weighted Mean Temperature over Europe for GNSS Applications

Abstract

Development of the so-called global navigation satellite system (GNSS) meteorology is based on the possibility of determining a precipitable water vapor (PWV) from a GNSS zenith wet delay (ZWD). Conversion of ZWD to the PWV requires application of water vapor weighted mean temperature ($T_m$) measurements, which can be done using a surface temperature ($T_s$) and its linear dependency to the $T_m$. In this study we analyzed up to 24 years (1994–2018) of data from 49 radio-sounding (RS) stations over Europe to determine reliable coefficients of the $T_m-T_s$ relationship. Their accuracy was verified using 109 RS stations. The analysis showed that for most of the stations, there are visible differences between coefficients estimated for the time of day and night. Consequently, the ETm4 model containing coefficients determined four times a day is presented. For hours other than the primary synoptic hours, linear interpolation was used. However, since this approach was not enough in some cases, we applied the dependence of $T_m-T_s$ coefficients on the time of day using a polynomial (ETmPoly model). This resulted in accuracy at the level of 2.8 ± 0.3 K. We also conducted an analysis of the impact of this model on the PWV GNSS. Analysis showed that differences in PWV reached 0.8 mm compared to other commonly used models.

1. Introduction

Earth’s atmosphere affects electromagnetic global navigation satellite system (GNSS) signals by delaying their propagation. While the ionosphere influence can be reduced using dual-frequency measurements, the delay caused by the troposphere has to be estimated during the precise positioning process. Tropospheric delay is usually expressed in the zenith direction (zenith tropospheric delay, ZTD) and it describes the total impact of neutral atmosphere on the GNSS signals. Therefore, ZTD contains delays caused by the hydrostatic (zenith hydrostatic delay, ZHD) and the wet (zenith wet delay, ZWD) part of the troposphere. Both ZHD and ZWD depend on the meteorological parameters and, consequently, changes in these parameters over time. ZHD is correlated with changes of air pressure and temperature, while ZWD depends on the amount of water vapor in the atmosphere.

The sensitivity of the ZWD parameter to the water vapor content variation has triggered the development of the so-called GNSS meteorology. It is based on the possibility of converting the ZWD to precipitable water vapor (PWV) [1], which describes the total content of water vapor in the air column above a GNSS antenna. Short-term monitoring of the water vapor variations is an important element in monitoring and studying of rapid changes of weather conditions [2,3,4]. On the other hand, estimation of the PWV long-term parameters can be used as part of investigation of the climate change [5,6,7]. Vertical PWV variations can be obtained via GNSS tomography [8]. Although this approach enables to estimate accurate humidity distribution through troposphere, it requires delivering of temperature profiles, using dense network of GNSS receivers and complex calculations. PWV can be also calculated based on the ZTD derived from GNSS processing and various external meteorological parameters. Temperature and air pressure profiles are used to estimate ZHD, which is necessary to obtain ZWD. In turn, profiles of temperature and water vapor pressure are required to calculate water vapor weighted mean temperature ($T_m$), which takes part in the conversion from ZWD to PWV. In situ measurements of the profiles can be done using radiosonde or microwave radiometers. Unfortunately, the vast majority of GNSS stations are not collocated with such instruments. Thus, it is not possible to obtain the required meteorological data. Therefore, it is necessary to apply some simplifications. ZHD can be accurately modeled using total air pressure at the GNSS antenna location [9], while $T_m$ can be obtained based on the linear relationship with the surface temperature ($T_s$): $T_m=a\cdot T_S+b$. This relationship (henceforth named $T_m-T_s$) was firstly proposed by Bevis et al. [1]. In their model, the $a$ (slope) and $b$ (intercept) coefficients were calculated on the basis of 8712 radiosonde profiles made on 13 U.S. stations over a two-year period. They determined the values of the coefficients at 0.72 and 70.2 for $a$ and $b$, respectively, with the root mean square error (RMSE) at the level of 4.74 K. Braun and Van Hove [10] used, based on a personal communication from Michael Bevis, revised version of these coefficients. Based on the profiles derived from nearly globally distributed radio-soundings stations, updated coefficients values were estimated: $a=0.668$ and $b=85.63$.

A linear regression analysis of the radiosonde profiles was also conducted by Mendes and Langley [11] who based on the data from stations located between 62°S to 83°N proposed that the following coefficients can be used: $a=0.789$ and $b=50.4$. Ross and Rosefeld [12] found that geographical localization may affect the reliability of the $T_m-T_s$ model and have proposed to use the nonlinear function ($T_m=3.402\times {10}^{-6}\cdot T_s^3+196.05$) for stations located at higher latitudes. Consequently, regional models have started to be developed. Here, Emardson et al. [13] can be mentioned, who have proposed to use a model which included information about time of the year for the Scandinavia region, or Solbrig et al. [14], who have used simple $T_m-T_s$ but estimated specifically for the territory of Germany. Despite the fact that they only included data from Germany, they obtained results ($a=0.77$, $b=54.7$) similar to those presented by Mendes and Langley on a global scale. The optimization of the $T_m-T_s$ for the specific area was also undertaken by, e.g., Liu et al. for Taiwan [15], Baltink et al. for the Netherlands [16], Bokoye et al. for Canada [17], Suresh et al. for India [18], Sappucci for Brazil [19], Mekik and Deniz for Turkey [20], Liu et al. for China [21], or Zhang et al. for Tahiti [22]. These works show the importance of developing accurate spatio-temporal models of $T_m$. Simultaneously, new global coefficients were also modeled. Based on the data from the numeric model, Schueler et al. [23] proposed three approaches: Harmonics, which took into account the day of year and seasonal variations of $T_m$, a linear $T_m-T_s$, and a combined periodic/linear assumption. Yao et al. [24] have proposed a station location (latitude, longitude, and elevation) and season-dependent coefficients for the $T_m-T_s$ linear relation. Lan et al. [25] proposed a gridded model with 2° $\times$ 2.5° window size to deliver the site-accurate $a$ and $b$ coefficients, while Yao et al. [26] developed a globally applicable global weighted mean temperature (GWMT) model which considered annual cycle, as well as semi-annual and diurnal variation of $T_m$ sourced from the Global Geodetic Observing System (GGOS).

All so far developed models have some advantages and limitations. Local models are best suited to the prevailing weather conditions occurring over a given region, while global models can be used on larger scale. On the other hand, in the case of global models, the coverage of radio-sounding (RS) stations is not uniform. This evokes the need to accept a compromise between model accuracy and global availability or the need to source data from weather models, not in situ observations. The numerical weather models, in turn, contain the average values of meteorological parameters, which can cause deterioration of the results. Until now, there is a lack of a regional model that could be used across the European continent with high accuracy. In this paper, we focus on a regional model of the $T_m-T_s$ relationship, developed exclusively on the basis of in situ measurements. In total, we used observations from 109 RS stations evenly distributed throughout Europe to deliver and validate reliable coefficients of the $T_m-T_s$. More details about the employed data are given in Section 2 together with a description of the methodology. The accuracy of the selected models is assessed in Section 3. The process of estimation of new coefficients is presented in Section 4, while the results are described and discussed in Section 5.

2. Data and Methodology

ZTD, together with the station coordinates, is a direct product of GNSS processing. It contains hydrostatic (ZHD) and wet components (ZWD):

$$ZTD=ZHD+ZWD.$$

(1)

ZWD is strictly correlated to the amount of water vapor in the atmosphere and can be converted to PWV by:

$$PWV= \Pi \left({\mathrm{T}}_m\right)\cdot \mathrm{ZWD},$$

(2)

where $\Pi$ is a dimensionless quantity dependent on ${\mathrm{T}}_m$:

$${\Pi }^{-1}={10}^{-8}\times \left(R_v\cdot \left(K_3/T_m+K_2^{\prime }\right)\right),$$

(3)

where $R_v$ is the specific gas constant of water vapor (461.5 J/kg K), and:

$$K_2^{\prime }=K_2-K_1\frac{M_w}{M_d}\approx 22.1 \left[k/hPa\right],$$

(4)

where $K_1, K_2, K_3$ denotes air refractivity parameters. Here, we applied the “best average” values estimated by Ruger [27], equal to 77.689 ± 0.0094 K/hPa, 71.295 ± 1.3 K/hPa, and 3.75463 × 10⁵ ± 0.0076 K²/hPa, respectively. $M_w=18.0151$ and $M_d=28.9644$ are the molar masses of wet and dry air expressed in g/mol.

We have stated that ${\mathrm{T}}_m$ can be obtained by vertical integration of the temperature in the atmosphere and requires measurements from a sounding station or a microwave radiometer, or data from an accurate and validated numerical weather model. The other way is to use the linear relationship with the ${\mathrm{T}}_s$. To estimate $T_m-T_s$ coefficients, it is necessary to know both values at the station location. In this study, observations from 109 RS stations located in Europe were used (Figure 1). Values of $T_s$ were derived directly from the measurements while $T_m$ were calculated by vertical integration of the temperature ($T$) and water vapor pressure ($P_w$) profiles:

$$T_m={{\displaystyle \int }}_{h_0}^{h_{max}}\frac{P_w/T}{P_w/T^2}dh,$$

(5)

where $h_0$ is a station altitude and $h_{max}$ denotes the maximum height of the atmospheric profile. We used only profiles made up to at least 8 km.

Thus obtained time series of $T_m$ and $T_s$ were divided into two groups. The first of them contained stations which have been working since 1994 to 2018, collected at least 80% of observations, and were free from visible shifts and significant gaps (henceforth named estimation group). This group consists of 49 stations and constitutes a basis for estimation new coefficients of $T_m-T_s$ (Figure 1, the red triangles) using linear regression. To ensure independent verification of our results, we selected the second group containing 60 RS stations (Figure 1, the light violet circles). These stations (henceforth named external group) have at least 2 years of observations and together with the estimation group were used as an external indicator of the accuracy of all coefficients shown in this paper. All observations were obtained from Earth System Research Laboratory (ESRL) database and were checked before analysis. Incomplete or too low (below 8 km) profiles as well as abnormal recorded values were removed from further processing.

In Section 5, we provide the impact of the new estimated coefficients on the PWV. Therefore, we had to extract ZWD from Equation (1). To this end, we used the Saastamoinen hydrostatic model [9] and surface meteorological parameters derived from the ERA5 model [28] and interpolated at the locations of the stations [29]. As we focus on PWV differences, this approach was acceptable.

3. Accuracy of the Existing ${\mathit{T}}_{\mathit{m}}\mathbf{-}{\mathit{T}}_{\mathit{s}}$ Coefficients

Here, we present the assessment of the accuracy of the $T_m-T_s$ coefficients that can be adopted in Europe. We took into consideration only these models which were estimated based on the in situ observations. In the described analysis, all 109 RS stations presented in Figure 1 were used. For each of them the $T_m$ values were estimated using Equation (5). They were compared to the $T_m$ values calculated from the relationship with $T_s$ and various coefficients. Obtained differences have allowed for an assessment which coefficients provide better accuracy and if this accuracy is evenly distributed over the entire area.

We decided to consider coefficients developed by Bevis et al. [1] (henceforth named Bevis coefficients), since they are still one of the most commonly used. Moreover, their revised versions [10] (BevisRev coefficients) as well as values presented by Mendes and Langley [11] (Mendes coefficients) and Solbrig [14] (Solbrig coefficients) were taken into account. Figure 2 presents the RMSE values for each of the analyzed station and set of coefficients. Using Bevis coefficients, we achieved mean RMSE value at the level of 3.1 ± 0.4 K. As it is seen in Figure 2a, the distribution of RMSE for the individual stations is not homogeneous. The greatest consistency can be found in the north-east area, where for most of the stations the RMSE value is below 3.0 K. This may result from the fact that these stations are located in the temperate oceanic climate zone for which the Atlantic Ocean has a smoothing effect. Similar results can be found when BevisRev coefficients were used (Figure 2b). In this case the mean RMSE was at the level 3.3 ± 0.4 K. These slightly poorer results are probably caused by the spatial distribution of stations used for coefficients estimation (regional in Bevis and global in BevisRev). Application of Mendes coefficients (Figure 2c) did not bring significant difference in results—the mean RMSE value was 3.1 ± 0.5 K. Interesting results were obtained for $T_m$ values estimated based on Solbrig coefficients (Figure 2d). Although the coefficients were estimated only for Germany, they can be used for stations located in almost all of Europe, especially the northern part. The obtained accuracy did not diverge significantly from others. The mean RMSE for this solution was 3.2 ± 0.6 K. In this case, the RMSE is affected by poor results obtained for stations located in the Mediterranean area. It is also worth to notice here, as opposed to the expectations, that the stations located in Germany are not characterized by higher accuracy. Summary of the provided results is given in

Table 1. Average root mean square error (RMSE) ($\overline{RMSE}$) of the analyzed coefficients calculated based on 109 RS stations.

${\mathit{T}}_{\mathit{m}}=\mathit{a}\cdot {\mathit{T}}_{\mathit{s}}+\mathit{b}$
Coefficients	a	b	$\overline{\mathit{R}\mathit{M}\mathit{S}\mathit{E}}$
Bevis	0.72	70.2	3.1 ± 0.4
BevisRev	0.668	85.63	3.3 ± 0.4
Mendes	0.789	50.4	3.1 ± 0.5
Solbrig	0.77	54.7	3.2 ± 0.6

.

Although part of the stations is in a very good agreement with the $T_m$ values obtained by the way of direct measurement, the limitations of these models are also clear. This motivated us to reanalyze RS data from Europe and to verify whether it is possible to estimate coefficients with higher accuracy.

4. Estimation of the New ${\mathit{T}}_{\mathit{m}}\mathbf{-}{\mathit{T}}_{\mathit{s}}$ Coefficients

In the first approach, the $T_m-T_s$ coefficients were estimated using the same method as in case of the above described studies. The $T_m$ and $T_s$ data from all stations belonging to the estimation group were analyzed in one common process, in which a linear regression was conducted to the whole dataset consisting of 967,277 records. As a result, the coefficients (named ETm) of the $T_m-T_s$ were estimated at 0.7440 ± 0.0004 and 62.84 ± 0.10, for $a$ and $b$, respectively, with RMSE equal 3.0 K.

The second approach was preceded by an investigation of the relationship between time-varying $T_m$ and $T_s$. The possible dependency between $T_m-T_s$ coefficients and the time of day were investigated, e.g., by Mekik and Deniz [19]. In their study, they proved that the differences in RMSE values obtained based on coefficients estimated during day and night are not significant, amounting to about 0.2 K, for RS stations located on the Turkish territory. To verify whether such dependency occur at European stations or not, we estimated $T_m-T_s$ coefficients for two separated synoptic hours: At 00:00 and 12:00 UTC. We found clear differences between the $T_m-T_s$ (Figure 3a–c) estimated during the night (higher RMSE values) and the day (lower RMSE values). In the case of some stations, the RMSE between night and day drops from 3.08 K to 2.19 K, 2.97 K to 2.56 K, and from 3.29 K to 2.41 K for 07761, 08430, and 14240 station, respectively. Only for a few RS stations, the obtained differences were less significant, e.g., 01415 for which RMSE dropped from 2.34 K to 2.06 K (Figure 3d).

The fact that the linear regression coefficients show dependencies on the time of day prompted us to include this fact into the estimation of the new $T_m-T_s$ model. In this approach, observations from the 00:00 and 12:00 UTC collected from all stations belonging to the estimation group were analyzed separately. We obtained the following coefficients for daytime: $a$ = 0.7430 ± 0.004 and $b$ = 61.84 ± 0.13, and for nighttime: $a$ = 0.8436 ± 0.0006 and $b$ = 35.88 ± 0.17. These coefficients, named ETm2, were estimated with RMSE of 2.46 K and 2.91 K, respectively, for the day and the night. It is seen that the obtained RMSE are lower than in case of ETm coefficients, which proves that ETm2 ensures a better fit for the linear model.

Many RS stations had observations performed also at additional hours, mostly at 06:00 and 18:00 UTC. This fact allowed for an estimation of linear regression coefficients (ETm4) also for these terms. As a result, we have obtained four coefficients of $T_m-T_s$ which are presented in

Table 2. Summary of the estimated $T_m-T_s$ coefficients.

${\mathit{T}}_{\mathit{m}}-{\mathit{T}}_{\mathit{s}} \mathbf{Coefficients}$
Name	Hours (UTC)	$\mathit{a}$	$\mathit{b}$	$\overline{\mathit{R}\mathit{M}\mathit{S}\mathit{E}}$	Number of Observations
ETm	all available	0.7440 ± 0.0004	62.84 ± 0.10	3.02	967 277
ETm2	00:00	0.8436 ± 0.0006	35.88 ± 0.17	2.91	391 401
	12:00	0.7430 ± 0.0004	61.84 ± 0.13	2.46	381 106
ETm4	00:00	0.8436 ± 0.0006	35.88 ± 0.17	2.91	391 401
	06:00	0.7997 ± 0.0014	48.07 ± 0.40	2.95	26 160
	12:00	0.7430 ± 0.0004	61.84 ± 0.13	2.46	381 106
	18:00	0.7478 ± 0.0011	61.00 ± 0.31	2.43	24 862

together with RMSE. It can be seen that the uncertainties of estimated coefficients at 06:00 and 18:00 have higher values. This is due to the fact of using fewer stations for calculation than on the other terms. Nonetheless, the differences between the coefficients are significantly higher than their uncertainties. Determination the coefficients for several terms in a day can increase the $T_m$ accuracy by making the simple $T_m-T_s$ formula time dependent. This will cause the coefficients to be interpolated at a specific time of a day.

5. Results and Discussion

In this section we provide results of validation of obtained coefficients. To this end, in situ measurements from all RS stations were used. Then, we considered the impact of our computed coefficients on PWV calculation at both RS and GNSS stations.

5.1. Accuracy of the New $T_m-T_s$ Coefficients

The accuracy of the presented $T_m-T_s$ coefficients was evaluated in the same way as previously described and expressed by RMSE values. As a reference the $T_m$ values obtained directly from RS observations were adopted. In Figure 4, the RMS values obtained for each RS station are presented. The triangles denote stations on the basis of which the coefficients were estimated. Furthermore, an external verification was done by the stations marked by the circles. For the coefficients estimated for different terms of a day (ETm2 and ETm4) we made linear interpolation to obtain values for a given hour of measurement.

The mean RMSE of the $T_m$ values calculated using ETm coefficients amounts 3.1 ± 0.5 K. As can be noticed (Figure 4a), the distribution of the RMS value for individual stations is close to those obtained for Bevis and Mendes coefficients. North-eastern Europe shows the higher consistency between estimated and measured $T_m$ value, while the Mediterranean area is affected by larger discrepancies between these values. Stations with the poorest accuracy are 14430 (Zadar–Zemunik, Croatia), 16144 (S. Pietro Capofiume Molin, Italy), and 11240 (Graz, Austria). A higher accuracy was obtained when the ETm2 coefficients were used. In this case, the RMSE was at the level of 2.8 ± 0.4 K. In contrast to the ETm, in ETm2 (Figure 4b) there is no inferior fit between measured and estimated parameters in the Mediterranean region. The use of the ETm4 coefficients resulted in the mean RMSE at the level of 2.8 ± 0.4 (Figure 4c). Although, at first sight, these results do not significantly vary from the ETm2 results, it is worth to notice that stations with the poorest (so far) accuracy have improved their results, especially station 11240 (Graz, Austria) (4.4 K to 4.1 K).

The application of linear interpolation to obtain values of coefficients at the specific epoch may introduce some uncertainties because in reality $T_m$ changes between four primary synoptic hours may not have a linear nature. Thus, we decided to express daily changes of $a$ and $b$ values using polynomial function. For this purpose, $T_m-T_s$ coefficients were estimated for each epoch when the total number of the RS observations was over 20,000. Then, we fitted polynomials of different degrees. Results are presented in Figure 5. Generally, all polynomials are well fitted in points where the observations were performed (Figure 5, black dots). Therefore, the average RMSE values calculated at all stations are very similar. The difference is visible in detailed comparison of the RS stations which performed measurements at other epochs, e.g., 11240 (Graz, Austria). For this station, the RMSE obtained using fifth degree polynomials (3.7 K) was better by about 0.1 K and 0.5 K then the results obtained using fourth and third degrees, respectively. Application of the sixth degree polynomial coefficients gives worse accuracy represented by the higher average RMSE value of about 0.2 ± 0.2 K (maximum 1.4 K) with comparison to the fifth degree polynomial. We assumed that those fifth degree coefficients allow accurate description of the time-varying $T_m-T_s$ coefficients (named ETmPoly):

$$\begin{gathered} a=-10.07t^5+23.95t^4-19.08t^3+5.998t^2-0.7914t+0.8436, \\ b=2985t^5-7200t^4+5882t^3-1923t^2+256.8t+35.87, \end{gathered}$$

(6)

where $t$ is a UTC time of the observation, expressed as a fraction of the day.

The results achieved using ETmPoly coefficients are characterized, in general, by an accuracy equal to 2.8 ± 0.3 K (Figure 4d). Most of the stations show similar results to those obtained with ETm4. The 91% of the stations are characterized by differences smaller than 0.1 K and 80% of the stations had differences smaller than 0.05 K. However, we noticed that the highest differences occurred, again, for 11240 station (Graz, Austria), which RMSE decreased from 4.1 K (ETm4) to 3.7 K (ETmPoly). This station has practically no observations at primary synoptic hours (00:00, 06:00, 12:00, and 18:00 UTC), and most of the data concerns hours between 02:00 and 04:00 UTC. Thus, although application of coefficients from the ETm4 model makes it possible to estimate $T_m$, the final accuracy suffers from the lack of observation at primary synoptic hours for which model ETm4 suits the best. Therefore, limitations of the linear interpolation of the ETm4 coefficients and benefits of introducing the non-linear ETmPoly model are easy to seen at this station.

5.2. Impact of the New Coefficients on the PWV Estimated at RS Stations

Here, an assessment of the impact of the new $T_m-T_s$ coefficients on the PWV estimation at RS stations is shown. The analysis was conducted for the ETmPoly, since these coefficients provide the highhest accuracy. PWV calculated using the Bevis coefficients are also shown as a comparison. Both series of PWV were calculated using ZWD and $T_s$ obtained at RS stations and Equations (2) and (3). Then, we compared them to the in situ measurements of PWV. Figure 6 presents the RMS for Bevis (left) and ETmPoly (right) coefficients for each of analyzed RS station. In general, the ETmPoly model is characterized by the higher accuracy in southern Europe, the slightly higher in central Europe, and similar in northern Europe, compared to Bevis model.

The above results were obtained for the entire period. Therefore, they represent the average differences and some remarkable details are not visible. In Figure 7, we present PWV differences for selected RS stations and for a representative time span. We selected stations from the areas where the differences between ETmPoly and Bevis were highest. It can be seen that characteristics of the presented variations differ. For station 15614 (Sofia, Bulgaria) a small bias, in the amount of 0.1 mm, is visible for the entire period. However, in the case of station 16144 (S. Pietro Capofiume Molin, Italy), the differences were determined by the season. In winter, the PWV obtained using ETmPoly and Bevis coefficients are almost the same, while during the summer months they are shifted by 0.2 mm. It is worth to notice that for both stations the ETmPoly solution gives more accurate results than the Bevis model.

5.3. Impact of the New Coefficients on the PWV Estimated at GNSS Stations

PWV can be calculated based on GNSS ZWD data which, in contrast to RS data, can be characterized by a much higher frequency of measurements. Here, we show an assessment of the impact of the ETmPoly coefficients on the PWV obtained from hourly ZWD values. Dataset contains ZTD values for 27 GNSS stations belonging to the EUREF Permanent Network (EPN). Tropospheric parameters were estimated using precise point positioning (PPP) method with global mapping function (GMF) [30]. Used data were described in detail in our previous study [29]. Hourly ZTD for each station was converted to the PWV according to the described methodology in Section 2.

Figure 8 shows PWV differences for selected stations. For the stations located at the southern part of the Europe, e.g., EBRE, SFER, or MEDI, the differences exceed ±0.5 mm and for the WTZR station reached 0.8 mm. For stations located at the northern part of the Europe, e.g., KIRU, the differences are smaller and contain usually in the range of ±0.2 mm. These discrepancies result from two reasons. Firstly, it is related to the accuracy of the model coefficients, which was discussed above. Secondly, it is associated with the resolution of the data. In case of GNSS technique, we usually deal with much higher data frequency than in the case of RS measurements. We assumed that the ETmPoly coefficients are linked to more accurate results than the Bevis coefficients. This results from the fact, that for the RS stations which collected observation on different terms than the 00:00, 06:00, 12:00, and 18:00 UTC, the ETmPoly allowed to get the highest accuracy of the ${\mathrm{T}}_m$ values, which was described on the basis of 11240 (Graz, Austria) station.

Like in the case of PWV estimated on the basis of RS data (Figure 7), also in the case of GNSS PWV, an increase in the discrepancy between ETMPoly and Bevis solutions took place during the summer months. To verify the spread of these differences depending on the season, we calculated the STD of the differences obtained during December/January/February (DJF) and June/July/August (JJA) seasons. The differences between GNSS PWV estimated using ETmPoly and Bevis are smaller during the DJF (Figure 9a) than during the JJA (Figure 9b) season. Moreover, the discrepancies between the results for individual stations are less significant in this case. Generally, the STDs for DJF are in a range from 0.02 mm to 0.08 mm, with slightly higher values for stations located along the coastal area, compared to the stations located at the continental part of Europe. In JJA, STDs of the differences are between 0.08 mm to 0.22 mm (WARE station) with slightly lower values for stations located at the northern part of Europe. The only exception of this is the GRAS station, which despite the location in the southern part of the continent is characterized by low STD value. This station, however, is also located at the highest altitude (1319 m) from all analyzed cases. Therefore, the smaller humidity fluctuations at this altitude, as well as at the northern part of Europe during the JJA and continental part of Europe during the DJF, supports the determination of similar $T_m$ values for both models.

6. Conclusions

This paper presents a comparative analysis of the existing $T_m-T_s$ linear coefficients together with the development of a new form. In total, 109 RS stations over Europe were used to access and validate the accuracy of the $T_m$ estimation. Firstly, the Bevis, BevisRev, Mendes, and Solbrig coefficients were examined by comparing $T_m$ time series to the in situ measurements derived from RS stations. This indicated that Bevis and Mendes have similar accuracy and spatial distribution of the obtained results. The mean RMSE were 3.1 ± 0.4 K and 3.1 ± 0.5 K, respectively. In both cases, the highest consistency between the estimated and the observed data was found in north-western Europe, while the worst results were seen for the south-east. Similar results for north-western stations were found for the Solbrig coefficients. In other cases, significantly worse results were achieved. The $T_m$ estimation using BevisRev coefficients, which is Bevis global revision, was characterized by a mean RMSE at the level of 3.3 ± 0.4 K. Although this model represents a similar distribution of the results, in general, it did not achieve RMSE values as low as in the case of previous models (especially in the north-east area). In addition, stations located at the highest latitudes were represented by much lower accuracy. This is probably a consequence of using low latitude stations and stations from the southern hemisphere for coefficient estimation.

The coefficients we propose were estimated based on 49 RS stations over that were continuously operating for at least 20 years. In total, almost 970,000 records were used to estimate the $T_m-T_s$ coefficients following three different approaches: (i) All observations were used simultaneously to estimate one set of coefficients (ETm), (ii) two set of coefficients were estimated based on a separate calculation of observations only from 0 and 12 hours (ETm2), and (iii) the coefficients were estimated at four (00:00, 06:00, 12:00, 18:00) synoptic hours (ETm4). The determination the coefficients several times a day allowed us to perform an interpolation between terms to increase the $T_m$ determination accuracy. We showed that an application of even a simple linear interpolation brings tangible benefits especially for the measurements not performed at four primary synoptic hours. This might be exemplified by differences in $T_m$ RMSE between ETm and ETm4 models. For 101 out of 109 RS stations, we noted a reduction in RMSE value. Its mean value was 0.3 K and the highest reached up to 0.8 K (16546, Piras, Italy).

As a complement to the above approaches, we decided to apply a five-degree polynomial fitting instead of a linear interpolation to estimate $a$ and $b$ coefficients for a specific epoch during daytime (ETmPoly). Although this model is based on the exact same set of data as ETm4, we noticed some differences between the obtained results. Despite the fact that the mean RMSE of $T_m$ for these two models is the same (2.8 K), the STD was reduced from 0.4 K to 0.3 K. The highest differences between ETm4 and ETmPoly were found for stations which have observations out of primary synoptic hours. Here, the 11240 (Graz, Austria) station, for which RMSE was reduced by 0.4 K, can be mentioned.

Each of the introduced models containing values of time-varying $a$ and $b$ coefficients represents higher accuracy than the Bevis model. For ETm2, ETm4, and ETmPoly we achieved accuracy at the level of 2.8 ± 0.4 K, 2.8 ± 0.4 K, and 2.8 ± 0.3 K, respectively, while in the case of Bevis model it was 3.1 ± 0.4. In addition to the RMSE, differences between the analyzed approaches can be found during estimation of the PWV value, especially based on GNSS data. We showed that although mean values of differences between, e.g., Bevis and ETmPoly coefficients may not be significant, applying these two models to the data with high temporal resolution results in discrepancies in PWV values up to 0.8 mm.

The most reliable results were obtained using ETmPoly coefficients. This is due to the fact that in addition to time-varying $a$ and $b$ coefficients, this model does not assume a linear character of the changes between main synoptic terms. Taking into account the daily variability of $T_m-T_s$ coefficients during the day by fitting the polynomial bring benefits and is especially important in GNSS meteorology, where time resolution of the observations is much higher than in case of RS.